# Newton method implemented manually does not converge

**URL:** <https://fenicsproject.discourse.group/t/newton-method-implemented-manually-does-not-converge/4141>\
**Category:** Uncategorized\
**Created:** [September 22, 2020, 8:46pm UTC](https://fenicsproject.discourse.group/t/newton-method-implemented-manually-does-not-converge/4141 "2020-09-22T20:46:00Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![sabyasachi](https://avatars.discourse-cdn.com/v4/letter/s/58956e/32.png) [@sabyasachi](https://fenicsproject.discourse.group/u/sabyasachi)\
**Post date:** [September 22, 2020, 8:46pm UTC](https://fenicsproject.discourse.group/t/newton-method-implemented-manually-does-not-converge/4141/1 "2020-09-22T20:46:00Z")

</div>

Hello,

I have written a code to solve a large deformation elasticity problem. I apply DirichletBC on the top while the bottom is fixed (along the Y direction). I am using Newton method manually to solve this, but it is not converging. Here is a minimal code:

```auto
from __future__ import print_function
from fenics import *
from ufl import *
import numpy as np 
from dolfin import *
print(DirichletBC. __doc__ )

# Parameters
size = 1.0; 
mu = 41.0
lamb = 87.2
T_ = 1.0 # final time
num_steps = 2 # number of time steps
dt = T_ / num_steps # time step size

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(size, size, size), 5,5,5) # refinement, refinement, refinement)

V = VectorFunctionSpace(mesh, 'P', 1)

#Boundary conditions
v_top= 1e-3 

u_top = Expression('v_top*t',
                 degree=1, v_top=v_top, t=0.0)

def boundary_bot(x, on_boundary):
    return on_boundary and near(x[1],0.0,tol)

def boundary_top(x, on_boundary):
    return on_boundary and near(x[1],size,tol)

bc_bottom = DirichletBC(V.sub(1), Constant(0.0), boundary_bot)

bc_top = DirichletBC(V.sub(1), u_top, boundary_top)

bcs=[bc_bottom,bc_top] 

# Define variational problem
u = Function(V)
u_ = TestFunction(V)
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
d = u.geometric_dimension() # space dimension
v = TrialFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 0, 0))
I = Identity(d) # Identity tensor

def C_int(d):
    I = np.eye(d)
    II = np.einsum("ik,jl->ijkl", I, I) + np.einsum("il,jk->ijkl", I, I)
    JJ = np.einsum("ij,kl->ijkl", I, I)
    return as_tensor(mu*II + lamb*JJ)

def get_F(u):
    F= I + grad(u)
    return F

def get_S(F):
    E = 0.5*(F.T*F - I) # Right Cauchy-Green tensor
    i,j,k,l = indices(4)
    S= as_tensor(C_int(d)[i,j,k,l]*E[k,l],(i,j))
    return S

def get_P(u):
    F=get_F(u)
    E = 0.5*(F.T*F - I)   
    S= get_S(F)
    P = F*S 
    return P 

# residual and Jacobian
res= inner(nabla_grad(u_),get_P(Du))*dx - dot(T, u_)*ds  

Jac= derivative(res,Du,v)

# Time-stepping
t = 0.0
tol=1e-8
Nitermax=10

for n in range(num_steps):
    t += dt
    print('t=',round(t,3),'T=',T_)
    u_top.t = t 
    
    Du.interpolate(Constant((0, 0, 0)))
    niter = 0

    while niter<Nitermax:
        A, Res = assemble_system(Jac, -res, bcs)
        solve(A, du.vector(), Res) 
        Du.assign(Du+du)

        a= assemble(res)
        for bc in bcs: 
            bc.apply(a)

        nRes = a.norm("l2")
        print(" Residual:", nRes)
        niter += 1
    
    u.assign(u+Du)

```

Ideally, the norm of the residual should decrease quickly inside the while loop. However, it is almost not changing at all. Here is the output:

```auto
t= 0.5 T= 1.0
    Residual: 0.003000596335349779
    Residual: 0.0030000046508762163
    Residual: 0.0030000046197302807
    Residual: 0.003000004621295162
    Residual: 0.003000004622873061
    Residual: 0.003000004624452816
    Residual: 0.003000004626034452
    Residual: 0.0030000046276179465
    Residual: 0.003000004629203315
    Residual: 0.0030000046307905636
t= 1.0 T= 1.0
    Residual: 0.0060047676568817715
    Residual: 0.006000037496833246
    Residual: 0.006000036995869351
    Residual: 0.0060000370208320025
    Residual: 0.0060000370462124746
    Residual: 0.006000037071653371
    Residual: 0.006000037097154581
    Residual: 0.006000037122716251
    Residual: 0.006000037148338622
    Residual: 0.00600003717402185

```

Can anyone please let me know what is wrong here.  
Thanks.

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [September 22, 2020, 9:01pm UTC](https://fenicsproject.discourse.group/t/newton-method-implemented-manually-does-not-converge/4141/2 "2020-09-22T21:01:35Z")

</div>

The problem comes from applying the Dirichlet bcs for each iteration corrections, [see also here](https://fenicsproject.discourse.group/t/custom-newton-solver-problem-with-dirichlet-conditions/2329)

---

<div class="post-metadata">

**Author:** ![sabyasachi](https://avatars.discourse-cdn.com/v4/letter/s/58956e/32.png) [@sabyasachi](https://fenicsproject.discourse.group/u/sabyasachi)\
**Post date:** [September 22, 2020, 9:39pm UTC](https://fenicsproject.discourse.group/t/newton-method-implemented-manually-does-not-converge/4141/3 "2020-09-22T21:39:39Z")

</div>

Thanks a lot, that indeed solves the problem.