Implementing the Newton Method per hand

Linear Algebra
3

Update: I got a working code now. The Form and its Jacobi had to be updated at the beggining of each loop. Apart from that I wasn’t properly applying the Boundary Conditions on the solution increment.
Although my code apprently does converge to the solution, the absolute and relative errors don’t match with and it’s slower than the built-in NewtonSolver.

Built-in NewtonSolver with residual as criterion:

 Newton iteration 0: r (abs) = 1.776e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
 Newton iteration 1: r (abs) = 9.035e-03 (tol = 1.000e-10) r (rel) = 5.087e-04 (tol = 1.000e-06)
 Newton iteration 2: r (abs) = 2.718e-04 (tol = 1.000e-10) r (rel) = 1.530e-05 (tol = 1.000e-06)
 Newton iteration 3: r (abs) = 1.267e-07 (tol = 1.000e-10) r (rel) = 7.130e-09 (tol = 1.000e-06)
 Newton solver finished in 3 iterations and 3 linear solver iterations.

My NewtonSolver with residual as criterion:

Newton iteration 0: r (abs) = 4.439e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 5.642e-01 (tol = 1.000e-10) r (rel) = 1.271e-01 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.339e-02 (tol = 1.000e-10) r (rel) = 2.374e-02 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 8.456e-06 (tol = 1.000e-10) r (rel) = 6.314e-04 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 5.085e-12 (tol = 1.000e-10) r (rel) = 6.014e-07 (tol = 1.000e-06)

Built-in NewtonSolver with incremental as criterion:

Newton iteration 0: r (abs) = 3.489e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 1.016e+00 (tol = 1.000e-10) r (rel) = 2.913e-02 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.866e-02 (tol = 1.000e-10) r (rel) = 5.349e-04 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 6.931e-06 (tol = 1.000e-10) r (rel) = 1.987e-07 (tol = 1.000e-06)
Newton solver finished in 4 iterations and 4 linear solver iterations.

My NewtonSolver with incremental as criterion:

Newton iteration 0: r (abs) = 4.524e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 9.949e+00 (tol = 1.000e-10) r (rel) = 2.199e-01 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 1.806e+00 (tol = 1.000e-10) r (rel) = 1.815e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 5.229e-02 (tol = 1.000e-10) r (rel) = 2.896e-02 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 4.304e-05 (tol = 1.000e-10) r (rel) = 8.232e-04 (tol = 1.000e-06)
Newton iteration 5: r (abs) = 3.050e-11 (tol = 1.000e-10) r (rel) = 7.087e-07 (tol = 1.000e-06)

If someone could help me find where does my Solver diverges from the built in one, I’d really appreciate it. I want to be sure my Solver is properly written and that it is not conveging out of luck/simplicity of the problem but it’s applicable to every solveable problem

Cheers,

from dolfin import *
import numpy as np
import sys

class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

def FFunc(u, v):
    return (inner((1 + u**2)*grad(u), grad(v))- f*v)*dx

def JFunc(F, u, du):
    return derivative(F, u, du)

def ConvergenceFunction(iteration, v, v0, abs_tol, rel_tol, convergence):
    absolute = v.norm("l2")
    if (iteration == 0):
        absolute0 = absolute
    else:
        absolute0 = v0
    relative = absolute/absolute0
    if absolute < abs_tol or relative < rel_tol:
        convergence = True
    return convergence, absolute, relative

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1)

v   = TestFunction(V)
du  = TrialFunction(V)
u   = Function(V)
du_ = Function(V)
u_  = Function(V)

f = Expression("x[0]*sin(x[1])", degree=2)

bc      = DirichletBC(V, Constant(1.0), DirichletBoundary())
bc_du   = DirichletBC(V, Constant(0.0), DirichletBoundary())

OwnNewton       = True
criterion       = "incremental"
if OwnNewton:
    u = interpolate(Constant(0.0), V)
    bc.apply(u.vector())
    relaxation               = 1.0
    iteration                = 0
    iteration_max            = 10
    absolute                 = 1.0
    absolute_tol             = 1.0E-10
    relative_tol             = 1.0E-6
    convergence              = False
    while iteration < iteration_max and convergence != True:
        F = FFunc(u, v)
        J = JFunc(F, u, du)
        A, b                                = assemble_system(J, -F, bc_du)
        solve(A, du_.vector(), b)
        u_.vector()[:]                       = u.vector() + relaxation * du_.vector()
        u.assign(u_)
        if (criterion == "residual"):
            b = assemble(-FFunc(u,v))
            bc_du.apply(b)
            convergence, absolute, relative = ConvergenceFunction(iteration, b  , absolute, absolute_tol, relative_tol, convergence)
        elif (criterion == "incremental"):
            convergence, absolute, relative = ConvergenceFunction(iteration, du_.vector(), absolute, absolute_tol, relative_tol, convergence)
        else:
            print("Convergence criterion not valid")
            sys.exit()
        print("Newton iteration %d: r (abs) = %.3e (tol = %.3e) r (rel) = %.3e (tol = %.3e)" % (
         iteration, absolute, absolute_tol, relative, relative_tol))
        iteration                           += 1
    file = File("Results/mod_nonlinear_poisson_own.pvd")
    file << u
else:
    solve(FFunc(u, v) == 0, u, bc, J = JFunc(FFunc(u, v), u, du), 
    solver_parameters={"newton_solver":{"relative_tolerance":1e-6, "absolute_tolerance":1e-10, "convergence_criterion":criterion}})
    file = File("Results/mod_nonlinear_poisson_builtin.pvd")
    file << u