Add omegaH_dt diagnostic to gyrokinetic species

gyrokinetic/apps/gk_species.c

@@ -494,6 +494,14 @@ gk_species_calc_integrated_mom_dynamic(gkyl_gyrokinetic_app* app, struct gk_spec
   }
   gkyl_dynvec_append(gks->integ_diag, tm, avals_global);
 
+  if (gks->info.omegaH_dt_diagnostic) {
+    double omegaH_dt_local = gk_species_omegaH_dt(app, gks, gks->f);
+    double omegaH_dt_global = DBL_MAX;
+    gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MIN, 1,
+      &omegaH_dt_local, &omegaH_dt_global);
+    gkyl_dynvec_append(gks->omegaH_dt, tm, &omegaH_dt_global);
+  }
+
   if (gks->info.time_rate_diagnostics) {
     // Reduce (sum) over whole domain, append to diagnostics.
     gkyl_array_accumulate(gks->fdot_mom_new, -1.0, gks->fdot_mom_old);
@@ -551,6 +559,37 @@ gk_species_write_integrated_mom_dynamic(gkyl_gyrokinetic_app *app, struct gk_spe
   }
   gkyl_dynvec_clear(gks->integ_diag);
   app->stat.n_diag_io += 1;
+
+  if (gks->info.omegaH_dt_diagnostic) {
+    if (rank == 0) {
+      const char *fmt = "%s-%s_omegaH_dt.gkyl";
+      int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name);
+      char fileNm[sz+1];
+      snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name);
+
+      if (gks->is_first_omegaH_dt_write_call) {
+        struct gkyl_msgpack_map_elem io_meta_omegaH_dt[] = {
+          { .key = "Description", .elem_type = GKYL_MP_STRING,
+            .cval = "Stable time step for the omega_H mode for this charged species." }
+        };
+        int io_meta_len[] = {gks->io_meta_basic_len, app->gk_geom->io_meta_basic_len, 1};
+        const struct gkyl_msgpack_map_elem* io_meta[] = {
+          gks->io_meta_basic, app->gk_geom->io_meta_basic, io_meta_omegaH_dt
+        };
+        struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union(
+          sizeof(io_meta_len)/sizeof(int), io_meta_len, io_meta);
+
+        gkyl_dynvec_write_wmeta(gks->omegaH_dt, fileNm, mt);
+        gks->is_first_omegaH_dt_write_call = false;
+        gkyl_msgpack_data_release(mt);
+      }
+      else {
+        gkyl_dynvec_awrite(gks->omegaH_dt, fileNm);
+      }
+    }
+    gkyl_dynvec_clear(gks->omegaH_dt);
+    app->stat.n_diag_io += 1;
+  }
 
   if (gks->info.time_rate_diagnostics) {
     if (rank == 0) {
@@ -722,6 +761,8 @@ gk_species_release_dynamic(const gkyl_gyrokinetic_app* app, const struct gk_spec
 
   // Release integrated diag memory.
   gkyl_dynvec_release(s->integ_diag);
+  if (s->info.omegaH_dt_diagnostic)
+    gkyl_dynvec_release(s->omegaH_dt);
   if (app->use_gpu) {
     gkyl_cu_free(s->red_integ_diag);
     gkyl_cu_free(s->red_integ_diag_global);
@@ -805,6 +846,11 @@ gk_species_init_dynamic(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app
   gks->integ_diag = gkyl_dynvec_new(GKYL_DOUBLE, gks->integ_moms.num_mom);
   gks->is_first_integ_write_call = true;
 
+  if (gks->info.omegaH_dt_diagnostic) {
+    gks->omegaH_dt = gkyl_dynvec_new(GKYL_DOUBLE, 1);
+    gks->is_first_omegaH_dt_write_call = true;
+  }
+
   // Allocate dynamic-vector to store Delta f integrated moments.
   if (gks->info.time_rate_diagnostics) {
     gks->fdot_mom_old = mkarr(app->use_gpu, gks->integ_moms.marr->ncomp, gks->integ_moms.marr->size);

gyrokinetic/apps/gkyl_gyrokinetic.h

@@ -452,6 +452,7 @@ struct gkyl_gyrokinetic_species {
   enum gkyl_distribution_moments integrated_diag_moments[12]; // List of integrated diagnostic moments.
   bool time_rate_diagnostics; // Whether to ouput df/dt diagnostics.
   bool write_omega_cfl; // Whether to ouput dt diagnostic for the CFL constraint.
+  bool omegaH_dt_diagnostic; // Whether to output the omega_H stable time step as a dynvector.
 
   struct gkyl_gyrokinetic_collisionless collisionless; // Collisionless terms.
 

gyrokinetic/apps/gkyl_gyrokinetic_multib.h

@@ -33,6 +33,7 @@ struct gkyl_gyrokinetic_multib_species {
   int num_integrated_diag_moments; // Number of integrated diagnostic moments.
   enum gkyl_distribution_moments integrated_diag_moments[12]; // List of integrated diagnostic moments.
   bool time_rate_diagnostics; // Whether to ouput df/dt diagnostics.
+  bool omegaH_dt_diagnostic; // Whether to output the omega_H stable time step as a dynvector.
 
   struct gkyl_phase_diagnostics_inp boundary_flux_diagnostics;
 

gyrokinetic/apps/gkyl_gyrokinetic_priv.h

@@ -1047,6 +1047,9 @@ struct gk_species {
   gkyl_dynvec integ_diag; // Integrated moments reduced across grid
   bool is_first_integ_write_call;  // Whether dynvec is being written for the first time.
 
+  gkyl_dynvec omegaH_dt; // omega_H stable time step reduced across MPI ranks.
+  bool is_first_omegaH_dt_write_call; // Whether omega_H dt dynvec is being written for the first time.
+
   struct gkyl_array *fdot_mom_old, *fdot_mom_new; // Moments of f_old and f_new.
   gkyl_dynvec fdot_integ_diag; // Integrated moments of Delta f=f_new - f_old..
   bool is_first_fdot_integ_write_call; // Whether dynvec is being written for the first time.

gyrokinetic/apps/gyrokinetic.c

@@ -3387,6 +3387,8 @@ gkyl_gyrokinetic_app_from_frame_species(gkyl_gyrokinetic_app *app, int sidx, int
   // Append to existing integrated diagnostics.
   app->is_first_dt_write_call = false;
   gk_s->is_first_integ_write_call = false;
+  if (gk_s->info.omegaH_dt_diagnostic)
+    gk_s->is_first_omegaH_dt_write_call = false;
   gk_s->is_first_L2norm_write_call = false;
   for (int b=0; b<gk_s->bflux.num_boundaries; ++b)
     gk_s->bflux.is_first_intmom_write_call[b] = false;

gyrokinetic/apps/gyrokinetic_lw.c

@@ -356,6 +356,7 @@ gyrokinetic_species_lw_new(lua_State *L)
   gk_species.charge = glua_tbl_get_number(L, "charge", 0.0);
   gk_species.mass = glua_tbl_get_number(L, "mass", 1.0);
   gk_species.polarization_density = glua_tbl_get_number(L, "polarizationDensity", 0.0);
+  gk_species.omegaH_dt_diagnostic = glua_tbl_get_bool(L, "omegaHDtDiagnostic", false);
 
   with_lua_tbl_tbl(L, "cells") {
     vdim = glua_objlen(L);

gyrokinetic/apps/gyrokinetic_multib.c

@@ -176,6 +176,7 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid,
       species_inp.integrated_diag_moments[n] = sp->integrated_diag_moments[n];
     }
     species_inp.time_rate_diagnostics = sp->time_rate_diagnostics;
+    species_inp.omegaH_dt_diagnostic = sp->omegaH_dt_diagnostic;
     species_inp.boundary_flux_diagnostics = sp->boundary_flux_diagnostics;
 
     // Choose proper block-specific species input.