# Correct definition of the jacobian for a nonlinearproblem

Hi, i’m defining by hand the jacobian J but there’s something wrong. The code is something like this:

V = VectorElement('Lagrange',triangle,1)
Q = FiniteElement('Lagrange',triangle,1)
element = Q*V
Z = FunctionSpace(mesh,element)
z = Function(Z)
dz = TrialFunction(Z)
(dη,du) = split(dz)
(η,u) = split(z)
(ϕ,w) = TestFunctions(Z)
z_n = Function(Z)
dz_n = Function(Z)
η_n, u_n = split(z_n)
dη_n, du_n = dz_n.split()

def norma(u):
return sqrt(u[0]**2+u[1]**2)

θ=0.5
ηc=θ*η+(1-θ)*η_n
Hc=ηc+h
U=θ*u+(1-θ)*u_n
dηc=θ*dη+(1-θ)*dη_n
dHc=dηc+h
dU=θ*du+(1-θ)*du_n

B = inner((η-η_n)/k,ϕ)*dx - inner(Hc*U,nabla_grad(ϕ))*dx + inner((u-u_n)/k,w)*dx \

η_i=project(Expression('2*ze*z0/L*(x[0]/L*cos(om1*t)-x[1]/L*sin(om1*t)-ze/2/L)', degree=2, ze=ze, z0=z0, L=L, om1=om1, t=0.0),Z.sub(0).collapse())
u_i=project(Expression(('-ze*om1*sin(om1*t)','-ze*om1*cos(om1*t)'), degree=2, ze=ze, om1=om1, t=0.0),Z.sub(1).collapse())
dη_i=project(Constant(0.0), Z.sub(0).collapse())
du_i=project(Constant((0.0,0.0)),Z.sub(1).collapse())
assign(z_n.sub(0), η_i)
assign(z_n.sub(1), u_i)
assign(dz_n.sub(0), dη_i)
assign(dz_n.sub(1), du_i)

t = dt
for m in range(num_steps):
η_D.t=t
u_D.t=t
bcη = DirichletBC(Z.sub(0), η_D, boundary)
bcu = DirichletBC(Z.sub(1), u_D, boundary)
bcs=[bcη, bcu]
problem = NonlinearVariationalProblem(B, z, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
η1, u1 = z.split()
assign(z_n.sub(0), η1)
assign(z_n.sub(1), u1)
t += dt


If i define:

J = derivative(B, z, dz)


it works perfectly, so i guess there is something wrong in my J or some space definitions, can you find what is it?
Thanks!

Hi,
I assume that terms like \eta_n are solutions of the previous time step, they will not appear in the Jacobian when you differentiate with respect to z.

That might be true, but still, taking away all of the \eta_n and u_n, i get this error in J:

ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ().


The error arises when i try to solve the problem, when i define J the code doesn’t complaint.

I did not check your full derivation but clearly you are not computing properly the directional derivatives, for instance you should have (if my understanding that h is a constant is correct)

dηc=θ*dη
dHc=dηc
dU=θ*du


You were right, it was just a mistake in the derivative computation