How can I catch divergence with the Newton method?


266
views
0
3 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 3 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 3 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 3 months ago by pf4d  

2 Answers


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

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()​
Cheers mate, and nicely demonstrated.
written 3 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 3 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 3 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 3 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 3 months ago by pf4d  
0
3 months ago by
pf4d  

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

Please login to add an answer/comment or follow this question.

Similar posts:
Search »