Design Review

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, therfore 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 neccisary 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.

$$ \frac{\partial f}{\partial t} =\rm{min}(1, \Delta t_{\rm CFL} / \Delta t_{\rm min}) [H, f ] + C + S $$

This slows down the advective Poisson bracket by the appropriate ratio when $\Delta t_{\mathrm{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.

Prototyping branch: https://github.com/ammarhakim/gkeyll/tree/slowdown-to-omega-cfl

Implementation Details

Key changes: To calculate $\rm{min}(1, \Delta t_{\rm CFL} / \Delta t_{\rm min})$, 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.

The array methods are determined within gk_species_collisionles.c inside gk_species_collisionless_fdot_scaling_enabled. We could create a seperate function for this method, but this works well.

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.
    },

Preliminary testing

Ammar from my dream wanted testing results to see if this is a good idea. 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,