13 Parallel Computing in CMPS
13.1 Process model
CMPS uses MPI domain decomposition with a dedicated driver/GUI rank and worker-owned numerical partitions. Rank 0 is the main process. The remaining ranks own the finite-volume cells, faces, solution arrays, distributed matrix/vector objects, and multigrid levels used by the solver.
This division has an important numerical consequence: rank 0 coordinates the calculation but does not assemble worker-owned CFD equations. Any algorithm added to the coupled solver must preserve the collective schedule expected by all workers.
13.2 Partitioning and mappings
The fine mesh is partitioned into worker-owned subdomains. CMPS maintains local/global cell and face mappings so that a worker can use compact local indices while the main process can still gather results or relate partition data to the original mesh.
For partition \(r\), write the local-to-global cell map as \[g=\mathcal{L}_r(i),\] with inverse map \[i=\mathcal{G}_r(g)\] for global cells owned or referenced by that partition. Similar mappings are maintained for faces.
Graph adjacency among partitions is derived from inter-partition connectivity. The graph communicator is then used for neighbour-oriented communication after partitioning.
13.3 Interface faces and halo data
A face that was an interior face on the undistributed mesh becomes a partition-interface face if its two cells are assigned to different workers. Each worker stores the local owner state and receives the neighbour state from the adjacent worker. Numerically, the face flux remains the same function \[\mathbf{F}_f=\mathbf{F}_f(\mathbf{q}_L,\mathbf{q}_R),\] but \(\mathbf{q}_R\) may be a received interface/halo value rather than a local owned-cell object.
Before an assembly that requires neighbour states or gradients, workers exchange the corresponding interface data. The ordering of state exchange, gradient exchange, residual assembly, and correction update must remain identical on all workers that participate in the communication graph.
13.4 Distributed coupled solve
Each worker executes the active coupled matrix driver on its local partition. Interior and boundary faces are assembled exactly as in serial execution, while interface data supply remote neighbour states. The distributed sparse matrix and vectors are then solved through the configured CMPS linear-solver path. CPU workers use the Intel MKL implementation, while GPU-enabled workers can use the NVIDIA CUDA implementation.
The global nonlinear progress is a collective property. If \(r_j^{(k)}\) denotes a local residual measure for worker \(k\), a global norm/reduction must be formed collectively, for example \[\lVert r_j\rVert_2 =\left(\sum_k\lVert r_j^{(k)}\rVert_2^2\right)^{1/2}.\] Rank 0 receives the reduced monitoring/status data and decides whether the synchronized solver loop continues, stops, or responds to a user command.
13.5 GPU acceleration
CMPS supports GPU-enabled execution of the scalable linear-algebra path in addition to CPU execution. The finite-volume domain decomposition and nonlinear assembly remain defined by the same MPI ownership model; GPU acceleration changes where supported matrix/vector and solver operations are executed, not the governing equations or the partition topology.
For a worker rank assigned to an NVIDIA GPU, the CMPS CUDA path keeps supported sparse linear-algebra operations on the device. The computationally intensive work of Krylov iteration and AMG—sparse matrix-vector products, block operations, vector updates, dot products/reductions, relaxation kernels, and multilevel transfer operations—can consequently execute on the GPU. Host/device transfers should be avoided inside the iteration because repeated movement of the coupled vectors can dominate the cost for otherwise fast CUDA kernels.
The block-coupled formulation is well suited to accelerator execution because each mesh adjacency carries a dense \(n_b\times n_b\) coupling block. The benefit nevertheless depends on problem size, block size, AMG hierarchy, device memory capacity, and how much of the selected preconditioner remains device-resident.
13.6 Multi-GPU and GPU-cluster execution
GPU-cluster execution combines the existing MPI decomposition with GPU-enabled worker ranks. Conceptually, for \(N_r\) worker ranks and \(N_g\) accelerators, each worker is bound to an available GPU according to the launch/runtime mapping. The CFD mesh remains partitioned by MPI exactly as in the CPU cluster path, \[\Omega=\bigcup_{r=1}^{N_r}\Omega_r,\] and each rank operates on its owned rows of the distributed coupled matrix.
The resulting execution has two distinct communication scales:
intra-rank device work: local sparse/block operations, vector kernels, and AMG smoothing/coarsening operations execute on the rank’s GPU when supported by the configured backend;
inter-rank communication: partition-interface values, Krylov reductions, and AMG coarse-level communication are exchanged through MPI across GPUs and across compute nodes.
Thus CMPS GPU-cluster support is not a separate CFD discretization. It is the distributed-memory solver architecture of Sections 13.1–13.4 with accelerator-backed local algebra. A multi-node run can therefore use one or more GPUs per node while retaining the same global block system and MPI partition interfaces.
For strong scaling, three costs become increasingly important: the surface-to-volume ratio of the CFD partitions, global reductions required by the selected Krylov method, and communication/coarse-grid work in AMG. AMG-preconditioned flexible Krylov methods are particularly relevant on large GPU clusters because local GPU throughput alone cannot remove global low-frequency error; a scalable coarse correction is required to prevent the Krylov iteration count from growing excessively with partition count.
13.7 GPU-aware correctness and performance requirements
A GPU-enabled distributed solve must preserve the same mathematical result as the CPU path to the tolerance of floating-point reduction ordering. In addition to the MPI requirements listed later in this chapter, the accelerator path should ensure that:
a worker uses the GPU assigned to that rank and does not accidentally contend with unrelated ranks for the same device unless oversubscription is intentional;
the selected CUDA matrix, vector, Krylov, and AMG kernels are compatible with the active block layout and precision;
halo exchanges and Krylov reductions operate on valid, synchronized data;
preconditioner setup and rebuild events do not leave stale device-side matrix data after a nonlinear Jacobian update;
fallback to a CPU-only operation is treated as a performance issue, not as a reason to change the residual or matrix definition; and
GPU memory use is budgeted for the fine matrix, Krylov basis vectors, and all active AMG levels.
13.8 Communication volume
For a partition, the dominant nearest-neighbour solution exchange scales approximately with \[V_{comm}\propto N_{if}\,n_b,\] where \(N_{if}\) is the number of interface faces/cells exchanged and \(n_b\) is the number of active variables. Consequently a good partition minimizes interface size while balancing owned-cell work. Adding species, turbulence, or particle equations increases both local block cost and interface payload.
This is why minimizing only the number of cells per partition is insufficient. Graph partition quality, surface-to-volume ratio, block size, linear solver communication, and load imbalance all influence parallel efficiency.
13.9 Multigrid ownership
The agglomeration/geometric multigrid hierarchy is worker-owned. Workers build and operate their local coarse levels while preserving the existing fine-grid partition and interface architecture. Rank 0 may gather diagnostics or visualization information, but it must not reconstruct an alternative coarse solver object that competes with the worker hierarchy.
Collective convergence/failure logic is particularly important in multigrid and FMG initialization: one worker cannot return early while others are still entering interface exchanges or reductions. A detected invalid state must be converted into a collective failure decision before the hierarchy is unwound.
13.10 Rank-0 GUI and data collection
The main process owns the GUI interaction and sends case/solver-option changes to the workers. For post-processing, solution data can be gathered from worker-local indexing to the global visualization representation. Because a full gather can become expensive for large three-dimensional problems, diagnostics should prefer reduced quantities or partition-aware visualization paths when the complete field is not required.
13.11 Correctness requirements for new parallel features
A new CMPS numerical feature is parallel-safe only if it satisfies all of the following:
worker ranks enter required MPI collectives in the same order;
rank 0 does not execute worker-owned matrix or mesh operations;
interface data are refreshed after state changes that invalidate halo values;
global convergence and failure decisions are collective;
local/global mappings remain valid after any hierarchy or case change;
rollback or cleanup does not delete fine-grid worker objects needed by the normal solver.
These conditions are part of the numerical algorithm, not merely implementation details, because violating them can change the residual seen by a partition or deadlock the distributed iteration.