### How can I catch divergence with the Newton method?

I am solving a nonlinear problem with

```
solver = fenics.NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['maximum_iterations'] = nlp_max_iterations
solver.parameters['newton_solver']['relative_tolerance'] = nlp_relative_tolerance
solver.parameters['newton_solver']['error_on_nonconvergence'] = False
iteration_count, converged = solver.solve()
```

For some problems I'm trying to solve, the Newton method will diverge rapidly, e.g.

(sorry I can't find a good way to format this with AllAnswered)

"

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Solving nonlinear variational problem.

Newton iteration 0: r (abs) = 1.069e+07 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 1: r (abs) = 1.622e+06 (tol = 1.000e-10) r (rel) = 1.517e-01 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 2: r (abs) = 1.327e+05 (tol = 1.000e-10) r (rel) = 1.242e-02 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 3: r (abs) = 1.617e+06 (tol = 1.000e-10) r (rel) = 1.513e-01 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 4: r (abs) = 2.822e+07 (tol = 1.000e-10) r (rel) = 2.640e+00 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 5: r (abs) = 5.849e+09 (tol = 1.000e-10) r (rel) = 5.472e+02 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 6: r (abs) = 1.140e+11 (tol = 1.000e-10) r (rel) = 1.067e+04 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 7: r (abs) = 1.903e+16 (tol = 1.000e-10) r (rel) = 1.780e+09 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

Newton iteration 8: r (abs) = 4.308e+23 (tol = 1.000e-10) r (rel) = 4.030e+16 (tol = 1.000e-04)

Solving linear system of size 8244 x 8244 (PETSc LU solver, umfpack).

PETSc Krylov solver starting to solve 8244 x 8244 system.

"

What I am missing is an option for a divergence tolerance. Something like this exists in SNES; but so far I'm less successful with the 'snes' method than with 'newton', and I don't want to make such a large change right now. So I would like to be able to write, e.g.

`solver.parameters['newton_solver']['divergence_tolerance'] = 1.e4`

Does such an option indeed not exist in fenics?

Thanks,

Alex

### 2 Answers

Adapt to fit your needs (based on nonlinear Poisson demo):

```
import matplotlib.pyplot as plt
from dolfin import *
class CustomProblem(NonlinearProblem):
def __init__(self, a, L, bcs):
NonlinearProblem.__init__(self)
self.a = a
self.L = L
self.bcs = bcs
def F(self, b, x):
assembler = SystemAssembler(self.a, self.L, self.bcs)
assembler.assemble(b, x)
def J(self, A, x):
assembler = SystemAssembler(self.a, self.L, self.bcs)
assembler.assemble(A)
class CustomSolver(NewtonSolver):
def converged(self, residual, problem, iteration):
# Some custom convergence criterion here
rnorm = residual.norm("l2")
if rnorm < 1e-6:
print("%d Residual norm ;) %.3e" % (iteration, rnorm))
return True
print("%d Residual norm :( %.3e" % (iteration, rnorm))
return False
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1)
g = Constant(1.0)
bc = DirichletBC(V, g, DirichletBoundary())
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])", degree=2)
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx
solver = CustomSolver()
problem = CustomProblem(derivative(F, u), F, bc)
solver.solve(problem, u.vector())
plot(u, title="Solution")
plt.show()
```

Thanks for the reply. I'm trying to implement this now.

How would I obtain the relative residual in this scope? Looking at the doc for NewtonSolver, I thought I should be able to do

```
def converged(self, residual, problem, iteration):
rnorm = self.relative_residual()
...
```

but this seems to always be nan (while it is a reasonable number when I use your code as is).

Now I'm just missing the relative_residual option mentioned in my previous comment. I'm not sure about etiquette here. Do I accept the answer now or after we've resolved the relative_residual question? I suppose that could be a new question entirely.

So right now I'm using the answer almost verbatum, with

```
class CustomNewtonSolver(fenics.NewtonSolver):
custom_parameters = {'divergence_threshold': nlp_divergence_threshold}
def converged(self, residual, problem, iteration):
rnorm = residual.norm("l2")
print("Newton iteration %d: r (abs) = %.3e (tol = %.3e)" % (iteration, rnorm,
self.parameters['absolute_tolerance']))
assert(rnorm < self.custom_parameters['divergence_threshold'])
if rnorm < self.parameters['absolute_tolerance']:
return True
return False
```

In some situations having this assertion is definitely better than what I had before; but what I really need is for the Newton solver to return with converged = False. This way my code can continue to handle the situation. I don't see any way to do this. Ideas? Thanks.

https://bitbucket.org/fenics-project/dolfin/src/96a53874966cd1cc9caea4a71f0a7e82a0bc591e/dolfin/nls/NewtonSolver.cpp?at=master&fileviewer=file-view-default

If you need, you can simply code the algorithm from scratch. Then you can just do whatever you want with your logic.

```
from dolfin import *
import numpy as np
mesh = IntervalMesh(100, 0,1)
V = FunctionSpace(mesh, "Lagrange",1) # Order 1 function space
dU = TrialFunction(V)
test_U = TestFunction(V)
U = Function(V)
left, right = compile_subdomains([
"(std::abs( x[0] ) < DOLFIN_EPS) && on_boundary",
"(std::abs( x[0]-1.0 ) < DOLFIN_EPS) && on_boundary"])
bcs = [
DirichletBC(V, 1, left),
DirichletBC(V, 0 , right),
]
k = 1.0*U+1.0
F = inner(k*grad(U), grad(test_U))*dx
Jac = derivative(F, U, dU)
U_init = Expression("1-x[0]")
U.interpolate(U_init)
if False:
problem = NonlinearVariationalProblem(F, U, bcs, Jac) #,form_compiler_parameters=ffc_options
solver = NonlinearVariationalSolver(problem)
solver.solve()
else:
a_tol, r_tol = 1e-7, 1e-10
bcs_du = homogenize(bcs)
U_inc = Function(V)
nIter = 0
eps = 1
while eps > 1e-10 and nIter < 10: # Newton iterations
nIter += 1
A, b = assemble_system(Jac, -F, bcs_du)
solve(A, U_inc.vector(), b) # Determine step direction
eps = np.linalg.norm(U_inc.vector().array(), ord = 2)
a = assemble(F)
for bc in bcs_du:
bc.apply(a)
print 'b.norm("l2") = ', b.norm('l2'), 'np.linalg.norm(a.array(), ord = 2) = ', np.linalg.norm(a.array(), ord = 2)
fnorm = b.norm('l2')
lmbda = 1.0 # step size, initially 1
U.vector()[:] += lmbda*U_inc.vector() # New u vector
print ' {0:2d} {1:3.2E} {2:5e}'.format(nIter, eps, fnorm)
plot(U)
interactive()
```

From our old discussion :

https://fenicsproject.org/qa/536/newton-method-programmed-manually

I wrote my own Newton solver that I was using for a while; but it is only serial and I have no desire to write a parallel (shared and distributed memory) solver.

It seems to me that the fenics Newton solver has all of the information needed to check for divergence. It just doesn't have an option to actually do it.