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.

[manauref]

Looking at the "Inclusion of Sheath Wall Potentials", why do you have to introduce augmented escape barriers (EBs), instead of defining the left and right EBs as

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

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

where

$$U_L(z,\mu)= \mu B(z)+q(\phi(z)-\phi(z_L)), $$

$$U_R(z,\mu)= \mu B(z)+q(\phi(z)-\phi(z_R)), $$

? If this is valid it seems to me like it simplifies the physical motivation/logic.

[Maxwell-Rosen]

I'm not sure I understand. It seems like the equations you have in the comment are just a rephrasing of the ones in the PR description. If anything, I think my formulation of this is clearer because I explicitly use the term $\phi^{bc}$ for the boundary condition, so it's a bit unclear from your formulation where the boundary condition enters.

I asked Codex about this, and it says the sheath potential modifications aren't currently implemented. I think it's misunderstanding, though, since $\phi^{bc}$ is zero for the sheath boundary conditions.

[Maxwell-Rosen]

I meant to have the latest, most powerful models inspect this PR. It found that the wall potential wasn't included in this implementation (although it's zero for the sheath BCs), so I had it implement that to be more comprehensive and physical

[manauref]

Detailed comments in the review. Overall this seems good, thanks for doing it Max. My only overall comments are:

1. It's clear that LLMs wrote this and several of its choices seem odd. I think the code would've been improved (and reviewer's time saved) with a more careful read/thought prior to PR submission.
2. Please explain how correctness is determined or understood. What is each unit test checking and how was the correct value obtained?
3. I think the equations/logic for escape barriers can be made clearer, as in my earlier comment, but we don't have to act on that now.

[manauref]

why are you editing the gk_species_damping module? I thought you weren't using this nor needed to modify it?

[Maxwell-Rosen]

I'm editing the damping module because there is code in there that utilizes the loss cone updater. Nobody is using this feature, so I deprecated it.

[manauref]

should this comparison use the derivative normalized to some sensible value instead? perhaps something like

dbmag_vals[i] / ( B at this location / dz ) < 1e-10

typically comparing a physical value (dB/dz in this case) to a single number tends to be less robust.

Same for all the dbmag comparisons below.

[Maxwell-Rosen]

It's just checking if Bmag is increasing or decreasing here. The difference is taken elsewhere

[manauref]

how is this reg test different than rt_gk_mirror_boltz_elc_poa_1x2v_p1.c‎? if it has kinetic elc, then name it so.

and are all these regression tests running in under 90 s?

[Maxwell-Rosen]

Deleted. I think this was doing ion-electron collisions, which is why it was a _2species

[manauref]

don't see why this is needed

[Maxwell-Rosen]

Not needed. deleted

[manauref]

don't think this is needed
if necessary see gkyl_const.h

[Maxwell-Rosen]

I agree. It might be doing this because there is a longstanding linting issue (or its just on my laptop) where it doesn't recognize M_PI.

[manauref]

note that this hardcodes a particular ordering to the nodes, which is not desirable. It's best to use ranges where that order is abstracted out and can more easily be changed

[Maxwell-Rosen]

Updated

[manauref]

add some prefix, like gk_lcm_, to all these static private functions so it's clear where they are and that they don't come from other headers

[manauref]

and add comments describing them, and their inputs when they are not absolutely obvious (cdim is an example of obvious)

[Maxwell-Rosen]

Done

[manauref]

don't know why this is needed. We've never had this before

[Maxwell-Rosen]

Removed

[manauref]

i'd prefer that the use of extern C follows the pattern used before and in other .cu files

[Maxwell-Rosen]

Changed