I am implementing a phase-field model for brittle fracture and have run into a classic numerical stability issue where the solver converges perfectly during the initial elastic loading but consistently fails the moment the crack is expected to initiate and propagate.
I’m confident the physics formulation is correct, but my monolithic solver setup cannot handle the severe ill-conditioning of the Jacobian at the onset of fracture. I’m seeking advice on the recommended best practices and robust solver strategies for this class of problem in FEniCSx.
I am intending to solve a fully coupled, monolithic solver. I started with a direct LU solver, which worked perfectly up until the moment of fracture. When I switched to an iterative solver using a fieldsplit preconditioner (gamg for the displacement field and lu for the phase-field), it also failed at the same point, with the linear solver stalling. Even trying to stabilize the elasticity block by attaching a nullspace didn’t help with convergence. I appreciate if you can guide me to any better debuging or manual setup to solve the divergence problem of solver (SNES Diverged (Reason -3 or -6)).
# 1. Imports and setup
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx.io import VTXWriter
from dolfinx import mesh, fem, default_scalar_type, la
from dolfinx.fem.petsc import NonlinearProblem
from petsc4py import PETSc
# 2. Simulation parameters
# Geometry
width, height = 1.0, 2.0
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array([0, 0]), np.array([width, height])],
[50, 100],
cell_type=mesh.CellType.triangle)
# Material properties
E_mat, nu_mat = 21.7e9, 0.25
E_ref = E_mat
Gc = 8e4 / E_ref
E_norm = E_mat / E_ref
# Phase-field parameters
epsilon = 0.04 # Regularization length
eta = 1e-6 # Residual stiffness
m_penalty = 1.0e3
phi_m = 0.98
l_2 = 0.01
# 3. Function spaces and fields
V_u = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
V_phi = fem.functionspace(domain, ("Lagrange", 1))
# Solution fields
u = fem.Function(V_u, name="Displacement")
phi = fem.Function(V_phi, name="PhaseField")
phi_from_last_step = fem.Function(V_phi, name="PhaseFieldHistory")
# Test/trial functions
u_test, phi_test = ufl.TestFunction(V_u), ufl.TestFunction(V_phi)
u_trial, phi_trial = ufl.TrialFunction(V_u), ufl.TrialFunction(V_phi)
# 4. Variational formulation: Energies
def W0_iso(u):
"""Isotropic elastic energy."""
mu = E_norm / (2 * (1 + nu_mat))
lmbda = E_norm * nu_mat / ((1 + nu_mat) * (1 - 2 * nu_mat))
eps = ufl.sym(ufl.grad(u))
return 0.5 * lmbda * ufl.tr(eps)**2 + mu * ufl.inner(eps, eps)
def gamma_surface_energy(phi):
"""Crack surface energy."""
grad_phi = ufl.grad(phi)
return Gc * (phi**2 / (2 * epsilon) + (epsilon / 2) * ufl.dot(grad_phi, grad_phi))
def H_l(x, x_m, l_scale):
return 0.5 * (1 + ufl.tanh((x - x_m) / l_scale))
def irreversibility_penalty(phi, phi_prev):
"""Penalty to prevent crack healing."""
return m_penalty * H_l(phi_prev, phi_m, l_2) * (phi - phi_prev)**2
def total_energy(u, phi, phi_prev):
"""Total potential energy."""
g_degradation = (1 - phi)**2
elastic = (g_degradation + eta) * W0_iso(u)
surface = gamma_surface_energy(phi)
penalty = irreversibility_penalty(phi, phi_prev)
return elastic + surface + penalty
# 5. Boundary & initial conditions
def initial_crack(x):
"""Initializes a horizontal crack in the middle of the plate."""
in_crack_y = np.abs(x[1] - height / 2) < 0.01
in_crack_x = np.abs(x[0] - width / 2) <= width / 4
return np.where(in_crack_x & in_crack_y, 1.0, 0.0)
phi.interpolate(initial_crack)
phi_from_last_step.interpolate(initial_crack)
# Tensile test boundary conditions
load_val = fem.Constant(domain, default_scalar_type(0.0))
# Fix bottom boundary
bottom_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, lambda x: np.isclose(x[1], 0))
bottom_dofs_u0 = fem.locate_dofs_topological(V_u.sub(0), domain.topology.dim - 1, bottom_facets)
bottom_dofs_u1 = fem.locate_dofs_topological(V_u.sub(1), domain.topology.dim - 1, bottom_facets)
bc_bottom_x = fem.dirichletbc(default_scalar_type(0), bottom_dofs_u0, V_u.sub(0))
bc_bottom_y = fem.dirichletbc(default_scalar_type(0), bottom_dofs_u1, V_u.sub(1))
# Fix top boundary in x, apply displacement in y
top_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, lambda x: np.isclose(x[1], height))
top_dofs_u0 = fem.locate_dofs_topological(V_u.sub(0), domain.topology.dim - 1, top_facets)
top_dofs_u1 = fem.locate_dofs_topological(V_u.sub(1), domain.topology.dim - 1, top_facets)
bc_top_x = fem.dirichletbc(default_scalar_type(0), top_dofs_u0, V_u.sub(0))
bc_top_y = fem.dirichletbc(load_val, top_dofs_u1, V_u.sub(1))
bcs_u = [bc_bottom_x, bc_bottom_y, bc_top_x, bc_top_y]
def build_nullspace_2D(V):
"""Build PETSc nullspace for 2D elasticity (translations and rotation)."""
basis_funcs = [fem.Function(V) for _ in range(3)]
# Translations
basis_funcs[0].sub(0).interpolate(lambda x: np.full_like(x[0], 1.0))
basis_funcs[1].sub(1).interpolate(lambda x: np.full_like(x[0], 1.0))
# Rotation
basis_funcs[2].sub(0).interpolate(lambda x: -x[1])
basis_funcs[2].sub(1).interpolate(lambda x: x[0])
basis_dolfinx = [f.x for f in basis_funcs]
la.orthonormalize(basis_dolfinx)
basis_petsc = [vec.petsc_vec for vec in basis_dolfinx]
return PETSc.NullSpace().create(vectors=basis_petsc)
# 6. Solver configuration
dx = ufl.Measure("dx", domain=domain)
# Residual and Jacobian
Pi = total_energy(u, phi, phi_from_last_step) * dx
Res_u = ufl.derivative(Pi, u, u_test)
Res_phi = ufl.derivative(Pi, phi, phi_test)
Res_phi_scaled = (1 / (Gc / epsilon)) * Res_phi
J_uu = ufl.derivative(Res_u, u, u_trial)
J_u_phi = ufl.derivative(Res_u, phi, phi_trial)
J_phi_u = ufl.derivative(Res_phi_scaled, u, u_trial)
J_phi_phi = ufl.derivative(Res_phi_scaled, phi, phi_trial)
F_forms = [Res_u, Res_phi_scaled]
J_forms = [[J_uu, J_u_phi], [J_phi_u, J_phi_phi]]
# PETSc options for fieldsplit solver
petsc_options = {
"snes_type": "newtonls", "snes_linesearch_type": "bt", "snes_monitor": None,
"snes_max_it": 50, "snes_rtol": 1e-8,
"ksp_type": "gmres", "pc_type": "fieldsplit",
"ksp_rtol": 1e-6,
"fieldsplit_0_ksp_type": "preonly", "fieldsplit_0_pc_type": "gamg",
"fieldsplit_1_ksp_type": "preonly", "fieldsplit_1_pc_type": "lu",
}
# Create nonlinear problem and attach nullspace
prob_uphi = NonlinearProblem(F_forms, [u, phi], bcs=bcs_u, J=J_forms,
petsc_options=petsc_options,
petsc_options_prefix="uphi_snes",
kind="nest")
solver = prob_uphi.solver
ksp = solver.getKSP()
pc = ksp.getPC()
pc.setUp() # Must be called before getting sub-KSPs
ksp_subs = pc.getFieldSplitSubKSP()
ns = build_nullspace_2D(V_u)
ksp_subs[0].getPC().getOperators()[0].setNearNullSpace(ns)
# 7. Main loading loop
T_total = 1.0e-2
delta_T = 1.0e-4
u_pull = 1.0 # Apply displacement upwards
num_steps = int(T_total / delta_T)
print("Starting simulation...")
# Setup output file
vtx_writer = VTXWriter(domain.comm, "solution.bp", [u, phi], engine="BP4")
vtx_writer.write(0.0)
for i in range(num_steps):
t = (i + 1) * delta_T
load_val.value = t * u_pull
if domain.comm.rank == 0:
print(f"\n--- Step {i+1}/{num_steps}, Applied Disp = {load_val.value:.3e} ---")
# Solve the system
prob_uphi.solve()
reason = solver.getConvergedReason()
iters = solver.getIterationNumber()
if domain.comm.rank == 0:
print(f"SNES converged in {iters} iterations with reason {reason}.")
if reason < 0:
if domain.comm.rank == 0:
print("!!! SNES Diverged. Aborting.")
break
# Make sure crack growth is irreversible
phi.x.array[:] = np.maximum(phi.x.array, phi_from_last_step.x.array)
phi.x.scatter_forward()
# Update history for the next step
phi_from_last_step.x.array[:] = phi.x.array
vtx_writer.write(t)
vtx_writer.close()
if domain.comm.rank == 0:
print("\nSimulation finished.")