# Performance Considerations for Large-Scale Phase-Field Simulations

> Discover performance considerations for large-scale phase-field simulations. PhaseFieldX achieves massive scalability using MPI-parallel execution and advanced solvers on distributed clusters.

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

---

**PhaseFieldX leverages MPI-parallel execution, blocked Newton solvers, and configurable PETSc KSP options to scale phase-field fracture simulations to tens of millions of degrees of freedom on distributed clusters.**

PhaseFieldX is built on **DolfinX** and **PETSc**, designed from the ground up for high-performance computing environments. Understanding the specific performance considerations for large-scale simulations requires examining how the library handles domain decomposition, linear algebra configurations, and parallel I/O strategies across its core solver modules.

## MPI-Parallel Architecture and Domain Decomposition

The foundation of PhaseFieldX's scalability lies in its strict adherence to MPI communicator patterns throughout the codebase. Every critical object—from meshes to function spaces—receives `msh.comm` and maintains data locality across ranks.

### Communicator Distribution

All parallel operations in [`solver_ener_variational.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/solver_ener_variational.py) utilize the mesh communicator to ensure each MPI rank works exclusively on its local sub-domain. When you create a mesh using `dolfinx.mesh.create_mesh`, DolfinX automatically partitions cells across available ranks. Subsequently, function spaces (`V_u`, `V_Φ`, `V_λ`) are constructed on these local sub-meshes, ensuring memory per rank scales proportionally with local degrees of freedom rather than global problem size.

### Scalable Global Reductions

Global quantities such as energies, norms, and reaction forces aggregate via `comm.allreduce` operations. In [`src/phasefieldx/Element/Phase_Field_Fracture/energy.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Element/Phase_Field_Fracture/energy.py), the `calculate_crack_surface_energy` and `compute_total_energies` functions perform local assembly followed by a single reduction call. This pattern avoids gathering large temporary arrays on rank 0, preventing out-of-memory crashes when meshes exceed millions of DOFs. Similarly, [`src/phasefieldx/norms.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/norms.py) implements parallel norm evaluation using logarithmic-cost reductions—O(log P) communication where P is the number of ranks.

## Solver Optimization Strategies

The default linear and nonlinear solver configurations in PhaseFieldX prioritize robustness, but specific tuning is required for large-scale performance.

### Blocked Newton Solver Design

The `RealSpaceNewtonSolver` class (found 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)) groups the displacement, phase-field, and Lagrange multiplier into a single block system. This blocked approach reduces the number of global solves required per load increment compared to segregated solvers. The residual and Jacobian assemblies utilize `dolfinx.fem.petsc.assemble_vector_block` and `dolfinx.fem.petsc.assemble_matrix_block`, avoiding monolithic matrix construction that would increase memory footprint and preconditioner complexity.

### PETSc KSP Configuration

By default, PhaseFieldX configures PETSc with direct LU factorization via MUMPS:

```python
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps"
}

```

These settings, located in [`solver_ener_variational.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/solver_ener_variational.py), provide robust convergence for moderate problem sizes but exhibit **O(N³)** memory growth. For simulations exceeding several million DOFs, switch to iterative Krylov methods:

```python
custom_options = {
    "ksp_type": "gmres",
    "ksp_max_it": 1000,
    "pc_type": "gamg",
    "pc_gamg_aggressive_coarsening": "1",
    "pc_gamg_type": "agg",
    "pc_gamg_threshold": "0.02"
}

pfx.Element.Phase_Field_Fracture.solver.solver_ener_variational.solve(
    ...,
    petsc_options=custom_options,
)

```

This GMRES + GAMG configuration scales roughly linearly in both memory and compute time, though it requires monitoring iteration counts via `ksp_max_it`.

### Adaptive Pseudo-Time Stepping

The adaptive load-stepping mechanism in [`solver_ener_variational.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/solver_ener_variational.py) dynamically adjusts the pseudo-time increment `dtau` based on Newton iteration counts. After each successful solve, the code increases `dtau` if convergence required two or fewer iterations, or decreases it using a golden-ratio factor of approximately 0.618 otherwise. This prevents excessive failed Newton steps and minimizes total solver calls for slowly evolving crack propagation.

## Memory Management and I/O Scalability

Efficient data output and logging are critical bottlenecks in large-scale distributed simulations.

### Parallel Output Formats

PhaseFieldX supports two distinct I/O strategies controlled via the `Data` configuration object:

- **XDMF with HDF5**: Activated when `Data.save_solution_xdmf = True`, this writes a single parallel file per simulation using the HDF5 format. This approach scales effectively on high-performance filesystems (Lustre, GPFS) and minimizes file descriptor exhaustion.
- **VTU (per-rank)**: Activated when `Data.save_solution_vtu = True`, this generates one file per MPI rank. While useful for post-processing in ParaView, this can overwhelm the filesystem when running on thousands of ranks.

Directory creation in [`src/phasefieldx/files.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/files.py) occurs exclusively on rank 0, followed by `comm.Barrier()` synchronization, eliminating race conditions during parallel startup.

### Memory-Light Post-Processing

Diagnostic functions throughout the codebase avoid memory aggregation. For instance, `Logger` utilities in [`src/phasefieldx/Logger/library_versions.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Logger/library_versions.py) write human-readable logs only from rank 0, and certain diagnostics (maximum displacement, reaction forces) in [`solver_ener_variational.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/solver_ener_variational.py) (lines 406–428) skip communication entirely when `comm.Get_size() == 1`, saving valuable bandwidth during large runs.

## Practical Implementation for Large Meshes

### Basic Large-Scale Driver Setup

The following pattern demonstrates proper initialization for a 100×100×100 hexahedral mesh:

```python
import dolfinx.mesh
import phasefieldx as pfx
import numpy as np

# Parallel mesh creation

mesh = dolfinx.mesh.create_box(
    dolfinx.mpi.MPI.COMM_WORLD,
    [np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0])],
    [100, 100, 100],
    cell_type=dolfinx.mesh.CellType.hexahedron,
)

# Function spaces on local sub-domains

V_u = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 2, (mesh.topology.dim,)))
V_phi = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))

# Data container with XDMF output enabled

class SimData:
    def __init__(self):
        self.Gc = 2.7e-3
        self.l = 0.01
        self.save_solution_xdmf = True  # Scalable parallel output

        self.save_solution_vtu = False  # Avoid per-rank file overhead

        self.results_folder_name = "large_sim"

data = SimData()

# Solver invocation

pfx.Element.Phase_Field_Fracture.solver.solver_ener_variational.solve(
    Data=data,
    msh=mesh,
    final_gamma=0.05,
    V_u=V_u,
    V_Φ=V_phi,
    dtau=5e-4,
    # petsc_options can be customized here for iterative solvers

)

```

### Profiling MPI Performance

Enable PETSc's built-in logging to identify load imbalances:

```bash
export PETSC_OPTIONS="-log_view"
mpirun -np 64 python driver.py

```

This outputs per-process assembly and solve times, revealing if specific ranks consistently lag in `dolfinx.fem.petsc.assemble_vector_block` operations—a symptom of mesh partition imbalance.

## Summary

- **MPI Distribution**: PhaseFieldX uses `msh.comm` throughout [`solver_ener_variational.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/solver_ener_variational.py) and energy modules to maintain data locality and perform logarithmic-cost global reductions.
- **Solver Configuration**: The default MUMPS direct solver provides robustness for moderate problems, but large-scale simulations require switching to `gmres` + `gamg` iterative methods to avoid **O(N³)** memory growth.
- **Block Assembly**: The `RealSpaceNewtonSolver` utilizes blocked vector and matrix assembly to reduce system size and improve preconditioner efficiency compared to monolithic approaches.
- **Adaptive Loading**: The golden-ratio based `dtau` adjustment minimizes Newton iteration counts during quasi-static crack evolution.
- **Scalable I/O**: XDMF/HDF5 output (`Data.save_solution_xdmf`) is strongly preferred over VTU for simulations running on hundreds or thousands of MPI ranks to prevent filesystem bottlenecks.

## Frequently Asked Questions

### How does PhaseFieldX handle mesh partitioning across MPI ranks?

When you invoke `dolfinx.mesh.create_mesh` with `dolfinx.mpi.MPI.COMM_WORLD`, DolfinX automatically distributes cells across ranks usingParMETIS or a similar partitioner. Function spaces created via `dolfinx.fem.FunctionSpace` operate exclusively on these local sub-meshes, ensuring that memory consumption scales linearly with local DOF count rather than global problem size. Global quantities are computed via `comm.allreduce` operations as demonstrated in [`src/phasefieldx/Element/Phase_Field_Fracture/energy.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Element/Phase_Field_Fracture/energy.py).

### When should I switch from direct to iterative solvers?

You should switch from the default `mumps` direct solver to iterative Krylov methods (e.g., `gmres` with `gamg` preconditioning) when your simulation exceeds approximately one million total degrees of freedom or when you observe memory exhaustion during factorization. Direct solvers exhibit **O(N³)** memory complexity, making them impractical for tens of millions of DOFs, whereas properly configured iterative solvers scale linearly in memory.

### What is the difference between XDMF and VTU output for large simulations?

XDMF format (`Data.save_solution_xdmf = True`) writes a single parallel HDF5 file accessible to all ranks, minimizing filesystem metadata operations and scaling to thousands of MPI processes. VTU format (`Data.save_solution_vtu = True`) generates one file per rank, which simplifies post-processing for small runs but can overwhelm parallel filesystems and exceed inode limits when using more than a few hundred ranks. For production large-scale runs, XDMF is strongly recommended.

### How does the adaptive dtau mechanism improve performance?

The adaptive pseudo-time step mechanism monitors Newton iteration counts after each load increment. If convergence requires two or fewer iterations, `dtau` increases to accelerate progress; otherwise, it decreases by a factor of approximately 0.618 to maintain stability. This prevents costly Newton solve failures in regions of rapid crack propagation while maximizing step size during stable elastic deformation, thereby reducing the total number of linear system solves required to reach the final load factor.