[PR] FLR effects with Padé approximants

[Antoinehoff]

This PR implements the proposed algorithm of DR #797.

Implementation details

The code change resides in the implementation of a new superlu/cusparse updater to apply a LHS matrix. This is the solution found to perform the invert_flr step, i.e.,

$$\phi = (1 - \rho_i^2 \nabla_\perp^2) \Phi_0$$

where $\Phi_0$ is the solution of the quasi-neutrality equation with FLR corrected charge density (see #797 for more details).
This new updater, gkyl_fem_poisson_perp_lhs_apply has the advantage to use the exact same representation as the gkyl_fem_poisson_perp_solve, which should ensure coherence in the discretization. A unit test verify the capability of the updater by solving the problem

$$(1 - \rho^2 \nabla_\perp) \phi = w,$$

using gkyl_fem_poisson_perp_solve where the RHS is built as

$$w = (1 - \rho^2 \nabla_\perp) g,$$

using the new updater gkyl_fem_poisson_perp_apply_lhs on an analytical function $g$.

Interface

We add a new entry gyroradius into the gkyl_gyrokinetic_flr table and add an instance of this structure in the field table. The new regression test, rt_gk_flr_d3d_iwl_2x2v_p1.c is testing the FLR feature. For the species, we pass

    .flr = {
      .type = GKYL_GK_FLR_PADE_CONST,
      .Tperp = ctx.Te0,
    },

where .Tperp will be used to compute the gyroradius (the magnetic field amplitude can be specified too but here it is taken from the app table). For the field, we pass

    .flr = {
      .type = GKYL_GK_FLR_PADE_CONST,
      .gyroradius = ctx.rho_i,
    },

where .gyroardius is the ion gyroradius that will be used to get $\phi$ from $\Phi_0$.
The type GKYL_GK_FLR_PADE_CONST means that we use only one reference gyroradius. The former version of the code used a sum over the gyroradii which is not correct but is still available using GKYL_GK_FLR_PADE_CONST_SUM.

Caveats and open questions

This implementation should conserve particle and energy (to be checked). One important limitation is that we cannot use varying boundary conditions, nor varying gyroradii, if we want to maintain conservation. While the former is not very common, the latter may limit the physical accuracy of our simulations.

This PR reintroduce cuSPARSE, that I heard was removed from Gkeyll a while ago. I don't see an interference on CPU nor GPU (cuda13 on Perlmutter). I could not use cuDSS only because it does not have a routine to apply a matrix. One thing that is not very logical in the PR is that the cuSPARSE use is in the cuDSS file. Should we create a new cuSPARSE file ?

Additional feature

By mistake I implemented the apply LHS updater using the gkyl_fem_poisson in moments. This updater also has a unit test and seems to work. I don't know if we want to keep it or remove it, because it is not used.

Results

TCV 2x2v

We use the TCV 2x2v regression test and increase the resolution to 24x16x12x8 to compare quickly the effect of the FLR model. We see an effect of the FLR operator on the potential.

Postgkyl command

pgkyl gk-load-quantity --quantity field --name gkeyll --species ion --frame : --path /pscratch/sd/a/ah1032/gkeyll_main/flr/tcv_2x2v/wk interpolate --interp 2 select --z1 0.0 --comp 0 collect plot --figure 0 --yshift -0.04 --zmax 300.0 --zmin -300.0 --no-legend --xlabel '' --ylabel '' --clabel '$\phi$ (V)' --title 'TCV 2x2v FLR' --subplot-xlabels '$t$ (s)' --subplot-ylabels '$R-R_{lcfs}$ (m)'grid --colormap RdBu_r

One main difference is that the initial poisson solve is non zero when we use FLR, even if $n_i^g=n_e^g$ at $t=0$. This may be physical because with FLR effects activated, the densities $n_s^g=\int dv f_s$ are passed through the FLR operator before computing the charge density. Since the FLR operator is dependent on the gyroradius, the resulting gyroaverage densities change between electron and ions which will create a non zero charge density at $t=0$.

Postgkyl command

pgkyl gk-load-quantity --quantity field --name gkeyll --species ion --frame 0 --path wk dg-local-poly select --z1 0.0 plot --figure 0 --xshift -0.04 --no-legend --xlabel '' --ylabel '' --title 'TCV 2x2v FLR $t=0$' --subplot-xlabels '$R-R_{lcfs}$ (m)' --subplot-ylabels '$\phi$ (V)'grid