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

> Explore the internal workings of the Newton solver in PhaseFieldX. Learn how it leverages PETSc for efficient nonlinear phase-field problem solutions, configuring tolerances, sub-solvers, and damping.

- Repository: [Miguel Castillón/phasefieldx](https://github.com/castillonmiguel/phasefieldx)
- Tags: internals
- Published: 2026-02-26

---

**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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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:

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

```

This construction (line 30 in [`newton.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/newton.py):

```python
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`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Element/Phase_Field_Fracture/solver/solver_ener_variational.py) (lines 78-86):

```python

# 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`](https://github.com/castillonmiguel/phasefieldx/blob/main/newton.py)) writes the current Newton and KSP settings to a standard Python logger:

```python
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:

```python
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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`](https://github.com/castillonmiguel/phasefieldx/blob/main/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`.