OpenFOAM Parallel Computing — Domain Decomposition and Scaling

#+CATEGORY: openfoam # :PROPERTIES: # :ID: uuid-openfoam-parallel-computing # :END:

OpenFOAM Parallel Computing :: Domain Decomposition and HPC Execution

This note documents OpenFOAM's parallel execution framework: domain decomposition methods, parallel solver execution, communication overhead, and scaling behaviour. Parallel computing is essential for production CFD — the mesh sizes required for accurate simulations (10⁶–10⁸ cells) are simply too large for single-core execution.

The theoretical underpinning is the conservation equation (Conservation Laws) — each subdomain solves the same discretised equation, and the interface fluxes must be exactly consistent to preserve global conservation. If the decomposition creates an imbalance at domain interfaces, the global conservation error increases, which directly affects solution accuracy.

**Table of Contents**

• Parallel Execution Model :: MPI-based domain decomposition
• Decomposition Methods :: scotch, hierarchical, ptscotch, metis
• decomposeParDict :: Configuration
• Parallel Solver Execution :: mpirun and HPC cluster setup
• Communication Overhead :: Ghost cells, halo exchange
• Scaling Behaviour :: Amdahl's law and strong/weak scaling
• Reconstruction :: reconstructPar and field averaging
• See Also
• References

Parallel Execution Model :: MPI-Based Domain Decomposition

OpenFOAM uses Message Passing Interface (MPI) for parallel execution. The computational domain is decomposed into subdomains, each assigned to one MPI rank (process). Each rank owns a subset of cells and solves the discretised equations independently. At each iteration/rank communicates boundary fluxes and field data with neighbouring ranks through halo (ghost) regions.

OpenFOAM's decomposition happens before the simulation starts — the decomposePar utility partitions the mesh into subdomains and writes separate data for each rank. The solver then runs in parallel with all ranks decomposed. Fluent and Star-CCM+ perform decomposition dynamically during the simulation, which provides better load balancing but with higher overhead.

Decomposition Methods :: scotch, hierarchical, METIS

**scotch (Scalable Community Toolkit for Optimization of Graph and Sparse Matrices)**

The default decomposition method. Scotch uses multi-level graph partitioning: the mesh graph is coarsened (cells merged) into a smaller graph, partitioned on the coarse level, then mapped back to the fine mesh. This produces high-quality decompositions with minimal interface area (which minimises communication).

scotch is installed with most OpenFOAM distributions. The library is libscotch.so (serial) and libptscotch.so (parallel). Scotch typically produces better load balancing than hierarchical decomposition, especially for non-uniform meshes where cell sizes vary significantly.

**ptscotch (Parallel Scotch)**

The parallel version of scotch. Use ptscotch when decomposing large meshes on a machine with many cores (> 100 cores for the decomposePar step itself). ptscotch distributes the decomposition computation across multiple cores, significantly reducing the decomposePar runtime.

**hierarchical (recursive bisection)**

Hierarchical decomposition uses recursive coordinate bisection (RCB) or recursive inertial bisection (RIB). RCB bisects the domain along the longest axis; RIB uses the inertia tensor of the cell distribution to find the optimal bisecting plane. Hierarchical is faster than scotch (simpler algorithm) but produces lower-quality decompositions — more cells at domain interfaces, poorer load balancing for non-uniform meshes.

Use hierarchical for simple, uniform meshes (structured hex from blockMesh) where the domain is already well-balanced. For unstructured meshes (snappyHexMesh, cfMesh), scotch or METIS produce better load balancing.

**METIS (multi-level graph partitioning)**

METIS is the original multi-level graph partitioning library. Scotch is inspired by METIS. The two often produce similar-quality decompositions for engineering meshes. METIS is faster on some problem types (structured or nearly-structured meshes) and slower on others (highly anisotropic meshes). Use either scotch or METIS depending on which produces better load balancing for your specific mesh — compare the interface area (number of faces at domain boundaries) for both.

decomposeParDict :: Parallel Decomposition Configuration

The decomposeParDict file in system/ controls how the mesh is decomposed:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      decomposePar;
}

method            scotch;              // decomposition method

// scotch or ptscotch
scotch
{
    cutType           scotchMesh;
    nDims             1;               // or: 2, 3; 1 = automatic
    coeffs
    {
        distWeight    0;              // importance of load balancing (0 = balance, 1 = interface)
        commWeight    1;              // importance of communication
    }
}

// METIS
// metis { ... }

// hierarchical (RCB/RIB)
// hierarchical { nSubDomains 4; coefficients (1 1 1); }

// manual
// manual
// {
//     coefficients
//     (
//         (  // Rank 0
//             0 0 0
//             20 20 20
//         )
//         (  // Rank 1
//             20 0 0
//             40 20 20
//         )
//         // ... more ranks
//     );
// }

// Number of subdomains (ranks)
// Set automatically based on number of MPI processes
// nSubDomains is NOT in decomposeParDict — the number of ranks is inferred from the number of subdirectories created

subDomains
(
    4 4 4   // nSubDomains x y z = 64 total
);

For structured hexahedral meshes, nDims = 3 (3D decomposition) is standard. For thin domains (e.g., boundary layer meshes where the wall-normal direction has far fewer cells than streamwise/spanwise), use nDims = 2 or manual decomposition to avoid splitting thin cells across domains.

Parallel Solver Execution :: mpirun and HPC Cluster Setup

Parallel solvers are invoked with mpirun (or mpiexec):

[#BEGIN_SRC bash # Local parallel execution (8 ranks on this machine) mpirun -np 8 simpleFoam -parallel

srun -n 64 pimpleFoam -parallel

mpiexec -n 128 pimpleFoam -parallel

#+END_SRC

The -parallel flag tells OpenFOAM to run in parallel mode. Each MPI rank owns a subdomain and writes to its own subdirectory in the case directory (0.0/, 0.1/, ..., 0.(N-1)/ for rank 0, 1, ..., N-1). Field data at write time is split across ranks and must be reconstructed for post-processing (see Reconstruction below).

HPC cluster execution: most HPC clusters use Slurm, PBS/Torque, or LSF job schedulers. The OpenFOAM installation on the cluster must be configured with:

. MPI library (OpenMPI, MPICH, or Intel MPI) — compiled with the same compiler used to build OpenFOAM . Shared filesystem (for reading input data, writing output data) . Node-local scratch disk (for temporary files, if needed)

Parallel I/O is a common bottleneck in HPC clusters: if all ranks write to a shared filesystem simultaneously, I/O performance degrades. The recommended practice is to use a node-local scratch disk for intermediate writes, then copy the final results to the shared filesystem at the end of the simulation.

Communication Overhead :: Ghost Cells and Halo Exchange

The communication overhead is the single largest factor limiting parallel efficiency. Every iteration of the solver requires halo exchange — the transfer of field data across domain boundaries:

. Pressure equation: all ranks must exchange pressure values at domain boundaries before solving the pressure Poisson equation. The pressure solver (GAMG) is inherently parallel, but each AMG cycle requires halo exchange. . /Velocity/Momentum/: flux at domain boundaries must be consistent across ranks. When one rank computes a flux through a shared face, the neighbouring rank must compute the exact same flux (conservation across domain boundaries). . Turbulence quantities: k, epsilon, omega at domain boundaries must be communicated. If a turbulence model uses gradients across boundaries, halo exchange must include both the field values and their gradients.

The communication volume scales as O(N^{2/3}) for 3D meshes with N total cells and uniform decomposition — as you add more ranks, the surface-area-to-volume ratio of each subdomain increases, meaning more communication per computation. This is the reason strong scaling (fixed mesh, more ranks) eventually breaks down: at some point, each rank has so few cells that communication dominates computation.

For a mesh with N = 10⁷ cells decomposed onto P = 1000 cores: each core has ~10⁴ cells. The interface area (number of cells at domain boundaries) is approximately O(N^{2/3}) = O(10^{14/3}) ≈ 10⁴·⁷ ≈ 50,000 cells at interfaces. Of these ~50,000, about half are domain-boundary cells — so each core has roughly ~5000 interface cells per core, which must communicate with their neighbours at every iteration. For GAMG (pressure solver), this halo exchange happens once per iteration — for the momentum equation, it happens once per iteration. The total communication is proportional to the number of interface cells times the number of fields communicated (U, p, k, omega, epsilon, T, species fractions).

Scaling Behaviour :: Amdahl's Law and Strong/Weak Scaling

OpenFOAM's scaling behaviour follows standard MPI parallel computing principles:

**Strong scaling** (fixed mesh, more cores):

For a 10⁶-cell mesh:

The efficiency drops because communication overhead grows as O(N^{2/3}) while computation scales as O(N/P). At P = 1000, communication dominates.

**Weak scaling** (10⁴ cells/core):

For a 10⁴–10⁵ cells/core mesh, scaling is nearly perfect up to ~1000 cores, because the communication-to-computation ratio remains constant.

**Amdahl's Law**: The serial fraction of the solver limits speedup. For OpenFOAM, the serial fraction is typically 5–10% (I/O, load balancing, GAMG coarse-grid correction with sequential coarsening). The theoretical maximum speedup is 1/0.05 = 20× for strong scaling.

For production runs, the recommended strategy is weak scaling: use 10⁴–10⁵ cells per core, with 100–500 cores total. This gives the best parallel efficiency with minimal communication overhead.

Reconstruction and Post-Processing :: reconstructPar

After parallel execution, field data is written to per-rank subdirectories (0.0/, 0.1/, ..., 0.(N-1)/). The reconstructPar utility combines these into a single coherent time directory:

[#BEGIN_SRC bash reconstructPar # reconstruct all times reconstructPar -time 0 # reconstruct only time 0 reconstructPar -latestTime # reconstruct only the latest time

#+END_SRC

reconstructPar creates a new reconstruction/ directory with subdirectories for each time (0/, 1/, ..., 100/). The reconstructed data can be visualised directly (paraFoam -builtin reconstruction/0) or analysed for field statistics.

For parallel post-processing, OpenFOAM provides:

. fieldValues functionObjects — run each on their own rank, then reconstructPar combines the results . parallel sampling — sample each rank locally, then reconstructPar combines the samples . parallel VTK export — foamToVTK -parallel writes per-rank VTK files, which ParaView can read in parallel (reducing I/O overhead compared to reconstructPar + serial VTK export)

For field averaging in parallel: use the fieldAverage functionObject with base time — OpenFOAM accumulates time-averaged statistics on each rank in parallel, then reconstructPar produces the global average from per-rank averages (weighted by cell volume).

fieldAverage
{
    type fieldAverage;
    libs ("libfieldFunctionObjects.so");
    fields
    (
        U
        {
            mean on;
            prime2Mean on;
            base time;
        }
        p
        {
            mean on;
            base time;
        }
    );
}

The prime2Mean option computes Reynolds-stress-like quantities =\overline{u'_i u'_j} = \overline{u_i u_j} - \bar{u}_i \bar{u}_j/, which are the turbulence statistics needed for LES/RANS comparison.

See Also :: Related Notes

. Conservation Laws — conservation across domain boundaries requires exact flux matching . FVM Discretisation — face fluxes must be consistent across subdomain boundaries . OpenFOAM Solver Selection — parallel solvers (each solver must be compiled with -DUSE_PARALLEL) . OpenFOAM Post-Processing — parallel VTK export, parallel functionObjects . OpenFOAM Multiphase Flows — parallel VOF (MULES requires halo exchange for alpha) . ANSYS vs OpenFOAM Comparison — parallel execution comparison

References

. OpenFOAM User Guide. ESI OpenCFD. (Parallel execution, decomposePar, reconstrucPar.) . OpenFOAM Programming Guide. ESI OpenCFD. (Parallel API: reduce, sum, gSum.) . Scotch documentation. INRIA. (Graph partitioning algorithm details.) . Amdahl, G.M. (1967). "Validity of the single-processor approach to achieving large scale computing capabilities." AFIPS Conference Proceedings. (Amdahl's law.) . HPC cluster documentation. (Slurm, PBS/Torque, LSF job scheduling.)