[ju1s04]

PR summary

This proposes optimising the PPM characteristic projection by reducing thread-local temporary storage and merging the left and right state projections. Previously, the wave integrals ($I_p$) and ($I_m$) were materialized as two NEIGN x NEIGN thread-local arrays.

Instead, this change caches the reconstructed parabolic edge values for each variable and computes the corresponding wave integrals at their point of use. The left and right characteristic projections are also combined into a single sweep over the waves, allowing wave-dependent quantities to be computed once and reused for both projections.

PR motivation

The default reconstructs a parabola for each variable ($n$), computes its left and right reference states, and evaluates the wave integrals $I_p$ and $I_m$ for every characteristic wave $ii$. Materialising in $I_p [7][7]$ and $I_m[7][7]$. Together they contain $98$ double-precision values and therefore require $784$ bytes of storage per thread. And is only later on the algorithm, that they are taken for each wave it projects Im and Ip to compute the right and left states respectively.

Instead of storing every $I_m[n][ii]$ and $I_p[n][ii]$, only the reconstructed parabolic edge values $s_m$ and $s_p$ are stored in the arrays $s_{mc}[7]$ and $s_{pc}[7]$, requiring 14 doubles, or 112 bytes per thread. Each pair of wave integrals is then recomputed as a scalar at its point of use within the characteristic projection by calling:

ppm_int_profile_single(smc[n], spc[n], qbar[n], lambda[ii], delta, Ipt, Imt);

This reduces temporary storage by a factor of $7$, while the number of function calls remains unchanged at $49$, one for each variable-wave combination, computing $I_m$ and $I_p$ only when needed.

The optimisation also merges the two projections, as shown in the picture below; instead of having two separate loops for looping waves $ii$ to compute the right and left states separately (lines $15-27$ and lines $30-42$), it is restructured such that in a single loop the right and left contributions are computed in a sweep (lines $11-27$ optimised).

This also merges the two projection loops. Instead of computing the right and left states in separate wave loops (lines $15-27$ and $30-42$), both contributions are accumulated in a single loop (lines $11-27$).

Within this loop, the wave-sign tests and slope weights, $t^-=(\lambda_{ii}\leq0)$, $t^+=(\lambda_{ii}\geq0)$, $w^-=\lambda_0-\lambda_{ii}$ and $w^+=\lambda_6-\lambda_{ii}$, are evaluated once per wave and reused for both states. Here, $\lambda_0$ and $\lambda_6$ are the minimum and maximum MHD characteristic speeds. The original accumulation order is preserved: variables are accumulated into $LdQ$ for each wave, followed by accumulation over waves for each projected state.

When evaluated using the 3D Orszag-Tang test on an NVIDIA A100 GPU, this single optimisation reduced the runtime by $40.9$ seconds compared to the original implementation. A consistent performance improvement was also observed for the MHD Rotor and MHD Blast Wave tests and is expected to generalise to other MHD tests.

PR checklist

• test suite needs to be run on this PR
• this PR will change answers in the test suite to more than roundoff level
• all newly-added functions have docstrings as per the coding conventions
• the CHANGES file has been updated, if appropriate
• if appropriate, this change is described in the docs