This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Newton method programmed manually

+5 votes

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

asked Jul 12, 2013 by mwelland FEniCS User (8,410 points)

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$

1 Answer

+2 votes
 
Best answer

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.

answered Jul 14, 2013 by mwelland FEniCS User (8,410 points)
selected Jul 14, 2013 by Jan Blechta

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?

...