How to Calculate Error Norms and Check Convergence in PhaseFieldX

PhaseFieldX provides dedicated modules in src/phasefieldx/norms.py and src/phasefieldx/errors_functions.py to compute L², H¹, and Lᵖ error norms, while convergence monitoring is handled by the Newton solver in src/phasefieldx/solvers/newton.py using PETSc SNES residual norms.

PhaseFieldX is an open-source finite element framework built on DOLFINx for simulating phase-field fracture and multi-physics problems. When validating numerical solutions against analytical references or debugging solver behavior, you need rigorous tools to calculate error norms and check convergence of nonlinear iterations.

Computing Error Norms in PhaseFieldX

Core Norm Implementations in norms.py

The foundation for all error calculations resides in src/phasefieldx/norms.py. These functions implement standard functional norms using DOLFINx forms with MPI-aware reduction, returning a single global float across all ranks.

• norm_L2(field, msh, dx=ufl.dx) computes $\sqrt{\int |field|^2 , dx}$
• norm_Lp(field, msh, p=2, dx=ufl.dx) computes $(\int |field|^p , dx)^{1/p}$
• norm_semiH1(field, msh, dx=ufl.dx) computes $\sqrt{\int |\nabla field|^2 , dx}$
• norm_H1(field, msh, dx=ufl.dx) computes $\sqrt{\int |field|^2 , dx + \int |\nabla field|^2 , dx}$

These are pure wrappers around DOLFINx scalar assembly, ensuring accurate global norms across distributed meshes.

High-Level Error Functions in errors_functions.py

For practical verification, src/phasefieldx/errors_functions.py provides convenience wrappers that accept two dolfinx.Function objects—your numerical solution uh and an exact/reference solution u_ex—and delegate to the norm utilities:

• compute_error_L2, compute_error_H1, compute_error_semiH1 call the corresponding norm on the difference field_a - field_b.
• eval_error_LP supports arbitrary $L^p$ norms.
• error_L2_higher_order_space builds a higher-order function space, interpolates both fields there, and integrates the squared difference—essential when the exact solution lives in a richer space than the FE approximation.
• eval_error_L2_normalized returns the relative L² error $|u_h - u_{ex}|{L^2} / |u{ex}|_{L^2}$ with safeguards against division-by-zero.

Monitoring Solver Convergence

Newton Solver Configuration in newton.py

Nonlinear solver convergence is managed by the NewtonSolver class in src/phasefieldx/solvers/newton.py, which wraps PETSc's SNES (Scalable Nonlinear Equations Solvers). The wrapper exposes three critical attributes for controlling convergence criteria:

• solver.rtol – relative tolerance for the residual norm.
• solver.atol – absolute tolerance for the residual norm.
• solver.convergence_criterion – defaults to "residual", meaning SNES reports convergence when the residual norm falls below the specified tolerances.

Setting solver.report = True enables verbose logging of the residual norm at each iteration via logger.info calls.

Residual Norm Logging and Inspection

During each SNES iteration, the Newton solver retrieves the residual norm via snes.getFunctionNorm(). This value is printed to the log when reporting is enabled, allowing you to track how quickly the nonlinear residual decreases. For phase-field fracture simulations, the residual history is often written to .conv files, which can be post-processed using scripts like examples/PhaseFieldFracture/plot_1718.py to generate convergence plots.

Complete Workflow Example

The following example demonstrates how to calculate error norms and check convergence in a complete PhaseFieldX simulation workflow:

import dolfinx.mesh, dolfinx.fem
import ufl, numpy as np
from mpi4py import MPI
from phasefieldx.errors_functions import compute_error_L2, compute_error_H1
from phasefieldx.norms import norm_L2, norm_H1
from phasefieldx.solvers.newton import NewtonSolver

# ----------------------------------------------------------------------

# 1️⃣  Set up mesh, function spaces and exact solution

# ----------------------------------------------------------------------

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))

# Numerical solution (already computed by a solver)

uh = dolfinx.fem.Function(V)

# ... fill `uh` via your simulation ...

# Exact solution (example: u(x,y) = sin(pi*x)*sin(pi*y))

exact_expr = ufl.sin(np.pi * ufl.x) * ufl.sin(np.pi * ufl.y)
u_ex = dolfinx.fem.Function(V)
u_ex.interpolate(dolfinx.fem.Expression(exact_expr, V.element.interpolation_points))

# ----------------------------------------------------------------------

# 2️⃣  Compute error norms

# ----------------------------------------------------------------------

err_L2 = compute_error_L2(uh, u_ex, mesh)          # → float

err_H1 = compute_error_H1(uh, u_ex, mesh)          # → float

rel_err = err_L2 / norm_L2(u_ex, mesh)             # manual relative error

print(f"L2 error = {err_L2:.3e}, H1 error = {err_H1:.3e}, Relative L2 = {rel_err:.3e}")

# ----------------------------------------------------------------------

# 3️⃣  Run Newton solver and monitor convergence

# ----------------------------------------------------------------------

# Assuming `problem` supplies residual and Jacobian forms

solver = NewtonSolver(problem)
solver.solver.rtol = 1e-8                           # tighter relative tolerance

solver.solver.atol = 1e-9                           # absolute tolerance

solver.solver.report = True                         # prints residual each iteration

solver.solve()                                      # logs "Residual norm u: …"

Summary

• Error norms in PhaseFieldX are computed via src/phasefieldx/norms.py, which provides MPI-aware implementations of L², Lᵖ, H¹-semi, and full H¹ norms.
• Error calculation between numerical and exact solutions is streamlined through src/phasefieldx/errors_functions.py, supporting absolute, relative, and higher-order space projections.
• Convergence monitoring is built into the NewtonSolver class in src/phasefieldx/solvers/newton.py, leveraging PETSc SNES residual norms with configurable relative and absolute tolerances.

Frequently Asked Questions

What norm types does PhaseFieldX support?

PhaseFieldX supports L², Lᵖ (for arbitrary p), H¹-semi, and full H¹ norms through the functions in src/phasefieldx/norms.py. These norms are computed using DOLFINx forms with MPI-aware reduction, ensuring accurate global values across distributed meshes.

How do I compute relative error norms?

Use the eval_error_L2_normalized function from src/phasefieldx/errors_functions.py, which returns $|u_h - u_{ex}|{L^2} / |u{ex}|_{L^2}$ with built-in safeguards against division-by-zero. Alternatively, manually divide the absolute error computed by compute_error_L2 by the norm of the exact solution using norm_L2(u_ex, mesh).

Can I monitor convergence without using the Newton solver?

While the built-in convergence monitoring is integrated into NewtonSolver via PETSc SNES, you can implement custom convergence checks by evaluating the residual norm directly within your own solver loop. The snes.getFunctionNorm() method (or equivalent DOLFINx/PETSc calls) provides the residual magnitude needed to assess convergence criteria manually.

Where are convergence histories stored in PhaseFieldX?

For phase-field fracture simulations, convergence histories are typically written to .conv files during the simulation. These files can be post-processed using scripts such as examples/PhaseFieldFracture/plot_1718.py to generate visualizations of residual norms versus iteration count or time steps.

{
  "mcpServers": {
    "instagit": {
      "command": "npx",
      "args": ["-y", "instagit@latest"]
    }
  }
}