Parallel Decomposition and Load Balance

Parallel decomposition divides a numerical-relativity problem among processes, threads, or devices. It changes execution order, communication, memory placement, and sometimes floating-point reduction, so it can affect both time to solution and bit-level results. Performance evidence is meaningful only for a fully specified physical problem and accuracy target. Strong and weak scaling, compute and communication time, imbalance, and reproducibility must be reported together; a raw core count or peak throughput is not a scientific result.

Helpful background. Formulations of General Relativity supplies the evolved systems, Ordinary Differential Equations supplies time stepping, Computational Relativity supplies reproducibility standards, and Probability and Statistical Inference supplies uncertainty summaries for variable timings.

Decompose work and data

Structured grids are commonly split into blocks with ghost exchanges. Spectral-element and discontinuous-Galerkin codes distribute elements and send interface traces. Particle, ray, or Monte Carlo work may require dynamic queues. The partition should minimize surface-to-volume communication while respecting memory limits and keeping work balanced.

Record the global mesh, local blocks or elements, refinement hierarchy, ghost width, communication graph, process and thread topology, affinity, device mapping, and solver collectives. A task graph can overlap communication with computation and expose irregular parallelism; Kidder et al. 2017, §4, §§5.1–5.3, Figures 4–6 and §6 describes this design in a discontinuous-Galerkin relativity code.

Load is not proportional to cell count when refinement levels subcycle, primitive recovery iterates variably, or horizon and radiation tasks are localized. Weight partitions with measured costs, then report compute, communication, idle, and imbalance times separately. Repartitioning has a migration cost and can perturb cache or reduction order; include it in the total.

Strong and weak scaling

For strong scaling, hold the complete physical problem fixed. With time $T_P$ on $P$ resources, define

$$S_P=\frac{T_{P_0}}{T_P},\qquad E_P=\frac{P_0T_{P_0}}{PT_P},$$

relative to an explicitly reported baseline $P_0$, not an inferred serial run that could not hold the problem. Saturation occurs when communication, synchronization, and insufficient local work stop reducing $T_P$.

For weak scaling, increase global problem size so work per resource is approximately fixed. State exactly which dimensions, AMR blocks, samples, or modes grow. Weak scaling is not comparable if the physics changes from a smooth wave to a shock-rich case or if refinement activity per resource changes. Report multiple repetitions, median and spread, warm-up policy, queue and filesystem effects, and whether initialization and output are included.

Kalinani et al. 2025, §§3.9–3.11, §4 and §§5.1–5.4, Tables 3–5 provides a modern GPU-accelerated GRMHD example with explicit code tests and scaling. Its reported machine and workload are contextual evidence, not a portable speed factor for another architecture.

Specified wave and scaling record

Use one physical problem for accuracy, strong scaling, and weak scaling:

$$\partial_tu+\partial_xu=0, \qquad u_{\rm ex}(t,x)=\sin\!\bigl[2\pi(x-t)\bigr].$$

The base domain is $[0,1]^3$, all faces are periodic, and $u(0,x)=\sin(2\pi x)$. Use fourth-order centered differences, RK4, and CFL $0.2$; run $N^3$ grids with $N=64,128,256$ to $t=1$. The primary reference is the analytic solution above. Compare the volume-weighted relative $L_2$ error $E_2$ and the phase of the independently computed $k_x=2\pi$ Fourier coefficient. Accept the $256^3$ solution only if $E_2\leq10^{-5}$, $|\Delta\phi|\leq10^{-5}\,\mathrm{rad}$, and the last two refinement ratios are consistent with fourth-order convergence to within $20\%$.

For strong scaling set $P_0=1$ node, hold the complete $256^3$ problem and output cadence fixed, and run on $P=1,2,4,8$ otherwise identical nodes. For weak scaling, tile $1,2,4,8$ unit cubes along $x$, retain the same periodic wave and $128^3$ cells per cube, and assign one cube per node. Pin ranks, threads, devices, and affinity per node. Use one warm-up followed by seven measured repetitions and report medians and interquartile ranges. The eight-node strong-scaling efficiency must be at least $0.50$, and the eight-node weak-scaling efficiency $T_1/T_8$ at least $0.80$. More importantly, every layout must retain $E_2\leq10^{-5}$ and $|\Delta\phi|\leq10^{-5}\,\mathrm{rad}$; its field difference from the corresponding one-node decomposition, after restricting to the same physical cube, must be at most $10^{-10}$ in relative $L_2$ norm. The named machine and software environment are part of the record, so these efficiency limits are test criteria rather than portable claims.

The weak-scaling tile needs its own accuracy qualification. For the stated fourth-order derivative, a Fourier mode with $\theta=2\pi/N$ sees the modified wave number

$$k_h=\frac{8\sin\theta-\sin(2\theta)}{6h}, \qquad h=1/N.$$

Applying the RK4 polynomial $R(z)=1+z+z^2/2+z^3/6+z^4/24$ for $5N$ steps gives amplitude $R(-ik_h\Delta t)^{5N}$ at $t=1$. Its phase error is approximately $1.94\times10^{-5}$ at $N=64$ and $1.22\times10^{-6}$ at $N=128$. These are direct evaluations of the discrete Fourier symbol, not timings or a three-dimensional implementation test. They explain why a $64^3$ tile cannot meet the declared $10^{-5}$ phase target merely by adding nodes.

Observed discrepancy: not available—specification only.

Preserve the scientific solution

Before comparing performance, demonstrate that every decomposition meets the same error target. Run an analytic wave or stationary black hole across several process counts and layouts. Compare constraints, conserved quantities, horizon invariants, and waveform phase. Non-associative reductions can shift last bits; the shift should remain below the declared numerical uncertainty and must not change convergence or failure decisions.

Use a production-representative benchmark with pinned resolution and output cadence. Also report a communication-focused microbenchmark only as diagnosis. Vary block size and rank placement, and identify the resource count at which efficiency loss outweighs wall-time gain. The validity domain is the named machine, software stack, compiler, network, and workload.

Environment and failure criteria

Pin machine and partition, node and device model, interconnect, compiler, MPI and accelerator runtimes, affinity, process/thread/device layout, problem configuration, input checksum, output cadence, profiling method, and code revision. A clean interface could be scale-run --config binary-grid.toml --layouts layouts.yaml --output run/, emitting timings by category, hardware counters where available, scientific diagnostics, repetition statistics, and environment metadata. It specifies a rerun record rather than claiming universal performance.

Reject a scaling claim with no real baseline, changed resolution or physics, omitted I/O that dominates production, single noisy timings, or speed obtained by loosening tolerances. Failure modes include load weights measured on an unrepresentative phase, global reductions at every substep, oversubscribed threads, hidden device transfers, and a decomposition whose different reduction order pushes an observable beyond its error budget.

Diagnostic

A run is four times faster on eight times as many GPUs, but waveform phase differs from the baseline by twice the accepted error because the accelerated layout uses a different reduction and looser primitive-recovery tolerance. Is the speedup valid for the scientific claim?

Solution

No. The two runs do not meet the same numerical target, so their wall times are not comparable as time to the validated result. Restore identical tolerances, quantify reduction-order sensitivity, repeat timings, and require phase, constraints, and conserved quantities to agree within the declared uncertainty. Then report speedup and efficiency relative to the actual baseline.

References

• Kalinani, J. V., Ji, L., Ennoggi, L., et al. (2025). “AsterX: A new open-source GPU-accelerated GRMHD code for dynamical spacetimes.” Classical and Quantum Gravity 42, 015010. arXiv.
• Kidder, L. E., Field, S. E., Foucart, F., et al. (2017). “SpECTRE: A task-based discontinuous Galerkin code for relativistic astrophysics.” Journal of Computational Physics 335, 84–114. arXiv.