Newton solver, nonlinear problem convergence

1

Hello, I am working on a problem where the norm of a stress tensor is involved in the weak form, the weak form is solved on a 2D rectangular domain, where the left and right boundaries have prescribed values for the vector u and the scalar P. I have been having trouble to achieve convergence, see following an example code:

from dolfin import *
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
solver.parameters["report"] = False


class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 0.0    # u0-x all 0
        values[1] = 0.0    # u0-y all 0
        values[2] = 0.0    # P0 all 0
    def value_shape(self):return (3,)

class Problem(NonlinearProblem):
    def __init__(self, L, a, bc_u_left, bc_u_right, bc_P_left, bc_P_right):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bc_u_left = bc_u_left
        self.bc_u_right = bc_u_right
        self.bc_P_left = bc_P_left
        self.bc_P_right = bc_P_right
    def F(self, b, x):
        assemble(self.L, tensor=b)
        self.bc_u_left.apply(b,x)
        self.bc_u_right.apply(b,x)
        self.bc_P_left.apply(b,x)
        self.bc_P_right.apply(b,x)

    def J(self, A, x):
        assemble(self.a, tensor=A)
        self.bc_u_left.apply(A)
        self.bc_u_right.apply(A)
        self.bc_P_left.apply(A)
        self.bc_P_right.apply(A)


mesh = RectangleMesh(Point(0.,0.),Point(1.0,1.0),20,20,diagonal="crossed")
P1 = VectorElement("CG", mesh.ufl_cell(), 2)
P2 = FiniteElement("CG", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, MixedElement( [P1, P2] ))

dSol  = TrialFunction(ME)
testfunc = TestFunction(ME)
v,w = split(testfunc)
Sol   = Function(ME)
du, dP = split(dSol)
u,  P  = split(Sol)
Sol_init = InitialConditions(degree=2)
Sol.interpolate(Sol_init)

def eps(u):return sym(nabla_grad(u))

# weak forms

sigma0 = 77.*eps(u) 
sigma  = sigma0 + 58.**div(u)*Identity(2)
L0 = inner( sigma, eps(v) )*dx

sigma0_norm = sqrt(inner(sigma0,sigma0))
L1 =   (0.06* (1.-exp(-300.*P)) + 0.05 - sigma0_norm)*w*dx

L  = L0 + L1
a = derivative(L, Sol, dSol)

# Prescribed Boundaries
def left(x, on_boundary):return near(x[0],0.)
def right(x, on_boundary):return near(x[0],1.)

g_P  = Constant(0.0)
bc_P_left  = DirichletBC(ME.sub(1), g_P, left)
bc_P_right = DirichletBC(ME.sub(1), g_P, right)

g_u  = Constant((0.0,0.0))
bc_u_left  = DirichletBC(ME.sub(0), g_u, left)

Ubc   = 0.0
u_bc  = Constant((Ubc,0.0))
bc_u_right  = DirichletBC(ME.sub(0), u_bc, right)

problem = Problem(L, a, bc_u_left, bc_u_right, bc_P_left, bc_P_right)

u_bc.assign(Constant((0.0001,0.0)))
Sol0.vector()[:] = Sol.vector()
solver.solve(problem, Sol.vector())

I think the convergence difficulty is mainly with the term sigma0_norm, I am not sure how to go from here, does anyone have some insight on this problem? Thank you!

Traceback (most recent call last):

File "z.py", line 92, in <module>

solver.solve(problem, Sol.vector())

RuntimeError:

*** -------------------------------------------------------------------------

*** DOLFIN encountered an error. If you are not able to resolve this issue

*** using the information listed below, you can ask for help at

***

*** fenics-support@googlegroups.com

***

*** Remember to include the error message listed below and, if possible,

*** include a *minimal* running example to reproduce the error.

***

*** -------------------------------------------------------------------------

*** Error: Unable to solve nonlinear system with NewtonSolver.

*** Reason: Newton solver did not converge because maximum number of iterations reached.

*** Where: This error was encountered inside NewtonSolver.cpp.

*** Process: 0

***

*** DOLFIN version: 2019.1.0

*** Git changeset: 15b823f2c45d3036b6af931284d0f8e3c77b6845

*** -------------------------------------------------------------------------
2

The following looks like a typo to me:

sigma  = sigma0 + 58.**div(u)*Identity(2)

I assume you meant to type 58.*div(u) instead of 58.**div(u). Also, it’s generally safer in conjunction with automated symbolic differentiation to regularize norms by a small amount, e.g.,

sigma0_norm = sqrt(inner(sigma0,sigma0) + DOLFIN_EPS)

to avoid dividing by zero in the Jacobian when inner(sigma0,sigma0) is exactly zero (as in the initial condition).

Another note: It’s generally a good idea to wrap numerical literals as Constants. That way, when you change their values, you don’t need to re-compile.

1 Like
3

Hello David,
Yes, the double multiplication was a typo, and very good catch on the division by zero!! looks like it is working now for the example code, I will implement this on my full code and see how it works! Thank you!!!

4

David,

I have two follow-up questions, maybe you have some ideas on how to deal with them:

  1. In my full code, I have a tensor term in my weak form that starts as a zero tensor, and gets updated every time step based on the solutions at each time step, see the following pseudo-code:
T = as_tensor([[0.0,0.0],[0.0,0.0]])
Weak form <- T
for t in ts:
    sol <- solve weak form
    update T with sol

I tried to use T.assign(...) to update but I received this error:

AttributeError: ‘Zero’ object has no attribute ‘assign’

  1. I receive this error regarding integration, should I be worried about it? or how should I use quadrature to fix it?

WARNING: The number of integration points for each cell will be: 196
Consider using the option ‘quadrature_degree’ to reduce the number of points

5

In my full code, I have a tensor term in my weak form that starts as a zero tensor, and gets updated every time step based on the solutions at each time step

For this case, you should make T a Function in some TensorFunctionSpace, the appropriate choice of which would depend on the details of the problem. A typical choice for history variables in solid mechanics would be a function space of type "Quadrature". There is an example of doing this with a scalar-valued function space in @bleyerj’s tutorial here. However, to use quadrature elements effectively, you want to integrate all terms with the same quadrature rule, which segues nicely into your next question.

I receive this error regarding integration, should I be worried about it? or how should I use quadrature to fix it?

For a well-designed formulation, over-integration shouldn’t lead to worse results, but it will do excessive computation during the assembly phase. You can manually set quadrature degree to some integer deg (i.e., the degree up to which quadrature is exact for polynomial integrands) by

dx = dx(metadata={"quadrature_degree":deg})

and likewise for ds or dS, if you are using those measures. (See the linked tutorial for other quadrature settings that are needed with quadrature function spaces.)

1 Like