Fokker-Planck Collisions for Vlasov Solver

gyrokinetic/zero/gkyl_prim_lbo_gyrokinetic_priv.h

@@ -31,7 +31,7 @@ static const gkyl_prim_lbo_gyrokinetic_kern_list ser_self_prim_kernels[] = {
   // 2x kernels
   { NULL, gyrokinetic_self_prim_moments_2x2v_ser_p1, gyrokinetic_self_prim_moments_2x2v_ser_p2 }, // 2
   // 3x kernels
-  { NULL, gyrokinetic_self_prim_moments_3x2v_ser_p1, gyrokinetic_self_prim_moments_3x2v_ser_p2 }, // 3
+  { NULL, gyrokinetic_self_prim_moments_3x2v_ser_p1, NULL }, // 3
 };
 
 // cross primitive moment kernel list
@@ -43,7 +43,7 @@ static const gkyl_prim_lbo_gyrokinetic_cross_kern_list ser_cross_prim_kernels[]
   // 2x kernels
   { NULL, gyrokinetic_cross_prim_moments_2x2v_ser_p1, gyrokinetic_cross_prim_moments_2x2v_ser_p2 }, // 2
   // 3x kernels
-  { NULL, gyrokinetic_cross_prim_moments_3x2v_ser_p1, gyrokinetic_cross_prim_moments_3x2v_ser_p2 }, // 3
+  { NULL, gyrokinetic_cross_prim_moments_3x2v_ser_p1, NULL }, // 3
 };
 
 struct prim_lbo_type_gyrokinetic {

vlasov/apps/gkyl_vlasov.h

@@ -45,6 +45,7 @@ struct gkyl_vlasov_projection {
 // Parameters for species collisions
 struct gkyl_vlasov_collisions {
   enum gkyl_collision_id collision_id; // type of collisions (see gkyl_eqn_type.h)
+  bool write_diagnostics; // Whether to write diagnostics out.
 
   void *ctx; // context for collision function
   // function for computing self-collision frequency
@@ -552,6 +553,16 @@ void gkyl_vlasov_app_write_species(gkyl_vlasov_app* app, int sidx, double tm, in
  */
 void gkyl_vlasov_app_write_species_lte(gkyl_vlasov_app* app, int sidx, double tm, int frame);
 
+/**
+ * Write species data to file for H, G, drag coeff, diff coeff.
+ * 
+ * @param app App object.
+ * @param sidx Index of species to initialize.
+ * @param tm Time-stamp
+ * @param frame Frame number
+ */
+void gkyl_vlasov_app_write_species_fpo(gkyl_vlasov_app* app, int sidx, double tm, int frame);
+
 /**
  * Write fluid species data to file. 
  * 

vlasov/apps/gkyl_vlasov_priv.h

@@ -34,11 +34,14 @@
 #include <gkyl_dg_canonical_pb.h>
 #include <gkyl_dg_canonical_pb_fluid.h>
 #include <gkyl_dg_euler.h>
+#include <gkyl_dg_fpo_vlasov_diff_coeff.h>
+#include <gkyl_dg_fpo_vlasov_drag_coeff.h>
 #include <gkyl_dg_maxwell.h>
 #include <gkyl_dg_updater_fluid.h>
 #include <gkyl_dg_updater_diffusion_fluid.h>
 #include <gkyl_dg_updater_diffusion_gen.h>
 #include <gkyl_dg_updater_lbo_vlasov.h>
+#include <gkyl_dg_updater_fpo_vlasov.h>
 #include <gkyl_dg_updater_rad_vlasov.h>
 #include <gkyl_dg_updater_moment.h>
 #include <gkyl_dg_updater_vlasov.h>
@@ -49,13 +52,17 @@
 #include <gkyl_dynvec.h>
 #include <gkyl_elem_type.h>
 #include <gkyl_eqn_type.h>
+#include <gkyl_fpo_vlasov_coeff_correct.h>
+#include <gkyl_fpo_vlasov_coeff_recovery.h>
+#include <gkyl_fpo_proj_maxwellian_pots_on_basis.h>
 #include <gkyl_eval_on_nodes.h>
 #include <gkyl_fem_poisson.h>
 #include <gkyl_ghost_surf_calc.h>
 #include <gkyl_hyper_dg.h>
 #include <gkyl_mom_bcorr_lbo_vlasov.h>
 #include <gkyl_mom_calc.h>
 #include <gkyl_mom_calc_bcorr.h>
+#include <gkyl_mom_fpo_vlasov.h>
 #include <gkyl_mom_vlasov.h>
 #include <gkyl_mom_vlasov_sr.h>
 #include <gkyl_null_pool.h>
@@ -205,6 +212,40 @@ struct vm_rad_drag {
   gkyl_dg_updater_rad_vlasov *rad_slvr; // radiation solver
 };
 
+struct vm_fpo_collisions {
+  bool write_diagnostics; // Whether to write diagnostics out.
+  struct gkyl_array *gamma; // FPO Gamma factor
+  struct gkyl_array *h, *g; // Rosenbluth potentials
+  struct gkyl_array *h_host, *g_host;
+
+  // Maxwellian potentials and derivatives on velocity space edges for boundary conditions
+  struct gkyl_array *h_surf, *g_surf;
+  struct gkyl_array *dhdv_surf, *dgdv_surf, *d2gdv2_surf; // dG/dv evaluated at transverse surfaces, dH/dv and d2G/dv2 at same surface
+
+  gkyl_proj_maxwellian_pots_on_basis *pot_slvr; // potential solver for Maxwellian potentials
+
+  struct vm_species_moment lte_moms; // calculator for LTE moments for potential calculation (n u_drift, T/m)
+  struct vm_species_moment moms; // calculator for moments (M0, M1i, M2)
+
+  struct gkyl_array *drag_coeff, *diff_coeff; // Drag and diffusion coefficients
+  struct gkyl_array *drag_coeff_host, *diff_coeff_host; // Drag and diffusion coefficients
+  struct gkyl_array *drag_coeff_surf, *diff_coeff_surf; // Drag and diffusion coefficient surface expansions at lower cell boundary
+  struct gkyl_array *sgn_drag_coeff_surf, *const_sgn_drag_coeff_surf; // Sign of drag coefficient at lower boundary of cell, and boolean for checking if sign(drag_coeff) is constant along boundary
+
+  gkyl_fpo_vlasov_coeff_recovery *coeff_recovery; // Struct for drag and diffusion coeff calculation
+
+  // Momentum and energy conservation corrections
+  struct gkyl_array *fpo_moms, *boundary_corrections; // Primitive moments and boundary corrections
+  struct gkyl_array *drag_diff_coeff_corrs; // Correction quantities added to drag and diffusion coefficients
+
+  struct gkyl_mom_calc *fpo_mom_calc; // FPO volume corrections calculator 
+  struct gkyl_mom_calc_bcorr *bcorr_calc; // FPO boundary corrections calculator
+  gkyl_fpo_coeff_correct *coeff_correct_calc; // FPO drag and diffusion coeff correction calculator
+
+  gkyl_dg_updater_collisions *coll_slvr; // collision solver
+  long offsets[36]; // Array of relative offsets for 3- and 6-cell stencils. Will probably be moved to hyper_dg
+};
+
 struct vm_boundary_fluxes {
   struct gkyl_rect_grid boundary_grid[2*GKYL_MAX_CDIM];
   struct gkyl_array *flux_arr[2*GKYL_MAX_CDIM];
@@ -429,6 +470,9 @@ struct vm_species {
     struct {
       struct vm_bgk_collisions bgk; // BGK collisions object
     };
+    struct {
+      struct vm_fpo_collisions fpo; // FPO collisions object
+    };
   }; 
 
   enum gkyl_radiation_id radiation_id; // type of radiation
@@ -1090,6 +1134,52 @@ void vm_species_bgk_rhs(gkyl_vlasov_app *app,
  */
 void vm_species_bgk_release(const struct gkyl_vlasov_app *app, const struct vm_bgk_collisions *bgk);
 
+/** vm_species_fpo API */
+
+/**
+ * Initialize species fpo collisions object.
+ *
+ * @param app Vlasov app object
+ * @param s Species object 
+ * @param fpo Species fpo object
+ * @param collides_with_fluid Boolean for if kinetic species collides with a fluid species
+ */
+void vm_species_fpo_init(struct gkyl_vlasov_app *app, struct vm_species *s,
+  struct vm_fpo_collisions *fpo);
+
+/**
+ * Compute drag and diffusion coefficients for update
+ *
+ * @param app Vlasov app object
+ * @param species Pointer to species
+ * @param fpo Pointer to fpo
+ * @param fin Input distribution function
+ */
+void vm_species_fpo_drag_diff_coeffs(gkyl_vlasov_app *app, const struct vm_species *s,
+  struct vm_fpo_collisions *fpo, const struct gkyl_array *fin);
+
+/**
+ * Compute RHS from fpo collisions
+ *
+ * @param app Vlasov app object
+ * @param species Pointer to species
+ * @param fpo Pointer to fpo
+ * @param fin Input distribution function
+ * @param rhs On output, the RHS from fpo
+ * @return Maximum stable time-step
+ */
+void vm_species_fpo_rhs(gkyl_vlasov_app *app,
+  const struct vm_species *species,
+  struct vm_fpo_collisions *fpo,
+  const struct gkyl_array *fin, struct gkyl_array *rhs);
+
+/**
+ * Release species fpo object.
+ *
+ * @param app Vlasov app object
+ * @param sm Species fpo object to release
+ */
+void vm_species_fpo_release(const struct gkyl_vlasov_app *app, const struct vm_fpo_collisions *fpo);
 /** vm_species_radiation API */
 
 /**

vlasov/apps/vlasov.c

@@ -610,6 +610,9 @@ gkyl_vlasov_app_write(gkyl_vlasov_app* app, double tm, int frame)
     if (app->species[i].info.output_f_lte) {
       gkyl_vlasov_app_write_species_lte(app, i, tm, frame);
     }
+    if ((app->species[i].collision_id == GKYL_FPO_COLLISIONS) && (app->species[i].fpo.write_diagnostics)) {
+      gkyl_vlasov_app_write_species_fpo(app, i, tm, frame);
+    }
   }
   for (int i=0; i<app->num_fluid_species; ++i) {
     gkyl_vlasov_app_write_fluid_species(app, i, tm, frame);
@@ -781,6 +784,59 @@ gkyl_vlasov_app_write_species_lte(gkyl_vlasov_app* app, int sidx, double tm, int
   vlasov_array_meta_release(mt); 
 }
 
+void
+gkyl_vlasov_app_write_species_fpo(gkyl_vlasov_app* app, int sidx, double tm, int frame)
+{
+  struct gkyl_msgpack_data *mt = vlasov_array_meta_new( (struct vlasov_output_meta) {
+      .frame = frame,
+      .stime = tm,
+      .poly_order = app->poly_order,
+      .basis_type = app->basis.id
+    }
+  );  
+
+  struct vm_species *vm_s = &app->species[sidx];
+
+  const char *fmt_h = "%s-%s_H_%d.gkyl";
+  const char *fmt_g = "%s-%s_G_%d.gkyl";
+  const char *fmt_drag = "%s-%s_drag_coeff_%d.gkyl";
+  const char *fmt_diff = "%s-%s_diff_coeff_%d.gkyl";
+  int sz = gkyl_calc_strlen(fmt_drag, app->name,vm_s->info.name, frame);
+  char fileNm[sz+1]; // ensures no buffer overflow
+
+  vm_species_fpo_drag_diff_coeffs(app, vm_s, &vm_s->fpo, vm_s->f);
+
+  if (app->use_gpu) {
+    // copy data from device to host before writing it out
+    gkyl_array_copy(vm_s->fpo.h_host, vm_s->fpo.h);
+    gkyl_array_copy(vm_s->fpo.g_host, vm_s->fpo.g);
+    gkyl_array_copy(vm_s->fpo.drag_coeff_host, vm_s->fpo.drag_coeff);
+    gkyl_array_copy(vm_s->fpo.diff_coeff_host, vm_s->fpo.diff_coeff);
+  }
+
+  // Write H
+  snprintf(fileNm, sizeof fileNm, fmt_h, app->name,vm_s->info.name, frame);
+  gkyl_comm_array_write(vm_s->comm, &vm_s->grid, &vm_s->local, 
+    mt, vm_s->fpo.h_host, fileNm);
+
+  // Write G
+  snprintf(fileNm, sizeof fileNm, fmt_g, app->name,vm_s->info.name, frame);
+  gkyl_comm_array_write(vm_s->comm, &vm_s->grid, &vm_s->local, 
+    mt, vm_s->fpo.g_host, fileNm);
+
+  // Write drag coefficient
+  snprintf(fileNm, sizeof fileNm, fmt_drag, app->name,vm_s->info.name, frame);
+  gkyl_comm_array_write(vm_s->comm, &vm_s->grid, &vm_s->local, 
+    mt, vm_s->fpo.drag_coeff_host, fileNm);
+
+  // Write diffusion tensor
+  snprintf(fileNm, sizeof fileNm, fmt_diff, app->name,vm_s->info.name, frame);
+  gkyl_comm_array_write(vm_s->comm, &vm_s->grid, &vm_s->local, 
+    mt, vm_s->fpo.diff_coeff_host, fileNm);
+
+  vlasov_array_meta_release(mt);     
+}
+
 void
 gkyl_vlasov_app_write_fluid_species(gkyl_vlasov_app* app, int sidx, double tm, int frame)
 {

vlasov/apps/vlasov_forward_euler.c

@@ -46,6 +46,10 @@ vlasov_forward_euler(gkyl_vlasov_app* app, double tcurr, double dt,
       vm_species_bgk_moms(app, &app->species[i], 
         &app->species[i].bgk, fin[i]);
     }
+    else if (app->species[i].collision_id == GKYL_FPO_COLLISIONS) {
+      vm_species_fpo_drag_diff_coeffs(app, &app->species[i], 
+        &app->species[i].fpo, fin[i]);
+    }
   }
 
   // compute necessary moments for cross-species collisions

vlasov/apps/vm_species.c

@@ -225,7 +225,7 @@ vm_species_init(struct gkyl_vm *vm, struct gkyl_vlasov_app *app, struct vm_speci
     s->sgn_alpha_surf = mkarr(app->use_gpu, sgn_alpha_surf_sz, s->local_ext.volume);
     s->const_sgn_alpha = mk_int_arr(app->use_gpu, (cdim + vdim), s->local_ext.volume);
 
-    // Pre-compute alpha_surf, sgn_alpha_surf, const_sgn_alpha, and cot_vec since they are time-independent
+    // Pre-compute alpha_surf, sgn_alpha_surf, and const_sgn_alpha since they are time-independent
     struct gkyl_dg_calc_canonical_pb_vars *calc_vars = gkyl_dg_calc_canonical_pb_vars_new(&s->grid, 
       &app->confBasis, &app->basis, app->use_gpu);
     gkyl_dg_calc_canonical_pb_vars_alpha_surf(calc_vars, &app->local, &s->local, &s->local_ext, s->hamil,
@@ -241,7 +241,7 @@ vm_species_init(struct gkyl_vm *vm, struct gkyl_vlasov_app *app, struct vm_speci
   }
   else {
     if (s->field_id == GKYL_FIELD_NULL || s->field_id == GKYL_FIELD_E_B) {
-      struct gkyl_dg_vlasov_auxfields aux_inp = {.field = s->qmem, .cot_vec = 0, 
+      struct gkyl_dg_vlasov_auxfields aux_inp = {.field = s->qmem, 
         .alpha_surf = 0, .sgn_alpha_surf = 0, .const_sgn_alpha = 0 };
       s->slvr = gkyl_dg_updater_vlasov_new(&s->grid, &app->confBasis, &app->basis, 
         &app->local, &s->local_vel, &s->local, is_zero_flux, s->model_id, s->field_id, &aux_inp, app->use_gpu);
@@ -327,6 +327,7 @@ vm_species_init(struct gkyl_vm *vm, struct gkyl_vlasov_app *app, struct vm_speci
   // initialize empty collision structs so inputs of structs are set to 0
   s->lbo = (struct vm_lbo_collisions) { };
   s->bgk = (struct vm_bgk_collisions) { };
+  s->fpo = (struct vm_fpo_collisions) { };
   if (s->info.output_f_lte){
     // Always have correct moments on for the f_lte output
     struct correct_all_moms_inp corr_inp = { .correct_all_moms = true, 
@@ -340,6 +341,9 @@ vm_species_init(struct gkyl_vm *vm, struct gkyl_vlasov_app *app, struct vm_speci
   else if (s->collision_id == GKYL_BGK_COLLISIONS) {
     vm_species_bgk_init(app, s, &s->bgk);
   }
+  else if (s->collision_id == GKYL_FPO_COLLISIONS) {
+    vm_species_fpo_init(app, s, &s->fpo);
+  }
 
   // determine radiation type to use in vlasov update
   s->radiation_id = s->info.radiation.radiation_id;
@@ -515,6 +519,9 @@ vm_species_rhs(gkyl_vlasov_app *app, struct vm_species *species,
     species->bgk.implicit_step = false;
     vm_species_bgk_rhs(app, species, &species->bgk, fin, rhs);
   }
+  else if (species->collision_id == GKYL_FPO_COLLISIONS) {
+    vm_species_fpo_rhs(app, species, &species->fpo, fin, rhs);
+  }
 
   if (species->calc_bflux) {
     vm_species_bflux_rhs(app, species, &species->bflux, fin, rhs);
@@ -810,6 +817,9 @@ vm_species_release(const gkyl_vlasov_app* app, const struct vm_species *s)
   else if (s->collision_id == GKYL_BGK_COLLISIONS) {
     vm_species_bgk_release(app, &s->bgk);
   }
+  else if (s->collision_id == GKYL_FPO_COLLISIONS) {
+    vm_species_fpo_release(app, &s->fpo);
+  }
 
   if (s->radiation_id == GKYL_VM_COMPTON_RADIATION) {
     vm_species_radiation_release(app, &s->rad);