Newton method programmed manually

Question

Hi all, This is probably a simple thing but I've been hitting my head against it. I need to program the newton method manually, and want to compare it to the nonlinearvariationalsolver. Both converge but I can't seem to match my residual norms. Here is my code:

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)
        fnorm = np.linalg.norm(a.array(), ord = 2)
        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()

Change the True / False above to switch methods. With the nonlinearvariationalsolver I get the output:

Solving nonlinear variational problem.
  Newton iteration 1: r (abs) = 5.898e-03 (tol = 1.000e-10) r (rel) = 5.927e-02 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 6.381e-06 (tol = 1.000e-10) r (rel) = 6.413e-05 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 7.288e-12 (tol = 1.000e-10) r (rel) = 7.325e-11 (tol = 1.000e-09)
  Newton solver finished in 3 iterations and 3 linear solver iterations.

With my code I get:

   1  6.24E-01  2.236057e+00
   2  1.56E-02  2.121959e+00
   3  1.14E-05  2.121320e+00
   4  7.20E-12  2.121320e+00

Can anyone point out my error?

Thanks,
Mike

Comments

In your second column there is a norm of the increment while NewtonSolver checks convergence using norm of the equation residual by default. You can change this behaviour by

solver.parameters['newton_solver']['convergence_criterion'] = 'incremental'

---

Thanks Jan, but I am interested in how to calculate the norm of the equation residual actually. I thought thats what I was doing with the commands:

a = assemble(F)
fnorm = np.linalg.norm(a.array(), ord = 2)

But it doesn't seem to match the output from the newton solver... Maybe I have teh wronge equation?

---

There has been a little update on how and when the norms are calculated in Dolfin's Newton method, cf. https://bitbucket.org/fenics-project/dolfin/commits/0345014b7cd09bff6bbd8f56665a2832b13ddebb. You may want to build from Git master.

---

Nico, this is barely significant to the problem. The problem actually is that mwelland's home-made Newton solver does not show convergence in terms of $||F||_{\ell^2}$. He gets still $||F||_{\ell^2}\approx 2$

Answer 1

Thanks for the comments guys, looks like my problem was not applying the bcs_du to the matrix 'a' before calculating the norm. I tried just b.norm (inspired from nschole's notes) and it matches the nonlinearvariationalsolver output.

Comments

Yeah, it was tricky for me too to realize that assemble(F) needs application of bcs_du because test space of Dirichlet problem F is constrained by bcs_du ;-)

---

can you clarify exactly what you did to make it work, i.e. code?

---

I guess

...

a = assemble(F)
bcs_du.apply(a) # this line added
fnorm = norm(a, 'l2') # better than numpy function
lmbda = 1.0     # step size, initially 1

...

---

actually, I just replaced the calculation of fnorm to
fnorm = b.norm('l2')
(which avoids having to assemble F again too!) Here is modified code with a bit more description:

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

Hope this helps!

---

perfect, thanks!

---

I'm having an issue with the solver when solving a system of equations :

from dolfin import *
from pylab  import plot, show, spy
import sys

bb       = bool(int(sys.argv[1]))

mesh     = IntervalMesh(4, 0, 1)  
Q        = FunctionSpace(mesh, "Lagrange",1)       # Order 1 function space
MQ       = Q * Q

dU       = TrialFunction(MQ)
dU1, dU1 = split(dU)
test_U   = TestFunction(MQ)
t1, t2   = split(test_U)
U        = Function(MQ)
u1, u2   = split(U)

left, right = compile_subdomains([
  "(std::abs( x[0] )     < DOLFIN_EPS) && on_boundary",
  "(std::abs( x[0]-1.0 ) < DOLFIN_EPS) && on_boundary"])

bcs = [DirichletBC(MQ.sub(0), 1, left), 
       DirichletBC(MQ.sub(0), 0, right),
       DirichletBC(MQ.sub(1), 0, left),
       DirichletBC(MQ.sub(1), 1, right)]

k   = 1.0*u1 + 1.0

vnorm    = sqrt(dot(u2, u2) + 1e-10)
cellh    = CellSize(mesh)
t2hat    = t2 + cellh/(2*vnorm)*dot(u2, t2.dx(0))

F1  = + u1 * u1.dx(0) * t1.dx(0) * dx \
      + u1.dx(0) * t1.dx(0) * dx
c   = 1e-9
F2  = c * u2 * t2hat * dx
#F2  = u2 * t2 * dx
F   = F1 + F2

u1_init = Expression("1-x[0]")
u2_init = Constant(1)

U_k = project(as_vector([u1_init, u2_init]), MQ) # project inital values
U.vector().set_local(U_k.vector().array())       # initalize u1, u2 in solution

Jac     = derivative(F, U, dU)

if bb:
  problem = NonlinearVariationalProblem(F, U, bcs, Jac)
  solver  = NonlinearVariationalSolver(problem)
  solver.solve()
else:
  atol, rtol = 1e-7, 1e-10                      # abs/rel tolerances
  lmbda      = 1.0                              # relaxation parameter    
  U_inc      = Function(MQ)                     # residual
  bcs_u      = homogenize(bcs)                  # residual is zero on boundary
  nIter      = 0                                # number of iterations
  residual   = 1                                # residual
  rel_res    = residual                         # relative residual
  maxIter    = 100                              # max iterations

  while residual > atol and rel_res > rtol and nIter < maxIter: 
    nIter  += 1                                 # increment interation
    A, b    = assemble_system(Jac, -F, bcs_u)   # assemble system
    solve(A, U_inc.vector(), b)                 # determine step direction
    rel_res = U_inc.vector().norm('l2')         # calculate norm

    # calculate absolute residual :
    a = assemble(F)
    for bc in bcs_u:
      bc.apply(a)
    residual = b.norm('l2')

    U.vector()[:] += lmbda*U_inc.vector()       # New u vector

    string = "Newton iteration %d: r (abs) = %.3e (tol = %.3e) " \
             +"r (rel) = %.3e (tol = %.3e)"
    print string % (nIter, residual, atol, rel_res, rtol)

xcrd = mesh.coordinates()[:,0]
plot(xcrd, project(u1, Q).vector().array())
plot(xcrd, project(u2, Q).vector().array())

show()

it appears that the u2 solution's boundary conditions are not being applied properly. Any ideas?