I am working on setting up a custom Newton Solver framework for a contact problem that I am solving, but wanted to validate it with a simpler problem first. I adopted the code from the Custom Newton Solver tutorial for a hyperelasticity problem that converges well with the built-in Newton solver. However, the update step that I am calculating is very small (~1e-6 even at the first step), and the solver takes far more steps than the built in solver. When it does finish (based on the magnitude of the step size), the solution is inconsistent with the built-in solver (and is illogical - no displacement except on single faces).
Are there particular changes that need to be made to the Custom Newton Solver tutorial to adapt it to vector-valued function problems?
I have included the code for the hyperelasticity problem below with both solvers, as well as the solutions that are produced. (All code is using dolfin-x in Docker)
import numpy as np
import ufl
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
L=20.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0],[L,1,1]],[20,5,5],mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain,("CG",2))
def left(x):
return np.isclose(x[0],0)
def right(x):
return np.isclose(x[0],L)
fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain,fdim,left)
right_facets = mesh.locate_entities_boundary(domain,fdim,right)
marked_facets = np.hstack([left_facets,right_facets])
marked_values = np.hstack([np.full(len(left_facets),1,dtype=np.int32),np.full(len(right_facets),2,dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain,fdim,marked_facets[sorted_facets],marked_values[sorted_facets])
u_bc = np.array((0,)*domain.geometry.dim,dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V,facet_tag.dim,facet_tag.indices[facet_tag.values==1])
right_dofs = fem.locate_dofs_topological(V,facet_tag.dim,facet_tag.indices[facet_tag.values==2])
bcs = [fem.dirichletbc(u_bc,left_dofs,V),fem.dirichletbc(u_bc,right_dofs,V)]
B = fem.Constant(domain,PETSc.ScalarType((0,0,-1)))
T = fem.Constant(domain,PETSc.ScalarType((0,0,0)))
v,u = ufl.TestFunction(V),fem.Function(V)
d = len(u)
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.49)
mu = fem.Constant(domain,E/(2*(1+nu)))
lmbda = fem.Constant(domain,E*nu/((1+nu)*(1-2*nu)))
psi = (mu/2)*(Ic-3) - mu*ufl.ln(J) + (lmbda/2)*(ufl.ln(J))**2
P = ufl.diff(psi,F)
metadata = {"quadrature_degree":4}
ds = ufl.Measure('ds',domain=domain,subdomain_data=facet_tag,metadata=metadata)
dx = ufl.Measure("dx",domain=domain,metadata=metadata)
PI = ufl.inner(P,ufl.grad(v))*dx - ufl.inner(B,v)*dx - ufl.inner(T,v)*ds(2)
##### Choose which solver to use (0: built-in solver, 1: custom Newton solver) #####
solver_tool = ['built_in','custom'][1]
if solver_tool=='built_in':
problem = fem.petsc.NonlinearProblem(PI,u,bcs=bcs)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)
solver.atol = 1e-12
solver.rtol = 1e-12
solver.convergence_criterion = "incremental"
elif solver_tool=='custom':
### Custom Solver ###
F_form = fem.form(PI) # form of the residual function
J = ufl.derivative(PI, u) # jacobian of the residual function
J_form = fem.form(J) # form of the jacobian
du = fem.Function(V) # increment of the solution
u_start = fem.Function(V) # holds previous solution for reverting in case of convergence failure
A = fem.petsc.create_matrix(J_form) # preallocating sparse jacobian matrix
L = fem.petsc.create_vector(F_form) # preallocating residual vector
solver = PETSc.KSP().create(domain.comm)
opts = PETSc.Options()
option_prefix = solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
solver.setFromOptions()
solver.setOperators(A)
# iteration parameters
relaxation_param = 1 # full Newton's Method (<1 for restraint)
max_it = int(1e3)
tol = 1e-8
### iteration adopted from Custom Newton Solver tutorial ###
def NewtonSolve(print_steps=True,print_solution=True):
i = 0 # number of iterations of the Newton solver
converged = False
while i<max_it:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, J_form, bcs=bcs)
A.assemble()
fem.petsc.assemble_vector(L, F_form)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [J_form], [bcs], x0=[u.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, bcs, u.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linear problem
solver.solve(L, du.vector)
du.x.scatter_forward()
# Update u_{i+1} = u_i + relaxation_param * delta x_i
u.x.array[:] += du.x.array[:]
i += 1
# Compute norm of update
correction_norm = du.vector.norm(0)
error_norm = L.norm(0)
if print_steps:
print(f" Iteration {i}: Correction norm {correction_norm}, Residual {error_norm}")
if correction_norm < tol:
converged = True
break
if print_solution:
if converged:
print(f"Solution reached in {i} iterations.")# (Residual norm {error_norm})")
else:
print(f"No solution found after {i} iterations. Revert to previous solution and adjust solver parameters.")
return converged, u_start, i
import pyvista
pyvista.set_jupyter_backend("ipygany")
grid = pyvista.UnstructuredGrid(*plot.create_vtk_mesh(domain,domain.topology.dim))
def plot_function(t,uh):
p = pyvista.Plotter()
topology,cells,geometry = plot.create_vtk_mesh(uh.function_space)
function_grid = pyvista.UnstructuredGrid(topology,cells,geometry)
var_name = f"u({t})"
values = np.zeros((geometry.shape[0],3))
values[:,:len(uh)] = uh.x.array.reshape(geometry.shape[0],len(uh))
function_grid[var_name] = values
function_grid.set_active_vectors(var_name)
warped = function_grid.warp_by_vector(var_name,factor=1)
actor = p.add_mesh(warped)
p.show_axes()
if not pyvista.OFF_SCREEN:
p.show()
else:
pyvista.start_xvfb()
figure_as_array = p.screenshot(f"Hyperelasticity_{t:.2f}.png")
p.remove_actor(actor)
from dolfinx import log
log.set_log_level(log.LogLevel.INFO)
print('Using',solver_tool,'Newton solver')
if solver_tool=='built_in':
num_its,converged = solver.solve(u)
elif solver_tool=='custom':
converged, u_start, num_its = NewtonSolve()
assert(converged)
print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")
plot_function(n,u)
Solution Comparison:
Any help would be greatly appreciated!