Manual newton iteration NOT converging in FEniCS

General
1

I am trying to write a manual newton iteration code in FEniCS for some problem, i.e., I want to supply my own Jacobian and residual information instead of using the derivative command. To test my algorithm in this context I am considering a simple linear poisson problem with homogeneous boundary data on a unit square domain. Then, writing u = u^0 + \delta u, I obtain the following to be solved:

(\nabla \delta u, \nabla v) = -(\nabla u^0, \nabla v) + (f, v)

In order to do this I have written the following minimal code:

from dolfin import *

# Data
N = 8; deg = 1; nr = 4
uex = Expression('sin(pi*x[0])*sin(pi*x[1])', degree = deg + 3)
f = Expression('2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1])', degree = deg + 3)
u0 = Constant(0)

# Domain
mesh = UnitSquareMesh(4)

# Function space
V = FunctionSpace(mesh, 'CG', 1)

# Boundary
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u0, boundary)    

u = TrialFunction(V)
v = TestFunction(V)
Eq1 = dot(grad(u), grad(v))*dx
Eq2 = f*v*dx

u0 = interpolate(Constant(0), V)

iter = 0; R = 1; tol = 1e-5

while (R > tol):
    J = Eq1 # Jacobian
    res = - action(Eq1, u0) + Eq2 # residual

    deltau = Function(V)
    solve(J == res, deltau, bc)

    error = (u0 - deltau)**2*dx
    R = sqrt(abs(assemble(error)))    

    iter += 1

    u0 += deltau

    print('|res| = %e, iteration %d \n' % (R, iter))

Error: 'Sum' object has no attribute 'ufl_function_space' File "/Users/jaitushar/Documents/myFEniCS/mySolvers/mySSNStAdSolver.py", line 33, in <module> res = - action(Eq1, u0) + Eq2

2 
from dolphin import *

Which library is this? It could be overriding what you expect from dolfin

Are you sure your formulation is correct? See for example Custom Newton solvers — FEniCSx tutorial.

1 Like
3

Hi nate. Sorry it was a typo, the library is dolfin only. I have updated it in the code

4

Your Newton algorithm does not seem correct. See

1 Like
5

After looking at the algorithm, I realised that I was not performing the following update $$u0 += \deltau$$. But now, after this update u0 becomes a variable of type <class \;\; 'ufl.algebra.Sum'>. and thus produces an error 'Sum' object has no attribute 'ufl_function_space' on the residual calculation step, specifically the part action(Eq1, u0). How to bypass this error

6

You want to assign data to the vector u0. i.e.
u0.vector()[:] += deltau.vector()[:]

7

Thank you, dockken!. It is updating now. I ran the code, but the residual remains stable i.e., the res that I am printing in my code remains |res| = 4.301287e-01, iteration 1 for all iterations and never crosses the tolerance value 1e-5, with u0

[0.         0.         0.         0.         0.46778505 0.
 0.         0.67186899 0.67186899 0.         0.         0.4823812
 0.95016624 0.4823812  0.         0.         0.67186899 0.67186899
 0.         0.         0.46778505 0.         0.         0.
 0.        ]

My Jacobian and residual forms are correct (which I have updated in the question also).

8

Please supply your updated code. Without it we would have to guess as to how it now looks. Please also add all output.

Also note that there are several examples on the forum on custom newton solvers, such as: Custom Newton Solver problem with Dirichlet conditions - #5 by javijv4 and

9

Updated code:

from dolfin import *

uex = Expression('sin(pi*x[0])*sin(pi*x[1])', degree = 2)
f = Expression('2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1])', degree = 2)
u0 = Constant(0)

mesh = UnitSquareMesh(4,4)

V = FunctionSpace(mesh, 'CG', 1)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u0, boundary)    

u = TrialFunction(V)
v = TestFunction(V)
Eq1 = dot(grad(u), grad(v))*dx
Eq2 = f*v*dx

u0 = interpolate(Constant(0), V)


iter = 0; R = 1; tol = 1e-5

while (R > tol) and (iter < 5):
    iter += 1

    J = Eq1
    res = - action(Eq1, u0) + Eq2

    deltau = Function(V)
    solve(J == res, deltau, bc)

    error = (u0 - deltau)**2*dx
    R = sqrt(abs(assemble(error)))    

    # assign data to the vector
    u0.vector()[:] += deltau.vector()[:]

    print(u0.vector().get_local())
    print('|res| = %e, iteration %d \n' % (R, iter))

OUTPUT:

[0.         0.         0.         0.         0.46664257 0.
 0.         0.67055373 0.67055373 0.         0.         0.48166361
 0.94830618 0.48166361 0.         0.         0.67055373 0.67055373
 0.         0.         0.46664257 0.         0.         0.
 0.        ]
|res| = 4.292924e-01, iteration 1 

[0.         0.         0.         0.         0.46664257 0.
 0.         0.67055373 0.67055373 0.         0.         0.48166361
 0.94830618 0.48166361 0.         0.         0.67055373 0.67055373
 0.         0.         0.46664257 0.         0.         0.
|res| = 4.292924e-01, iteration 2

Thank you, I will read the forum suggested by you.

10

It worked with the following code.

from dolfin import *

uex = Expression('sin(pi*x[0])*sin(pi*x[1])', degree = 2)
f = Expression('2*pow(pi, 2)*sin(pi*x[0])*sin(pi*x[1])', degree = 2)
u0 = Constant(0)

mesh = UnitSquareMesh(4,4)

V = FunctionSpace(mesh, 'CG', 1)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u0, boundary)    

u = TrialFunction(V)
v = TestFunction(V)
Eq1 = dot(grad(u), grad(v))*dx
Eq2 = f*v*dx

u0 = interpolate(Constant(0), V)


iter = 0; R = 1; tol = 1e-5

while (R > tol) and (iter < 5):
    iter += 1

    J = Eq1
    res = - action(Eq1, u0) + Eq2

    deltau = Function(V)
    solve(J == res, deltau, bc)

    error = (deltau)**2*dx
    R = sqrt(abs(assemble(error)))    

    # assign data to the vector
    u0.vector()[:] += deltau.vector()[:]

    print('|res| = %e, iteration %d \n' % (R, iter))

Turns out I was calculating R in the code the wrong way. Now it converges with 2 iterations.