From 958b04a4aa08b9d27f701dd3261481f6c58334f6 Mon Sep 17 00:00:00 2001 From: Maxwell-Rosen Date: Fri, 28 Aug 2026 18:11:02 -0400 Subject: [PATCH 1/3] Refactor biased wall functionality and introduce wall potential management for species --- gyrokinetic/apps/gk_field.c | 6 - gyrokinetic/apps/gk_field_biased_wall.c | 96 ------------ gyrokinetic/apps/gk_species.c | 8 +- gyrokinetic/apps/gk_species_fdot_multiplier.c | 22 ++- gyrokinetic/apps/gk_species_wall_potential.c | 74 ++++++++++ gyrokinetic/apps/gkyl_gk_field_priv.h | 34 +---- gyrokinetic/apps/gkyl_gyrokinetic.h | 13 +- gyrokinetic/apps/gkyl_gyrokinetic_multib.h | 9 -- gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 57 +++++--- gyrokinetic/apps/gyrokinetic.c | 12 +- gyrokinetic/apps/gyrokinetic_multib.c | 27 +--- .../unit/ctest_gk_species_wall_potential.c | 138 ++++++++++++++++++ 12 files changed, 280 insertions(+), 216 deletions(-) delete mode 100644 gyrokinetic/apps/gk_field_biased_wall.c create mode 100644 gyrokinetic/apps/gk_species_wall_potential.c create mode 100644 gyrokinetic/unit/ctest_gk_species_wall_potential.c diff --git a/gyrokinetic/apps/gk_field.c b/gyrokinetic/apps/gk_field.c index 8e81b86f67..108c1c82e8 100644 --- a/gyrokinetic/apps/gk_field.c +++ b/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); } diff --git a/gyrokinetic/apps/gk_field_biased_wall.c b/gyrokinetic/apps/gk_field_biased_wall.c deleted file mode 100644 index 3658b29a71..0000000000 --- a/gyrokinetic/apps/gk_field_biased_wall.c +++ /dev/null @@ -1,96 +0,0 @@ -#include -#include -#include -#include - -void -gk_field_biased_wall_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f) -{ - f->phi_wall_lo = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); - f->has_phi_wall_lo = false; - f->phi_wall_lo_evolve = false; - if (f->info.phi_wall_lo) - { - f->has_phi_wall_lo = true; - if (f->info.phi_wall_lo_evolve) { - f->phi_wall_lo_evolve = f->info.phi_wall_lo_evolve; - } - - f->phi_wall_lo_host = f->phi_wall_lo; - if (app->use_gpu) { - f->phi_wall_lo_host = mkarr(false, f->phi_wall_lo->ncomp, f->phi_wall_lo->size); - } - - f->phi_wall_lo_proj = gkyl_eval_on_nodes_new(&app->grid, &app->basis, - 1, f->info.phi_wall_lo, f->info.phi_wall_lo_ctx); - - // Compute phi_wall_lo at t = 0 - gkyl_eval_on_nodes_advance(f->phi_wall_lo_proj, 0.0, &app->local_ext, f->phi_wall_lo_host); - if (app->use_gpu) { // note: phi_wall_lo_host is same as phi_wall_lo when not on GPUs - gkyl_array_copy(f->phi_wall_lo, f->phi_wall_lo_host); - } - } - - // Set up biased upper wall (same size as electrostatic potential), by default is 0.0 - f->phi_wall_up = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume); - f->has_phi_wall_up = false; - f->phi_wall_up_evolve = false; - if (f->info.phi_wall_up) - { - f->has_phi_wall_up = true; - if (f->info.phi_wall_up_evolve) { - f->phi_wall_up_evolve = f->info.phi_wall_up_evolve; - } - - f->phi_wall_up_host = f->phi_wall_up; - if (app->use_gpu) { - f->phi_wall_up_host = mkarr(false, f->phi_wall_up->ncomp, f->phi_wall_up->size); - } - - f->phi_wall_up_proj = gkyl_eval_on_nodes_new(&app->grid, &app->basis, - 1, f->info.phi_wall_up, f->info.phi_wall_up_ctx); - - // Compute phi_wall_up at t = 0. - gkyl_eval_on_nodes_advance(f->phi_wall_up_proj, 0.0, &app->local_ext, f->phi_wall_up_host); - if (app->use_gpu) { // Note: phi_wall_up_host is same as phi_wall_up when not on GPUs. - gkyl_array_copy(f->phi_wall_up, f->phi_wall_up_host); - } - } -} - -void -gk_field_calc_phi_wall(gkyl_gyrokinetic_app *app, struct gk_field *field, double tm) -{ - if (field->has_phi_wall_lo && field->phi_wall_lo_evolve) { - gkyl_eval_on_nodes_advance(field->phi_wall_lo_proj, tm, &app->local_ext, field->phi_wall_lo_host); - if (app->use_gpu) { - gkyl_array_copy(field->phi_wall_lo, field->phi_wall_lo_host); - } - } - if (field->has_phi_wall_up && field->phi_wall_up_evolve) { - gkyl_eval_on_nodes_advance(field->phi_wall_up_proj, tm, &app->local_ext, field->phi_wall_up_host); - if (app->use_gpu) { - gkyl_array_copy(field->phi_wall_up, field->phi_wall_up_host); - } - } -} - -void -gk_field_biased_wall_release(const struct gkyl_gyrokinetic_app *app, struct gk_field *f) -{ - gkyl_array_release(f->phi_wall_lo); - if (f->has_phi_wall_lo) { - gkyl_eval_on_nodes_release(f->phi_wall_lo_proj); - if (app->use_gpu) { - gkyl_array_release(f->phi_wall_lo_host); - } - } - - gkyl_array_release(f->phi_wall_up); - if (f->has_phi_wall_up) { - gkyl_eval_on_nodes_release(f->phi_wall_up_proj); - if (app->use_gpu) { - gkyl_array_release(f->phi_wall_up_host); - } - } -} diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index 2b93be63cf..67c505057e 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -211,7 +211,7 @@ 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: 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); @@ -233,7 +233,7 @@ 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: 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); @@ -1497,6 +1497,8 @@ 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; } } + + gk_species_wall_potential_init(app, gks); // Species properties metadata. struct gkyl_msgpack_map_elem io_meta_sprop[] = { @@ -2036,6 +2038,8 @@ gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *gks gk_species_fdot_multiplier_release(app, &gks->fdot_mult); + gk_species_wall_potential_release(app, gks); + gk_species_lbo_release(app, &gks->lbo); gk_species_bgk_release(app, &gks->bgk); diff --git a/gyrokinetic/apps/gk_species_fdot_multiplier.c b/gyrokinetic/apps/gk_species_fdot_multiplier.c index 7992e1f6d1..cc9d74a13d 100644 --- a/gyrokinetic/apps/gk_species_fdot_multiplier.c +++ b/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) { diff --git a/gyrokinetic/apps/gk_species_wall_potential.c b/gyrokinetic/apps/gk_species_wall_potential.c new file mode 100644 index 0000000000..befc2c082e --- /dev/null +++ b/gyrokinetic/apps/gk_species_wall_potential.c @@ -0,0 +1,74 @@ +#include +#include + +static void +wall_potential_advance(gkyl_gyrokinetic_app *app, + struct gk_species_wall_potential *wall, double tm) +{ + if (wall->projector) { + 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); + } +} + +static void +wall_potential_init(gkyl_gyrokinetic_app *app, + const struct gkyl_gyrokinetic_bc *bc, struct gk_species_wall_potential *wall) +{ + *wall = (struct gk_species_wall_potential) { }; + 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_potential_advance(app, wall, 0.0); + } +} + +static void +wall_potential_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); + } +} + +void +gk_species_wall_potential_init(gkyl_gyrokinetic_app *app, struct gk_species *species) +{ + int par_dir = app->cdim-1; + wall_potential_init(app, &species->lower_bc[par_dir], &species->phi_wall_lo); + wall_potential_init(app, &species->upper_bc[par_dir], &species->phi_wall_up); +} + +void +gk_species_wall_potential_advance(gkyl_gyrokinetic_app *app, + struct gk_species *species, double tm) +{ + wall_potential_advance(app, &species->phi_wall_lo, tm); + wall_potential_advance(app, &species->phi_wall_up, tm); +} + +void +gk_species_wall_potential_release(const gkyl_gyrokinetic_app *app, + const struct gk_species *species) +{ + wall_potential_release(app, &species->phi_wall_lo); + wall_potential_release(app, &species->phi_wall_up); +} diff --git a/gyrokinetic/apps/gkyl_gk_field_priv.h b/gyrokinetic/apps/gkyl_gk_field_priv.h index da502d8cd9..eee5b8281b 100644 --- a/gyrokinetic/apps/gkyl_gk_field_priv.h +++ b/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); diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h index 65c8da144b..31938b0e63 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic.h @@ -213,7 +213,8 @@ struct gkyl_gyrokinetic_bc { enum gkyl_edge_loc edge; // Which edge this BC is for. enum gkyl_gyrokinetic_bc_type type; // BC type flag. double value[3]; // Meaning depends on type. - void (*aux_profile)(double t, const double *xn, double *fout, void *ctx); // Auxiliary function (e.g. wall potential). + // Auxiliary profile. For a species sheath BC, this is the material-wall potential. + void (*aux_profile)(double t, const double *xn, double *fout, void *ctx); void *aux_ctx; // Context for aux_profile. struct gkyl_gyrokinetic_projection projection; // Projection object input (e.g. for FIXED_FUNC). struct gkyl_gyrokinetic_emission_inp emission; @@ -569,16 +570,6 @@ struct gkyl_gyrokinetic_field { void (*init_field_profile)(double t, const double *xn, double *out, void *ctx); void *init_field_profile_ctx; - void *phi_wall_lo_ctx; // context for biased wall potential on lower wall - // pointer to biased wall potential on lower wall function - void (*phi_wall_lo)(double t, const double *xn, double *phi_wall_lo_out, void *ctx); - bool phi_wall_lo_evolve; // set to true if biased wall potential on lower wall function is time dependent - - void *phi_wall_up_ctx; // context for biased wall potential on upper wall - // pointer to biased wall potential on upper wall function - void (*phi_wall_up)(double t, const double *xn, double *phi_wall_up_out, void *ctx); - bool phi_wall_up_evolve; // set to true if biased wall potential on upper wall function is time dependent - struct gkyl_poisson_bias_line_list *bias_line_list; // Biased lines constraining the solution. }; diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_multib.h b/gyrokinetic/apps/gkyl_gyrokinetic_multib.h index e47288d6f1..3706894eab 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_multib.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_multib.h @@ -110,15 +110,6 @@ struct gkyl_gyrokinetic_multib_field_pb { bool time_rate_diagnostics; // Writes the time rate of change of field energy. - void *phi_wall_lo_ctx; // context for biased wall potential on lower wall - // pointer to biased wall potential on lower wall function - void (*phi_wall_lo)(double t, const double *xn, double *phi_wall_lo_out, void *ctx); - bool phi_wall_lo_evolve; // set to true if biased wall potential on lower wall function is time dependent - - void *phi_wall_up_ctx; // context for biased wall potential on upper wall - // pointer to biased wall potential on upper wall function - void (*phi_wall_up)(double t, const double *xn, double *phi_wall_up_out, void *ctx); - bool phi_wall_up_evolve; // set to true if biased wall potential on upper wall function is time dependent }; // Field input diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 4621a1f499..7d557e423d 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -999,6 +999,13 @@ struct gk_positivity { struct gk_positivity *pos); }; +// Projected material-wall potential associated with a species sheath BC. +struct gk_species_wall_potential { + struct gkyl_array *phi; + struct gkyl_array *phi_host; + gkyl_eval_on_nodes *projector; +}; + // Species data. struct gk_species { struct gkyl_gyrokinetic_species info; // Input data. @@ -1061,6 +1068,8 @@ struct gk_species { // Boundary conditions on lower/upper edges in each direction. struct gkyl_gyrokinetic_bc lower_bc[GKYL_MAX_CDIM], upper_bc[GKYL_MAX_CDIM]; + // Wall potentials supplied by the parallel sheath BCs. + struct gk_species_wall_potential phi_wall_lo, phi_wall_up; // gyrokinetic sheath boundary conditions struct gkyl_bc_sheath_gyrokinetic *bc_sheath_lo; struct gkyl_bc_sheath_gyrokinetic *bc_sheath_up; @@ -1367,18 +1376,6 @@ struct gk_field { gkyl_dynvec integ_energy_dot; // d/dt of integrated energy components. bool is_first_energy_dot_write_call; // flag for d(energy)/dt dynvec written first time - bool has_phi_wall_lo; // flag to indicate there is biased wall potential on lower wall - bool phi_wall_lo_evolve; // flag to indicate biased wall potential on lower wall is time dependent - struct gkyl_array *phi_wall_lo; // biased wall potential on lower wall - struct gkyl_array *phi_wall_lo_host; // host copy for use in IO and projecting - gkyl_eval_on_nodes *phi_wall_lo_proj; // projector for biased wall potential on lower wall - - bool has_phi_wall_up; // flag to indicate there is biased wall potential on upper wall - bool phi_wall_up_evolve; // flag to indicate biased wall potential on upper wall is time dependent - struct gkyl_array *phi_wall_up; // biased wall potential on upper wall - struct gkyl_array *phi_wall_up_host; // host copy for use in IO and projecting - gkyl_eval_on_nodes *phi_wall_up_proj; // projector for biased wall potential on upper wall - // Pointer to function that computes the time rate of change of the energy. void (*calc_energy_dt_func)(gkyl_gyrokinetic_app *app, const struct gk_field *field, double dt, double *energy_reduced); @@ -3339,6 +3336,33 @@ gk_species_calc_int_mom_dt(gkyl_gyrokinetic_app* app, struct gk_species *gks, do */ void gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *s); +/** + * Initialize the wall potentials owned by a species' parallel sheath BCs. + * + * @param app Gyrokinetic app object. + * @param species Species object. + */ +void gk_species_wall_potential_init(gkyl_gyrokinetic_app *app, struct gk_species *species); + +/** + * Evaluate a species' sheath wall-potential profiles at the requested time. + * + * @param app Gyrokinetic app object. + * @param species Species object. + * @param tm Time at which to evaluate the profiles. + */ +void gk_species_wall_potential_advance(gkyl_gyrokinetic_app *app, + struct gk_species *species, double tm); + +/** + * Release a species' sheath wall-potential resources. + * + * @param app Gyrokinetic app object. + * @param species Species object. + */ +void gk_species_wall_potential_release(const gkyl_gyrokinetic_app *app, + const struct gk_species *species); + /** gk_neut_species_moment API */ /** @@ -4134,15 +4158,6 @@ void gk_neut_species_release(const gkyl_gyrokinetic_app* app, const struct gk_ne */ struct gk_field* gk_field_new(struct gkyl_gk *gk, struct gkyl_gyrokinetic_app *app); -/** - * Compute biased wall potentials. - * - * @param app gyrokinetic app object. - * @param field Pointer to field. - * @param tm Time to compute biased wall potentials at. - */ -void gk_field_calc_phi_wall(gkyl_gyrokinetic_app *app, struct gk_field *field, double tm); - /** * Accumulate charge density for Poisson solve. * diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c index 526aeebc67..c18c9bf6f2 100644 --- a/gyrokinetic/apps/gyrokinetic.c +++ b/gyrokinetic/apps/gyrokinetic.c @@ -690,12 +690,6 @@ gyrokinetic_calc_field_enabled(gkyl_gyrokinetic_app* app, double tcurr, // Compute electrostatic potential from gyrokinetic Poisson's equation. gk_field_accumulate_rho_c(app, app->field, fin, bflux); - // Compute biased wall potential if present and time-dependent. - // Note: biased wall potential use eval_on_nodes. - // so does copy to GPU every call if app->use_gpu = true. - if (app->field->phi_wall_lo_evolve || app->field->phi_wall_up_evolve) - gk_field_calc_phi_wall(app, app->field, tcurr); - // Solve the field equation. gk_field_rhs(app, app->field); app->stat.field_tm += gkyl_time_diff_now_sec(wtm); @@ -982,6 +976,9 @@ void gyrokinetic_calc_field(gkyl_gyrokinetic_app* app, double tcurr, const struct gkyl_array *fin[], struct gkyl_array **bflux[]) { + for (int i=0; inum_species; ++i) + gk_species_wall_potential_advance(app, &app->species[i], tcurr); + app->calc_field_func(app, tcurr, fin, bflux); } @@ -1061,6 +1058,9 @@ gkyl_gyrokinetic_app_apply_ic(gkyl_gyrokinetic_app* app, double t0) for (int i=0; inum_neut_species; ++i) gkyl_gyrokinetic_app_apply_ic_cross_neut_species(app, i, t0); + for (int i=0; inum_species; ++i) + gk_species_wall_potential_advance(app, &app->species[i], t0); + // Compute the fields and apply BCs. struct gkyl_array *distf[app->num_species]; struct gkyl_array **bflux[app->num_species]; diff --git a/gyrokinetic/apps/gyrokinetic_multib.c b/gyrokinetic/apps/gyrokinetic_multib.c index de0352dee2..cfccffc254 100644 --- a/gyrokinetic/apps/gyrokinetic_multib.c +++ b/gyrokinetic/apps/gyrokinetic_multib.c @@ -193,10 +193,8 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid, species_inp.polarization_density = sp_pb->polarization_density; // By default, skip BCs altogether. - for (int i=0; i<2*GKYL_MAX_CDIM; i++) { + for (int i=0; i<2*GKYL_MAX_CDIM; i++) species_inp.bcs[i].type = GKYL_BC_GK_SKIP; - species_inp.bcs[i].type = GKYL_BC_GK_SKIP; - } int pardir = cdim-1; int num_below = gkyl_multib_conn_get_num_connected(mbapp->block_topo, @@ -285,10 +283,8 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid, neut_species_inp.source = nsp_pb->source; // By default, skip BCs altogether. - for (int i=0; i<2*GKYL_MAX_CDIM; i++) { + for (int i=0; i<2*GKYL_MAX_CDIM; i++) neut_species_inp.bcs[i].type = GKYL_BC_GK_SKIP; - neut_species_inp.bcs[i].type = GKYL_BC_GK_SKIP; - } // Set species physical BCs. int bc_count_nsp[num_blocks]; @@ -353,15 +349,6 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid, const struct gkyl_gyrokinetic_multib_field_pb *fld_pb = &fld->blocks[0]; // Choose proper block-specific field input. - if (!fld->duplicate_across_blocks) { - for (int i=0; iblocks[i].block_id) { - const struct gkyl_gyrokinetic_multib_field_pb *fld_pb = &fld->blocks[i]; - break; - } - } - } - if (!fld->duplicate_across_blocks) { for (int i=0; iblocks[i].block_id) { @@ -374,15 +361,6 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid, field_inp.polarization_bmag = fld_pb->polarization_bmag ? fld_pb->polarization_bmag : mbapp->bmag_ref; field_inp.kperpSq = fld_pb->kperpSq; field_inp.time_rate_diagnostics = fld_pb->time_rate_diagnostics; - - field_inp.phi_wall_lo_ctx = fld_pb->phi_wall_lo_ctx; - field_inp.phi_wall_lo = fld_pb->phi_wall_lo; - field_inp.phi_wall_lo_evolve = fld_pb->phi_wall_lo_evolve; - - field_inp.phi_wall_up_ctx = fld_pb->phi_wall_up_ctx; - field_inp.phi_wall_up = fld_pb->phi_wall_up; - field_inp.phi_wall_up_evolve = fld_pb->phi_wall_up_evolve; - // Copy field input into app input. memcpy(&app_inp.field, &field_inp, sizeof(struct gkyl_gyrokinetic_field)); @@ -932,6 +910,7 @@ gyrokinetic_multib_apply_bc(struct gkyl_gyrokinetic_multib_app* app, double tcur int li_charged = b * app->num_species; int li_neut = b * app->num_neut_species; for (int i=0; inum_species; ++i) { + gk_species_wall_potential_advance(sbapp, &sbapp->species[i], tcurr); gk_species_apply_bc(sbapp, &sbapp->species[i], distf[li_charged+i]); } for (int i=0; inum_neut_species; ++i) { diff --git a/gyrokinetic/unit/ctest_gk_species_wall_potential.c b/gyrokinetic/unit/ctest_gk_species_wall_potential.c new file mode 100644 index 0000000000..8944f5be72 --- /dev/null +++ b/gyrokinetic/unit/ctest_gk_species_wall_potential.c @@ -0,0 +1,138 @@ +#include + +#include + +#include + +struct wall_profile_ctx { + double offset; + double time_slope; + double space_slope; +}; + +static void +wall_profile(double t, const double *xn, double *fout, void *ctx) +{ + const struct wall_profile_ctx *profile = ctx; + fout[0] = profile->offset + profile->time_slope*t + profile->space_slope*xn[0]; +} + +static gkyl_gyrokinetic_app +make_app(void) +{ + gkyl_gyrokinetic_app app = { + .cdim = 1, + .use_gpu = false, + }; + double lower[] = { -1.0 }, upper[] = { 1.0 }; + int cells[] = { 4 }, ghost[] = { 1 }; + + gkyl_rect_grid_init(&app.grid, app.cdim, lower, upper, cells); + gkyl_cart_modal_serendip(&app.basis, app.cdim, 1); + gkyl_create_grid_ranges(&app.grid, ghost, &app.local_ext, &app.local); + + return app; +} + +static void +check_profile(const gkyl_gyrokinetic_app *app, const struct gkyl_array *phi, + const struct wall_profile_ctx *profile, double tm) +{ + struct gkyl_range_iter iter; + gkyl_range_iter_init(&iter, &app->local_ext); + while (gkyl_range_iter_next(&iter)) { + long lidx = gkyl_range_idx(&app->local_ext, iter.idx); + const double *phi_c = gkyl_array_cfetch(phi, lidx); + double xc[GKYL_MAX_DIM]; + gkyl_rect_grid_cell_center(&app->grid, iter.idx, xc); + + const double eta[] = { -1.0, 0.0, 1.0 }; + for (int i=0; i<3; ++i) { + double x = xc[0] + 0.5*app->grid.dx[0]*eta[i]; + double expected = profile->offset + profile->time_slope*tm + profile->space_slope*x; + TEST_CHECK(fabs(app->basis.eval_expand(&eta[i], phi_c)-expected) < 1e-12); + } + } +} + +static void +test_species_wall_profiles(void) +{ + gkyl_gyrokinetic_app app = make_app(); + struct wall_profile_ctx lower_ctx = { + .offset = 1.5, + .time_slope = 2.0, + .space_slope = -0.25, + }; + struct wall_profile_ctx upper_ctx = { + .offset = -0.5, + .time_slope = -1.0, + .space_slope = 0.75, + }; + struct gk_species species = { + .lower_bc[0] = { + .type = GKYL_BC_GK_SPECIES_SHEATH, + .aux_profile = wall_profile, + .aux_ctx = &lower_ctx, + }, + .upper_bc[0] = { + .type = GKYL_BC_GK_SPECIES_SHEATH, + .aux_profile = wall_profile, + .aux_ctx = &upper_ctx, + }, + }; + + gk_species_wall_potential_init(&app, &species); + TEST_ASSERT(species.phi_wall_lo.phi != 0); + TEST_ASSERT(species.phi_wall_up.phi != 0); + TEST_ASSERT(species.phi_wall_lo.projector != 0); + TEST_ASSERT(species.phi_wall_up.projector != 0); + check_profile(&app, species.phi_wall_lo.phi, &lower_ctx, 0.0); + check_profile(&app, species.phi_wall_up.phi, &upper_ctx, 0.0); + + double tm = 0.625; + gk_species_wall_potential_advance(&app, &species, tm); + check_profile(&app, species.phi_wall_lo.phi, &lower_ctx, tm); + check_profile(&app, species.phi_wall_up.phi, &upper_ctx, tm); + + gk_species_wall_potential_release(&app, &species); +} + +static void +test_grounded_and_non_sheath_walls(void) +{ + gkyl_gyrokinetic_app app = make_app(); + struct gk_species species = { + .lower_bc[0] = { + .type = GKYL_BC_GK_SPECIES_SHEATH, + }, + .upper_bc[0] = { + .type = GKYL_BC_GK_SPECIES_REFLECT, + }, + }; + + gk_species_wall_potential_init(&app, &species); + TEST_ASSERT(species.phi_wall_lo.phi != 0); + TEST_CHECK(species.phi_wall_lo.phi_host == species.phi_wall_lo.phi); + TEST_CHECK(species.phi_wall_lo.projector == 0); + TEST_CHECK(species.phi_wall_up.phi == 0); + TEST_CHECK(species.phi_wall_up.phi_host == 0); + TEST_CHECK(species.phi_wall_up.projector == 0); + + struct gkyl_range_iter iter; + gkyl_range_iter_init(&iter, &app.local_ext); + while (gkyl_range_iter_next(&iter)) { + long lidx = gkyl_range_idx(&app.local_ext, iter.idx); + const double *phi_c = gkyl_array_cfetch(species.phi_wall_lo.phi, lidx); + for (int k=0; k Date: Mon, 31 Aug 2026 13:38:35 -0400 Subject: [PATCH 2/3] Refactor wall potential management in species: introduce phi_wall functions and update boundary condition applications --- gyrokinetic/apps/gk_species.c | 20 +++++-- ...wall_potential.c => gk_species_phi_wall.c} | 56 ++++++------------- gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 30 +++++----- gyrokinetic/apps/gyrokinetic.c | 12 +--- gyrokinetic/apps/gyrokinetic_multib.c | 14 ++++- .../unit/ctest_gk_species_wall_potential.c | 16 ++++-- 6 files changed, 72 insertions(+), 76 deletions(-) rename gyrokinetic/apps/{gk_species_wall_potential.c => gk_species_phi_wall.c} (52%) diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index 67c505057e..0e674aebe0 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/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,6 +211,7 @@ 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, species->phi_wall_lo.phi, f, &app->local); break; @@ -232,6 +234,7 @@ 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, species->phi_wall_up.phi, f, &app->local); break; @@ -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 } @@ -1498,7 +1502,9 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st } } - gk_species_wall_potential_init(app, gks); + int wall_dir = app->cdim-1; + gk_species_phi_wall_init(app, &gks->lower_bc[wall_dir], &gks->phi_wall_lo); + gk_species_phi_wall_init(app, &gks->upper_bc[wall_dir], &gks->phi_wall_up); // Species properties metadata. struct gkyl_msgpack_map_elem io_meta_sprop[] = { @@ -1937,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 @@ -2038,7 +2045,8 @@ gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *gks gk_species_fdot_multiplier_release(app, &gks->fdot_mult); - gk_species_wall_potential_release(app, gks); + gk_species_phi_wall_release(app, &gks->phi_wall_lo); + gk_species_phi_wall_release(app, &gks->phi_wall_up); gk_species_lbo_release(app, &gks->lbo); diff --git a/gyrokinetic/apps/gk_species_wall_potential.c b/gyrokinetic/apps/gk_species_phi_wall.c similarity index 52% rename from gyrokinetic/apps/gk_species_wall_potential.c rename to gyrokinetic/apps/gk_species_phi_wall.c index befc2c082e..c67e1d4c12 100644 --- a/gyrokinetic/apps/gk_species_wall_potential.c +++ b/gyrokinetic/apps/gk_species_phi_wall.c @@ -1,19 +1,8 @@ #include #include -static void -wall_potential_advance(gkyl_gyrokinetic_app *app, - struct gk_species_wall_potential *wall, double tm) -{ - if (wall->projector) { - 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); - } -} - -static void -wall_potential_init(gkyl_gyrokinetic_app *app, +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) { }; @@ -30,12 +19,23 @@ wall_potential_init(gkyl_gyrokinetic_app *app, wall->projector = gkyl_eval_on_nodes_new(&app->grid, &app->basis, 1, bc->aux_profile, bc->aux_ctx); - wall_potential_advance(app, wall, 0.0); + gk_species_phi_wall_advance(app, wall, 0.0); } } -static void -wall_potential_release(const gkyl_gyrokinetic_app *app, +void +gk_species_phi_wall_advance(gkyl_gyrokinetic_app *app, + const struct gk_species_wall_potential *wall, double tm) +{ + if (wall->projector) { + 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_release(const gkyl_gyrokinetic_app *app, const struct gk_species_wall_potential *wall) { if (!wall->phi) @@ -48,27 +48,3 @@ wall_potential_release(const gkyl_gyrokinetic_app *app, gkyl_array_release(wall->phi_host); } } - -void -gk_species_wall_potential_init(gkyl_gyrokinetic_app *app, struct gk_species *species) -{ - int par_dir = app->cdim-1; - wall_potential_init(app, &species->lower_bc[par_dir], &species->phi_wall_lo); - wall_potential_init(app, &species->upper_bc[par_dir], &species->phi_wall_up); -} - -void -gk_species_wall_potential_advance(gkyl_gyrokinetic_app *app, - struct gk_species *species, double tm) -{ - wall_potential_advance(app, &species->phi_wall_lo, tm); - wall_potential_advance(app, &species->phi_wall_up, tm); -} - -void -gk_species_wall_potential_release(const gkyl_gyrokinetic_app *app, - const struct gk_species *species) -{ - wall_potential_release(app, &species->phi_wall_lo); - wall_potential_release(app, &species->phi_wall_up); -} diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 7d557e423d..105e2fb0c4 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -1135,7 +1135,7 @@ struct gk_species { double (*rhs_implicit_func)(gkyl_gyrokinetic_app *app, struct gk_species *species, const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms, double dt); void (*bc_func)(gkyl_gyrokinetic_app *app, const struct gk_species *species, - struct gkyl_array *f); + double tm, struct gkyl_array *f); void (*release_func)(const gkyl_gyrokinetic_app* app, const struct gk_species *s); void (*step_f_func)(struct gkyl_array* out, double dt, const struct gkyl_array* inp); void (*combine_func)(struct gkyl_array *out, double c1, @@ -3240,9 +3240,11 @@ void gk_species_apply_pos_shift(gkyl_gyrokinetic_app* app, struct gk_species *gk * * @param app gyrokinetic app object. * @param species Pointer to species. + * @param tm Time at which to apply the BCs. * @param f Field to apply BCs. */ -void gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f); +void gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species, + double tm, struct gkyl_array *f); /** * Fill stat object in app with collision timers. @@ -3337,31 +3339,33 @@ gk_species_calc_int_mom_dt(gkyl_gyrokinetic_app* app, struct gk_species *gks, do void gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *s); /** - * Initialize the wall potentials owned by a species' parallel sheath BCs. + * Initialize a wall potential owned by a species sheath BC. * * @param app Gyrokinetic app object. - * @param species Species object. + * @param bc Species sheath BC input. + * @param wall Wall-potential object. */ -void gk_species_wall_potential_init(gkyl_gyrokinetic_app *app, struct gk_species *species); +void gk_species_phi_wall_init(gkyl_gyrokinetic_app *app, + const struct gkyl_gyrokinetic_bc *bc, struct gk_species_wall_potential *wall); /** - * Evaluate a species' sheath wall-potential profiles at the requested time. + * Evaluate a species sheath wall-potential profile at the requested time. * * @param app Gyrokinetic app object. - * @param species Species object. + * @param wall Wall-potential object. * @param tm Time at which to evaluate the profiles. */ -void gk_species_wall_potential_advance(gkyl_gyrokinetic_app *app, - struct gk_species *species, double tm); +void gk_species_phi_wall_advance(gkyl_gyrokinetic_app *app, + const struct gk_species_wall_potential *wall, double tm); /** - * Release a species' sheath wall-potential resources. + * Release a species sheath wall-potential resource. * * @param app Gyrokinetic app object. - * @param species Species object. + * @param wall Wall-potential object. */ -void gk_species_wall_potential_release(const gkyl_gyrokinetic_app *app, - const struct gk_species *species); +void gk_species_phi_wall_release(const gkyl_gyrokinetic_app *app, + const struct gk_species_wall_potential *wall); /** gk_neut_species_moment API */ diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c index c18c9bf6f2..4edc456b48 100644 --- a/gyrokinetic/apps/gyrokinetic.c +++ b/gyrokinetic/apps/gyrokinetic.c @@ -976,9 +976,6 @@ void gyrokinetic_calc_field(gkyl_gyrokinetic_app* app, double tcurr, const struct gkyl_array *fin[], struct gkyl_array **bflux[]) { - for (int i=0; inum_species; ++i) - gk_species_wall_potential_advance(app, &app->species[i], tcurr); - app->calc_field_func(app, tcurr, fin, bflux); } @@ -996,7 +993,7 @@ gyrokinetic_calc_field_and_apply_bc(gkyl_gyrokinetic_app* app, double tcurr, // Apply boundary conditions. struct timespec wst = gkyl_wall_clock(); for (int i=0; inum_species; ++i) { - gk_species_apply_bc(app, &app->species[i], distf[i]); + gk_species_apply_bc(app, &app->species[i], tcurr, distf[i]); } for (int i=0; inum_neut_species; ++i) { gk_neut_species_apply_bc(app, &app->neut_species[i], distf_neut[i]); @@ -1058,9 +1055,6 @@ gkyl_gyrokinetic_app_apply_ic(gkyl_gyrokinetic_app* app, double t0) for (int i=0; inum_neut_species; ++i) gkyl_gyrokinetic_app_apply_ic_cross_neut_species(app, i, t0); - for (int i=0; inum_species; ++i) - gk_species_wall_potential_advance(app, &app->species[i], t0); - // Compute the fields and apply BCs. struct gkyl_array *distf[app->num_species]; struct gkyl_array **bflux[app->num_species]; @@ -1115,7 +1109,7 @@ gkyl_gyrokinetic_app_apply_ic(gkyl_gyrokinetic_app* app, double t0) // Apply boundary conditions. for (int i=0; inum_species; ++i) { - gk_species_apply_bc(app, &app->species[i], distf[i]); + gk_species_apply_bc(app, &app->species[i], t0, distf[i]); } for (int i=0; inum_neut_species; ++i) { if (!app->neut_species[i].info.is_static) { @@ -3506,7 +3500,7 @@ gkyl_gyrokinetic_app_read_from_frame(gkyl_gyrokinetic_app *app, int frame) // Apply boundary conditions. for (int i=0; inum_species; ++i) { - gk_species_apply_bc(app, &app->species[i], distf[i]); + gk_species_apply_bc(app, &app->species[i], rstat.stime, distf[i]); } for (int i=0; inum_neut_species; ++i) { gk_neut_species_apply_bc(app, &app->neut_species[i], distf_neut[i]); diff --git a/gyrokinetic/apps/gyrokinetic_multib.c b/gyrokinetic/apps/gyrokinetic_multib.c index cfccffc254..792032fb92 100644 --- a/gyrokinetic/apps/gyrokinetic_multib.c +++ b/gyrokinetic/apps/gyrokinetic_multib.c @@ -335,7 +335,7 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid, field_inp.poisson_bcs[i].value[k] = -1.0e3; } } - + for (int i=0; inum_physical_bcs; i++) { if (bid == fld->bcs[i].bidx) { struct gkyl_gyrokinetic_bc *bc_curr = gk_fetch_bc_with_dir_edge(field_inp.poisson_bcs, 2*cdim, fld->bcs[i].dir, fld->bcs[i].edge); @@ -349,6 +349,15 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid, const struct gkyl_gyrokinetic_multib_field_pb *fld_pb = &fld->blocks[0]; // Choose proper block-specific field input. + if (!fld->duplicate_across_blocks) { + for (int i=0; iblocks[i].block_id) { + const struct gkyl_gyrokinetic_multib_field_pb *fld_pb = &fld->blocks[i]; + break; + } + } + } + if (!fld->duplicate_across_blocks) { for (int i=0; iblocks[i].block_id) { @@ -910,8 +919,7 @@ gyrokinetic_multib_apply_bc(struct gkyl_gyrokinetic_multib_app* app, double tcur int li_charged = b * app->num_species; int li_neut = b * app->num_neut_species; for (int i=0; inum_species; ++i) { - gk_species_wall_potential_advance(sbapp, &sbapp->species[i], tcurr); - gk_species_apply_bc(sbapp, &sbapp->species[i], distf[li_charged+i]); + gk_species_apply_bc(sbapp, &sbapp->species[i], tcurr, distf[li_charged+i]); } for (int i=0; inum_neut_species; ++i) { gk_neut_species_apply_bc(sbapp, &sbapp->neut_species[i], distf_neut[li_neut+i]); diff --git a/gyrokinetic/unit/ctest_gk_species_wall_potential.c b/gyrokinetic/unit/ctest_gk_species_wall_potential.c index 8944f5be72..498f3827e3 100644 --- a/gyrokinetic/unit/ctest_gk_species_wall_potential.c +++ b/gyrokinetic/unit/ctest_gk_species_wall_potential.c @@ -82,7 +82,8 @@ test_species_wall_profiles(void) }, }; - gk_species_wall_potential_init(&app, &species); + gk_species_phi_wall_init(&app, &species.lower_bc[0], &species.phi_wall_lo); + gk_species_phi_wall_init(&app, &species.upper_bc[0], &species.phi_wall_up); TEST_ASSERT(species.phi_wall_lo.phi != 0); TEST_ASSERT(species.phi_wall_up.phi != 0); TEST_ASSERT(species.phi_wall_lo.projector != 0); @@ -91,11 +92,14 @@ test_species_wall_profiles(void) check_profile(&app, species.phi_wall_up.phi, &upper_ctx, 0.0); double tm = 0.625; - gk_species_wall_potential_advance(&app, &species, tm); + gk_species_phi_wall_advance(&app, &species.phi_wall_lo, tm); check_profile(&app, species.phi_wall_lo.phi, &lower_ctx, tm); + check_profile(&app, species.phi_wall_up.phi, &upper_ctx, 0.0); + gk_species_phi_wall_advance(&app, &species.phi_wall_up, tm); check_profile(&app, species.phi_wall_up.phi, &upper_ctx, tm); - gk_species_wall_potential_release(&app, &species); + gk_species_phi_wall_release(&app, &species.phi_wall_lo); + gk_species_phi_wall_release(&app, &species.phi_wall_up); } static void @@ -111,7 +115,8 @@ test_grounded_and_non_sheath_walls(void) }, }; - gk_species_wall_potential_init(&app, &species); + gk_species_phi_wall_init(&app, &species.lower_bc[0], &species.phi_wall_lo); + gk_species_phi_wall_init(&app, &species.upper_bc[0], &species.phi_wall_up); TEST_ASSERT(species.phi_wall_lo.phi != 0); TEST_CHECK(species.phi_wall_lo.phi_host == species.phi_wall_lo.phi); TEST_CHECK(species.phi_wall_lo.projector == 0); @@ -128,7 +133,8 @@ test_grounded_and_non_sheath_walls(void) TEST_CHECK(phi_c[k] == 0.0); } - gk_species_wall_potential_release(&app, &species); + gk_species_phi_wall_release(&app, &species.phi_wall_lo); + gk_species_phi_wall_release(&app, &species.phi_wall_up); } TEST_LIST = { From 3a97aaa01994fe51ba59b95bcbd702ada8ae0500 Mon Sep 17 00:00:00 2001 From: Maxwell-Rosen Date: Mon, 31 Aug 2026 14:46:21 -0400 Subject: [PATCH 3/3] Refactor wall potential management: implement advance functions and update initialization in species --- gyrokinetic/apps/gk_species.c | 11 +++----- gyrokinetic/apps/gk_species_phi_wall.c | 26 ++++++++++++++----- gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 2 ++ .../unit/ctest_gk_species_wall_potential.c | 7 +++++ 4 files changed, 33 insertions(+), 13 deletions(-) diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index 0e674aebe0..1c1b68d641 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -681,6 +681,7 @@ gk_species_release_dynamic(const gkyl_gyrokinetic_app* app, const struct gk_spec for (int d=0; dcdim; ++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); @@ -694,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); @@ -881,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); @@ -941,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); @@ -1502,10 +1506,6 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st } } - int wall_dir = app->cdim-1; - gk_species_phi_wall_init(app, &gks->lower_bc[wall_dir], &gks->phi_wall_lo); - gk_species_phi_wall_init(app, &gks->upper_bc[wall_dir], &gks->phi_wall_up); - // Species properties metadata. struct gkyl_msgpack_map_elem io_meta_sprop[] = { { .key = "mass", .elem_type = GKYL_MP_DOUBLE, .dval = gks->info.mass }, @@ -2045,9 +2045,6 @@ gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *gks gk_species_fdot_multiplier_release(app, &gks->fdot_mult); - gk_species_phi_wall_release(app, &gks->phi_wall_lo); - gk_species_phi_wall_release(app, &gks->phi_wall_up); - gk_species_lbo_release(app, &gks->lbo); gk_species_bgk_release(app, &gks->bgk); diff --git a/gyrokinetic/apps/gk_species_phi_wall.c b/gyrokinetic/apps/gk_species_phi_wall.c index c67e1d4c12..c5f211947b 100644 --- a/gyrokinetic/apps/gk_species_phi_wall.c +++ b/gyrokinetic/apps/gk_species_phi_wall.c @@ -1,11 +1,28 @@ #include #include +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) { }; + *wall = (struct gk_species_wall_potential) { + .advance_func = gk_species_phi_wall_advance_disabled, + }; if (bc->type != GKYL_BC_GK_SPECIES_SHEATH) return; @@ -19,6 +36,7 @@ gk_species_phi_wall_init(gkyl_gyrokinetic_app *app, 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); } } @@ -27,11 +45,7 @@ void gk_species_phi_wall_advance(gkyl_gyrokinetic_app *app, const struct gk_species_wall_potential *wall, double tm) { - if (wall->projector) { - 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); - } + wall->advance_func(app, wall, tm); } void diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 105e2fb0c4..fb15362ffe 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -1004,6 +1004,8 @@ struct gk_species_wall_potential { struct gkyl_array *phi; struct gkyl_array *phi_host; gkyl_eval_on_nodes *projector; + void (*advance_func)(gkyl_gyrokinetic_app *app, + const struct gk_species_wall_potential *wall, double tm); }; // Species data. diff --git a/gyrokinetic/unit/ctest_gk_species_wall_potential.c b/gyrokinetic/unit/ctest_gk_species_wall_potential.c index 498f3827e3..7b10f8abdb 100644 --- a/gyrokinetic/unit/ctest_gk_species_wall_potential.c +++ b/gyrokinetic/unit/ctest_gk_species_wall_potential.c @@ -88,6 +88,8 @@ test_species_wall_profiles(void) TEST_ASSERT(species.phi_wall_up.phi != 0); TEST_ASSERT(species.phi_wall_lo.projector != 0); TEST_ASSERT(species.phi_wall_up.projector != 0); + TEST_ASSERT(species.phi_wall_lo.advance_func != 0); + TEST_ASSERT(species.phi_wall_up.advance_func != 0); check_profile(&app, species.phi_wall_lo.phi, &lower_ctx, 0.0); check_profile(&app, species.phi_wall_up.phi, &upper_ctx, 0.0); @@ -120,9 +122,14 @@ test_grounded_and_non_sheath_walls(void) TEST_ASSERT(species.phi_wall_lo.phi != 0); TEST_CHECK(species.phi_wall_lo.phi_host == species.phi_wall_lo.phi); TEST_CHECK(species.phi_wall_lo.projector == 0); + TEST_ASSERT(species.phi_wall_lo.advance_func != 0); TEST_CHECK(species.phi_wall_up.phi == 0); TEST_CHECK(species.phi_wall_up.phi_host == 0); TEST_CHECK(species.phi_wall_up.projector == 0); + TEST_ASSERT(species.phi_wall_up.advance_func != 0); + + gk_species_phi_wall_advance(&app, &species.phi_wall_lo, 1.0); + gk_species_phi_wall_advance(&app, &species.phi_wall_up, 1.0); struct gkyl_range_iter iter; gkyl_range_iter_init(&iter, &app.local_ext);