Floor cross species thermal velocities at zero

gyrokinetic/apps/gk_species_lbo.c

@@ -107,43 +107,55 @@ gklbo_cross_moms_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gk
   // Empty method.
 }
 
+// Compute unscaled cross primitive moments (u_par_sr, vtsq_sr) for cross-collision species coll_idx.
+// Result is written into lbo->nu_boundary_corrections, which is used as scratch.
+static void
+gklbo_calc_cross_prim_moms(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_lbo_collisions *lbo, int coll_idx)
+{
+  lbo->alpha_E_func(app, gks, lbo, coll_idx);
+
+  // Multiply moments and boundary corrections by cross nu.
+  for (int d=0; d<3; d++)
+    gkyl_dg_mul_op(&app->basis, d, lbo->nu_moms, d, lbo->moms.marr, 0, lbo->cross_nu[coll_idx]);
+  for (int d=0; d<2; d++)
+    gkyl_dg_mul_op(&app->basis, d, lbo->nu_boundary_corrections, d, lbo->boundary_corrections, 0, lbo->cross_nu[coll_idx]);
+
+  // Compute cross primitive moments.
+  // Recycle the boundary_corrections array because we don't need those anymore.
+  gkyl_prim_lbo_cross_calc_advance(lbo->cross_calc, &app->local, lbo->alpha_E,
+    gks->info.mass, lbo->nu_moms, lbo->prim_moms,
+    lbo->other_m[coll_idx], lbo->collide_with[coll_idx]->lbo.moms.marr, lbo->other_prim_moms[coll_idx],
+    lbo->nu_boundary_corrections, lbo->cross_nu[coll_idx], lbo->nu_boundary_corrections);
+
+    struct gkyl_array *cross_prim_moms = lbo->nu_boundary_corrections;
+
+    // Floor the cell-average of vtSq_{sr} to zero, setting the rest of the coefficients to zero as well
+    gkyl_array_set_offset_range(lbo->cross_vtsq, 1.0, cross_prim_moms, 1*app->basis.num_basis, &app->local);
+    gkyl_array_max_by_cell_per_cell_avg_range(lbo->cross_vtsq, lbo->cross_vtsq_floor, &app->local);
+    gkyl_array_set_offset_range(cross_prim_moms, 1.0, lbo->cross_vtsq, 1*app->basis.num_basis, &app->local);
+}
+
 static void
 gklbo_cross_moms_enabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_lbo_collisions *lbo)
 {
-  // Compute primitive moments for cross-species collisions.
   struct timespec wst = gkyl_wall_clock();
-
-  for (int i=0; i<lbo->num_cross_collisions; ++i) {
 
-    // Compute the cross-species collision frequency.
+  for (int i=0; i<lbo->num_cross_collisions; ++i) {
     lbo->cross_nu_func(app, gks, lbo, i);
+    gklbo_calc_cross_prim_moms(app, gks, lbo, i);
 
-    // Compute alpha_E.
-    lbo->alpha_E_func(app, gks, lbo, i);
-
-    // Multiply moments and boundary corrections by cross nu.
-    for (int d=0; d<3; d++)
-      gkyl_dg_mul_op(&app->basis, d, lbo->nu_moms, d, lbo->moms.marr, 0, lbo->cross_nu[i]);
-    for (int d=0; d<2; d++)
-      gkyl_dg_mul_op(&app->basis, d, lbo->nu_boundary_corrections, d, lbo->boundary_corrections, 0, lbo->cross_nu[i]);
-
-    // Compute cross primitive moments.
-    // Recycle the boundary_corrections array because we don't need those anymore.
+    // Scale upar_{sr} and vtSq_{sr} by nu_{sr} and accumulate.
     struct gkyl_array *cross_prim_moms = lbo->nu_boundary_corrections;
-    gkyl_prim_lbo_cross_calc_advance(lbo->cross_calc, &app->local, lbo->alpha_E, 
-      gks->info.mass, lbo->nu_moms, lbo->prim_moms,
-      lbo->other_m[i], lbo->collide_with[i]->lbo.moms.marr, lbo->other_prim_moms[i],
-      lbo->nu_boundary_corrections, lbo->cross_nu[i], cross_prim_moms);
 
     // Scale upar_{sr} and vtSq_{sr} by nu_{sr}.
     for (int d=0; d<2; d++)
       gkyl_dg_mul_op(&app->basis, d, cross_prim_moms, d, cross_prim_moms, 0, lbo->cross_nu[i]);
 
     gkyl_array_accumulate(lbo->nu_prim_moms, 1.0, cross_prim_moms);
-
   }
-  app->stat.species_coll_mom_tm += gkyl_time_diff_now_sec(wst);    
+  app->stat.species_coll_mom_tm += gkyl_time_diff_now_sec(wst);
 }
 
 static void
@@ -194,6 +206,14 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
   const struct gkyl_msgpack_map_elem* io_meta_nu_prim[] = {gks->io_meta_conf, app->gk_geom->io_meta_basic, desc_nu_prim};
   struct gkyl_msgpack_data *mt_nu_prim = gkyl_msgpack_create_union(sizeof(io_meta_nu_prim_len)/sizeof(int), io_meta_nu_prim_len, io_meta_nu_prim);
 
+  struct gkyl_msgpack_map_elem desc_nu_cross[] = {
+    { .key = "Description", .elem_type = GKYL_MP_STRING,
+      .cval = "Cross-species drift velocity and thermal speed squared, or cross-species collision frequency." }
+  };
+  int io_meta_nu_cross_len[] = {app->io_meta_dg_len, app->gk_geom->io_meta_basic_len, 1};
+  const struct gkyl_msgpack_map_elem* io_meta_nu_cross[] = {app->io_meta_dg, app->gk_geom->io_meta_basic, desc_nu_cross};
+  struct gkyl_msgpack_data *mt_nu_cross = gkyl_msgpack_create_union(sizeof(io_meta_nu_cross_len)/sizeof(int), io_meta_nu_cross_len, io_meta_nu_cross);
+
   // Write out nu_sum and nu_prim_moms.
   const char *fmt = "%s-%s_lbo_nu_sum_%d.gkyl";
   int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name, frame);
@@ -215,8 +235,68 @@ gklbo_write_mom_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, doubl
   gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt_nu_prim, gks->lbo.nu_prim_moms_host, fileNm_nu_prim);
   app->stat.n_diag_io += 2;
 
+  // Per-cross-species diagnostics: unscaled cross primitive moments and cross-nu.
+  if (gks->lbo.write_diagnostics) {
+    // Self primitive moments (u_par_s, vtsq_s) before any cross-species contribution.
+    const char *fmt_prim = "%s-%s_lbo_self_prim_moms_%d.gkyl";
+    int sz_prim = gkyl_calc_strlen(fmt_prim, app->name, gks->info.name, frame);
+    char fileNm_prim[sz_prim+1];
+    snprintf(fileNm_prim, sizeof fileNm_prim, fmt_prim, app->name, gks->info.name, frame);
+    struct gkyl_array *prim_moms_io;
+    if (app->use_gpu) {
+      gkyl_array_copy(gks->lbo.nu_prim_moms_host, gks->lbo.prim_moms);
+      prim_moms_io = gks->lbo.nu_prim_moms_host;
+    } else {
+      prim_moms_io = gks->lbo.prim_moms;
+    }
+    struct gkyl_msgpack_map_elem desc_prim[] = {
+      { .key = "Description", .elem_type = GKYL_MP_STRING,
+        .cval = "Self drift velocity and thermal speed squared, before cross-species contributions." }
+    };
+    int io_meta_prim_len[] = {app->io_meta_dg_len, app->gk_geom->io_meta_basic_len, 1};
+    const struct gkyl_msgpack_map_elem* io_meta_prim[] = {app->io_meta_dg, app->gk_geom->io_meta_basic, desc_prim};
+    struct gkyl_msgpack_data *mt_prim = gkyl_msgpack_create_union(sizeof(io_meta_prim_len)/sizeof(int), io_meta_prim_len, io_meta_prim);
+    gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt_prim, prim_moms_io, fileNm_prim);
+    gkyl_msgpack_data_release(mt_prim);
+    app->stat.n_diag_io += 1;
+    for (int i=0; i<gks->lbo.num_cross_collisions; ++i) {
+      const char *other_name = gks->lbo.collide_with[i]->info.name;
+
+      const char *fmt_cprim = "%s-%s_lbo_cross_%s_prim_moms_%d.gkyl";
+      int sz_cprim = gkyl_calc_strlen(fmt_cprim, app->name, gks->info.name, other_name, frame);
+      char fileNm_cprim[sz_cprim+1];
+      snprintf(fileNm_cprim, sizeof fileNm_cprim, fmt_cprim, app->name, gks->info.name, other_name, frame);
+
+      const char *fmt_cnu = "%s-%s_lbo_cross_%s_nu_%d.gkyl";
+      int sz_cnu = gkyl_calc_strlen(fmt_cnu, app->name, gks->info.name, other_name, frame);
+      char fileNm_cnu[sz_cnu+1];
+      snprintf(fileNm_cnu, sizeof fileNm_cnu, fmt_cnu, app->name, gks->info.name, other_name, frame);
+
+      gklbo_calc_cross_prim_moms(app, gks, &gks->lbo, i);
+      struct gkyl_array *cross_prim_moms_scratch = gks->lbo.nu_boundary_corrections;
+
+      struct gkyl_array *cross_prim_moms_io, *cross_nu_io;
+      if (app->use_gpu) {
+        gkyl_array_copy(gks->lbo.nu_prim_moms_host, cross_prim_moms_scratch);
+        cross_prim_moms_io = gks->lbo.nu_prim_moms_host;
+        gkyl_array_copy(gks->lbo.nu_sum_host, gks->lbo.cross_nu[i]);
+        cross_nu_io = gks->lbo.nu_sum_host;
+      }
+      else {
+        cross_prim_moms_io = cross_prim_moms_scratch;
+        cross_nu_io = gks->lbo.cross_nu[i];
+      }
+
+      gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt_nu_cross, cross_prim_moms_io, fileNm_cprim);
+      gkyl_comm_array_write(app->comm, &app->grid, &app->local, mt_nu_cross, cross_nu_io, fileNm_cnu);
+      app->stat.n_diag_io += 2;
+    }
+  }
+
   gkyl_msgpack_data_release(mt_nu_sum);
   gkyl_msgpack_data_release(mt_nu_prim);
+  gkyl_msgpack_data_release(mt_nu_cross);
+
   app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wtm);
 }
 
@@ -373,6 +453,12 @@ gk_species_lbo_cross_init(struct gkyl_gyrokinetic_app *app, struct gk_species *g
 
       // Morse's alpha_E.
       lbo->alpha_E = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+
+      // Scratch space for flooring the vtSq component of the cross-prim moments.
+      lbo->cross_vtsq = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+      lbo->cross_vtsq_floor = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+      gkyl_array_clear(lbo->cross_vtsq_floor, 0.0); // Floor value of zero.
+
       for (int i=0; i<lbo->num_cross_collisions; ++i) {
         // Cross-species collision frequency, nu_sr.
         lbo->cross_nu[i] = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
@@ -509,6 +595,8 @@ gk_species_lbo_release(const struct gkyl_gyrokinetic_app *app, const struct gk_l
         gkyl_array_release(lbo->cross_nu[i]);
 
       gkyl_array_release(lbo->alpha_E);
+      gkyl_array_release(lbo->cross_vtsq);
+      gkyl_array_release(lbo->cross_vtsq_floor);
     }
 
     gkyl_dg_updater_lbo_gyrokinetic_release(lbo->coll_slvr);

gyrokinetic/apps/gkyl_gyrokinetic_priv.h

@@ -322,10 +322,10 @@ struct gk_collisionless {
     struct gk_collisionless *gkcls, double tm, int frame);
 };
 
-struct gk_lbo_collisions {  
+struct gk_lbo_collisions {
   enum gkyl_collision_id collision_id; // type of collisions
   bool write_diagnostics; // Whether to write diagnostics out.
-  bool not_in_dfdt; // Whether to not add collision contribution to df/dt 
+  bool not_in_dfdt; // Whether to not add collision contribution to df/dt
 
   struct gkyl_array *self_nu; // Self-collision frequency.
   struct gkyl_array *boundary_corrections; // LBO boundary corrections.
@@ -355,6 +355,8 @@ struct gk_lbo_collisions {
   struct gkyl_array *cross_nu[GKYL_MAX_SPECIES]; // Cross-species collision frequencies.
   struct gkyl_array *cross_nu_prim_moms; // Weak multiplication of collision frequency and primitive moments.
   struct gkyl_array *alpha_E; // Morse's alpha_E factor.
+  struct gkyl_array *cross_vtsq; // Scratch space holding the vtSq component of the cross-prim moments.
+  struct gkyl_array *cross_vtsq_floor; // Zeroed array used to floor cross vtSq cell-average to >= 0.
   gkyl_prim_lbo_cross_calc *cross_calc; // LBO cross-primitive moment calculator
 
   struct gk_species_moment moms; // Moments needed in LBO (M0, M1, M2).