From 5932f2c4ec08871a662c3a1e6eebf0f545e2589f Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Mon, 11 May 2026 13:00:17 -0400
Subject: [PATCH 01/11] Add test with showing negative vtsquared in the ion lbo
 moments when using ion-electron collisions in 2x

---
 gyrokinetic/creg/rt_gk_mirror_failure.c | 481 ++++++++++++++++++++++++
 1 file changed, 481 insertions(+)
 create mode 100644 gyrokinetic/creg/rt_gk_mirror_failure.c

diff --git a/gyrokinetic/creg/rt_gk_mirror_failure.c b/gyrokinetic/creg/rt_gk_mirror_failure.c
new file mode 100644
index 0000000000..1de04ededc
--- /dev/null
+++ b/gyrokinetic/creg/rt_gk_mirror_failure.c
@@ -0,0 +1,481 @@
+#include <math.h>
+#include <stdio.h>
+#include <time.h>
+
+#include <gkyl_alloc.h>
+#include <gkyl_const.h>
+#include <gkyl_eqn_type.h>
+#include <gkyl_fem_poisson_bctype.h>
+#include <gkyl_gyrokinetic.h>
+#include <gkyl_gyrokinetic_run.h>
+#include <gkyl_math.h>
+
+#include <rt_arg_parse.h>
+
+// Define the context of the simulation. This is basically all the globals
+struct gk_mirror_ctx
+{
+  int cdim, vdim; // Dimensionality.
+  // Plasma parameters
+  double mi;
+  double qi;
+  double me;
+  double qe;
+  double Te0;
+  double n0;
+  double B_p; // Hardcoded magnetic field at midplane 0.53 for double lorentzian
+  double Bmag_midp; // Magnetic field magnitude at midplane (Z=0)
+  double beta;
+  double tau;
+  double Ti0;
+  double nuFrac;
+  // Ion-ion collision freq.
+  double logLambdaIon;
+  double nuIon;
+  double vti, vte;
+  double RatZeq0; // Radius of the field line at Z=0.
+  double kperp; // Perpendicular wavenumber for ES gyrokinetics.
+  // Axial coordinate Z extents. Endure that Z=0 is not on
+  double z_min;
+  double z_max;
+  double psi_eval;
+  double psi_max;
+  double psi_min;
+  double theta_eval;
+  double theta_min;
+  double theta_max;
+  // Physics parameters at mirror throat
+  double vpar_max_ion;
+  double mu_max_ion;
+  double vpar_max_elc;
+  double mu_max_elc;
+
+  int Npsi;
+  int Ntheta;
+  int Nz;
+  int Nvpar;
+  int Nmu;
+  int Nvpar_elc;
+  int Nmu_elc;
+  int cells[GKYL_MAX_DIM]; // Number of cells in all directions.
+  int poly_order;
+
+  // Source parameters
+  double ion_source_amplitude;
+  double ion_source_sigma;
+  double ion_source_temp;
+
+  double t_end; // End time.
+  int num_frames; // Number of output frames.
+  int num_phases; // Number of phases.
+  double write_phase_freq; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
+  double int_diag_calc_freq; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
+  double dt_failure_tol; // Minimum allowable fraction of initial time-step.
+  int num_failures_max; // Maximum allowable number of consecutive small time-steps.
+
+  // Geometry parameters for Lorentzian mirror
+  double mcB;     // Magnetic field parameter
+  double mcB_inner; // Magnetic field parameter for inner mirror
+  double mcB_outer; // Magnetic field parameter for outer mirror
+  double gamma;   // Width parameter for Lorentzian profile
+  double gamma_inner; // Width parameter for Lorentzian profile for inner mirror
+  double gamma_outer; // Width parameter for Lorentzian profile for outer mirror
+  double Z_m;     // Mirror throat location
+  double L_center; // Length of central region
+  double L_plugs; // Length of plugs
+  double Z_min;   // Minimum Z coordinate
+  double Z_max;   // Maximum Z coordinate
+  double psi_in;  // Working variable for psi integration
+  double z_in;    // Working variable for z integration
+};
+
+
+// Evaluate initial conditions
+// I think ICs should be constructed to match the expander rather than the center
+// as the center is quickly filled in the OAP
+void
+initial_density(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  fout[0] = 1e17;
+}
+
+void
+initial_upar(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  double z = xn[1];
+  double c_s = 7 * sqrt(app->Te0/app->mi);
+  if (fabs(z) <= app->Z_m)
+  {
+    fout[0] = 0.0;
+  }
+  else
+  {
+    fout[0] = fabs(z) / z * c_s * tanh(4 * (app->Z_max - app->Z_m) * fabs(fabs(z) - app->Z_m));
+  }
+}
+
+void
+eval_zero(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  fout[0] = 0.0;
+}
+
+void
+initial_temp_ion(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  fout[0] = app->Ti0/10.0;
+}
+
+void
+initial_temp_elc(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  fout[0] = app->Te0;
+}
+
+void mapc2p_vel_ion(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  double vpar_max_ion = app->vpar_max_ion;
+  double mu_max_ion = app->mu_max_ion;
+
+  double cvpar = vc[0], cmu = vc[1];
+  double b = 1.4;
+  vp[0] = vpar_max_ion*tan(cvpar*b)/tan(b);
+  vp[1] = mu_max_ion*pow(cmu,2);  // Cubic map in mu.
+}
+
+void mapc2p_vel_elc(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  double mu_max_elc = app->mu_max_elc;
+  double vpar_max_elc = app->vpar_max_elc;
+  double cvpar = vc[0], cmu = vc[1];
+  vp[0] = vpar_max_elc*cvpar;
+  vp[1] = mu_max_elc*pow(cmu,4);  // Cubic map in mu.
+}
+
+struct gk_mirror_ctx
+create_ctx(void)
+{
+  int cdim = 2, vdim = 2; // Dimensionality.
+  int poly_order = 1;
+
+  // Universal constant parameters.
+  double eps0 = GKYL_EPSILON0;
+  double mu0 = GKYL_MU0; // Not sure if this is right
+  double eV = GKYL_ELEMENTARY_CHARGE;
+  double mp = GKYL_PROTON_MASS; // ion mass
+  double me = GKYL_ELECTRON_MASS;
+  double qi = eV;  // ion charge
+  double qe = -eV; // electron charge
+
+  // Plasma parameters.
+  double mi = 2.014 * mp;
+  double Te0 = 940 * eV;
+  double n0 = 3e19;
+  double B_p = 0.53; // Bmag at z=0
+  double beta = 0.4;
+  double tau = pow(B_p, 2.) * beta / (2.0 * mu0 * n0 * Te0) - 1.;
+  double Ti0 = tau * Te0;
+
+  // Ion-ion collision freq.
+  double nuFrac = 1.0;
+  double logLambdaIon = 6.6 - 0.5 * log(n0 / 1e20) + 1.5 * log(Ti0 / eV);
+  double nuIon = nuFrac * logLambdaIon * pow(eV, 4.) * n0 /
+                 (12 * pow(M_PI, 3. / 2.) * pow(eps0, 2.) * sqrt(mi) * pow(Ti0, 3. / 2.));
+
+  // Thermal speeds.
+  double vti = sqrt(Ti0 / mi);
+  double vte = sqrt(Te0 / me);
+
+  // Grid parameters
+  double vpar_max_ion = 16 * vti;
+  double mu_max_ion = mi * pow(3. * vti, 2.) / (2. * B_p);
+  double vpar_max_elc = 4 * vte;
+  double mu_max_elc = me * pow(4. * vte, 2.) / (2. * B_p);
+
+  int Nz = 32;
+  int Npsi = 2;
+  int Nvpar = 16;
+  int Nmu = 8;
+  int Nvpar_elc = 8;
+  int Nmu_elc = 8;
+
+  // Geometry parameters.
+  double RatZeq0 = 0.10; // Radius of the field line at Z=0.
+  double Z_min = -2.5;
+  double Z_max =  2.5;
+  double mcB = 3.691260;
+  double gamma = 0.226381;
+  double Z_m = 0.98;
+
+  // Calculate phase structure
+  double write_phase_freq = 1; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
+  double int_diag_calc_freq = 100; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
+  double dt_failure_tol = 1.0e-4; // Minimum allowable fraction of initial time-step.
+  int num_failures_max = 20; // Maximum allowable number of consecutive small time-steps.
+
+  struct gk_mirror_ctx ctx = {
+    .cdim = cdim,
+    .vdim = vdim,
+    .poly_order = poly_order,
+
+    .mi = mi,
+    .qi = qi,
+    .me = me,
+    .qe = qe,
+    .Te0 = Te0,
+    .n0 = n0,
+    .B_p = B_p,
+    .beta = beta,
+    .tau = tau,
+    .Ti0 = Ti0,
+
+    .nuFrac = nuFrac,
+    .logLambdaIon = logLambdaIon,
+    .nuIon = nuIon,
+    .vti = vti,
+    .vte = vte,
+    .RatZeq0 = RatZeq0,
+    .vpar_max_ion = vpar_max_ion,
+    .mu_max_ion = mu_max_ion,
+    .vpar_max_elc = vpar_max_elc,
+    .mu_max_elc = mu_max_elc,
+
+    .Npsi = Npsi,
+    .Nz = Nz,
+    .Nvpar = Nvpar,
+    .Nmu = Nmu,
+    .Nvpar_elc = Nvpar_elc,
+    .Nmu_elc = Nmu_elc,
+    .cells = {Npsi, Nz, Nvpar, Nmu},
+    .Z_min = Z_min,
+    .Z_max = Z_max,
+
+    .write_phase_freq = write_phase_freq,
+    .int_diag_calc_freq = int_diag_calc_freq,
+    .dt_failure_tol = dt_failure_tol,
+    .num_failures_max = num_failures_max,
+
+    .mcB = mcB,
+    .gamma = gamma,
+    .Z_m = Z_m,
+  };
+  
+  // Populate a couple more values in the context.
+  ctx.psi_max = 3e-3;
+  ctx.psi_min = 1e-5;
+  ctx.psi_eval = (ctx.psi_max + ctx.psi_min) / 2.0;
+
+  ctx.z_min    = -2.5;
+  ctx.z_max    = 2.5;
+
+  return ctx;
+}
+
+int main(int argc, char **argv)
+{
+  struct gkyl_app_args app_args = parse_app_args(argc, argv);
+
+#ifdef GKYL_HAVE_MPI
+  if (app_args.use_mpi) MPI_Init(&argc, &argv);
+#endif
+
+  if (app_args.trace_mem) {
+    gkyl_cu_dev_mem_debug_set(true);
+    gkyl_mem_debug_set(true);
+  }
+
+  struct gk_mirror_ctx ctx = create_ctx(); // Context for init functions.
+
+  int rank = 0;
+#ifdef GKYL_HAVE_MPI
+  if (app_args.use_mpi) MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+#endif
+  // if (rank == 0)
+  //   print_ctx(&ctx);
+
+  int cells_x[ctx.cdim], cells_v[ctx.vdim];
+  for (int d=0; d<ctx.cdim; d++)
+    cells_x[d] = APP_ARGS_CHOOSE(app_args.xcells[d], ctx.cells[d]);
+  for (int d=0; d<ctx.vdim; d++)
+    cells_v[d] = APP_ARGS_CHOOSE(app_args.vcells[d], ctx.cells[ctx.cdim+d]);
+
+  // Construct communicator for use in app.
+  struct gkyl_comm *comm = gkyl_gyrokinetic_comms_new(app_args.use_mpi, app_args.use_gpu, stderr);
+
+  struct gkyl_gyrokinetic_species ion = {
+    .name = "ion",
+    .charge = ctx.qi,
+    .mass = ctx.mi,
+    .vdim = ctx.vdim,
+    .lower = {-1.0, 0.0},
+    .upper = { 1.0, 1.0},
+    .cells = { cells_v[0], cells_v[1]},
+    .polarization_density = ctx.n0,
+
+    .projection = {
+      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
+      .density = initial_density,
+      .ctx_density = &ctx,
+      .upar = initial_upar,
+      .ctx_upar = &ctx,
+      .temp = initial_temp_ion,
+      .ctx_temp = &ctx,
+    },
+
+    .mapc2p = {
+      .mapping = mapc2p_vel_ion,
+      .ctx = &ctx,
+    },
+
+    .collisionless = {
+      .type = GKYL_GK_COLLISIONLESS_ES,
+      .write_diagnostics = true,
+    },
+
+    .collisions = {
+      .collision_id = GKYL_LBO_COLLISIONS,
+      .den_ref = ctx.n0,
+      .temp_ref = ctx.Ti0,
+      .num_cross_collisions = 1,
+      .collide_with = { "elc" },
+      .write_diagnostics = true,
+    },
+
+    .write_omega_cfl = true,
+    .num_diag_moments = 8,
+    .diag_moments = {GKYL_F_MOMENT_BIMAXWELLIAN, GKYL_F_MOMENT_M0, GKYL_F_MOMENT_M1, GKYL_F_MOMENT_M2, GKYL_F_MOMENT_M2PAR, GKYL_F_MOMENT_M2PERP, GKYL_F_MOMENT_M3PAR, GKYL_F_MOMENT_M3PERP },
+    .num_integrated_diag_moments = 1,
+    .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP },
+    .time_rate_diagnostics = true,
+
+    .boundary_flux_diagnostics = {
+      .num_integrated_diag_moments = 1,
+      .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP},
+    },
+  };
+
+
+  struct gkyl_gyrokinetic_species elc = {
+    .name = "elc",
+    .charge = ctx.qe,
+    .mass = ctx.me,
+    .vdim = ctx.vdim,
+    .lower = {-1.0, 0.0},
+    .upper = { 1.0, 1.0},
+    .cells = { ctx.Nvpar_elc, ctx.Nmu_elc },
+
+    .polarization_density = ctx.n0,
+
+    .mapc2p = {
+      .mapping = mapc2p_vel_elc,
+      .ctx = &ctx,
+    },
+
+    .projection = {
+      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
+      .density = initial_density,
+      .ctx_density = &ctx,
+      .upar = eval_zero,
+      .ctx_upar = &ctx,
+      .temp = initial_temp_elc,
+      .ctx_temp = &ctx,
+      .correct_all_moms = true,
+    },
+
+    .collisions =  {
+      .collision_id = GKYL_LBO_COLLISIONS,
+      .den_ref = ctx.n0,
+      .temp_ref = ctx.Te0,
+      .num_cross_collisions = 1,
+      .collide_with = { "ion" },
+      .write_diagnostics = true,
+      .not_in_dfdt = true,
+    },
+
+    .num_diag_moments = 1,
+    .diag_moments = {GKYL_F_MOMENT_MAXWELLIAN},
+  };
+
+  struct gkyl_gyrokinetic_field field = {
+    .gkfield_id = GKYL_GK_FIELD_BOLTZMANN,
+    .electron_mass = ctx.me,
+    .electron_charge = ctx.qe,
+    .electron_temp = ctx.Te0,
+    .is_static = false,
+  };
+
+  struct gkyl_mirror_geo_grid_inp grid_inp = {
+    .filename_psi = "gyrokinetic/data/eqdsk/lorentzian_R10.geqdsk_psi.gkyl", // psi file to use
+    .rclose = 0.2, // closest R to region of interest
+    .zmin = -2.5,  // Z of lower boundary
+    .zmax =  2.5,  // Z of upper boundary
+    .include_axis = false, // Include R=0 axis in grid
+    .fl_coord = GKYL_GEOMETRY_MIRROR_GRID_GEN_SQRT_PSI_CART_Z, // coordinate system for psi grid
+  };
+
+  struct gkyl_gk app_inp = {  // GK app
+    .name = "gk_lorentzian_mirror",
+    .cdim = ctx.cdim,
+    .upper = {ctx.psi_max, ctx.Z_max},
+    .lower = {ctx.psi_min, ctx.Z_min},
+    .cells = { cells_x[0], cells_x[1] },
+    .poly_order = ctx.poly_order,
+    .basis_type = app_args.basis_type,
+
+    .geometry = {
+      .geometry_id = GKYL_GEOMETRY_MIRROR,
+      .world = {0.0},
+      .mirror_grid_info = grid_inp,
+    },
+
+    .num_periodic_dir = 0,
+    .periodic_dirs = {},
+
+    .num_species = 2,
+    .species = {ion, elc},
+
+    .field = field,
+
+    .parallelism = {
+      .use_gpu = app_args.use_gpu,
+      .cuts = { app_args.cuts[0], app_args.cuts[1] },
+      .comm = comm,
+    },
+  };
+  struct gkyl_gyrokinetic_run_inp run_inp = {
+    .app_inp = app_inp, 
+    .time_stepping = {
+      .t_end = 5e-9,
+      .num_frames = 50,
+      .write_phase_freq = 1,
+      .int_diag_calc_num = 50*100,
+      .dt_failure_tol = ctx.dt_failure_tol,
+      .num_failures_max = ctx.num_failures_max,
+      .is_restart = false,
+      .restart_frame = 0,
+      .num_steps = app_args.num_steps,
+    },
+    .print_verbosity = {
+      .disable_timings = true,
+      .estimate_completion_time = true,
+      .enabled = true,
+      .frequency = 1,
+    },
+  };
+
+  gkyl_gyrokinetic_run_simulation(&run_inp);
+
+  gkyl_gyrokinetic_comms_release(comm);
+  
+#ifdef GKYL_HAVE_MPI
+  if (app_args.use_mpi)
+    MPI_Finalize();
+#endif
+  return 0;
+}

From d9808a785e1af1dd553a4a29a62e0c6a7fe05ea8 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Tue, 12 May 2026 15:08:09 -0400
Subject: [PATCH 02/11] Add support for per-cross-species diagnostics in LBO
 collisions

---
 gyrokinetic/apps/gk_species_lbo.c        | 112 ++++++++++++++++++-----
 gyrokinetic/apps/gkyl_gyrokinetic.h      |   1 +
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h |   7 +-
 gyrokinetic/creg/rt_gk_mirror_failure.c  |   1 +
 4 files changed, 93 insertions(+), 28 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_lbo.c b/gyrokinetic/apps/gk_species_lbo.c
index a2b1e038a3..b9d85d5e0c 100644
--- a/gyrokinetic/apps/gk_species_lbo.c
+++ b/gyrokinetic/apps/gk_species_lbo.c
@@ -107,43 +107,40 @@ gklbo_cross_moms_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gk
   // Empty method.
 }
 
+// Compute unscaled cross primitive moments (u_par_sr, vtsq_sr) for cross-collision species coll_idx.
+// Result is written into lbo->nu_boundary_corrections, which is used as scratch.
+static void
+gklbo_calc_cross_prim_moms(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_lbo_collisions *lbo, int coll_idx)
+{
+  lbo->alpha_E_func(app, gks, lbo, coll_idx);
+  for (int d=0; d<3; d++)
+    gkyl_dg_mul_op(app->basis, d, lbo->nu_moms, d, lbo->moms.marr, 0, lbo->cross_nu[coll_idx]);
+  for (int d=0; d<2; d++)
+    gkyl_dg_mul_op(app->basis, d, lbo->nu_boundary_corrections, d, lbo->boundary_corrections, 0, lbo->cross_nu[coll_idx]);
+  gkyl_prim_lbo_cross_calc_advance(lbo->cross_calc, &app->local, lbo->alpha_E,
+    gks->info.mass, lbo->nu_moms, lbo->prim_moms,
+    lbo->other_m[coll_idx], lbo->collide_with[coll_idx]->lbo.moms.marr, lbo->other_prim_moms[coll_idx],
+    lbo->nu_boundary_corrections, lbo->cross_nu[coll_idx], lbo->nu_boundary_corrections);
+}
+
 static void
 gklbo_cross_moms_enabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_lbo_collisions *lbo)
 {
-  // Compute primitive moments for cross-species collisions.
   struct timespec wst = gkyl_wall_clock();
-  
-  for (int i=0; i<lbo->num_cross_collisions; ++i) {
 
-    // Compute the cross-species collision frequency.
+  for (int i=0; i<lbo->num_cross_collisions; ++i) {
     lbo->cross_nu_func(app, gks, lbo, i);
+    gklbo_calc_cross_prim_moms(app, gks, lbo, i);
 
-    // Compute alpha_E.
-    lbo->alpha_E_func(app, gks, lbo, i);
-
-    // Multiply moments and boundary corrections by cross nu.
-    for (int d=0; d<3; d++)
-      gkyl_dg_mul_op(app->basis, d, lbo->nu_moms, d, lbo->moms.marr, 0, lbo->cross_nu[i]);
-    for (int d=0; d<2; d++)
-      gkyl_dg_mul_op(app->basis, d, lbo->nu_boundary_corrections, d, lbo->boundary_corrections, 0, lbo->cross_nu[i]);
-
-    // Compute cross primitive moments.
-    // Recycle the boundary_corrections array because we don't need those anymore.
+    // Scale upar_{sr} and vtSq_{sr} by nu_{sr} and accumulate.
     struct gkyl_array *cross_prim_moms = lbo->nu_boundary_corrections;
-    gkyl_prim_lbo_cross_calc_advance(lbo->cross_calc, &app->local, lbo->alpha_E, 
-      gks->info.mass, lbo->nu_moms, lbo->prim_moms,
-      lbo->other_m[i], lbo->collide_with[i]->lbo.moms.marr, lbo->other_prim_moms[i],
-      lbo->nu_boundary_corrections, lbo->cross_nu[i], cross_prim_moms);
-
-    // Scale upar_{sr} and vtSq_{sr} by nu_{sr}.
     for (int d=0; d<2; d++)
       gkyl_dg_mul_op(app->basis, d, cross_prim_moms, d, cross_prim_moms, 0, lbo->cross_nu[i]);
-
     gkyl_array_accumulate(lbo->nu_prim_moms, 1.0, cross_prim_moms);
-
   }
-  app->stat.species_coll_mom_tm += gkyl_time_diff_now_sec(wst);    
+  app->stat.species_coll_mom_tm += gkyl_time_diff_now_sec(wst);
 }
 
 static void
@@ -204,7 +201,55 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
   gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt, gks->lbo.nu_prim_moms_host, fileNm_nu_prim);
   app->stat.n_diag_io += 2;
 
-  gkyl_msgpack_data_release(mt); 
+  // Per-cross-species diagnostics: unscaled cross primitive moments and cross-nu.
+  if (gks->lbo.write_cross_diagnostics) {
+    // Self primitive moments (u_par_s, vtsq_s) before any cross-species contribution.
+    const char *fmt_prim = "%s-%s_lbo_self_prim_moms_%d.gkyl";
+    int sz_prim = gkyl_calc_strlen(fmt_prim, app->name, gks->info.name, frame);
+    char fileNm_prim[sz_prim+1];
+    snprintf(fileNm_prim, sizeof fileNm_prim, fmt_prim, app->name, gks->info.name, frame);
+    struct gkyl_array *prim_moms_io;
+    if (app->use_gpu) {
+      gkyl_array_copy(gks->lbo.nu_prim_moms_host, gks->lbo.prim_moms);
+      prim_moms_io = gks->lbo.nu_prim_moms_host;
+    } else {
+      prim_moms_io = gks->lbo.prim_moms;
+    }
+    gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt, prim_moms_io, fileNm_prim);
+    app->stat.n_diag_io += 1;
+    for (int i=0; i<gks->lbo.num_cross_collisions; ++i) {
+      const char *other_name = gks->lbo.collide_with[i]->info.name;
+
+      const char *fmt_cprim = "%s-%s_lbo_cross_%s_prim_moms_%d.gkyl";
+      int sz_cprim = gkyl_calc_strlen(fmt_cprim, app->name, gks->info.name, other_name, frame);
+      char fileNm_cprim[sz_cprim+1];
+      snprintf(fileNm_cprim, sizeof fileNm_cprim, fmt_cprim, app->name, gks->info.name, other_name, frame);
+
+      const char *fmt_cnu = "%s-%s_lbo_cross_%s_nu_%d.gkyl";
+      int sz_cnu = gkyl_calc_strlen(fmt_cnu, app->name, gks->info.name, other_name, frame);
+      char fileNm_cnu[sz_cnu+1];
+      snprintf(fileNm_cnu, sizeof fileNm_cnu, fmt_cnu, app->name, gks->info.name, other_name, frame);
+
+      gklbo_calc_cross_prim_moms(app, gks, &gks->lbo, i);
+      struct gkyl_array *cross_prim_moms_scratch = gks->lbo.nu_boundary_corrections;
+
+      struct gkyl_array *cross_prim_moms_io;
+      if (app->use_gpu) {
+        gkyl_array_copy(gks->lbo.nu_prim_moms_host, cross_prim_moms_scratch);
+        cross_prim_moms_io = gks->lbo.nu_prim_moms_host;
+        gkyl_array_copy(gks->lbo.cross_nu_host[i], gks->lbo.cross_nu[i]);
+      }
+      else {
+        cross_prim_moms_io = cross_prim_moms_scratch;
+      }
+
+      gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt, cross_prim_moms_io, fileNm_cprim);
+      gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt, gks->lbo.cross_nu_host[i], fileNm_cnu);
+      app->stat.n_diag_io += 2;
+    }
+  }
+
+  gkyl_msgpack_data_release(mt);
   app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wtm);
 }
 
@@ -213,6 +258,7 @@ gk_species_lbo_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks, st
 {
   lbo->collision_id = gks->info.collisions.collision_id;
   lbo->write_diagnostics = gks->info.collisions.write_diagnostics;
+  lbo->write_cross_diagnostics = gks->info.collisions.write_cross_diagnostics;
 
   // Empty methods.
   lbo->moms_func = gklbo_moms_disabled;
@@ -368,6 +414,17 @@ gk_species_lbo_cross_init(struct gkyl_gyrokinetic_app *app, struct gk_species *g
         lbo->other_prim_moms[i] = lbo->collide_with[i]->lbo.prim_moms;
       }
 
+      if (lbo->write_cross_diagnostics) {
+        for (int i=0; i<lbo->num_cross_collisions; ++i) {
+          if (app->use_gpu) {
+            lbo->cross_nu_host[i] = mkarr(false, lbo->cross_nu[i]->ncomp, lbo->cross_nu[i]->size);
+          }
+          else {
+            lbo->cross_nu_host[i] = lbo->cross_nu[i];
+          }
+        }
+      }
+
       double nu_frac = gks->info.collisions.nu_frac ? gks->info.collisions.nu_frac : 1.0;
 
       // Compute the time-independent part of alpha_E.
@@ -497,6 +554,11 @@ gk_species_lbo_release(const struct gkyl_gyrokinetic_app *app, const struct gk_l
         gkyl_array_release(lbo->cross_nu[i]);
 
       gkyl_array_release(lbo->alpha_E);
+
+      if (lbo->write_cross_diagnostics && app->use_gpu) {
+        for (int i=0; i<lbo->num_cross_collisions; ++i)
+          gkyl_array_release(lbo->cross_nu_host[i]);
+      }
     }
 
     gkyl_dg_updater_lbo_gyrokinetic_release(lbo->coll_slvr);
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 77d82317ca..ca6eefb085 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -106,6 +106,7 @@ struct gkyl_gyrokinetic_collisionless {
 struct gkyl_gyrokinetic_collisions {
   enum gkyl_collision_id collision_id; // type of collisions (see gkyl_eqn_type.h)
   bool write_diagnostics; // Whether to output diagnostics.
+  bool write_cross_diagnostics; // Whether to output per-cross-species diagnostics (prim_moms_sr, nu_sr).
   bool not_in_dfdt; // If true, the collision operator will not be added to df/dt.
     // Used to ignore the collisional updates of this species, while updating cross-species collisions.
 
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 0d26fe4b6b..54e16770fa 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -321,10 +321,11 @@ struct gk_collisionless {
     struct gk_collisionless *gkcls, double tm, int frame);
 };
 
-struct gk_lbo_collisions {  
+struct gk_lbo_collisions {
   enum gkyl_collision_id collision_id; // type of collisions
   bool write_diagnostics; // Whether to write diagnostics out.
-  bool not_in_dfdt; // Whether to not add collision contribution to df/dt 
+  bool write_cross_diagnostics; // Whether to write per-cross-species diagnostics (prim_moms_sr, nu_sr).
+  bool not_in_dfdt; // Whether to not add collision contribution to df/dt
 
   struct gkyl_array *self_nu; // Self-collision frequency.
   struct gkyl_array *boundary_corrections; // LBO boundary corrections.
@@ -350,8 +351,8 @@ struct gk_lbo_collisions {
   double delta_sr; // Free parameter in relationship between alpha_E and nu_sr.
   double other_m[GKYL_MAX_SPECIES]; // Masses of species colliding with.
   struct gkyl_array *other_prim_moms[GKYL_MAX_SPECIES]; // Self-primitive moments of species colliding with.
-  struct gkyl_array *cross_prim_moms[GKYL_MAX_SPECIES]; // Cross-primitive moments.
   struct gkyl_array *cross_nu[GKYL_MAX_SPECIES]; // Cross-species collision frequencies.
+  struct gkyl_array *cross_nu_host[GKYL_MAX_SPECIES]; // Host arrays for I/O.
   struct gkyl_array *cross_nu_prim_moms; // Weak multiplication of collision frequency and primitive moments.
   struct gkyl_array *alpha_E; // Morse's alpha_E factor.
   gkyl_prim_lbo_cross_calc *cross_calc; // LBO cross-primitive moment calculator
diff --git a/gyrokinetic/creg/rt_gk_mirror_failure.c b/gyrokinetic/creg/rt_gk_mirror_failure.c
index 1de04ededc..586966e308 100644
--- a/gyrokinetic/creg/rt_gk_mirror_failure.c
+++ b/gyrokinetic/creg/rt_gk_mirror_failure.c
@@ -345,6 +345,7 @@ int main(int argc, char **argv)
       .num_cross_collisions = 1,
       .collide_with = { "elc" },
       .write_diagnostics = true,
+      .write_cross_diagnostics = true,
     },
 
     .write_omega_cfl = true,

From 68befa637543031144f0378d84f7c49a689cfbee Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Wed, 13 May 2026 10:15:50 -0400
Subject: [PATCH 03/11] Enhance cross primitive moments calculation in LBO
 collisions and remove unused mirror failure file

---
 gyrokinetic/apps/gk_species_lbo.c       |   5 +
 gyrokinetic/creg/rt_gk_mirror_failure.c | 482 ------------------------
 2 files changed, 5 insertions(+), 482 deletions(-)
 delete mode 100644 gyrokinetic/creg/rt_gk_mirror_failure.c

diff --git a/gyrokinetic/apps/gk_species_lbo.c b/gyrokinetic/apps/gk_species_lbo.c
index b9d85d5e0c..b47f84031e 100644
--- a/gyrokinetic/apps/gk_species_lbo.c
+++ b/gyrokinetic/apps/gk_species_lbo.c
@@ -114,10 +114,15 @@ gklbo_calc_cross_prim_moms(gkyl_gyrokinetic_app *app, const struct gk_species *g
   struct gk_lbo_collisions *lbo, int coll_idx)
 {
   lbo->alpha_E_func(app, gks, lbo, coll_idx);
+
+  // Multiply moments and boundary corrections by cross nu.
   for (int d=0; d<3; d++)
     gkyl_dg_mul_op(app->basis, d, lbo->nu_moms, d, lbo->moms.marr, 0, lbo->cross_nu[coll_idx]);
   for (int d=0; d<2; d++)
     gkyl_dg_mul_op(app->basis, d, lbo->nu_boundary_corrections, d, lbo->boundary_corrections, 0, lbo->cross_nu[coll_idx]);
+    
+  // Compute cross primitive moments.
+  // Recycle the boundary_corrections array because we don't need those anymore.
   gkyl_prim_lbo_cross_calc_advance(lbo->cross_calc, &app->local, lbo->alpha_E,
     gks->info.mass, lbo->nu_moms, lbo->prim_moms,
     lbo->other_m[coll_idx], lbo->collide_with[coll_idx]->lbo.moms.marr, lbo->other_prim_moms[coll_idx],
diff --git a/gyrokinetic/creg/rt_gk_mirror_failure.c b/gyrokinetic/creg/rt_gk_mirror_failure.c
deleted file mode 100644
index 586966e308..0000000000
--- a/gyrokinetic/creg/rt_gk_mirror_failure.c
+++ /dev/null
@@ -1,482 +0,0 @@
-#include <math.h>
-#include <stdio.h>
-#include <time.h>
-
-#include <gkyl_alloc.h>
-#include <gkyl_const.h>
-#include <gkyl_eqn_type.h>
-#include <gkyl_fem_poisson_bctype.h>
-#include <gkyl_gyrokinetic.h>
-#include <gkyl_gyrokinetic_run.h>
-#include <gkyl_math.h>
-
-#include <rt_arg_parse.h>
-
-// Define the context of the simulation. This is basically all the globals
-struct gk_mirror_ctx
-{
-  int cdim, vdim; // Dimensionality.
-  // Plasma parameters
-  double mi;
-  double qi;
-  double me;
-  double qe;
-  double Te0;
-  double n0;
-  double B_p; // Hardcoded magnetic field at midplane 0.53 for double lorentzian
-  double Bmag_midp; // Magnetic field magnitude at midplane (Z=0)
-  double beta;
-  double tau;
-  double Ti0;
-  double nuFrac;
-  // Ion-ion collision freq.
-  double logLambdaIon;
-  double nuIon;
-  double vti, vte;
-  double RatZeq0; // Radius of the field line at Z=0.
-  double kperp; // Perpendicular wavenumber for ES gyrokinetics.
-  // Axial coordinate Z extents. Endure that Z=0 is not on
-  double z_min;
-  double z_max;
-  double psi_eval;
-  double psi_max;
-  double psi_min;
-  double theta_eval;
-  double theta_min;
-  double theta_max;
-  // Physics parameters at mirror throat
-  double vpar_max_ion;
-  double mu_max_ion;
-  double vpar_max_elc;
-  double mu_max_elc;
-
-  int Npsi;
-  int Ntheta;
-  int Nz;
-  int Nvpar;
-  int Nmu;
-  int Nvpar_elc;
-  int Nmu_elc;
-  int cells[GKYL_MAX_DIM]; // Number of cells in all directions.
-  int poly_order;
-
-  // Source parameters
-  double ion_source_amplitude;
-  double ion_source_sigma;
-  double ion_source_temp;
-
-  double t_end; // End time.
-  int num_frames; // Number of output frames.
-  int num_phases; // Number of phases.
-  double write_phase_freq; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
-  double int_diag_calc_freq; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
-  double dt_failure_tol; // Minimum allowable fraction of initial time-step.
-  int num_failures_max; // Maximum allowable number of consecutive small time-steps.
-
-  // Geometry parameters for Lorentzian mirror
-  double mcB;     // Magnetic field parameter
-  double mcB_inner; // Magnetic field parameter for inner mirror
-  double mcB_outer; // Magnetic field parameter for outer mirror
-  double gamma;   // Width parameter for Lorentzian profile
-  double gamma_inner; // Width parameter for Lorentzian profile for inner mirror
-  double gamma_outer; // Width parameter for Lorentzian profile for outer mirror
-  double Z_m;     // Mirror throat location
-  double L_center; // Length of central region
-  double L_plugs; // Length of plugs
-  double Z_min;   // Minimum Z coordinate
-  double Z_max;   // Maximum Z coordinate
-  double psi_in;  // Working variable for psi integration
-  double z_in;    // Working variable for z integration
-};
-
-
-// Evaluate initial conditions
-// I think ICs should be constructed to match the expander rather than the center
-// as the center is quickly filled in the OAP
-void
-initial_density(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  fout[0] = 1e17;
-}
-
-void
-initial_upar(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  double z = xn[1];
-  double c_s = 7 * sqrt(app->Te0/app->mi);
-  if (fabs(z) <= app->Z_m)
-  {
-    fout[0] = 0.0;
-  }
-  else
-  {
-    fout[0] = fabs(z) / z * c_s * tanh(4 * (app->Z_max - app->Z_m) * fabs(fabs(z) - app->Z_m));
-  }
-}
-
-void
-eval_zero(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  fout[0] = 0.0;
-}
-
-void
-initial_temp_ion(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  fout[0] = app->Ti0/10.0;
-}
-
-void
-initial_temp_elc(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  fout[0] = app->Te0;
-}
-
-void mapc2p_vel_ion(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  double vpar_max_ion = app->vpar_max_ion;
-  double mu_max_ion = app->mu_max_ion;
-
-  double cvpar = vc[0], cmu = vc[1];
-  double b = 1.4;
-  vp[0] = vpar_max_ion*tan(cvpar*b)/tan(b);
-  vp[1] = mu_max_ion*pow(cmu,2);  // Cubic map in mu.
-}
-
-void mapc2p_vel_elc(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  double mu_max_elc = app->mu_max_elc;
-  double vpar_max_elc = app->vpar_max_elc;
-  double cvpar = vc[0], cmu = vc[1];
-  vp[0] = vpar_max_elc*cvpar;
-  vp[1] = mu_max_elc*pow(cmu,4);  // Cubic map in mu.
-}
-
-struct gk_mirror_ctx
-create_ctx(void)
-{
-  int cdim = 2, vdim = 2; // Dimensionality.
-  int poly_order = 1;
-
-  // Universal constant parameters.
-  double eps0 = GKYL_EPSILON0;
-  double mu0 = GKYL_MU0; // Not sure if this is right
-  double eV = GKYL_ELEMENTARY_CHARGE;
-  double mp = GKYL_PROTON_MASS; // ion mass
-  double me = GKYL_ELECTRON_MASS;
-  double qi = eV;  // ion charge
-  double qe = -eV; // electron charge
-
-  // Plasma parameters.
-  double mi = 2.014 * mp;
-  double Te0 = 940 * eV;
-  double n0 = 3e19;
-  double B_p = 0.53; // Bmag at z=0
-  double beta = 0.4;
-  double tau = pow(B_p, 2.) * beta / (2.0 * mu0 * n0 * Te0) - 1.;
-  double Ti0 = tau * Te0;
-
-  // Ion-ion collision freq.
-  double nuFrac = 1.0;
-  double logLambdaIon = 6.6 - 0.5 * log(n0 / 1e20) + 1.5 * log(Ti0 / eV);
-  double nuIon = nuFrac * logLambdaIon * pow(eV, 4.) * n0 /
-                 (12 * pow(M_PI, 3. / 2.) * pow(eps0, 2.) * sqrt(mi) * pow(Ti0, 3. / 2.));
-
-  // Thermal speeds.
-  double vti = sqrt(Ti0 / mi);
-  double vte = sqrt(Te0 / me);
-
-  // Grid parameters
-  double vpar_max_ion = 16 * vti;
-  double mu_max_ion = mi * pow(3. * vti, 2.) / (2. * B_p);
-  double vpar_max_elc = 4 * vte;
-  double mu_max_elc = me * pow(4. * vte, 2.) / (2. * B_p);
-
-  int Nz = 32;
-  int Npsi = 2;
-  int Nvpar = 16;
-  int Nmu = 8;
-  int Nvpar_elc = 8;
-  int Nmu_elc = 8;
-
-  // Geometry parameters.
-  double RatZeq0 = 0.10; // Radius of the field line at Z=0.
-  double Z_min = -2.5;
-  double Z_max =  2.5;
-  double mcB = 3.691260;
-  double gamma = 0.226381;
-  double Z_m = 0.98;
-
-  // Calculate phase structure
-  double write_phase_freq = 1; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
-  double int_diag_calc_freq = 100; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
-  double dt_failure_tol = 1.0e-4; // Minimum allowable fraction of initial time-step.
-  int num_failures_max = 20; // Maximum allowable number of consecutive small time-steps.
-
-  struct gk_mirror_ctx ctx = {
-    .cdim = cdim,
-    .vdim = vdim,
-    .poly_order = poly_order,
-
-    .mi = mi,
-    .qi = qi,
-    .me = me,
-    .qe = qe,
-    .Te0 = Te0,
-    .n0 = n0,
-    .B_p = B_p,
-    .beta = beta,
-    .tau = tau,
-    .Ti0 = Ti0,
-
-    .nuFrac = nuFrac,
-    .logLambdaIon = logLambdaIon,
-    .nuIon = nuIon,
-    .vti = vti,
-    .vte = vte,
-    .RatZeq0 = RatZeq0,
-    .vpar_max_ion = vpar_max_ion,
-    .mu_max_ion = mu_max_ion,
-    .vpar_max_elc = vpar_max_elc,
-    .mu_max_elc = mu_max_elc,
-
-    .Npsi = Npsi,
-    .Nz = Nz,
-    .Nvpar = Nvpar,
-    .Nmu = Nmu,
-    .Nvpar_elc = Nvpar_elc,
-    .Nmu_elc = Nmu_elc,
-    .cells = {Npsi, Nz, Nvpar, Nmu},
-    .Z_min = Z_min,
-    .Z_max = Z_max,
-
-    .write_phase_freq = write_phase_freq,
-    .int_diag_calc_freq = int_diag_calc_freq,
-    .dt_failure_tol = dt_failure_tol,
-    .num_failures_max = num_failures_max,
-
-    .mcB = mcB,
-    .gamma = gamma,
-    .Z_m = Z_m,
-  };
-  
-  // Populate a couple more values in the context.
-  ctx.psi_max = 3e-3;
-  ctx.psi_min = 1e-5;
-  ctx.psi_eval = (ctx.psi_max + ctx.psi_min) / 2.0;
-
-  ctx.z_min    = -2.5;
-  ctx.z_max    = 2.5;
-
-  return ctx;
-}
-
-int main(int argc, char **argv)
-{
-  struct gkyl_app_args app_args = parse_app_args(argc, argv);
-
-#ifdef GKYL_HAVE_MPI
-  if (app_args.use_mpi) MPI_Init(&argc, &argv);
-#endif
-
-  if (app_args.trace_mem) {
-    gkyl_cu_dev_mem_debug_set(true);
-    gkyl_mem_debug_set(true);
-  }
-
-  struct gk_mirror_ctx ctx = create_ctx(); // Context for init functions.
-
-  int rank = 0;
-#ifdef GKYL_HAVE_MPI
-  if (app_args.use_mpi) MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-#endif
-  // if (rank == 0)
-  //   print_ctx(&ctx);
-
-  int cells_x[ctx.cdim], cells_v[ctx.vdim];
-  for (int d=0; d<ctx.cdim; d++)
-    cells_x[d] = APP_ARGS_CHOOSE(app_args.xcells[d], ctx.cells[d]);
-  for (int d=0; d<ctx.vdim; d++)
-    cells_v[d] = APP_ARGS_CHOOSE(app_args.vcells[d], ctx.cells[ctx.cdim+d]);
-
-  // Construct communicator for use in app.
-  struct gkyl_comm *comm = gkyl_gyrokinetic_comms_new(app_args.use_mpi, app_args.use_gpu, stderr);
-
-  struct gkyl_gyrokinetic_species ion = {
-    .name = "ion",
-    .charge = ctx.qi,
-    .mass = ctx.mi,
-    .vdim = ctx.vdim,
-    .lower = {-1.0, 0.0},
-    .upper = { 1.0, 1.0},
-    .cells = { cells_v[0], cells_v[1]},
-    .polarization_density = ctx.n0,
-
-    .projection = {
-      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
-      .density = initial_density,
-      .ctx_density = &ctx,
-      .upar = initial_upar,
-      .ctx_upar = &ctx,
-      .temp = initial_temp_ion,
-      .ctx_temp = &ctx,
-    },
-
-    .mapc2p = {
-      .mapping = mapc2p_vel_ion,
-      .ctx = &ctx,
-    },
-
-    .collisionless = {
-      .type = GKYL_GK_COLLISIONLESS_ES,
-      .write_diagnostics = true,
-    },
-
-    .collisions = {
-      .collision_id = GKYL_LBO_COLLISIONS,
-      .den_ref = ctx.n0,
-      .temp_ref = ctx.Ti0,
-      .num_cross_collisions = 1,
-      .collide_with = { "elc" },
-      .write_diagnostics = true,
-      .write_cross_diagnostics = true,
-    },
-
-    .write_omega_cfl = true,
-    .num_diag_moments = 8,
-    .diag_moments = {GKYL_F_MOMENT_BIMAXWELLIAN, GKYL_F_MOMENT_M0, GKYL_F_MOMENT_M1, GKYL_F_MOMENT_M2, GKYL_F_MOMENT_M2PAR, GKYL_F_MOMENT_M2PERP, GKYL_F_MOMENT_M3PAR, GKYL_F_MOMENT_M3PERP },
-    .num_integrated_diag_moments = 1,
-    .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP },
-    .time_rate_diagnostics = true,
-
-    .boundary_flux_diagnostics = {
-      .num_integrated_diag_moments = 1,
-      .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP},
-    },
-  };
-
-
-  struct gkyl_gyrokinetic_species elc = {
-    .name = "elc",
-    .charge = ctx.qe,
-    .mass = ctx.me,
-    .vdim = ctx.vdim,
-    .lower = {-1.0, 0.0},
-    .upper = { 1.0, 1.0},
-    .cells = { ctx.Nvpar_elc, ctx.Nmu_elc },
-
-    .polarization_density = ctx.n0,
-
-    .mapc2p = {
-      .mapping = mapc2p_vel_elc,
-      .ctx = &ctx,
-    },
-
-    .projection = {
-      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
-      .density = initial_density,
-      .ctx_density = &ctx,
-      .upar = eval_zero,
-      .ctx_upar = &ctx,
-      .temp = initial_temp_elc,
-      .ctx_temp = &ctx,
-      .correct_all_moms = true,
-    },
-
-    .collisions =  {
-      .collision_id = GKYL_LBO_COLLISIONS,
-      .den_ref = ctx.n0,
-      .temp_ref = ctx.Te0,
-      .num_cross_collisions = 1,
-      .collide_with = { "ion" },
-      .write_diagnostics = true,
-      .not_in_dfdt = true,
-    },
-
-    .num_diag_moments = 1,
-    .diag_moments = {GKYL_F_MOMENT_MAXWELLIAN},
-  };
-
-  struct gkyl_gyrokinetic_field field = {
-    .gkfield_id = GKYL_GK_FIELD_BOLTZMANN,
-    .electron_mass = ctx.me,
-    .electron_charge = ctx.qe,
-    .electron_temp = ctx.Te0,
-    .is_static = false,
-  };
-
-  struct gkyl_mirror_geo_grid_inp grid_inp = {
-    .filename_psi = "gyrokinetic/data/eqdsk/lorentzian_R10.geqdsk_psi.gkyl", // psi file to use
-    .rclose = 0.2, // closest R to region of interest
-    .zmin = -2.5,  // Z of lower boundary
-    .zmax =  2.5,  // Z of upper boundary
-    .include_axis = false, // Include R=0 axis in grid
-    .fl_coord = GKYL_GEOMETRY_MIRROR_GRID_GEN_SQRT_PSI_CART_Z, // coordinate system for psi grid
-  };
-
-  struct gkyl_gk app_inp = {  // GK app
-    .name = "gk_lorentzian_mirror",
-    .cdim = ctx.cdim,
-    .upper = {ctx.psi_max, ctx.Z_max},
-    .lower = {ctx.psi_min, ctx.Z_min},
-    .cells = { cells_x[0], cells_x[1] },
-    .poly_order = ctx.poly_order,
-    .basis_type = app_args.basis_type,
-
-    .geometry = {
-      .geometry_id = GKYL_GEOMETRY_MIRROR,
-      .world = {0.0},
-      .mirror_grid_info = grid_inp,
-    },
-
-    .num_periodic_dir = 0,
-    .periodic_dirs = {},
-
-    .num_species = 2,
-    .species = {ion, elc},
-
-    .field = field,
-
-    .parallelism = {
-      .use_gpu = app_args.use_gpu,
-      .cuts = { app_args.cuts[0], app_args.cuts[1] },
-      .comm = comm,
-    },
-  };
-  struct gkyl_gyrokinetic_run_inp run_inp = {
-    .app_inp = app_inp, 
-    .time_stepping = {
-      .t_end = 5e-9,
-      .num_frames = 50,
-      .write_phase_freq = 1,
-      .int_diag_calc_num = 50*100,
-      .dt_failure_tol = ctx.dt_failure_tol,
-      .num_failures_max = ctx.num_failures_max,
-      .is_restart = false,
-      .restart_frame = 0,
-      .num_steps = app_args.num_steps,
-    },
-    .print_verbosity = {
-      .disable_timings = true,
-      .estimate_completion_time = true,
-      .enabled = true,
-      .frequency = 1,
-    },
-  };
-
-  gkyl_gyrokinetic_run_simulation(&run_inp);
-
-  gkyl_gyrokinetic_comms_release(comm);
-  
-#ifdef GKYL_HAVE_MPI
-  if (app_args.use_mpi)
-    MPI_Finalize();
-#endif
-  return 0;
-}

From f7a48880fb5575139ad61045b6cafa92c3ee2309 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Wed, 13 May 2026 10:20:31 -0400
Subject: [PATCH 04/11] Add cross-primitive moments array to gk_lbo_collisions
 struct

---
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 54e16770fa..dc9562223b 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -351,6 +351,7 @@ struct gk_lbo_collisions {
   double delta_sr; // Free parameter in relationship between alpha_E and nu_sr.
   double other_m[GKYL_MAX_SPECIES]; // Masses of species colliding with.
   struct gkyl_array *other_prim_moms[GKYL_MAX_SPECIES]; // Self-primitive moments of species colliding with.
+  struct gkyl_array *cross_prim_moms[GKYL_MAX_SPECIES]; // Cross-primitive moments.
   struct gkyl_array *cross_nu[GKYL_MAX_SPECIES]; // Cross-species collision frequencies.
   struct gkyl_array *cross_nu_host[GKYL_MAX_SPECIES]; // Host arrays for I/O.
   struct gkyl_array *cross_nu_prim_moms; // Weak multiplication of collision frequency and primitive moments.

From 7d539c3610436f600de21f1f561568889bea834f Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Wed, 13 May 2026 10:23:29 -0400
Subject: [PATCH 05/11] Refactor cross-nu handling in gklbo_write_mom_enabled
 for improved I/O efficiency

---
 gyrokinetic/apps/gk_species_lbo.c        | 23 ++++-------------------
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h |  1 -
 2 files changed, 4 insertions(+), 20 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_lbo.c b/gyrokinetic/apps/gk_species_lbo.c
index b47f84031e..31b8d38c2a 100644
--- a/gyrokinetic/apps/gk_species_lbo.c
+++ b/gyrokinetic/apps/gk_species_lbo.c
@@ -238,18 +238,19 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
       gklbo_calc_cross_prim_moms(app, gks, &gks->lbo, i);
       struct gkyl_array *cross_prim_moms_scratch = gks->lbo.nu_boundary_corrections;
 
-      struct gkyl_array *cross_prim_moms_io;
+      struct gkyl_array *cross_prim_moms_io, *cross_nu_io;
       if (app->use_gpu) {
         gkyl_array_copy(gks->lbo.nu_prim_moms_host, cross_prim_moms_scratch);
         cross_prim_moms_io = gks->lbo.nu_prim_moms_host;
-        gkyl_array_copy(gks->lbo.cross_nu_host[i], gks->lbo.cross_nu[i]);
+        gkyl_array_copy(gks->lbo.nu_sum_host, gks->lbo.cross_nu[i]);
       }
       else {
         cross_prim_moms_io = cross_prim_moms_scratch;
+        cross_nu_io = gks->lbo.cross_nu[i];
       }
 
       gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt, cross_prim_moms_io, fileNm_cprim);
-      gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt, gks->lbo.cross_nu_host[i], fileNm_cnu);
+      gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt, cross_nu_io, fileNm_cnu);
       app->stat.n_diag_io += 2;
     }
   }
@@ -419,17 +420,6 @@ gk_species_lbo_cross_init(struct gkyl_gyrokinetic_app *app, struct gk_species *g
         lbo->other_prim_moms[i] = lbo->collide_with[i]->lbo.prim_moms;
       }
 
-      if (lbo->write_cross_diagnostics) {
-        for (int i=0; i<lbo->num_cross_collisions; ++i) {
-          if (app->use_gpu) {
-            lbo->cross_nu_host[i] = mkarr(false, lbo->cross_nu[i]->ncomp, lbo->cross_nu[i]->size);
-          }
-          else {
-            lbo->cross_nu_host[i] = lbo->cross_nu[i];
-          }
-        }
-      }
-
       double nu_frac = gks->info.collisions.nu_frac ? gks->info.collisions.nu_frac : 1.0;
 
       // Compute the time-independent part of alpha_E.
@@ -559,11 +549,6 @@ gk_species_lbo_release(const struct gkyl_gyrokinetic_app *app, const struct gk_l
         gkyl_array_release(lbo->cross_nu[i]);
 
       gkyl_array_release(lbo->alpha_E);
-
-      if (lbo->write_cross_diagnostics && app->use_gpu) {
-        for (int i=0; i<lbo->num_cross_collisions; ++i)
-          gkyl_array_release(lbo->cross_nu_host[i]);
-      }
     }
 
     gkyl_dg_updater_lbo_gyrokinetic_release(lbo->coll_slvr);
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index dc9562223b..d42bae4a13 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -353,7 +353,6 @@ struct gk_lbo_collisions {
   struct gkyl_array *other_prim_moms[GKYL_MAX_SPECIES]; // Self-primitive moments of species colliding with.
   struct gkyl_array *cross_prim_moms[GKYL_MAX_SPECIES]; // Cross-primitive moments.
   struct gkyl_array *cross_nu[GKYL_MAX_SPECIES]; // Cross-species collision frequencies.
-  struct gkyl_array *cross_nu_host[GKYL_MAX_SPECIES]; // Host arrays for I/O.
   struct gkyl_array *cross_nu_prim_moms; // Weak multiplication of collision frequency and primitive moments.
   struct gkyl_array *alpha_E; // Morse's alpha_E factor.
   gkyl_prim_lbo_cross_calc *cross_calc; // LBO cross-primitive moment calculator

From f6ad28323f4cf81a8505457b406693211dda14d9 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Wed, 13 May 2026 10:36:37 -0400
Subject: [PATCH 06/11] Add test again

---
 gyrokinetic/creg/rt_gk_mirror_failure.c | 482 ++++++++++++++++++++++++
 1 file changed, 482 insertions(+)
 create mode 100644 gyrokinetic/creg/rt_gk_mirror_failure.c

diff --git a/gyrokinetic/creg/rt_gk_mirror_failure.c b/gyrokinetic/creg/rt_gk_mirror_failure.c
new file mode 100644
index 0000000000..586966e308
--- /dev/null
+++ b/gyrokinetic/creg/rt_gk_mirror_failure.c
@@ -0,0 +1,482 @@
+#include <math.h>
+#include <stdio.h>
+#include <time.h>
+
+#include <gkyl_alloc.h>
+#include <gkyl_const.h>
+#include <gkyl_eqn_type.h>
+#include <gkyl_fem_poisson_bctype.h>
+#include <gkyl_gyrokinetic.h>
+#include <gkyl_gyrokinetic_run.h>
+#include <gkyl_math.h>
+
+#include <rt_arg_parse.h>
+
+// Define the context of the simulation. This is basically all the globals
+struct gk_mirror_ctx
+{
+  int cdim, vdim; // Dimensionality.
+  // Plasma parameters
+  double mi;
+  double qi;
+  double me;
+  double qe;
+  double Te0;
+  double n0;
+  double B_p; // Hardcoded magnetic field at midplane 0.53 for double lorentzian
+  double Bmag_midp; // Magnetic field magnitude at midplane (Z=0)
+  double beta;
+  double tau;
+  double Ti0;
+  double nuFrac;
+  // Ion-ion collision freq.
+  double logLambdaIon;
+  double nuIon;
+  double vti, vte;
+  double RatZeq0; // Radius of the field line at Z=0.
+  double kperp; // Perpendicular wavenumber for ES gyrokinetics.
+  // Axial coordinate Z extents. Endure that Z=0 is not on
+  double z_min;
+  double z_max;
+  double psi_eval;
+  double psi_max;
+  double psi_min;
+  double theta_eval;
+  double theta_min;
+  double theta_max;
+  // Physics parameters at mirror throat
+  double vpar_max_ion;
+  double mu_max_ion;
+  double vpar_max_elc;
+  double mu_max_elc;
+
+  int Npsi;
+  int Ntheta;
+  int Nz;
+  int Nvpar;
+  int Nmu;
+  int Nvpar_elc;
+  int Nmu_elc;
+  int cells[GKYL_MAX_DIM]; // Number of cells in all directions.
+  int poly_order;
+
+  // Source parameters
+  double ion_source_amplitude;
+  double ion_source_sigma;
+  double ion_source_temp;
+
+  double t_end; // End time.
+  int num_frames; // Number of output frames.
+  int num_phases; // Number of phases.
+  double write_phase_freq; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
+  double int_diag_calc_freq; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
+  double dt_failure_tol; // Minimum allowable fraction of initial time-step.
+  int num_failures_max; // Maximum allowable number of consecutive small time-steps.
+
+  // Geometry parameters for Lorentzian mirror
+  double mcB;     // Magnetic field parameter
+  double mcB_inner; // Magnetic field parameter for inner mirror
+  double mcB_outer; // Magnetic field parameter for outer mirror
+  double gamma;   // Width parameter for Lorentzian profile
+  double gamma_inner; // Width parameter for Lorentzian profile for inner mirror
+  double gamma_outer; // Width parameter for Lorentzian profile for outer mirror
+  double Z_m;     // Mirror throat location
+  double L_center; // Length of central region
+  double L_plugs; // Length of plugs
+  double Z_min;   // Minimum Z coordinate
+  double Z_max;   // Maximum Z coordinate
+  double psi_in;  // Working variable for psi integration
+  double z_in;    // Working variable for z integration
+};
+
+
+// Evaluate initial conditions
+// I think ICs should be constructed to match the expander rather than the center
+// as the center is quickly filled in the OAP
+void
+initial_density(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  fout[0] = 1e17;
+}
+
+void
+initial_upar(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  double z = xn[1];
+  double c_s = 7 * sqrt(app->Te0/app->mi);
+  if (fabs(z) <= app->Z_m)
+  {
+    fout[0] = 0.0;
+  }
+  else
+  {
+    fout[0] = fabs(z) / z * c_s * tanh(4 * (app->Z_max - app->Z_m) * fabs(fabs(z) - app->Z_m));
+  }
+}
+
+void
+eval_zero(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  fout[0] = 0.0;
+}
+
+void
+initial_temp_ion(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  fout[0] = app->Ti0/10.0;
+}
+
+void
+initial_temp_elc(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  fout[0] = app->Te0;
+}
+
+void mapc2p_vel_ion(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  double vpar_max_ion = app->vpar_max_ion;
+  double mu_max_ion = app->mu_max_ion;
+
+  double cvpar = vc[0], cmu = vc[1];
+  double b = 1.4;
+  vp[0] = vpar_max_ion*tan(cvpar*b)/tan(b);
+  vp[1] = mu_max_ion*pow(cmu,2);  // Cubic map in mu.
+}
+
+void mapc2p_vel_elc(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
+{
+  struct gk_mirror_ctx *app = ctx;
+  double mu_max_elc = app->mu_max_elc;
+  double vpar_max_elc = app->vpar_max_elc;
+  double cvpar = vc[0], cmu = vc[1];
+  vp[0] = vpar_max_elc*cvpar;
+  vp[1] = mu_max_elc*pow(cmu,4);  // Cubic map in mu.
+}
+
+struct gk_mirror_ctx
+create_ctx(void)
+{
+  int cdim = 2, vdim = 2; // Dimensionality.
+  int poly_order = 1;
+
+  // Universal constant parameters.
+  double eps0 = GKYL_EPSILON0;
+  double mu0 = GKYL_MU0; // Not sure if this is right
+  double eV = GKYL_ELEMENTARY_CHARGE;
+  double mp = GKYL_PROTON_MASS; // ion mass
+  double me = GKYL_ELECTRON_MASS;
+  double qi = eV;  // ion charge
+  double qe = -eV; // electron charge
+
+  // Plasma parameters.
+  double mi = 2.014 * mp;
+  double Te0 = 940 * eV;
+  double n0 = 3e19;
+  double B_p = 0.53; // Bmag at z=0
+  double beta = 0.4;
+  double tau = pow(B_p, 2.) * beta / (2.0 * mu0 * n0 * Te0) - 1.;
+  double Ti0 = tau * Te0;
+
+  // Ion-ion collision freq.
+  double nuFrac = 1.0;
+  double logLambdaIon = 6.6 - 0.5 * log(n0 / 1e20) + 1.5 * log(Ti0 / eV);
+  double nuIon = nuFrac * logLambdaIon * pow(eV, 4.) * n0 /
+                 (12 * pow(M_PI, 3. / 2.) * pow(eps0, 2.) * sqrt(mi) * pow(Ti0, 3. / 2.));
+
+  // Thermal speeds.
+  double vti = sqrt(Ti0 / mi);
+  double vte = sqrt(Te0 / me);
+
+  // Grid parameters
+  double vpar_max_ion = 16 * vti;
+  double mu_max_ion = mi * pow(3. * vti, 2.) / (2. * B_p);
+  double vpar_max_elc = 4 * vte;
+  double mu_max_elc = me * pow(4. * vte, 2.) / (2. * B_p);
+
+  int Nz = 32;
+  int Npsi = 2;
+  int Nvpar = 16;
+  int Nmu = 8;
+  int Nvpar_elc = 8;
+  int Nmu_elc = 8;
+
+  // Geometry parameters.
+  double RatZeq0 = 0.10; // Radius of the field line at Z=0.
+  double Z_min = -2.5;
+  double Z_max =  2.5;
+  double mcB = 3.691260;
+  double gamma = 0.226381;
+  double Z_m = 0.98;
+
+  // Calculate phase structure
+  double write_phase_freq = 1; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
+  double int_diag_calc_freq = 100; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
+  double dt_failure_tol = 1.0e-4; // Minimum allowable fraction of initial time-step.
+  int num_failures_max = 20; // Maximum allowable number of consecutive small time-steps.
+
+  struct gk_mirror_ctx ctx = {
+    .cdim = cdim,
+    .vdim = vdim,
+    .poly_order = poly_order,
+
+    .mi = mi,
+    .qi = qi,
+    .me = me,
+    .qe = qe,
+    .Te0 = Te0,
+    .n0 = n0,
+    .B_p = B_p,
+    .beta = beta,
+    .tau = tau,
+    .Ti0 = Ti0,
+
+    .nuFrac = nuFrac,
+    .logLambdaIon = logLambdaIon,
+    .nuIon = nuIon,
+    .vti = vti,
+    .vte = vte,
+    .RatZeq0 = RatZeq0,
+    .vpar_max_ion = vpar_max_ion,
+    .mu_max_ion = mu_max_ion,
+    .vpar_max_elc = vpar_max_elc,
+    .mu_max_elc = mu_max_elc,
+
+    .Npsi = Npsi,
+    .Nz = Nz,
+    .Nvpar = Nvpar,
+    .Nmu = Nmu,
+    .Nvpar_elc = Nvpar_elc,
+    .Nmu_elc = Nmu_elc,
+    .cells = {Npsi, Nz, Nvpar, Nmu},
+    .Z_min = Z_min,
+    .Z_max = Z_max,
+
+    .write_phase_freq = write_phase_freq,
+    .int_diag_calc_freq = int_diag_calc_freq,
+    .dt_failure_tol = dt_failure_tol,
+    .num_failures_max = num_failures_max,
+
+    .mcB = mcB,
+    .gamma = gamma,
+    .Z_m = Z_m,
+  };
+  
+  // Populate a couple more values in the context.
+  ctx.psi_max = 3e-3;
+  ctx.psi_min = 1e-5;
+  ctx.psi_eval = (ctx.psi_max + ctx.psi_min) / 2.0;
+
+  ctx.z_min    = -2.5;
+  ctx.z_max    = 2.5;
+
+  return ctx;
+}
+
+int main(int argc, char **argv)
+{
+  struct gkyl_app_args app_args = parse_app_args(argc, argv);
+
+#ifdef GKYL_HAVE_MPI
+  if (app_args.use_mpi) MPI_Init(&argc, &argv);
+#endif
+
+  if (app_args.trace_mem) {
+    gkyl_cu_dev_mem_debug_set(true);
+    gkyl_mem_debug_set(true);
+  }
+
+  struct gk_mirror_ctx ctx = create_ctx(); // Context for init functions.
+
+  int rank = 0;
+#ifdef GKYL_HAVE_MPI
+  if (app_args.use_mpi) MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+#endif
+  // if (rank == 0)
+  //   print_ctx(&ctx);
+
+  int cells_x[ctx.cdim], cells_v[ctx.vdim];
+  for (int d=0; d<ctx.cdim; d++)
+    cells_x[d] = APP_ARGS_CHOOSE(app_args.xcells[d], ctx.cells[d]);
+  for (int d=0; d<ctx.vdim; d++)
+    cells_v[d] = APP_ARGS_CHOOSE(app_args.vcells[d], ctx.cells[ctx.cdim+d]);
+
+  // Construct communicator for use in app.
+  struct gkyl_comm *comm = gkyl_gyrokinetic_comms_new(app_args.use_mpi, app_args.use_gpu, stderr);
+
+  struct gkyl_gyrokinetic_species ion = {
+    .name = "ion",
+    .charge = ctx.qi,
+    .mass = ctx.mi,
+    .vdim = ctx.vdim,
+    .lower = {-1.0, 0.0},
+    .upper = { 1.0, 1.0},
+    .cells = { cells_v[0], cells_v[1]},
+    .polarization_density = ctx.n0,
+
+    .projection = {
+      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
+      .density = initial_density,
+      .ctx_density = &ctx,
+      .upar = initial_upar,
+      .ctx_upar = &ctx,
+      .temp = initial_temp_ion,
+      .ctx_temp = &ctx,
+    },
+
+    .mapc2p = {
+      .mapping = mapc2p_vel_ion,
+      .ctx = &ctx,
+    },
+
+    .collisionless = {
+      .type = GKYL_GK_COLLISIONLESS_ES,
+      .write_diagnostics = true,
+    },
+
+    .collisions = {
+      .collision_id = GKYL_LBO_COLLISIONS,
+      .den_ref = ctx.n0,
+      .temp_ref = ctx.Ti0,
+      .num_cross_collisions = 1,
+      .collide_with = { "elc" },
+      .write_diagnostics = true,
+      .write_cross_diagnostics = true,
+    },
+
+    .write_omega_cfl = true,
+    .num_diag_moments = 8,
+    .diag_moments = {GKYL_F_MOMENT_BIMAXWELLIAN, GKYL_F_MOMENT_M0, GKYL_F_MOMENT_M1, GKYL_F_MOMENT_M2, GKYL_F_MOMENT_M2PAR, GKYL_F_MOMENT_M2PERP, GKYL_F_MOMENT_M3PAR, GKYL_F_MOMENT_M3PERP },
+    .num_integrated_diag_moments = 1,
+    .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP },
+    .time_rate_diagnostics = true,
+
+    .boundary_flux_diagnostics = {
+      .num_integrated_diag_moments = 1,
+      .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP},
+    },
+  };
+
+
+  struct gkyl_gyrokinetic_species elc = {
+    .name = "elc",
+    .charge = ctx.qe,
+    .mass = ctx.me,
+    .vdim = ctx.vdim,
+    .lower = {-1.0, 0.0},
+    .upper = { 1.0, 1.0},
+    .cells = { ctx.Nvpar_elc, ctx.Nmu_elc },
+
+    .polarization_density = ctx.n0,
+
+    .mapc2p = {
+      .mapping = mapc2p_vel_elc,
+      .ctx = &ctx,
+    },
+
+    .projection = {
+      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
+      .density = initial_density,
+      .ctx_density = &ctx,
+      .upar = eval_zero,
+      .ctx_upar = &ctx,
+      .temp = initial_temp_elc,
+      .ctx_temp = &ctx,
+      .correct_all_moms = true,
+    },
+
+    .collisions =  {
+      .collision_id = GKYL_LBO_COLLISIONS,
+      .den_ref = ctx.n0,
+      .temp_ref = ctx.Te0,
+      .num_cross_collisions = 1,
+      .collide_with = { "ion" },
+      .write_diagnostics = true,
+      .not_in_dfdt = true,
+    },
+
+    .num_diag_moments = 1,
+    .diag_moments = {GKYL_F_MOMENT_MAXWELLIAN},
+  };
+
+  struct gkyl_gyrokinetic_field field = {
+    .gkfield_id = GKYL_GK_FIELD_BOLTZMANN,
+    .electron_mass = ctx.me,
+    .electron_charge = ctx.qe,
+    .electron_temp = ctx.Te0,
+    .is_static = false,
+  };
+
+  struct gkyl_mirror_geo_grid_inp grid_inp = {
+    .filename_psi = "gyrokinetic/data/eqdsk/lorentzian_R10.geqdsk_psi.gkyl", // psi file to use
+    .rclose = 0.2, // closest R to region of interest
+    .zmin = -2.5,  // Z of lower boundary
+    .zmax =  2.5,  // Z of upper boundary
+    .include_axis = false, // Include R=0 axis in grid
+    .fl_coord = GKYL_GEOMETRY_MIRROR_GRID_GEN_SQRT_PSI_CART_Z, // coordinate system for psi grid
+  };
+
+  struct gkyl_gk app_inp = {  // GK app
+    .name = "gk_lorentzian_mirror",
+    .cdim = ctx.cdim,
+    .upper = {ctx.psi_max, ctx.Z_max},
+    .lower = {ctx.psi_min, ctx.Z_min},
+    .cells = { cells_x[0], cells_x[1] },
+    .poly_order = ctx.poly_order,
+    .basis_type = app_args.basis_type,
+
+    .geometry = {
+      .geometry_id = GKYL_GEOMETRY_MIRROR,
+      .world = {0.0},
+      .mirror_grid_info = grid_inp,
+    },
+
+    .num_periodic_dir = 0,
+    .periodic_dirs = {},
+
+    .num_species = 2,
+    .species = {ion, elc},
+
+    .field = field,
+
+    .parallelism = {
+      .use_gpu = app_args.use_gpu,
+      .cuts = { app_args.cuts[0], app_args.cuts[1] },
+      .comm = comm,
+    },
+  };
+  struct gkyl_gyrokinetic_run_inp run_inp = {
+    .app_inp = app_inp, 
+    .time_stepping = {
+      .t_end = 5e-9,
+      .num_frames = 50,
+      .write_phase_freq = 1,
+      .int_diag_calc_num = 50*100,
+      .dt_failure_tol = ctx.dt_failure_tol,
+      .num_failures_max = ctx.num_failures_max,
+      .is_restart = false,
+      .restart_frame = 0,
+      .num_steps = app_args.num_steps,
+    },
+    .print_verbosity = {
+      .disable_timings = true,
+      .estimate_completion_time = true,
+      .enabled = true,
+      .frequency = 1,
+    },
+  };
+
+  gkyl_gyrokinetic_run_simulation(&run_inp);
+
+  gkyl_gyrokinetic_comms_release(comm);
+  
+#ifdef GKYL_HAVE_MPI
+  if (app_args.use_mpi)
+    MPI_Finalize();
+#endif
+  return 0;
+}

From 4b71ce549ec4e86d9337b1e62fd893b381ea5205 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sun, 17 May 2026 10:38:43 -0400
Subject: [PATCH 07/11] Refactor diagnostics handling by renaming
 write_cross_diagnostics to write_diagnostics

---
 gyrokinetic/apps/gk_species_lbo.c        | 3 +--
 gyrokinetic/apps/gkyl_gyrokinetic.h      | 1 -
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 1 -
 3 files changed, 1 insertion(+), 4 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_lbo.c b/gyrokinetic/apps/gk_species_lbo.c
index 31b8d38c2a..b1f398bbb3 100644
--- a/gyrokinetic/apps/gk_species_lbo.c
+++ b/gyrokinetic/apps/gk_species_lbo.c
@@ -207,7 +207,7 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
   app->stat.n_diag_io += 2;
 
   // Per-cross-species diagnostics: unscaled cross primitive moments and cross-nu.
-  if (gks->lbo.write_cross_diagnostics) {
+  if (gks->lbo.write_diagnostics) {
     // Self primitive moments (u_par_s, vtsq_s) before any cross-species contribution.
     const char *fmt_prim = "%s-%s_lbo_self_prim_moms_%d.gkyl";
     int sz_prim = gkyl_calc_strlen(fmt_prim, app->name, gks->info.name, frame);
@@ -264,7 +264,6 @@ gk_species_lbo_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks, st
 {
   lbo->collision_id = gks->info.collisions.collision_id;
   lbo->write_diagnostics = gks->info.collisions.write_diagnostics;
-  lbo->write_cross_diagnostics = gks->info.collisions.write_cross_diagnostics;
 
   // Empty methods.
   lbo->moms_func = gklbo_moms_disabled;
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index ca6eefb085..77d82317ca 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -106,7 +106,6 @@ struct gkyl_gyrokinetic_collisionless {
 struct gkyl_gyrokinetic_collisions {
   enum gkyl_collision_id collision_id; // type of collisions (see gkyl_eqn_type.h)
   bool write_diagnostics; // Whether to output diagnostics.
-  bool write_cross_diagnostics; // Whether to output per-cross-species diagnostics (prim_moms_sr, nu_sr).
   bool not_in_dfdt; // If true, the collision operator will not be added to df/dt.
     // Used to ignore the collisional updates of this species, while updating cross-species collisions.
 
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index d42bae4a13..bf79444f2e 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -324,7 +324,6 @@ struct gk_collisionless {
 struct gk_lbo_collisions {
   enum gkyl_collision_id collision_id; // type of collisions
   bool write_diagnostics; // Whether to write diagnostics out.
-  bool write_cross_diagnostics; // Whether to write per-cross-species diagnostics (prim_moms_sr, nu_sr).
   bool not_in_dfdt; // Whether to not add collision contribution to df/dt
 
   struct gkyl_array *self_nu; // Self-collision frequency.

From d08714799d889429994238e62ec1afe960d3c641 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 13 Jun 2026 18:32:08 -0700
Subject: [PATCH 08/11] Claude found the bug in 30 seconds that was crashing my
 simulations. This took me hours to debug and claude found it in 30 seconds,
 but I wouldn't have been able to focus claude in the right areas without the
 debugging.

---
 gyrokinetic/apps/gk_species_lbo.c | 1 +
 1 file changed, 1 insertion(+)

diff --git a/gyrokinetic/apps/gk_species_lbo.c b/gyrokinetic/apps/gk_species_lbo.c
index b1f398bbb3..abdf402406 100644
--- a/gyrokinetic/apps/gk_species_lbo.c
+++ b/gyrokinetic/apps/gk_species_lbo.c
@@ -243,6 +243,7 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
         gkyl_array_copy(gks->lbo.nu_prim_moms_host, cross_prim_moms_scratch);
         cross_prim_moms_io = gks->lbo.nu_prim_moms_host;
         gkyl_array_copy(gks->lbo.nu_sum_host, gks->lbo.cross_nu[i]);
+        cross_nu_io = gks->lbo.nu_sum_host;
       }
       else {
         cross_prim_moms_io = cross_prim_moms_scratch;

From 346365a1d09ab367e381dda31ba4800981591f0e Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Wed, 17 Jun 2026 21:04:07 -0700
Subject: [PATCH 09/11] Apply a floor on vtsq cell average as a limiter against
 it going negative. Cell average being zero means that there should be no
 cross species collisions, reducing to the self-species collisions

---
 gyrokinetic/apps/gk_species_lbo.c        | 27 +++++++++++++++++++++++-
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h |  2 ++
 gyrokinetic/creg/rt_gk_mirror_failure.c  |  1 -
 3 files changed, 28 insertions(+), 2 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_lbo.c b/gyrokinetic/apps/gk_species_lbo.c
index c9c38ba2a6..6cfb4e27e7 100644
--- a/gyrokinetic/apps/gk_species_lbo.c
+++ b/gyrokinetic/apps/gk_species_lbo.c
@@ -127,6 +127,13 @@ gklbo_calc_cross_prim_moms(gkyl_gyrokinetic_app *app, const struct gk_species *g
     gks->info.mass, lbo->nu_moms, lbo->prim_moms,
     lbo->other_m[coll_idx], lbo->collide_with[coll_idx]->lbo.moms.marr, lbo->other_prim_moms[coll_idx],
     lbo->nu_boundary_corrections, lbo->cross_nu[coll_idx], lbo->nu_boundary_corrections);
+
+    struct gkyl_array *cross_prim_moms = lbo->nu_boundary_corrections;
+
+    // Floor the cell-average of vtSq_{sr} to zero, setting the rest of the coefficients to zero as well
+    gkyl_array_set_offset_range(lbo->cross_vtsq, 1.0, cross_prim_moms, 1*app->basis.num_basis, &app->local);
+    gkyl_array_max_by_cell_per_cell_avg_range(lbo->cross_vtsq, lbo->cross_vtsq_floor, &app->local);
+    gkyl_array_set_offset_range(cross_prim_moms, 1.0, lbo->cross_vtsq, 1*app->basis.num_basis, &app->local);
 }
 
 static void
@@ -141,6 +148,8 @@ gklbo_cross_moms_enabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks
 
     // Scale upar_{sr} and vtSq_{sr} by nu_{sr} and accumulate.
     struct gkyl_array *cross_prim_moms = lbo->nu_boundary_corrections;
+
+    // Scale upar_{sr} and vtSq_{sr} by nu_{sr}.
     for (int d=0; d<2; d++)
       gkyl_dg_mul_op(app->basis, d, cross_prim_moms, d, cross_prim_moms, 0, lbo->cross_nu[i]);
     gkyl_array_accumulate(lbo->nu_prim_moms, 1.0, cross_prim_moms);
@@ -239,7 +248,15 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
     } else {
       prim_moms_io = gks->lbo.prim_moms;
     }
-    gkyl_comm_array_write(app->comm, &app->grid, &app->local, desc_nu_prim, prim_moms_io, fileNm_prim);
+    struct gkyl_msgpack_map_elem desc_prim[] = {
+      { .key = "Description", .elem_type = GKYL_MP_STRING,
+        .cval = "Self drift velocity and thermal speed squared, before cross-species contributions." }
+    };
+    int io_meta_prim_len[] = {app->io_meta_grid_len, app->gk_geom->io_meta_basic_len, 1};
+    const struct gkyl_msgpack_map_elem* io_meta_prim[] = {app->io_meta_grid, app->gk_geom->io_meta_basic, desc_prim};
+    struct gkyl_msgpack_data *mt_prim = gkyl_msgpack_create_union(sizeof(io_meta_prim_len)/sizeof(int), io_meta_prim_len, io_meta_prim);
+    gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt_prim, prim_moms_io, fileNm_prim);
+    gkyl_msgpack_data_release(mt_prim);
     app->stat.n_diag_io += 1;
     for (int i=0; i<gks->lbo.num_cross_collisions; ++i) {
       const char *other_name = gks->lbo.collide_with[i]->info.name;
@@ -435,6 +452,12 @@ gk_species_lbo_cross_init(struct gkyl_gyrokinetic_app *app, struct gk_species *g
 
       // Morse's alpha_E.
       lbo->alpha_E = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+
+      // Scratch space for flooring the vtSq component of the cross-prim moments.
+      lbo->cross_vtsq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+      lbo->cross_vtsq_floor = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+      gkyl_array_clear(lbo->cross_vtsq_floor, 0.0); // Floor value of zero.
+
       for (int i=0; i<lbo->num_cross_collisions; ++i) {
         // Cross-species collision frequency, nu_sr.
         lbo->cross_nu[i] = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
@@ -571,6 +594,8 @@ gk_species_lbo_release(const struct gkyl_gyrokinetic_app *app, const struct gk_l
         gkyl_array_release(lbo->cross_nu[i]);
 
       gkyl_array_release(lbo->alpha_E);
+      gkyl_array_release(lbo->cross_vtsq);
+      gkyl_array_release(lbo->cross_vtsq_floor);
     }
 
     gkyl_dg_updater_lbo_gyrokinetic_release(lbo->coll_slvr);
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 8731f21136..e62aeac99d 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -355,6 +355,8 @@ struct gk_lbo_collisions {
   struct gkyl_array *cross_nu[GKYL_MAX_SPECIES]; // Cross-species collision frequencies.
   struct gkyl_array *cross_nu_prim_moms; // Weak multiplication of collision frequency and primitive moments.
   struct gkyl_array *alpha_E; // Morse's alpha_E factor.
+  struct gkyl_array *cross_vtsq; // Scratch space holding the vtSq component of the cross-prim moments.
+  struct gkyl_array *cross_vtsq_floor; // Zeroed array used to floor cross vtSq cell-average to >= 0.
   gkyl_prim_lbo_cross_calc *cross_calc; // LBO cross-primitive moment calculator
   
   struct gk_species_moment moms; // Moments needed in LBO (M0, M1, M2).
diff --git a/gyrokinetic/creg/rt_gk_mirror_failure.c b/gyrokinetic/creg/rt_gk_mirror_failure.c
index 586966e308..1de04ededc 100644
--- a/gyrokinetic/creg/rt_gk_mirror_failure.c
+++ b/gyrokinetic/creg/rt_gk_mirror_failure.c
@@ -345,7 +345,6 @@ int main(int argc, char **argv)
       .num_cross_collisions = 1,
       .collide_with = { "elc" },
       .write_diagnostics = true,
-      .write_cross_diagnostics = true,
     },
 
     .write_omega_cfl = true,

From 261e0e5bd7599b94607fa68386f2ba8b23895f34 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Jul 2026 16:00:35 -0700
Subject: [PATCH 10/11] Remove rt_gk_mirror_failure.c file

---
 gyrokinetic/creg/rt_gk_mirror_failure.c | 481 ------------------------
 1 file changed, 481 deletions(-)
 delete mode 100644 gyrokinetic/creg/rt_gk_mirror_failure.c

diff --git a/gyrokinetic/creg/rt_gk_mirror_failure.c b/gyrokinetic/creg/rt_gk_mirror_failure.c
deleted file mode 100644
index 1de04ededc..0000000000
--- a/gyrokinetic/creg/rt_gk_mirror_failure.c
+++ /dev/null
@@ -1,481 +0,0 @@
-#include <math.h>
-#include <stdio.h>
-#include <time.h>
-
-#include <gkyl_alloc.h>
-#include <gkyl_const.h>
-#include <gkyl_eqn_type.h>
-#include <gkyl_fem_poisson_bctype.h>
-#include <gkyl_gyrokinetic.h>
-#include <gkyl_gyrokinetic_run.h>
-#include <gkyl_math.h>
-
-#include <rt_arg_parse.h>
-
-// Define the context of the simulation. This is basically all the globals
-struct gk_mirror_ctx
-{
-  int cdim, vdim; // Dimensionality.
-  // Plasma parameters
-  double mi;
-  double qi;
-  double me;
-  double qe;
-  double Te0;
-  double n0;
-  double B_p; // Hardcoded magnetic field at midplane 0.53 for double lorentzian
-  double Bmag_midp; // Magnetic field magnitude at midplane (Z=0)
-  double beta;
-  double tau;
-  double Ti0;
-  double nuFrac;
-  // Ion-ion collision freq.
-  double logLambdaIon;
-  double nuIon;
-  double vti, vte;
-  double RatZeq0; // Radius of the field line at Z=0.
-  double kperp; // Perpendicular wavenumber for ES gyrokinetics.
-  // Axial coordinate Z extents. Endure that Z=0 is not on
-  double z_min;
-  double z_max;
-  double psi_eval;
-  double psi_max;
-  double psi_min;
-  double theta_eval;
-  double theta_min;
-  double theta_max;
-  // Physics parameters at mirror throat
-  double vpar_max_ion;
-  double mu_max_ion;
-  double vpar_max_elc;
-  double mu_max_elc;
-
-  int Npsi;
-  int Ntheta;
-  int Nz;
-  int Nvpar;
-  int Nmu;
-  int Nvpar_elc;
-  int Nmu_elc;
-  int cells[GKYL_MAX_DIM]; // Number of cells in all directions.
-  int poly_order;
-
-  // Source parameters
-  double ion_source_amplitude;
-  double ion_source_sigma;
-  double ion_source_temp;
-
-  double t_end; // End time.
-  int num_frames; // Number of output frames.
-  int num_phases; // Number of phases.
-  double write_phase_freq; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
-  double int_diag_calc_freq; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
-  double dt_failure_tol; // Minimum allowable fraction of initial time-step.
-  int num_failures_max; // Maximum allowable number of consecutive small time-steps.
-
-  // Geometry parameters for Lorentzian mirror
-  double mcB;     // Magnetic field parameter
-  double mcB_inner; // Magnetic field parameter for inner mirror
-  double mcB_outer; // Magnetic field parameter for outer mirror
-  double gamma;   // Width parameter for Lorentzian profile
-  double gamma_inner; // Width parameter for Lorentzian profile for inner mirror
-  double gamma_outer; // Width parameter for Lorentzian profile for outer mirror
-  double Z_m;     // Mirror throat location
-  double L_center; // Length of central region
-  double L_plugs; // Length of plugs
-  double Z_min;   // Minimum Z coordinate
-  double Z_max;   // Maximum Z coordinate
-  double psi_in;  // Working variable for psi integration
-  double z_in;    // Working variable for z integration
-};
-
-
-// Evaluate initial conditions
-// I think ICs should be constructed to match the expander rather than the center
-// as the center is quickly filled in the OAP
-void
-initial_density(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  fout[0] = 1e17;
-}
-
-void
-initial_upar(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  double z = xn[1];
-  double c_s = 7 * sqrt(app->Te0/app->mi);
-  if (fabs(z) <= app->Z_m)
-  {
-    fout[0] = 0.0;
-  }
-  else
-  {
-    fout[0] = fabs(z) / z * c_s * tanh(4 * (app->Z_max - app->Z_m) * fabs(fabs(z) - app->Z_m));
-  }
-}
-
-void
-eval_zero(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  fout[0] = 0.0;
-}
-
-void
-initial_temp_ion(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  fout[0] = app->Ti0/10.0;
-}
-
-void
-initial_temp_elc(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  fout[0] = app->Te0;
-}
-
-void mapc2p_vel_ion(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  double vpar_max_ion = app->vpar_max_ion;
-  double mu_max_ion = app->mu_max_ion;
-
-  double cvpar = vc[0], cmu = vc[1];
-  double b = 1.4;
-  vp[0] = vpar_max_ion*tan(cvpar*b)/tan(b);
-  vp[1] = mu_max_ion*pow(cmu,2);  // Cubic map in mu.
-}
-
-void mapc2p_vel_elc(double t, const double *vc, double* GKYL_RESTRICT vp, void *ctx)
-{
-  struct gk_mirror_ctx *app = ctx;
-  double mu_max_elc = app->mu_max_elc;
-  double vpar_max_elc = app->vpar_max_elc;
-  double cvpar = vc[0], cmu = vc[1];
-  vp[0] = vpar_max_elc*cvpar;
-  vp[1] = mu_max_elc*pow(cmu,4);  // Cubic map in mu.
-}
-
-struct gk_mirror_ctx
-create_ctx(void)
-{
-  int cdim = 2, vdim = 2; // Dimensionality.
-  int poly_order = 1;
-
-  // Universal constant parameters.
-  double eps0 = GKYL_EPSILON0;
-  double mu0 = GKYL_MU0; // Not sure if this is right
-  double eV = GKYL_ELEMENTARY_CHARGE;
-  double mp = GKYL_PROTON_MASS; // ion mass
-  double me = GKYL_ELECTRON_MASS;
-  double qi = eV;  // ion charge
-  double qe = -eV; // electron charge
-
-  // Plasma parameters.
-  double mi = 2.014 * mp;
-  double Te0 = 940 * eV;
-  double n0 = 3e19;
-  double B_p = 0.53; // Bmag at z=0
-  double beta = 0.4;
-  double tau = pow(B_p, 2.) * beta / (2.0 * mu0 * n0 * Te0) - 1.;
-  double Ti0 = tau * Te0;
-
-  // Ion-ion collision freq.
-  double nuFrac = 1.0;
-  double logLambdaIon = 6.6 - 0.5 * log(n0 / 1e20) + 1.5 * log(Ti0 / eV);
-  double nuIon = nuFrac * logLambdaIon * pow(eV, 4.) * n0 /
-                 (12 * pow(M_PI, 3. / 2.) * pow(eps0, 2.) * sqrt(mi) * pow(Ti0, 3. / 2.));
-
-  // Thermal speeds.
-  double vti = sqrt(Ti0 / mi);
-  double vte = sqrt(Te0 / me);
-
-  // Grid parameters
-  double vpar_max_ion = 16 * vti;
-  double mu_max_ion = mi * pow(3. * vti, 2.) / (2. * B_p);
-  double vpar_max_elc = 4 * vte;
-  double mu_max_elc = me * pow(4. * vte, 2.) / (2. * B_p);
-
-  int Nz = 32;
-  int Npsi = 2;
-  int Nvpar = 16;
-  int Nmu = 8;
-  int Nvpar_elc = 8;
-  int Nmu_elc = 8;
-
-  // Geometry parameters.
-  double RatZeq0 = 0.10; // Radius of the field line at Z=0.
-  double Z_min = -2.5;
-  double Z_max =  2.5;
-  double mcB = 3.691260;
-  double gamma = 0.226381;
-  double Z_m = 0.98;
-
-  // Calculate phase structure
-  double write_phase_freq = 1; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
-  double int_diag_calc_freq = 100; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
-  double dt_failure_tol = 1.0e-4; // Minimum allowable fraction of initial time-step.
-  int num_failures_max = 20; // Maximum allowable number of consecutive small time-steps.
-
-  struct gk_mirror_ctx ctx = {
-    .cdim = cdim,
-    .vdim = vdim,
-    .poly_order = poly_order,
-
-    .mi = mi,
-    .qi = qi,
-    .me = me,
-    .qe = qe,
-    .Te0 = Te0,
-    .n0 = n0,
-    .B_p = B_p,
-    .beta = beta,
-    .tau = tau,
-    .Ti0 = Ti0,
-
-    .nuFrac = nuFrac,
-    .logLambdaIon = logLambdaIon,
-    .nuIon = nuIon,
-    .vti = vti,
-    .vte = vte,
-    .RatZeq0 = RatZeq0,
-    .vpar_max_ion = vpar_max_ion,
-    .mu_max_ion = mu_max_ion,
-    .vpar_max_elc = vpar_max_elc,
-    .mu_max_elc = mu_max_elc,
-
-    .Npsi = Npsi,
-    .Nz = Nz,
-    .Nvpar = Nvpar,
-    .Nmu = Nmu,
-    .Nvpar_elc = Nvpar_elc,
-    .Nmu_elc = Nmu_elc,
-    .cells = {Npsi, Nz, Nvpar, Nmu},
-    .Z_min = Z_min,
-    .Z_max = Z_max,
-
-    .write_phase_freq = write_phase_freq,
-    .int_diag_calc_freq = int_diag_calc_freq,
-    .dt_failure_tol = dt_failure_tol,
-    .num_failures_max = num_failures_max,
-
-    .mcB = mcB,
-    .gamma = gamma,
-    .Z_m = Z_m,
-  };
-  
-  // Populate a couple more values in the context.
-  ctx.psi_max = 3e-3;
-  ctx.psi_min = 1e-5;
-  ctx.psi_eval = (ctx.psi_max + ctx.psi_min) / 2.0;
-
-  ctx.z_min    = -2.5;
-  ctx.z_max    = 2.5;
-
-  return ctx;
-}
-
-int main(int argc, char **argv)
-{
-  struct gkyl_app_args app_args = parse_app_args(argc, argv);
-
-#ifdef GKYL_HAVE_MPI
-  if (app_args.use_mpi) MPI_Init(&argc, &argv);
-#endif
-
-  if (app_args.trace_mem) {
-    gkyl_cu_dev_mem_debug_set(true);
-    gkyl_mem_debug_set(true);
-  }
-
-  struct gk_mirror_ctx ctx = create_ctx(); // Context for init functions.
-
-  int rank = 0;
-#ifdef GKYL_HAVE_MPI
-  if (app_args.use_mpi) MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-#endif
-  // if (rank == 0)
-  //   print_ctx(&ctx);
-
-  int cells_x[ctx.cdim], cells_v[ctx.vdim];
-  for (int d=0; d<ctx.cdim; d++)
-    cells_x[d] = APP_ARGS_CHOOSE(app_args.xcells[d], ctx.cells[d]);
-  for (int d=0; d<ctx.vdim; d++)
-    cells_v[d] = APP_ARGS_CHOOSE(app_args.vcells[d], ctx.cells[ctx.cdim+d]);
-
-  // Construct communicator for use in app.
-  struct gkyl_comm *comm = gkyl_gyrokinetic_comms_new(app_args.use_mpi, app_args.use_gpu, stderr);
-
-  struct gkyl_gyrokinetic_species ion = {
-    .name = "ion",
-    .charge = ctx.qi,
-    .mass = ctx.mi,
-    .vdim = ctx.vdim,
-    .lower = {-1.0, 0.0},
-    .upper = { 1.0, 1.0},
-    .cells = { cells_v[0], cells_v[1]},
-    .polarization_density = ctx.n0,
-
-    .projection = {
-      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
-      .density = initial_density,
-      .ctx_density = &ctx,
-      .upar = initial_upar,
-      .ctx_upar = &ctx,
-      .temp = initial_temp_ion,
-      .ctx_temp = &ctx,
-    },
-
-    .mapc2p = {
-      .mapping = mapc2p_vel_ion,
-      .ctx = &ctx,
-    },
-
-    .collisionless = {
-      .type = GKYL_GK_COLLISIONLESS_ES,
-      .write_diagnostics = true,
-    },
-
-    .collisions = {
-      .collision_id = GKYL_LBO_COLLISIONS,
-      .den_ref = ctx.n0,
-      .temp_ref = ctx.Ti0,
-      .num_cross_collisions = 1,
-      .collide_with = { "elc" },
-      .write_diagnostics = true,
-    },
-
-    .write_omega_cfl = true,
-    .num_diag_moments = 8,
-    .diag_moments = {GKYL_F_MOMENT_BIMAXWELLIAN, GKYL_F_MOMENT_M0, GKYL_F_MOMENT_M1, GKYL_F_MOMENT_M2, GKYL_F_MOMENT_M2PAR, GKYL_F_MOMENT_M2PERP, GKYL_F_MOMENT_M3PAR, GKYL_F_MOMENT_M3PERP },
-    .num_integrated_diag_moments = 1,
-    .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP },
-    .time_rate_diagnostics = true,
-
-    .boundary_flux_diagnostics = {
-      .num_integrated_diag_moments = 1,
-      .integrated_diag_moments = { GKYL_F_MOMENT_M0M1M2PARM2PERP},
-    },
-  };
-
-
-  struct gkyl_gyrokinetic_species elc = {
-    .name = "elc",
-    .charge = ctx.qe,
-    .mass = ctx.me,
-    .vdim = ctx.vdim,
-    .lower = {-1.0, 0.0},
-    .upper = { 1.0, 1.0},
-    .cells = { ctx.Nvpar_elc, ctx.Nmu_elc },
-
-    .polarization_density = ctx.n0,
-
-    .mapc2p = {
-      .mapping = mapc2p_vel_elc,
-      .ctx = &ctx,
-    },
-
-    .projection = {
-      .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM,
-      .density = initial_density,
-      .ctx_density = &ctx,
-      .upar = eval_zero,
-      .ctx_upar = &ctx,
-      .temp = initial_temp_elc,
-      .ctx_temp = &ctx,
-      .correct_all_moms = true,
-    },
-
-    .collisions =  {
-      .collision_id = GKYL_LBO_COLLISIONS,
-      .den_ref = ctx.n0,
-      .temp_ref = ctx.Te0,
-      .num_cross_collisions = 1,
-      .collide_with = { "ion" },
-      .write_diagnostics = true,
-      .not_in_dfdt = true,
-    },
-
-    .num_diag_moments = 1,
-    .diag_moments = {GKYL_F_MOMENT_MAXWELLIAN},
-  };
-
-  struct gkyl_gyrokinetic_field field = {
-    .gkfield_id = GKYL_GK_FIELD_BOLTZMANN,
-    .electron_mass = ctx.me,
-    .electron_charge = ctx.qe,
-    .electron_temp = ctx.Te0,
-    .is_static = false,
-  };
-
-  struct gkyl_mirror_geo_grid_inp grid_inp = {
-    .filename_psi = "gyrokinetic/data/eqdsk/lorentzian_R10.geqdsk_psi.gkyl", // psi file to use
-    .rclose = 0.2, // closest R to region of interest
-    .zmin = -2.5,  // Z of lower boundary
-    .zmax =  2.5,  // Z of upper boundary
-    .include_axis = false, // Include R=0 axis in grid
-    .fl_coord = GKYL_GEOMETRY_MIRROR_GRID_GEN_SQRT_PSI_CART_Z, // coordinate system for psi grid
-  };
-
-  struct gkyl_gk app_inp = {  // GK app
-    .name = "gk_lorentzian_mirror",
-    .cdim = ctx.cdim,
-    .upper = {ctx.psi_max, ctx.Z_max},
-    .lower = {ctx.psi_min, ctx.Z_min},
-    .cells = { cells_x[0], cells_x[1] },
-    .poly_order = ctx.poly_order,
-    .basis_type = app_args.basis_type,
-
-    .geometry = {
-      .geometry_id = GKYL_GEOMETRY_MIRROR,
-      .world = {0.0},
-      .mirror_grid_info = grid_inp,
-    },
-
-    .num_periodic_dir = 0,
-    .periodic_dirs = {},
-
-    .num_species = 2,
-    .species = {ion, elc},
-
-    .field = field,
-
-    .parallelism = {
-      .use_gpu = app_args.use_gpu,
-      .cuts = { app_args.cuts[0], app_args.cuts[1] },
-      .comm = comm,
-    },
-  };
-  struct gkyl_gyrokinetic_run_inp run_inp = {
-    .app_inp = app_inp, 
-    .time_stepping = {
-      .t_end = 5e-9,
-      .num_frames = 50,
-      .write_phase_freq = 1,
-      .int_diag_calc_num = 50*100,
-      .dt_failure_tol = ctx.dt_failure_tol,
-      .num_failures_max = ctx.num_failures_max,
-      .is_restart = false,
-      .restart_frame = 0,
-      .num_steps = app_args.num_steps,
-    },
-    .print_verbosity = {
-      .disable_timings = true,
-      .estimate_completion_time = true,
-      .enabled = true,
-      .frequency = 1,
-    },
-  };
-
-  gkyl_gyrokinetic_run_simulation(&run_inp);
-
-  gkyl_gyrokinetic_comms_release(comm);
-  
-#ifdef GKYL_HAVE_MPI
-  if (app_args.use_mpi)
-    MPI_Finalize();
-#endif
-  return 0;
-}

From 04df9e9486f14c7bc5d9d2ce30e7c762eb11459c Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Jul 2026 16:36:02 -0700
Subject: [PATCH 11/11] Fix indexing for io_meta_dg in gklbo_write_mom_enabled
 function

---
 gyrokinetic/apps/gk_species_lbo.c | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_lbo.c b/gyrokinetic/apps/gk_species_lbo.c
index 1edd89f9aa..ee509ec02b 100644
--- a/gyrokinetic/apps/gk_species_lbo.c
+++ b/gyrokinetic/apps/gk_species_lbo.c
@@ -152,7 +152,7 @@ gklbo_cross_moms_enabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks
     // Scale upar_{sr} and vtSq_{sr} by nu_{sr}.
     for (int d=0; d<2; d++)
       gkyl_dg_mul_op(&app->basis, d, cross_prim_moms, d, cross_prim_moms, 0, lbo->cross_nu[i]);
-      
+
     gkyl_array_accumulate(lbo->nu_prim_moms, 1.0, cross_prim_moms);
   }
   app->stat.species_coll_mom_tm += gkyl_time_diff_now_sec(wst);
@@ -210,8 +210,8 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
     { .key = "Description", .elem_type = GKYL_MP_STRING,
       .cval = "Cross-species drift velocity and thermal speed squared, or cross-species collision frequency." }
   };
-  int io_meta_nu_cross_len[] = {app->io_meta_grid_len, app->gk_geom->io_meta_basic_len, 1};
-  const struct gkyl_msgpack_map_elem* io_meta_nu_cross[] = {app->io_meta_grid, app->gk_geom->io_meta_basic, desc_nu_cross};
+  int io_meta_nu_cross_len[] = {app->io_meta_dg_len, app->gk_geom->io_meta_basic_len, 1};
+  const struct gkyl_msgpack_map_elem* io_meta_nu_cross[] = {app->io_meta_dg, app->gk_geom->io_meta_basic, desc_nu_cross};
   struct gkyl_msgpack_data *mt_nu_cross = gkyl_msgpack_create_union(sizeof(io_meta_nu_cross_len)/sizeof(int), io_meta_nu_cross_len, io_meta_nu_cross);
 
   // Write out nu_sum and nu_prim_moms.
@@ -253,8 +253,8 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
       { .key = "Description", .elem_type = GKYL_MP_STRING,
         .cval = "Self drift velocity and thermal speed squared, before cross-species contributions." }
     };
-    int io_meta_prim_len[] = {app->io_meta_grid_len, app->gk_geom->io_meta_basic_len, 1};
-    const struct gkyl_msgpack_map_elem* io_meta_prim[] = {app->io_meta_grid, app->gk_geom->io_meta_basic, desc_prim};
+    int io_meta_prim_len[] = {app->io_meta_dg_len, app->gk_geom->io_meta_basic_len, 1};
+    const struct gkyl_msgpack_map_elem* io_meta_prim[] = {app->io_meta_dg, app->gk_geom->io_meta_basic, desc_prim};
     struct gkyl_msgpack_data *mt_prim = gkyl_msgpack_create_union(sizeof(io_meta_prim_len)/sizeof(int), io_meta_prim_len, io_meta_prim);
     gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt_prim, prim_moms_io, fileNm_prim);
     gkyl_msgpack_data_release(mt_prim);