# How to Run Parallel Simulations Efficiently with PhaseFieldX

> Learn to run parallel simulations efficiently with PhaseFieldX. Discover how PhaseFieldX automatically handles domain decomposition, parallel assembly, and distributed I/O using its FEniCSx/PETSc backend.

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

---

**Launch your script with `mpirun -n <N> python script.py` and pass `mpi4py.MPI.COMM_WORLD` to DOLFINx mesh constructors; PhaseFieldX automatically handles domain decomposition, parallel assembly, and distributed I/O through its FEniCSx/PETSc backend.**

PhaseFieldX is an open-source finite element framework built on top of **FEniCSx**, designed for phase-field modeling. Running parallel simulations efficiently requires minimal code changes because the library delegates all distributed memory operations to its underlying high-performance computing stack. This guide explains the MPI workflow, optimization strategies, and key implementation details found in the `castillonmiguel/phasefieldx` repository.

## How Parallel Execution Works in PhaseFieldX

PhaseFieldX leverages **MPI** (Message Passing Interface) through **mpi4py** to enable distributed memory parallelism. The architecture consists of three tightly coupled layers that operate transparently once you provide the communicator to mesh creation functions.

### Domain Decomposition

When you create a mesh using `dolfinx.mesh.create_*` functions and pass `mpi4py.MPI.COMM_WORLD`, FEniCSx automatically partitions the domain across available ranks. The partitioner uses **ParMETIS** or **SCOTCH** to balance element counts while minimizing inter-process edge cuts. Each rank owns a distinct subset of elements and maintains ghost cells for boundary communication.

In [`examples/PhaseField/plot_2003.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/examples/PhaseField/plot_2003.py), the mesh constructor receives the global communicator:

```python
msh = dolfinx.mesh.create_rectangle(
    mpi4py.MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([lx, ly])],
    [divx, divy],
    cell_type=dolfinx.mesh.CellType.quadrilateral,
)

```

### Parallel Linear and Nonlinear Solvers

PhaseFieldX uses **PETSc** via the `dolfinx.nls.petsc.NewtonSolver` class defined in [`src/phasefieldx/solvers/newton.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/solvers/newton.py). The solver inherits the MPI communicator from the mesh and distributes the linear system automatically across ranks. You can tune solver options through the `PetscOptions` object attached to the KSP (Krylov Subspace) object, selecting methods like **GMRES** or **CG** with preconditioners such as **GAMG** (Generalized Algebraic Multigrid).

### Distributed I/O and Post-Processing

Each MPI rank writes its local solution data to separate `.vtu` files (VTK unstructured grid). Rank 0 generates a master `.pvtu` (parallel VTU) file that references all partitions, enabling seamless visualization in ParaView. The `AllResults` class and post-processing utilities in [`src/phasefieldx/PostProcessing/ReferenceResult.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/PostProcessing/ReferenceResult.py) handle this automatically based on flags in your `Input` configuration.

## Step-by-Step Parallel Workflow

Follow these steps to execute scalable simulations:

1. **Launch with MPI** – Start the program using `mpirun -n <N> python <script>.py`, where `<N>` specifies the number of processes.

2. **Create Distributed Mesh** – Pass `mpi4py.MPI.COMM_WORLD` to DOLFINx mesh constructors. The library handles partitioning via ParMETIS/SCOTCH.

3. **Parallel Assembly** – Local element contributions are assembled on each rank. Ghost cell values synchronize through hidden MPI `Allreduce` and `Allgather` operations inside DOLFINx.

4. **Solve** – The Newton solver builds a PETSc KSP object on the same communicator. PETSc performs parallel preconditioning and Krylov iterations.

5. **Write Output** – After convergence, each rank writes its slice using XDMF/VTU APIs. Rank 0 writes the lightweight `.pvtu` master file for visualization.

## Code Example: Running a Parallel Phase-Field Simulation

The following minimal script demonstrates a static linear phase-field problem running in parallel. It requires no explicit MPI calls beyond passing the communicator to the mesh constructor.

```python

# parallel_demo.py

import numpy as np, dolfinx, mpi4py, os
from phasefieldx.Element.Phase_Field.Input import Input
from phasefieldx.Element.Phase_Field.solver.solver import solve
from phasefieldx.Boundary.boundary_conditions import bc_phi, get_ds_bound_from_marker

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

# 1️⃣  Mesh – distributed across all MPI ranks

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

divx, divy = 80, 40
lx, ly = 1.0, 0.5
msh = dolfinx.mesh.create_rectangle(
    mpi4py.MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([lx, ly])],
    [divx, divy],
    cell_type=dolfinx.mesh.CellType.quadrilateral,
)

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

# 2️⃣  Boundary identification

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

def bottom(x):
    return np.logical_and(np.isclose(x[1], 0), np.less(x[0], 0.5))

fdim = msh.topology.dim - 1
bottom_facets = dolfinx.mesh.locate_entities_boundary(msh, fdim, bottom)
ds_bottom = get_ds_bound_from_marker(bottom_facets, msh, fdim)

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

# 3️⃣  Function space & Dirichlet bc

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

V_phi = dolfinx.fem.functionspace(msh, ("Lagrange", 1))
bc = bc_phi(bottom_facets, V_phi, fdim, value=1.0)

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

# 4️⃣  Simulation data container

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

Data = Input(l=0.25, save_solution_vtu=True, results_folder_name="parallel_demo")

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

# 5️⃣  Call the generic Phase‑Field solver

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

solve(
    Data,
    msh,
    final_time=1.0,
    V_phi=V_phi,
    bcs_list_phi=[bc],
    update_boundary_conditions=None,
    update_loading=None,
    ds_list=[[ds_bottom, "bottom"]],
    dt=1.0,
    quadrature_degree=2,
)

```

Execute on eight processes:

```bash
mpirun -n 8 python parallel_demo.py

```

After completion, `parallel_demo/paraview-solutions_vtu/` contains:
- `phasefieldx_p0_000000.vtu` through `phasefieldx_p7_000000.vtu` (one per rank)
- `phasefieldx000000.pvtu` (master file for ParaView)

## Optimizing Performance for Large-Scale Runs

Efficient parallel execution requires balancing computational load and minimizing communication overhead. Configure these parameters based on your cluster architecture and problem size.

### Mesh Resolution per Rank

Maintain a moderate number of elements per MPI process. Too few elements leads to communication-dominated execution where MPI overhead exceeds computation. Adjust `divx` and `divy` parameters in your script to ensure each rank holds sufficient work.

### PETSc Solver Configuration

Select appropriate Krylov methods and preconditioners in [`src/phasefieldx/solvers/newton.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/solvers/newton.py). For large 2D and 3D phase-field problems, configure:

```python
-ksp_type gmres -pc_type gamg

```

These options enable Generalized Minimal Residual iteration with Algebraic Multigrid preconditioning, which scales well for elliptic phase-field operators. Set these via command line arguments or modify `self.opts` in the Newton solver class.

### Performance Profiling

Enable PETSc’s built-in logging to identify bottlenecks:

```bash
mpirun -n 8 python parallel_demo.py -log_view

```

This prints detailed timing breakdowns per rank, distinguishing between assembly and solver phases.

### I/O Frequency Reduction

Minimize filesystem traffic by writing VTU files only at selected timesteps. Control this behavior through the `save_solution_vtu` flag in the `Input` data class ([`src/phasefieldx/Element/Phase_Field/Input.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Element/Phase_Field/Input.py)).

## Key Source Files for Parallel Implementation

Understanding these files helps customize parallel behavior:

- **[`examples/PhaseField/plot_2003.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/examples/PhaseField/plot_2003.py)** – Full-featured MPI example showing mesh creation, boundary setup, and domain decomposition visualization.

- **[`src/phasefieldx/solvers/newton.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/solvers/newton.py)** – Wrapper around PETSc Newton solver exposing KSP/PC options and convergence tolerances.

- **[`src/phasefieldx/Element/Phase_Field/solver/solver.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Element/Phase_Field/solver/solver.py)** – Generic solver building residuals and invoking `NewtonSolver`.

- **[`src/phasefieldx/Boundary/boundary_conditions.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/Boundary/boundary_conditions.py)** – Utilities for Dirichlet/Neumann conditions on distributed meshes.

- **[`src/phasefieldx/PostProcessing/ReferenceResult.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/PostProcessing/ReferenceResult.py)** – Handles per-rank VTU writing and master PVTU file generation.

## Summary

- **PhaseFieldX parallelizes via FEniCSx/PETSc** – You only need to pass `mpi4py.MPI.COMM_WORLD` to mesh constructors.
- **Launch with `mpirun`** – No code changes required between serial and parallel execution.
- **Automatic domain decomposition** – ParMETIS/SCOTCH partition meshes; ghost cells handle communication.
- **Tune PETSc options** – Use GMRES + GAMG for large problems; set via command line or [`newton.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/newton.py).
- **Distributed I/O** – Each rank writes `.vtu` files; open the `.pvtu` master file in ParaView.

## Frequently Asked Questions

### Do I need to modify my code to run in parallel?

No. PhaseFieldX scripts remain identical between serial and parallel execution. The only requirement is passing `mpi4py.MPI.COMM_WORLD` to DOLFINx mesh creation functions. The underlying FEniCSx library handles partitioning, assembly, and solver distribution automatically.

### How do I visualize parallel output?

Open the `.pvtu` (Parallel VTU) file in ParaView. This lightweight XML file references all individual `.vtu` pieces written by each MPI rank. ParaView automatically assembles the distributed data for visualization without manual file merging.

### What is the optimal number of MPI ranks?

The ideal count depends on problem size and cluster architecture. Aim for at least several hundred elements per rank to maintain compute/communication balance. Use PETSc's `-log_view` flag to diagnose performance degradation from excessive communication overhead when scaling beyond optimal rank counts.

### Which PETSc preconditioner works best for phase-field problems?

For large 2D and 3D phase-field simulations, **GAMG** (Generalized Algebraic Multigrid) provides excellent scalability when paired with **GMRES** or **CG** Krylov solvers. Configure these via command line (`-ksp_type gmres -pc_type gamg`) or by modifying the `PetscOptions` object in [`src/phasefieldx/solvers/newton.py`](https://github.com/castillonmiguel/phasefieldx/blob/main/src/phasefieldx/solvers/newton.py).