Non-convergence issue in the Newton solver while solving the hydrogen diffusion equation in a phase-field model for hydrogen-assisted cracking under tensile loading

I am implementing a phase-field model for hydrogen-assisted cracking under tensile load in FEniCS. The model includes displacement, phase-field, and hydrogen diffusion equations. However, I am encountering a Newton solver non-convergence issue while solving the hydrogen diffusion equation. The error message is:
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
o address this, I have:

  • Used MUMPS as the linear solver.
  • Set strict absolute/relative tolerances.
  • Reduced the time step dynamically when the residual is too large.
  • Initialized the hydrogen concentration (C.interpolate(Constant(0.1))) before solving.

Despite these steps, the solver still fails after several iterations. What strategies or solver settings can I try to improve convergence?

Here is the code snippet.

# Hydrogen diffusion (nonlinear problem with stress-driven term)
E_C = (( (C - C_old) / Constant(deltaT) ) * q * dx
       + D_eff * dot(grad(C), grad(q)) * dx
       - (D_eff * VH * C) / (R * T) * dot(grad(sigma_H), grad(q)) * dx
       + 1e-3 * dot(grad(C), grad(q)) * dx)  # Preconditioning term

# Jacobian for nonlinear solver
J_C = derivative(E_C, C)

# Define nonlinear problem and solver
problem_C = NonlinearVariationalProblem(E_C, C, bc_C, J_C)
solver_C = NonlinearVariationalSolver(problem_C)

# Solver parameters
solver_C.parameters["newton_solver"]["linear_solver"] = "mumps"
solver_C.parameters["newton_solver"]["absolute_tolerance"] = 1e-10
solver_C.parameters["newton_solver"]["relative_tolerance"] = 1e-8
solver_C.parameters["newton_solver"]["maximum_iterations"] = 100  # Increase if needed
solver_C.parameters["newton_solver"]["relaxation_parameter"] = 0.5  # Reduce step size
solver_C.parameters["newton_solver"]["error_on_nonconvergence"] = False  # Prevent crash

# Time-stepping loop
for step in range(num_steps):
    t += deltaT
    tension.t = t

    solver_disp.solve()
    u_old.assign(u)

    Hold.assign(project(psi(u), V))
    phi_old.assign(phi)
    solver_phi.solve()

    # Check residual before solving for C
    residual_C = assemble(E_C)
    residual_norm_C = norm(residual_C, 'l2')

    if residual_norm_C > 1e6:
        print(f"Warning: Large residual ({residual_norm_C}) in hydrogen diffusion solver at step {step}, reducing Δt.")
        deltaT *= 0.5  # Dynamically reduce time step if divergence occurs
        if deltaT < 1e-6:
            print("Time step too small, stopping simulation.")
            break

    # Ensure C starts from a reasonable value
    C.interpolate(Constant(0.1))
    C_old.assign(C)

    solver_C.solve()

Could you maybe provide a full MWE so that it can be run by other people?
But just looking at this, I’m wondering, shouldn’t you not update your C_old when your Newton solver doesn’t converge and you reduce the time step?

1 Like

Thanks for your response!
When the Newton solver fails to converge and we reduce the time step (Δt), we should NOT update C_old. Otherwise, we are moving forward with a bad solution, making convergence even harder.
Please correct me if i am wrong.

Here is my code for your kind perusal.
Thanks in advance.

from dolfin import *
from mshr import *
import numpy as np

# Define and refine the mesh
mesh = RectangleMesh(Point(-0.5, -0.5), Point(0.5, 0.5), 200, 200, "crossed")

# Refine mesh around the initial horizontal crack (y=0, x <=0)
class CrackDomain(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1]) < 0.01 and x[0] <= 0.0  # Horizontal crack at y=0

crack_domain = CrackDomain()
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)
crack_domain.mark(cell_markers, True)
mesh = refine(mesh, cell_markers)

# Function spaces
V = FunctionSpace(mesh, 'CG', 1)  # Phase-field
W = VectorFunctionSpace(mesh, 'CG', 1)  # Displacement
Q = FunctionSpace(mesh, 'CG', 1)  # Hydrogen concentration

# Define trial and test functions
du, v = TrialFunction(W), TestFunction(W)
dphi, w = TrialFunction(V), TestFunction(V)
dC, q = TrialFunction(Q), TestFunction(Q)

# Functions
u, u_old = Function(W), Function(W)
phi, phi_old = Function(V), Function(V)
C, C_old = Function(Q), Function(Q)
Hold = Function(V)

# Material parameters
E = 210e9  # Young's modulus (Pa)
nu = 0.3   # Poisson's ratio
lmbda = E * nu / ((1 + nu) * (1 - 2*nu))  # Lamé parameter λ
mu = E / (2 * (1 + nu))  # Lamé parameter μ
Gc_0 = 2.7  # Initial fracture energy (N/m)
l = 0.001   # Phase-field length scale (m)
D0 = 1e-9   # Baseline diffusion coefficient (m²/s)
VH = 2e-6   # Partial molar volume (m³/mol)
R = 8.314   # Gas constant (J/(mol·K))
T = 300     # Temperature (K)
chi = 0.89  # Hydrogen embrittlement coefficient
Delta_g0 = 30000  # Gibbs free energy difference (J/mol)

# Boundary conditions
bcbot = DirichletBC(W, Constant((0.0, 0.0)), "near(x[1], -0.5) && on_boundary")
tension = Expression(("0.0", "u_r * t / total_time"), u_r=0.007, t=0.0, total_time=3.0, degree=1)
bctop = DirichletBC(W, tension, "near(x[1], 0.5) && on_boundary")
bc_u = [bcbot, bctop]

bc_C = DirichletBC(Q, Constant(0.1), "on_boundary")

class Crack(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1]) < 1e-3 and x[0] <= 0.0

crack_surface = Crack()
bc_phi = DirichletBC(V, Constant(1.0), crack_surface)

# Initial crack
class InitialCrack(UserExpression):
    def eval(self, values, x):
        if abs(x[1]) < 1e-3 and x[0] <= 0.0:
            values[0] = 1.0
        else:
            values[0] = 0.0
    def value_shape(self):
        return ()

phi.interpolate(InitialCrack())

# Constitutive relations
def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return 2 * mu * epsilon(u) + lmbda * tr(epsilon(u)) * Identity(len(u))

def psi(u):
    return 0.5 * lmbda * (tr(epsilon(u)))**2 + mu * inner(epsilon(u), epsilon(u))

sigma_H = tr(sigma(u))

theta = C / (C + exp(-Delta_g0 / (R * T)))
Gc = Gc_0 * (1 - chi * theta)
D_eff = (1 - chi * theta) * D0

# Variational forms
# Displacement
E_du = ((1.0 - phi_old**2) * inner(sigma(du), grad(v)) - inner(Constant((0.0, 0.0)), v)) * dx

# Phase-field
def H(u_old, u_new, Hold):
    psi_new = psi(u_new)
    return conditional(gt(psi_new, Hold), psi_new, Hold)

F_phi = (Gc * l * inner(grad(phi), grad(w)) + (Gc / l) * phi * w - 2 * (1 - phi_old) * H(u_old, u, Hold) * w) * dx
J_phi = derivative(F_phi, phi)

# Hydrogen diffusion
F_C = (( (C - C_old)/Constant(0.01) ) * q * dx
       + D_eff * dot(grad(C), grad(q)) * dx
       - (D_eff * VH * C) / (R * T) * dot(grad(sigma_H), grad(q)) * dx )
J_C = derivative(F_C, C)

# Solvers
problem_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), u, bc_u)
solver_disp = LinearVariationalSolver(problem_disp)

problem_phi = NonlinearVariationalProblem(F_phi, phi, bc_phi, J=J_phi)
solver_phi = NonlinearVariationalSolver(problem_phi)
prm_phi = solver_phi.parameters
prm_phi['newton_solver']['linear_solver'] = 'gmres'
prm_phi['newton_solver']['preconditioner'] = 'ilu'
prm_phi['newton_solver']['absolute_tolerance'] = 1e-10
prm_phi['newton_solver']['relative_tolerance'] = 1e-8

problem_C = NonlinearVariationalProblem(F_C, C, bc_C, J=J_C)
solver_C = NonlinearVariationalSolver(problem_C)
prm_C = solver_C.parameters
prm_C['newton_solver']['linear_solver'] = 'gmres'
prm_C['newton_solver']['preconditioner'] = 'ilu'
prm_C['newton_solver']['absolute_tolerance'] = 1e-10
prm_C['newton_solver']['relative_tolerance'] = 1e-8

# Time-stepping parameters
t = 0.0
deltaT = 0.001
num_steps = 3000

# Output
conc_file = XDMFFile("results.xdmf")
conc_file.parameters["flush_output"] = True

# Time loop
max_attempts = 5
successful_step = False
current_deltaT = deltaT
for step in range(num_steps):
    t += current_deltaT
    tension.t = t
    
    attempts = 0
    while attempts < max_attempts:
        try:
            # Solve displacement using u_old (previous step's solution)
            solver_disp.solve()
            
            # Update history variable using u and u_old
            Hold.assign(project(H(u_old, u, Hold), V))
            
            # Solve phase-field using phi_old (previous) and u (current)
            solver_phi.solve()
            
            # Solve hydrogen diffusion using C_old (previous) and current variables
            solver_C.solve()
            
            # Update old variables for next time step
            u_old.assign(u)
            phi_old.assign(phi)
            C_old.assign(C)
            
            successful_step = True
            break  # Break out of the retry loop if all solves succeed
        except RuntimeError as e:
            print(f"Solver failed at time {t} with deltaT {current_deltaT}. Retrying with reduced time step.")
            current_deltaT *= 0.5  # Reduce time step by half
            attempts += 1
    
    if not successful_step:
        print(f"Failed to converge after {max_attempts} attempts. Terminating simulation.")
        break

    # Write output
    C.rename("HydrogenConcentration", "label")
    phi.rename("PhaseField", "label")
    conc_file.write(C, t)
    conc_file.write(phi, t)

conc_file.close()
print("Simulation completed.")