Newton Solver Does Not Converge in Mooney-Rivlin Hyperelasticity Example - Maximum Iterations Reached

Errors
1

Hello,

I am working on a Mooney-Rivlin hyperelasticity problem using FEniCSx and encountering convergence issues when solving the nonlinear problem with the Newton solver.

I have a simple code setup for the problem, which defines a unit square mesh, applies Dirichlet boundary conditions, and uses a Mooney-Rivlin strain energy model. However, I am unable to get the Newton solver to converge.

Code:

import dolfinx
from mpi4py import MPI
import basix
import numpy as np
import ufl
from dolfinx.nls.petsc import NewtonSolver

# Create a unit square mesh
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 4, 4, dolfinx.mesh.CellType.quadrilateral)
gdim = mesh.geometry.dim

# Define the displacement element for finite elements
element_disp = basix.ufl.element("Lagrange", mesh.topology.cell_name(), degree=1, shape=(gdim,))
V = dolfinx.fem.functionspace(mesh, element_disp)
δ_u = ufl.TestFunction(V)
u = dolfinx.fem.Function(V, name="Displacement")

# Boundary conditions
def left(x): return np.isclose(x[0], 0)
def bottom(x): return np.isclose(x[1], 0)
def top(x): return np.isclose(x[1], 1)
def right(x): return np.isclose(x[0], 1)

fdim = mesh.topology.dim - 1
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, left)
right_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, right)
bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom)
top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top)

marked_facets = np.hstack([left_facets, right_facets, bottom_facets, top_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2), np.full_like(bottom_facets, 3), np.full_like(top_facets, 4)])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

# Define measures for integration
metadata = {"quadrature_degree": 4}
dx = ufl.Measure("dx", domain=mesh, metadata=metadata)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag, metadata=metadata)

# Boundary conditions
fixed_bottom_dofs = dolfinx.fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(3))
fixed_right_dofs = dolfinx.fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(2))
fixed_bottom = dolfinx.fem.dirichletbc(0.0, fixed_bottom_dofs, V.sub(1))
fixed_right = dolfinx.fem.dirichletbc(0.0, fixed_right_dofs, V.sub(0))

u_load = 1.0
top_boundary_dofs = dolfinx.fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(4))
top_displacement = dolfinx.fem.dirichletbc(u_load, top_boundary_dofs, V.sub(1))

# Set boundary conditions
bc_u = [fixed_bottom, fixed_right, top_displacement]

# Material parameters for Mooney-Rivlin model
K = 21.01
C_01 = 7.003
C_10 = 2.617

dim = len(u)
I = ufl.Identity(dim)
F = ufl.variable(I + ufl.grad(u))
J = ufl.det(F)

# Strain energy function for Mooney-Rivlin model
def invariants(u):
    C = ufl.variable(ufl.dot(F.T, F))
    I1 = ufl.variable(ufl.tr(C))
    I2 = ufl.variable(0.5 * (I1**2 - ufl.tr(ufl.dot(C, C))))
    I1_bar = ufl.variable((J**(-2/3)) * I1)
    I2_bar = ufl.variable((J**(-4/3)) * I2)
    return {'C': C, 'I1_bar': I1_bar, 'I2_bar': I2_bar}

def mooney_rivlin_strain_energy(u):
    invariant = invariants(u)
    I1_bar, I2_bar = invariant['I1_bar'], invariant['I2_bar']
    psi_iso = C_10 * (I1_bar - 3) + C_01 * ((I2_bar**(3/2)) - 3**(3/2))
    psi_vol = (K/2) * ((J-1)**2)
    H = ufl.conditional(ufl.ge(J,0), 1, 0)
    psi_vol_plus = H * psi_vol
    psi_mr_plus = (psi_iso + psi_vol_plus)
    return psi_mr_plus

def stress(u):
    psi_mr_plus = mooney_rivlin_strain_energy(u)
    P_plus = ufl.diff(psi_mr_plus, F)
    return P_plus

P = stress(u)
weak_form_disp = ufl.inner(ufl.grad(δ_u), P) * dx

# Set up and solve the nonlinear problem
problem = dolfinx.fem.petsc.NonlinearProblem(weak_form_disp, u, bcs=bc_u)
solver = NewtonSolver(mesh.comm, problem)
solver.atol = 1e-8  
solver.rtol = 1e-8  
solver.convergence_criterion = "incremental"

# Solve the problem
num_its, converged = solver.solve(u)
assert (converged)
u.x.scatter_forward()

Error Message:

num_its, converged = solver.solve(u)

  File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/nls/petsc.py:47 in solve
    n, converged = super().solve(u.x.petsc_vec)

RuntimeError: Newton solver did not converge because maximum number of iterations reached

I have tried adjusting the tolerance and the solver parameters, but the problem still does not converge. Could anyone help me understand why the Newton solver is not converging, or suggest potential modifications to improve convergence?

Thank you for your help!

2

I don’t have the bandwidth to try your problem, but a first and straightforward avenue of investigation would be a continuation method.

One could start by using a value of u_load = 0.01 and incrementing over a number of discretized steps to 1.0 using the previous solution as the initial guess for each subsequent solve to see if you can achieve your desired setup.

Hyperelasticity problems are notorious for their multitude of solutions for seemingly simple setups. Defcon exemplifies this beautifully in, for example, its hyperelasticity demo.

3

To debug this, one should start with adding:

dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

which will give you info about the convergence:

[2024-12-03 20:50:59.115] [info] Column ghost size increased from 0 to 0
[2024-12-03 20:50:59.117] [info] PETSc Krylov solver starting to solve system.
[2024-12-03 20:50:59.118] [info] PETSc Krylov solver starting to solve system.
[2024-12-03 20:50:59.119] [info] Newton iteration 2: r (abs) = -nan (tol = 1e-08), r (rel) = -nan (tol = 1e-08)

Then, next I would change
solver.convergence_criterion = "incremental"
to
solver.convergence_criterion = "residual"
which gives:

[2024-12-03 20:51:46.479] [info] Column ghost size increased from 0 to 0
[2024-12-03 20:51:46.480] [info] Newton iteration 0: r (abs) = 35.740443828512376 (tol = 1e-08), r (rel) = inf (tol = 1e-08)
[2024-12-03 20:51:46.480] [info] PETSc Krylov solver starting to solve system.
[2024-12-03 20:51:46.481] [info] Newton iteration 1: r (abs) = -nan (tol = 1e-08), r (rel) = -nan (tol = 1e-08)

Next, I would modify the krylov solver for the underlying system, i.e.

solver.convergence_criterion = "residual"
solver.error_on_nonconvergence = False
from petsc4py import PETSc
ksp = solver.krylov_solver
options = {"ksp_type": "ksp", "pc_type": "lu", "pc_factor_mat_solver_type": "superlu",
           "ksp_error_if_not_converged": None, "ksp_monitor": None}
opts = PETSc.Options()
for key, val in options.items():
    opts.setValue(key, val)
ksp.setFromOptions()

However, I am not sure the variational form is correct, as

looks strange to me.
However, hyperelastic models is not my field of expertise.
There are other posts regarding Mooney-Rivlin at: Hyperelastic beam
which seems to use a very different variational formulation than you.
Similarly, see: 2D Hyper-elasticity — FEniCS solid tutorial 2019.0.2 documentation

1 Like