Wrong behaviour on the second time the system is solved

Question

Hello everybody, I have to deal with a really weird behaviour in the resolution of a non linear mixed functions problem.

When I try to call the solve() method on the solver, the first time I get the right solution. The second time, I should expect the same residual and hence only a single newton iteration, but I get instead that the system can't converge.

What is going wrong?

import dolfin

# Creates a mesh and define the function spaces
mesh = dolfin.UnitIntervalMesh(10000)

V = dolfin.FunctionSpace(mesh, 'CG', 1)
V_n = dolfin.FunctionSpace(mesh, 'CG', 1)
W = dolfin.MixedFunctionSpace([V,V_n])

# Defines the boundary conditions
def left_boundary(x, on_boundary):
    return abs(x[0] - 0) < dolfin.DOLFIN_EPS and on_boundary

def right_boundary(x, on_boundary):
    return abs(x[0] - 1) < dolfin.DOLFIN_EPS and on_boundary

left_value = dolfin.Constant(0.)
right_value = dolfin.Constant(1.)

bcs_phi = [
    dolfin.DirichletBC(W.sub(0), left_value, left_boundary),
    dolfin.DirichletBC(W.sub(0), right_value, right_boundary),
]
left_value_n = dolfin.Constant(0.)
right_value_n = dolfin.Constant(0.)
bcs_n = [
    dolfin.DirichletBC(W.sub(1), left_value_n, left_boundary),
    dolfin.DirichletBC(W.sub(1), right_value_n, right_boundary),
]
bcs = bcs_phi + bcs_n

# Defines the variational problem
solution = dolfin.Function(W)
phi, phi_n = dolfin.split(solution)
v, v_n = dolfin.TestFunction(W)

g_1 = dolfin.inner(dolfin.grad(phi), dolfin.grad(v)) * dolfin.dx
g_2 = dolfin.inner(dolfin.grad(phi_n), dolfin.grad(v_n)) * dolfin.dx

F = g_1 + g_2
d_solution = dolfin.TrialFunction(W)
J = dolfin.derivative(F, solution, d_solution)
problem = dolfin.NonlinearVariationalProblem(F, solution, bcs, J)
solver = dolfin.NonlinearVariationalSolver(problem)

# Works correctly
solver.solve()

# If called the second time, the solve() doesn't work
solver.solve()

This is the result

Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.000e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 1.399e-10 (tol = 1.000e-10) r (rel) = 1.399e-10 (tol = 1.000e-09)
  Newton solver finished in 1 iterations and 1 linear solver iterations.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.399e-10 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 1.496e-10 (tol = 1.000e-10) r (rel) = 1.069e+00 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 1.338e-10 (tol = 1.000e-10) r (rel) = 9.564e-01 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 1.334e-10 (tol = 1.000e-10) r (rel) = 9.537e-01 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 1.373e-10 (tol = 1.000e-10) r (rel) = 9.818e-01 (tol = 1.000e-09)
  Newton iteration 5: r (abs) = 1.370e-10 (tol = 1.000e-10) r (rel) = 9.796e-01 (tol = 1.000e-09)
  ...
  Newton iteration 49: r (abs) = 1.315e-10 (tol = 1.000e-10) r (rel) = 9.399e-01 (tol = 1.000e-09)
  Newton iteration 50: r (abs) = 1.286e-10 (tol = 1.000e-10) r (rel) = 9.192e-01 (tol = 1.000e-09)
Traceback (most recent call last):
  File "try.py", line 50, in <module>
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 1.6.0
*** Git changeset:  9e4783019af77769cd36820da32e58f498cbeaa5
*** -------------------------------------------------------------------------

umberto · Answer 1 · Oct 12, 2015

Hi,

The Newton method is behaving exactly as expected.

The first time you call the Newton method you are using a zero initial guess.
The Newton methods converges after one iteration because the norm of the residual dropped of 9 order of magnitude with respect to the initial residual norm:
r_0 = 1.000e+00 , r_1 = 1.399e-10.

The second time you call the newton solver you provide the solution as initial guess, so that the norm of the residual is already extremely small. Since Newton uses iterative methods to solve the Jacobian matrix and since you are already extremely close to the solution, Newton stagnates.

Also, the problem you implemented is linear....

Wrong behaviour on the second time the system is solved

1 Answer