### Custom Newton solver for nonlinear problem and AttributeError: 'list' object has no attribute 'homogenize'

312
views
0
8 months ago by
I am solving a 2D nonlinear elasticity problem using the default NonlinearVariationalSolver and the custom Newton solver. This is an example to myself for how to setup a custom Newton solver to solve the nonlinear problem, similar to the previous post. However, I have two problems with my code.

First, there is an error with applying boundary condition to the incremental du:

AttributeError: 'list' object has no attribute 'homogenize'

This error is related with the code:

bcs_du = bcs.homogenize()

where bcs is a list of DirichletBC.

Second, if I manually specify the BCs of du, instead of using bcs_du = bcs.homogenize(), i.e.,
bcs_du = [ DirichletBC(V.sub(0), 0, left),
DirichletBC(V.sub(1), 0, bot),
DirichletBC(V.sub(1), 0, top) ]
The error disappears. However, the result from the the custom Newton solver is different from the default NonlinearVariationalSolver.

The minimum working code is as follows. Note that currently the code is dealing with a constant body force B(u) = [0 , 1], I am using a nonlinear solver to solve a linear elasticity problem. A nonlinear problem appears if the body force vector B(u) indeed depends on u, e.g., B(u) = [u.sub(0)^2, u.sub(1)^2 ]. The nonlinear case is actually what I want to do. But first I need to fix the problems mentioned before.

from dolfin import *
import os, shutil, math, sys, sympy
import numpy as np

# result folder
savedir = "Result/Example"
if os.path.isdir(savedir):
    shutil.rmtree(savedir)

# Create mesh and define function space
mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 2.0), 25, 50, "right/left")

d = mesh.topology().dim()
e1 = [Constant([1.,0.]),Constant((1.,0.,0.))][d-2]
e2 = [Constant([0.,1.]),Constant((0.,1.,0.))][d-2]

# Mark boundary subdomains
class left(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[0], 0.)
class bot(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[1], 0.)
class top(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[1], 2.)
left, bot, top = left(), bot(), top()
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
bot.mark(boundaries, 2)
top.mark(boundaries, 3)
ds = Measure("ds", subdomain_data=boundaries)

# Define functions
V = VectorFunctionSpace(mesh, "Lagrange", 1)
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration

# Define Dirichlet boundary
bc_l = DirichletBC(V.sub(0), 0, left)
bc_b = DirichletBC(V.sub(1), 0, bot)
bc_t = DirichletBC(V.sub(1), 0.2, top)
bcs = [ bc_l, bc_b, bc_t ]

# Elasticity parameters
mu, lmbda = 1.0, 1.0

# Strain and stress
def eps(u):
     return sym(grad(u))

def sigma(u):
     return 2.0*mu*(eps(u)) + lmbda*tr(eps(u))*Identity(d)

def B(u):
     uy = u.sub(1)
     return as_vector( [ 0.0, 1 ] )

T = Constant((0.0, 0.0)) # Traction force on the boundary

file_u = File(savedir+"/u.pvd");
file_s = File(savedir+"/stress.pvd")

# Total potential energy
Pi = 0.5*inner(sigma(u),eps(u))*dx - dot(B(u), u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
Jac = derivative(F, u, du)

u_init = Expression(("0", "0"), degree=1)
u.interpolate(u_init)

abs_tol, rel_tol = 1e-10, 1e-8
abs_res, rel_res = 1.0, 1.0
max_iter = 100

# bcs_du = [ DirichletBC(V.sub(0), 0, left),
#                   DirichletBC(V.sub(1), 0, bot),
#                   DirichletBC(V.sub(1), 0, top) ]
bcs_du = bcs.homogenize()
u_inc = Function(V)
nIter = 0

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

if flag:    print ("Default solver")
     problem = NonlinearVariationalProblem(F, u, bcs, Jac)
     solver = NonlinearVariationalSolver(problem)
     solver.solve()
else:    print ("Custom solver")
     while abs_res > abs_tol and nIter < max_iter:  # Newton iterations
        nIter += 1
         abs_res_old = abs_res
         A, b = assemble_system(Jac, -F, bcs_du)
         solve(A, u_inc.vector(), b) # Determine step direction
         u_inc_norm = np.linalg.norm(u_inc.vector().array(), ord = 2)
         a = assemble(F)
         for bc in bcs_du:
             bc.apply(a)
         abs_res = b.norm('l2')
         rel_res = abs_res/abs_res_old
         alpha = 1.0 # step size, initially 1
         u.vector()[:] += alpha*u_inc.vector() # New u vector
         string = " Newton iteration %d: r (abs) = %.3e (tol = %.3e) " \
         +"r (rel) = %.3e (tol = %.3e)"
         print string % (nIter, abs_res, abs_tol, rel_res, rel_tol)
file_u << u
Community: FEniCS Project

2
8 months ago by
bc0 = [DirichletBC(bc) for bc in bcs]       # List of copies of your original bcs
[bc.homogenize() for bc in bc0]             # homogenize every item (DirichletBC object)​

The list bc0 contains the homogenized versions of your bcs.

An idea as to why your residuals do not match:
The relative residual in the FEniCS NewtonSolver method is the residual in the current iteration divided by the residual in the zeroth iteration, before anything is solved. So basically the residual, your initial guess produces. If I read your code correctly you define your relative residual as the ratio of the current residual and the residual from the previous iteration.

[ [ [ [ [ [ [ [ i agree ] ] ] ] ] ] ] ]
written 8 months ago by pf4d
I found the mismatch of the result from the two solver disappears if I delete the non-zero boundary condition bc_t on the top boundary. They match perfectly. But if I add the bc_t boundary, the result of the Newton solve is not correct. I have no idea why this happens. Could you help me with this?
written 8 months ago by Wilhelm
1
It looks like you need to apply the boundary condition to the initial guess before entering Newton.
written 8 months ago by pf4d
1
Beat me to it... damn :)
written 8 months ago by klunkean
1
[ [ [ [ [ [ [ [ [ [ [ [ i know right?! ] ] ] ] ] ] ] ] ] ] ] ] ]
written 8 months ago by pf4d
I get what you mean. The initial guess should satisfy the boundary condition.
written 8 months ago by Wilhelm
1
well, the boundary data in the examples where you cut your code:

https://fenicsproject.org/qa/536/newton-method-programmed-manually

are all homogeneous, but it should have been included anyway because the use of homogenize really doesn't make any sense without it.

You have to apply the boundary condition to the unknown before entering newton, because the direction of decent is set to zero along Dirichlet boundaries, meaning that the value of the unknown will not change as it is incremented down the gradient.
written 8 months ago by pf4d
yeah. now it matches perfectly. cheers.

Thank you, pf4d.
written 8 months ago by Wilhelm
Thank you for your comment, klunken. Your suggestion works for my first question. The following two statements are equivalent to each other, right?
 bcs_du = [ DirichletBC(V.sub(0), 0, left),
    DirichletBC(V.sub(1), 0, bot),
    DirichletBC(V.sub(1), 0, top) ]

and
bcs_du = [DirichletBC(bc) for bc in bcs] # List of copies of your original bcs
[bc.homogenize() for bc in bcs_du] # homogenize every item
.

Regarding to the second question, the criterion of relative residual is not used for the convergence check in my code, although I computed it ("incorrectly". Thank you for pointing out how the relative residual is computed in Fenics, I do not know that before). I apply the homogenized boundary conditions to du, as you indicated. However, the result of the Newton solver is still not correct. So I think the mismatach of the two solvers are not due to the relative residual convergence check.

I found the mismatch of the result from the two solver disappears if I delete the non-zero boundary condition bc_t on the top boundary. They match perfectly. But if I add the bc_t boundary, the result of the Newton solve is not correct. I have no idea why this happens. Could you help me with this?
written 8 months ago by Wilhelm
1
I can't find a line where you apply your actual (non-homogenous) Dirichlet boundary conditions to your solution vector u and I think there's your problem.
First I would delete the lines
u_init = Expression(("0", "0"), degree=1)
u.interpolate(u_init)​

since the initialization of u already writes zeros too every entry.

[bc.apply(u.vector()) for bc in bcs]

- where bcs are your original non-homogenous Dirichlet boundary conditions - somewhere before your iteration loop.

The only reason you have to set homogenous Dirichlet boundary conditions on your increment is that you suppose that your initial guess already satisfies all Dirichlet boundary conditions. When the initial guess gets incremented in a Newton iteration, you'll want to add a zero in every entry where Dirichlet boundary conditions are prescribed otherwise you would change the values that originally satisfied the BC - thus homogenous DirichletBC on the increment. But for this to work, your initial u vector has to satisfy the Dirichlet BC in the first place, hence the line above.

written 8 months ago by klunkean
Thank you, klunkean.

I change my initial guess to
u_init = Expression(("0", "0.1*x[1]"), degree=1)
u.interpolate(u_init)

to meet the boundary conditions.
written 8 months ago by Wilhelm
Oh and your code statements are equivalent only for your special case obviously, but in that special case they return the same results.
written 8 months ago by klunkean
what do you mean by "your code statements are equivalent only for your special case obviously"?
written 8 months ago by Wilhelm
bcs_du = [ DirichletBC(V.sub(0), 0, left),
DirichletBC(V.sub(1), 0, bot),
DirichletBC(V.sub(1), 0, top) ]

is an explicitly given special case, tailored to the non-homogenous BC

# Define Dirichlet boundary
bc_l = DirichletBC(V.sub(0), 0, left)
bc_b = DirichletBC(V.sub(1), 0, bot)
bc_t = DirichletBC(V.sub(1), 0.2, top)
bcs = [ bc_l, bc_b, bc_t ]

whereas the expression I suggested works for any case of Dirichlet BC. Thus I'd prefer the second variant so I can simply change the non-homogenous BC and the homogenized version adapts accordingly.

written 8 months ago by klunkean
Absolutely. You are right.
written 8 months ago by Wilhelm