How the Newton Solver Works Internally Within PhaseFieldX: PETSc Configuration and Implementation

The Newton solver in PhaseFieldX is a thin wrapper around dolfinx's PETSc-based NewtonSolver that configures convergence tolerances, Krylov sub-solvers, and damping parameters to solve nonlinear phase-field problems.

PhaseFieldX is an open-source finite element framework for phase-field modeling of fracture, elasticity, and Allen-Cahn dynamics. At its core, the library relies on a robust Newton solver to handle the nonlinear systems arising from these physics. This article examines how the Newton solver works internally within PhaseFieldX, tracing the implementation from the PETSc wrapper in src/phasefieldx/solvers/newton.py to its usage in variational fracture solvers.

Newton Solver Architecture and PETSc Integration

The NewtonSolver class in src/phasefieldx/solvers/newton.py does not implement Newton's method from scratch. Instead, it wraps dolfinx's dolfinx.nls.petsc.NewtonSolver, which itself interfaces with PETSc's SNES (Scalable Nonlinear Equations Solvers) framework.

At initialization, the wrapper instantiates the underlying PETSc solver:

self.solver = dolfinx.nls.petsc.NewtonSolver(
    mpi4py.MPI.COMM_WORLD, problem)

This construction (line 30 in newton.py) accepts an MPI communicator and a dolfinx.fem.petsc.NewtonSolverNonlinearProblem object. By using MPI.COMM_WORLD, the solver operates in parallel by default, distributing the Jacobian and residual assemblies across processes.

Configuring Convergence Criteria and Solver Parameters

Immediately after construction, the wrapper applies a consistent set of convergence parameters suitable for phase-field problems. These settings are hardcoded in the __init__ method (lines 34-40):

Setting	Value	Meaning
max_it	500	Maximum Newton iterations per load step
rtol	1e-8	Relative tolerance on the Newton residual
atol	1e-9	Absolute tolerance on the Newton residual
convergence_criterion	"residual"	Use residual norm to decide convergence
report	True	Print convergence report each step
relaxation_parameter	1.0	Full Newton step (no damping)

The relative tolerance (rtol) of 1e-8 and absolute tolerance (atol) of 1e-9 provide tight convergence for the stiff nonlinearities typical in fracture mechanics. The convergence criterion set to "residual" uses the norm of the residual vector to determine convergence, as opposed to "incremental" which would check the solution update norm. The relaxation parameter of 1.0 indicates no line search or damping is applied by default, though users can adjust this for highly nonlinear problems.

Krylov Sub-Solver (KSP) Configuration

Each Newton iteration requires solving a linear system involving the Jacobian matrix. PhaseFieldX exposes the underlying Krylov solver (KSP) from PETSc for this purpose. The configuration occurs in lines 42-51 of newton.py:

self.ksp = self.solver.krylov_solver
self.opts = petsc4py.PETSc.Options()
self.option_prefix = self.ksp.getOptionsPrefix()

# Example user overrides (commented in source):

# self.opts[f"{self.option_prefix}ksp_type"] = "gmres"

# self.opts[f"{self.option_prefix}pc_type"] = "gamg"

self.ksp.setFromOptions()
self.ksp.setConvergenceHistory()

The wrapper retrieves the krylov_solver attribute from the dolfinx Newton solver, which returns a petsc4py.PETSc.KSP object. It then creates a PETSc options database to allow runtime configuration. The commented lines in the source demonstrate how users could override the default linear solver to use GMRES with a GAMG (Geometric Algebraic Multigrid) preconditioner, though by default the solver relies on PETSc's default choices (typically CG with Jacobi or ILU for symmetric systems).

Calling setFromOptions() applies any command-line PETSc options (e.g., -ksp_type gmres), while setConvergenceHistory() enables storage of the residual history for debugging.

Solving Nonlinear Problems: The Execution Flow

When solving a concrete phase-field problem, the workflow follows a three-step pattern demonstrated in src/phasefieldx/Element/Phase_Field_Fracture/solver/solver_ener_variational.py (lines 78-86):

# 1. Build the nonlinear problem object

problem = dolfinx.fem.petsc.NewtonSolverNonlinearProblem(
    residual_form, Φ, bcs=bc_list_phi)

# 2. Wrap it with PhaseFieldX's NewtonSolver

solver_phi = NewtonSolver(problem)

# 3. Execute the Newton iteration

phi_iterations, converged = solver_phi.solver.solve(Φ)

The underlying PETSc Newton solver performs the standard Newton-Krylov loop:

1. Assemble the residual F at the current iterate.
2. Assemble the Jacobian J (or apply a matrix-free action).
3. Solve J·Δx = -F using the configured KSP.
4. Update the iterate x ← x + λ·Δx where λ is the relaxation parameter.
5. Test convergence using the residual norm against rtol and atol.

The loop repeats until convergence or until the maximum iteration count is reached. The solution is stored in-place in the function Φ, and the method returns the iteration count and a boolean convergence flag.

Logging and Solver Inspection

PhaseFieldX provides utilities to inspect the solver configuration at runtime. The save_log_info method (lines 63-71 in newton.py) writes the current Newton and KSP settings to a standard Python logger:

def save_log_info(self, logger):
    logger.info(f"NewtonSolver: max_it={self.solver.max_it}, "
                f"rtol={self.solver.rtol}, atol={self.solver.atol}, "
                f"convergence_criterion={self.solver.convergence_criterion}")

Additionally, the __str__ method (lines 84-91) constructs a human-readable configuration string that queries the PETSc options database for the current KSP type and preconditioner:

def __str__(self):
    string = f"NewtonSolver: max_it={self.solver.max_it}, rtol={self.solver.rtol}, ..."
    string += f"\n  KSP type: {self.opts.get(self.option_prefix + 'ksp_type', 'default')}"
    string += f"\n  PC type: {self.opts.get(self.option_prefix + 'pc_type', 'default')}"
    return string

These inspection tools are invaluable for debugging convergence issues or verifying that command-line PETSc options have been correctly applied.

Integration with Blocked Solvers for Multiphysics

While the NewtonSolver wrapper handles single-field problems, PhaseFieldX also supports multi-field coupled problems (e.g., displacement and phase-field) through blocked Newton solvers. The variational fracture solver in solver_ener_variational.py demonstrates this pattern: it uses scifem.BlockedNewtonSolver for the coupled system, which handles its own residual and Jacobian assembly for the block system but still delegates linear solves to PETSc KSP objects configured via the same petsc_options dictionary (lines 62-65).

For pure phase-field sub-problems within these multiphysics simulations, the simple NewtonSolver wrapper is used directly. This modular design allows users to mix monolithic and staggered solution strategies while maintaining consistent linear solver configuration across the codebase.

Summary

• PhaseFieldX implements the Newton solver as a thin wrapper around dolfinx's PETSc-based NewtonSolver in src/phasefieldx/solvers/newton.py.
• Default convergence settings include rtol=1e-8, atol=1e-9, and max_it=500, using the residual norm as the convergence criterion.
• Linear solves within each Newton iteration use configurable Krylov (KSP) solvers from PETSc, accessible via solver.krylov_solver and customizable through PETSc options.
• The execution flow involves creating a NewtonSolverNonlinearProblem, wrapping it with NewtonSolver, and calling solver.solve() to perform the Newton-Krylov loop.
• For multiphysics problems, PhaseFieldX integrates with blocked Newton solvers while maintaining the same PETSc configuration interface for linear solves.

Frequently Asked Questions

How does PhaseFieldX handle convergence failures in the Newton solver?

If the Newton solver fails to converge within the maximum iteration count (default 500) or diverges, the solve() method returns converged=False along with the iteration count. The calling solver in PhaseFieldX typically checks this flag and can raise an exception or adjust the load step size in quasi-static simulations. Users can inspect the residual history via solver.krylov_solver.getConvergenceHistory() or enable report=True to see iteration-by-iteration diagnostics printed to stdout.

Can I change the linear solver from the default CG to GMRES?

Yes. The Krylov solver (KSP) is fully configurable via PETSc options. You can either set options programmatically through the petsc4py.PETSc.Options() object accessible via solver.opts, or pass command-line arguments when running your script. For example, to use GMRES with a Jacobi preconditioner, you would pass -ksp_type gmres -pc_type jacobi to the PETSc options, or set solver.opts[f"{solver.option_prefix}ksp_type"] = "gmres" before calling solver.ksp.setFromOptions().

What is the difference between the standard NewtonSolver and the blocked Newton solver in PhaseFieldX?

The standard NewtonSolver in src/phasefieldx/solvers/newton.py is designed for single-field nonlinear problems (e.g., a pure phase-field equation). It wraps dolfinx's PETSc Newton solver directly. The blocked Newton solver (typically scifem.BlockedNewtonSolver used in solver_ener_variational.py) is designed for coupled multi-field problems where multiple unknowns (displacement, phase-field, Lagrange multipliers) are solved simultaneously in a monolithic system. While the blocked solver handles its own residual and Jacobian assembly for the block matrix structure, it still delegates the linear solves within each Newton iteration to PETSc KSP objects configured via the same options interface.

How do I adjust the relaxation parameter for line search in the Newton solver?

The relaxation parameter (damping factor) is exposed via solver.solver.relaxation_parameter and defaults to 1.0 (full Newton step). To enable damping, set this value between 0 and 1 before calling solve(). For example, solver.solver.relaxation_parameter = 0.5 applies half of the computed Newton update at each iteration. Note that this applies a fixed damping factor; for adaptive line search algorithms (e.g., backtracking), you would need to configure PETSc's SNES line search options directly via the PETSc options database using prefixes like -snes_linesearch_type.

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