Newton Solver did not converge

mathematics
1

I am having problems with a nonlinear equation and getting it solved in the program. All the previous errors were relating to lines ~75-82, because that is where are weak form is defined. Are there suggestions to help the solver converge, or ways to format/type our code? This is also my first time posting in a help forum, please let me know if I did anything incorrectly. Anything would be appreciated, thank you!

from dolfin import *
#mesh = Mesh('mesh.xml')
#mesh = Mesh('untitled.msh')
from mshr import *
L = 1
P = 0.01
N = 150
domain = Rectangle(Point(0.0,0.0), Point(L,L)) - Rectangle(Point(0.25,0.5-(P/2)), Point(0.75, 0.5+(P/2)))
mesh = generate_mesh(domain,N)

import matplotlib.pyplot as plt
#plot(mesh)

V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
dp, p, q = TrialFunction(V), TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
Gc = 2.7
#Original
#l = 0.015
#Displacement boundary condition
E = 400e3
l = 0.01
nu = 0.2 #poisson ratio
lmbda = (nu * E) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
sf = 100 #maximum/fracture strength
psic = (sf**2)/(2*E)
m = (3*Gc)/(4*l*psic)

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))

def psi(u):
    return 0.5 * (lmbda + mu) * (0.5 * (tr(epsilon(u))+abs(tr(epsilon(u)))))**2 + mu * inner(dev(epsilon(u)),dev(epsilon(u)))

def H(uold,unew,Hold):
    return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)

def g(phi):
    return ((1 - phi)**2) / (((1 - phi)**2) + m * phi * (1 - phi))

def dg(phi):
    return (1 - phi)**2*(m*phi - m*(1 - phi) - 2*phi + 2)/(m*phi*(1 - phi) + (1 - phi)**2)**2 + (2*phi - 2)/(m*phi*(1 - phi) + (1 - phi)**2)

top = CompiledSubDomain("near(x[1], 1.0) && on_boundary")

bot = CompiledSubDomain("near(x[1], 0.0) && on_boundary")

def Crack(x):
    return 0.5-1e-3 <(x[1])< 0.5+1e-03 and x[0] <= 0.5

load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]
bc_phi = [DirichletBC(V, Constant(0.0), top)]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)

n = FacetNormal(mesh)
unew, uold = Function(W), Function(W)

pnew, pold, Hold = Function(V), Function(V), Function(V)

E_du = (g(pold))*inner(grad(v),sigma(u))*dx
E_phi = (0.5*(3.0/4.0)*Gc*l*inner(grad(p),grad(q)) + ((3*Gc)/(4*l))*q + H(uold,unew,Hold)*inner(q,dg(p)))*dx
F_phi = action(E_phi, pnew)
J_phi = derivative(F_phi, pnew, p)
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = NonlinearVariationalProblem(F_phi, pnew, bc_phi, J_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = NonlinearVariationalSolver(p_phi)
t = 0

u_r = 0.007
deltaT = 0.1
tol = 1e-5
conc_f = File ("./QuasiBritCohResults/phi.pvd")
fname = open('QBCForcevsDisp.txt', 'w')
while t<=2.0:
    t += deltaT
    if t >=0.7:
        deltaT = 0.01
    load.t=t*u_r
    iter = 0
    err = 1
    while err > tol:
        iter += 1
        solver_disp.solve()
        solver_phi.solve()
        err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
        err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
        err = max(err_u,err_phi)
        uold.assign(unew)
        pold.assign(pnew)
        Hold.assign(project(psi(unew), WW))
        if err < tol:
            print ('Iterations:', iter, ', Total time', t)
            if round(t*1e4) % 10 == 0:
                conc_f << pnew
                Traction = dot(sigma(unew),n)
                fy = Traction[1]*ds(1)
                fname.write(str(t*u_r) + "\t")
                fname.write(str(assemble(fy)) + "\n")

fname.close()
print ('Simulation completed')

*** 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.2.0.dev0
*** Git changeset: unknown

2

After defining the nonlinear solver object solver_phi, you can configure its parameters to make it more robust. The following seems to work for the current problem:

solver_phi = NonlinearVariationalSolver(p_phi)

# Access parameters object for solver.
param = solver_phi.parameters

# Use PETSc's implementation of a black-box nonlinear
# algebraic solver.  ("newton" uses the default
# implementation of Newton's method.)
stype = "snes" #"newton"
param["nonlinear_solver"] = stype
sparam = param[stype+"_solver"]

# Set the line search algorithm to back-tracking;
# the default, "basic", is just standard Newton
# iteration.
sparam["line_search"] = "bt" #"basic"

# This is sufficient for most practical problems.
sparam["relative_tolerance"] = 1e-3

To print out all the available options for a parameter set param, you can run info(param,True).