From 04d0050b25a6ca6668eb7b57c901bc3254ed91e6 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Fri, 27 Mar 2026 19:03:30 -0400
Subject: [PATCH 01/25] Implement low-pass filter damping in gyrokinetic
 species model

- Add low-pass filter damping type to gkyl_gyrokinetic_damping enum.
- Initialize filtered distribution function in gk_species_apply_ic.
- Implement damping initialization and calculation functions for low-pass filter.
- Update SSP RK3 integration to advance filtered distributions.
---
 gyrokinetic/apps/gk_species.c                 |  3 +
 gyrokinetic/apps/gk_species_damping.c         | 62 +++++++++++++++++++
 gyrokinetic/apps/gkyl_gyrokinetic.h           |  1 +
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h      | 28 +++++++++
 gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c | 42 +++++++++++++
 5 files changed, 136 insertions(+)

diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c
index 3345748d97..f6c454c34e 100644
--- a/gyrokinetic/apps/gk_species.c
+++ b/gyrokinetic/apps/gk_species.c
@@ -1732,6 +1732,9 @@ gk_species_apply_ic(gkyl_gyrokinetic_app *app, struct gk_species *gks, double t0
   if (gks->info.init_from_file.type == 0)
     gk_species_projection_calc(app, gks, &gks->proj_init, gks->f, t0);
 
+  // Initialize fbar for low-pass filter damping
+  gk_species_damping_init_fbar_from_f(gks, &gks->damping, gks->f);
+
   // We are pre-computing source for now as it is time-independent.
   gk_species_source_calc(app, gks, &gks->src, gks->lte.f_lte, t0);
 }
diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index cbc907ddfd..b1855234ea 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -4,6 +4,18 @@
 #include <gkyl_alloc.h>
 #include <gkyl_dg_basis_ops.h>
 
+void
+gk_species_damping_init_fbar_from_f(const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *f)
+{
+  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+    // Initialize the filtered distribution to match the initial distribution
+    gkyl_array_set(damp->fbar, 1.0, f);
+    gkyl_array_set(damp->fbar1, 1.0, f);
+    gkyl_array_set(damp->fbarnew, 1.0, f);
+  }
+}
+
 void
 gk_species_damping_write_disabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
 {
@@ -194,6 +206,19 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       // Multiply by the user's scaling profile.
       gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
     }
+    else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+      // Allocate filtered distribution function array
+      damp->fbar = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
+      damp->fbar1 = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
+      damp->fbarnew = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
+      damp->fbar_host = damp->fbar;
+      if (app->use_gpu)
+        damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
+
+      // Initialize fbar from the projection of the initial distribution
+      // (will be set from the initial f in the main app loop)
+      gkyl_array_set(damp->fbar, 0.0);
+    }
 
     // Set function pointers chosen at runtime.
     if (damp->evolve) {
@@ -205,6 +230,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   }
 }
 
+
 void
 gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
   const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
@@ -234,6 +260,22 @@ gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *g
       gkyl_array_accumulate(rhs, -1.0, f_buffer);
 
     }
+    else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+      // Match fbar state to the RK stage of fin.
+      const struct gkyl_array *fbar_in = damp->fbar;
+      if (fin == gks->f1)
+        fbar_in = damp->fbar1;
+      else if (fin == gks->fnew)
+        fbar_in = damp->fbarnew;
+
+      // Compute f - fbar and scale by the damping rate: rate * (f - fbar)
+      gkyl_array_set(f_buffer, 1.0, fin);      // f_buffer = f
+      gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = f - fbar
+      gkyl_array_scale_by_cell(f_buffer, damp->rate);    // f_buffer = rate * (f - fbar)
+      
+      // Add damping term to RHS: df/dt -= rate * (f - fbar)
+      gkyl_array_accumulate(rhs, -1.0, f_buffer);
+    }
 
     // Add the frequency to the CFL frequency.
     gkyl_array_accumulate(cflrate, 1.0, damp->rate);
@@ -242,12 +284,25 @@ gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *g
   }
 }
 
+
 void
 gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
 {
   gks->damping.write_func(app, gks, tm, frame);
 }
 
+void
+gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+    // rhs_fbar = rate * (f - fbar)
+    gkyl_array_set(rhs_fbar, 1.0, fin);
+    gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
+    gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
+  }
+}
+
 void
 gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp)
 {
@@ -275,5 +330,12 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
       gkyl_loss_cone_mask_gyrokinetic_release(damp->lcm_proj_op);
       gkyl_array_release(damp->scale_prof);
     }
+    else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+      gkyl_array_release(damp->fbar);
+      gkyl_array_release(damp->fbar1);
+      gkyl_array_release(damp->fbarnew);
+      if (app->use_gpu)
+        gkyl_array_release(damp->fbar_host);
+    }
   }
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 0a890ef220..bc846ea381 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -318,6 +318,7 @@ enum gkyl_gyrokinetic_damping_type {
   GKYL_GK_DAMPING_NONE = 0,
   GKYL_GK_DAMPING_USER_INPUT,
   GKYL_GK_DAMPING_LOSS_CONE,
+  GKYL_GK_DAMPING_LOW_PASS_FILTER,
 };
 
 struct gkyl_gyrokinetic_damping {
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index ae971282b2..e8d62ae9c1 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -808,6 +808,10 @@ struct gk_damping {
   double *bmag_max_coord; // Location of bmag_max.
   double *phi_m, *phi_m_global; // Electrostatic potential at bmag_max.
   struct gkyl_array *scale_prof; // Conf-space scaling factor profile.
+  // Low-pass filter variables
+  struct gkyl_array *fbar; // Filtered/averaged distribution function.
+  struct gkyl_array *fbar1, *fbarnew; // SSPRK3 stage arrays for filtered distribution.
+  struct gkyl_array *fbar_host; // Host copy of fbar for use in IO.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
 };
@@ -2771,6 +2775,18 @@ void gk_species_source_release(const struct gkyl_gyrokinetic_app *app, const str
  */
 void gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks, struct gk_damping *damp);
 
+/**
+ * Initialize filtered distribution function from initial distribution.
+ * Should be called after initial conditions are applied to f.
+ *
+ * @param gks Pointer to species.
+ * @param damp Species damping object.
+ * @param f Initial distribution function.
+ */
+void gk_species_damping_init_fbar_from_f(const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *f);
+
+
 /**
  * Compute species applied source term.
  *
@@ -2797,6 +2813,18 @@ void gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_speci
  */
 void gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
 
+/**
+ * Compute filtered distribution RHS for low-pass filter damping:
+ *   d(fbar)/dt = rate * (f - fbar).
+ *
+ * @param damp Species damping object.
+ * @param fin Current distribution function f.
+ * @param fbar_in Current filtered distribution function fbar.
+ * @param rhs_fbar RHS buffer for fbar.
+ */
+void gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
+
 /**
  * Release species damping object.
  *
diff --git a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
index 5a76659633..d6befeb874 100644
--- a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
+++ b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
@@ -1,5 +1,17 @@
 #include <gkyl_gyrokinetic_priv.h>
 
+static void
+gk_species_damping_step_fbar_ssp_rk3(struct gk_species *gks, const struct gkyl_array *fin,
+  const struct gkyl_array *fbar_in, struct gkyl_array *fbar_out, double dt)
+{
+  if (gks->damping.type != GKYL_GK_DAMPING_LOW_PASS_FILTER)
+    return;
+
+  // Reuse existing species scratch buffer for rhs_fbar.
+  gk_species_damping_calc_fbar_rhs(&gks->damping, fin, fbar_in, gks->lte.f_lte);
+  gk_species_step_f(gks, fbar_out, dt, fbar_in);
+}
+
 static void
 gyrokinetic_forward_euler(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   const struct gkyl_array *fin[], struct gkyl_array *fout[], 
@@ -90,6 +102,13 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
         dt = st.dt_actual;
 
+        // Advance filtered distributions with the same stage-1 FE update.
+        for (int i=0; i<app->num_species; ++i) {
+          struct gk_species *gks = &app->species[i];
+          gk_species_damping_step_fbar_ssp_rk3(gks, gks->f,
+            gks->damping.fbar, gks->damping.fbar1, dt);
+        }
+
         // Subtract boundary flux f from f1 so that we only step boundary
         // fluxes during a given time step, not over all time. And so that the
         // boundary flux in f is kept in case a later RK stage fails.
@@ -141,6 +160,13 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
         gyrokinetic_forward_euler(app, tcurr+dt, dt, fin, fout, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
+        // Advance filtered distributions with the same stage-2 FE update.
+        for (int i=0; i<app->num_species; ++i) {
+          struct gk_species *gks = &app->species[i];
+          gk_species_damping_step_fbar_ssp_rk3(gks, gks->f1,
+            gks->damping.fbar1, gks->damping.fbarnew, st.dt_actual);
+        }
+
         if (st.dt_actual < dt) {
 
           // Recalculate the field.
@@ -168,6 +194,10 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           for (int i=0; i<app->num_species; ++i) {
             struct gk_species *gks = &app->species[i];
             gk_species_combine(gks, gks->f1, 3.0/4.0, gks->f, 1.0/4.0, gks->fnew, &gks->local_ext);
+            if (gks->damping.type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+              gk_species_combine(gks, gks->damping.fbar1, 3.0/4.0, gks->damping.fbar,
+                1.0/4.0, gks->damping.fbarnew, &gks->local_ext);
+            }
             gk_species_bflux_set(app, &gks->bflux, gks->bflux.f1, 1.0/4.0, gks->bflux.fnew);
           }
           for (int i=0; i<app->num_neut_species; ++i) {
@@ -213,6 +243,13 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
         gyrokinetic_forward_euler(app, tcurr+dt/2, dt, fin, fout, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
+        // Advance filtered distributions with the same stage-3 FE update.
+        for (int i=0; i<app->num_species; ++i) {
+          struct gk_species *gks = &app->species[i];
+          gk_species_damping_step_fbar_ssp_rk3(gks, gks->f1,
+            gks->damping.fbar1, gks->damping.fbarnew, st.dt_actual);
+        }
+
         if (st.dt_actual < dt) {
           // Recalculate the field.
           for (int i=0; i<app->num_species; ++i) {
@@ -242,6 +279,11 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
             // Step f.
             gk_species_combine(gks, gks->f1, 1.0/3.0, gks->f, 2.0/3.0, gks->fnew, &gks->local_ext);
             gk_species_copy_range(gks, gks->f, gks->f1, &gks->local_ext);
+            if (gks->damping.type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+              gk_species_combine(gks, gks->damping.fbar1, 1.0/3.0, gks->damping.fbar,
+                2.0/3.0, gks->damping.fbarnew, &gks->local_ext);
+              gk_species_copy_range(gks, gks->damping.fbar, gks->damping.fbar1, &gks->local_ext);
+            }
             // Step boundary fluxes.
             gk_species_bflux_set(app, &gks->bflux, gks->bflux.f, 2.0/3.0, gks->bflux.fnew);
             gk_species_bflux_calc_voltime_integrated_mom(app, gks, &gks->bflux, tcurr);

From b99a98cca981eef8700d74ca15791e9bbc58c1a6 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Fri, 27 Mar 2026 19:05:37 -0400
Subject: [PATCH 02/25] Add reset function for species damping in gyrokinetic
 app

---
 gyrokinetic/apps/gk_species_damping.c    | 13 +++++++++++++
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 12 ++++++++++++
 2 files changed, 25 insertions(+)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index b1855234ea..1eac356900 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -339,3 +339,16 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
     }
   }
 }
+
+void
+gk_species_damping_reset(gkyl_gyrokinetic_app* app, double tm, struct gk_species *gks,
+  struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp)
+{
+  (void) tm;
+
+  gk_species_damping_release(app, damp);
+
+  gks->info.damping = damp_inp;
+  gk_species_damping_init(app, gks, damp);
+  gk_species_damping_init_fbar_from_f(gks, damp, gks->f);
+}
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index e8d62ae9c1..c95238d3d3 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -2833,6 +2833,18 @@ void gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
  */
 void gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp);
 
+/**
+ * Reset species damping object.
+ *
+ * @param app gyrokinetic app object.
+ * @param tm Time.
+ * @param gks Species object.
+ * @param damp Species damping object.
+ * @param damp_inp New damping input.
+ */
+void gk_species_damping_reset(gkyl_gyrokinetic_app* app, double tm, struct gk_species *gks,
+  struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp);
+
 /** gk_species_fdot_multiplier API */
 
 /**

From 41e4a2726855f4148718bc6d40cb3a6abf78bb62 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Fri, 27 Mar 2026 19:15:55 -0400
Subject: [PATCH 03/25] Refactor damping functions for low-pass filter
 implementation in gyrokinetic app

---
 gyrokinetic/apps/gk_species.c            |   2 +-
 gyrokinetic/apps/gk_species_damping.c    | 197 ++++++++++++++++-------
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h |   9 +-
 3 files changed, 146 insertions(+), 62 deletions(-)

diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c
index f6c454c34e..d8c042ff98 100644
--- a/gyrokinetic/apps/gk_species.c
+++ b/gyrokinetic/apps/gk_species.c
@@ -1733,7 +1733,7 @@ gk_species_apply_ic(gkyl_gyrokinetic_app *app, struct gk_species *gks, double t0
     gk_species_projection_calc(app, gks, &gks->proj_init, gks->f, t0);
 
   // Initialize fbar for low-pass filter damping
-  gk_species_damping_init_fbar_from_f(gks, &gks->damping, gks->f);
+  gk_species_damping_set_fbar_to_f(gks, &gks->damping, gks->f);
 
   // We are pre-computing source for now as it is time-independent.
   gk_species_source_calc(app, gks, &gks->src, gks->lte.f_lte, t0);
diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 1eac356900..1e8bf3cb58 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -4,16 +4,44 @@
 #include <gkyl_alloc.h>
 #include <gkyl_dg_basis_ops.h>
 
+void gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
+  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
+void gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
+  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
+void gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
+  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
+void gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
+  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
+void gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
+void gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
+
+
 void
-gk_species_damping_init_fbar_from_f(const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk_damping *damp, 
   const struct gkyl_array *f)
 {
-  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-    // Initialize the filtered distribution to match the initial distribution
-    gkyl_array_set(damp->fbar, 1.0, f);
-    gkyl_array_set(damp->fbar1, 1.0, f);
-    gkyl_array_set(damp->fbarnew, 1.0, f);
-  }
+  gkyl_array_set(damp->fbar, 1.0, f);
+  gkyl_array_set(damp->fbar1, 1.0, f);
+  gkyl_array_set(damp->fbarnew, 1.0, f);
+}
+
+void
+gk_species_damping_set_fbar_to_f_disabled(const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *f)
+{
+}
+
+void
+gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *f)
+{
+  damp->set_fbar_to_f_func(gks, damp, f);
 }
 
 void
@@ -83,6 +111,9 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
   // Default function pointers.
   damp->write_func = gk_species_damping_write_disabled;
+  damp->advance_func = gk_species_damping_advance_disabled;
+  damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_disabled;
+  damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_disabled;
 
   if (damp->type) {
     // Allocate rate array.
@@ -113,6 +144,8 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
       if (num_quad == 1)
         gkyl_array_scale_range(damp->rate, 1.0/pow(sqrt(2.0),gks->grid.ndim), &gks->local);
+
+      damp->advance_func = gk_species_damping_advance_user_input;
     }
     else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
       damp->evolve = true; // Since the loss cone boundary is proportional to phi(t).
@@ -139,14 +172,14 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       if (app->use_gpu) {
         damp->bmag_max = gkyl_cu_malloc(sizeof(double));
         damp->bmag_max_coord = gkyl_cu_malloc(app->cdim*sizeof(double));
-	gkyl_cu_memcpy(damp->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
-	gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(damp->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
       }
       else {
         damp->bmag_max = gkyl_malloc(sizeof(double));
         damp->bmag_max_coord = gkyl_malloc(app->cdim*sizeof(double));
-	memcpy(damp->bmag_max, &bmag_max_global, sizeof(double));
-	memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double));
+        memcpy(damp->bmag_max, &bmag_max_global, sizeof(double));
+        memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double));
       }
 
       // Electrostatic potential at bmag_max_coord.
@@ -205,6 +238,8 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
         app->field->phi_smooth, damp->phi_m_global, damp->rate);
       // Multiply by the user's scaling profile.
       gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
+
+      damp->advance_func = gk_species_damping_advance_loss_cone;
     }
     else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
       // Allocate filtered distribution function array
@@ -218,6 +253,10 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       // Initialize fbar from the projection of the initial distribution
       // (will be set from the initial f in the main app loop)
       gkyl_array_set(damp->fbar, 0.0);
+
+      damp->advance_func = gk_species_damping_advance_low_pass_filter;
+      damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_enabled;
+      damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_enabled;
     }
 
     // Set function pointers chosen at runtime.
@@ -230,58 +269,85 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   }
 }
 
+void
+gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+}
+
 
 void
-gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
   const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
-  if (damp->type) {
-    struct timespec wst = gkyl_wall_clock();
-    if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
-      gkyl_array_set(f_buffer, 1.0, fin);
-      gkyl_array_scale_by_cell(f_buffer, damp->rate);
-      gkyl_array_accumulate(rhs, -1.0, f_buffer);
-    }
-    else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
-      // Find the potential at the mirror throat.
-      gkyl_dg_basis_ops_eval_array_at_coord_comp(phi, damp->bmag_max_coord,
-        app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
-      gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
+  gkyl_array_set(f_buffer, 1.0, fin);
+  gkyl_array_scale_by_cell(f_buffer, damp->rate);
 
-      // Project the loss cone mask.
-      gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
-        phi, damp->phi_m_global, damp->rate);
+  // Add damping to f and the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
 
-      // Assemble the damping term -scale_prof * mask * f.
-      gkyl_array_set(f_buffer, 1.0, fin);
-      gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
-      gkyl_array_scale_by_cell(f_buffer, damp->rate);
-      gkyl_array_accumulate(rhs, -1.0, f_buffer);
+void
+gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  // Find the potential at the mirror throat.
+  gkyl_dg_basis_ops_eval_array_at_coord_comp(phi, damp->bmag_max_coord,
+    app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
+  gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
+
+  // Project the loss cone mask.
+  gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
+    phi, damp->phi_m_global, damp->rate);
+
+  // Assemble the damping term -scale_prof * mask * f.
+  gkyl_array_set(f_buffer, 1.0, fin);
+  gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
+  gkyl_array_scale_by_cell(f_buffer, damp->rate);
+
+  // Add damping to f and the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
 
-    }
-    else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-      // Match fbar state to the RK stage of fin.
-      const struct gkyl_array *fbar_in = damp->fbar;
-      if (fin == gks->f1)
-        fbar_in = damp->fbar1;
-      else if (fin == gks->fnew)
-        fbar_in = damp->fbarnew;
-
-      // Compute f - fbar and scale by the damping rate: rate * (f - fbar)
-      gkyl_array_set(f_buffer, 1.0, fin);      // f_buffer = f
-      gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = f - fbar
-      gkyl_array_scale_by_cell(f_buffer, damp->rate);    // f_buffer = rate * (f - fbar)
-      
-      // Add damping term to RHS: df/dt -= rate * (f - fbar)
-      gkyl_array_accumulate(rhs, -1.0, f_buffer);
-    }
+void
+gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  // Match fbar state to the RK stage of fin.
+  const struct gkyl_array *fbar_in = damp->fbar;
+  if (fin == gks->f1)
+    fbar_in = damp->fbar1;
+  else if (fin == gks->fnew)
+    fbar_in = damp->fbarnew;
+
+  // Compute f - fbar and scale by the damping rate: rate * (f - fbar)
+  gkyl_array_set(f_buffer, 1.0, fin);      // f_buffer = f
+  gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = f - fbar
+  gkyl_array_scale_by_cell(f_buffer, damp->rate);    // f_buffer = rate * (f - fbar)
+  
+  // Add damping term to RHS: df/dt -= rate * (f - fbar)
+  // Add to the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
 
-    // Add the frequency to the CFL frequency.
-    gkyl_array_accumulate(cflrate, 1.0, damp->rate);
 
-    app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
-  }
+void
+gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  struct timespec wst = gkyl_wall_clock();
+
+  damp->advance_func(app, gks, damp, phi, fin, f_buffer, rhs, cflrate);
+  
+  app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
 }
 
 
@@ -291,16 +357,27 @@ gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks, doub
   gks->damping.write_func(app, gks, tm, frame);
 }
 
+void
+gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+}
+
+void
+gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  // rhs_fbar = rate * (f - fbar)
+  gkyl_array_set(rhs_fbar, 1.0, fin);
+  gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
+  gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
+}
+
 void
 gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
 {
-  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-    // rhs_fbar = rate * (f - fbar)
-    gkyl_array_set(rhs_fbar, 1.0, fin);
-    gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
-    gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
-  }
+  damp->calc_fbar_rhs_func(damp, fin, fbar_in, rhs_fbar);
 }
 
 void
@@ -350,5 +427,5 @@ gk_species_damping_reset(gkyl_gyrokinetic_app* app, double tm, struct gk_species
 
   gks->info.damping = damp_inp;
   gk_species_damping_init(app, gks, damp);
-  gk_species_damping_init_fbar_from_f(gks, damp, gks->f);
+  gk_species_damping_set_fbar_to_f(gks, damp, gks->f);
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index c95238d3d3..486715a358 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -814,6 +814,13 @@ struct gk_damping {
   struct gkyl_array *fbar_host; // Host copy of fbar for use in IO.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
+  void (*advance_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+    const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+    struct gkyl_array *rhs, struct gkyl_array *cflrate);
+  void (*set_fbar_to_f_func)(const struct gk_species *gks, struct gk_damping *damp,
+    const struct gkyl_array *f);
+  void (*calc_fbar_rhs_func)(const struct gk_damping *damp,
+    const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
 };
 
 struct gk_fdot_multiplier {
@@ -2783,7 +2790,7 @@ void gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species
  * @param damp Species damping object.
  * @param f Initial distribution function.
  */
-void gk_species_damping_init_fbar_from_f(const struct gk_species *gks, struct gk_damping *damp, 
+void gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping *damp, 
   const struct gkyl_array *f);
 
 

From a768b4f15748533b4dd8fab9c68a0d335c6651c1 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Fri, 27 Mar 2026 19:16:16 -0400
Subject: [PATCH 04/25] Uncrustify format

---
 gyrokinetic/apps/gk_species_damping.c | 135 +++++++++++++++-----------
 1 file changed, 76 insertions(+), 59 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 1e8bf3cb58..8d45d053ff 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -13,7 +13,8 @@ void gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const stru
 void gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
   struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
-void gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+void gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app,
+  const struct gk_species *gks,
   struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
   struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
 void gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
@@ -21,9 +22,8 @@ void gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
 void gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
 
-
 void
-gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk_damping *damp,
   const struct gkyl_array *f)
 {
   gkyl_array_set(damp->fbar, 1.0, f);
@@ -32,25 +32,27 @@ gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk
 }
 
 void
-gk_species_damping_set_fbar_to_f_disabled(const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_set_fbar_to_f_disabled(const struct gk_species *gks, struct gk_damping *damp,
   const struct gkyl_array *f)
 {
 }
 
 void
-gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping *damp,
   const struct gkyl_array *f)
 {
   damp->set_fbar_to_f_func(gks, damp, f);
 }
 
 void
-gk_species_damping_write_disabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_write_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
+  int frame)
 {
 }
 
 void
-gk_species_damping_write_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_write_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
+  int frame)
 {
   struct timespec wst = gkyl_wall_clock();
   // DG metadata for damping rate.
@@ -58,19 +60,21 @@ gk_species_damping_write_enabled(gkyl_gyrokinetic_app* app, struct gk_species *g
     { .key = "poly_order", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = 0 },
     { .key = "basis_type", .elem_type = GKYL_MP_STRING, .cval = "serendipity" },
   };
-  int mpe_drate_len = sizeof(mpe_drate)/sizeof(mpe_drate[0]);
+  int mpe_drate_len = sizeof(mpe_drate) / sizeof(mpe_drate[0]);
   // Update app basic metada with time/frame.
   gkyl_msgpack_map_elem_set_double(app->io_meta_basic_len, app->io_meta_basic, "time", tm);
   gkyl_msgpack_map_elem_set_uint(app->io_meta_basic_len, app->io_meta_basic, "frame", frame);
   // Package metadata.
-  int io_meta_len[] = {app->io_meta_basic_len, mpe_drate_len, app->gk_geom->io_meta_len};
-  const struct gkyl_msgpack_map_elem* io_meta[] = {app->io_meta_basic, mpe_drate, app->gk_geom->io_meta};
-  struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union(sizeof(io_meta_len)/sizeof(int), io_meta_len, io_meta);
+  int io_meta_len[] = { app->io_meta_basic_len, mpe_drate_len, app->gk_geom->io_meta_len };
+  const struct gkyl_msgpack_map_elem *io_meta[] = { app->io_meta_basic, mpe_drate,
+                                                    app->gk_geom->io_meta };
+  struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union(sizeof(io_meta_len) / sizeof(int),
+    io_meta_len, io_meta);
 
   // Write out the damping rate.
   const char *fmt = "%s-%s_damping_rate_%d.gkyl";
   int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name, frame);
-  char fileNm[sz+1]; // ensures no buffer overflow
+  char fileNm[sz + 1]; // ensures no buffer overflow
   snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name, frame);
 
   // Copy data from device to host before writing it out.
@@ -80,12 +84,13 @@ gk_species_damping_write_enabled(gkyl_gyrokinetic_app* app, struct gk_species *g
   gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt, gks->damping.rate_host, fileNm);
   app->stat.n_io += 1;
 
-  gkyl_msgpack_data_release(mt); 
+  gkyl_msgpack_data_release(mt);
   app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wst);
 }
 
 void
-gk_species_damping_write_init_only(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_write_init_only(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
+  int frame)
 {
   gk_species_damping_write_enabled(app, gks, tm, frame);
   gks->damping.write_func = gk_species_damping_write_disabled;
@@ -117,33 +122,33 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
   if (damp->type) {
     // Allocate rate array.
-    damp->rate = mkarr(app->use_gpu, num_quad==1? 1 : gks->basis.num_basis, gks->local_ext.volume);
+    damp->rate = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
+      gks->local_ext.volume);
     damp->rate_host = damp->rate;
     if (app->use_gpu)
-      damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size); 
+      damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size);
 
     if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
       struct gk_proj_on_basis_c2p_func_ctx proj_on_basis_c2p_ctx; // c2p function context.
       proj_on_basis_c2p_ctx.cdim = app->cdim;
       proj_on_basis_c2p_ctx.vdim = gks->local_vel.ndim;
       proj_on_basis_c2p_ctx.vel_map = gks->vel_map;
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
-          .grid = &gks->grid,
-          .basis = &gks->basis,
-          .num_quad = num_quad,
-          .num_ret_vals = 1,
-          .eval = gks->info.damping.rate_profile,
-          .ctx = gks->info.damping.rate_profile_ctx,
-          .c2p_func = proj_on_basis_c2p_phase_func,
-          .c2p_func_ctx = &proj_on_basis_c2p_ctx,
-        }
-      );
+      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_inew(&(struct gkyl_proj_on_basis_inp) {
+        .grid = &gks->grid,
+        .basis = &gks->basis,
+        .num_quad = num_quad,
+        .num_ret_vals = 1,
+        .eval = gks->info.damping.rate_profile,
+        .ctx = gks->info.damping.rate_profile_ctx,
+        .c2p_func = proj_on_basis_c2p_phase_func,
+        .c2p_func_ctx = &proj_on_basis_c2p_ctx,
+      });
       gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, damp->rate_host);
       gkyl_proj_on_basis_release(projup);
       gkyl_array_copy(damp->rate, damp->rate_host);
 
       if (num_quad == 1)
-        gkyl_array_scale_range(damp->rate, 1.0/pow(sqrt(2.0),gks->grid.ndim), &gks->local);
+        gkyl_array_scale_range(damp->rate, 1.0 / pow(sqrt(2.0), gks->grid.ndim), &gks->local);
 
       damp->advance_func = gk_species_damping_advance_user_input;
     }
@@ -154,32 +159,38 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       // NOTE: if the same max bmag occurs at multiple locations,
       // bmag_max_coord may have different values on different MPI processes.
       double bmag_max_coord_ho[GKYL_MAX_CDIM];
-      double bmag_max_ho = gkyl_gk_geometry_reduce_arg_bmag(app->gk_geom, GKYL_MAX, bmag_max_coord_ho);
+      double bmag_max_ho = gkyl_gk_geometry_reduce_arg_bmag(app->gk_geom, GKYL_MAX,
+        bmag_max_coord_ho);
       double bmag_max_local = bmag_max_ho;
       double bmag_max_global;
-      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, &bmag_max_local, &bmag_max_global);
+      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, &bmag_max_local,
+        &bmag_max_global);
       double bmag_max_coord_local[app->cdim], bmag_max_coord_global[app->cdim];
       if (fabs(bmag_max_ho - bmag_max_global) < 1e-16) {
-        for (int d=0; d<app->cdim; d++)
+        for (int d = 0; d < app->cdim; d++) {
           bmag_max_coord_local[d] = bmag_max_coord_ho[d];
+        }
       }
       else {
-        for (int d=0; d<app->cdim; d++)
+        for (int d = 0; d < app->cdim; d++) {
           bmag_max_coord_local[d] = -DBL_MAX;
+        }
       }
-      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, app->cdim, bmag_max_coord_local, bmag_max_coord_global);
+      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, app->cdim, bmag_max_coord_local,
+        bmag_max_coord_global);
 
       if (app->use_gpu) {
         damp->bmag_max = gkyl_cu_malloc(sizeof(double));
-        damp->bmag_max_coord = gkyl_cu_malloc(app->cdim*sizeof(double));
+        damp->bmag_max_coord = gkyl_cu_malloc(app->cdim * sizeof(double));
         gkyl_cu_memcpy(damp->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
-        gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double),
+          GKYL_CU_MEMCPY_H2D);
       }
       else {
         damp->bmag_max = gkyl_malloc(sizeof(double));
-        damp->bmag_max_coord = gkyl_malloc(app->cdim*sizeof(double));
+        damp->bmag_max_coord = gkyl_malloc(app->cdim * sizeof(double));
         memcpy(damp->bmag_max, &bmag_max_global, sizeof(double));
-        memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double));
+        memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double));
       }
 
       // Electrostatic potential at bmag_max_coord.
@@ -197,9 +208,9 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
         .phase_grid = &gks->grid,
         .conf_basis = &app->basis,
         .phase_basis = &gks->basis,
-        .conf_range =  &app->local,
+        .conf_range = &app->local,
         .conf_range_ext = &app->local_ext,
-        .vel_range = &gks->local_vel, 
+        .vel_range = &gks->local_vel,
         .vel_map = gks->vel_map,
         .bmag = app->gk_geom->geo_int.bmag,
         .bmag_max = damp->bmag_max,
@@ -209,21 +220,25 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
         .num_quad = num_quad,
         .use_gpu = app->use_gpu,
       };
-      damp->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew( &inp_proj );
+      damp->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew(&inp_proj);
 
       // Project the conf-space rate profile provided.
-      struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-      struct gkyl_array *scale_prof_high_order_ho = app->use_gpu? mkarr(false, scale_prof_high_order->ncomp, scale_prof_high_order->size)
+      struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, gks->basis.num_basis,
+        gks->local_ext.volume);
+      struct gkyl_array *scale_prof_high_order_ho = app->use_gpu? mkarr(false,
+        scale_prof_high_order->ncomp, scale_prof_high_order->size)
                                                      : gkyl_array_acquire(scale_prof_high_order);
-      
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1, 
+
+      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1,
         gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
       gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
       gkyl_proj_on_basis_release(projup);
       gkyl_array_copy(scale_prof_high_order, scale_prof_high_order_ho);
 
-      damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis, gks->local_ext.volume);
-      gkyl_array_set_offset(damp->scale_prof, pow(sqrt(2.0),gks->grid.ndim), scale_prof_high_order, 0);
+      damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
+        gks->local_ext.volume);
+      gkyl_array_set_offset(damp->scale_prof, pow(sqrt(2.0), gks->grid.ndim), scale_prof_high_order,
+        0);
 
       gkyl_array_release(scale_prof_high_order_ho);
       gkyl_array_release(scale_prof_high_order);
@@ -270,15 +285,16 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 }
 
 void
-gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
   const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
 }
 
-
 void
-gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
   const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
@@ -291,7 +307,8 @@ gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk
 }
 
 void
-gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
   const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
@@ -315,7 +332,8 @@ gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_
 }
 
 void
-gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
   const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
@@ -330,29 +348,28 @@ gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const stru
   gkyl_array_set(f_buffer, 1.0, fin);      // f_buffer = f
   gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = f - fbar
   gkyl_array_scale_by_cell(f_buffer, damp->rate);    // f_buffer = rate * (f - fbar)
-  
+
   // Add damping term to RHS: df/dt -= rate * (f - fbar)
   // Add to the CFL frequency.
   gkyl_array_accumulate(rhs, -1.0, f_buffer);
   gkyl_array_accumulate(cflrate, 1.0, damp->rate);
 }
 
-
 void
-gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
   const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
   struct timespec wst = gkyl_wall_clock();
 
   damp->advance_func(app, gks, damp, phi, fin, f_buffer, rhs, cflrate);
-  
+
   app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
 }
 
-
 void
-gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_write(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm, int frame)
 {
   gks->damping.write_func(app, gks, tm, frame);
 }
@@ -418,10 +435,10 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
 }
 
 void
-gk_species_damping_reset(gkyl_gyrokinetic_app* app, double tm, struct gk_species *gks,
+gk_species_damping_reset(gkyl_gyrokinetic_app *app, double tm, struct gk_species *gks,
   struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp)
 {
-  (void) tm;
+  (void)tm;
 
   gk_species_damping_release(app, damp);
 

From 107f2cd32aa6e36a6847354b91955fd777036f11 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Fri, 27 Mar 2026 19:17:51 -0400
Subject: [PATCH 05/25] Refactor damping functions in gk_species_damping.c for
 improved clarity and functionality

---
 gyrokinetic/apps/gk_species_damping.c | 245 ++++++++++++--------------
 1 file changed, 114 insertions(+), 131 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 8d45d053ff..bd296585e8 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -4,24 +4,6 @@
 #include <gkyl_alloc.h>
 #include <gkyl_dg_basis_ops.h>
 
-void gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
-  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
-void gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
-  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
-void gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
-  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
-void gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app,
-  const struct gk_species *gks,
-  struct gk_damping *damp, const struct gkyl_array *phi, const struct gkyl_array *fin,
-  struct gkyl_array *f_buffer, struct gkyl_array *rhs, struct gkyl_array *cflrate);
-void gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
-  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
-void gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
-  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
-
 void
 gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk_damping *damp,
   const struct gkyl_array *f)
@@ -104,6 +86,120 @@ proj_on_basis_c2p_phase_func(const double *xcomp, double *xphys, void *ctx)
   gkyl_velocity_map_eval_c2p(c2p_ctx->vel_map, &xcomp[cdim], &xphys[cdim]);
 }
 
+
+void
+gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+}
+
+void
+gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  gkyl_array_set(f_buffer, 1.0, fin);
+  gkyl_array_scale_by_cell(f_buffer, damp->rate);
+
+  // Add damping to f and the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
+
+void
+gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  // Find the potential at the mirror throat.
+  gkyl_dg_basis_ops_eval_array_at_coord_comp(phi, damp->bmag_max_coord,
+    app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
+  gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
+
+  // Project the loss cone mask.
+  gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
+    phi, damp->phi_m_global, damp->rate);
+
+  // Assemble the damping term -scale_prof * mask * f.
+  gkyl_array_set(f_buffer, 1.0, fin);
+  gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
+  gkyl_array_scale_by_cell(f_buffer, damp->rate);
+
+  // Add damping to f and the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
+
+void
+gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  // Match fbar state to the RK stage of fin.
+  const struct gkyl_array *fbar_in = damp->fbar;
+  if (fin == gks->f1)
+    fbar_in = damp->fbar1;
+  else if (fin == gks->fnew)
+    fbar_in = damp->fbarnew;
+
+  // Compute f - fbar and scale by the damping rate: rate * (f - fbar)
+  gkyl_array_set(f_buffer, 1.0, fin);      // f_buffer = f
+  gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = f - fbar
+  gkyl_array_scale_by_cell(f_buffer, damp->rate);    // f_buffer = rate * (f - fbar)
+
+  // Add damping term to RHS: df/dt -= rate * (f - fbar)
+  // Add to the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
+
+void
+gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  struct timespec wst = gkyl_wall_clock();
+
+  damp->advance_func(app, gks, damp, phi, fin, f_buffer, rhs, cflrate);
+
+  app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
+}
+
+void
+gk_species_damping_write(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm, int frame)
+{
+  gks->damping.write_func(app, gks, tm, frame);
+}
+
+void
+gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+}
+
+void
+gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  // rhs_fbar = rate * (f - fbar)
+  gkyl_array_set(rhs_fbar, 1.0, fin);
+  gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
+  gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
+}
+
+void
+gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  damp->calc_fbar_rhs_func(damp, fin, fbar_in, rhs_fbar);
+}
+
 void
 gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
   struct gk_damping *damp)
@@ -284,119 +380,6 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   }
 }
 
-void
-gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
-  struct gkyl_array *rhs, struct gkyl_array *cflrate)
-{
-}
-
-void
-gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
-  struct gkyl_array *rhs, struct gkyl_array *cflrate)
-{
-  gkyl_array_set(f_buffer, 1.0, fin);
-  gkyl_array_scale_by_cell(f_buffer, damp->rate);
-
-  // Add damping to f and the CFL frequency.
-  gkyl_array_accumulate(rhs, -1.0, f_buffer);
-  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
-}
-
-void
-gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
-  struct gkyl_array *rhs, struct gkyl_array *cflrate)
-{
-  // Find the potential at the mirror throat.
-  gkyl_dg_basis_ops_eval_array_at_coord_comp(phi, damp->bmag_max_coord,
-    app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
-  gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
-
-  // Project the loss cone mask.
-  gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
-    phi, damp->phi_m_global, damp->rate);
-
-  // Assemble the damping term -scale_prof * mask * f.
-  gkyl_array_set(f_buffer, 1.0, fin);
-  gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
-  gkyl_array_scale_by_cell(f_buffer, damp->rate);
-
-  // Add damping to f and the CFL frequency.
-  gkyl_array_accumulate(rhs, -1.0, f_buffer);
-  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
-}
-
-void
-gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
-  struct gkyl_array *rhs, struct gkyl_array *cflrate)
-{
-  // Match fbar state to the RK stage of fin.
-  const struct gkyl_array *fbar_in = damp->fbar;
-  if (fin == gks->f1)
-    fbar_in = damp->fbar1;
-  else if (fin == gks->fnew)
-    fbar_in = damp->fbarnew;
-
-  // Compute f - fbar and scale by the damping rate: rate * (f - fbar)
-  gkyl_array_set(f_buffer, 1.0, fin);      // f_buffer = f
-  gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = f - fbar
-  gkyl_array_scale_by_cell(f_buffer, damp->rate);    // f_buffer = rate * (f - fbar)
-
-  // Add damping term to RHS: df/dt -= rate * (f - fbar)
-  // Add to the CFL frequency.
-  gkyl_array_accumulate(rhs, -1.0, f_buffer);
-  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
-}
-
-void
-gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
-  struct gkyl_array *rhs, struct gkyl_array *cflrate)
-{
-  struct timespec wst = gkyl_wall_clock();
-
-  damp->advance_func(app, gks, damp, phi, fin, f_buffer, rhs, cflrate);
-
-  app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
-}
-
-void
-gk_species_damping_write(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm, int frame)
-{
-  gks->damping.write_func(app, gks, tm, frame);
-}
-
-void
-gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
-  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
-{
-}
-
-void
-gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
-  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
-{
-  // rhs_fbar = rate * (f - fbar)
-  gkyl_array_set(rhs_fbar, 1.0, fin);
-  gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
-  gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
-}
-
-void
-gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
-  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
-{
-  damp->calc_fbar_rhs_func(damp, fin, fbar_in, rhs_fbar);
-}
-
 void
 gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp)
 {

From 3f3be48effa66ccf2faea0a41bd185b17a9dabda Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Fri, 27 Mar 2026 19:23:01 -0400
Subject: [PATCH 06/25] Refactor gk_species_damping_step_fbar_ssp_rk3 for
 improved clarity and functionality

---
 gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
index d6befeb874..6d3dce5447 100644
--- a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
+++ b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
@@ -4,6 +4,17 @@ static void
 gk_species_damping_step_fbar_ssp_rk3(struct gk_species *gks, const struct gkyl_array *fin,
   const struct gkyl_array *fbar_in, struct gkyl_array *fbar_out, double dt)
 {
+  // Okay so this is not good
+  // There is significant amounts of refactoring needed to be done here
+  // This if statement is no good
+  // We should have these functions point to pointers
+  // We have to be careful that we do not step uninitilized data
+  // It would be so much better if we can combine this into stepping like a species
+  // We should't have to initilize a seperate species, but that's almost what this is doing. Damping is creating a new distribution function which must be stepped.
+  // The stepping should happen inside the SSPRK3 algorithm
+  // I'm not sure how sundials will change all this logic, so we should wait for that.
+  // The current implementation here is functional, but it's not designed well
+  
   if (gks->damping.type != GKYL_GK_DAMPING_LOW_PASS_FILTER)
     return;
 

From f35880e8c31fd3369b2d29bc3d1387ae599fdd0f Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 28 Mar 2026 12:14:04 -0400
Subject: [PATCH 07/25] Implement low-pass filter damping functions and
 integrate into gyrokinetic update process

---
 gyrokinetic/apps/gk_species_damping.c         | 162 +++++++++++++++---
 gyrokinetic/apps/gkyl_gyrokinetic.h           |   6 +
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h      |  42 +++++
 gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c |  75 ++------
 4 files changed, 203 insertions(+), 82 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index bd296585e8..9b192d46d0 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -86,6 +86,35 @@ proj_on_basis_c2p_phase_func(const double *xcomp, double *xphys, void *ctx)
   gkyl_velocity_map_eval_c2p(c2p_ctx->vel_map, &xcomp[cdim], &xphys[cdim]);
 }
 
+static void
+gk_species_damping_project_phase_rate(const struct gkyl_gyrokinetic_app *app,
+  const struct gk_species *gks, int num_quad,
+  void (*rate_profile)(double t, const double *xn, double *fout, void *ctx),
+  void *rate_profile_ctx, struct gkyl_array *rate_host, struct gkyl_array *rate)
+{
+  struct gk_proj_on_basis_c2p_func_ctx proj_on_basis_c2p_ctx; // c2p function context.
+  proj_on_basis_c2p_ctx.cdim = app->cdim;
+  proj_on_basis_c2p_ctx.vdim = gks->local_vel.ndim;
+  proj_on_basis_c2p_ctx.vel_map = gks->vel_map;
+
+  gkyl_proj_on_basis *projup = gkyl_proj_on_basis_inew(&(struct gkyl_proj_on_basis_inp) {
+    .grid = &gks->grid,
+    .basis = &gks->basis,
+    .num_quad = num_quad,
+    .num_ret_vals = 1,
+    .eval = rate_profile,
+    .ctx = rate_profile_ctx,
+    .c2p_func = proj_on_basis_c2p_phase_func,
+    .c2p_func_ctx = &proj_on_basis_c2p_ctx,
+  });
+  gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, rate_host);
+  gkyl_proj_on_basis_release(projup);
+  gkyl_array_copy(rate, rate_host);
+
+  if (num_quad == 1)
+    gkyl_array_scale_range(rate, 1.0 / pow(sqrt(2.0), gks->grid.ndim), &gks->local);
+}
+
 
 void
 gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
@@ -200,12 +229,81 @@ gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
   damp->calc_fbar_rhs_func(damp, fin, fbar_in, rhs_fbar);
 }
 
+static void
+gk_species_damping_forward_euler_disabled(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt)
+{
+}
+
+static void
+gk_species_damping_forward_euler_enabled(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt)
+{
+  // Reuse species scratch storage to assemble fbar RHS for this FE substep.
+  gk_species_damping_calc_fbar_rhs(&gks->damping, fin, fbar_in, gks->lte.f_lte);
+  gk_species_step_f(gks, fbar_out, dt, fbar_in);
+}
+
+void
+gk_species_damping_forward_euler(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt)
+{
+  gks->damping.forward_euler_func(gks, fin, fbar_in, fbar_out, dt);
+}
+
+static void
+gk_species_damping_combine_disabled(struct gk_species *gks, struct gkyl_array *fout,
+  double c1, const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng)
+{
+}
+
+static void
+gk_species_damping_combine_enabled(struct gk_species *gks, struct gkyl_array *fout,
+  double c1, const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng)
+{
+  gk_species_combine(gks, fout, c1, f1, c2, f2, rng);
+}
+
+void
+gk_species_damping_combine(struct gk_species *gks, struct gkyl_array *fout, double c1,
+  const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng)
+{
+  gks->damping.combine_func(gks, fout, c1, f1, c2, f2, rng);
+}
+
+static void
+gk_species_damping_copy_range_disabled(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range)
+{
+}
+
+static void
+gk_species_damping_copy_range_enabled(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range)
+{
+  gk_species_copy_range(gks, fout, fin, range);
+}
+
+void
+gk_species_damping_copy_range(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range)
+{
+  gks->damping.copy_func(gks, fout, fin, range);
+}
+
 void
 gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
   struct gk_damping *damp)
 {
   damp->type = gks->info.damping.type;
   damp->evolve = false; // Whether the rate is time dependent.
+  const double damping_const = gks->info.damping.damping_const;
 
   int num_quad = gks->info.damping.num_quad? gks->info.damping.num_quad : 1; // Default is a p=0 mask.
   assert(num_quad == 1); // MF 2025/06/11: Limited to this for now.
@@ -215,6 +313,9 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->advance_func = gk_species_damping_advance_disabled;
   damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_disabled;
   damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_disabled;
+  damp->forward_euler_func = gk_species_damping_forward_euler_disabled;
+  damp->combine_func = gk_species_damping_combine_disabled;
+  damp->copy_func = gk_species_damping_copy_range_disabled;
 
   if (damp->type) {
     // Allocate rate array.
@@ -225,26 +326,15 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size);
 
     if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
-      struct gk_proj_on_basis_c2p_func_ctx proj_on_basis_c2p_ctx; // c2p function context.
-      proj_on_basis_c2p_ctx.cdim = app->cdim;
-      proj_on_basis_c2p_ctx.vdim = gks->local_vel.ndim;
-      proj_on_basis_c2p_ctx.vel_map = gks->vel_map;
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_inew(&(struct gkyl_proj_on_basis_inp) {
-        .grid = &gks->grid,
-        .basis = &gks->basis,
-        .num_quad = num_quad,
-        .num_ret_vals = 1,
-        .eval = gks->info.damping.rate_profile,
-        .ctx = gks->info.damping.rate_profile_ctx,
-        .c2p_func = proj_on_basis_c2p_phase_func,
-        .c2p_func_ctx = &proj_on_basis_c2p_ctx,
-      });
-      gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, damp->rate_host);
-      gkyl_proj_on_basis_release(projup);
-      gkyl_array_copy(damp->rate, damp->rate_host);
-
-      if (num_quad == 1)
-        gkyl_array_scale_range(damp->rate, 1.0 / pow(sqrt(2.0), gks->grid.ndim), &gks->local);
+      // USER_INPUT supports either a projected profile, or a uniform constant if no profile is provided.
+      if (gks->info.damping.rate_profile) {
+        gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
+          gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
+        gkyl_array_scale(damp->rate, damping_const == 0.0 ? 1.0 : damping_const);
+      }
+      else {
+        gkyl_array_clear(damp->rate, damping_const);
+      }
 
       damp->advance_func = gk_species_damping_advance_user_input;
     }
@@ -325,11 +415,18 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
         scale_prof_high_order->ncomp, scale_prof_high_order->size)
                                                      : gkyl_array_acquire(scale_prof_high_order);
 
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1,
-        gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
-      gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
-      gkyl_proj_on_basis_release(projup);
-      gkyl_array_copy(scale_prof_high_order, scale_prof_high_order_ho);
+      if (gks->info.damping.rate_profile) {
+        gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1,
+          gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
+        gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
+        gkyl_proj_on_basis_release(projup);
+        gkyl_array_copy(scale_prof_high_order, scale_prof_high_order_ho);
+      }
+      else {
+        gkyl_array_clear(scale_prof_high_order, 1.0);
+      }
+
+      gkyl_array_scale(scale_prof_high_order, damping_const == 0.0 ? 1.0 : damping_const);
 
       damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
         gks->local_ext.volume);
@@ -353,6 +450,16 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       damp->advance_func = gk_species_damping_advance_loss_cone;
     }
     else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+      // LOW_PASS_FILTER supports either a projected phase-space profile or a uniform constant.
+      if (gks->info.damping.rate_profile) {
+        gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
+          gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
+        gkyl_array_scale(damp->rate, damping_const == 0.0 ? 1.0 : damping_const);
+      }
+      else {
+        gkyl_array_clear(damp->rate, damping_const);
+      }
+
       // Allocate filtered distribution function array
       damp->fbar = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
       damp->fbar1 = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
@@ -363,11 +470,14 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
       // Initialize fbar from the projection of the initial distribution
       // (will be set from the initial f in the main app loop)
-      gkyl_array_set(damp->fbar, 0.0);
+      gkyl_array_clear(damp->fbar, 0.0);
 
       damp->advance_func = gk_species_damping_advance_low_pass_filter;
       damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_enabled;
       damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_enabled;
+      damp->forward_euler_func = gk_species_damping_forward_euler_enabled;
+      damp->combine_func = gk_species_damping_combine_enabled;
+      damp->copy_func = gk_species_damping_copy_range_enabled;
     }
 
     // Set function pointers chosen at runtime.
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index bc846ea381..c9d447c9e6 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -329,6 +329,12 @@ struct gkyl_gyrokinetic_damping {
   //   - I_loss(z) * scale_factor * scale_profile(z), where I_loss(z) is =1 in the loss
   //     cone and 0 in the confined region (type = GKYL_PROPORTIONAL_TERM_LOSS_CONE).
   enum gkyl_gyrokinetic_damping_type type;
+  // Optional constant damping factor.
+  // Semantics:
+  //  - if rate_profile is NULL: use damping_const as the full rate profile.
+  //  - if rate_profile is provided: effective rate is damping_const*rate_profile,
+  //    where damping_const defaults to 1.0 when left as 0.0.
+  double damping_const;
   void (*rate_profile)(double t, const double *xn, double *fout, void *ctx);
   void *rate_profile_ctx; // Context for rate_profile function.
   int num_quad; // Number of quadrature points in each direction to use in projecting the rate.
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 486715a358..d2bbc9fa5d 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -821,6 +821,13 @@ struct gk_damping {
     const struct gkyl_array *f);
   void (*calc_fbar_rhs_func)(const struct gk_damping *damp,
     const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
+  void (*forward_euler_func)(struct gk_species *gks, const struct gkyl_array *fin,
+    const struct gkyl_array *fbar_in, struct gkyl_array *fbar_out, double dt);
+  void (*combine_func)(struct gk_species *gks, struct gkyl_array *fout, double c1,
+    const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+    const struct gkyl_range *rng);
+  void (*copy_func)(struct gk_species *gks, struct gkyl_array *fout,
+    const struct gkyl_array *fin, const struct gkyl_range *range);
 };
 
 struct gk_fdot_multiplier {
@@ -2832,6 +2839,41 @@ void gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks,
 void gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
 
+/**
+ * Take a forward-Euler substep for the filtered distribution.
+ *
+ * This is a runtime-dispatched no-op unless the damping type enables
+ * filtered-state evolution (currently low-pass filter).
+ *
+ * @param gks Species object.
+ * @param fin Current distribution function f.
+ * @param fbar_in Current filtered distribution state.
+ * @param fbar_out Output filtered distribution state after FE substep.
+ * @param dt Time-step for substep.
+ */
+void gk_species_damping_forward_euler(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt);
+
+/**
+ * Combine filtered distributions like gk_species_combine.
+ *
+ * This is runtime-dispatched and a no-op unless the damping type enables
+ * filtered-state evolution (currently low-pass filter).
+ */
+void gk_species_damping_combine(struct gk_species *gks, struct gkyl_array *fout, double c1,
+  const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng);
+
+/**
+ * Copy filtered distribution ranges like gk_species_copy_range.
+ *
+ * This is runtime-dispatched and a no-op unless the damping type enables
+ * filtered-state evolution (currently low-pass filter).
+ */
+void gk_species_damping_copy_range(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range);
+
 /**
  * Release species damping object.
  *
diff --git a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
index 6d3dce5447..22ec1cf2bc 100644
--- a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
+++ b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
@@ -1,31 +1,9 @@
 #include <gkyl_gyrokinetic_priv.h>
 
-static void
-gk_species_damping_step_fbar_ssp_rk3(struct gk_species *gks, const struct gkyl_array *fin,
-  const struct gkyl_array *fbar_in, struct gkyl_array *fbar_out, double dt)
-{
-  // Okay so this is not good
-  // There is significant amounts of refactoring needed to be done here
-  // This if statement is no good
-  // We should have these functions point to pointers
-  // We have to be careful that we do not step uninitilized data
-  // It would be so much better if we can combine this into stepping like a species
-  // We should't have to initilize a seperate species, but that's almost what this is doing. Damping is creating a new distribution function which must be stepped.
-  // The stepping should happen inside the SSPRK3 algorithm
-  // I'm not sure how sundials will change all this logic, so we should wait for that.
-  // The current implementation here is functional, but it's not designed well
-  
-  if (gks->damping.type != GKYL_GK_DAMPING_LOW_PASS_FILTER)
-    return;
-
-  // Reuse existing species scratch buffer for rhs_fbar.
-  gk_species_damping_calc_fbar_rhs(&gks->damping, fin, fbar_in, gks->lte.f_lte);
-  gk_species_step_f(gks, fbar_out, dt, fbar_in);
-}
-
 static void
 gyrokinetic_forward_euler(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   const struct gkyl_array *fin[], struct gkyl_array *fout[], 
+  const struct gkyl_array *fbar_in[], struct gkyl_array *fbar_out[],
   struct gkyl_array **bflux_in[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], 
   struct gkyl_array **bflux_in_neut[], struct gkyl_array **bflux_out_neut[], 
@@ -49,6 +27,7 @@ gyrokinetic_forward_euler(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   for (int i=0; i<app->num_species; ++i) {
     struct gk_species *gks = &app->species[i];
     gk_species_step_f(gks, fout[i], dta, fin[i]);
+    gk_species_damping_forward_euler(gks, fin[i], fbar_in[i], fbar_out[i], dta);
     gk_species_bflux_accumulate(app, &gks->bflux, bflux_out[i], 1.0, bflux_in[i]);
   }
   for (int i=0; i<app->num_neut_species; ++i) {
@@ -69,6 +48,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
   // from the actual time-step.
   const struct gkyl_array *fin[app->num_species];
   struct gkyl_array *fout[app->num_species];
+  const struct gkyl_array *fbar_in[app->num_species];
+  struct gkyl_array *fbar_out[app->num_species];
   struct gkyl_array **bflux_in[app->num_species];
   struct gkyl_array **bflux_out[app->num_species];
 
@@ -90,6 +71,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           struct gk_species *gks = &app->species[i];
           fin[i] = gks->f;
           fout[i] = gks->f1;
+          fbar_in[i] = gks->damping.fbar;
+          fbar_out[i] = gks->damping.fbar1;
           // Boundary fluxes.
           bflux_in[i] = gks->bflux.f;
           bflux_out[i] = gks->bflux.f1;
@@ -109,17 +92,10 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           gk_species_source_adapt(app, gks, &gks->src, gks->lte.f_lte, bflux_in, tcurr);
         }
 
-        gyrokinetic_forward_euler(app, tcurr, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_forward_euler(app, tcurr, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
         dt = st.dt_actual;
 
-        // Advance filtered distributions with the same stage-1 FE update.
-        for (int i=0; i<app->num_species; ++i) {
-          struct gk_species *gks = &app->species[i];
-          gk_species_damping_step_fbar_ssp_rk3(gks, gks->f,
-            gks->damping.fbar, gks->damping.fbar1, dt);
-        }
-
         // Subtract boundary flux f from f1 so that we only step boundary
         // fluxes during a given time step, not over all time. And so that the
         // boundary flux in f is kept in case a later RK stage fails.
@@ -155,6 +131,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           struct gk_species *gks = &app->species[i];
           fin[i] = gks->f1;
           fout[i] = gks->fnew;
+          fbar_in[i] = gks->damping.fbar1;
+          fbar_out[i] = gks->damping.fbarnew;
           // Boundary fluxes.
           bflux_in[i] = gks->bflux.f1;
           bflux_out[i] = gks->bflux.fnew;
@@ -168,16 +146,9 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           bflux_out_neut[i] = gkns->bflux.fnew;
         }
 
-        gyrokinetic_forward_euler(app, tcurr+dt, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_forward_euler(app, tcurr+dt, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
-        // Advance filtered distributions with the same stage-2 FE update.
-        for (int i=0; i<app->num_species; ++i) {
-          struct gk_species *gks = &app->species[i];
-          gk_species_damping_step_fbar_ssp_rk3(gks, gks->f1,
-            gks->damping.fbar1, gks->damping.fbarnew, st.dt_actual);
-        }
-
         if (st.dt_actual < dt) {
 
           // Recalculate the field.
@@ -205,10 +176,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           for (int i=0; i<app->num_species; ++i) {
             struct gk_species *gks = &app->species[i];
             gk_species_combine(gks, gks->f1, 3.0/4.0, gks->f, 1.0/4.0, gks->fnew, &gks->local_ext);
-            if (gks->damping.type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-              gk_species_combine(gks, gks->damping.fbar1, 3.0/4.0, gks->damping.fbar,
-                1.0/4.0, gks->damping.fbarnew, &gks->local_ext);
-            }
+            gk_species_damping_combine(gks, gks->damping.fbar1, 3.0/4.0, gks->damping.fbar,
+              1.0/4.0, gks->damping.fbarnew, &gks->local_ext);
             gk_species_bflux_set(app, &gks->bflux, gks->bflux.f1, 1.0/4.0, gks->bflux.fnew);
           }
           for (int i=0; i<app->num_neut_species; ++i) {
@@ -238,6 +207,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           struct gk_species *gks = &app->species[i];
           fin[i] = gks->f1;
           fout[i] = gks->fnew;
+          fbar_in[i] = gks->damping.fbar1;
+          fbar_out[i] = gks->damping.fbarnew;
           // Boundary fluxes.
           bflux_in[i] = gks->bflux.f1;
           bflux_out[i] = gks->bflux.fnew;
@@ -251,16 +222,9 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           bflux_out_neut[i] = gkns->bflux.fnew;
         }
 
-        gyrokinetic_forward_euler(app, tcurr+dt/2, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_forward_euler(app, tcurr+dt/2, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
-        // Advance filtered distributions with the same stage-3 FE update.
-        for (int i=0; i<app->num_species; ++i) {
-          struct gk_species *gks = &app->species[i];
-          gk_species_damping_step_fbar_ssp_rk3(gks, gks->f1,
-            gks->damping.fbar1, gks->damping.fbarnew, st.dt_actual);
-        }
-
         if (st.dt_actual < dt) {
           // Recalculate the field.
           for (int i=0; i<app->num_species; ++i) {
@@ -290,11 +254,10 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
             // Step f.
             gk_species_combine(gks, gks->f1, 1.0/3.0, gks->f, 2.0/3.0, gks->fnew, &gks->local_ext);
             gk_species_copy_range(gks, gks->f, gks->f1, &gks->local_ext);
-            if (gks->damping.type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-              gk_species_combine(gks, gks->damping.fbar1, 1.0/3.0, gks->damping.fbar,
-                2.0/3.0, gks->damping.fbarnew, &gks->local_ext);
-              gk_species_copy_range(gks, gks->damping.fbar, gks->damping.fbar1, &gks->local_ext);
-            }
+            gk_species_damping_combine(gks, gks->damping.fbar, 1.0/3.0, gks->damping.fbar1,
+              2.0/3.0, gks->damping.fbarnew, &gks->local_ext);
+            gk_species_damping_copy_range(gks, gks->damping.fbar1, gks->damping.fbar,
+              &gks->local_ext);
             // Step boundary fluxes.
             gk_species_bflux_set(app, &gks->bflux, gks->bflux.f, 2.0/3.0, gks->bflux.fnew);
             gk_species_bflux_calc_voltime_integrated_mom(app, gks, &gks->bflux, tcurr);

From e3bce23234d9053da3cda44ffa333577407f70ad Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 28 Mar 2026 12:17:16 -0400
Subject: [PATCH 08/25] Uncrustify format

---
 gyrokinetic/apps/gk_species_damping.c | 1 -
 1 file changed, 1 deletion(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 9b192d46d0..43a48d3368 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -115,7 +115,6 @@ gk_species_damping_project_phase_rate(const struct gkyl_gyrokinetic_app *app,
     gkyl_array_scale_range(rate, 1.0 / pow(sqrt(2.0), gks->grid.ndim), &gks->local);
 }
 
-
 void
 gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,

From 135599f97155e8d11c980204be8954a086c853a7 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 28 Mar 2026 12:21:54 -0400
Subject: [PATCH 09/25] Add comments for damping helpers and low-pass filter
 functions

---
 gyrokinetic/apps/gk_species_damping.c | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 43a48d3368..6d4e906281 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -4,6 +4,8 @@
 #include <gkyl_alloc.h>
 #include <gkyl_dg_basis_ops.h>
 
+// Damping state synchronization helpers.
+
 void
 gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk_damping *damp,
   const struct gkyl_array *f)
@@ -26,6 +28,8 @@ gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping
   damp->set_fbar_to_f_func(gks, damp, f);
 }
 
+// Damping diagnostics write helpers.
+
 void
 gk_species_damping_write_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
   int frame)
@@ -78,6 +82,8 @@ gk_species_damping_write_init_only(gkyl_gyrokinetic_app *app, struct gk_species
   gks->damping.write_func = gk_species_damping_write_disabled;
 }
 
+// Damping rate projection helpers.
+
 static void
 proj_on_basis_c2p_phase_func(const double *xcomp, double *xphys, void *ctx)
 {
@@ -115,6 +121,8 @@ gk_species_damping_project_phase_rate(const struct gkyl_gyrokinetic_app *app,
     gkyl_array_scale_range(rate, 1.0 / pow(sqrt(2.0), gks->grid.ndim), &gks->local);
 }
 
+// Damping RHS assembly dispatch helpers.
+
 void
 gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
@@ -205,6 +213,8 @@ gk_species_damping_write(gkyl_gyrokinetic_app *app, struct gk_species *gks, doub
   gks->damping.write_func(app, gks, tm, frame);
 }
 
+// Low-pass filter fbar RHS helpers.
+
 void
 gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
@@ -228,6 +238,8 @@ gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
   damp->calc_fbar_rhs_func(damp, fin, fbar_in, rhs_fbar);
 }
 
+// Low-pass filter stage update helpers.
+
 static void
 gk_species_damping_forward_euler_disabled(struct gk_species *gks,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
@@ -253,6 +265,8 @@ gk_species_damping_forward_euler(struct gk_species *gks,
   gks->damping.forward_euler_func(gks, fin, fbar_in, fbar_out, dt);
 }
 
+// Low-pass filter RK stage combine/copy helpers.
+
 static void
 gk_species_damping_combine_disabled(struct gk_species *gks, struct gkyl_array *fout,
   double c1, const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
@@ -296,6 +310,8 @@ gk_species_damping_copy_range(struct gk_species *gks, struct gkyl_array *fout,
   gks->damping.copy_func(gks, fout, fin, range);
 }
 
+// Damping object lifecycle.
+
 void
 gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
   struct gk_damping *damp)
@@ -526,6 +542,8 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
   }
 }
 
+// Damping runtime reset.
+
 void
 gk_species_damping_reset(gkyl_gyrokinetic_app *app, double tm, struct gk_species *gks,
   struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp)

From 98492c3c999cff489b26049315710e540faf88fb Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 28 Mar 2026 12:24:51 -0400
Subject: [PATCH 10/25] Update damping parameters to use rate_const instead of
 damping_const for clarity and consistency

---
 gyrokinetic/apps/gk_species_damping.c | 12 ++++++------
 gyrokinetic/apps/gkyl_gyrokinetic.h   |  8 ++++----
 gyrokinetic/creg/rt_gk_wham_1x2v_p1.c |  5 +++++
 3 files changed, 15 insertions(+), 10 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 6d4e906281..41eab39f77 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -318,7 +318,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 {
   damp->type = gks->info.damping.type;
   damp->evolve = false; // Whether the rate is time dependent.
-  const double damping_const = gks->info.damping.damping_const;
+  const double rate_const = gks->info.damping.rate_const;
 
   int num_quad = gks->info.damping.num_quad? gks->info.damping.num_quad : 1; // Default is a p=0 mask.
   assert(num_quad == 1); // MF 2025/06/11: Limited to this for now.
@@ -345,10 +345,10 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       if (gks->info.damping.rate_profile) {
         gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
           gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
-        gkyl_array_scale(damp->rate, damping_const == 0.0 ? 1.0 : damping_const);
+        gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
       }
       else {
-        gkyl_array_clear(damp->rate, damping_const);
+        gkyl_array_clear(damp->rate, rate_const);
       }
 
       damp->advance_func = gk_species_damping_advance_user_input;
@@ -441,7 +441,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
         gkyl_array_clear(scale_prof_high_order, 1.0);
       }
 
-      gkyl_array_scale(scale_prof_high_order, damping_const == 0.0 ? 1.0 : damping_const);
+      gkyl_array_scale(scale_prof_high_order, rate_const == 0.0 ? 1.0 : rate_const);
 
       damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
         gks->local_ext.volume);
@@ -469,10 +469,10 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       if (gks->info.damping.rate_profile) {
         gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
           gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
-        gkyl_array_scale(damp->rate, damping_const == 0.0 ? 1.0 : damping_const);
+        gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
       }
       else {
-        gkyl_array_clear(damp->rate, damping_const);
+        gkyl_array_clear(damp->rate, rate_const);
       }
 
       // Allocate filtered distribution function array
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index c9d447c9e6..6e02a15cc3 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -331,10 +331,10 @@ struct gkyl_gyrokinetic_damping {
   enum gkyl_gyrokinetic_damping_type type;
   // Optional constant damping factor.
   // Semantics:
-  //  - if rate_profile is NULL: use damping_const as the full rate profile.
-  //  - if rate_profile is provided: effective rate is damping_const*rate_profile,
-  //    where damping_const defaults to 1.0 when left as 0.0.
-  double damping_const;
+  //  - if rate_profile is NULL: use rate_const as the full rate profile.
+  //  - if rate_profile is provided: effective rate is rate_const*rate_profile,
+  //    where rate_const defaults to 1.0 when left as 0.0.
+  double rate_const;
   void (*rate_profile)(double t, const double *xn, double *fout, void *ctx);
   void *rate_profile_ctx; // Context for rate_profile function.
   int num_quad; // Number of quadrature points in each direction to use in projecting the rate.
diff --git a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
index b7462eb3c4..ab9391aca2 100644
--- a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
@@ -608,6 +608,11 @@ int main(int argc, char **argv)
       .write_diagnostics = true,
     },
 
+    .damping = {
+      .type = GKYL_GK_DAMPING_LOW_PASS_FILTER,
+      .rate_const = 1e-6,
+    },
+
     .source = {
       .source_id = GKYL_PROJ_SOURCE,
       .num_sources = 1,

From ec79063290d76df07fbaf5b37cb0fa9eb6754b76 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 28 Mar 2026 12:38:13 -0400
Subject: [PATCH 11/25] Add low-pass filter fbar write functions and update
 damping struct for diagnostics

---
 gyrokinetic/apps/gk_species_damping.c    | 51 ++++++++++++++++++++++--
 gyrokinetic/apps/gkyl_gyrokinetic.h      |  1 +
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h |  2 +
 3 files changed, 50 insertions(+), 4 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 41eab39f77..f20b115443 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -36,6 +36,41 @@ gk_species_damping_write_disabled(gkyl_gyrokinetic_app *app, struct gk_species *
 {
 }
 
+static void
+gk_species_damping_write_fbar_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+}
+
+static void
+gk_species_damping_write_fbar_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+  // Metadata from app, species, and geometry (same pattern as species f write).
+  gkyl_msgpack_map_elem_set_double(app->io_meta_basic_len, app->io_meta_basic, "time", tm);
+  gkyl_msgpack_map_elem_set_uint(app->io_meta_basic_len, app->io_meta_basic, "frame", frame);
+  int io_meta_fbar_len[] = {app->io_meta_basic_len, gks->io_meta_len, app->gk_geom->io_meta_len};
+  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = {app->io_meta_basic, gks->io_meta,
+    app->gk_geom->io_meta};
+  struct gkyl_msgpack_data *mt_fbar = gkyl_msgpack_create_union(
+    sizeof(io_meta_fbar_len) / sizeof(int), io_meta_fbar_len, io_meta_fbar);
+
+  const char *fmt_fbar = "%s-%s_fbar_%d.gkyl";
+  int sz_fbar = gkyl_calc_strlen(fmt_fbar, app->name, gks->info.name, frame);
+  char fileNm_fbar[sz_fbar + 1]; // ensures no buffer overflow
+  snprintf(fileNm_fbar, sizeof fileNm_fbar, fmt_fbar, app->name, gks->info.name, frame);
+
+  // Copy fbar from device to host before writing it out.
+  if (app->use_gpu)
+    gkyl_array_copy(gks->damping.fbar_host, gks->damping.fbar);
+
+  gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt_fbar, gks->damping.fbar_host,
+    fileNm_fbar);
+  app->stat.n_io += 1;
+
+  gkyl_msgpack_data_release(mt_fbar);
+}
+
 void
 gk_species_damping_write_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
   int frame)
@@ -70,6 +105,8 @@ gk_species_damping_write_enabled(gkyl_gyrokinetic_app *app, struct gk_species *g
   gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt, gks->damping.rate_host, fileNm);
   app->stat.n_io += 1;
 
+  gks->damping.write_fbar_func(app, gks, tm, frame);
+
   gkyl_msgpack_data_release(mt);
   app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wst);
 }
@@ -318,6 +355,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 {
   damp->type = gks->info.damping.type;
   damp->evolve = false; // Whether the rate is time dependent.
+  damp->write_fbar = gks->info.damping.write_fbar;
   const double rate_const = gks->info.damping.rate_const;
 
   int num_quad = gks->info.damping.num_quad? gks->info.damping.num_quad : 1; // Default is a p=0 mask.
@@ -325,6 +363,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
   // Default function pointers.
   damp->write_func = gk_species_damping_write_disabled;
+  damp->write_fbar_func = gk_species_damping_write_fbar_disabled;
   damp->advance_func = gk_species_damping_advance_disabled;
   damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_disabled;
   damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_disabled;
@@ -479,9 +518,13 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       damp->fbar = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
       damp->fbar1 = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
       damp->fbarnew = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-      damp->fbar_host = damp->fbar;
-      if (app->use_gpu)
-        damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
+      damp->fbar_host = 0;
+      if (damp->write_fbar) {
+        damp->write_fbar_func = gk_species_damping_write_fbar_enabled;
+        damp->fbar_host = damp->fbar;
+        if (app->use_gpu)
+          damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
+      }
 
       // Initialize fbar from the projection of the initial distribution
       // (will be set from the initial f in the main app loop)
@@ -536,7 +579,7 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
       gkyl_array_release(damp->fbar);
       gkyl_array_release(damp->fbar1);
       gkyl_array_release(damp->fbarnew);
-      if (app->use_gpu)
+      if (app->use_gpu && damp->write_fbar)
         gkyl_array_release(damp->fbar_host);
     }
   }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 6e02a15cc3..e88ca74a26 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -338,6 +338,7 @@ struct gkyl_gyrokinetic_damping {
   void (*rate_profile)(double t, const double *xn, double *fout, void *ctx);
   void *rate_profile_ctx; // Context for rate_profile function.
   int num_quad; // Number of quadrature points in each direction to use in projecting the rate.
+  bool write_fbar; // For low-pass filter damping, write fbar diagnostics each output frame.
 };
 
 enum gkyl_gyrokinetic_fdot_multiplier_type {
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index d2bbc9fa5d..8cc8ce2172 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -801,6 +801,7 @@ struct gk_source {
 struct gk_damping {
   enum gkyl_gyrokinetic_damping_type type; // Type of damping term.
   bool evolve; // Whether the source is time dependent.
+  bool write_fbar; // Whether to write low-pass-filter fbar diagnostics.
   struct gkyl_array *rate; // Damping rate.
   struct gkyl_array *rate_host; // Host copy for use in IO and projecting.
   struct gkyl_loss_cone_mask_gyrokinetic *lcm_proj_op; // Operator that projects the loss cone mask.
@@ -814,6 +815,7 @@ struct gk_damping {
   struct gkyl_array *fbar_host; // Host copy of fbar for use in IO.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
+  void (*write_fbar_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*advance_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
     const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
     struct gkyl_array *rhs, struct gkyl_array *cflrate);

From a799125b8c0aca7e8258b5a1c357f4a378eb7172 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 28 Mar 2026 12:38:51 -0400
Subject: [PATCH 12/25] Remove unused parameter from gk_species_damping_reset
 function

---
 gyrokinetic/apps/gk_species_damping.c | 2 --
 1 file changed, 2 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index f20b115443..aabb695f3c 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -591,8 +591,6 @@ void
 gk_species_damping_reset(gkyl_gyrokinetic_app *app, double tm, struct gk_species *gks,
   struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp)
 {
-  (void)tm;
-
   gk_species_damping_release(app, damp);
 
   gks->info.damping = damp_inp;

From 740d6301b271ea3bcbd6b1c705046a00d25cb987 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 28 Mar 2026 17:16:12 -0400
Subject: [PATCH 13/25] Set evolve flag for low-pass filter damping to enable
 time evolution of fbar

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

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index aabb695f3c..7f187ab83f 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -505,6 +505,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
     }
     else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
       // LOW_PASS_FILTER supports either a projected phase-space profile or a uniform constant.
+      damp->evolve = true; // Since fbar must evolve in time.
       if (gks->info.damping.rate_profile) {
         gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
           gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);

From 40a87ecbfc530fdaaca3e52168f68f513eff5bca Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sun, 29 Mar 2026 19:57:42 -0400
Subject: [PATCH 14/25] Add write_rate flag to damping struct and implement
 rate diagnostics in low-pass filter

---
 gyrokinetic/apps/gk_species_damping.c         | 406 +++++++++---------
 gyrokinetic/apps/gkyl_gyrokinetic.h           |   1 +
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h      |   2 +
 .../rt_gk_mirror_boltz_elc_damped_1x2v_p1.c   |   1 +
 gyrokinetic/creg/rt_gk_wham_1x2v_p1.c         |   1 +
 5 files changed, 218 insertions(+), 193 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 7f187ab83f..bee0b7c161 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -36,6 +36,48 @@ gk_species_damping_write_disabled(gkyl_gyrokinetic_app *app, struct gk_species *
 {
 }
 
+static void
+gk_species_damping_write_rate_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+}
+
+static void
+gk_species_damping_write_rate_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+  // DG metadata for damping rate.
+  struct gkyl_msgpack_map_elem mpe_drate[] = {
+    { .key = "poly_order", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = 0 },
+    { .key = "basis_type", .elem_type = GKYL_MP_STRING, .cval = "serendipity" },
+  };
+  int mpe_drate_len = sizeof(mpe_drate) / sizeof(mpe_drate[0]);
+  // Update app basic metadata with time/frame.
+  gkyl_msgpack_map_elem_set_double(app->io_meta_basic_len, app->io_meta_basic, "time", tm);
+  gkyl_msgpack_map_elem_set_uint(app->io_meta_basic_len, app->io_meta_basic, "frame", frame);
+  // Package metadata.
+  int io_meta_len[] = { app->io_meta_basic_len, mpe_drate_len, app->gk_geom->io_meta_len };
+  const struct gkyl_msgpack_map_elem *io_meta[] = { app->io_meta_basic, mpe_drate,
+                                                    app->gk_geom->io_meta };
+  struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union(sizeof(io_meta_len) / sizeof(int),
+    io_meta_len, io_meta);
+
+  // Write out the damping rate.
+  const char *fmt = "%s-%s_damping_rate_%d.gkyl";
+  int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name, frame);
+  char fileNm[sz + 1]; // ensures no buffer overflow
+  snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name, frame);
+
+  // Copy data from device to host before writing it out.
+  if (app->use_gpu)
+    gkyl_array_copy(gks->damping.rate_host, gks->damping.rate);
+
+  gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt, gks->damping.rate_host, fileNm);
+  app->stat.n_io += 1;
+
+  gkyl_msgpack_data_release(mt);
+}
+
 static void
 gk_species_damping_write_fbar_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
   double tm, int frame)
@@ -49,9 +91,9 @@ gk_species_damping_write_fbar_enabled(gkyl_gyrokinetic_app *app, struct gk_speci
   // Metadata from app, species, and geometry (same pattern as species f write).
   gkyl_msgpack_map_elem_set_double(app->io_meta_basic_len, app->io_meta_basic, "time", tm);
   gkyl_msgpack_map_elem_set_uint(app->io_meta_basic_len, app->io_meta_basic, "frame", frame);
-  int io_meta_fbar_len[] = {app->io_meta_basic_len, gks->io_meta_len, app->gk_geom->io_meta_len};
-  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = {app->io_meta_basic, gks->io_meta,
-    app->gk_geom->io_meta};
+  int io_meta_fbar_len[] = { app->io_meta_basic_len, gks->io_meta_len, app->gk_geom->io_meta_len };
+  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = { app->io_meta_basic, gks->io_meta,
+                                                         app->gk_geom->io_meta };
   struct gkyl_msgpack_data *mt_fbar = gkyl_msgpack_create_union(
     sizeof(io_meta_fbar_len) / sizeof(int), io_meta_fbar_len, io_meta_fbar);
 
@@ -76,38 +118,9 @@ gk_species_damping_write_enabled(gkyl_gyrokinetic_app *app, struct gk_species *g
   int frame)
 {
   struct timespec wst = gkyl_wall_clock();
-  // DG metadata for damping rate.
-  struct gkyl_msgpack_map_elem mpe_drate[] = {
-    { .key = "poly_order", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = 0 },
-    { .key = "basis_type", .elem_type = GKYL_MP_STRING, .cval = "serendipity" },
-  };
-  int mpe_drate_len = sizeof(mpe_drate) / sizeof(mpe_drate[0]);
-  // Update app basic metada with time/frame.
-  gkyl_msgpack_map_elem_set_double(app->io_meta_basic_len, app->io_meta_basic, "time", tm);
-  gkyl_msgpack_map_elem_set_uint(app->io_meta_basic_len, app->io_meta_basic, "frame", frame);
-  // Package metadata.
-  int io_meta_len[] = { app->io_meta_basic_len, mpe_drate_len, app->gk_geom->io_meta_len };
-  const struct gkyl_msgpack_map_elem *io_meta[] = { app->io_meta_basic, mpe_drate,
-                                                    app->gk_geom->io_meta };
-  struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union(sizeof(io_meta_len) / sizeof(int),
-    io_meta_len, io_meta);
-
-  // Write out the damping rate.
-  const char *fmt = "%s-%s_damping_rate_%d.gkyl";
-  int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name, frame);
-  char fileNm[sz + 1]; // ensures no buffer overflow
-  snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name, frame);
-
-  // Copy data from device to host before writing it out.
-  if (app->use_gpu)
-    gkyl_array_copy(gks->damping.rate_host, gks->damping.rate);
-
-  gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt, gks->damping.rate_host, fileNm);
-  app->stat.n_io += 1;
+  gks->damping.write_rate_func(app, gks, tm, frame);
 
   gks->damping.write_fbar_func(app, gks, tm, frame);
-
-  gkyl_msgpack_data_release(mt);
   app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wst);
 }
 
@@ -355,6 +368,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 {
   damp->type = gks->info.damping.type;
   damp->evolve = false; // Whether the rate is time dependent.
+  damp->write_rate = gks->info.damping.write_rate;
   damp->write_fbar = gks->info.damping.write_fbar;
   const double rate_const = gks->info.damping.rate_const;
 
@@ -363,6 +377,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
   // Default function pointers.
   damp->write_func = gk_species_damping_write_disabled;
+  damp->write_rate_func = gk_species_damping_write_rate_disabled;
   damp->write_fbar_func = gk_species_damping_write_fbar_disabled;
   damp->advance_func = gk_species_damping_advance_disabled;
   damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_disabled;
@@ -371,181 +386,186 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->combine_func = gk_species_damping_combine_disabled;
   damp->copy_func = gk_species_damping_copy_range_disabled;
 
-  if (damp->type) {
-    // Allocate rate array.
-    damp->rate = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
-      gks->local_ext.volume);
-    damp->rate_host = damp->rate;
-    if (app->use_gpu)
-      damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size);
+  if (!damp->type) {
+    return;
+  }
 
-    if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
-      // USER_INPUT supports either a projected profile, or a uniform constant if no profile is provided.
-      if (gks->info.damping.rate_profile) {
-        gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
-          gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
-        gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
-      }
-      else {
-        gkyl_array_clear(damp->rate, rate_const);
-      }
+  if (damp->write_rate)
+    damp->write_rate_func = gk_species_damping_write_rate_enabled;
 
-      damp->advance_func = gk_species_damping_advance_user_input;
+  // Allocate rate array.
+  damp->rate = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
+    gks->local_ext.volume);
+  damp->rate_host = damp->rate;
+  if (app->use_gpu)
+    damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size);
+
+  if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
+    // USER_INPUT supports either a projected profile, or a uniform constant if no profile is provided.
+    if (gks->info.damping.rate_profile) {
+      gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
+        gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
+      gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
+    }
+    else {
+      gkyl_array_clear(damp->rate, rate_const);
     }
-    else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
-      damp->evolve = true; // Since the loss cone boundary is proportional to phi(t).
-
-      // Maximum bmag and its location.
-      // NOTE: if the same max bmag occurs at multiple locations,
-      // bmag_max_coord may have different values on different MPI processes.
-      double bmag_max_coord_ho[GKYL_MAX_CDIM];
-      double bmag_max_ho = gkyl_gk_geometry_reduce_arg_bmag(app->gk_geom, GKYL_MAX,
-        bmag_max_coord_ho);
-      double bmag_max_local = bmag_max_ho;
-      double bmag_max_global;
-      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, &bmag_max_local,
-        &bmag_max_global);
-      double bmag_max_coord_local[app->cdim], bmag_max_coord_global[app->cdim];
-      if (fabs(bmag_max_ho - bmag_max_global) < 1e-16) {
-        for (int d = 0; d < app->cdim; d++) {
-          bmag_max_coord_local[d] = bmag_max_coord_ho[d];
-        }
-      }
-      else {
-        for (int d = 0; d < app->cdim; d++) {
-          bmag_max_coord_local[d] = -DBL_MAX;
-        }
-      }
-      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, app->cdim, bmag_max_coord_local,
-        bmag_max_coord_global);
-
-      if (app->use_gpu) {
-        damp->bmag_max = gkyl_cu_malloc(sizeof(double));
-        damp->bmag_max_coord = gkyl_cu_malloc(app->cdim * sizeof(double));
-        gkyl_cu_memcpy(damp->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
-        gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double),
-          GKYL_CU_MEMCPY_H2D);
-      }
-      else {
-        damp->bmag_max = gkyl_malloc(sizeof(double));
-        damp->bmag_max_coord = gkyl_malloc(app->cdim * sizeof(double));
-        memcpy(damp->bmag_max, &bmag_max_global, sizeof(double));
-        memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double));
-      }
-
-      // Electrostatic potential at bmag_max_coord.
-      if (app->use_gpu) {
-        damp->phi_m = gkyl_cu_malloc(sizeof(double));
-        damp->phi_m_global = gkyl_cu_malloc(sizeof(double));
-      }
-      else {
-        damp->phi_m = gkyl_malloc(sizeof(double));
-        damp->phi_m_global = gkyl_malloc(sizeof(double));
-      }
 
-      // Operator that projects the loss cone mask.
-      struct gkyl_loss_cone_mask_gyrokinetic_inp inp_proj = {
-        .phase_grid = &gks->grid,
-        .conf_basis = &app->basis,
-        .phase_basis = &gks->basis,
-        .conf_range = &app->local,
-        .conf_range_ext = &app->local_ext,
-        .vel_range = &gks->local_vel,
-        .vel_map = gks->vel_map,
-        .bmag = app->gk_geom->geo_int.bmag,
-        .bmag_max = damp->bmag_max,
-        .bmag_max_loc = damp->bmag_max_coord,
-        .mass = gks->info.mass,
-        .charge = gks->info.charge,
-        .num_quad = num_quad,
-        .use_gpu = app->use_gpu,
-      };
-      damp->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew(&inp_proj);
-
-      // Project the conf-space rate profile provided.
-      struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, gks->basis.num_basis,
-        gks->local_ext.volume);
-      struct gkyl_array *scale_prof_high_order_ho = app->use_gpu? mkarr(false,
-        scale_prof_high_order->ncomp, scale_prof_high_order->size)
-                                                     : gkyl_array_acquire(scale_prof_high_order);
-
-      if (gks->info.damping.rate_profile) {
-        gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1,
-          gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
-        gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
-        gkyl_proj_on_basis_release(projup);
-        gkyl_array_copy(scale_prof_high_order, scale_prof_high_order_ho);
+    damp->advance_func = gk_species_damping_advance_user_input;
+  }
+  else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
+    damp->evolve = true; // Since the loss cone boundary is proportional to phi(t).
+
+    // Maximum bmag and its location.
+    // NOTE: if the same max bmag occurs at multiple locations,
+    // bmag_max_coord may have different values on different MPI processes.
+    double bmag_max_coord_ho[GKYL_MAX_CDIM];
+    double bmag_max_ho = gkyl_gk_geometry_reduce_arg_bmag(app->gk_geom, GKYL_MAX,
+      bmag_max_coord_ho);
+    double bmag_max_local = bmag_max_ho;
+    double bmag_max_global;
+    gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, &bmag_max_local,
+      &bmag_max_global);
+    double bmag_max_coord_local[app->cdim], bmag_max_coord_global[app->cdim];
+    if (fabs(bmag_max_ho - bmag_max_global) < 1e-16) {
+      for (int d = 0; d < app->cdim; d++) {
+        bmag_max_coord_local[d] = bmag_max_coord_ho[d];
       }
-      else {
-        gkyl_array_clear(scale_prof_high_order, 1.0);
+    }
+    else {
+      for (int d = 0; d < app->cdim; d++) {
+        bmag_max_coord_local[d] = -DBL_MAX;
       }
+    }
+    gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, app->cdim, bmag_max_coord_local,
+      bmag_max_coord_global);
+
+    if (app->use_gpu) {
+      damp->bmag_max = gkyl_cu_malloc(sizeof(double));
+      damp->bmag_max_coord = gkyl_cu_malloc(app->cdim * sizeof(double));
+      gkyl_cu_memcpy(damp->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
+      gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double),
+        GKYL_CU_MEMCPY_H2D);
+    }
+    else {
+      damp->bmag_max = gkyl_malloc(sizeof(double));
+      damp->bmag_max_coord = gkyl_malloc(app->cdim * sizeof(double));
+      memcpy(damp->bmag_max, &bmag_max_global, sizeof(double));
+      memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double));
+    }
 
-      gkyl_array_scale(scale_prof_high_order, rate_const == 0.0 ? 1.0 : rate_const);
-
-      damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
-        gks->local_ext.volume);
-      gkyl_array_set_offset(damp->scale_prof, pow(sqrt(2.0), gks->grid.ndim), scale_prof_high_order,
-        0);
+    // Electrostatic potential at bmag_max_coord.
+    if (app->use_gpu) {
+      damp->phi_m = gkyl_cu_malloc(sizeof(double));
+      damp->phi_m_global = gkyl_cu_malloc(sizeof(double));
+    }
+    else {
+      damp->phi_m = gkyl_malloc(sizeof(double));
+      damp->phi_m_global = gkyl_malloc(sizeof(double));
+    }
 
-      gkyl_array_release(scale_prof_high_order_ho);
-      gkyl_array_release(scale_prof_high_order);
+    // Operator that projects the loss cone mask.
+    struct gkyl_loss_cone_mask_gyrokinetic_inp inp_proj = {
+      .phase_grid = &gks->grid,
+      .conf_basis = &app->basis,
+      .phase_basis = &gks->basis,
+      .conf_range = &app->local,
+      .conf_range_ext = &app->local_ext,
+      .vel_range = &gks->local_vel,
+      .vel_map = gks->vel_map,
+      .bmag = app->gk_geom->geo_int.bmag,
+      .bmag_max = damp->bmag_max,
+      .bmag_max_loc = damp->bmag_max_coord,
+      .mass = gks->info.mass,
+      .charge = gks->info.charge,
+      .num_quad = num_quad,
+      .use_gpu = app->use_gpu,
+    };
+    damp->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew(&inp_proj);
+
+    // Project the conf-space rate profile provided.
+    struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, gks->basis.num_basis,
+      gks->local_ext.volume);
+    struct gkyl_array *scale_prof_high_order_ho = app->use_gpu? mkarr(false,
+      scale_prof_high_order->ncomp, scale_prof_high_order->size)
+                                                    : gkyl_array_acquire(scale_prof_high_order);
+
+    if (gks->info.damping.rate_profile) {
+      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1,
+        gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
+      gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
+      gkyl_proj_on_basis_release(projup);
+      gkyl_array_copy(scale_prof_high_order, scale_prof_high_order_ho);
+    }
+    else {
+      gkyl_array_clear(scale_prof_high_order, 1.0);
+    }
 
-      // Compute the initial damping rate (assuming phi=0 because phi hasn't been computed).
-      // Find the potential at the mirror throat.
-      gkyl_dg_basis_ops_eval_array_at_coord_comp(app->field->phi_smooth, damp->bmag_max_coord,
-        app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
-      gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
-      // Project the loss cone mask.
-      gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
-        app->field->phi_smooth, damp->phi_m_global, damp->rate);
-      // Multiply by the user's scaling profile.
-      gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
+    gkyl_array_scale(scale_prof_high_order, rate_const == 0.0 ? 1.0 : rate_const);
 
-      damp->advance_func = gk_species_damping_advance_loss_cone;
+    damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
+      gks->local_ext.volume);
+    gkyl_array_set_offset(damp->scale_prof, pow(sqrt(2.0), gks->grid.ndim), scale_prof_high_order,
+      0);
+
+    gkyl_array_release(scale_prof_high_order_ho);
+    gkyl_array_release(scale_prof_high_order);
+
+    // Compute the initial damping rate (assuming phi=0 because phi hasn't been computed).
+    // Find the potential at the mirror throat.
+    gkyl_dg_basis_ops_eval_array_at_coord_comp(app->field->phi_smooth, damp->bmag_max_coord,
+      app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
+    gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
+    // Project the loss cone mask.
+    gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
+      app->field->phi_smooth, damp->phi_m_global, damp->rate);
+    // Multiply by the user's scaling profile.
+    gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
+
+    damp->advance_func = gk_species_damping_advance_loss_cone;
+  }
+  else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+    // LOW_PASS_FILTER supports either a projected phase-space profile or a uniform constant.
+    damp->evolve = true; // Since fbar must evolve in time.
+    if (gks->info.damping.rate_profile) {
+      gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
+        gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
+      gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
+    }
+    else {
+      gkyl_array_clear(damp->rate, rate_const);
     }
-    else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-      // LOW_PASS_FILTER supports either a projected phase-space profile or a uniform constant.
-      damp->evolve = true; // Since fbar must evolve in time.
-      if (gks->info.damping.rate_profile) {
-        gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
-          gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
-        gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
-      }
-      else {
-        gkyl_array_clear(damp->rate, rate_const);
-      }
 
-      // Allocate filtered distribution function array
-      damp->fbar = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-      damp->fbar1 = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-      damp->fbarnew = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-      damp->fbar_host = 0;
-      if (damp->write_fbar) {
-        damp->write_fbar_func = gk_species_damping_write_fbar_enabled;
-        damp->fbar_host = damp->fbar;
-        if (app->use_gpu)
-          damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
-      }
+    // Allocate filtered distribution function array
+    damp->fbar = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
+    damp->fbar1 = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
+    damp->fbarnew = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
+    damp->fbar_host = 0;
+    if (damp->write_fbar) {
+      damp->write_fbar_func = gk_species_damping_write_fbar_enabled;
+      damp->fbar_host = damp->fbar;
+      if (app->use_gpu)
+        damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
+    }
 
-      // Initialize fbar from the projection of the initial distribution
-      // (will be set from the initial f in the main app loop)
-      gkyl_array_clear(damp->fbar, 0.0);
+    // Initialize fbar from the projection of the initial distribution
+    // (will be set from the initial f in the main app loop)
+    gkyl_array_clear(damp->fbar, 0.0);
 
-      damp->advance_func = gk_species_damping_advance_low_pass_filter;
-      damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_enabled;
-      damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_enabled;
-      damp->forward_euler_func = gk_species_damping_forward_euler_enabled;
-      damp->combine_func = gk_species_damping_combine_enabled;
-      damp->copy_func = gk_species_damping_copy_range_enabled;
-    }
+    damp->advance_func = gk_species_damping_advance_low_pass_filter;
+    damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_enabled;
+    damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_enabled;
+    damp->forward_euler_func = gk_species_damping_forward_euler_enabled;
+    damp->combine_func = gk_species_damping_combine_enabled;
+    damp->copy_func = gk_species_damping_copy_range_enabled;
+  }
 
-    // Set function pointers chosen at runtime.
-    if (damp->evolve) {
-      damp->write_func = gk_species_damping_write_enabled;
-    }
-    else {
-      damp->write_func = gk_species_damping_write_init_only;
-    }
+  // Set function pointers chosen at runtime.
+  if (damp->evolve) {
+    damp->write_func = gk_species_damping_write_enabled;
+  }
+  else {
+    damp->write_func = gk_species_damping_write_init_only;
   }
 }
 
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index e88ca74a26..60c6b26f22 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -338,6 +338,7 @@ struct gkyl_gyrokinetic_damping {
   void (*rate_profile)(double t, const double *xn, double *fout, void *ctx);
   void *rate_profile_ctx; // Context for rate_profile function.
   int num_quad; // Number of quadrature points in each direction to use in projecting the rate.
+  bool write_rate; // Whether to write damping-rate diagnostics.
   bool write_fbar; // For low-pass filter damping, write fbar diagnostics each output frame.
 };
 
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 8cc8ce2172..ca2967986a 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -801,6 +801,7 @@ struct gk_source {
 struct gk_damping {
   enum gkyl_gyrokinetic_damping_type type; // Type of damping term.
   bool evolve; // Whether the source is time dependent.
+  bool write_rate; // Whether to write damping-rate diagnostics.
   bool write_fbar; // Whether to write low-pass-filter fbar diagnostics.
   struct gkyl_array *rate; // Damping rate.
   struct gkyl_array *rate_host; // Host copy for use in IO and projecting.
@@ -815,6 +816,7 @@ struct gk_damping {
   struct gkyl_array *fbar_host; // Host copy of fbar for use in IO.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
+  void (*write_rate_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*write_fbar_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*advance_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
     const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
diff --git a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
index d6f94e1324..4729b5274e 100644
--- a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
@@ -710,6 +710,7 @@ int main(int argc, char **argv)
       .type = GKYL_GK_DAMPING_LOSS_CONE,
       .rate_profile = loss_cone_damping_rate_scaling,
       .rate_profile_ctx = &ctx,
+      .write_rate = true,
     },
 
     .bcs = {
diff --git a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
index ab9391aca2..c42332086f 100644
--- a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
@@ -611,6 +611,7 @@ int main(int argc, char **argv)
     .damping = {
       .type = GKYL_GK_DAMPING_LOW_PASS_FILTER,
       .rate_const = 1e-6,
+      .write_rate = true,
     },
 
     .source = {

From d4be6d140a64921a44284103b8720c6badbc75ee Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 12:08:56 -0400
Subject: [PATCH 15/25] Change the fbar distribution (low passed) that we
 update to P0 to save a ton of memory and computation on the updates

---
 gyrokinetic/apps/gk_species_damping.c         | 66 +++++++++++--------
 .../rt_gk_mirror_boltz_elc_damped_1x2v_p1.c   |  3 -
 gyrokinetic/creg/rt_gk_wham_1x2v_p1.c         |  3 +-
 3 files changed, 42 insertions(+), 30 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index bee0b7c161..eb16106e62 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -10,9 +10,9 @@ void
 gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk_damping *damp,
   const struct gkyl_array *f)
 {
-  gkyl_array_set(damp->fbar, 1.0, f);
-  gkyl_array_set(damp->fbar1, 1.0, f);
-  gkyl_array_set(damp->fbarnew, 1.0, f);
+  gkyl_array_set_offset(damp->fbar, 1.0, f, 0);
+  gkyl_array_set_offset(damp->fbar1, 1.0, f, 0);
+  gkyl_array_set_offset(damp->fbarnew, 1.0, f, 0);
 }
 
 void
@@ -88,16 +88,28 @@ static void
 gk_species_damping_write_fbar_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
   double tm, int frame)
 {
-  // Metadata from app, species, and geometry (same pattern as species f write).
+  // Metadata from app/species/geometry. fbar is always stored as a p0 field.
   gkyl_msgpack_map_elem_set_double(app->io_meta_basic_len, app->io_meta_basic, "time", tm);
   gkyl_msgpack_map_elem_set_uint(app->io_meta_basic_len, app->io_meta_basic, "frame", frame);
-  int io_meta_fbar_len[] = { app->io_meta_basic_len, gks->io_meta_len, app->gk_geom->io_meta_len };
-  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = { app->io_meta_basic, gks->io_meta,
-                                                         app->gk_geom->io_meta };
+
+  struct gkyl_msgpack_map_elem mpe_fbar_p0[] = {
+    { .key = "poly_order", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = 0 },
+    { .key = "basis_type", .elem_type = GKYL_MP_STRING, .cval = "serendipity" },
+  };
+  int io_meta_fbar_len[] = {
+    app->io_meta_basic_len,
+    (int)(sizeof(mpe_fbar_p0) / sizeof(mpe_fbar_p0[0])),
+    app->gk_geom->io_meta_len
+  };
+  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = {
+    app->io_meta_basic,
+    mpe_fbar_p0,
+    app->gk_geom->io_meta
+  };
   struct gkyl_msgpack_data *mt_fbar = gkyl_msgpack_create_union(
     sizeof(io_meta_fbar_len) / sizeof(int), io_meta_fbar_len, io_meta_fbar);
 
-  const char *fmt_fbar = "%s-%s_fbar_%d.gkyl";
+  const char *fmt_fbar = "%s-%s_damping_fbar_%d.gkyl";
   int sz_fbar = gkyl_calc_strlen(fmt_fbar, app->name, gks->info.name, frame);
   char fileNm_fbar[sz_fbar + 1]; // ensures no buffer overflow
   snprintf(fileNm_fbar, sizeof fileNm_fbar, fmt_fbar, app->name, gks->info.name, frame);
@@ -153,9 +165,12 @@ gk_species_damping_project_phase_rate(const struct gkyl_gyrokinetic_app *app,
   proj_on_basis_c2p_ctx.vdim = gks->local_vel.ndim;
   proj_on_basis_c2p_ctx.vel_map = gks->vel_map;
 
+  struct gkyl_basis basis_mult;
+  gkyl_cart_modal_serendip(&basis_mult, gks->basis.ndim, 0);
+
   gkyl_proj_on_basis *projup = gkyl_proj_on_basis_inew(&(struct gkyl_proj_on_basis_inp) {
     .grid = &gks->grid,
-    .basis = &gks->basis,
+    .basis = &basis_mult,
     .num_quad = num_quad,
     .num_ret_vals = 1,
     .eval = rate_profile,
@@ -233,10 +248,10 @@ gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const stru
   else if (fin == gks->fnew)
     fbar_in = damp->fbarnew;
 
-  // Compute f - fbar and scale by the damping rate: rate * (f - fbar)
-  gkyl_array_set(f_buffer, 1.0, fin);      // f_buffer = f
-  gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = f - fbar
-  gkyl_array_scale_by_cell(f_buffer, damp->rate);    // f_buffer = rate * (f - fbar)
+  // Compute f - fbar, where fbar is stored as a cellwise-constant p0 field.
+  gkyl_array_set(f_buffer, 1.0, fin); // f_buffer = fin
+  gkyl_array_accumulate_offset(f_buffer, -1.0, fbar_in, 0); // f_buffer = fin - fbar
+  gkyl_array_scale_by_cell(f_buffer, damp->rate); // f_buffer = rate * (fin - fbar)
 
   // Add damping term to RHS: df/dt -= rate * (f - fbar)
   // Add to the CFL frequency.
@@ -276,7 +291,7 @@ gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
 {
   // rhs_fbar = rate * (f - fbar)
-  gkyl_array_set(rhs_fbar, 1.0, fin);
+  gkyl_array_set_offset(rhs_fbar, 1.0, fin, 0);
   gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
   gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
 }
@@ -302,8 +317,7 @@ gk_species_damping_forward_euler_enabled(struct gk_species *gks,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
   struct gkyl_array *fbar_out, double dt)
 {
-  // Reuse species scratch storage to assemble fbar RHS for this FE substep.
-  gk_species_damping_calc_fbar_rhs(&gks->damping, fin, fbar_in, gks->lte.f_lte);
+  gk_species_damping_calc_fbar_rhs(&gks->damping, fin, fbar_in, fbar_out);
   gk_species_step_f(gks, fbar_out, dt, fbar_in);
 }
 
@@ -394,8 +408,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
     damp->write_rate_func = gk_species_damping_write_rate_enabled;
 
   // Allocate rate array.
-  damp->rate = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
-    gks->local_ext.volume);
+  damp->rate = mkarr(app->use_gpu, 1, gks->local_ext.volume);
   damp->rate_host = damp->rate;
   if (app->use_gpu)
     damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size);
@@ -484,14 +497,16 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
     damp->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew(&inp_proj);
 
     // Project the conf-space rate profile provided.
-    struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, gks->basis.num_basis,
+    struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, 1,
       gks->local_ext.volume);
     struct gkyl_array *scale_prof_high_order_ho = app->use_gpu? mkarr(false,
       scale_prof_high_order->ncomp, scale_prof_high_order->size)
                                                     : gkyl_array_acquire(scale_prof_high_order);
 
     if (gks->info.damping.rate_profile) {
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1,
+      struct gkyl_basis basis_mult;
+      gkyl_cart_modal_serendip(&basis_mult, gks->basis.ndim, 0);
+      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &basis_mult, num_quad, 1,
         gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
       gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
       gkyl_proj_on_basis_release(projup);
@@ -503,8 +518,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
     gkyl_array_scale(scale_prof_high_order, rate_const == 0.0 ? 1.0 : rate_const);
 
-    damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis,
-      gks->local_ext.volume);
+    damp->scale_prof = mkarr(app->use_gpu, 1, gks->local_ext.volume);
     gkyl_array_set_offset(damp->scale_prof, pow(sqrt(2.0), gks->grid.ndim), scale_prof_high_order,
       0);
 
@@ -536,10 +550,10 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       gkyl_array_clear(damp->rate, rate_const);
     }
 
-    // Allocate filtered distribution function array
-    damp->fbar = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-    damp->fbar1 = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-    damp->fbarnew = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
+    // Allocate filtered distribution function array as p0.
+    damp->fbar = mkarr(app->use_gpu, 1, gks->local_ext.volume);
+    damp->fbar1 = mkarr(app->use_gpu, 1, gks->local_ext.volume);
+    damp->fbarnew = mkarr(app->use_gpu, 1, gks->local_ext.volume);
     damp->fbar_host = 0;
     if (damp->write_fbar) {
       damp->write_fbar_func = gk_species_damping_write_fbar_enabled;
diff --git a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
index 4729b5274e..a50ff1162d 100644
--- a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
@@ -704,9 +704,6 @@ int main(int argc, char **argv)
     },
 
     .damping = {
-//      .type = GKYL_GK_DAMPING_USER_INPUT,
-//      .rate_profile = loss_cone_damping_rate_profile,
-//      .rate_profile_ctx = &ctx,
       .type = GKYL_GK_DAMPING_LOSS_CONE,
       .rate_profile = loss_cone_damping_rate_scaling,
       .rate_profile_ctx = &ctx,
diff --git a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
index c42332086f..664643e299 100644
--- a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
@@ -465,7 +465,7 @@ create_ctx(void)
   int poly_order = 1;
 
   double t_end = 1e-8;
-  int num_frames = 10;
+  int num_frames = 1;
   double write_phase_freq = 0.2; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
   int int_diag_calc_num = num_frames*100;
   double dt_failure_tol = 1.0e-4; // Minimum allowable fraction of initial time-step.
@@ -612,6 +612,7 @@ int main(int argc, char **argv)
       .type = GKYL_GK_DAMPING_LOW_PASS_FILTER,
       .rate_const = 1e-6,
       .write_rate = true,
+      .write_fbar = true,
     },
 
     .source = {

From df74b33697976bade7f2924012c49feee13c8604 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 12:23:23 -0400
Subject: [PATCH 16/25] There was some really crazy spacing issues that look
 like it might come from someone using Vim in the
 gyrokinetic_multib_update_ssp_rk3. I fixed the spacing issues. I also found
 an issue with the implementation where the fbar was not being passed
 appropriately from ssprk3 down to the damping modules. Syntax is updated

---
 gyrokinetic/apps/gk_species.c                 | 13 ++--
 gyrokinetic/apps/gk_species_damping.c         | 24 +++----
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h      | 18 +++--
 gyrokinetic/apps/gyrokinetic.c                |  5 +-
 .../apps/gyrokinetic_multib_update_ssp_rk3.c  | 70 ++++++++++++-------
 gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c |  3 +-
 6 files changed, 80 insertions(+), 53 deletions(-)

diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c
index d8c042ff98..0bf3791865 100644
--- a/gyrokinetic/apps/gk_species.c
+++ b/gyrokinetic/apps/gk_species.c
@@ -75,7 +75,8 @@ gk_species_omegaH_dt(gkyl_gyrokinetic_app *app, struct gk_species *gks, const st
 
 static double
 gk_species_rhs_dynamic(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
 {
   // Gyroaverage the potential if needed.
   species->gyroaverage(app, species, app->field->phi_smooth, species->gyro_phi);
@@ -88,7 +89,7 @@ gk_species_rhs_dynamic(gkyl_gyrokinetic_app *app, struct gk_species *species,
 
   // Damping term.
   gk_species_damping_advance(app, species, &species->damping, app->field->phi_smooth, fin,
-    species->lte.f_lte, rhs, species->cflrate);
+    fbar_in, species->lte.f_lte, rhs, species->cflrate);
 
   // LBO Collisions.
   gk_species_lbo_rhs(app, species, &species->lbo, fin, rhs);
@@ -174,7 +175,8 @@ gk_species_rhs_implicit_dynamic(gkyl_gyrokinetic_app *app, struct gk_species *sp
 
 static double
 gk_species_rhs_static(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
 {
   double omega_cfl = 1/DBL_MAX;
   return app->cfl/omega_cfl;
@@ -1782,9 +1784,10 @@ gk_species_apply_ic_cross(gkyl_gyrokinetic_app *app, struct gk_species *gks_self
 
 double
 gk_species_rhs(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
 {
-  return species->rhs_func(app, species, fin, rhs, bflux_moms);
+  return species->rhs_func(app, species, fin, fbar_in, rhs, bflux_moms);
 }
 
 double
diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index eb16106e62..347566a6d5 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -191,7 +191,8 @@ gk_species_damping_project_phase_rate(const struct gkyl_gyrokinetic_app *app,
 void
 gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
 }
@@ -199,7 +200,8 @@ gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_s
 void
 gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
   gkyl_array_set(f_buffer, 1.0, fin);
@@ -213,7 +215,8 @@ gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk
 void
 gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
   // Find the potential at the mirror throat.
@@ -238,16 +241,10 @@ gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_
 void
 gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
-  // Match fbar state to the RK stage of fin.
-  const struct gkyl_array *fbar_in = damp->fbar;
-  if (fin == gks->f1)
-    fbar_in = damp->fbar1;
-  else if (fin == gks->fnew)
-    fbar_in = damp->fbarnew;
-
   // Compute f - fbar, where fbar is stored as a cellwise-constant p0 field.
   gkyl_array_set(f_buffer, 1.0, fin); // f_buffer = fin
   gkyl_array_accumulate_offset(f_buffer, -1.0, fbar_in, 0); // f_buffer = fin - fbar
@@ -262,12 +259,13 @@ gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const stru
 void
 gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate)
 {
   struct timespec wst = gkyl_wall_clock();
 
-  damp->advance_func(app, gks, damp, phi, fin, f_buffer, rhs, cflrate);
+  damp->advance_func(app, gks, damp, phi, fin, fbar_in, f_buffer, rhs, cflrate);
 
   app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index ca2967986a..2ec502bfe6 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -819,7 +819,8 @@ struct gk_damping {
   void (*write_rate_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*write_fbar_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*advance_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
-    const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+    const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+    struct gkyl_array *f_buffer,
     struct gkyl_array *rhs, struct gkyl_array *cflrate);
   void (*set_fbar_to_f_func)(const struct gk_species *gks, struct gk_damping *damp,
     const struct gkyl_array *f);
@@ -1049,7 +1050,8 @@ struct gk_species {
 
   // Pointer to various functions selected at runtime.
   double (*rhs_func)(gkyl_gyrokinetic_app *app, struct gk_species *species,
-    const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
+    const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+    struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
   double (*rhs_implicit_func)(gkyl_gyrokinetic_app *app, struct gk_species *species,
     const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms, double dt);
   void (*bc_func)(gkyl_gyrokinetic_app *app, const struct gk_species *species,
@@ -2813,12 +2815,14 @@ void gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_da
  * @param damp Species damping object.
  * @param phi Current electrostatic potential.
  * @param fin Current distribution function.
+ * @param fbar_in Filtered distribution corresponding to fin RK stage.
  * @param f_buffer Phase-space buffer.
  * @param rhs df/dt damping term gets added to.
  * @param cflrate CFL frequency in phase space.
  */
 void gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate);
 
 /**
@@ -3091,12 +3095,14 @@ void gk_species_apply_ic_cross(gkyl_gyrokinetic_app *app, struct gk_species *spe
  * @param app gyrokinetic app object.
  * @param species Pointer to species.
  * @param fin Input distribution function.
+ * @param fbar_in Input filtered distribution for the current RK stage.
  * @param rhs On output, the RHS from the species object.
  * @param bflux_moms Output boundary flux moments.
  * @return Maximum stable time-step.
  */
 double gk_species_rhs(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
 
 /**
  * Compute the *implicit* RHS from species distribution function
@@ -4056,6 +4062,7 @@ void gyrokinetic_calc_field_and_apply_bc(gkyl_gyrokinetic_app* app, double tcurr
  * @param tcurr Current simulation time.
  * @param dt Suggested time step.
  * @param fin Input array of charged-species distribution functions.
+ * @param fbar_in Input array of charged-species filtered distributions.
  * @param fout Output array of charged-species distribution functions.
  * @param bflux_out Output array of charged-species boundary fluxes.
  * @param fin_neut Input array of neutral-species distribution functions.
@@ -4064,7 +4071,8 @@ void gyrokinetic_calc_field_and_apply_bc(gkyl_gyrokinetic_app* app, double tcurr
  * @param st Time stepping status object.
  */
 void gyrokinetic_rhs(gkyl_gyrokinetic_app* app, double tcurr, double dt,
-  const struct gkyl_array *fin[], struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
+  const struct gkyl_array *fin[], const struct gkyl_array *fbar_in[],
+  struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], struct gkyl_array **bflux_out_neut[],
   struct gkyl_update_status *st); 
 
diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c
index b728aaef51..6e3bd13f58 100644
--- a/gyrokinetic/apps/gyrokinetic.c
+++ b/gyrokinetic/apps/gyrokinetic.c
@@ -1801,7 +1801,8 @@ gkyl_gyrokinetic_app_write(gkyl_gyrokinetic_app* app, double tm, int frame)
 
 void
 gyrokinetic_rhs(gkyl_gyrokinetic_app* app, double tcurr, double dt,
-  const struct gkyl_array *fin[], struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
+  const struct gkyl_array *fin[], const struct gkyl_array *fbar_in[],
+  struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], struct gkyl_array **bflux_out_neut[], 
   struct gkyl_update_status *st)
 {
@@ -1846,7 +1847,7 @@ gyrokinetic_rhs(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   // Compute df/dt (not including sources).
   for (int i=0; i<app->num_species; ++i) {
     struct gk_species *gk_s = &app->species[i];
-    double dt1 = gk_species_rhs(app, gk_s, fin[i], fout[i], bflux_out[i]);
+    double dt1 = gk_species_rhs(app, gk_s, fin[i], fbar_in[i], fout[i], bflux_out[i]);
     dtmin = fmin(dtmin, dt1);
   }
   for (int i=0; i<app->num_neut_species; ++i) {
diff --git a/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c b/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c
index 8f079ac882..fa2e1a8839 100644
--- a/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c
+++ b/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c
@@ -2,7 +2,8 @@
 
 static void
 gyrokinetic_multib_forward_euler(struct gkyl_gyrokinetic_multib_app* app, double tcurr, double dt,
-  const struct gkyl_array *fin[], struct gkyl_array *fout[], 
+  const struct gkyl_array *fin[], struct gkyl_array *fout[],
+  const struct gkyl_array *fbar_in[], struct gkyl_array *fbar_out[],
   struct gkyl_array **bflux_in[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], 
   struct gkyl_array **bflux_in_neut[], struct gkyl_array **bflux_out_neut[], 
@@ -22,7 +23,7 @@ gyrokinetic_multib_forward_euler(struct gkyl_gyrokinetic_multib_app* app, double
   for (int b=0; b<app->num_local_blocks; ++b) {
     int li_charged = b * app->num_species;
     int li_neut = b * app->num_neut_species;
-    gyrokinetic_rhs(app->singleb_apps[b], tcurr, dt, &fin[li_charged], &fout[li_charged],
+    gyrokinetic_rhs(app->singleb_apps[b], tcurr, dt, &fin[li_charged], &fbar_in[li_charged], &fout[li_charged],
       &bflux_out[li_charged], &fin_neut[li_neut], &fout_neut[li_neut], &bflux_out_neut[li_neut], st);
     dtmin = fmin(dtmin, st->dt_actual);
   }
@@ -44,6 +45,7 @@ gyrokinetic_multib_forward_euler(struct gkyl_gyrokinetic_multib_app* app, double
     for (int i=0; i<app->num_species; ++i) {
       struct gk_species *gks = &sbapp->species[i];
       gk_species_step_f(gks, fout[li_charged+i], dta, fin[li_charged+i]);
+      gk_species_damping_forward_euler(gks, fin[li_charged+i], fbar_in[li_charged+i], fbar_out[li_charged+i], dta);
       gk_species_bflux_step_f(sbapp, &gks->bflux, bflux_out[li_charged+i], 1.0, bflux_in[li_charged+i]);
     }
     for (int i=0; i<app->num_neut_species; ++i) {
@@ -69,6 +71,8 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
 
   const struct gkyl_array *fin[ns_charged * nblocks_local];
   struct gkyl_array *fout[ns_charged * nblocks_local];
+  const struct gkyl_array *fbar_in[ns_charged * nblocks_local];
+  struct gkyl_array *fbar_out[ns_charged * nblocks_local];
   struct gkyl_array **bflux_in[ns_charged * nblocks_local];
   struct gkyl_array **bflux_out[ns_charged * nblocks_local];
 
@@ -94,6 +98,8 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
             struct gk_species *gks = &sbapp->species[i];
             fin[li_charged+i] = gks->f;
             fout[li_charged+i] = gks->f1;
+            fbar_in[li_charged+i] = gks->damping.fbar;
+            fbar_out[li_charged+i] = gks->damping.fbar1;
             // Boundary fluxes.
             bflux_in[li_charged+i] = gks->bflux.f;
             bflux_out[li_charged+i] = gks->bflux.f1;
@@ -101,14 +107,14 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           for (int i=0; i<ns_neut; ++i) {
             struct gk_neut_species *gkns = &sbapp->neut_species[i];
             fin_neut[li_neut+i] = gkns->f;
-	    fout_neut[li_neut+i] = gkns->f1;
+            fout_neut[li_neut+i] = gkns->f1;
             // Boundary fluxes.
             bflux_in_neut[li_neut+i] = gkns->bflux.f;
             bflux_out_neut[li_neut+i] = gkns->bflux.f1;
           }
         }
 
-        gyrokinetic_multib_forward_euler(app, tcurr, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_multib_forward_euler(app, tcurr, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
         dt = st.dt_actual;
 
@@ -121,11 +127,11 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           int li_charged = b * ns_charged;
           int li_neut = b * ns_neut;
           for (int i=0; i<ns_charged; ++i) {
-	    struct gk_species *gks = &sbapp->species[i]; 
+            struct gk_species *gks = &sbapp->species[i]; 
             gk_species_bflux_accumulate(sbapp, &gks->bflux, bflux_out[li_charged+i], -1.0, bflux_in[li_charged+i]);
-	  }
+          }
           for (int i=0; i<ns_neut; ++i) {
-	    struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
+            struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
             gk_neut_species_bflux_accumulate(sbapp, &gkns->bflux, bflux_out_neut[li_neut+i], -1.0, bflux_in_neut[li_neut+i]);
           }
         }
@@ -159,21 +165,23 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
             struct gk_species *gks = &sbapp->species[i];
             fin[li_charged+i] = gks->f1;
             fout[li_charged+i] = gks->fnew;
+            fbar_in[li_charged+i] = gks->damping.fbar1;
+            fbar_out[li_charged+i] = gks->damping.fbarnew;
             // Boundary fluxes.
             bflux_in[li_charged+i] = gks->bflux.f1;
             bflux_out[li_charged+i] = gks->bflux.fnew;
           }
           for (int i=0; i<ns_neut; ++i) {
             struct gk_neut_species *gkns = &sbapp->neut_species[i];
-	    fin_neut[li_neut+i] = gkns->f1;
-	    fout_neut[li_neut+i] = gkns->fnew;
+            fin_neut[li_neut+i] = gkns->f1;
+            fout_neut[li_neut+i] = gkns->fnew;
             // Boundary fluxes.
             bflux_in_neut[li_neut+i] = gkns->bflux.f1;
             bflux_out_neut[li_neut+i] = gkns->bflux.fnew;
           }
         }
 
-        gyrokinetic_multib_forward_euler(app, tcurr+dt, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_multib_forward_euler(app, tcurr+dt, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
         if (st.dt_actual < dt) {
@@ -207,13 +215,15 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           for (int b=0; b<nblocks_local; ++b) {
             struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
             for (int i=0; i<ns_charged; ++i) {
-	      struct gk_species *gks = &sbapp->species[i]; 
-	      gk_species_combine(gks, gks->f1, 3.0/4.0, gks->f, 1.0/4.0, gks->fnew, &gks->local_ext);
+              struct gk_species *gks = &sbapp->species[i]; 
+              gk_species_combine(gks, gks->f1, 3.0/4.0, gks->f, 1.0/4.0, gks->fnew, &gks->local_ext);
+              gk_species_damping_combine(gks, gks->damping.fbar1, 3.0/4.0, gks->damping.fbar,
+                1.0/4.0, gks->damping.fbarnew, &gks->local_ext);
               gk_species_bflux_set(sbapp, &gks->bflux, gks->bflux.f1, 1.0/4.0, gks->bflux.fnew);
-	    }
+            }
             for (int i=0; i<ns_neut; ++i) {
-	      struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
-	      gk_neut_species_combine(gkns, gkns->f1, 3.0/4.0, gkns->f, 1.0/4.0, gkns->fnew, &gkns->local_ext);
+              struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
+              gk_neut_species_combine(gkns, gkns->f1, 3.0/4.0, gkns->f, 1.0/4.0, gkns->fnew, &gkns->local_ext);
               gk_neut_species_bflux_set(sbapp, &gkns->bflux, gkns->bflux.f1, 1.0/4.0, gkns->bflux.fnew);
             }
           }
@@ -245,24 +255,26 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           int li_charged = b * ns_charged;
           int li_neut = b * ns_neut;
           for (int i=0; i<ns_charged; ++i) {
-	    struct gk_species *gks = &sbapp->species[i]; 
+            struct gk_species *gks = &sbapp->species[i]; 
             fin[li_charged+i] = gks->f1;
             fout[li_charged+i] = gks->fnew;
+            fbar_in[li_charged+i] = gks->damping.fbar1;
+            fbar_out[li_charged+i] = gks->damping.fbarnew;
             // Boundary fluxes.
             bflux_in[li_charged+i] = gks->bflux.f1;
             bflux_out[li_charged+i] = gks->bflux.fnew;
           }
           for (int i=0; i<ns_neut; ++i) {
             struct gk_neut_species *gkns = &sbapp->neut_species[i];
-	    fin_neut[li_neut+i] = sbapp->neut_species[i].f1;
-	    fout_neut[li_neut+i] = sbapp->neut_species[i].fnew;
+            fin_neut[li_neut+i] = sbapp->neut_species[i].f1;
+            fout_neut[li_neut+i] = sbapp->neut_species[i].fnew;
             // Boundary fluxes.
             bflux_in_neut[li_neut+i] = gkns->bflux.f1;
             bflux_out_neut[li_neut+i] = gkns->bflux.fnew;
           }
         }
 
-        gyrokinetic_multib_forward_euler(app, tcurr+dt/2, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_multib_forward_euler(app, tcurr+dt/2, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
         if (st.dt_actual < dt) {
@@ -296,18 +308,22 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           for (int b=0; b<nblocks_local; ++b) {
             struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
             for (int i=0; i<ns_charged; ++i) {
-	      struct gk_species *gks = &sbapp->species[i]; 
+              struct gk_species *gks = &sbapp->species[i]; 
               // Step f.
-	      gk_species_combine(gks, gks->f1, 1.0/3.0, gks->f, 2.0/3.0, gks->fnew, &gks->local_ext);
-	      gk_species_copy_range(gks, gks->f, gks->f1, &gks->local_ext);
+              gk_species_combine(gks, gks->f1, 1.0/3.0, gks->f, 2.0/3.0, gks->fnew, &gks->local_ext);
+              gk_species_copy_range(gks, gks->f, gks->f1, &gks->local_ext);
+              gk_species_damping_combine(gks, gks->damping.fbar, 1.0/3.0, gks->damping.fbar1,
+                2.0/3.0, gks->damping.fbarnew, &gks->local_ext);
+              gk_species_damping_copy_range(gks, gks->damping.fbar1, gks->damping.fbar,
+                &gks->local_ext);
               // Step boundary fluxes.
               gk_species_bflux_set(sbapp, &gks->bflux, gks->bflux.f, 2.0/3.0, gks->bflux.fnew);
               gk_species_bflux_calc_voltime_integrated_mom(sbapp, gks, &gks->bflux, tcurr);
-	    }
+            }
             for (int i=0; i<ns_neut; ++i) {
-	      struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
-	      gk_neut_species_combine(gkns, gkns->f1, 1.0/3.0, gkns->f, 2.0/3.0, gkns->fnew, &gkns->local_ext);
-	      gk_neut_species_copy_range(gkns, gkns->f, gkns->f1, &gkns->local_ext);
+              struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
+              gk_neut_species_combine(gkns, gkns->f1, 1.0/3.0, gkns->f, 2.0/3.0, gkns->fnew, &gkns->local_ext);
+              gk_neut_species_copy_range(gkns, gkns->f, gkns->f1, &gkns->local_ext);
               // Step boundary fluxes.
               gk_neut_species_bflux_set(sbapp, &gkns->bflux, gkns->bflux.f, 2.0/3.0, gkns->bflux.fnew);
               gk_neut_species_bflux_calc_voltime_integrated_mom(sbapp, gkns, &gkns->bflux, tcurr);
@@ -350,7 +366,7 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
 
           for (int b=0; b<nblocks_local; ++b) {
             struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
-	    for (int i=0; i<ns_charged; ++i) {
+            for (int i=0; i<ns_charged; ++i) {
               struct gk_species *gks = &sbapp->species[i];
               // Compute moment of f_new to compute moment of df/dt.
               // Need to do it after the fields are updated.
diff --git a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
index 22ec1cf2bc..1fb640da1c 100644
--- a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
+++ b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
@@ -19,7 +19,8 @@ gyrokinetic_forward_euler(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   app->stat.nfeuler += 1;
 
   // Compute the time rate of change of the distributions, df/dt.
-  gyrokinetic_rhs(app, tcurr, dt, fin, fout, bflux_out, fin_neut, fout_neut, bflux_out_neut, st);
+  gyrokinetic_rhs(app, tcurr, dt, fin, fbar_in, fout, bflux_out, fin_neut, fout_neut,
+    bflux_out_neut, st);
 
   struct timespec wst = gkyl_wall_clock();
   // Complete update of distribution functions.

From d1a1d43967e2bfec0fd16d1e8a44638e0f02dd39 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 13:10:15 -0400
Subject: [PATCH 17/25] Depricate damping with the loss cone multiplier. This
 will cause merge conflicts with another PR and nobody is using this feature,
 as far as I know

---
 gyrokinetic/apps/gk_species_damping.c         | 156 ------------------
 gyrokinetic/apps/gk_species_fdot_multiplier.c |   4 +-
 gyrokinetic/apps/gkyl_gyrokinetic.h           |  12 +-
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h      |   5 -
 .../rt_gk_mirror_boltz_elc_damped_1x2v_p1.c   |   2 +-
 5 files changed, 6 insertions(+), 173 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 347566a6d5..82eae86f0d 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -1,8 +1,6 @@
 #include <assert.h>
 #include <gkyl_gyrokinetic_priv.h>
-#include <gkyl_loss_cone_mask_gyrokinetic.h>
 #include <gkyl_alloc.h>
-#include <gkyl_dg_basis_ops.h>
 
 // Damping state synchronization helpers.
 
@@ -212,32 +210,6 @@ gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk
   gkyl_array_accumulate(cflrate, 1.0, damp->rate);
 }
 
-void
-gk_species_damping_advance_loss_cone(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
-  struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
-  struct gkyl_array *f_buffer,
-  struct gkyl_array *rhs, struct gkyl_array *cflrate)
-{
-  // Find the potential at the mirror throat.
-  gkyl_dg_basis_ops_eval_array_at_coord_comp(phi, damp->bmag_max_coord,
-    app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
-  gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
-
-  // Project the loss cone mask.
-  gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
-    phi, damp->phi_m_global, damp->rate);
-
-  // Assemble the damping term -scale_prof * mask * f.
-  gkyl_array_set(f_buffer, 1.0, fin);
-  gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
-  gkyl_array_scale_by_cell(f_buffer, damp->rate);
-
-  // Add damping to f and the CFL frequency.
-  gkyl_array_accumulate(rhs, -1.0, f_buffer);
-  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
-}
-
 void
 gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
@@ -424,118 +396,6 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 
     damp->advance_func = gk_species_damping_advance_user_input;
   }
-  else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
-    damp->evolve = true; // Since the loss cone boundary is proportional to phi(t).
-
-    // Maximum bmag and its location.
-    // NOTE: if the same max bmag occurs at multiple locations,
-    // bmag_max_coord may have different values on different MPI processes.
-    double bmag_max_coord_ho[GKYL_MAX_CDIM];
-    double bmag_max_ho = gkyl_gk_geometry_reduce_arg_bmag(app->gk_geom, GKYL_MAX,
-      bmag_max_coord_ho);
-    double bmag_max_local = bmag_max_ho;
-    double bmag_max_global;
-    gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, &bmag_max_local,
-      &bmag_max_global);
-    double bmag_max_coord_local[app->cdim], bmag_max_coord_global[app->cdim];
-    if (fabs(bmag_max_ho - bmag_max_global) < 1e-16) {
-      for (int d = 0; d < app->cdim; d++) {
-        bmag_max_coord_local[d] = bmag_max_coord_ho[d];
-      }
-    }
-    else {
-      for (int d = 0; d < app->cdim; d++) {
-        bmag_max_coord_local[d] = -DBL_MAX;
-      }
-    }
-    gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, app->cdim, bmag_max_coord_local,
-      bmag_max_coord_global);
-
-    if (app->use_gpu) {
-      damp->bmag_max = gkyl_cu_malloc(sizeof(double));
-      damp->bmag_max_coord = gkyl_cu_malloc(app->cdim * sizeof(double));
-      gkyl_cu_memcpy(damp->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
-      gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double),
-        GKYL_CU_MEMCPY_H2D);
-    }
-    else {
-      damp->bmag_max = gkyl_malloc(sizeof(double));
-      damp->bmag_max_coord = gkyl_malloc(app->cdim * sizeof(double));
-      memcpy(damp->bmag_max, &bmag_max_global, sizeof(double));
-      memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim * sizeof(double));
-    }
-
-    // Electrostatic potential at bmag_max_coord.
-    if (app->use_gpu) {
-      damp->phi_m = gkyl_cu_malloc(sizeof(double));
-      damp->phi_m_global = gkyl_cu_malloc(sizeof(double));
-    }
-    else {
-      damp->phi_m = gkyl_malloc(sizeof(double));
-      damp->phi_m_global = gkyl_malloc(sizeof(double));
-    }
-
-    // Operator that projects the loss cone mask.
-    struct gkyl_loss_cone_mask_gyrokinetic_inp inp_proj = {
-      .phase_grid = &gks->grid,
-      .conf_basis = &app->basis,
-      .phase_basis = &gks->basis,
-      .conf_range = &app->local,
-      .conf_range_ext = &app->local_ext,
-      .vel_range = &gks->local_vel,
-      .vel_map = gks->vel_map,
-      .bmag = app->gk_geom->geo_int.bmag,
-      .bmag_max = damp->bmag_max,
-      .bmag_max_loc = damp->bmag_max_coord,
-      .mass = gks->info.mass,
-      .charge = gks->info.charge,
-      .num_quad = num_quad,
-      .use_gpu = app->use_gpu,
-    };
-    damp->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew(&inp_proj);
-
-    // Project the conf-space rate profile provided.
-    struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, 1,
-      gks->local_ext.volume);
-    struct gkyl_array *scale_prof_high_order_ho = app->use_gpu? mkarr(false,
-      scale_prof_high_order->ncomp, scale_prof_high_order->size)
-                                                    : gkyl_array_acquire(scale_prof_high_order);
-
-    if (gks->info.damping.rate_profile) {
-      struct gkyl_basis basis_mult;
-      gkyl_cart_modal_serendip(&basis_mult, gks->basis.ndim, 0);
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &basis_mult, num_quad, 1,
-        gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
-      gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
-      gkyl_proj_on_basis_release(projup);
-      gkyl_array_copy(scale_prof_high_order, scale_prof_high_order_ho);
-    }
-    else {
-      gkyl_array_clear(scale_prof_high_order, 1.0);
-    }
-
-    gkyl_array_scale(scale_prof_high_order, rate_const == 0.0 ? 1.0 : rate_const);
-
-    damp->scale_prof = mkarr(app->use_gpu, 1, gks->local_ext.volume);
-    gkyl_array_set_offset(damp->scale_prof, pow(sqrt(2.0), gks->grid.ndim), scale_prof_high_order,
-      0);
-
-    gkyl_array_release(scale_prof_high_order_ho);
-    gkyl_array_release(scale_prof_high_order);
-
-    // Compute the initial damping rate (assuming phi=0 because phi hasn't been computed).
-    // Find the potential at the mirror throat.
-    gkyl_dg_basis_ops_eval_array_at_coord_comp(app->field->phi_smooth, damp->bmag_max_coord,
-      app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
-    gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
-    // Project the loss cone mask.
-    gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
-      app->field->phi_smooth, damp->phi_m_global, damp->rate);
-    // Multiply by the user's scaling profile.
-    gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
-
-    damp->advance_func = gk_species_damping_advance_loss_cone;
-  }
   else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
     // LOW_PASS_FILTER supports either a projected phase-space profile or a uniform constant.
     damp->evolve = true; // Since fbar must evolve in time.
@@ -592,22 +452,6 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
     if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
       // Nothing to release.
     }
-    else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
-      if (app->use_gpu) {
-        gkyl_cu_free(damp->bmag_max);
-        gkyl_cu_free(damp->bmag_max_coord);
-        gkyl_cu_free(damp->phi_m);
-        gkyl_cu_free(damp->phi_m_global);
-      }
-      else {
-        gkyl_free(damp->bmag_max);
-        gkyl_free(damp->bmag_max_coord);
-        gkyl_free(damp->phi_m);
-        gkyl_free(damp->phi_m_global);
-      }
-      gkyl_loss_cone_mask_gyrokinetic_release(damp->lcm_proj_op);
-      gkyl_array_release(damp->scale_prof);
-    }
     else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
       gkyl_array_release(damp->fbar);
       gkyl_array_release(damp->fbar1);
diff --git a/gyrokinetic/apps/gk_species_fdot_multiplier.c b/gyrokinetic/apps/gk_species_fdot_multiplier.c
index 0a983aeaa6..78a697e669 100644
--- a/gyrokinetic/apps/gk_species_fdot_multiplier.c
+++ b/gyrokinetic/apps/gk_species_fdot_multiplier.c
@@ -277,10 +277,10 @@ gk_species_fdot_multiplier_release(const struct gkyl_gyrokinetic_app *app, const
     if (fdmul->write_diagnostics)
       gkyl_array_release(fdmul->multiplier_host);
 
-    if (fdmul->type == GKYL_GK_DAMPING_USER_INPUT) {
+    if (fdmul->type == GKYL_GK_FDOT_MULTIPLIER_USER_INPUT) {
       // Nothing to release.
     }
-    else if (fdmul->type == GKYL_GK_DAMPING_LOSS_CONE) {
+    else if (fdmul->type == GKYL_GK_FDOT_MULTIPLIER_LOSS_CONE) {
       if (app->use_gpu) {
         gkyl_cu_free(fdmul->bmag_max);
         gkyl_cu_free(fdmul->bmag_max_coord);
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 60c6b26f22..0e4da30825 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -317,7 +317,6 @@ struct gkyl_gyrokinetic_positivity {
 enum gkyl_gyrokinetic_damping_type {
   GKYL_GK_DAMPING_NONE = 0,
   GKYL_GK_DAMPING_USER_INPUT,
-  GKYL_GK_DAMPING_LOSS_CONE,
   GKYL_GK_DAMPING_LOW_PASS_FILTER,
 };
 
@@ -325,15 +324,10 @@ struct gkyl_gyrokinetic_damping {
   // Add a damping term to the RHS of the gyrokinetic equation
   //   df/dt = - rate(z) * f
   // with the function rate(z) being:
-  //   - a function given by the user (type = GKYL_PROPORTIONAL_TERM_USER_INPUT).
-  //   - I_loss(z) * scale_factor * scale_profile(z), where I_loss(z) is =1 in the loss
-  //     cone and 0 in the confined region (type = GKYL_PROPORTIONAL_TERM_LOSS_CONE).
+  //   - a function given by the user (type = GKYL_GK_DAMPING_USER_INPUT).
+  //   Alternativly, a low-pass filter can be applied (type = GKYL_GK_DAMPING_LOW_PASS_FILTER), which projects f to a lower-dimensional function fbar(z) and applies the damping as df/dt = - rate(z) * (f - fbar).
   enum gkyl_gyrokinetic_damping_type type;
-  // Optional constant damping factor.
-  // Semantics:
-  //  - if rate_profile is NULL: use rate_const as the full rate profile.
-  //  - if rate_profile is provided: effective rate is rate_const*rate_profile,
-  //    where rate_const defaults to 1.0 when left as 0.0.
+
   double rate_const;
   void (*rate_profile)(double t, const double *xn, double *fout, void *ctx);
   void *rate_profile_ctx; // Context for rate_profile function.
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 2ec502bfe6..d24e6d8324 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -805,11 +805,6 @@ struct gk_damping {
   bool write_fbar; // Whether to write low-pass-filter fbar diagnostics.
   struct gkyl_array *rate; // Damping rate.
   struct gkyl_array *rate_host; // Host copy for use in IO and projecting.
-  struct gkyl_loss_cone_mask_gyrokinetic *lcm_proj_op; // Operator that projects the loss cone mask.
-  double *bmag_max; // Maximum magnetic field amplitude.
-  double *bmag_max_coord; // Location of bmag_max.
-  double *phi_m, *phi_m_global; // Electrostatic potential at bmag_max.
-  struct gkyl_array *scale_prof; // Conf-space scaling factor profile.
   // Low-pass filter variables
   struct gkyl_array *fbar; // Filtered/averaged distribution function.
   struct gkyl_array *fbar1, *fbarnew; // SSPRK3 stage arrays for filtered distribution.
diff --git a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
index a50ff1162d..91128d814e 100644
--- a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
@@ -704,7 +704,7 @@ int main(int argc, char **argv)
     },
 
     .damping = {
-      .type = GKYL_GK_DAMPING_LOSS_CONE,
+      .type = GKYL_GK_DAMPING_USER_INPUT,
       .rate_profile = loss_cone_damping_rate_scaling,
       .rate_profile_ctx = &ctx,
       .write_rate = true,

From bd9a6170d25a2ab98e1ac01a06e12ae0025c6c78 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 15:31:08 -0400
Subject: [PATCH 18/25] Bring back whether fbar is gkhyb basis or p0. I get a
 feeling that having fbar having gkhyb basis will help with rapid oscillations
 in gradients of f, driving negativity

---
 gyrokinetic/apps/gk_species_damping.c    | 199 ++++++++++++++++-------
 gyrokinetic/apps/gkyl_gyrokinetic.h      |   1 +
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h |   1 +
 gyrokinetic/creg/rt_gk_wham_1x2v_p1.c    |   3 +-
 4 files changed, 147 insertions(+), 57 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 82eae86f0d..9e8e7dbf9b 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -4,8 +4,9 @@
 
 // Damping state synchronization helpers.
 
-void
-gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk_damping *damp,
+static void
+gk_species_damping_set_fbar_to_f_cellwise_const(const struct gk_species *gks,
+  struct gk_damping *damp,
   const struct gkyl_array *f)
 {
   gkyl_array_set_offset(damp->fbar, 1.0, f, 0);
@@ -13,6 +14,15 @@ gk_species_damping_set_fbar_to_f_enabled(const struct gk_species *gks, struct gk
   gkyl_array_set_offset(damp->fbarnew, 1.0, f, 0);
 }
 
+static void
+gk_species_damping_set_fbar_to_f_same_basis(const struct gk_species *gks, struct gk_damping *damp,
+  const struct gkyl_array *f)
+{
+  gkyl_array_set(damp->fbar, 1.0, f);
+  gkyl_array_set(damp->fbar1, 1.0, f);
+  gkyl_array_set(damp->fbarnew, 1.0, f);
+}
+
 void
 gk_species_damping_set_fbar_to_f_disabled(const struct gk_species *gks, struct gk_damping *damp,
   const struct gkyl_array *f)
@@ -83,7 +93,7 @@ gk_species_damping_write_fbar_disabled(gkyl_gyrokinetic_app *app, struct gk_spec
 }
 
 static void
-gk_species_damping_write_fbar_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+gk_species_damping_write_fbar_cellwise_const(gkyl_gyrokinetic_app *app, struct gk_species *gks,
   double tm, int frame)
 {
   // Metadata from app/species/geometry. fbar is always stored as a p0 field.
@@ -123,6 +133,43 @@ gk_species_damping_write_fbar_enabled(gkyl_gyrokinetic_app *app, struct gk_speci
   gkyl_msgpack_data_release(mt_fbar);
 }
 
+static void
+gk_species_damping_write_fbar_same_basis(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+  // Metadata from app/species/geometry using the species phase-space basis.
+  gkyl_msgpack_map_elem_set_double(app->io_meta_basic_len, app->io_meta_basic, "time", tm);
+  gkyl_msgpack_map_elem_set_uint(app->io_meta_basic_len, app->io_meta_basic, "frame", frame);
+
+  int io_meta_fbar_len[] = {
+    app->io_meta_basic_len,
+    gks->io_meta_len,
+    app->gk_geom->io_meta_len
+  };
+  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = {
+    app->io_meta_basic,
+    gks->io_meta,
+    app->gk_geom->io_meta
+  };
+  struct gkyl_msgpack_data *mt_fbar = gkyl_msgpack_create_union(
+    sizeof(io_meta_fbar_len) / sizeof(int), io_meta_fbar_len, io_meta_fbar);
+
+  const char *fmt_fbar = "%s-%s_damping_fbar_%d.gkyl";
+  int sz_fbar = gkyl_calc_strlen(fmt_fbar, app->name, gks->info.name, frame);
+  char fileNm_fbar[sz_fbar + 1]; // ensures no buffer overflow
+  snprintf(fileNm_fbar, sizeof fileNm_fbar, fmt_fbar, app->name, gks->info.name, frame);
+
+  // Copy fbar from device to host before writing it out.
+  if (app->use_gpu)
+    gkyl_array_copy(gks->damping.fbar_host, gks->damping.fbar);
+
+  gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt_fbar, gks->damping.fbar_host,
+    fileNm_fbar);
+  app->stat.n_io += 1;
+
+  gkyl_msgpack_data_release(mt_fbar);
+}
+
 void
 gk_species_damping_write_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
   int frame)
@@ -210,8 +257,9 @@ gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk
   gkyl_array_accumulate(cflrate, 1.0, damp->rate);
 }
 
-void
-gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+static void
+gk_species_damping_advance_low_pass_filter_cellwise_const(gkyl_gyrokinetic_app *app,
+  const struct gk_species *gks,
   struct gk_damping *damp,
   const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
   struct gkyl_array *f_buffer,
@@ -228,6 +276,25 @@ gk_species_damping_advance_low_pass_filter(gkyl_gyrokinetic_app *app, const stru
   gkyl_array_accumulate(cflrate, 1.0, damp->rate);
 }
 
+static void
+gk_species_damping_advance_low_pass_filter_same_basis(gkyl_gyrokinetic_app *app,
+  const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  // Compute f - fbar, where fbar uses the same basis as f.
+  gkyl_array_set(f_buffer, 1.0, fin); // f_buffer = fin
+  gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = fin - fbar
+  gkyl_array_scale_by_cell(f_buffer, damp->rate); // f_buffer = rate * (fin - fbar)
+
+  // Add damping term to RHS: df/dt -= rate * (f - fbar)
+  // Add to the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
+
 void
 gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_damping *damp,
@@ -256,8 +323,8 @@ gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
 {
 }
 
-void
-gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
+static void
+gk_species_damping_calc_fbar_rhs_cellwise_const(const struct gk_damping *damp,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
 {
   // rhs_fbar = rate * (f - fbar)
@@ -266,6 +333,16 @@ gk_species_damping_calc_fbar_rhs_enabled(const struct gk_damping *damp,
   gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
 }
 
+static void
+gk_species_damping_calc_fbar_rhs_same_basis(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  // rhs_fbar = rate * (f - fbar)
+  gkyl_array_set(rhs_fbar, 1.0, fin);
+  gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
+  gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
+}
+
 void
 gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
   const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
@@ -354,6 +431,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->evolve = false; // Whether the rate is time dependent.
   damp->write_rate = gks->info.damping.write_rate;
   damp->write_fbar = gks->info.damping.write_fbar;
+  damp->cellwise_const = gks->info.damping.cellwise_const;
   const double rate_const = gks->info.damping.rate_const;
 
   int num_quad = gks->info.damping.num_quad? gks->info.damping.num_quad : 1; // Default is a p=0 mask.
@@ -370,70 +448,79 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->combine_func = gk_species_damping_combine_disabled;
   damp->copy_func = gk_species_damping_copy_range_disabled;
 
+  // Default pointers for fbar arrays
+  damp->fbar = 0;
+  damp->fbar1 = 0;
+  damp->fbarnew = 0;
+  damp->fbar_host = 0;
+
   if (!damp->type) {
     return;
   }
 
-  if (damp->write_rate)
+  if (damp->write_rate) {
     damp->write_rate_func = gk_species_damping_write_rate_enabled;
+  }
 
   // Allocate rate array.
   damp->rate = mkarr(app->use_gpu, 1, gks->local_ext.volume);
   damp->rate_host = damp->rate;
-  if (app->use_gpu)
+  if (app->use_gpu) {
     damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size);
+  }
+
+  if (gks->info.damping.rate_profile) {
+    gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
+      gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
+    gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
+  }
+  else {
+    gkyl_array_clear(damp->rate, rate_const);
+  }
 
   if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
     // USER_INPUT supports either a projected profile, or a uniform constant if no profile is provided.
-    if (gks->info.damping.rate_profile) {
-      gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
-        gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
-      gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
-    }
-    else {
-      gkyl_array_clear(damp->rate, rate_const);
-    }
-
     damp->advance_func = gk_species_damping_advance_user_input;
   }
   else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
     // LOW_PASS_FILTER supports either a projected phase-space profile or a uniform constant.
-    damp->evolve = true; // Since fbar must evolve in time.
-    if (gks->info.damping.rate_profile) {
-      gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
-        gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
-      gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
-    }
-    else {
-      gkyl_array_clear(damp->rate, rate_const);
-    }
+    damp->evolve = true;
+    
+    // Allocate filtered distribution function array in the requested basis.
+    int fbar_ncomp = damp->cellwise_const ? 1 : gks->basis.num_basis;
+    damp->fbar = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
+    damp->fbar1 = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
+    damp->fbarnew = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
+    damp->fbar_host = damp->fbar;
+    gkyl_array_clear(damp->fbar, 0.0);
 
-    // Allocate filtered distribution function array as p0.
-    damp->fbar = mkarr(app->use_gpu, 1, gks->local_ext.volume);
-    damp->fbar1 = mkarr(app->use_gpu, 1, gks->local_ext.volume);
-    damp->fbarnew = mkarr(app->use_gpu, 1, gks->local_ext.volume);
-    damp->fbar_host = 0;
     if (damp->write_fbar) {
-      damp->write_fbar_func = gk_species_damping_write_fbar_enabled;
-      damp->fbar_host = damp->fbar;
+      if (damp->cellwise_const) {
+        damp->write_fbar_func = gk_species_damping_write_fbar_cellwise_const;
+      }
+      else {
+        damp->write_fbar_func = gk_species_damping_write_fbar_same_basis;
+      }
       if (app->use_gpu)
         damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
     }
 
-    // Initialize fbar from the projection of the initial distribution
-    // (will be set from the initial f in the main app loop)
-    gkyl_array_clear(damp->fbar, 0.0);
-
-    damp->advance_func = gk_species_damping_advance_low_pass_filter;
-    damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_enabled;
-    damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_enabled;
+    if (damp->cellwise_const) {
+      damp->advance_func = gk_species_damping_advance_low_pass_filter_cellwise_const;
+      damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_cellwise_const;
+      damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_cellwise_const;
+    }
+    else {
+      damp->advance_func = gk_species_damping_advance_low_pass_filter_same_basis;
+      damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_same_basis;
+      damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_same_basis;
+    }
     damp->forward_euler_func = gk_species_damping_forward_euler_enabled;
     damp->combine_func = gk_species_damping_combine_enabled;
     damp->copy_func = gk_species_damping_copy_range_enabled;
   }
 
-  // Set function pointers chosen at runtime.
-  if (damp->evolve) {
+  if (damp->evolve && damp->write_rate) {
     damp->write_func = gk_species_damping_write_enabled;
   }
   else {
@@ -444,20 +531,20 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 void
 gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp)
 {
-  if (damp->type) {
-    gkyl_array_release(damp->rate);
-    if (app->use_gpu)
-      gkyl_array_release(damp->rate_host);
+  if (!damp->type) {
+    return;
+  }
+  gkyl_array_release(damp->rate);
+  if (app->use_gpu) {
+    gkyl_array_release(damp->rate_host);
+  }
 
-    if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
-      // Nothing to release.
-    }
-    else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-      gkyl_array_release(damp->fbar);
-      gkyl_array_release(damp->fbar1);
-      gkyl_array_release(damp->fbarnew);
-      if (app->use_gpu && damp->write_fbar)
-        gkyl_array_release(damp->fbar_host);
+  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+    gkyl_array_release(damp->fbar);
+    gkyl_array_release(damp->fbar1);
+    gkyl_array_release(damp->fbarnew);
+    if (app->use_gpu && damp->write_fbar) {
+      gkyl_array_release(damp->fbar_host);
     }
   }
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 0e4da30825..388501ae5e 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -334,6 +334,7 @@ struct gkyl_gyrokinetic_damping {
   int num_quad; // Number of quadrature points in each direction to use in projecting the rate.
   bool write_rate; // Whether to write damping-rate diagnostics.
   bool write_fbar; // For low-pass filter damping, write fbar diagnostics each output frame.
+  bool cellwise_const; // For low-pass filter damping, whether fbar uses a single value per cell.
 };
 
 enum gkyl_gyrokinetic_fdot_multiplier_type {
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index d24e6d8324..edb7265fec 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -803,6 +803,7 @@ struct gk_damping {
   bool evolve; // Whether the source is time dependent.
   bool write_rate; // Whether to write damping-rate diagnostics.
   bool write_fbar; // Whether to write low-pass-filter fbar diagnostics.
+  bool cellwise_const; // Whether low-pass-filter fbar is stored as one value per cell.
   struct gkyl_array *rate; // Damping rate.
   struct gkyl_array *rate_host; // Host copy for use in IO and projecting.
   // Low-pass filter variables
diff --git a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
index 664643e299..f7f489ff17 100644
--- a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
@@ -613,8 +613,9 @@ int main(int argc, char **argv)
       .rate_const = 1e-6,
       .write_rate = true,
       .write_fbar = true,
+      .cellwise_const = true,
     },
-
+    
     .source = {
       .source_id = GKYL_PROJ_SOURCE,
       .num_sources = 1,

From 0b064240ffb874bf4dbe95d6f6e893dc2de36c11 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 15:45:48 -0400
Subject: [PATCH 19/25] Add the appropriate reset method to gyrokinetic.c

---
 gyrokinetic/apps/gyrokinetic.c | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c
index 6e3bd13f58..67f5ec886f 100644
--- a/gyrokinetic/apps/gyrokinetic.c
+++ b/gyrokinetic/apps/gyrokinetic.c
@@ -3072,6 +3072,14 @@ gkyl_gyrokinetic_app_reset_species_positivity(gkyl_gyrokinetic_app* app, double
   gk_species_positivity_reset(app, tm, gks, &gks->positivity, pos_inp);
 }
 
+void
+gkyl_gyrokinetic_app_reset_species_damping(gkyl_gyrokinetic_app* app, double tm,
+  const char *species_name, struct gkyl_gyrokinetic_damping damping_inp)
+{
+  struct gk_species *gks = gk_find_species(app, species_name);
+  gk_species_damping_reset(app, tm, gks, &gks->damping, damping_inp);
+}
+
 void
 gkyl_gyrokinetic_app_reset_field(gkyl_gyrokinetic_app* app, double tm,
   struct gkyl_gyrokinetic_field field_inp)

From 26b3917649b599bef1c2cff9366469c225eeda8a Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 16:04:17 -0400
Subject: [PATCH 20/25] Make sure the timings are where they should be. Update
 the header to the function call for reset

---
 gyrokinetic/apps/gk_species_damping.c |  1 -
 gyrokinetic/apps/gkyl_gyrokinetic.h   | 12 ++++++++++++
 2 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 9e8e7dbf9b..e5123fb1e7 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -176,7 +176,6 @@ gk_species_damping_write_enabled(gkyl_gyrokinetic_app *app, struct gk_species *g
 {
   struct timespec wst = gkyl_wall_clock();
   gks->damping.write_rate_func(app, gks, tm, frame);
-
   gks->damping.write_fbar_func(app, gks, tm, frame);
   app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wst);
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 388501ae5e..5d6ccf4ec7 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -1347,6 +1347,18 @@ void gkyl_gyrokinetic_app_reset_species_collisionless(gkyl_gyrokinetic_app* app,
 void gkyl_gyrokinetic_app_reset_species_positivity(gkyl_gyrokinetic_app* app, double tm,
   const char *species_name, struct gkyl_gyrokinetic_positivity pos_inp);
 
+/**
+ * Reset the damping for a given species.
+ *
+ * @param app App object.
+ * @param tm Time-stamp.
+ * @param species_name Name of the species to reset.
+ * @param damping_inp Input parameters for damping terms.
+ */
+void
+gkyl_gyrokinetic_app_reset_species_damping(gkyl_gyrokinetic_app* app, double tm,
+  const char *species_name, struct gkyl_gyrokinetic_damping damping_inp);
+
 /**
  * Reset the field solver.
  *

From 6e876a2eed389f758fbb3e8e93acb4b4f309146d Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 20:40:12 -0400
Subject: [PATCH 21/25] Change the flag for wite init only dependent on evolve.
 This was a mistake on my part

---
 gyrokinetic/apps/gk_species_damping.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index e5123fb1e7..a965220207 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -519,7 +519,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
     damp->copy_func = gk_species_damping_copy_range_enabled;
   }
 
-  if (damp->evolve && damp->write_rate) {
+  if (damp->evolve) {
     damp->write_func = gk_species_damping_write_enabled;
   }
   else {

From 44bd9c21f927f072dd6cf7cbdb588e99655db2ea Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 22:52:13 -0400
Subject: [PATCH 22/25] Implement low-pass filter damping enhancements:
 preserve fbar state across resets and add configuration options

This was pretty difficult to figure out and make sure it's valgrind clean. The damping_release method needed to be circumvented, posing memory sanitization issues. This is valgrind clean with and without the flag. rt_gk_sheath_1x2v is also valgrind clean
---
 gyrokinetic/apps/gk_species.c                 |  4 ++
 gyrokinetic/apps/gk_species_damping.c         | 56 ++++++++++++++++++-
 gyrokinetic/apps/gkyl_gyrokinetic.h           |  1 +
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h      |  3 +
 .../creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c | 19 ++++++-
 5 files changed, 78 insertions(+), 5 deletions(-)

diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c
index 0bf3791865..4981d55614 100644
--- a/gyrokinetic/apps/gk_species.c
+++ b/gyrokinetic/apps/gk_species.c
@@ -1505,6 +1505,7 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st
   gk_species_anomalous_diff_init(app, gks, &gks->anom_diff);
 
   // Damping term -nu*f on RHS.
+  gks->damping = (struct gk_damping) { .fbar_backup_initialized = false };
   gk_species_damping_init(app, gks, &gks->damping);
 
   // Function multiplying df/dt.
@@ -1914,6 +1915,9 @@ gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *gks
   gk_species_source_release(app, &gks->src);
 
   gk_species_damping_release(app, &gks->damping);
+  if (gks->damping.fbar_backup_initialized) {
+    gkyl_array_release(gks->damping.fbar_backup);
+  }
 
   gk_species_fdot_multiplier_release(app, &gks->fdot_mult);
 
diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index a965220207..d6ccee72f1 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -431,6 +431,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->write_rate = gks->info.damping.write_rate;
   damp->write_fbar = gks->info.damping.write_fbar;
   damp->cellwise_const = gks->info.damping.cellwise_const;
+  damp->do_not_reset_fbar = gks->info.damping.do_not_reset_fbar;
   const double rate_const = gks->info.damping.rate_const;
 
   int num_quad = gks->info.damping.num_quad? gks->info.damping.num_quad : 1; // Default is a p=0 mask.
@@ -452,6 +453,9 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->fbar1 = 0;
   damp->fbarnew = 0;
   damp->fbar_host = 0;
+  if (!damp->fbar_backup_initialized) {
+    damp->fbar_backup = 0;
+  }
 
   if (!damp->type) {
     return;
@@ -500,8 +504,9 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       else {
         damp->write_fbar_func = gk_species_damping_write_fbar_same_basis;
       }
-      if (app->use_gpu)
+      if (app->use_gpu) {
         damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
+      }
     }
 
     if (damp->cellwise_const) {
@@ -549,14 +554,61 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
 }
 
 // Damping runtime reset.
+static bool
+gk_species_damping_can_restore_fbar(const struct gk_species *gks,
+  const struct gkyl_gyrokinetic_damping *damp_inp, const struct gkyl_array *fbar_backup)
+{
+  // Determines whether the provided backup of fbar state is compatible with the requested damping configuration, and can be used to preserve fbar state across a reset.
+  if (damp_inp->type != GKYL_GK_DAMPING_LOW_PASS_FILTER || !fbar_backup
+    || !damp_inp->do_not_reset_fbar) {
+    return false;
+  }
+
+  int fbar_ncomp = damp_inp->cellwise_const ? 1 : gks->basis.num_basis;
+  return fbar_backup->ncomp == fbar_ncomp && fbar_backup->size == gks->local_ext.volume;
+}
 
 void
 gk_species_damping_reset(gkyl_gyrokinetic_app *app, double tm, struct gk_species *gks,
   struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp)
 {
+  bool preserve_fbar = false;
+  struct gkyl_array *fbar_backup = 0;
+  bool keep_lpf_state = damp_inp.type == GKYL_GK_DAMPING_LOW_PASS_FILTER
+    && damp_inp.do_not_reset_fbar;
+
+  if (keep_lpf_state) {
+    // Prefer latest LPF state if currently active, otherwise keep parked backup.
+    if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER && damp->fbar) {
+      fbar_backup = gkyl_array_acquire(damp->fbar);
+    }
+    else if (damp->fbar_backup_initialized) {
+      fbar_backup = damp->fbar_backup;
+      damp->fbar_backup = 0;
+      damp->fbar_backup_initialized = false;
+    }
+    preserve_fbar = gk_species_damping_can_restore_fbar(gks, &damp_inp, fbar_backup);
+  }
+  else if (damp->fbar_backup_initialized) {
+    // User disabled preservation: drop any parked state.
+    gkyl_array_release(damp->fbar_backup);
+    damp->fbar_backup = 0;
+    damp->fbar_backup_initialized = false;
+  }
+
   gk_species_damping_release(app, damp);
 
+  damp->fbar_backup = fbar_backup;
+  damp->fbar_backup_initialized = fbar_backup != 0;
   gks->info.damping = damp_inp;
   gk_species_damping_init(app, gks, damp);
-  gk_species_damping_set_fbar_to_f(gks, damp, gks->f);
+
+  gk_species_damping_set_fbar_to_f(gks, damp, preserve_fbar ? damp->fbar_backup : gks->f);
+
+  // If backup has been consumed to initialize LPF state, drop parked handle.
+  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER && damp->fbar_backup_initialized) {
+    gkyl_array_release(damp->fbar_backup);
+    damp->fbar_backup = 0;
+    damp->fbar_backup_initialized = false;
+  }
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 5d6ccf4ec7..14ab2b33fd 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -335,6 +335,7 @@ struct gkyl_gyrokinetic_damping {
   bool write_rate; // Whether to write damping-rate diagnostics.
   bool write_fbar; // For low-pass filter damping, write fbar diagnostics each output frame.
   bool cellwise_const; // For low-pass filter damping, whether fbar uses a single value per cell.
+  bool do_not_reset_fbar; // On runtime reset, preserve the existing low-pass-filter fbar instead of reinitializing from f.
 };
 
 enum gkyl_gyrokinetic_fdot_multiplier_type {
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index edb7265fec..413547129a 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -804,12 +804,15 @@ struct gk_damping {
   bool write_rate; // Whether to write damping-rate diagnostics.
   bool write_fbar; // Whether to write low-pass-filter fbar diagnostics.
   bool cellwise_const; // Whether low-pass-filter fbar is stored as one value per cell.
+  bool do_not_reset_fbar; // Whether runtime reset should preserve existing fbar.
   struct gkyl_array *rate; // Damping rate.
   struct gkyl_array *rate_host; // Host copy for use in IO and projecting.
   // Low-pass filter variables
   struct gkyl_array *fbar; // Filtered/averaged distribution function.
   struct gkyl_array *fbar1, *fbarnew; // SSPRK3 stage arrays for filtered distribution.
   struct gkyl_array *fbar_host; // Host copy of fbar for use in IO.
+  struct gkyl_array *fbar_backup; // Parked backup retained across resets while LPF is inactive.
+  bool fbar_backup_initialized; // Whether fbar_backup holds a valid parked handle.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*write_rate_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
diff --git a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
index 87c13efe32..8a18f1a128 100644
--- a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
@@ -27,6 +27,7 @@ struct gk_poa_phase_params {
   bool is_static_field; // Whether to evolve the field.
   bool is_positivity_enabled; // Whether positivity is enabled.
   enum gkyl_gyrokinetic_fdot_multiplier_type fdot_mult_type; // Type of df/dt multipler.
+  enum gkyl_gyrokinetic_damping_type damping_type; // Type of damping.
 };
 
 // Define the context of the simulation. This is basically all the globals
@@ -576,6 +577,7 @@ create_ctx(void)
     poa_phases[2*i].is_static_field = is_static_field_oap;
     poa_phases[2*i].fdot_mult_type = fdot_mult_type_oap;
     poa_phases[2*i].is_positivity_enabled = is_positivity_enabled_oap;
+    poa_phases[2*i].damping_type = GKYL_GK_DAMPING_LOW_PASS_FILTER;
 
     // FDPs.
     poa_phases[2*i+1].phase = GK_POA_FDP;
@@ -585,6 +587,7 @@ create_ctx(void)
     poa_phases[2*i+1].is_static_field = is_static_field_fdp;
     poa_phases[2*i+1].fdot_mult_type = fdot_mult_type_fdp;
     poa_phases[2*i+1].is_positivity_enabled = is_positivity_enabled_fdp;
+    poa_phases[2*i+1].damping_type = GKYL_GK_DAMPING_NONE;
   }
   // Add an extra, longer FDP.
   poa_phases[num_phases-1].phase = GK_POA_FDP;
@@ -594,6 +597,7 @@ create_ctx(void)
   poa_phases[num_phases-1].is_static_field = is_static_field_fdp;
   poa_phases[num_phases-1].fdot_mult_type = fdot_mult_type_fdp;
   poa_phases[num_phases-1].is_positivity_enabled = is_positivity_enabled_fdp;
+  poa_phases[num_phases-1].damping_type = GKYL_GK_DAMPING_LOW_PASS_FILTER;
 
   double write_phase_freq = 0.5; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
   double int_diag_calc_freq = 5; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
@@ -777,10 +781,19 @@ void run_phase(gkyl_gyrokinetic_app* app, struct gk_mirror_ctx *ctx, double num_
     .type = pparams->is_positivity_enabled? GKYL_GK_POSITIVITY_SHIFT : GKYL_GK_POSITIVITY_NONE,
     .write_diagnostics = pparams->is_positivity_enabled,
   };
+  struct gkyl_gyrokinetic_damping damping_inp = {
+    .type = pparams->damping_type,
+    .cellwise_const = false,
+    .rate_const = 1/5e-6,
+    .write_fbar = true,
+    .write_rate = true,
+    .do_not_reset_fbar = true,
+  };
 
   gkyl_gyrokinetic_app_reset_species_collisionless(app, t_curr, "ion", collisionless_inp);
   gkyl_gyrokinetic_app_reset_species_fdot_multiplier(app, t_curr, "ion", fdot_mult_inp);
   gkyl_gyrokinetic_app_reset_species_positivity(app, t_curr, "ion", positivity_inp);
+  gkyl_gyrokinetic_app_reset_species_damping(app, t_curr, "ion", damping_inp);
   gkyl_gyrokinetic_app_reset_field(app, t_curr, field_inp);
 
   // Compute initial guess of maximum stable time-step.
@@ -908,7 +921,7 @@ int main(int argc, char **argv)
       .num_sources = 1,
       .projection[0] = {
         .proj_id = GKYL_PROJ_MAXWELLIAN_PRIM, 
-	.density = eval_density_ion_source,
+        .density = eval_density_ion_source,
         .upar = eval_upar_ion_source,
         .temp = eval_temp_ion_source,
         .ctx_density = &ctx,
@@ -1038,8 +1051,8 @@ int main(int argc, char **argv)
     write_data(&trig_write_conf, &trig_write_phase, app, tfs.t_curr, true);
   }
 
-  if (app_args.num_steps != INT_MAX)
-    phase_idx_end = 1;
+  // if (app_args.num_steps != INT_MAX)
+  //   phase_idx_end = 1;
 
   // Loop over number of number of phases;
   for (int pit=phase_idx_init; pit<phase_idx_end; pit++) {

From 6d13d26217b1f6ea5ea5f6c65bf71ceb33414bcf Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sat, 4 Apr 2026 22:52:27 -0400
Subject: [PATCH 23/25] Push changes to the regression test too

---
 gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
index 8a18f1a128..0dae9173ba 100644
--- a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
@@ -577,7 +577,7 @@ create_ctx(void)
     poa_phases[2*i].is_static_field = is_static_field_oap;
     poa_phases[2*i].fdot_mult_type = fdot_mult_type_oap;
     poa_phases[2*i].is_positivity_enabled = is_positivity_enabled_oap;
-    poa_phases[2*i].damping_type = GKYL_GK_DAMPING_LOW_PASS_FILTER;
+    poa_phases[2*i].damping_type = GKYL_GK_DAMPING_NONE;
 
     // FDPs.
     poa_phases[2*i+1].phase = GK_POA_FDP;
@@ -587,7 +587,7 @@ create_ctx(void)
     poa_phases[2*i+1].is_static_field = is_static_field_fdp;
     poa_phases[2*i+1].fdot_mult_type = fdot_mult_type_fdp;
     poa_phases[2*i+1].is_positivity_enabled = is_positivity_enabled_fdp;
-    poa_phases[2*i+1].damping_type = GKYL_GK_DAMPING_NONE;
+    poa_phases[2*i+1].damping_type = GKYL_GK_DAMPING_LOW_PASS_FILTER;
   }
   // Add an extra, longer FDP.
   poa_phases[num_phases-1].phase = GK_POA_FDP;

From 4b22dc1fd627da8fd22fc2f080020e5696d99b42 Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Wed, 8 Apr 2026 16:49:28 -0400
Subject: [PATCH 24/25] Clean up the reset of the fbar across phases. Valgrind
 clean

---
 gyrokinetic/apps/gk_species.c            |  5 +-
 gyrokinetic/apps/gk_species_damping.c    | 73 +++++++++---------------
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h |  3 +-
 3 files changed, 28 insertions(+), 53 deletions(-)

diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c
index 4981d55614..8be5ea819a 100644
--- a/gyrokinetic/apps/gk_species.c
+++ b/gyrokinetic/apps/gk_species.c
@@ -1505,7 +1505,7 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st
   gk_species_anomalous_diff_init(app, gks, &gks->anom_diff);
 
   // Damping term -nu*f on RHS.
-  gks->damping = (struct gk_damping) { .fbar_backup_initialized = false };
+  gks->damping = (struct gk_damping) { };
   gk_species_damping_init(app, gks, &gks->damping);
 
   // Function multiplying df/dt.
@@ -1915,9 +1915,6 @@ gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *gks
   gk_species_source_release(app, &gks->src);
 
   gk_species_damping_release(app, &gks->damping);
-  if (gks->damping.fbar_backup_initialized) {
-    gkyl_array_release(gks->damping.fbar_backup);
-  }
 
   gk_species_fdot_multiplier_release(app, &gks->fdot_mult);
 
diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index d6ccee72f1..70973dcc31 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -453,9 +453,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->fbar1 = 0;
   damp->fbarnew = 0;
   damp->fbar_host = 0;
-  if (!damp->fbar_backup_initialized) {
-    damp->fbar_backup = 0;
-  }
+  damp->fbar_initialized = false;
 
   if (!damp->type) {
     return;
@@ -495,6 +493,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
     damp->fbar1 = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
     damp->fbarnew = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
     damp->fbar_host = damp->fbar;
+    damp->fbar_initialized = true;
     gkyl_array_clear(damp->fbar, 0.0);
 
     if (damp->write_fbar) {
@@ -535,6 +534,11 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
 void
 gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp)
 {
+  // Release remembered fbar state even when damping is currently disabled.
+  if (damp->fbar_initialized) {
+    gkyl_array_release(damp->fbar);
+  }
+
   if (!damp->type) {
     return;
   }
@@ -542,9 +546,7 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
   if (app->use_gpu) {
     gkyl_array_release(damp->rate_host);
   }
-
   if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
-    gkyl_array_release(damp->fbar);
     gkyl_array_release(damp->fbar1);
     gkyl_array_release(damp->fbarnew);
     if (app->use_gpu && damp->write_fbar) {
@@ -553,62 +555,39 @@ gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct
   }
 }
 
-// Damping runtime reset.
-static bool
-gk_species_damping_can_restore_fbar(const struct gk_species *gks,
-  const struct gkyl_gyrokinetic_damping *damp_inp, const struct gkyl_array *fbar_backup)
-{
-  // Determines whether the provided backup of fbar state is compatible with the requested damping configuration, and can be used to preserve fbar state across a reset.
-  if (damp_inp->type != GKYL_GK_DAMPING_LOW_PASS_FILTER || !fbar_backup
-    || !damp_inp->do_not_reset_fbar) {
-    return false;
-  }
-
-  int fbar_ncomp = damp_inp->cellwise_const ? 1 : gks->basis.num_basis;
-  return fbar_backup->ncomp == fbar_ncomp && fbar_backup->size == gks->local_ext.volume;
-}
-
 void
 gk_species_damping_reset(gkyl_gyrokinetic_app *app, double tm, struct gk_species *gks,
   struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp)
 {
-  bool preserve_fbar = false;
+  bool has_fbar_backup = false;
   struct gkyl_array *fbar_backup = 0;
-  bool keep_lpf_state = damp_inp.type == GKYL_GK_DAMPING_LOW_PASS_FILTER
-    && damp_inp.do_not_reset_fbar;
 
-  if (keep_lpf_state) {
-    // Prefer latest LPF state if currently active, otherwise keep parked backup.
-    if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER && damp->fbar) {
-      fbar_backup = gkyl_array_acquire(damp->fbar);
-    }
-    else if (damp->fbar_backup_initialized) {
-      fbar_backup = damp->fbar_backup;
-      damp->fbar_backup = 0;
-      damp->fbar_backup_initialized = false;
+  // Create a reset-local backup if requested and there is restorable state.
+  if (damp_inp.do_not_reset_fbar) {
+    if (damp->fbar_initialized) {
+      has_fbar_backup = true;
+      fbar_backup = mkarr(app->use_gpu, damp->fbar->ncomp, damp->fbar->size);
+      gkyl_array_copy(fbar_backup, damp->fbar);
     }
-    preserve_fbar = gk_species_damping_can_restore_fbar(gks, &damp_inp, fbar_backup);
-  }
-  else if (damp->fbar_backup_initialized) {
-    // User disabled preservation: drop any parked state.
-    gkyl_array_release(damp->fbar_backup);
-    damp->fbar_backup = 0;
-    damp->fbar_backup_initialized = false;
   }
 
   gk_species_damping_release(app, damp);
 
-  damp->fbar_backup = fbar_backup;
-  damp->fbar_backup_initialized = fbar_backup != 0;
   gks->info.damping = damp_inp;
   gk_species_damping_init(app, gks, damp);
 
-  gk_species_damping_set_fbar_to_f(gks, damp, preserve_fbar ? damp->fbar_backup : gks->f);
+  gk_species_damping_set_fbar_to_f(gks, damp, gks->f);
 
-  // If backup has been consumed to initialize LPF state, drop parked handle.
-  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER && damp->fbar_backup_initialized) {
-    gkyl_array_release(damp->fbar_backup);
-    damp->fbar_backup = 0;
-    damp->fbar_backup_initialized = false;
+  if (has_fbar_backup) {
+    if (damp->fbar_initialized) {
+      gkyl_array_copy(damp->fbar, fbar_backup);
+      gkyl_array_copy(damp->fbar1, fbar_backup);
+      gkyl_array_copy(damp->fbarnew, fbar_backup);
+    }
+    else {
+      damp->fbar = gkyl_array_acquire(fbar_backup);
+      damp->fbar_initialized = true;
+    }
+    gkyl_array_release(fbar_backup);
   }
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 413547129a..9be93d65bb 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -808,11 +808,10 @@ struct gk_damping {
   struct gkyl_array *rate; // Damping rate.
   struct gkyl_array *rate_host; // Host copy for use in IO and projecting.
   // Low-pass filter variables
+  bool fbar_initialized; // Whether fbar is allocated/owned and should be released.
   struct gkyl_array *fbar; // Filtered/averaged distribution function.
   struct gkyl_array *fbar1, *fbarnew; // SSPRK3 stage arrays for filtered distribution.
   struct gkyl_array *fbar_host; // Host copy of fbar for use in IO.
-  struct gkyl_array *fbar_backup; // Parked backup retained across resets while LPF is inactive.
-  bool fbar_backup_initialized; // Whether fbar_backup holds a valid parked handle.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*write_rate_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);

From 7d80ccaeffd18ebc61a5b8d936029cb6775f9f9e Mon Sep 17 00:00:00 2001
From: Maxwell-Rosen <mrquell@gmail.com>
Date: Sun, 6 Sep 2026 10:47:23 -0400
Subject: [PATCH 25/25] =?UTF-8?q?=E2=80=BA=20When=20simulations=20are=20re?=
 =?UTF-8?q?started,=20the=20low-pass=20filter=20fbar=20should=20be=20resta?=
 =?UTF-8?q?rted=20as=20well=20from=20the=20file.=20Implement=20this?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

• Implemented fbar restart support.

  - Reads the matching damping_fbar frame file.
  - Synchronizes CPU/GPU state and SSP-RK stage arrays.
  - Supports single- and multi-block restarts.
  - Properly propagates missing/corrupt restart-file errors.
  - Falls back to restarted f when fbar output is disabled.

  Verified with reduced WHAM restart and Valgrind: 0 errors, 0 leaks. Build compiled and linked; only the sandbox-blocked external ADAS install copy failed.
---
 gyrokinetic/apps/gk_species_damping.c    | 34 ++++++++++++++++++++++++
 gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 17 ++++++++++++
 gyrokinetic/apps/gyrokinetic.c           | 19 ++++++++++++-
 gyrokinetic/apps/gyrokinetic_multib.c    | 10 ++++++-
 4 files changed, 78 insertions(+), 2 deletions(-)

diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 40639a3bd9..753da55ec2 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -36,6 +36,38 @@ gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping
   damp->set_fbar_to_f_func(gks, damp, f);
 }
 
+static enum gkyl_array_rio_status
+gk_species_damping_read_fbar_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  const char *fname)
+{
+  return GKYL_ARRAY_RIO_SUCCESS;
+}
+
+static enum gkyl_array_rio_status
+gk_species_damping_read_fbar_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  const char *fname)
+{
+  struct gk_damping *damp = &gks->damping;
+  enum gkyl_array_rio_status status =
+    gkyl_comm_array_read(gks->comm, &gks->grid, &gks->local, damp->fbar_host, fname);
+
+  if (status == GKYL_ARRAY_RIO_SUCCESS) {
+    if (app->use_gpu)
+      gkyl_array_copy(damp->fbar, damp->fbar_host);
+    gkyl_array_copy(damp->fbar1, damp->fbar);
+    gkyl_array_copy(damp->fbarnew, damp->fbar);
+  }
+
+  return status;
+}
+
+enum gkyl_array_rio_status
+gk_species_damping_read_fbar(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  const char *fname)
+{
+  return gks->damping.read_fbar_func(app, gks, fname);
+}
+
 // Damping diagnostics write helpers.
 
 void
@@ -441,6 +473,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
   damp->write_func = gk_species_damping_write_disabled;
   damp->write_rate_func = gk_species_damping_write_rate_disabled;
   damp->write_fbar_func = gk_species_damping_write_fbar_disabled;
+  damp->read_fbar_func = gk_species_damping_read_fbar_disabled;
   damp->advance_func = gk_species_damping_advance_disabled;
   damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_disabled;
   damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_disabled;
@@ -497,6 +530,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
     gkyl_array_clear(damp->fbar, 0.0);
 
     if (damp->write_fbar) {
+      damp->read_fbar_func = gk_species_damping_read_fbar_enabled;
       if (damp->cellwise_const) {
         damp->write_fbar_func = gk_species_damping_write_fbar_cellwise_const;
       }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 3044dddb28..5dafb412fe 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -856,6 +856,8 @@ struct gk_damping {
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*write_rate_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*write_fbar_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
+  enum gkyl_array_rio_status (*read_fbar_func)(gkyl_gyrokinetic_app* app,
+    struct gk_species *gks, const char *fname);
   void (*advance_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
     const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
     struct gkyl_array *f_buffer,
@@ -2973,6 +2975,21 @@ void gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_speci
  */
 void gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
 
+/**
+ * Read the filtered distribution function used by low-pass-filter damping.
+ *
+ * This is a runtime-dispatched no-op unless fbar output is enabled. On a
+ * successful read, all filtered-distribution RK stage arrays are synchronized
+ * to the restarted state.
+ *
+ * @param app Gyrokinetic app object.
+ * @param gks Species object.
+ * @param fname Name of the fbar restart file.
+ * @return Status of the file read.
+ */
+enum gkyl_array_rio_status gk_species_damping_read_fbar(gkyl_gyrokinetic_app *app,
+  struct gk_species *gks, const char *fname);
+
 /**
  * Compute filtered distribution RHS for low-pass filter damping:
  *   d(fbar)/dt = rate * (f - fbar).
diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c
index bf7cb60a7b..c00491364c 100644
--- a/gyrokinetic/apps/gyrokinetic.c
+++ b/gyrokinetic/apps/gyrokinetic.c
@@ -3332,6 +3332,8 @@ gkyl_gyrokinetic_app_from_file_species(gkyl_gyrokinetic_app *app, int sidx,
       gkyl_array_copy(gk_s->f, gk_s->f_host);
 
     if (rstat.io_status == GKYL_ARRAY_RIO_SUCCESS) {
+      // Use f as a valid fallback filtered state when no fbar file is configured.
+      gk_species_damping_set_fbar_to_f(gk_s, &gk_s->damping, gk_s->f);
       gk_species_source_calc(app, gk_s, &gk_s->src, gk_s->lte.f_lte, 0.0);
       // Read volume and time integrated boundary flux diagnostics.
       gk_species_bflux_read_voltime_integrated_mom(app, gk_s, &gk_s->bflux);
@@ -3385,6 +3387,13 @@ gkyl_gyrokinetic_app_from_frame_species(gkyl_gyrokinetic_app *app, int sidx, int
   struct gkyl_app_restart_status rstat = gkyl_gyrokinetic_app_from_file_species(app, sidx, fileNm.str);
   cstr_drop(&fileNm);
 
+  if (rstat.io_status == GKYL_ARRAY_RIO_SUCCESS) {
+    cstr fbarFileNm = cstr_from_fmt("%s-%s_damping_fbar_%d.gkyl", app->name,
+      gk_s->info.name, frame);
+    rstat.io_status = gk_species_damping_read_fbar(app, gk_s, fbarFileNm.str);
+    cstr_drop(&fbarFileNm);
+  }
+
   // Append to existing integrated diagnostics.
   app->is_first_dt_write_call = false;
   gk_s->is_first_integ_write_call = false;
@@ -3432,13 +3441,19 @@ gkyl_gyrokinetic_app_from_frame_neut_species(gkyl_gyrokinetic_app *app, int sidx
 struct gkyl_app_restart_status
 gkyl_gyrokinetic_app_read_from_frame(gkyl_gyrokinetic_app *app, int frame)
 {
-  struct gkyl_app_restart_status rstat;
+  struct gkyl_app_restart_status rstat = {
+    .io_status = GKYL_ARRAY_RIO_SUCCESS,
+    .frame = frame,
+    .stime = 0.0,
+  };
   for (int i=0; i<app->num_neut_species; i++) {
     if (app->neut_species[i].info.is_static) {
       gk_neut_species_apply_ic(app, &app->neut_species[i], 0.0);
     }
     else {
       rstat = gkyl_gyrokinetic_app_from_frame_neut_species(app, i, frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
   }
   for (int i=0; i<app->num_species; i++) {
@@ -3447,6 +3462,8 @@ gkyl_gyrokinetic_app_read_from_frame(gkyl_gyrokinetic_app *app, int frame)
     }
     else {
       rstat = gkyl_gyrokinetic_app_from_frame_species(app, i, frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
   }
 
diff --git a/gyrokinetic/apps/gyrokinetic_multib.c b/gyrokinetic/apps/gyrokinetic_multib.c
index de0352dee2..99a1558bd5 100644
--- a/gyrokinetic/apps/gyrokinetic_multib.c
+++ b/gyrokinetic/apps/gyrokinetic_multib.c
@@ -1067,7 +1067,11 @@ gkyl_gyrokinetic_multib_app_from_file_neut_species(gkyl_gyrokinetic_multib_app *
 struct gkyl_app_restart_status
 gkyl_gyrokinetic_multib_app_read_from_frame(gkyl_gyrokinetic_multib_app *app, int frame)
 {
-  struct gkyl_app_restart_status rstat;
+  struct gkyl_app_restart_status rstat = {
+    .io_status = GKYL_ARRAY_RIO_SUCCESS,
+    .frame = frame,
+    .stime = 0.0,
+  };
   for (int b=0; b<app->num_local_blocks; ++b) {
     struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
     for (int i=0; i<app->num_neut_species; i++) {
@@ -1076,9 +1080,13 @@ gkyl_gyrokinetic_multib_app_read_from_frame(gkyl_gyrokinetic_multib_app *app, in
         neut_frame = 0;
       }
       rstat = gkyl_gyrokinetic_app_from_frame_neut_species(sbapp, i, neut_frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
     for (int i=0; i<app->num_species; i++) {
       rstat = gkyl_gyrokinetic_app_from_frame_species(sbapp, i, frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
   }