Performance Considerations for Large-Scale Phase-Field Simulations

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 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, 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 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) 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:

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

These settings, located in 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:

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 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 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 write human-readable logs only from rank 0, and certain diagnostics (maximum displacement, reaction forces) in 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:

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:

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 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.

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.

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