# Newton solver converge initially, then diverge and blow up

Hello everyone,
I am trying to modify yield condition by tweaking @bleyerj 's tutorial Elasto-plastic analysis of a 2D von Mises material.
Given von-mises yield condition :

I want to change it to :

where :-

I have tried to do so in the following way (MWE) :-

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

domain = mshr.Rectangle(Point(0,0), Point(1.,1.))
mesh = mshr.generate_mesh(domain, 30)
plot(mesh);

E = 10e3
nu = 0.33
sig0C = 60.
sig0T = 45.
m = sig0C/sig0T

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
Et = E/100.
H = E*Et/(E-Et)

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")

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)
sk = sig_elas - s
i111 = tr(sk)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = (1./2.*m)*((m - 1)*i111+ (m+1)*sig_eq) - sig0T - 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)

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

class bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., 1e-8)

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

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
bottom().mark(boundaries, 1)
top().mark(boundaries, 2)

bc_bottom = DirichletBC(V.sub(1), Constant((0.)), boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
bc_u = [bc_bottom]

def F_ext(v):

a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)

Nitermax, tol = 200, 1e-8
Nincr = 1000
results = np.zeros((Nincr+1, 2))
A, Res = assemble_system(a_Newton, res, bc_u)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > 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_Newton, res, 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))

k = p_avg.vector().max()

if k >= 0.02:
break
``````

Initially, it converges nicely but as soon as it reaches the yield condition (i.e. plastic strain will become non-zero in next iteration), it diverges.

and solution to my problem according to above might be to try easing relaxation parameter.
Therefore, I also tried custom newton solver for the same but ran into Arity mismatch :-
What I tried :-

``````# Define a Nonlinear problem to assemble the residual and Jacobian
class Problem(NonlinearProblem):
def __init__(self, a, L):
self.a = a
self.L = L
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.L, tensor=b)

def J(self, A, x):
assemble(self.a, tensor=A)

# Use a Newton solver with custom damping parameter
class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())

def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)

PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps")

self.linear_solver().set_from_options()

def update_solution(self, x, dx, relaxation_parameter, nonlinear_problem, iteration):
tau = 1.0
theta = min(sqrt(2.0*tau/norm(dx, norm_type="l2", mesh=V.mesh())), 1.0)
info("Newton damping parameter: %.3e" % theta)
x.axpy(-theta, dx)

F = inner(eps(v), sigma_tang(eps(u_)))*dxm + inner(eps(u_), as_3D_tensor(sig))*dxm - F_ext(u_)
J = derivative(F, u)
problem = Problem(J, F)
solver = CustomSolver()
solver.solve(problem, u.vector())
``````

Error :-

``````---------------------------------------------------------------------------
ArityMismatch                             Traceback (most recent call last)
/tmp/ipykernel_720/716830184.py in <module>
26 problem = Problem(J, F)
27 solver = CustomSolver()
---> 28 solver.solve(problem, u.vector())

/tmp/ipykernel_720/2688388353.py in F(self, b, x)
7
8     def F(self, b, x):
----> 9         assemble(self.L, tensor=b)
10
11     def J(self, A, x):

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
200         dolfin_form = form
201     else:
--> 202         dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
203
204     # Create tensor

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
58         return form
59     elif isinstance(form, ufl.Form):
---> 60         return Form(form,
61                     form_compiler_parameters=form_compiler_parameters,
62                     function_spaces=function_spaces)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py in __init__(self, form, **kwargs)
41             form_compiler_parameters = {"external_include_dirs": dolfin_pc["include_dirs"]}
42
---> 43         ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
44                            mpi_comm=mesh.mpi_comm())
45         ufc_form = cpp.fem.make_ufc_form(ufc_form[0])

/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in mpi_jit(*args, **kwargs)
48         # Just call JIT compiler when running in serial
49         if MPI.size(mpi_comm) == 1:
---> 50             return local_jit(*args, **kwargs)
51
52         # Default status (0 == ok, 1 == fail)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in ffc_jit(ufl_form, form_compiler_parameters)
98     p.update(dict(parameters["form_compiler"]))
99     p.update(form_compiler_parameters or {})
--> 100     return ffc.jit(ufl_form, parameters=p)
101
102

/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit(ufl_object, parameters, indirect)
215
216     # Inspect cache and generate+build if necessary
--> 217     module = jit_build(ufl_object, module_name, parameters)
218
219     # Raise exception on failure to build or import module

/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit_build(ufl_object, module_name, parameters)
128
129     # Carry out jit compilation, calling jit_generate only if needed
--> 130     module, signature = dijitso.jit(jitable=ufl_object,
131                                     name=module_name,
132                                     params=params,

/usr/lib/python3/dist-packages/dijitso/jit.py in jit(jitable, name, params, generate, send, receive, wait)
163         elif generate:
164             # 1) Generate source code
--> 165             header, source, dependencies = generate(jitable, name, signature, params["generator"])
166             # Ensure we got unicode from generate

/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit_generate(ufl_object, module_name, signature, parameters)
63         compile_object = compile_coordinate_mapping
64
---> 65     code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
66             prefix=module_name, parameters=parameters, jit=True)
67

/usr/lib/python3/dist-packages/ffc/compiler.py in compile_form(forms, object_names, prefix, parameters, jit)
140                  prefix="Form", parameters=None, jit=False):
141     """This function generates UFC code for a given UFL form or list of UFL forms."""
--> 142     return compile_ufl_objects(forms, "form", object_names,
143                                prefix, parameters, jit)
144

/usr/lib/python3/dist-packages/ffc/compiler.py in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
183     # Stage 1: analysis
184     cpu_time = time()
--> 185     analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
186     _print_timing(1, time() - cpu_time)
187

/usr/lib/python3/dist-packages/ffc/analysis.py in analyze_ufl_objects(ufl_objects, kind, parameters)
87
88         # Analyze forms
---> 89         form_datas = tuple(_analyze_form(form, parameters)
90                            for form in forms)
91

/usr/lib/python3/dist-packages/ffc/analysis.py in <genexpr>(.0)
87
88         # Analyze forms
---> 89         form_datas = tuple(_analyze_form(form, parameters)
90                            for form in forms)
91

/usr/lib/python3/dist-packages/ffc/analysis.py in _analyze_form(form, parameters)
181             error("Failed to import tsfc.ufl_utils.compute_form_data when asked "
182                   "for tsfc representation.")
--> 183         form_data = tsfc_compute_form_data(form)

~/.local/lib/python3.10/site-packages/tsfc/ufl_utils.py in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_restrictions, do_estimate_degrees)
47     other form compilers based on TSFC.
48     """
---> 49     fd = ufl_compute_form_data(
50         form,
51         do_apply_function_pullbacks=do_apply_function_pullbacks,

/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
405     preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
406
--> 407     check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
408
409     # TODO: This member is used by unit tests, change the tests to

/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py in check_form_arity(form, arguments, complex_mode)
175 def check_form_arity(form, arguments, complex_mode=False):
176     for itg in form.integrals():
--> 177         check_integrand_arity(itg.integrand(), arguments, complex_mode)

/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py in check_integrand_arity(expr, arguments, complex_mode)
157                              key=lambda x: (x.number(), x.part())))
158     rules = ArityChecker(arguments)
--> 159     arg_tuples = map_expr_dag(rules, expr, compress=False)
160     args = tuple(a[0] for a in arg_tuples)
161     if args != arguments:

/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dag(function, expression, compress, vcache, rcache)
34     Return the result of the final function call.
35     """
---> 36     result, = map_expr_dags(function, [expression], compress=compress,
37                             vcache=vcache,
38                             rcache=rcache)

/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dags(function, expressions, compress, vcache, rcache)
97                 r = handlers[v._ufl_typecode_](v)
98             else:
---> 99                 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
100
101             # Optionally check if r is in rcache, a memory optimization

/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py in sum(self, o, a, b)
46     def sum(self, o, a, b):
47         if a != b:
---> 48             raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
49         return a
50

ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('v_0', 'v_1').
``````

Kindly suggest some pointers on how to go about the same. Any thoughts are most welcome.
Thanks.

ArityMismatch: Adding expressions with non-matching form arguments (‘v_0’,) vs (‘v_0’, ‘v_1’).

That typically means that you have the sum of two addends in which one contains only a test function and the other both a test function and a trial function. Check that `self.L` does not contain any trial function: it shouldn’t, as you would need to use the solution itself in place of any trial function.

Thanks for the reply and solution. I was using TrialFunction in F. Changing it to solution solves the problem. I have one last question. If I want to solve this nonlinear problem iteratively, how can I do it like as it is given in tutorial above mentioned (Please refer to above post) ?
For example - Tutorial’s newton-raphson loop is as follows -

``````Nitermax, tol = 200, 1e-8
Nincr = 1000
results = np.zeros((Nincr+1, 2))
A, Res = assemble_system(a_Newton, res, bc_u)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > 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_Newton, res, 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))
``````

My attempt -

``````for t in range(5):
print(t)

u.interpolate(Constant((0.1, 0.1)))
Du.interpolate(Constant((0, 0)))
problem = Problem(J, F)
solver = CustomSolver()
solver.solve(problem, u.vector())

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)

u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
t += 1
``````

But this is leading to -

``````0
Newton iteration 0: r (abs) = 1.674e-09 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton damping parameter: 1.000e+00
Newton iteration 1: r (abs) = 4.963e-11 (tol = 1.000e-10) r (rel) = 2.966e-02 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
1
Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_1605/1390456215.py in <module>
7     problem = Problem(J, F)
8     solver = CustomSolver()
----> 9     solver.solve(problem, u.vector())
10
11     Du.assign(Du+du)

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
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------
``````

Please do let me know if I need to make another topic for this.
Thanks for the help.

It’s difficult for me to comment, because I have no experience in elasto-plasticity.

However, NaN typically happens in case of meaningless floating point operations, like having a zero at the denominator of a fraction or (when using real numbers, as legacy FEniCS does) trying to take the square root of a negative number.

From your snippet, I would look at

``````    u.interpolate(Constant((0.1, 0.1)))
Du.interpolate(Constant((0, 0)))
``````

: do you really need to reset them at every time, or can you assign these values at time equal to zero and for later times use the converged ones at the previous time? I do understand that Jeremy is indeed resetting `Du` at every time, but I don’t see that being done for `u`.

Thanks for the suggestion. I tried it but 2nd iteration onwards when I interpolate u with previous solution. It produces nan values. I saw post on the same like PETSc Krylov Solver RuntimeError with fine mesh where nan values mean that my initial guess is solution.
Here is MWE with your suggestion -

``````from dolfin import *
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
import mshr

domain = mshr.Rectangle(Point(0,0), Point(1.,1.))
mesh = mshr.generate_mesh(domain, 30)

E = 10e3
nu = 0.33
sig0C = 60.
sig0T = 45.
m = sig0C/sig0T

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
Et = E/100.
H = E*Et/(E-Et)

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")
w = Function(V, name="Previous 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")

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)
sk = sig_elas - s
i111 = tr(sk)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = (1/(2.*m))*((m - 1)*i111+ (m+1)*sig_eq) - sig0T - 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)

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

class Problem(NonlinearProblem):
def __init__(self, a, L):
self.a = a
self.L = L
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.L, tensor=b)

def J(self, A, x):
assemble(self.a, tensor=A)

class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())

def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)

PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "umfpack")

self.linear_solver().set_from_options()

def update_solution(self, x, dx, relaxation_parameter, nonlinear_problem, iteration):
tau = 1.0
theta = min(sqrt(2.0*tau/norm(dx, norm_type="l2", mesh=V.mesh())), 1.0)
info("Newton damping parameter: %.3e" % theta)
x.axpy(-theta, dx)

class bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., 1e-8)

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

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
bottom().mark(boundaries, 1)
top().mark(boundaries, 2)

bc_bottom = DirichletBC(V.sub(1), Constant((0.)), boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
bc_u = [bc_bottom]

def F_ext(v):

F = inner(eps(u), sigma_tang(eps(u_)))*dxm + inner(eps(u_), as_3D_tensor(sig))*dxm - F_ext(u_)
J = derivative(F, u)

u.interpolate(Expression(('x[0]/1000.', 'x[1]/1000.'), degree = 1))

for t in range(5):

if t > 0:
u.interpolate(w)

Du.interpolate(Constant((0, 0)))
problem = Problem(J, F)
solver = CustomSolver()
solver.solve(problem, u.vector())

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)

w.assign(u)
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
t += 1
``````

Error -

``````0
Newton iteration 0: r (abs) = 5.794e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton damping parameter: 1.000e+00
Newton iteration 1: r (abs) = 7.321e-13 (tol = 1.000e-10) r (rel) = 1.263e-13 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
1
Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 2: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 3: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 4: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 5: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 6: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 7: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 8: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 9: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 10: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 11: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 12: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 13: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 14: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 15: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 16: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 17: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 18: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 19: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 20: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 21: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 22: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 23: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 24: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 25: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 26: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 27: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 28: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 29: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 30: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 31: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 32: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 33: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 34: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 35: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 36: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 37: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 38: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 39: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 40: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 41: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 42: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 43: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 44: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 45: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 46: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 47: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 48: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 49: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 50: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_2248/761677441.py in <module>
9     problem = Problem(J, F)
10     solver = CustomSolver()
---> 11     solver.solve(problem, u.vector())
12
13     Du.assign(Du+du)

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

For now, I would like to read How to approximate the error without exact solution? and see what can be done best from my side.
Thanks.

You mean, `w` did not have NaNs, but `u` after interpolation does? You can easily verify that by printing `u.vector().get_local()`, and similarly for `w`.

BTW, how did you install FEniCS to run this code?
I am unable to reproduce locally because I don’t have around a FEniCS installation with `tsfc`, and the `firedrake` installation I have surely ships a `tsfc` installation which is too new (for instance, it tries to import `ufl` instead of `ufl_legacy`).

This is the output when I’m printing values of `u` and `w` (decreased mesh density for sake of convenience) -

``````t:  0
u:  [0.001      0.         0.001      0.00025    0.00075    0.
0.000875   0.000125   0.000875   0.         0.001      0.000125
0.00070833 0.00029167 0.0005     0.         0.00072917 0.00014583
0.00060417 0.00014583 0.000625   0.         0.00085417 0.00027083
0.001      0.0005     0.00085417 0.00039583 0.001      0.000375
0.00075    0.000625   0.0005     0.0005     0.000375   0.00025
0.0004375  0.000375   0.00054167 0.00027083 0.0004375  0.000125
0.00025    0.         0.0003125  0.000125   0.000375   0.
0.00072917 0.00045833 0.000875   0.0005625  0.001      0.00075
0.000875   0.0006875  0.001      0.000625   0.00060417 0.00039583
0.000625   0.0005625  0.00053472 0.00074306 0.00029167 0.00070833
0.00041319 0.00072569 0.00051736 0.00062153 0.00039583 0.00060417
0.00024641 0.00046336 0.0003107  0.00035668 0.0003732  0.00048168
0.00026904 0.00058585 0.00064236 0.00068403 0.00075    0.001
0.00075    0.0008125  0.000875   0.000875   0.001      0.001
0.000875   0.001      0.001      0.000875   0.00064236 0.00087153
0.         0.00025    0.0001232  0.00035668 0.0001875  0.00025
0.000125   0.000125   0.         0.         0.000125   0.
0.         0.000125   0.0005     0.001      0.00051736 0.00087153
0.00039583 0.00085417 0.00025    0.001      0.         0.00075
0.         0.0005     0.0001232  0.00048168 0.         0.000375
0.000625   0.001      0.00014583 0.00072917 0.00014583 0.00060417
0.         0.000625   0.00027083 0.00085417 0.000125   0.000875
0.000375   0.001      0.         0.001      0.         0.000875
0.000125   0.001     ]
w:  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0.]
Newton iteration 0: r (abs) = 1.628e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
t:  1
u:  [0.00122984 0.00082161 0.00127491 0.00082161 0.00122984 0.00086668
0.00125238 0.00084415 0.00122984 0.00084415 0.00125238 0.00082161
0.00128243 0.0008742  0.00122984 0.00091176 0.00125613 0.00087044
0.00125613 0.00089298 0.00122984 0.00088922 0.00127867 0.0008479
0.00131999 0.00082161 0.00130121 0.0008479  0.00129745 0.00082161
0.00134252 0.00086668 0.00131999 0.00091176 0.00127491 0.0009343
0.00129745 0.00092303 0.00127867 0.00090425 0.00125238 0.00092303
0.00122984 0.00095683 0.00125238 0.00094556 0.00122984 0.0009343
0.00131247 0.00087044 0.00133126 0.00084415 0.00136506 0.00082161
0.00135379 0.00084415 0.00134252 0.00082161 0.00130121 0.00089298
0.00133126 0.00088922 0.00136381 0.0009055  0.00135755 0.00094932
0.00136068 0.00092741 0.0013419  0.00090863 0.00133877 0.00093054
0.00131338 0.00095748 0.00129415 0.00094589 0.00131668 0.00093462
0.00133546 0.0009534  0.00135317 0.00088609 0.00141013 0.00086668
0.00137633 0.00086668 0.0013876  0.00084415 0.00141013 0.00082161
0.00141013 0.00084415 0.0013876  0.00082161 0.00138697 0.00088609
0.00127491 0.00100191 0.00129415 0.00097969 0.00127491 0.0009681
0.00125238 0.00097937 0.00122984 0.00100191 0.00122984 0.00097937
0.00125238 0.00100191 0.00141013 0.00091176 0.00138697 0.00090863
0.00138384 0.00093054 0.00141013 0.00095683 0.00136506 0.00100191
0.00131999 0.00100191 0.00131668 0.00097969 0.00129745 0.00100191
0.00141013 0.00088922 0.0013613  0.00097561 0.00133877 0.00097561
0.00134252 0.00100191 0.00138384 0.00095308 0.0013876  0.00097937
0.00141013 0.0009343  0.00141013 0.00100191 0.0013876  0.00100191
0.00141013 0.00097937]
w:  [0.00122984 0.00082161 0.00127491 0.00082161 0.00122984 0.00086668
0.00125238 0.00084415 0.00122984 0.00084415 0.00125238 0.00082161
0.00128243 0.0008742  0.00122984 0.00091176 0.00125613 0.00087044
0.00125613 0.00089298 0.00122984 0.00088922 0.00127867 0.0008479
0.00131999 0.00082161 0.00130121 0.0008479  0.00129745 0.00082161
0.00134252 0.00086668 0.00131999 0.00091176 0.00127491 0.0009343
0.00129745 0.00092303 0.00127867 0.00090425 0.00125238 0.00092303
0.00122984 0.00095683 0.00125238 0.00094556 0.00122984 0.0009343
0.00131247 0.00087044 0.00133126 0.00084415 0.00136506 0.00082161
0.00135379 0.00084415 0.00134252 0.00082161 0.00130121 0.00089298
0.00133126 0.00088922 0.00136381 0.0009055  0.00135755 0.00094932
0.00136068 0.00092741 0.0013419  0.00090863 0.00133877 0.00093054
0.00131338 0.00095748 0.00129415 0.00094589 0.00131668 0.00093462
0.00133546 0.0009534  0.00135317 0.00088609 0.00141013 0.00086668
0.00137633 0.00086668 0.0013876  0.00084415 0.00141013 0.00082161
0.00141013 0.00084415 0.0013876  0.00082161 0.00138697 0.00088609
0.00127491 0.00100191 0.00129415 0.00097969 0.00127491 0.0009681
0.00125238 0.00097937 0.00122984 0.00100191 0.00122984 0.00097937
0.00125238 0.00100191 0.00141013 0.00091176 0.00138697 0.00090863
0.00138384 0.00093054 0.00141013 0.00095683 0.00136506 0.00100191
0.00131999 0.00100191 0.00131668 0.00097969 0.00129745 0.00100191
0.00141013 0.00088922 0.0013613  0.00097561 0.00133877 0.00097561
0.00134252 0.00100191 0.00138384 0.00095308 0.0013876  0.00097937
0.00141013 0.0009343  0.00141013 0.00100191 0.0013876  0.00100191
0.00141013 0.00097937]
Newton damping parameter: 1.000e+00
Newton iteration 1: r (abs) = 7.765e-14 (tol = 1.000e-10) r (rel) = 4.771e-15 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 2: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 3: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 4: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 5: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 6: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 7: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 8: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 9: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 10: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 11: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 12: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 13: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 14: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 15: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 16: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 17: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 18: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 19: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 20: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 21: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 22: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 23: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 24: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 25: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 26: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 27: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 28: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 29: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 30: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 31: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 32: r (abs) = -nan (tol = 1.000e-10) r (rel) =
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_2338/2221793411.py in <module>
13     problem = Problem(J, F)
14     solver = CustomSolver()
---> 15     solver.solve(problem, u.vector())
16
17     Du.assign(Du+du)

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

-nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 33: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 34: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 35: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 36: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 37: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 38: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 39: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 40: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 41: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 42: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 43: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 44: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 45: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 46: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 47: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 48: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 49: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
Newton iteration 50: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
``````

I installed WSL then on Ubuntu terminal executed -

``````sudo apt-get install software-properties-common
sudo apt-get update
sudo apt-get install fenics
``````

Later I added TSFC by following -

I won’t be able to install `tsfc`, because downgrading to `2018.1.0` is not really an option I can afford in my current setup.

Still, I have a question just looking at your output

``````t:  1
u:  [...]
w:  [...]
Newton damping parameter: 1.000e+00
Newton iteration 1: r (abs) = 7.765e-14 (tol = 1.000e-10) r (rel) = 4.771e-15 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
# ----------> WHY DOES IT START BACK FROM ITERATION 0 WITHOUT PRINTING t? <---------- #
Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton damping parameter: nan
``````

Thanks for pointing that out. I didn’t notice that. I don’t know why it is happening ? I’m using below loop.

``````for t in range(5):
print('t: ', t)
print('u: ', u.vector().get_local())
print('w: ', w.vector().get_local())

if t > 0:
u.interpolate(w)

Du.interpolate(Constant((0, 0)))
problem = Problem(J, F)
solver = CustomSolver()
solver.solve(problem, u.vector())

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)

w.assign(u)
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
t += 1
``````

One question :- If I remove `F_ext(u_)` and try to build displacement-controlled model, can it solve the diverging issue ?

``````F = inner(eps(u), sigma_tang(eps(u_)))*dxm + inner(eps(u_), as_3D_tensor(sig))*dxm - F_ext(u_)
``````