# Iterative solver diverging

Hi, I am trying to build a displacement-controlled version of this tutorial Elasto-plastic analysis of a 2D von Mises material in which the author has used external force.
I tried to follow these posts Plasticity in fenics and Custom Newton Solver problem with Dirichlet conditions where @bleyerj states that to make a displacement controlled model, we need to keep 2 things in mind -

1. Add a DirichletBC with imposed value equal to the displacement increment over the time step and not the total displacement.
2. Only the first iteration should have a non-zero DirichletBC all other iterations must be computed with zero DirichletBC, otherwise you will accumulate imposed displacements over all iterations and will not converge.

Following both advices, I modified the loop like this -

``````Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
results = np.zeros((Nincr+1, 2))
A, Res = assemble_system(a, L, bc_u)
for bc in bc_u:
bc.apply(du.vector())
bc.homogenize()
nRes0 = Res.norm("l2")
nRes = nRes0
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = assemble_system(a, L, bc_u)
nRes = Res.norm("l2")
print("    Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
results[i+1, :] = (u(0.5, 1.)[1], t)
``````

But the solver diverges, I think I’ve not understood both the points due to which there is divergence. Therefore, I request to kindly give me suggestions on the same.
Thanks.

PS: Here is my MWE-

``````from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = "tsfc"
import warnings

# elastic parameters
E = Constant(70e3 )
nu = 0.3
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
Et = E/100.  # tangent modulus
H = E*Et/(E-Et)  # hardening modulus
sig0 = Constant(250.)

from mshr import *
domain= Rectangle(Point(0., 0.), Point(1., 1.))
mesh = generate_mesh(domain, 20)

deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
W = FunctionSpace(mesh, We)
W0 = FunctionSpace(mesh, W0e)

sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")

# Boundary conditions
class top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1., 1e-8)
class bot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., 1e-8)

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)

top().mark(boundaries, 1)
bot().mark(boundaries, 2)

bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 2)

bc_u = [bc_bottom, bc_top]

def eps(v):
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])

def sigma(eps_el):
return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el

def as_3D_tensor(X):
return as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])

ppos = lambda x: (x+abs(x))/2.
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
dp = ppos(f_elas)/(3*mu+H)
n_elas = s/sig_eq*ppos(f_elas)/f_elas
beta = 3*mu*dp/sig_eq
new_sig = sig_elas-beta*s
return as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp

def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)

a = inner(eps(v), sigma_tang(eps(u_)))*dxm
L = -inner(eps(u_), as_3D_tensor(sig))*dxm

def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return

Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
results = np.zeros((Nincr+1, 2))

A, Res = assemble_system(a, L, bc_u)
for bc in bc_u:
bc.apply(du.vector())
bc.homogenize()
nRes0 = Res.norm("l2")
nRes = nRes0
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = assemble_system(a, L, bc_u)
nRes = Res.norm("l2")
print("    Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
results[i+1, :] = (u(0.5, 1.)[1], t)
``````

Hi, so I tried to use this logic -

``````at first iteration -> nonzero DBC
second iteration onwards -> zero DBC -> homogenize them
``````

In the code, I made the following changes but still it diverges -

``````Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure
Nincr = 20
results = np.zeros((Nincr+1, 2))

A, Res = assemble_system(a, L, bc_u)
nRes0 = Res.norm("l2")
nRes = nRes0
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = assemble_system(a, L, bc_u)
if niter >= 1:
for bc in bc_u:
bc.apply(Du.vector())
bc.homogenize()

nRes = Res.norm("l2")
print("    Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
results[i+1, :] = (u(0.5, 1.)[1], t)
``````

I don’t know what mistake I’m making, any pointers about the same would be highly appreciable.
Thanks.

Okay, so I tried to make a simpler version of the problem to track down the issue and changed my boundary conditions to -

``````# Boundary conditions
class top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1., 1e-8)
class bot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., 1e-8)

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)

top().mark(boundaries, 1)
bot().mark(boundaries, 2)

bc_top = DirichletBC(V.sub(1), Constant(1.), boundaries, 1)
bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 2)

bc_u = [bc_bottom, bc_top]

``````

Modified the loop to -

``````Nitermax, tol = 200, 1e-6 # parameters of the Newton-Raphson procedure
Nincr = 20
results = np.zeros((Nincr+1, 2))

A, Res = assemble_system(a, L, bc_u)
err0 = Res.norm("l2")
err1 = err0
print(err1/ err0)
print("Increment:", str(i + 1))

## 1st iteration with nonzero DirichletBC conditions
Du.interpolate(Constant((0, 0.1)))
niter = 0
while err1/err0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "umfpack")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = assemble_system(a, L, bc_u)
err1 = Res.norm("l2")
print(err1/err0)

## 2nd iteration onwards => zero DirichletBC conditions
for bc in bc_u:
bc.homogenize()

niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
results[i+1, :] = (u(0.5, 1.)[1], t)
``````

Now the issue is I am getting NaN from 3rd iteration, I also tried looking up Default absolute tolerance and relative tolerance and specially what to do in nan values case, it says to change the inital guess to a non-zero value but I’m not quite sure how to do that, therefore I’m stuck.
Any help would be highly appreciable.
Thanks.