Get EM back 🚧

core/zero/gkyl_eqn_type.h

@@ -34,11 +34,14 @@ enum gkyl_eqn_type {
 
 // Identifiers for specific gyrokinetic field object types.
 enum gkyl_gkfield_id {
-  GKYL_GK_FIELD_ES = 0, // Electrostatic GK. This is default.
+  GKYL_GK_FIELD = 0, // This is default.
   GKYL_GK_FIELD_BOLTZMANN = 1, // GK Boltzmann, isothermal electrons, phi = phi_sheath + (T_e/e)*ln(n_i/n_is).
   GKYL_GK_FIELD_ADIABATIC = 2, // GK field with an adiabatic species.
-  GKYL_GK_FIELD_ES_IWL = 3, // Inner-wall limited ES.
-  GKYL_GK_FIELD_EM = 4, // Electromagnetic GK.
+  GKYL_GK_FIELD_IWL = 3, // Inner-wall limited.
+  GKYL_GK_FIELD_ES = 4, // Electrostatic field.
+  GKYL_GK_FIELD_ES_IWL = 5, // Electromagnetic field with inner-wall limited BCs.
+  GKYL_GK_FIELD_EM = 6, // To force electromagnetic field solve (even if no species is EM).
+  GKYL_GK_FIELD_EM_IWL = 7, // As above but for IWL.
 };
 
 // Identifiers for specific field object types.
@@ -177,7 +180,6 @@ enum gkyl_gk_anomalous_diff_id {
 enum gkyl_gk_collisionless_type {
   GKYL_GK_COLLISIONLESS_NONE = 0, // Turns off collisionless terms.
   GKYL_GK_COLLISIONLESS_ES, // Electrostatic.
-  GKYL_GK_COLLISIONLESS_ES_NO_BY, // Electrostatic, w/o toroidal field (b_y=0).
   GKYL_GK_COLLISIONLESS_EM_BPERP, // Electromagnetic with B_perp fluctuations.
   GKYL_GK_COLLISIONLESS_EM_BPAR, // Electromagnetic with B_par fluctuations.
   GKYL_GK_COLLISIONLESS_EM, // Electromagnetic with B_perp and B_par fluctuations.

gyrokinetic/apps/gk_field.c

@@ -68,9 +68,19 @@ gk_field_calc_energy_dt_active(gkyl_gyrokinetic_app *app, const struct gk_field
   app->stat.phidot_tm += gkyl_time_diff_now_sec(wst);
 }
 
+static void
+gk_field_calc_apar_energy_dt_active(gkyl_gyrokinetic_app *app, const struct gk_field *field, double dt, double *energy_reduced)
+{
+  struct timespec wst = gkyl_wall_clock();
+  gkyl_array_integrate_advance(field->calc_em_energy, field->apar, 
+    1.0/dt, field->apar_energy_fac, &app->local, &app->local, energy_reduced);
+  app->stat.phidot_tm += gkyl_time_diff_now_sec(wst);
+}
+
 static void
 gk_field_calc_energy_dt_none(gkyl_gyrokinetic_app *app, const struct gk_field *field, double dt, double *energy_reduced)
 {
+  // Do nothing.
 }
 
 static void
@@ -122,6 +132,67 @@ gk_field_calc_energy_enabled(struct gkyl_gyrokinetic_app *app, const struct gk_f
 
     gkyl_dynvec_append(field->integ_energy_dot, tm, energy_dot_global);
   }
+
+  if (field->is_em) {
+    // Calculate Aparallel energy.
+    gkyl_array_integrate_advance(field->calc_em_energy, field->apar, 
+      1.0, field->apar_energy_fac, &app->local, &app->local, field->em_energy_red);
+
+    gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_SUM, 1, field->em_energy_red, field->em_energy_red_global);
+
+    double energy_global[1] = { 0.0 };
+    if (app->use_gpu)
+      gkyl_cu_memcpy(energy_global, field->em_energy_red_global, sizeof(double[1]), GKYL_CU_MEMCPY_D2H);
+    else
+      energy_global[0] = field->em_energy_red_global[0];
+
+    if (app->cdim == 1)
+      energy_global[0] *= field->apar_energy_fac_1d; 
+
+    gkyl_dynvec_append(field->integ_apar_energy, tm, energy_global);
+
+    // Calculate d(Aparallel)/dt energy.
+    gkyl_array_integrate_advance(field->calc_em_energy, field->apardot, 
+      1.0, field->apar_energy_fac, &app->local, &app->local, field->em_energy_red);
+
+    gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_SUM, 1, field->em_energy_red, field->em_energy_red_global);
+
+    double energy_dot_global[1] = { 0.0 };
+    if (app->use_gpu)
+      gkyl_cu_memcpy(energy_dot_global, field->em_energy_red_global, sizeof(double[1]), GKYL_CU_MEMCPY_D2H);
+    else
+      energy_dot_global[0] = field->em_energy_red_global[0];
+
+    // if (app->cdim == 1)
+    //   energy_dot_global[0] *= field->apar_energy_fac_1d;
+
+    gkyl_dynvec_append(field->integ_apardot_energy, tm, energy_dot_global);
+
+    if (field->info.time_rate_diagnostics) {
+      gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_SUM, 1, field->apar_energy_red_old, field->em_energy_red_global);
+      double energy_dot_global_old[1] = { 0.0 };
+      if (app->use_gpu)
+        gkyl_cu_memcpy(energy_dot_global_old, field->em_energy_red_global, sizeof(double[1]), GKYL_CU_MEMCPY_D2H);
+      else
+        energy_dot_global_old[0] = field->em_energy_red_global[0];
+      if (app->cdim == 1)
+        energy_dot_global_old[0] *= field->apar_energy_fac_1d; 
+
+      gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_SUM, 1, field->apar_energy_red_new, field->em_energy_red_global);
+      double energy_dot_global_new[1] = { 0.0 };
+      if (app->use_gpu)
+        gkyl_cu_memcpy(energy_dot_global_new, field->em_energy_red_global, sizeof(double[1]), GKYL_CU_MEMCPY_D2H);
+      else
+        energy_dot_global_new[0] = field->em_energy_red_global[0];
+      if (app->cdim == 1)
+        energy_dot_global_new[0] *= field->apar_energy_fac_1d; 
+
+      double energy_dot_global[1] = { 0.0 };
+      energy_dot_global[0] = energy_dot_global_new[0] - energy_dot_global_old[0];
+
+      gkyl_dynvec_append(field->integ_apar_energy_dot, tm, energy_dot_global);
+    }
+  }
 }
 
 static void
@@ -130,6 +201,52 @@ gk_field_calc_energy_disabled(struct gkyl_gyrokinetic_app *app, const struct gk_
   // Do nothing.
 }
 
+static void
+gk_field_combine_enabled(struct gkyl_array *out, double c1,
+  const struct gkyl_array *arr1, double c2, const struct gkyl_array *arr2,
+  const struct gkyl_range *rng)
+{
+  gkyl_array_accumulate_range(gkyl_array_set_range(out, c1, arr1, rng),
+    c2, arr2, rng);
+}
+
+static void
+gk_field_combine_disabled(struct gkyl_array *out, double c1,
+  const struct gkyl_array *arr1, double c2, const struct gkyl_array *arr2,
+  const struct gkyl_range *rng)
+{
+  // Do nothing.
+}
+
+static void
+gk_field_step_apar_enabled(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array* out, double dt,
+  const struct gkyl_array* inp)
+{
+  // Apar^{n+1} = Apar^{n} + dt*dApar/dt
+  gkyl_array_accumulate(gkyl_array_scale(out, dt), 1.0, inp);
+}
+
+static void
+gk_field_step_apar_disabled(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array* out, double dt,
+  const struct gkyl_array* inp)
+{
+  // Do nothing.
+}
+
+static void
+gk_field_copy_range_enabled(struct gkyl_array *out,
+  const struct gkyl_array *inp, const struct gkyl_range *range)
+{
+  gkyl_array_copy_range(out, inp, range);
+}
+
+static void 
+gk_field_copy_range_disabled(struct gkyl_array *out,
+  const struct gkyl_array *inp, const struct gkyl_range *range)
+{
+  // Do nothing.
+}
+
 static void
 gk_field_time_rate_diags_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f)
 {
@@ -147,6 +264,22 @@ gk_field_time_rate_diags_new(struct gkyl_gyrokinetic_app *app, struct gk_field *
   }
   f->integ_energy_dot = gkyl_dynvec_new(GKYL_DOUBLE, 1);
   f->is_first_energy_dot_write_call = true;
+  if (f->is_em){
+    f->calc_apar_energy_dt_func = gk_field_calc_apar_energy_dt_active;
+    if (app->use_gpu) {
+      f->apar_energy_red_new = gkyl_cu_malloc(sizeof(double[1]));
+      f->apar_energy_red_old = gkyl_cu_malloc(sizeof(double[1]));
+      gkyl_cu_memset(f->apar_energy_red_new, 0, sizeof(double[1]));
+      gkyl_cu_memset(f->apar_energy_red_old, 0, sizeof(double[1]));
+    } else {
+      f->apar_energy_red_new = gkyl_malloc(sizeof(double[1]));
+      f->apar_energy_red_old = gkyl_malloc(sizeof(double[1]));
+      memset(f->apar_energy_red_new, 0, sizeof(double[1]));
+      memset(f->apar_energy_red_old, 0, sizeof(double[1]));
+    }
+    f->integ_apar_energy_dot = gkyl_dynvec_new(GKYL_DOUBLE, 1);
+  }
+  f->is_first_apar_energy_dot_write_call = true;
 }
 
 static void
@@ -165,15 +298,22 @@ gk_field_energy_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f)
   f->is_first_energy_write_call = true;
 
   f->calc_energy_func = gk_field_calc_energy_enabled;
+
   f->calc_energy_dt_func = gk_field_calc_energy_dt_none;
-
-  if (f->info.time_rate_diagnostics) {
+  f->calc_apar_energy_dt_func = gk_field_calc_energy_dt_none;
+  if (f->info.time_rate_diagnostics)
     gk_field_time_rate_diags_new(app, f);
-  }
 
   // Factors for ES energy.
   f->es_energy_fac = mkarr(app->use_gpu, (2*(app->cdim/3)+1)*app->basis.num_basis, app->local_ext.volume);
   f->es_energy_fac_1d = 0.0;
+
+  if (f->is_em) {
+    f->integ_apar_energy = gkyl_dynvec_new(GKYL_DOUBLE, 1);
+    f->integ_apardot_energy = gkyl_dynvec_new(GKYL_DOUBLE, 1);
+    f->apar_energy_fac = mkarr(app->use_gpu, (2*(app->cdim/3)+1)*app->basis.num_basis, app->local_ext.volume);
+    f->apar_energy_fac_1d = 0;
+  }
 }
 
 static void
@@ -188,6 +328,16 @@ gk_field_time_rate_diags_release(const struct gkyl_gyrokinetic_app *app, struct
     gkyl_free(f->em_energy_red_old);
   }
   gkyl_dynvec_release(f->integ_energy_dot);
+  if (f->is_em) {
+    if (app->use_gpu) {
+      gkyl_cu_free(f->apar_energy_red_new);
+      gkyl_cu_free(f->apar_energy_red_old);
+    } else {
+      gkyl_free(f->apar_energy_red_new);
+      gkyl_free(f->apar_energy_red_old);
+    }
+    gkyl_dynvec_release(f->integ_apar_energy_dot);
+  }
 }
 
 static void
@@ -208,6 +358,12 @@ gk_field_energy_release(const struct gkyl_gyrokinetic_app *app, struct gk_field
   }
 
   gkyl_array_release(f->es_energy_fac);
+
+  if (f->is_em) {
+    gkyl_dynvec_release(f->integ_apar_energy);
+    gkyl_dynvec_release(f->integ_apardot_energy);
+    gkyl_array_release(f->apar_energy_fac);
+  }
 }
 
 // Related to enforcing parallel Vlasov boundary conditions
@@ -226,6 +382,25 @@ gk_field_new(struct gkyl_gk *gk, struct gkyl_gyrokinetic_app *app)
   f->info = gk->field;
 
   f->gkfield_id = f->info.gkfield_id ? f->info.gkfield_id : GKYL_GK_FIELD_ES;
+
+  f->is_em = f->info.gkfield_id == GKYL_GK_FIELD_EM || f->info.gkfield_id == GKYL_GK_FIELD_EM_IWL;
+  // Ensure that if any species is electromagnetic, the field is electromagnetic.
+  for (int i=0; i<app->num_species; ++i) {
+    struct gk_species *s = &app->species[i];
+    f->is_em = f->is_em || (s->info.collisionless.type == GKYL_GK_COLLISIONLESS_EM) || 
+                      (s->info.collisionless.type == GKYL_GK_COLLISIONLESS_EM_BPERP) ||
+                      (s->info.collisionless.type == GKYL_GK_COLLISIONLESS_EM_BPAR);
+  }
+  if (f->is_em) {
+    assert(f->info.mu0 > 0.0);
+    f->step_apar = gk_field_step_apar_enabled;
+    f->em_combine_func = gk_field_combine_enabled;
+    f->em_copy_func = gk_field_copy_range_enabled;
+  } else {
+    f->step_apar = gk_field_step_apar_disabled;
+    f->em_combine_func = gk_field_combine_disabled;
+    f->em_copy_func = gk_field_copy_range_disabled;
+  }
 
   f->calc_init_field = !f->info.zero_init_field;
   f->update_field = !f->info.is_static;
@@ -269,6 +444,46 @@ gk_field_rhs(gkyl_gyrokinetic_app *app, struct gk_field *field)
   app->stat.field_phi_solve_tm += gkyl_time_diff_now_sec(wst);
 }
 
+void
+gk_field_em_rhs(gkyl_gyrokinetic_app *app, struct gk_field *field, const struct gkyl_array *f_in[],  struct gkyl_array *rhs_in[])
+{
+  // Compute Apardot through Ohm's law.
+  struct timespec wst = gkyl_wall_clock();
+  field->em_rhs_func(app, field, f_in, rhs_in);
+  app->stat.field_tm += gkyl_time_diff_now_sec(wst);
+}
+
+void gk_field_calc_apar_ic(gkyl_gyrokinetic_app *app, struct gk_field *field)
+{
+  struct timespec wst = gkyl_wall_clock();
+  field->ampere_solve(app, field);
+  // Smooth Apar after solving Ampere's law.
+  gk_field_fem_projection_par(app, field, field->apar, field->apar);
+  app->stat.field_apar_solve_tm += gkyl_time_diff_now_sec(wst);
+}
+
+void gk_field_step_apar(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array* out, double a, const struct gkyl_array* inp)
+{
+  struct timespec wst = gkyl_wall_clock();
+  field->step_apar(app, field, out, a, inp);
+  app->stat.field_apar_solve_tm += gkyl_time_diff_now_sec(wst);
+}
+
+void
+gk_field_em_combine(struct gk_field *field, struct gkyl_array *out, double c1,
+  const struct gkyl_array *arr1, double c2, const struct gkyl_array *arr2,
+  const struct gkyl_range *rng)
+{
+  field->em_combine_func(out, c1, arr1, c2, arr2, rng);
+}
+
+void
+gk_field_em_copy_range(struct gk_field *field, struct gkyl_array *out,
+  const struct gkyl_array *inp, const struct gkyl_range *range)
+{
+  field->em_copy_func(out, inp, range);
+}
+
 void
 gk_field_calc_energy(gkyl_gyrokinetic_app *app, double tm, const struct gk_field *field)
 {
@@ -283,6 +498,12 @@ gk_field_calc_energy_dt(gkyl_gyrokinetic_app *app, const struct gk_field *field,
   field->calc_energy_dt_func(app, field, dt, energy_reduced);
 }
 
+void
+gk_field_calc_apar_energy_dt(gkyl_gyrokinetic_app *app, const struct gk_field *field, double dt, double *energy_reduced)
+{
+  field->calc_apar_energy_dt_func(app, field, dt, energy_reduced);
+}
+
 void gk_field_accumulate_rho_c_adiabatic(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gk_species *s, struct gkyl_array **bflux)
 {
   // Gyroaverage the density if needed.
@@ -316,6 +537,54 @@ gk_field_accumulate_rho_c(gkyl_gyrokinetic_app *app, struct gk_field *field,
   app->stat.field_phi_rhs_tm += gkyl_time_diff_now_sec(wst);
 }
 
+void
+gk_field_accumulate_current_dens(gkyl_gyrokinetic_app *app, struct gk_field *field, 
+  const struct gkyl_array *fin[])
+{
+  struct timespec wst = gkyl_wall_clock();
+  gkyl_array_clear(field->currentDens, 0.0);
+  for (int i=0; i<app->num_species; ++i) {
+    struct gk_species *s = &app->species[i];
+    gk_species_moment_calc(&s->m1, s->local, app->local, fin[i]);
+    gkyl_array_accumulate_range(field->currentDens, s->info.charge, s->m1.marr, &app->local);
+  } 
+  app->stat.field_apar_rhs_tm += gkyl_time_diff_now_sec(wst);
+}
+
+void
+gk_field_accumulate_current_dens_dot(gkyl_gyrokinetic_app *app, struct gk_field *field, 
+  struct gkyl_array *rhs_in[])
+{
+  struct timespec wst = gkyl_wall_clock();
+  gkyl_array_clear(field->currentDensdot, 0.0);
+  for (int i=0; i<app->num_species; ++i) {
+    struct gk_species *s = &app->species[i];
+
+    gk_species_moment_calc(&s->m1, s->local, app->local, rhs_in[i]);
+
+    gkyl_array_accumulate_range(field->currentDensdot, s->info.charge, s->m1.marr, &app->local);
+  } 
+  app->stat.field_apar_rhs_tm += gkyl_time_diff_now_sec(wst);
+}
+
+void
+gk_field_accumulate_ohms_kSq(gkyl_gyrokinetic_app *app, struct gk_field *field, 
+  const struct gkyl_array *fin[])
+{
+  struct timespec wst = gkyl_wall_clock();
+  gkyl_array_clear(field->dApartdtSlvr_kSq, 0.0);
+  for (int i=0; i<app->num_species; ++i) {
+    struct gk_species *s = &app->species[i];
+    gk_species_moment_calc(&s->m0, s->local, app->local, fin[i]);
+    // Gyroaverage the density if needed.
+    s->gyroaverage(app, s, s->m0.marr, s->m0_gyroavg);
+    // Use m0 to update also update the kSq matrix for the Ohm's law solver.
+    double fac = -s->info.charge*s->info.charge/s->info.mass;
+    gkyl_array_accumulate_range(field->dApartdtSlvr_kSq, fac, s->m0_gyroavg, &app->local);
+  }
+  app->stat.field_apar_rhs_tm += gkyl_time_diff_now_sec(wst);
+}
+
 void
 gk_field_fem_projection_par(gkyl_gyrokinetic_app *app, struct gk_field *field, struct gkyl_array *arr_dg, struct gkyl_array *arr_fem)
 {