Correct definition of the jacobian for a nonlinearproblem

variational formulation
1

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 \
    +inner(dot(U,grad(U))+f*R*U-1/Hc*(c_d*norma(U)*U),w)*dx - g*ηc*div(w)*dx + α_h*inner(grad(U),grad(w))*dx
J = inner((dη-dη_n)/k,ϕ)*dx - inner(dηc*U,grad(ϕ))*dx - inner(dHc*dU,grad(ϕ))*dx + inner((du-du_n)/k,w)*dx \
    +inner(dot(dU,grad(U)) + dot(U,grad(dU)) + f*R*dU + c_d/Hc*(dot(U,dU)/norma(U)*U+norma(U)*dU),w)*dx \
    -g*dηc*div(w)*dx + α_h*inner(grad(dU),grad(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!

2

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.

3

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.

4

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
5

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