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 :
for t in loads[20:25]:
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)
#loading
load_min, load_max, nsteps = 0.0, 0.3, 40
dt = load_max/nsteps
loads = np.linspace(load_min, load_max,nsteps)
# 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).
eps = ufl.sym( ufl.grad(u))
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))
v,eps_test = dolfin.TestFunction(V), ufl.sym( ufl.grad(v))
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)
for t in loads[20:25]:
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
Thanks for your help !
Regards,
Lamia