Hello everyone,
For some further applications, I need to access to the assembled system at each Newton iteration of a nonlinear hyperelasticity problem. Therefore I build my own Newton Solver using petsc4py:
def NewtonSolver(G, u, bcs, dG):
tol = 1e-10
error = 1.
it = 0
while error > tol and it < 10:
# Assemble system
dg, g = d.assemble_system(dG, G, bcs)
# Manual BC
#g = d.assemble(G)
#dg = d.assemble(dG)
# Apply BC
#for bc in bcs:
#bc.apply(g)
#bc.apply(dg)
error = g.norm('l2')
print('Iteration:', str(it), '; Error: %3.3e' % error)
if error < tol:
break
# PETSc solver
dGp = d.as_backend_type(dg).mat()
x, b = dGp.getVecs()
b.setValues(range(g.size()), g)
ksp = PETSc.KSP().create()
ksp.setType('cg')
ksp.setOperators(dGp)
ksp.solve(b, x)
# Update
u.vector()[:] = u.vector()[:] - x
it += 1
And I used it to solve a simple example of a square with zero displacement at the left boundary and a traction or a prescribed displacement at the right boundary (the following code is for the second example, but it can be easily changed to the first one):
#%% Mesh
E = 250.
nu = 0.3
mu = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
mesh = d.UnitSquareMesh(3,3)
class C2(d.SubDomain):
def inside(self,x,on_boundary):
return d.near(x[0],1)
class C4(d.SubDomain):
def inside(self,x,on_boundary):
return d.near(x[0],0)
c2 = C2()
c4 = C4()
#%% Function Spaces and Functions
V = d.VectorFunctionSpace(mesh, 'CG', 1)
u = d.Function(V)
du = d.TrialFunction(V)
delta_u = d.TestFunction(V)
#%% Boundary Conditions
boundary_markers = d.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
c2.mark(boundary_markers, 2)
c4.mark(boundary_markers, 4)
d.File('bc.pvd') << boundary_markers
u0 = d.Constant((0., 0.))
u1 = d.Constant((0.1, 0.))
bc4 = d.DirichletBC(V, u0, boundary_markers, 4)
bc2 = d.DirichletBC(V, u1, boundary_markers, 2)
ds = d.Measure("ds", domain = V.mesh(), subdomain_data=boundary_markers)
bcs = [bc2, bc4]
#%% NL problem
# Volume
dim = len(u)
I = d.Identity(dim) # Identity tensor
F = I + d.grad(u) # Deformation gradient
J = d.det(F)
FiT = d.inv(F.T)
PK1 = mu*F + (lmbda*d.ln(J) - mu)*FiT
gradu = d.grad(delta_u)
T = d.Constant((30.,0.))
G = d.inner(PK1, gradu)*d.dx #- d.dot(T, delta_u)*ds(2)
dG = d.derivative(G, u)
I am comparing the results with the NonlinearVariationalSolver:
#%% Solve
u.rename('u','u')
# Fenics solver
problem = d.NonlinearVariationalProblem(G, u, bcs, dG)
solver = d.NonlinearVariationalSolver(problem)
solver.solve()
d.File('u0.pvd') << u
# Custom solver
NewtonSolver(G, u, bcs, dG)
d.File('u1.pvd') << u
The problem with traction at the right end works very well, although the residual norm for the last steps is not exactly the same as in the Fenics solver:
NonlinearVariationalSolver:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.581e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.563e+00 (tol = 1.000e-10) r (rel) = 9.883e-02 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 1.169e-02 (tol = 1.000e-10) r (rel) = 7.395e-04 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 1.215e-06 (tol = 1.000e-10) r (rel) = 7.685e-08 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 7.951e-14 (tol = 1.000e-10) r (rel) = 5.028e-15 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
Custom Solver:
Iteration: 0 ; Error: 1.581e+01
Iteration: 1 ; Error: 1.563e+00
Iteration: 2 ; Error: 1.169e-02
Iteration: 3 ; Error: 1.235e-06
Iteration: 4 ; Error: 6.880e-12
Now, my problem is when I put a non-zero Dirichlet condition at the right end. The default fenics solver works perfect:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 5.418e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 4.638e-01 (tol = 1.000e-10) r (rel) = 8.560e-03 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 5.118e-04 (tol = 1.000e-10) r (rel) = 9.445e-06 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 7.032e-10 (tol = 1.000e-10) r (rel) = 1.298e-11 (tol = 1.000e-09)
Newton solver finished in 3 iterations and 3 linear solver iterations.
But my custom solver fails:
Iteration: 0 ; Error: 5.418e+01
Iteration: 1 ; Error: 6.535e+01
Iteration: 2 ; Error: 8.172e+01
Iteration: 3 ; Error: 1.067e+02
Iteration: 4 ; Error: 1.474e+02
Iteration: 5 ; Error: 2.192e+02
I think my solver fails because I am missing something on how boundary conditions should be handled. In fact, in the update is something wrong because at each iteration the vector x is 0.1 at the prescribed dofs of the right end which makes that u gets bigger and bigger at each iteration. Before I was using
g = d.assemble(G)
dg = d.assemble(dG)
# Apply BC
for bc in bcs:
bc.apply(g)
bc.apply(dg)
instead of dg, g = d.assemble_system(dG, G, bcs), but that is even worst:
Iteration: 0 ; Error: 2.000e-01
Iteration: 1 ; Error: nan
So, I had two questions, the first one is if someone has an idea of what I am doing wrong and how to fix it.
The second one is why is different to use assemble_system instead of assemble and then apply in this case. I have been doing some digging in this last issue and assemble_system maintains the assembled matrix symmetric. Can someone explain to me how assemble_system does that?
Thank you in advance and sorry for the long post, but there were a lot of things.