Evaluating the residual and Jacobian from a NonlinearVariationalProblem

Question

Hello,
How do you obtain the residual and jacobian from a NonlinearVariationalProblem? I realize there are functions residual_form and jacobian_form that return the forms, and that I can then use the action and assemble_system functions to create the jacobian matrix and residual vector. However, below is an example that makes me doubt that this is the correct way to do this. The following code solves a Monge-Ampere PDE using the built-in Newton's method. I then attempt to evaluate the residual at the solution and print out the L2 norm of the residual:

#Define the weak residual F
(Sxx, Sxy, Syy, u) = TrialFunction(MixedV)
(muxx, muxy, muyy, v) = TestFunction(MixedV)

F = inner(Sxx,muxx)*dx + 2*inner(Sxy,muxy)*dx + inner(Syy,muyy)*dx;
F += inner(Dx(u,0), Dx(muxx,0))*dx + inner(Dx(u,0), Dx(muxy,1))*dx;
F += inner(Dx(u,1), Dx(muxy,0))*dx + inner(Dx(u,1), Dx(muyy,1))*dx;

F += ( inner(Dx(Sxx,0), Dx(v,0)) + inner(Dx(Sxy,0), Dx(v,1)))*dx;
F += ( inner(Dx(Sxy,1), Dx(v,0)) + inner(Dx(Syy,1), Dx(v,1)))*dx;

F += (((Sxx*Syy - Sxy*Sxy)-(1 + (Dx(u,0)**2 + Dx(u,1)**2))**(2)) * K)*v*dx;

F -= (-gy*muxy*ds(1) + gx*muxy*ds(2) + gy*muxy*ds(3) - gx*muxy*ds(4));

# Solve the problem
initial = Function(MixedV)
R = action(F,initial);
DR = derivative(R, initial);
problem = NonlinearVariationalProblem(R,initial,bc,DR);
solver = NonlinearVariationalSolver(problem);
solver.solve();

#Evaluate the residual
Fw = action(F,initial);
A,b = assemble_system(DR,Fw,bc);
print "My computed residual = ", b.norm('l2')

The output for this program is the following, with set_log_level set high enough to output Newton iteration residual calculations:

Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.149e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.535e-03 (tol = 1.000e-10) r (rel) = 2.207e-04 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 4.556e-09 (tol = 1.000e-10) r (rel) = 3.965e-10 (tol = 1.000e-09)
Newton solver finished in 2 iterations and 2 linear solver iterations.
My computed residual =  29.4316515057

So the built-in Newton stops because the residual is 4.556e-09 but my calculated residual is 29.43165. My question then is am I calculating the residual correctly, and if not, what is the correct way?

Answer 1

Edit: Maybe my answer wasn't very clear. Edited.
The line

A,b = assemble_system(DR,Fw,bc);

imposes boundary conditions on A and b. So you take the norm of a vector with Dirichlet-BCs for the system A*x = b incorporated. I don't know what your BCs are, but you should impose zero boundary conditions. The system you assembled is in terms of the increment, and given an initial guess satisfying your Dirichlet BCs, zero BCs for the Newton update are justified.
Setting

bc.homogenize()
b = assemble(Fw)
bc.apply(b)
print('residual: {0}'.format(b.norm('l2')))

should give you the same result. Also make sure you to set parameters["newton_solver"]["convergence_criterion"] = "residual".

Comments

I agree, but in looking at this tutorial
http://fenicsproject.org/documentation/tutorial/nonlinear.html

it looks like when setting up a NonlinearVariationalProblem, you should use the problem boundary conditions, not homogeneous ones.

---

But does it work?
It's true that you pass the problem's bcs to NonlinearVariationProblem. Wondered about this, too, so I checked out the code :) During assembly (SystemAssembler.cpp -> ll. 250), the initial solution is substracted from the bc values before they are applied. Like doing

b = assemble(F) 
bc.apply(b, x_init.vector())