[DR] Refactor of SR Vlasov for energy conservation and improved robustness of correction algorithm

[JunoRavin]

This design review encompasses the two key improvements contained in the branch: https://github.com/ammarhakim/gkylzero/tree/dg_10m_sr_refactor

1. In the attached pdf I outline the proof but the essential idea is that before, the SR Vlasov solve was not energy conserving. The original solver was formulated to separately compute gamma = sqrt(1 + p^2) (the particle Lorentz boost factor) and p/gamma (the relativistic velocity) independently. However, our work on Hamiltonian systems demonstrates that these two quantities are not independent and to properly have an energy conserving scheme one must

• First construct gamma to be continuous in the DG basis (using eval_on_nodes)
• Anywhere p/gamma is needed (streaming term, Lorentz force, and the current density moment for coupling to Maxwell's equations), this expansion must be computed as grad_p gamma.

I have verified that with this construction our scheme is now energy conserving using the relativistic two-stream instability, and we obtain third order accuracy in the convergence of the energy due to the SSP RK3 method employed for time-stepping. See attached figure from the 1x1v regression test.

2. We were previously having difficulty with the overall robustness of the LTE moment calculation and correction algorithm in relativity, principally due to the fact that in regions of strong gradients the drift velocity may be super-luminal (even if the cell average is less than the speed of light). Constraining the drift velocity to be well-behaved was even harder when trying to compute |V_drift|^2 for computing the Lorentz boost factor and the rest-frame density. This difficulty can be clearly seen in the main branch at the time of this issue (07/09/2024) where the Vlasov SR BGK Sod Shock test fails for Nx > 128 due to the initial large flow generated by the pressure gradient.

To solve this issue, I have changed the relativistic solver to utilize the spatial component of the bulk four-velocity in the LTE moment calculator and correction algorithm. The spatial component of the bulk four-velocity is guaranteed to be well behaved as it does not have to be bounded from above, and computing the needed Lorentz boost factor from this quantity is also guaranteed to be positive definite. To achieve this refactor the steps are now:

• Compute V_drift as M1i/M0 using weak division (as before).
• Check if (1 - |V_drift|^2) > 0 at Gauss-Legendre quadrature points; if it is, utilized the basis_sqrt function which evaluates the square root at Gauss-Legendre points and projects onto the modal basis to construct the basis expansion of 1/GammaV = sqrt(1 - |V_drift|^2)
• If (1 - |V_drift|^2) < 0 at any quadrature point instead evaluate V_drift at p=1 Gauss-Lobatto points (the vertices of a cell) and compute (1 - |V_drift|^2) at all the vertices of the cell.
• If (1 - |V_drift|^2) < 0 at a Gauss-Lobatto point, floor the value to 1.0e-16 (Lorentz boost factor of 10^8, a reasonable upper limit for our applications).
• Reconstruct the p=1 modal representation from these Gauss-Lobatto points using the nodal to modal transformation.
• Utilize the basis_sqrt function which evaluates the square root at Gauss-Legendre points and projects onto the modal basis, but now returns p=1 basis expansion of 1/GammaV = sqrt(1 - |V_drift|^2). Note: because we are using p=1 and evaluating at the vertices, this projection is guaranteed to be positive everywhere in the interior of the cell.
• Find the rest-frame density n using weak multiplication n = (1/GammaV) M0.
• Obtain the relativistic bulk four-velocity with weak division: u_i = M1i/n.

• Utilize this relativistic bulk four-velocity to obtain needed quantities for the pressure computation and LTE distribution, e.g., GammaV = sqrt(1 + |u_i|^2). Note that the pressure velocity moment is now modified to

  gamma GammaV^2 - 2 GammaV (u_i . p) + 1/gamma ((u_i . p)^2 - 1)

where u_i is the four-velocity GammaV*V_drift_i, p is the particle velocity, and gamma = sqrt(1 + p^2) is the particle Lorentz boost factor.

I have also rearranged things such that the SR moment calculator does not take a significant number of auxiliary variables anymore, only gamma (since you only ever need gamma or you compute the derivative of gamma inside the kernel) and instead moved the pressure computation to the "sr_vars" updater (https://github.com/ammarhakim/gkylzero/blob/dg_10m_sr_refactor/zero/dg_calc_sr_vars.c). This "sr_vars" updater is principally used by the LTE Moments routine to get the rest-frame density and rest-frame pressure, but it also provides the necessary routines for getting the initial continuous projections of gamma and 1/gamma.

This new formulation now can run the Vlasov SR BGK Sod Shock at higher resolution and has further generalized the LTE moments to be computable in 2x and 3x simulations.

SR_Vlasov_Notes.pdf

[JonathanGorard]

I'm recording this comment here regarding a conversation @JunoRavin and I just had: I think it would be interesting to see (if only as an error quantification scheme, e.g. in a unit test or something) what the timelike components of the bulk four-velocity look like, when computed from the spatial components that are now being used in the LTE moment and correction routines. As far as I can tell from an initial glance, that seems like the primary thing that might potentially go wrong (even though I understand that it's already far more stable than it used to be) with the algorithm as described.

[JunoRavin]

Design review merged with https://github.com/ammarhakim/gkylzero/pull/429