Refactor phi wall boundary conditions into species table

gyrokinetic/apps/gk_field.c

@@ -306,9 +306,6 @@ gk_field_new(struct gkyl_gk *gk, struct gkyl_gyrokinetic_app *app)
     }
   }
 
-  // Initialize biased walls.
-  gk_field_biased_wall_new(app, f);
-
   return f;
 }
 
@@ -432,8 +429,5 @@ gk_field_release(const gkyl_gyrokinetic_app* app, struct gk_field *f)
   // Release energy diagnostics.
   gk_field_energy_release(app, f);
 
-  // Release biased walls.
-  gk_field_biased_wall_release(app, f);
-
   gkyl_free(f);
 }

gyrokinetic/apps/gk_species.c

@@ -197,7 +197,8 @@ gk_species_rhs_implicit_static(gkyl_gyrokinetic_app *app, struct gk_species *spe
 }
 
 static void
-gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f)
+gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *species,
+  double tm, struct gkyl_array *f)
 {
   struct timespec wst = gkyl_wall_clock();
 
@@ -210,8 +211,9 @@ gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *
 
       switch (species->lower_bc[d].type) {
         case GKYL_BC_GK_SPECIES_SHEATH:
+          gk_species_phi_wall_advance(app, &species->phi_wall_lo, tm);
           gkyl_bc_sheath_gyrokinetic_advance(species->bc_sheath_lo, app->field->phi_smooth, 
-            app->field->phi_wall_lo, f, &app->local);
+            species->phi_wall_lo.phi, f, &app->local);
           break;
         case GKYL_BC_GK_SPECIES_TWISTSHIFT:
           gkyl_bc_twistshift_advance(species->bc_ts_lo, f, f);
@@ -232,8 +234,9 @@ gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *
 
       switch (species->upper_bc[d].type) {
         case GKYL_BC_GK_SPECIES_SHEATH:
+          gk_species_phi_wall_advance(app, &species->phi_wall_up, tm);
           gkyl_bc_sheath_gyrokinetic_advance(species->bc_sheath_up, app->field->phi_smooth, 
-            app->field->phi_wall_up, f, &app->local);
+            species->phi_wall_up.phi, f, &app->local);
           break;
         case GKYL_BC_GK_SPECIES_TWISTSHIFT:
           gkyl_bc_twistshift_advance(species->bc_ts_up, f, f);
@@ -266,7 +269,8 @@ gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *
 }
 
 static void
-gk_species_apply_bc_static(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f)
+gk_species_apply_bc_static(gkyl_gyrokinetic_app *app, const struct gk_species *species,
+  double tm, struct gkyl_array *f)
 {
   // do nothing
 }
@@ -677,6 +681,7 @@ gk_species_release_dynamic(const gkyl_gyrokinetic_app* app, const struct gk_spec
   for (int d=0; d<app->cdim; ++d) {
     if (s->lower_bc[d].type == GKYL_BC_GK_SPECIES_SHEATH) {
       gkyl_bc_sheath_gyrokinetic_release(s->bc_sheath_lo);
+      gk_species_phi_wall_release(app, &s->phi_wall_lo);
     }
     else if (s->lower_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) { 
       gkyl_bc_twistshift_release(s->bc_ts_lo);
@@ -690,6 +695,7 @@ gk_species_release_dynamic(const gkyl_gyrokinetic_app* app, const struct gk_spec
 
     if (s->upper_bc[d].type == GKYL_BC_GK_SPECIES_SHEATH) {
       gkyl_bc_sheath_gyrokinetic_release(s->bc_sheath_up);
+      gk_species_phi_wall_release(app, &s->phi_wall_up);
     }
     else if (s->upper_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) {
       gkyl_bc_twistshift_release(s->bc_ts_up);
@@ -877,6 +883,7 @@ gk_species_init_dynamic(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app
       struct gkyl_range *sol_ghost = gk_app_inp->geometry.has_LCFS? &gks->local_lower_ghost_par_sol : &gks->local_lower_ghost[d];
       gks->bc_sheath_lo = gkyl_bc_sheath_gyrokinetic_new(d, GKYL_LOWER_EDGE, gks->basis_on_dev, 
         sol_skin, sol_ghost, gks->vel_map, cdim, 2.0*(gks->info.charge/gks->info.mass), app->use_gpu);
+      gk_species_phi_wall_init(app, &gks->lower_bc[d], &gks->phi_wall_lo);
     }
     else if (gks->lower_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) {
       assert(cdim == 3);
@@ -937,6 +944,7 @@ gk_species_init_dynamic(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app
       struct gkyl_range *sol_ghost = gk_app_inp->geometry.has_LCFS? &gks->local_upper_ghost_par_sol : &gks->local_upper_ghost[d];
       gks->bc_sheath_up = gkyl_bc_sheath_gyrokinetic_new(d, GKYL_UPPER_EDGE, gks->basis_on_dev, 
         sol_skin, sol_ghost, gks->vel_map, cdim, 2.0*(gks->info.charge/gks->info.mass), app->use_gpu);
+      gk_species_phi_wall_init(app, &gks->upper_bc[d], &gks->phi_wall_up);
     }
     else if (gks->upper_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) {
       assert(cdim == 3);
@@ -1497,7 +1505,7 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st
       gks->upper_bc[d].type = GKYL_BC_GK_SPECIES_PERIODIC;
     }
   }
- 
+
   // Species properties metadata.
   struct gkyl_msgpack_map_elem io_meta_sprop[] = {
     { .key = "mass", .elem_type = GKYL_MP_DOUBLE, .dval = gks->info.mass },
@@ -1935,9 +1943,10 @@ gk_species_copy_range(struct gk_species *species, struct gkyl_array *out,
 }
 
 void
-gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f)
+gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species,
+  double tm, struct gkyl_array *f)
 {
-  species->bc_func(app, species, f);
+  species->bc_func(app, species, tm, f);
 }
 
 void

gyrokinetic/apps/gk_species_fdot_multiplier.c

@@ -144,11 +144,11 @@ gk_species_fdot_multiplier_advance_loss_cone_mult(gkyl_gyrokinetic_app *app,
   int zdim = app->cdim - 1;
   if (gks->lower_bc[zdim].type == GKYL_BC_GK_SPECIES_SHEATH) {
     gkyl_comm_array_allgather(app->comm, &app->local, &app->global,
-      app->field->phi_wall_lo, fdmul->phi_wall_lo_global);
+      gks->phi_wall_lo.phi, fdmul->phi_wall_lo_global);
   }
   if (gks->upper_bc[zdim].type == GKYL_BC_GK_SPECIES_SHEATH) {
     gkyl_comm_array_allgather(app->comm, &app->local, &app->global,
-      app->field->phi_wall_up, fdmul->phi_wall_up_global);
+      gks->phi_wall_up.phi, fdmul->phi_wall_up_global);
   }
   gkyl_loss_cone_mask_gyrokinetic_advance(fdmul->lcm_proj_op, &gks->local, &app->global,
     fdmul->bmag_global, fdmul->phi_global, fdmul->phi_wall_lo_global,
@@ -338,6 +338,8 @@ gk_species_fdot_multiplier_init_comp(gkyl_gyrokinetic_app *app, struct gk_specie
     }
     else if (fdmul->type == GKYL_GK_FDOT_MULTIPLIER_LOSS_CONE) {
       // Operator that projects the loss cone mask.
+      fdmul->phi_wall_lo_global = 0;
+      fdmul->phi_wall_up_global = 0;
       int zdim = app->cdim - 1;
       struct gkyl_loss_cone_mask_gyrokinetic_inp inp_proj = {
         .conf_basis = &app->basis,
@@ -354,10 +356,12 @@ gk_species_fdot_multiplier_init_comp(gkyl_gyrokinetic_app *app, struct gk_specie
       fdmul->bmag_global = mkarr(app->use_gpu, app->gk_geom->geo_corn.bmag->ncomp,
         app->global_ext.volume);
       fdmul->phi_global = mkarr(app->use_gpu, app->basis.num_basis, app->global_ext.volume);
-      fdmul->phi_wall_lo_global = mkarr(app->use_gpu, app->basis.num_basis,
-        app->global_ext.volume);
-      fdmul->phi_wall_up_global = mkarr(app->use_gpu, app->basis.num_basis,
-        app->global_ext.volume);
+      if (inp_proj.lower_boundary == GKYL_LOSS_CONE_BC_SHEATH)
+        fdmul->phi_wall_lo_global = mkarr(app->use_gpu, app->basis.num_basis,
+          app->global_ext.volume);
+      if (inp_proj.upper_boundary == GKYL_LOSS_CONE_BC_SHEATH)
+        fdmul->phi_wall_up_global = mkarr(app->use_gpu, app->basis.num_basis,
+          app->global_ext.volume);
 
       gkyl_comm_array_allgather(app->comm, &app->local, &app->global, app->gk_geom->geo_corn.bmag,
         fdmul->bmag_global);
@@ -540,8 +544,10 @@ gk_species_fdot_multiplier_release_comp(const struct gkyl_gyrokinetic_app *app,
     gkyl_array_release(fdmul->buffer);
     gkyl_array_release(fdmul->bmag_global);
     gkyl_array_release(fdmul->phi_global);
-    gkyl_array_release(fdmul->phi_wall_lo_global);
-    gkyl_array_release(fdmul->phi_wall_up_global);
+    if (fdmul->phi_wall_lo_global)
+      gkyl_array_release(fdmul->phi_wall_lo_global);
+    if (fdmul->phi_wall_up_global)
+      gkyl_array_release(fdmul->phi_wall_up_global);
     gkyl_loss_cone_mask_gyrokinetic_release(fdmul->lcm_proj_op);
   }
   else if (fdmul->type == GKYL_GK_FDOT_MULTIPLIER_CONSTANT) {

gyrokinetic/apps/gk_species_phi_wall.c

@@ -0,0 +1,64 @@
+#include <gkyl_array_ops.h>
+#include <gkyl_gyrokinetic_priv.h>
+
+static void
+gk_species_phi_wall_advance_disabled(gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall, double tm)
+{
+}
+
+static void
+gk_species_phi_wall_advance_enabled(gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall, double tm)
+{
+  gkyl_eval_on_nodes_advance(wall->projector, tm, &app->local_ext, wall->phi_host);
+  if (app->use_gpu)
+    gkyl_array_copy(wall->phi, wall->phi_host);
+}
+
+void
+gk_species_phi_wall_init(gkyl_gyrokinetic_app *app,
+  const struct gkyl_gyrokinetic_bc *bc, struct gk_species_wall_potential *wall)
+{
+  *wall = (struct gk_species_wall_potential) {
+    .advance_func = gk_species_phi_wall_advance_disabled,
+  };
+  if (bc->type != GKYL_BC_GK_SPECIES_SHEATH)
+    return;
+
+  wall->phi = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+  gkyl_array_clear(wall->phi, 0.0);
+
+  wall->phi_host = wall->phi;
+  if (bc->aux_profile) {
+    if (app->use_gpu)
+      wall->phi_host = mkarr(false, wall->phi->ncomp, wall->phi->size);
+
+    wall->projector = gkyl_eval_on_nodes_new(&app->grid, &app->basis,
+      1, bc->aux_profile, bc->aux_ctx);
+    wall->advance_func = gk_species_phi_wall_advance_enabled;
+    gk_species_phi_wall_advance(app, wall, 0.0);
+  }
+}
+
+void
+gk_species_phi_wall_advance(gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall, double tm)
+{
+  wall->advance_func(app, wall, tm);
+}
+
+void
+gk_species_phi_wall_release(const gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall)
+{
+  if (!wall->phi)
+    return;
+
+  gkyl_array_release(wall->phi);
+  if (wall->projector) {
+    gkyl_eval_on_nodes_release(wall->projector);
+    if (app->use_gpu)
+      gkyl_array_release(wall->phi_host);
+  }
+}

gyrokinetic/apps/gkyl_gk_field_priv.h

@@ -10,7 +10,7 @@
  * This header contains forward declarations for internal field-related
  * functions used by the gyrokinetic application. These functions handle
  * FEM projections, field solves, boundary conditions, energy diagnostics,
- * FLR corrections, polarization potentials, and biased wall configurations.
+ * FLR corrections and polarization potentials.
  */
 
 /** FEM Initialization Functions **/
@@ -147,35 +147,3 @@ gk_field_invert_flr(gkyl_gyrokinetic_app *app, struct gk_field *field,
 void
 gk_field_invert_flr_none(gkyl_gyrokinetic_app *app, struct gk_field *field,
   struct gkyl_array *phi);
-
-/** Biased Wall Functions **/
-
-/**
- * Initialize biased wall boundary condition objects.
- * Sets up machinery for applying time-dependent wall potentials.
- *
- * @param app Gyrokinetic application object.
- * @param f Field object to initialize biased wall for.
- */
-void
-gk_field_biased_wall_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f);
-
-/**
- * Release biased wall resources.
- *
- * @param app Gyrokinetic application object.
- * @param f Field object whose biased wall resources are to be released.
- */
-void
-gk_field_biased_wall_release(const struct gkyl_gyrokinetic_app *app, struct gk_field *f);
-
-/**
- * Calculate and apply the wall potential at the current simulation time.
- * Updates phi at the wall boundaries based on the biased wall configuration.
- *
- * @param app Gyrokinetic application object.
- * @param field Field object containing the potential.
- * @param tm Current simulation time.
- */
-void
-gk_field_calc_phi_wall(gkyl_gyrokinetic_app *app, struct gk_field *field, double tm);