# Newton solver (manually computed) issue in a time loop

Dear Community,

I am trying to use the Newton Solver coded manually (I use the code posted in a previous post) but I get some trouble when I use it in a time loop.
I applied a pure shear on the top of a unit square.

When I test the code at one time step, the Newton solver works very well and I checked on paraview.

U_app.t = dt*20
NewtonSolver(residual, u, bcs, jacobian)
print("u.vector",u.vector().get_local() )

But when I test the code in a time loop I get some trouble. The first step computed is good and at the following steps, the displacement begins zeros :

U_app.t = t
print("t : ",t)
print("U_app.t : ",U_app.t, "t : ",t)
NewtonSolver(residual, u, bcs, jacobian)
print("u.vector",u.vector().get_local() )
vtkfileU      << u

For me, the Newton solver takes the previous solution to initialize the Jacobian and the residual at the next step. I probably miss something in the loop, but I don’t manage to find the mistake.

This the code :

import dolfin
from petsc4py import PETSc
import ufl
import numpy as np
import matplotlib.pyplot as plt
import os.path

dolfin.parameters["form_compiler"]["cpp_optimize"] = True
dolfin.parameters["form_compiler"]["representation"] = "uflacs"

# Material Parameters
Y, nu = 1e3, 0.3 # Young modulus and Poisson ratio
mu = dolfin.Constant(Y/(2*(1 + nu)))
lmbda = dolfin.Constant(Y*nu/((1 + nu)*(1 - 2*nu)))
K0 = dolfin.Constant(lmbda+2./3.*mu)

# Mesh & BC
Lx, Ly = 1.0, 1.0 # Dimensions of the beam
nx, ny = 2, 2 # Number of elements in each direction
mesh = dolfin.RectangleMesh(dolfin.Point(0,0),
dolfin.Point(Lx,Ly),
nx, ny)
left = dolfin.CompiledSubDomain("near(x[0],0) && on_boundary")
right = dolfin.CompiledSubDomain("near(x[0],Lx) && on_boundary", Lx=Lx)
bottom = dolfin.CompiledSubDomain("near(x[1],0) && on_boundary")
top = dolfin.CompiledSubDomain("near(x[1],Ly) && on_boundary", Ly=Ly)

boundary_markers = dolfin.MeshFunction("size_t", mesh, 1, 0)
boundary_indices = {"top": 1, "bottom": 2,  "right": 3, "left": 4}
top.mark(boundary_markers, boundary_indices["top"])
bottom.mark(boundary_markers, boundary_indices["bottom"])
right.mark(boundary_markers, boundary_indices["right"])
left.mark(boundary_markers, boundary_indices["left"])
ds = dolfin.ds(domain=mesh,subdomain_data=boundary_markers)
dx = dolfin.dx(domain=mesh)

# Space function
V_element = dolfin.VectorElement("CG", mesh.ufl_cell(), degree=2)
V = dolfin.FunctionSpace(mesh, V_element)
u = dolfin.Function(V, name="u")

# strain energy of a nonnlinear elastic material (deviatoric part).
plus_tr_eps = 0.5*( ufl.tr(eps)+abs( ufl.tr(eps)))
moins_tr_eps = 0.5*( ufl.tr(eps)-abs( ufl.tr(eps)))
sigma =  0.5* K0*(plus_tr_eps)*ufl.Identity(2) + mu * ufl.dev(eps) + 0.5* K0*(moins_tr_eps)*ufl.Identity(2)

tn = dolfin.Constant((0.0, 0.0))
residual = ufl.inner( sigma, eps_test ) * dx - ufl.dot( tn, v ) * ds # residual
jacobian = ufl.derivative(residual, u, dolfin.TrialFunction(V)) # jacobian

U_app = dolfin.Expression(("t","0.0"), t=0, degree=0)
bcs = [dolfin.DirichletBC( V , dolfin.Constant((0.,0.)), bottom),
dolfin.DirichletBC( V , U_app, top ),
dolfin.DirichletBC( V.sub(1) ,  dolfin.Constant(0.),right),
dolfin.DirichletBC( V.sub(1) ,  dolfin.Constant(0.), left)]

def NewtonSolver(residual, u, bcs, jacobian):
tol = 1e-8
error = 1.
it = 0

for bc in bcs:
bc.apply(u.vector())
bc.homogenize()

while error > tol and it < 8:
# Assemble
dg, g = dolfin.assemble_system(jacobian, -residual, bcs)
print('Jac_assemble :', dg.array())
print('Res_assemble :', g.get_local())

# PETSc solver
jacobian_p = dolfin.as_backend_type(dg).mat()
du, b = jacobian_p.getVecs()
b.setValues(range(g.size()), g)
ksp = PETSc.KSP().create()
ksp.setType('cg')
ksp.setOperators(jacobian_p)
ksp.solve(b, du)
#print('type(du) ;', type(du))
print('du.array ;', du.array)
error = g.norm('l2')
print('Iteration:', str(it), '; Error (g.norm(l2)): %3.3e' % error)
if error < tol:
break

# Update
u.vector()[:] = u.vector()[:] + du
it += 1

savedir        = "RESULT_test/"
vtkfileU       = dolfin.File(savedir+'displacement.pvd')

# Test at one time step t = load_max
u.interpolate(dolfin.Constant((0., 0.)))
# U_app.t = load_max # 0.3mm
# NewtonSolver(residual, u, bcs, jacobian)

U_app.t = t
print("t : ",t)
print("U_app.t : ",U_app.t)
NewtonSolver(residual, u, bcs, jacobian)
#newton_solver(u, monitor=simple_monitor )
print("u.vector",u.vector().get_local() )
vtkfileU      << u

Regards,

Lamia

Hello !

I managed to detect the error. It seems that at the second timestep the boundary conditions are not applied. I printed the u.vector() before entering in the iterative loop and after solve(A, du, b) with u = u +du.
But the code does something I don’t understand.
In the NewtonSolver function, before entering in the loop, there is :

for bc in bcs:
bc.apply(u.vector())
bc.homogenize() # otherwise it diverges

to initialize the displacement with BCs but it seems to work only at the first time step computed. But at the followings time step, any boundary conditions are applied.

This is the result printed when I start the time loop at

U_app.t = t
print("U_app.t : ",U_app.t)
NewtonSolver(residual, u, bcs, jacobian)
vtkfileU      << u

results :

U_app.t :  0.15384615384615383
u.vector before solve [0.15384615 0.         0.         0.         0.15384615 0.
0.         0.         0.15384615 0.         0.         0.
0.         0.         0.         0.         0.         0.        ]
du.array before solve :  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
du.array ; [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00  7.69230769e-02 -8.28330460e-17
0.00000000e+00  0.00000000e+00  7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00  7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00]
Iteration: 0 ; Error (g.norm(l2)): 4.831e+01
u.vector [ 1.53846154e-01  0.00000000e+00  0.00000000e+00  0.00000000e+00
1.53846154e-01  0.00000000e+00  7.69230769e-02 -8.28330460e-17
1.53846154e-01  0.00000000e+00  7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00  7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00]
du.array before solve :  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
du.array ; [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00  5.09813804e-17 -3.20120398e-20
0.00000000e+00  0.00000000e+00 -1.17392502e-16  0.00000000e+00
0.00000000e+00  0.00000000e+00  1.68442432e-16  0.00000000e+00
0.00000000e+00  0.00000000e+00]
Iteration: 1 ; Error (g.norm(l2)): 2.276e-13
U_app.t :  0.16153846153846152
u.vector before solve [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00  7.69230769e-02 -8.28330460e-17
0.00000000e+00  0.00000000e+00  7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00  7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00]
du.array before solve :  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
du.array ; [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00 -7.69230769e-02  8.10983225e-17
0.00000000e+00  0.00000000e+00 -7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00 -7.69230769e-02  0.00000000e+00
0.00000000e+00  0.00000000e+00]
Iteration: 0 ; Error (g.norm(l2)): 4.831e+01
u.vector [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00  1.38777878e-17 -1.73472348e-18
0.00000000e+00  0.00000000e+00 -6.93889390e-17  0.00000000e+00
0.00000000e+00  0.00000000e+00  8.32667268e-17  0.00000000e+00
0.00000000e+00  0.00000000e+00]
du.array before solve :  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
du.array ; [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
0.00000000e+00  0.00000000e+00 -1.38777878e-17  1.73472348e-18
0.00000000e+00  0.00000000e+00  6.93889390e-17  0.00000000e+00
0.00000000e+00  0.00000000e+00 -8.32667268e-17  0.00000000e+00
0.00000000e+00  0.00000000e+00]
Iteration: 1 ; Error (g.norm(l2)): 1.173e-13

Do you have an idea of what happen at the second timestep ?

Thanks for the help !

Lamia