Implementing the Newton Method per hand

Hi everyone,

I’m working on a non linear singular problem, where I have to build and apply the Nullspace in every iteration and solve it with a PETscKrylovSolver so I have to programm the Newton Method per hand. First I’m trying to adapt the nonlinear poisson demo with a NewtonSolver to get a working algorithm and then use the PETscKrylovSolver but my assembled system doesn’t seem to update in each iteration. When using a residual criterion the first absolute error coincides with the dema, while using a incremental criterion doesn’t. Appreciate any help!

Cheers

Extra Question: How do I get the NonLinearVariationalSolver to report the iteration number, absolute and relative errors?

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

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

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

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

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

F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx
J = derivative(F, u, du)

def ConvergenceFunction(iteration, v, v0, abs_tol, rel_tol, convergence):
    absolute = norm(v, "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

OwnNewton       = False
criterion       = "incremental"
if OwnNewton:
    u = interpolate(Constant(0.0), V)
    relaxation               = 1.0
    iteration                = 0
    iteration_max            = 5
    absolute                 = 1.0
    absolute_tol             = 1.0E-6
    relative_tol             = 1.0E-6
    convergence              = False
    while iteration < iteration_max and convergence != True:
        A, b                                = assemble_system(J, -F, bc)
        solve(A, du_.vector(), b)
        if (criterion == "residual"):
            convergence, absolute, relative = ConvergenceFunction(iteration, b, absolute, absolute_tol, relative_tol, convergence)
        elif (criterion == "incremental"):
            convergence, absolute, relative = ConvergenceFunction(iteration, du_, absolute, absolute_tol, relative_tol, convergence)
        else:
            print("Convergence criterion not valid")
            sys.exit()
        #norm_u0                             = norm(u, "l2")
        #norm_du                             = norm(du_, "l2")
        u.vector()[:]                       += relaxation * du_.vector()
        #norm_u                              = norm(u, "l2")
        print("Newton iteration %d: r (abs) = %.3e (tol = %.3e) r (rel) = %.3e (tol = %.3e)" % (
         iteration, absolute, absolute_tol, relative, relative_tol))
        #print("Newton iteration %d: norm_u0 = %.3e, norm_du = %.3e, norm_u = %.3e" % (iteration, norm_u0, norm_du, norm_u))
        iteration                           += 1
else:
    solve(F == 0, u, bc, J = J, 
        solver_parameters={"newton_solver":{"relative_tolerance":1e-6, "convergence_criterion":criterion}})
    #problem = NonlinearVariationalProblem(F, u, bc, J)
    #solver  = NonlinearVariationalSolver(problem)
    #solver.parameters["newton_solver"]["convergence_criterion"] = criterion
    #solver.parameters["newton_solver"]["relative_tolerance"] = 1e-6
    #solver.parameters["newton_solver"]["report"] = True
    #solver.solve
file = File("Results/mod_nonlinear_poisson.pvd")
file << u

I can offer help on the iteration number.

That solver.solve() call that you have commented out returns a tuple, (number of iterations as int, solver did converge as bool).

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

I know you’ve already gone through all this work, but why not just use the inbuilt NewtonSolver (or PETScSNESSolver) and attach the nullspace to the matrices using petsc4py?

See, for example, here. You can attach the nullspaces at the CustomSolver.solver_setup() stage.

Thanks for the suggestion! I’m not sure how to implement it though. I’m fairly inadvance when it comes to fenics, and programming for that matter, and I haven’t worked with petsc4py before. Could you maybe elaborate a bit more?

What I tried was attaching the Nullspace to the Operator A in the solver_setup as you suggested but then I need to orthogonalize b too, which is not called on that method. My idea was to attach it after the assembly with the method J, but my newton iterations then diverge.

def build_nullspace(V):
	x = Function(V).vector()
	nullspace_basis = [x.copy() for i in range(6)]
	V.sub(0).dofmap().set(nullspace_basis[0], 1.0)
	V.sub(1).dofmap().set(nullspace_basis[1], 1.0)
	V.sub(2).dofmap().set(nullspace_basis[2], 1.0)
	V.sub(0).set_x(nullspace_basis[3], -1.0, 1)
	V.sub(1).set_x(nullspace_basis[3],  1.0, 0)
	V.sub(0).set_x(nullspace_basis[4],  1.0, 2)
	V.sub(2).set_x(nullspace_basis[4], -1.0, 0)
	V.sub(1).set_x(nullspace_basis[5], -1.0, 2)
	V.sub(2).set_x(nullspace_basis[5],  1.0, 1)
	for x in nullspace_basis:
		x.apply("insert")
	basis = VectorSpaceBasis(nullspace_basis)
	basis.orthonormalize()
	return basis
class Problem(NonlinearProblem):
    def __init__(self, a, L, bcs):
        self.L      = L
        self.a      = a
        self.bcs    = bcs
        NonlinearProblem.__init__(self)
    def F(self, b, x):
        assemble(self.L, tensor = b)
        for bc in self.bcs:
            bc.apply(b, x)
        NS = build_nullspace(V)
        NS.orthogonalize(b)
    def J(self, A, x):
        assemble(self.a, tensor = A)
        for bc in self.bcs:
            bc.apply(A)
class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                                PETScKrylovSolver(),
                                PETScFactory.instance())
    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)
        NS = build_nullspace(V)
        as_backend_type(A).set_nullspace(NS)
        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("pc_type", "ilu")
        PETScOptions.set("ksp_atol", 1.0e-10)
        PETScOptions.set("ksp_rtol", 1.0e-8)
        self.linear_solver().set_from_options()

A question for better understanding: When creating an instance of the NonlinearProblem class, what it is passed and then applied are the Boundary Conditions of the solution and not the Boundary Conditions of the solution increment. From what I undestand the boundary conditions for the system

J\left( u_\textrm{h}, v\right) \delta U = -F \left( u_\textrm{h}, v\right)

have to be that of the solution increment. Does the NewtonSolver zeroes the given Boundary Conditions when it calls my instance of NonlinearProblem?

On a second note I’m still curious as to why my NewtonSolver diverges from the built-in one.