How can I catch divergence with the Newton method?

493
views
0
9 months ago by

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

Community: FEniCS Project
Not sure about official parameters or methods, but if you implemented your own newton solver, which is relatively straight-forward, you can set any divergence criteria you like along with associated logical steps.
written 9 months ago by pf4d

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.

written 9 months ago by Alexander G. Zimmerman
You do not require any extra code to implement in parallel with fenics; this is done automatically for you.
written 9 months ago by pf4d

7
9 months ago by
As that bloke Evan said, you can do it yourself. Also don't forget to use a relaxation parameter for highly nonlinear problems.

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)

solver = CustomSolver()
problem = CustomProblem(derivative(F, u), F, bc)
solver.solve(problem, u.vector())

plot(u, title="Solution")
plt.show()​
Cheers mate, and nicely demonstrated.
written 9 months ago by pf4d

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).

written 9 months ago by Alexander G. Zimmerman
This works just as expected. Thanks!

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.
written 9 months ago by Alexander G. Zimmerman

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.

written 8 months ago by Alexander G. Zimmerman
Check the source code for details:

https://bitbucket.org/fenics-project/dolfin/src/96a53874966cd1cc9caea4a71f0a7e82a0bc591e/dolfin/nls/NewtonSolver.cpp?at=master&fileviewer=file-view-default
written 8 months ago by pf4d
0
8 months ago by

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

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