GK: loss cone updater refactored for arbitrary geometry

[Maxwell-Rosen]

Hamiltonian-based Classification of Loss and Trapped Orbits in Mirror Geometry

This PR was generated with the assistance of GitHub Copilot for the unit tests and the updater functions. Most of it is LLM-generated, but I've read it all and understand it.

Consider a one-dimensional coordinate $z\in[z_L,z_R]$ aligned with a magnetic field line, with magnetic-field magnitude $B(z)$ and electrostatic potential $\phi(z)$. For a species of mass $m_s$ and charge $q_s$, the guiding-center Hamiltonian in $(z,v_\parallel,\mu)$ coordinates is

$$ H(z,v_\parallel,\mu)=\frac{1}{2}m_s v_\parallel^2+\mu B(z)+q_s\phi(z), $$

where $v_\parallel$ is the parallel velocity and $\mu$ is the magnetic moment. In the absence of collisions, $H$ and $\mu$ are invariants of motion, and particle trajectories lie along contours of constant $H$ in $(z, v_\parallel)$.
Turning points satisfy $v_\parallel=0$, hence

$$H=U(z,\mu), \qquad U(z,\mu)= \mu B(z)+q\phi(z), $$

where $U$ is the Yushmanov potential.

Prior work simply identified the trapped/passing boundary based on the
loss cone criterion

$$\mu = \frac{\frac{1}{2}m_s v_\parallel^2 + q_s\left(\phi(z) - \phi_m\right)}{B_m - B(z)}.$$

In that implementation, $B_m,\phi_m$ are determined at the peaks of $B$. However, the effective potential experienced by particles is $U$, which combines both $B$ and $\phi$. This can be essential for identifying the trapped region of phase space induced by angled neutral-beam injection of ions at $45^\circ$, where an off-axis density peak occurs. Furthermore, this implementation only uses a single loss criterion per region, when multiple can exist, for instance, in tandem mirrors. A Hamiltonian-based region-identification algorithm is more accurate and flexible because it directly addresses the Yushmanov potential.

Two-Wall Escape Criterion

One issue with a Hamiltonian-based identifier is that $H$ is non-unique. For instance, consider a symmetric simple mirror without an electric potential, just two peaks in B. The magnetic field has fourfold degeneracy at the half-maximum. The points outside the peaks in $B$ should be labeled passing, while those between peaks of $B$ should be trapped.

A trajectory at fixed $(H,\mu)$ can escape through the left wall at $z_L$ if it can traverse the interval $[z_L,z]$ without violating $H\ge U$. For a phase-space point $(z,v_\parallel,\mu)$, define the escape barrier to the left wall as

$$EB_L(z,\mu)=\max_{s\in[z_L,z]}U(s,\mu), $$

and the escape barrier to the right wall as

$$EB_R(z,\mu)=\max_{s\in[z,z_R]}U(s,\mu). $$

The minimum escape barrier required for escape through at least one wall is

$$EB(z,\mu)=\min\left(EB_L(z,\mu),EB_R(z,\mu)\right). $$

Defining the signed loss function

$$\mathcal{F}(z,v_\parallel,\mu)=H(z,v_\parallel,\mu)-EB(z,\mu). $$

$\mathcal{F}\ge 0$ means that there is no energy barrier between the wall and this point, indicating a passing orbit. If $\mathcal{F}<0$, then the escape barrier is at a higher energy than this point, so it is trapped. So the orbiting mask is

$$\chi_{\rm orbit}(z, v_\parallel,\mu) = \mathrm{bool}(\mathcal{F}(z,v_\parallel,\mu)<0)$$

Inclusion of Sheath Wall Potentials

The field $\phi(z)$ does not directly represent the sheath drop at the wall, but instead it is assumed that $\phi=0$ in Gkeyll's sheath boundary conditions. Let $\phi_L^{\mathrm{bc}}$ and $\phi_R^{\mathrm{bc}}$ denote prescribed wall potentials at $z_L$ and $z_R$, respectively, and let $\phi(z_L)$ and $\phi(z_R)$ be the values implied by the loaded field data near the walls. These are left general, but taken as zero in the simulations. The effective wall potentials are taken as

$$\begin{aligned} U_L^{\mathrm{wall}}(\mu) &=\max\left(\mu B(z_L)+q\phi(z_L),\mu B(z_L)+q\phi_L^{\mathrm{bc}}\right), \\ U_R^{\mathrm{wall}}(\mu) &=\max\left(\mu B(z_R)+q\phi(z_R),\mu B(z_R)+q\phi_R^{\mathrm{bc}}\right). \end{aligned}$$

The corresponding augmented path barriers are

$$\begin{aligned} \widetilde{EB}_L(z,\mu)&=\max\left(\mathcal{B}_L(z,\mu),U_L^{\mathrm{wall}}(\mu)\right), \\ \widetilde{EB}_R(z,\mu)&=\max\left(\mathcal{B}_R(z,\mu),U_R^{\mathrm{wall}}(\mu)\right), \end{aligned}$$

with an augmented escape barrier

$$\widetilde{EB}(z,\mu)=\min\left(\widetilde{EB}_L(z,\mu),\widetilde{EB}_R(z,\mu)\right). $$

Replacing $EB$ by $\widetilde{EB}$
yields a classifier that includes sheath reflection consistently.

Integration into the Gkeyll Workflow

Within Gkeyll, one first interpolates the fields to nodal values, computes $EB$, calculates the loss barrier, and then, if any corners of the cell are not orbiting, the cell is not evolved.

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: (x (yes), blank (no))

• 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

Here is a comparison in a beam simulation with the current loss mask and the old one.
This is an R=32 mirror case with Boltzmann electrons, showing the difference in the fdot multiplier. Red means these cells are now masked with the new updater, while blue means cells that were masked, but are not with this version. In other words, this is new_mask - old_mask. We see the population of small mu particles which are trapped due to the off-axis beam injection, which is what we expect.

Here is another view, this time in pyvista

Here is a simulation run with the new loss cone mask. I'm re-running a beam R=32 case, which was used for the paper Gyrokinetic equilibria of high temperature superconducting magnetic mirrors

This did decrease the time-step during the OAP to 3.4e-7 versus 2.6e-6 seconds. Perhaps a more coarse mesh in mu would help.