Orbit average phase for tandem mirrors

gyrokinetic/apps/gk_species.c

@@ -136,6 +136,7 @@ gk_species_rhs_dynamic(gkyl_gyrokinetic_app *app, struct gk_species *species,
 
   // Enforce the omega_H constraint on dt.
   double dt_omegaH = gk_species_omegaH_dt(app, species, fin);
+  gk_species_fdot_multiplier_advance_times_omegaH(app, species, &species->fdot_mult, &dt_omegaH);
   dt_out = fmin(dt_out, dt_omegaH);
 
   app->stat.species_omega_cfl_tm += gkyl_time_diff_now_sec(tm);

gyrokinetic/apps/gk_species_damping.c

@@ -121,14 +121,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.
@@ -183,7 +183,7 @@ gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks
       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);
+        app->field->phi_smooth, damp->phi_m_global, damp->phi_m_global, damp->rate);
       // Multiply by the user's scaling profile.
       gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
     }
@@ -218,7 +218,7 @@ gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *g
 
       // 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);
+        phi, damp->phi_m_global, damp->phi_m_global, damp->rate);
 
       // Assemble the damping term -scale_prof * mask * f.
       gkyl_array_set(f_buffer, 1.0, fin);

gyrokinetic/apps/gk_species_fdot_multiplier.c

@@ -53,6 +53,14 @@ gk_species_fdot_multiplier_advance_mult(gkyl_gyrokinetic_app *app, const struct
   gkyl_array_scale_by_cell(out, fdmul->multiplier);
 }
 
+void
+gk_species_fdot_multiplier_advance_omegaH_mult(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_fdot_multiplier *fdmul, double *out)
+{
+  // Multiply out by the multplier.
+  out[0] = out[0] / gks->collisionless.scale_fac;
+}
+
 void
 gk_species_fdot_multiplier_advance_loss_cone_mult(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_fdot_multiplier *fdmul, const struct gkyl_array *phi, struct gkyl_array *out)
@@ -62,9 +70,18 @@ gk_species_fdot_multiplier_advance_loss_cone_mult(gkyl_gyrokinetic_app *app, con
     app->basis_on_dev, &app->grid, &app->local, fdmul->phi_m);
   gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, fdmul->phi_m, fdmul->phi_m_global);
 
-  // Project the loss cone mask.
+  // For tandem mirrors, calculate the inner potential. For non-tandem, duplicate phi_m for safety.
+  if (fdmul->is_tandem) {
+    gkyl_dg_basis_ops_eval_array_at_coord_comp(phi, fdmul->bmag_tandem_coord,
+      app->basis_on_dev, &app->grid, &app->local, fdmul->phi_tandem);
+    gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, fdmul->phi_tandem, fdmul->phi_tandem_global);
   gkyl_loss_cone_mask_gyrokinetic_advance(fdmul->lcm_proj_op, &gks->local, &app->local,
-    phi, fdmul->phi_m_global, fdmul->multiplier);
+      phi, fdmul->phi_m_global, fdmul->phi_tandem_global, fdmul->multiplier);
+  }
+  else {
+    gkyl_loss_cone_mask_gyrokinetic_advance(fdmul->lcm_proj_op, &gks->local, &app->local,
+      phi, fdmul->phi_m_global, fdmul->phi_m_global, fdmul->multiplier);
+  }
 
   // Multiply out by the multplier.
   gkyl_array_scale_by_cell(out, fdmul->multiplier);
@@ -76,6 +93,13 @@ gk_species_fdot_multiplier_advance_disabled(gkyl_gyrokinetic_app *app, const str
 {
 }
 
+void
+gk_species_fdot_multiplier_advance_omegaH_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_fdot_multiplier *fdmul, double *out)
+{
+}
+
+
 static void
 proj_on_basis_c2p_phase_func(const double *xcomp, double *xphys, void *ctx)
 {
@@ -102,6 +126,7 @@ gk_species_fdot_multiplier_init(struct gkyl_gyrokinetic_app *app, struct gk_spec
   // Default function pointers.
   fdmul->write_func = gk_species_fdot_multiplier_write_disabled;
   fdmul->advance_times_cfl_func = gk_species_fdot_multiplier_advance_disabled;
+  fdmul->advance_times_omegaH_func = gk_species_fdot_multiplier_advance_omegaH_disabled;
   fdmul->advance_times_rate_func = gk_species_fdot_multiplier_advance_disabled;
 
   if (fdmul->type) {
@@ -144,6 +169,7 @@ gk_species_fdot_multiplier_init(struct gkyl_gyrokinetic_app *app, struct gk_spec
       gkyl_array_copy(fdmul->multiplier, fdmul->multiplier_host);
 
       fdmul->advance_times_cfl_func = gk_species_fdot_multiplier_advance_mult;
+      fdmul->advance_times_omegaH_func = gk_species_fdot_multiplier_advance_omegaH_mult;
       fdmul->advance_times_rate_func = gk_species_fdot_multiplier_advance_mult;
       if (fdmul->write_diagnostics)
         fdmul->write_func = gk_species_fdot_multiplier_write_init_only;
@@ -158,47 +184,117 @@ gk_species_fdot_multiplier_init(struct gkyl_gyrokinetic_app *app, struct gk_spec
       enum gkyl_quad_type qtype = GKYL_GAUSS_LOBATTO_QUAD;
       int num_quad = gks->basis.poly_order+1; // This can be p+1 or 1. Must be
                                               // at leat p+1 for Gauss-Lobatto.
+      int cdim = app->cdim;
+
+      // Find the peaks of B magnitude using the position map infrastructure
+      struct gkyl_position_map_inp pmap_info = {
+        .id = GKYL_PMAP_CONSTANT_DB_NUMERIC,
+        .map_strength = 0.0,
+      };
+      struct gkyl_position_map *gpm = gkyl_position_map_new( pmap_info, app->grid, \
+        app->local, app->local_ext, app->global, app->global_ext, app->basis);
+
+      gkyl_position_map_set_bmag(gpm, app->comm, app->gk_geom->geo_int.bmag);
+      gkyl_position_map_deflated_find_B_extrema(gpm, app->grid, app->global);
+
+      // Extract variables from position map context
+      double *bmag_extrema = gpm->constB_ctx->bmag_extrema;
+      double *theta_extrema = gpm->constB_ctx->theta_extrema;
+      int num_extrema = gpm->constB_ctx->num_extrema;
+
+      bool is_symmetric, is_tandem;
+      if (gkyl_compare_double(app->grid.lower[cdim-1], -app->grid.upper[cdim-1], 1e-14)) {
+        is_symmetric = true;
+      }
+      else if (gkyl_compare_double(app->grid.lower[cdim-1], 0.0, 1e-14)){
+        is_symmetric = false;
+      }
+      else {
+        assert(false); // Needs either the lower bound at 0 or symmetric grid
+      }
 
-      // 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];
+      if ( (is_symmetric && num_extrema == 5) || (!is_symmetric && num_extrema == 3) ) {
+        is_tandem = false;
+      }
+      else if ((is_symmetric && num_extrema == 9) || (!is_symmetric && num_extrema == 5)) {
+        is_tandem = true;
       }
       else {
-        for (int d=0; d<app->cdim; d++)
-          bmag_max_coord_local[d] = -DBL_MAX;
+        assert(false); // Unsupported number of extrema for loss-cone multiplier
+      }
+
+      // Identify wall values
+      double bmag_wall = bmag_extrema[num_extrema - 1];
+      double theta_wall = theta_extrema[num_extrema - 1];
+      // Identify walls (endpoints) and find the actual mirror point
+      double bmag_mirror = bmag_extrema[num_extrema - 2];
+      double theta_mirror = theta_extrema[num_extrema - 2];
+      // Find the inner peak for tandem mirrors, or duplicate for non-tandem
+      double bmag_tandem, theta_tandem;
+      if (!is_symmetric && !is_tandem) {
+        bmag_tandem = bmag_mirror;
+        theta_tandem = theta_mirror;
+      }
+      else {
+        bmag_tandem = bmag_extrema[num_extrema - 4];
+        theta_tandem = theta_extrema[num_extrema - 4];
+      }
+
+      double coord_mirror[GKYL_MAX_CDIM];
+      double coord_wall[GKYL_MAX_CDIM];
+      double coord_tandem[GKYL_MAX_CDIM];
+      if (app->cdim==1){
+        coord_mirror[0] = theta_mirror;
+        coord_wall[0] = theta_wall;
+        coord_tandem[0] = theta_tandem;
+      }
+      else{
+        assert(false); // only implemented for 1D config space
       }
-      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, app->cdim, bmag_max_coord_local, bmag_max_coord_global);
 
       if (app->use_gpu) {
         fdmul->bmag_max = gkyl_cu_malloc(sizeof(double));
         fdmul->bmag_max_coord = gkyl_cu_malloc(app->cdim*sizeof(double));
-	gkyl_cu_memcpy(fdmul->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
-	gkyl_cu_memcpy(fdmul->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
+        fdmul->bmag_wall = gkyl_cu_malloc(sizeof(double));
+        fdmul->bmag_wall_coord = gkyl_cu_malloc(app->cdim*sizeof(double));
+        fdmul->bmag_tandem = gkyl_cu_malloc(sizeof(double));
+        fdmul->bmag_tandem_coord = gkyl_cu_malloc(app->cdim*sizeof(double));
+
+        gkyl_cu_memcpy(fdmul->bmag_max, &bmag_mirror, sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(fdmul->bmag_max_coord, &coord_mirror, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(fdmul->bmag_wall, &bmag_wall, sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(fdmul->bmag_wall_coord, &coord_wall, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(fdmul->bmag_tandem, &bmag_tandem, sizeof(double), GKYL_CU_MEMCPY_H2D);
+        gkyl_cu_memcpy(fdmul->bmag_tandem_coord, &coord_tandem, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
       }
       else {
         fdmul->bmag_max = gkyl_malloc(sizeof(double));
         fdmul->bmag_max_coord = gkyl_malloc(app->cdim*sizeof(double));
-	memcpy(fdmul->bmag_max, &bmag_max_global, sizeof(double));
-	memcpy(fdmul->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double));
+        fdmul->bmag_wall = gkyl_malloc(sizeof(double));
+        fdmul->bmag_wall_coord = gkyl_malloc(app->cdim*sizeof(double));
+        fdmul->bmag_tandem = gkyl_malloc(sizeof(double));
+        fdmul->bmag_tandem_coord = gkyl_malloc(app->cdim*sizeof(double));
+
+        memcpy(fdmul->bmag_max, &bmag_mirror, sizeof(double));
+        memcpy(fdmul->bmag_max_coord, &coord_mirror, app->cdim*sizeof(double));
+        memcpy(fdmul->bmag_wall, &bmag_wall, sizeof(double));
+        memcpy(fdmul->bmag_wall_coord, &coord_wall, app->cdim*sizeof(double));
+        memcpy(fdmul->bmag_tandem, &bmag_tandem, sizeof(double));
+        memcpy(fdmul->bmag_tandem_coord, &coord_tandem, app->cdim*sizeof(double));
       }
 
       // Electrostatic potential at bmag_max_coord.
       if (app->use_gpu) {
         fdmul->phi_m = gkyl_cu_malloc(sizeof(double));
         fdmul->phi_m_global = gkyl_cu_malloc(sizeof(double));
+        fdmul->phi_tandem = gkyl_cu_malloc(sizeof(double));
+        fdmul->phi_tandem_global = gkyl_cu_malloc(sizeof(double));
       }
       else {
         fdmul->phi_m = gkyl_malloc(sizeof(double));
         fdmul->phi_m_global = gkyl_malloc(sizeof(double));
+        fdmul->phi_tandem = gkyl_malloc(sizeof(double));
+        fdmul->phi_tandem_global = gkyl_malloc(sizeof(double));
       }
 
       // Operator that projects the loss cone mask.
@@ -213,6 +309,11 @@ gk_species_fdot_multiplier_init(struct gkyl_gyrokinetic_app *app, struct gk_spec
         .bmag = app->gk_geom->geo_int.bmag,
         .bmag_max = fdmul->bmag_max,
         .bmag_max_loc = fdmul->bmag_max_coord,
+        .bmag_wall = fdmul->bmag_wall,
+        .bmag_wall_loc = fdmul->bmag_wall_coord,
+        .bmag_tandem = fdmul->bmag_tandem,
+        .bmag_tandem_loc = fdmul->bmag_tandem_coord,
+        .is_tandem = is_tandem,
         .mass = gks->info.mass,
         .charge = gks->info.charge,
         .qtype = qtype,
@@ -225,11 +326,13 @@ gk_species_fdot_multiplier_init(struct gkyl_gyrokinetic_app *app, struct gk_spec
       fdmul->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew( &inp_proj );
 
       fdmul->advance_times_cfl_func = gk_species_fdot_multiplier_advance_loss_cone_mult;
+      fdmul->advance_times_omegaH_func = gk_species_fdot_multiplier_advance_omegaH_mult;
       fdmul->advance_times_rate_func = gk_species_fdot_multiplier_advance_mult;
       if (fdmul->write_diagnostics)
         fdmul->write_func = gk_species_fdot_multiplier_write_enabled;
       else
         gkyl_array_release(fdmul->multiplier_host);
+      gkyl_position_map_release(gpm);
     }
   }
 }
@@ -244,6 +347,17 @@ gk_species_fdot_multiplier_advance_times_cfl(gkyl_gyrokinetic_app *app, const st
 
   app->stat.species_fdot_mult_tm += gkyl_time_diff_now_sec(wst);
 }
+
+void
+gk_species_fdot_multiplier_advance_times_omegaH(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_fdot_multiplier *fdmul, double *out)
+{
+  struct timespec wst = gkyl_wall_clock();
+
+  fdmul->advance_times_omegaH_func(app, gks, fdmul, out);
+
+  app->stat.species_fdot_mult_tm += gkyl_time_diff_now_sec(wst);
+}
 
 void
 gk_species_fdot_multiplier_advance_times_rate(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
@@ -278,14 +392,26 @@ gk_species_fdot_multiplier_release(const struct gkyl_gyrokinetic_app *app, const
       if (app->use_gpu) {
         gkyl_cu_free(fdmul->bmag_max);
         gkyl_cu_free(fdmul->bmag_max_coord);
+        gkyl_cu_free(fdmul->bmag_wall);
+        gkyl_cu_free(fdmul->bmag_wall_coord);
+        gkyl_cu_free(fdmul->bmag_tandem);
+        gkyl_cu_free(fdmul->bmag_tandem_coord);
         gkyl_cu_free(fdmul->phi_m);
         gkyl_cu_free(fdmul->phi_m_global);
+        gkyl_cu_free(fdmul->phi_tandem);
+        gkyl_cu_free(fdmul->phi_tandem_global);
       }
       else {
         gkyl_free(fdmul->bmag_max);
         gkyl_free(fdmul->bmag_max_coord);
+        gkyl_free(fdmul->bmag_wall);
+        gkyl_free(fdmul->bmag_wall_coord);
+        gkyl_free(fdmul->bmag_tandem);
+        gkyl_free(fdmul->bmag_tandem_coord);
         gkyl_free(fdmul->phi_m);
         gkyl_free(fdmul->phi_m_global);
+        gkyl_free(fdmul->phi_tandem);
+        gkyl_free(fdmul->phi_tandem_global);
       }
       gkyl_loss_cone_mask_gyrokinetic_release(fdmul->lcm_proj_op);
     }

gyrokinetic/apps/gkyl_gyrokinetic_priv.h

@@ -828,19 +828,27 @@ struct gk_fdot_multiplier {
   enum gkyl_gyrokinetic_fdot_multiplier_type type; // Type of multiplicative function term.
   bool write_diagnostics; // Whether to write diagnostics out.
   bool evolve; // Whether the multiplicative function is time dependent.
+  bool is_tandem; //  Whether we are doing a tandem mirror.
   struct gkyl_array *multiplier; // Damping rate.
   struct gkyl_array *multiplier_host; // Host copy for use in IO and projecting.
   struct gk_proj_on_basis_c2p_func_ctx proj_on_basis_c2p_ctx; // c2p function context.
   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 *bmag_wall; // Magnetic field at the wall.
+  double *bmag_wall_coord; // Location of bmag_wall.
+  double *bmag_tandem; // Magnetic field at the tandem mirror (for 7-extrema case).
+  double *bmag_tandem_coord; // Location of bmag_tandem.
   double *phi_m, *phi_m_global; // Electrostatic potential at bmag_max.
+  double *phi_tandem, *phi_tandem_global; // Electrostatic potential at bmag_tandem.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
   void (*advance_times_rate_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
     struct gk_fdot_multiplier *fdmul, const struct gkyl_array *phi, struct gkyl_array *out);
   void (*advance_times_cfl_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
     struct gk_fdot_multiplier *fdmul, const struct gkyl_array *phi, struct gkyl_array *out);
+  void (*advance_times_omegaH_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_fdot_multiplier *fdmul, double *out);
 };
 
 struct gk_heating {
@@ -2841,13 +2849,23 @@ void gk_species_fdot_multiplier_init(struct gkyl_gyrokinetic_app *app, struct gk
  * @param gks Species object.
  * @param fdmul Species df/dt multiplier object.
  * @param phi Current electrostatic potential.
- * @param fin Current distribution function.
- * @param f_buffer Phase-space buffer.
  * @param out CFL rate to multiply.
  */
 void gk_species_fdot_multiplier_advance_times_cfl(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
   struct gk_fdot_multiplier *fdmul, const struct gkyl_array *phi, struct gkyl_array *out);
 
+
+/**
+ * Multiply the omegaH rate.
+ *
+ * @param app gyrokinetic app object.
+ * @param gks Species object.
+ * @param fdmul Species df/dt multiplier object.
+ * @param out omegaH rate to multiply.
+ */
+void gk_species_fdot_multiplier_advance_times_omegaH(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_fdot_multiplier *fdmul, double *out);
+
 /**
  * Multiply df/dt.
  *