[DR] Incorrect calculation of boundary fluxes in GK solver

[manauref]

In the GK solver we compute boundary fluxes for at least two purposes:

1. Using them in the Boltzmann electron model, which obtains the sheath potential from ambipolarity at the sheath entrance: Gamma_e = Gamma_i (here Gamma_s is the particle flux of species s at the sheath entrance).
2. Computing fluxes due to recycling of ions, where the neutral flux into our domain is a function of the ion flux out of our domain.

At least for purpose 1 we were using the presently boundary flux, which is really the z surface term of the GK scaled, e.g. in 1x2v:

\mathrm{bflux}_k  = \frac{2}{\Delta z} \cdot \int_{-1}^{1} \mathrm{d}v_\parallel\,\mathrm{d}\mu~\left(\psi_k \dot{z}\cdot J\cdot B\cdot f\right)_{z=\pm 1}

where:

• $\Delta z$: cell length along z.
• $\psi_k$: phase basis function.
• $\dot{z}$: advection speed in the z direction.
• $J$: conf-space jacobian
• $B$: magnetic field amplitude
• $f$: distribution Later, to compute Gamma, integrate bflux over vpar-mu space, with the appropriate dimensional factors.

The above formula is incorrect if we want the physical particle flux through a surface (or in 1D, through a point). See Dhaeseleer's "Flux coordinates and magnetic field structure", equation 2.5.49a. We are missing a factor of the magnitude of the contravariant unit vector in z. The correct formula should be

\mathrm{bflux}_k  = \frac{2}{\Delta z} \cdot \int_{-1}^{1} \mathrm{d}v_\parallel\,\mathrm{d}\mu~\left(\psi_k \dot{z}\cdot J\cdot |\nabla z|\cdot B\cdot f\right)_{z=\pm 1}

One of the consequences of this error, which @Maxwell-Rosen has observed, is that in mirror sims with nonuniform z and Boltzmann electrons the sheath potential comes out different, compared to uniform z sims, and the difference is on the order of $\ln(\vert\nabla z\vert{\mathrm{uniform}}/\vert\nabla z{\mathrm{nonuniform}}\vert )$ .

Because of this extra factor of $|\nabla z|$ we can't use the collisionless GK boundary surf kernels, and we'll instead have to create new kernels for these calculations. I propose to change the bflux app and updater so that it uses such new kernels.

[JunoRavin]

Isn't grad(z) just the Jacobian arising from the non-uniform coordinates? Why isn't the Jacobian there the total Jacobian (Jacobian from field line following coordinates * Jacobian from further coordinate transformation)?

It seems like if we were missing this grad(z), the flux would be wrong at every interface in a non-uniform grid simulation.

[manauref]

grad{z} is not the Jacobian from the nonuniform coordinates, but it's related to it. The J in the above equations contains the Jacobian of the nonuniform coordinates, but that's not enough. We need to incorporate grad{z}.

The "flux" here is the physical particle flux, which we need for the Boltzmann electron model. The flux at interior surfaces is a computational flux, and is not quite the same and doesn't require this grad{z} factor, principally because our weak form is an inner product in computational coordinates, not in physical coordinates.