[DR] Pseudo-orbit average (POA) multiscale algorithm for magnetic mirrors

[Maxwell-Rosen]

Design Review

Overview

Description: This DR aims to describe the implementation of the pseudo-orbit average (POA) algorithm. The idea is to have a full dynamics phase (FDP) and an orbit-averaged phase (OAP). We alternate between the two phases to reach steady state.

Motivation: The magnetic mirror simulations are slow and we want to make them faster. We have explored a few mutli-scale algorithms and this one works well to capture the dynamics of the expander and the center of the mirror.

Impact: Many files will be affected in this change.

• A new timestepper will be added to the gyrokinetic/apps. This time stepper will be a superset of the SSPRK3, where this will pass parameters to the SSPRK3 stepper to toggle on and off certain parameters.
• Consideration must be taken when implementing the time stepper that we will aim to do a multi-species telescopic integraation scheme, freezing certain species at certain times.
• An updater is added to gyrokinetic/zero for determining the loss cone mask
• gk_species.c is modified to multiply the collisionless terms by a scalar number $\alpha$, so that the equation evolved is like $\partial_t f = \alpha \left\lbrace H, ~ f\right\rbrace + C + S + D$
• An app is added to multiply the whole GK equation by a multiplicative function $\mathcal{M}(x,v)$.
• An app is added to add a damping term to the RHS of the GK equation: $\partial_t f = \dots - \nu(x,v) f$.

Proposed Design

Key Design Elements: To the user, this should present as an enum option in the first layer of the input struct `.time_stepper = { .id = GK_PSEUDO_ORBIT_AVERAGE'. An option for SSP_RK3 should be added as well.
This should be in a struct where we can specify a number $\alpha$. In theory, $\alpha$ should be allowed to be different for each species in the simulation (thinking ahead for multi-species integration). The structure to specify the mask to apply during the OAP should be in this struct and depends on the species.

The way frames are output needs to be modified for this time stepper. It makes sense that the user would want N frames per phase, so each OAP gives N frames, and each FDP gives N frames. This goes against what we are doing in the input file for frame generation, so we must embed the input file loop inside the time stepper. From a software engineering perspective, this is very tricky, and I'm not sure what the appropriate design is for this.

Interaction with Existing Code: SSPRK3 will be the default option and it should be the default. Existing code should work out-of-the-box.

Prototyping Efforts (if applicable): Branch gk_dam_scale_hamil-prototyping

Sample Code:

  struct gkyl_gk app_inp = { 
    .time_stepper = {
      .id = GK_PSEUDO_ORBIT_AVERAGE,
      .frames_per_phase = 5,
      .scale_factor[2] = {alpha_ion, alpha_elc}, // Ions are index 0, electrons are index 1
      .damping[0] = { // The gk_damping_loss cone should distinguish where the trapped/passing boundaries are based on the charge of the species.
        .type = GKYL_GK_DAMPING_LOSS_CONE,
      }
      .damping[1] = { 
        .type = GKYL_GK_DAMPING_LOSS_CONE,
      }
    },
}

Design Review Checklist

Approval Criteria:

• Is the design coherent and consistent with Gkeyll's existing architecture?
• Does the design align with user-facing API standards?
• Do functionalities overlap with existing features?
• Are the algorithms robust and appropriate for the intended purpose?

Additional Notes

More work that needs to be done before this branch is merged into main:

1. The location of the mirror throat is hardcoded right here https://github.com/ammarhakim/gkeyll/blob/gk_damp_scale_hamil-merge/gyrokinetic/apps/gk_species_fdot_multiplier.c#L165. To fix this, we can utilize code from the position map, perhaps moving it into a general method inside calc_bmag, to calculate the mirror throat automatically. This would be a straightforward task; however, we want this to generalize to 2x so that each field line has its own independent Bmax and location. I'm not sure what this structure should look like. Perhaps we find the extrema of the magnetic field at the nodes, then interpolate between.

2. We must identify the loss boundary with more accuracy. Scaling with a factor of 2 is not good enough. We can identify the "loss" cells through the corners of the cells, i.e., the nodal values, rather than the cell-centered coordinate. This may improve accuracy.

3. Calculate the loss boundary for electrons in the expander region. There is a confined region of electrons here that traces a circle in phase space. Following the loss criterion arguments, you arrive at an equation

$$ \mu_{max} = \frac{m_s}{2 B(x)} \left( \frac{2 e \phi} {m_s} - v_{||}^2 \right) $$

This equation comes from thinking about

$$ m v^2 /2 < e \phi, $$

then separate the velocity into parallel and perpendicular, then solve for $\mu$.

4. Put all this work into a complete time-stepper app to complement SSP-RK3. At the input file level, I think a user should just specify $\alpha$ and identify the time stepper ID. This new time stepper will modulate between the orbit-averaged phase and the full dynamics phase. I would expect an input like "frames_per_phase = f" where it prints f frames per OAP or FDP phase. This contradicts the current loop structure in the input files, so I'm not sure how to implement this.
5. For beams and localized sources, we must compute the orbit-averaged source for the OAP. If we do not do this, the results will be terrible. We must average the source by integrating along orbit contours and taking the average flux of the source. We must compute

$$\mathring{S}(\vec x,\vec v) = \frac{\oint S(\vec\ell_x,\vec\ell_v)d\vec{\ell}}{\oint d\vec{\ell}}. $$

The orbit contours are

$$ d\ell_z = \frac{\partial \mathcal{H}(z,v_\parallel)}{\partial v_\parallel},\\ d\ell_{v_\parallel} = - \frac{\partial \mathcal{H}(z,v_\parallel)}{\partial z}. $$

Greg Hammett suggested an alternative. We could transform the source into $(E, \mu)$ coordinates and the orbit average will be much easier to perform.

Links

Additional Resources: Drafts of the paper that implements this scheme in some toy models.

Application_of_a_multiscale_pseudo_orbit_averaging_algorithm_for_partial_differential_equations_to_gyrokinetic_simulations_of_magnetic_mirrors.pdf
JCP__Part_1__A_multiscale_pseudo_orbit_averaging_algorithm_for_partial_differential_equations (1).pdf

[Maxwell-Rosen]

@manauref and I discussed that the input file structure should not look like this. Due to the interlinking with the POA scheme and the IO triggers, the POA scheme should be a modification to the while loop that is used in the input file. The while loop will toggle certain features of the code and determine whether to evolve the orbit-averaged equations or the full dynamics equations.

[Maxwell-Rosen]

Implemented in #934