GK sheath BC: adding a mu dependence to parallel cut velocity

[Antoinehoff]

In this PR, we implement the ability to apply a conducting sheath boundary condition (BC) with a $\mu$-dependent parallel velocity cut (see DR #957 for more details). This PR also includes and improves the sheath BC unit test presented in PR #956. The kernels are generated with the branch of this Gkylcas PR

Design choice

The main design decision is to encode the $\mu$ dependence in a non-dimensional factor $\alpha(\mu)$ such that,

$$v_{\parallel \text{cut}}^2(\mu) = \alpha(\mu) v_{\parallel \text{cut} 0}^2$$

where $v_{\parallel \text{cut} 0}$ is the currently implemented cutting parallel velocity. This allows to preserve the existing code structure, leaving the cut velocity evaluation at the kernel level.

Implementation details

The main changes are mainly centered around the sheath BC updater gkyl_bc_sheath_gyrokinetic. A new DG array is added to the updater structure,

  struct gkyl_array *alpha_mu; // factor for mu dependent vcut in sheath BC.

to store the $\alpha(\mu)$ factor as a 2D $(v_\parallel,\mu)$ GK hybrid DG array. The use of a 2D representation simplifies the kernel generation by making the nodal evaluation more straightforward.

The updater creator, gkyl_bc_sheath_gyrokinetic_new, sets alpha_mu = 1 by default, which yields the current state of the sheath BC implementation.
One can manipulate alpha_mu with the new function gkyl_bc_sheath_gyrokinetic_set_alpha_mu to copy an existing array into alpha_mu.

Demonstration

We run gyrokinetic/unit/ctest_bc_twistshift.c which sets up Maxwellian distribution functions in 1x2v, 2x2v, and 3x2v domains. This test also sets a $\mu$-dependent parallel cut velocity as

void
eval_func_alpha_mu(double t, const double *xn, double* GKYL_RESTRICT fout, void *ctx)
{
  double vpar = xn[0], mu = xn[1];

  double Lmu = 1.0; // Characteristic scale length in mu direction.
  double alpha_mu_0 = 1.0;

  double alpha_mu = alpha_mu_0 * (1 + exp(-mu/Lmu));
  fout[0] =  pow(alpha_mu, 2);
}

which adds an exponential decaying term to the cutting velocity. This function is projected on a DG GK hybrid 2v basis:

  gkyl_proj_on_basis *projalphamu = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
      .grid = &grid_vel->grid_vel,
      .basis = basis_vel,
      .num_ret_vals = 1,
      .eval = eval_func_alpha_mu,
      .ctx = NULL,
    }
  );
  gkyl_proj_on_basis_advance(projalphamu, 0.0, &grid_vel->local_vel, alpha_mu);
  gkyl_proj_on_basis_release(projalphamu);

where the grid_vel is the same as the one that define the bc_sheath updater and the basis_vel is obtained using a new bc_sheath interface

gkyl_bc_sheath_gyrokinetic_get_alpha_mu_basis(bcsheath, &basis_vel);

We then copy our alpha_mu array to the bc_sheath updater using

gkyl_bc_sheath_gyrokinetic_set_alpha_mu(bcsheath, alpha_mu);

We run and plot the resulting distribution functions using this script: bc_sheath_test.sh, which produces figures like:

Note that we varied the Maxwellian parallel flow depending on the dimensionality, specifically $u_\parallel=0$ for 1x2v (bottom row), $u_\parallel=v_{th}$ for 2x2v (middle row), and $u_\parallel=-v_{th}$ for 3x2v (top row).