Time dilation of collisionless {H,f} to control CFL

[Maxwell-Rosen]

Feature

Summary

Late on a Friday night, I had a research meeting with Ammar, Mana and Jimmy in my sleep. It was mostly with Ammar, but Jimmy and Mana were there. I explained to them the issues we had in G2 when we implemented force softening. In the 2023 paper, we published the results below.


The issue is that the force softening introduced higher temperatures in the expander. Expander heating is something we are incredibly concerned about right now because of the positivity hack, so this scheme is not ideal. I explained in my dream that the issue was that we are modifying the Hamiltonian, therefore modifying the orbit contours, and the orbit contours are essential for particle confinement. We really should not modify the Hamiltoninan for an appropriate scheme.

Instead, I was reminded of the work in this article (“Time-Dilation Methods for Extreme Multiscale Timestepping Problems” https://arxiv.org/html/2510.09756). They essentially discuss slowing down time near the extreme regions of the problem, effectively what we have done for the orbit average phase of the pseudo orbit averaging integrator. This method preserves the orbit contours, masks the CFL, and will speed up our simulations. Ammar was doubtful this would work and dismissed my idea, however I was very confident that this would work, so he wanted to see some evidence. I awoke to program in the middle of the night, teeming with excitement and motivation.

I realize another issue is that we used an overly harsh mask in the 2023 paper. It is not necessary to create a completely logical mask, going all the way from 1 to 0. Instead, the mask should optimally choose its strength based on the CFL frequency and rate. It should modify the problem the bare minimum in order to get the optimal CFL rate. That lead me to this modification of the GK vlasov equation.

Here is a derivation of this feature. It follows a similar derivation to Werner (2018) "Speeding Up Simulations By Slowing Down Particles: Speed-Limited Particle-In-Cell Simulation". Let’s start with a phase-space continuity equation of the form

$$ \frac{\partial f}{\partial t} = [\mathcal{H},f] + C(f) + S. $$

To dilate time, consider a slowed-down distribution function

$$ g(\vec x, \vec v, t) = \beta(\vec x, \vec v, t) f(\vec x, \vec v, t) $$

where $\beta(\vec x, \vec v, t)$ is the amount which time is dilated, $\beta \in (0,1]$. Simulations thus far in Gkeyll have been evolving $f(\vec x, \vec v, t)$, but we aim to evolve $g(\vec x, \vec v, t)$. Substituting $g/\beta$ into the time derivative of the phase-space continuity equation yields

$$ \frac{\partial(g/\beta)}{\partial t} = [\mathcal{H},f] + C(f) + S. $$

Applying the product rule, we get

$$ \frac{1}{\beta} \frac{\partial g}{\partial t} -\frac{g}{\beta^2} \frac{\partial \beta}{\partial t} = [\mathcal{H},f] + C(f) + S $$

Multiplying by $\beta$, this can be expressed in terms of the time derivative of $g$.

$$ \frac{\partial g}{\partial t} = \frac{g}{\beta} \frac{\partial \beta}{\partial t} + \beta \left([\mathcal{H},f] + C(f) + S \right). $$

This equation correctly slows down the advection, collisions, and source terms, which relaxes the time-step.
However, the resulting equation is not closed unless $\beta$ is specified independently, to which two options are presented.
One option is to approximate the time variation of $\beta$ as slow, setting $\partial_t\beta \approx 0$, or to treat $\beta$ as time-independent.
A second approach is to set $\beta$ as a state function, either of $g$ or the CFL limits imposed by equation \eqref{eq: Full vlasov equation toy models}. This means equation \eqref{eq: slowed down continuity equation with dbeta dt} is not determining the evolution of $\beta$ in time, but instead using $\partial_t \beta$ to compute, therefore an Euler step (or higher order RK methods) backwards in time can be made to compute $\partial_t \beta = \left( \beta_n - \beta_{n-1} \right) / \Delta t$.
The $\partial_t \beta$ term can be integrated exactly, avoiding additional timestep restrictions.
The time-derivative term involving $\beta$ is purely local in phase space and does not introduce a CFL constraint; the timestep remains controlled by the rescaled characteristic velocities.

The option which is implemented is to set

$$ \beta =\rm{min}(1, \Delta t_{\rm CFL} / \Delta t_{\Omega_H} ) $$

This slows down the advective Poisson bracket by the appropriate ratio when $ \Delta t_{ \rm CFL} < \Delta t_{ \Omega_H } $. This numerical approach to the mask is also more flexible, as the code or user does not need to know about the mirror loss cone, therefore it is a more general method of time-dilation in inaccessible regions of phase space.

Issue link: I will open a DR as well about this idea

Implementation Details

Key changes: To calculate $ \rm{min}( 1, \Delta t_{ \Omega_H } / \Delta t_{ \rm CFL}) $, some modifications are needed inside the array_ops and dg_bin_ops methods.

We need an inversion operation to invert a P0 array. This should use the dg_inv method, which is already written. However, there is no kernel for the P0 array. Division is tricky with DG arrays, so it should not be an array_op, even though this is a P0 array. I handwrote a P0 kernel for this method because the Maxima generators do not work for P0. It just does 1/f[0], which is a straightforward operation and I really don't want to make the operation 1/f[0] more complex then neccisary.

This is followed by a $\rm min()$ operation, so that also needed to be implemented. I decided on the location for this function inside array_ops for a few reasons. One, this is an operation on a single array. Second, it is a simple operation, not requiring a basis. Third, I really don't want to think about how this generalizes to a P1 array and implement a full kernel operation. I see the issue that this method is only intended for P0 arrays, however I think that is the purpose of these array_ops methods. It is for single array manipulation, independent of basis. I really don't want to make this method implementation more complex then neccisary and do P1,P2 kernels and write Maxima scripts when the operation is just `fmin( coeff, array)' element wise.

Dependencies:None

Automated testing: New unit tests are added to ctest_array_ops.c. A regression test is used to test the mask of the CFL.

Example Use

    .collisionless = {
      .type = GKYL_GK_COLLISIONLESS_ES,
      .cfl_dt_min_value = 1e-8, // User specified minimum dt
      .write_diagnostics = true,
      .cfl_dt_min_omegaH = true, // To limit the CFL based on the omegaH mode.
    },

Community Standards

• Documentation has been updated.
• My code follows the project's coding guidelines.
• Changes to layer/zero should have a unit test, e.g., core/zero.

Testing:

• I added a regression test to test this feature.
• I added this feature to an existing regression test.
• I added a unit test to test this feature.
• Ran make check and unit tests all pass.
• I ran the code with make valcheck, and it is clean.
• I ran the code through computer-sanitizer on GPU, and it is clean.
• I ran a few regression tests to ensure no apparent errors.
• Tested and works on CPU.
• Tested and works on multi-CPU.
• Tested and works on GPU.
• Tested and works on multi-GPU.

Additional Notes

I've been testing some simulations with this new feature. I will get some better results that are better direct comparsion. It speeds the FDP up by about 20x, which is huge. Preliminary results show no difference in the moments and a much better result thatn the skip_cell feature.

This figure compares two simulations, one with the time dilation feature, the other without it. This is a non-uniform z simulation with 216 cells using the positivity hack. The simulation executed 10 POA cycles, with a 20 μs FDP and 2 s OAP. The simulation ends around 20 s. The simulations compare nearly identically, except for minor deviations in Tpar in the mirror throat. The time dilation saves 20x during the FDP and the simulation goes from 22 hours to around 5-6 hours. It's a big cost saving,

I also wanted to compare a high resolution simulation Nz=800 without the positivity hack. The simulations executed 10 POA cycles with a 20 μs FDP and 2 s OAP. The results agree very well, although they needed to use a higher $\mu$ resolution of 32 cells instead of 16.

Here is the percent error for each quantity above, showing less than 1% deviation for all quantities with a 2x speedup.