# Newton solver converge initially, then diverge and blow up

**URL:** <https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315>\
**Category:** Legacy DOLFIN\
**Created:** [January 13, 2024, 11:56pm UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315 "2024-01-13T23:56:51Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![violetus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/violetus/32/5122_2.png) [@violetus](https://fenicsproject.discourse.group/u/violetus)\
**Post date:** [January 13, 2024, 11:56pm UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/1 "2024-01-13T23:56:51Z")

</div>

Hello everyone,  
I am trying to modify yield condition by tweaking @bleyerj 's tutorial [Elasto-plastic analysis of a 2D von Mises material](https://comet-fenics.readthedocs.io/en/latest/demo/2D_plasticity/vonMises_plasticity.py.html).  
Given von-mises yield condition :  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/b86ce64b8d360b8c769c2d58255168f6085b5f0a.png)

I want to change it to :  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/9de0bf57277d0155a9be6e794fe4391d0e96f784.png)

where :-

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/a/a0a7abd752bb95c6ca20227e326d45fc670b79b5.png)

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

```auto
from dolfin import *
parameters["form_compiler"]["representation"] = 'quadrature'
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
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)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
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):
    e = sym(grad(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)

metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

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)

loading = Expression("t", t=0, degree=1)
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):
    return loading*dot(n, v)*ds(2)

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
load_steps = np.linspace(0, 3000., Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):    
    loading.t = t
    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.  
I read @nate’s post on

> [@Default absolute tolerance and relative tolerance](https://fenicsproject.discourse.group/t/default-absolute-tolerance-and-relative-tolerance/3829/4):
>
> I write the following in the context of FE modelling rather than numerical analysis. I.e. this should be interpreted as tips and tricks. Perhaps someone else can contribute something more rigorous, or their own experience. Do you know what to expect the residual of your system to be to machine precision? Adjust your absolute tolerance accordingly. E.g: A Poisson problem with a diffusion coefficient of 1 may yield a residual to machine precision around 10^{-12}. Linear elasticity with a Young’…

and solution to my problem according to above might be to try easing relaxation parameter.  
Therefore, I also tried [custom newton solver](https://bitbucket.org/nate-sime/dolfin_dg/src/6eef8d491b5f299ce8b21c4f0aa4eb0231947ec4/demo/dolfin/compressible_navier_stokes_naca0012/dg_naca0012_2d.py#lines-51:55) for the same but ran into Arity mismatch :-  
What I tried :-

```auto
# 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 :-

```auto
---------------------------------------------------------------------------
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
    167 header = as_unicode(header)

/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)
    184 elif r == "quadrature":
    185 # quadrature representation

~/.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.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [January 14, 2024, 6:18am UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/2 "2024-01-14T06:18:36Z")

</div>

> 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.

---

<div class="post-metadata">

**Author:** ![violetus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/violetus/32/5122_2.png) [@violetus](https://fenicsproject.discourse.group/u/violetus)\
**Post date:** [January 14, 2024, 7:44am UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/3 "2024-01-14T07:44:58Z")

</div>

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 -

```auto
Nitermax, tol = 200, 1e-8  
Nincr = 1000
load_steps = np.linspace(0, 3000., Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):    
    loading.t = t
    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 -

```auto
for t in range(5):
    print(t)
    loading.t = 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 -

```auto
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
***
*** 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 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.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [January 14, 2024, 9:00am UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/4 "2024-01-14T09:00:55Z")

</div>

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

```auto
    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`.

---

<div class="post-metadata">

**Author:** ![violetus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/violetus/32/5122_2.png) [@violetus](https://fenicsproject.discourse.group/u/violetus)\
**Post date:** [January 14, 2024, 10:07am UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/5 "2024-01-14T10:07:38Z")

</div>

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](https://fenicsproject.discourse.group/t/petsc-krylov-solver-runtimeerror-with-fine-mesh/3418) where nan values mean that my initial guess is solution.  
Here is MWE with your suggestion -

```auto
from dolfin import *
parameters["form_compiler"]["representation"] = 'quadrature'
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
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)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
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):
    e = sym(grad(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)

metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)

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)

loading = Expression("t", t=0, degree=1)
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):
    return loading*dot(n, v)*ds(2)

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):
    loading.t = t
    
    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 -

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

```

For now, I would like to read [How to approximate the error without exact solution?](https://fenicsproject.discourse.group/t/how-to-approximate-the-error-without-exact-solution/2732) and see what can be done best from my side.  
Thanks.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [January 14, 2024, 10:58am UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/6 "2024-01-14T10:58:58Z")

</div>

> [@violetus](#):
>
> when I interpolate u with previous solution. It produces nan values

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`).

---

<div class="post-metadata">

**Author:** ![violetus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/violetus/32/5122_2.png) [@violetus](https://fenicsproject.discourse.group/u/violetus)\
**Post date:** [January 14, 2024, 12:28pm UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/7 "2024-01-14T12:28:59Z")

</div>

> [@francesco-ballarin](#):
>
> 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`.

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

```auto
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
***
*** 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.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)

```

> [@francesco-ballarin](#):
>
> BTW, how did you install FEniCS to run this code?

I installed WSL then on Ubuntu terminal executed -

```auto
sudo apt-get install software-properties-common
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt-get update
sudo apt-get install fenics

```

Later I added TSFC by following -

> [@Quadrature Representation](https://fenicsproject.discourse.group/t/quadrature-representation/2025/2):
>
> I’ve found the new form compiler [TSFC](https://doi.org/10.1137/17M1130642) to be more robust with quadrature elements. It’s not included in FEniCS by default, but it can be installed manually using pip3, pip3 install git+https://github.com/blechta/tsfc.git@2018.1.0 pip3 install git+https://github.com/blechta/COFFEE.git@2018.1.0 pip3 install git+https://github.com/blechta/FInAT.git@2018.1.0 pip3 install singledispatch networkx pulp then used in FEniCS Python scripts much like 'quadrature' representation: parameters["form\_compile…

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [January 15, 2024, 3:24pm UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/8 "2024-01-15T15:24:39Z")

</div>

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

```auto
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

```

---

<div class="post-metadata">

**Author:** ![violetus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/violetus/32/5122_2.png) [@violetus](https://fenicsproject.discourse.group/u/violetus)\
**Post date:** [January 15, 2024, 4:50pm UTC](https://fenicsproject.discourse.group/t/newton-solver-converge-initially-then-diverge-and-blow-up/13315/9 "2024-01-15T16:50:45Z")

</div>

> [@francesco-ballarin](#):
>
> ```auto
> # ----------> WHY DOES IT START BACK FROM ITERATION 0 WITHOUT PRINTING t? <---------- #
> 
> ```

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

```auto
for t in range(5):
    print('t: ', t)
    print('u: ', u.vector().get_local())
    print('w: ', w.vector().get_local())
    loading.t = t
    
    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 ?

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

```