From a35ddcdc3d63b96088387747947a63354db6d1cc Mon Sep 17 00:00:00 2001 From: manauref Date: Mon, 16 Mar 2026 09:48:26 -0700 Subject: [PATCH 1/3] Group density, momentum and energy terms in EIRENE coupling source. Add integrated diagnostics of the whole term. Not yet tested. --- gyrokinetic/apps/gk_species_source_bgk.c | 213 +++++++++++++++-------- gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 28 +++ 2 files changed, 166 insertions(+), 75 deletions(-) diff --git a/gyrokinetic/apps/gk_species_source_bgk.c b/gyrokinetic/apps/gk_species_source_bgk.c index 464e43e58..2312095a4 100644 --- a/gyrokinetic/apps/gk_species_source_bgk.c +++ b/gyrokinetic/apps/gk_species_source_bgk.c @@ -62,15 +62,14 @@ gk_species_source_bgk_rhs_default_enabled(gkyl_gyrokinetic_app *app, struct gk_s } static void -gk_species_source_bgk_rhs_density(gkyl_gyrokinetic_app *app, struct gk_species *species, - struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *rhs) +gks_src_bgk_rhs_accumulate_density_maxwellian(gkyl_gyrokinetic_app *app, struct gk_species *species, + struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *out) { struct timespec wst = gkyl_wall_clock(); // Compute Maxwellian moments (n, u_par, T/m). gk_species_moment_calc(&species->lte.moms, species->local, app->local, fin); - gkyl_dg_div_op_range(species->lte.moms.mem_geo, app->basis, 0, species->lte.moms.marr, - 0, species->lte.moms.marr, 0, app->gk_geom->geo_int.jacobgeo, &app->local); + gk_species_moment_diag_jacobgeo_div(app, &species->lte.moms, species->lte.moms.marr, species->lte.moms.marr); // Set a cap gkyl_array_set_offset(src->Jrate_cap, src->damping_factor, species->lte.moms.marr, 0*app->basis.num_basis); @@ -82,40 +81,22 @@ gk_species_source_bgk_rhs_density(gkyl_gyrokinetic_app *app, struct gk_species * // Set the LTE moments for projection and project gkyl_array_set_offset(species->lte.moms.marr, 1.0, src->Jrate_mom, 0*app->basis.num_basis); gk_species_lte_from_moms(app, species, &species->lte, species->lte.moms.marr); - // Multiply the Maxwellian by Jrate. - gkyl_dg_mul_conf_phase_op_range(&app->basis, &species->basis, src->Jrate_fmax, - src->Jrate, species->lte.f_lte, &app->local, &species->local); - // Assemble the BGK-like term and add it to rhs. - gkyl_bgk_collisions_advance(src->bgk_op, &app->local, &species->local, - src->rate, src->Jrate_fmax, fin, src->implicit_step, src->dt_implicit, rhs, species->cflrate); - app->stat.species_source_bgk_tm += gkyl_time_diff_now_sec(wst); + // Add this Maxwellian onto the average of Maxwellians. + gkyl_array_accumulate_range(out, 1.0/3.0, species->lte.f_lte, &species->local); - // Temporaty hack, write out the M0 moments - //gkyl_array_accumulate(src->Jrate_fmax, -1.0/src->coupling_time, fin); - //for (int m=0; minfo.num_diag_moments; ++m) { - // gk_species_moment_calc(&species->moms[m], species->local, app->local, src->Jrate_fmax); - // gkyl_dg_div_op_range(species->moms[m].mem_geo, app->basis, 0, species->moms[m].marr, 0, species->moms[m].marr, 0, app->gk_geom->geo_int.jacobgeo, &app->local); - // if (app->use_gpu) - // gkyl_array_copy(species->moms[m].marr_host, species->moms[m].marr); - // const char *fmt = "%s-%s_%s_nsource.gkyl"; - // int sz = gkyl_calc_strlen(fmt, app->name, species->info.name, gkyl_distribution_moments_strs[species->info.diag_moments[m]]); - // char fileNm[sz+1]; // ensures no buffer overflow - // snprintf(fileNm, sizeof fileNm, fmt, app->name, species->info.name, gkyl_distribution_moments_strs[species->info.diag_moments[m]]); - // gkyl_comm_array_write(app->comm, &app->grid, &app->local, 0, species->moms[m].marr_host, fileNm); - //} + app->stat.species_source_bgk_tm += gkyl_time_diff_now_sec(wst); } static void -gk_species_source_bgk_rhs_momentum(gkyl_gyrokinetic_app *app, struct gk_species *species, - struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *rhs) +gks_src_bgk_rhs_accumulate_momentum_maxwellian(gkyl_gyrokinetic_app *app, struct gk_species *species, + struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *out) { struct timespec wst = gkyl_wall_clock(); // Compute Maxwellian moments (n, u_par, T/m). gk_species_moment_calc(&species->lte.moms, species->local, app->local, fin); - gkyl_dg_div_op_range(species->lte.moms.mem_geo, app->basis, 0, species->lte.moms.marr, - 0, species->lte.moms.marr, 0, app->gk_geom->geo_int.jacobgeo, &app->local); + gk_species_moment_diag_jacobgeo_div(app, &species->lte.moms, species->lte.moms.marr, species->lte.moms.marr); // Divide M1dot by the rate and density then add on upar gkyl_dg_div_op_range(species->lte.moms.mem_geo, app->basis, 0, src->Jrate_mom, 0, src->M1dot, 0, src->rate, &app->local); @@ -124,40 +105,22 @@ gk_species_source_bgk_rhs_momentum(gkyl_gyrokinetic_app *app, struct gk_species // Set the LTE moments for projection and project gkyl_array_set_offset(species->lte.moms.marr, 1.0, src->Jrate_mom, 1*app->basis.num_basis); gk_species_lte_from_moms(app, species, &species->lte, species->lte.moms.marr); - // Multiply the Maxwellian by Jrate. - gkyl_dg_mul_conf_phase_op_range(&app->basis, &species->basis, src->Jrate_fmax, - src->Jrate, species->lte.f_lte, &app->local, &species->local); - // Assemble the BGK-like term and add it to rhs. - gkyl_bgk_collisions_advance(src->bgk_op, &app->local, &species->local, - src->rate, src->Jrate_fmax, fin, src->implicit_step, src->dt_implicit, rhs, species->cflrate); - app->stat.species_source_bgk_tm += gkyl_time_diff_now_sec(wst); + // Add this Maxwellian onto the average of Maxwellians. + gkyl_array_accumulate_range(out, 1.0/3.0, species->lte.f_lte, &species->local); - // Temporaty hack, write out the M2 moment - //gkyl_array_accumulate(src->Jrate_fmax, -1.0/src->coupling_time, fin); // Testing a factor of 2 - //for (int m=0; minfo.num_diag_moments; ++m) { - // gk_species_moment_calc(&species->moms[m], species->local, app->local, src->Jrate_fmax); - // gkyl_dg_div_op_range(species->moms[m].mem_geo, app->basis, 0, species->moms[m].marr, 0, species->moms[m].marr, 0, app->gk_geom->geo_int.jacobgeo, &app->local); - // if (app->use_gpu) - // gkyl_array_copy(species->moms[m].marr_host, species->moms[m].marr); - // const char *fmt = "%s-%s_%s_msource.gkyl"; - // int sz = gkyl_calc_strlen(fmt, app->name, species->info.name, gkyl_distribution_moments_strs[species->info.diag_moments[m]]); - // char fileNm[sz+1]; // ensures no buffer overflow - // snprintf(fileNm, sizeof fileNm, fmt, app->name, species->info.name, gkyl_distribution_moments_strs[species->info.diag_moments[m]]); - // gkyl_comm_array_write(app->comm, &app->grid, &app->local, 0, species->moms[m].marr_host, fileNm); - //} + app->stat.species_source_bgk_tm += gkyl_time_diff_now_sec(wst); } static void -gk_species_source_bgk_rhs_energy(gkyl_gyrokinetic_app *app, struct gk_species *species, - struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *rhs) +gks_src_bgk_rhs_accumulate_energy_maxwellian(gkyl_gyrokinetic_app *app, struct gk_species *species, + struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *out) { struct timespec wst = gkyl_wall_clock(); // Compute Maxwellian moments (n, u_par, T/m). gk_species_moment_calc(&species->lte.moms, species->local, app->local, fin); - gkyl_dg_div_op_range(species->lte.moms.mem_geo, app->basis, 0, species->lte.moms.marr, - 0, species->lte.moms.marr, 0, app->gk_geom->geo_int.jacobgeo, &app->local); + gk_species_moment_diag_jacobgeo_div(app, &species->lte.moms, species->lte.moms.marr, species->lte.moms.marr); // Set a cap gkyl_array_set_offset(src->Jrate_cap, src->damping_factor*species->info.mass, species->lte.moms.marr, 2*app->basis.num_basis); @@ -174,37 +137,31 @@ gk_species_source_bgk_rhs_energy(gkyl_gyrokinetic_app *app, struct gk_species *s // Set the LTE moments for projection and project gkyl_array_set_offset(species->lte.moms.marr, 1.0, src->Jrate_mom, 2*app->basis.num_basis); gk_species_lte_from_moms(app, species, &species->lte, species->lte.moms.marr); - // Multiply the Maxwellian by Jrate. - gkyl_dg_mul_conf_phase_op_range(&app->basis, &species->basis, src->Jrate_fmax, - src->Jrate, species->lte.f_lte, &app->local, &species->local); - // Assemble the BGK-like term and add it to rhs. - gkyl_bgk_collisions_advance(src->bgk_op, &app->local, &species->local, - src->rate, src->Jrate_fmax, fin, src->implicit_step, src->dt_implicit, rhs, species->cflrate); - app->stat.species_source_bgk_tm += gkyl_time_diff_now_sec(wst); + // Add this Maxwellian onto the average of Maxwellians. + gkyl_array_accumulate_range(out, 1.0/3.0, species->lte.f_lte, &species->local); - // Temporaty hack, write out the M2 moment - //gkyl_array_accumulate(src->Jrate_fmax, -1.0/src->coupling_time, fin); // Testing a factor of 2 - //for (int m=0; minfo.num_diag_moments; ++m) { - // gk_species_moment_calc(&species->moms[m], species->local, app->local, src->Jrate_fmax); - // gkyl_dg_div_op_range(species->moms[m].mem_geo, app->basis, 0, species->moms[m].marr, 0, species->moms[m].marr, 0, app->gk_geom->geo_int.jacobgeo, &app->local); - // if (app->use_gpu) - // gkyl_array_copy(species->moms[m].marr_host, species->moms[m].marr); - // const char *fmt = "%s-%s_%s_hsource.gkyl"; - // int sz = gkyl_calc_strlen(fmt, app->name, species->info.name, gkyl_distribution_moments_strs[species->info.diag_moments[m]]); - // char fileNm[sz+1]; // ensures no buffer overflow - // snprintf(fileNm, sizeof fileNm, fmt, app->name, species->info.name, gkyl_distribution_moments_strs[species->info.diag_moments[m]]); - // gkyl_comm_array_write(app->comm, &app->grid, &app->local, 0, species->moms[m].marr_host, fileNm); - //} + app->stat.species_source_bgk_tm += gkyl_time_diff_now_sec(wst); } static void gk_species_source_bgk_rhs_external_enabled(gkyl_gyrokinetic_app *app, struct gk_species *species, struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *rhs) { - gk_species_source_bgk_rhs_density(app, species, src, fin, rhs); - gk_species_source_bgk_rhs_momentum(app, species, src, fin, rhs); - gk_species_source_bgk_rhs_energy(app, species, src, fin, rhs); + gkyl_array_clear(src->Jrate_fmax, 0.0); + + gks_src_bgk_rhs_accumulate_density_maxwellian(app, species, src, fin, src->Jrate_fmax); + gks_src_bgk_rhs_accumulate_momentum_maxwellian(app, species, src, fin, src->Jrate_fmax); + gks_src_bgk_rhs_accumulate_energy_maxwellian(app, species, src, fin, src->Jrate_fmax); + + // Multiply the Maxwellian by Jrate. + gkyl_dg_mul_conf_phase_op_range(&app->basis, &species->basis, species->lte.f_lte, + src->Jrate, src->Jrate_fmax, &app->local, &species->local); + // Assemble the BGK-like term and add it to rhs. + gkyl_array_clear(src->Jrate_fmax, 0.0); + gkyl_bgk_collisions_advance(src->bgk_op, &app->local, &species->local, src->rate, species->lte.f_lte, + fin, src->implicit_step, src->dt_implicit, src->Jrate_fmax, species->cflrate); + gkyl_array_accumulate(rhs, 1.0, src->Jrate_fmax); } static void @@ -270,6 +227,80 @@ gk_species_source_bgk_write_diags_external_enabled(gkyl_gyrokinetic_app* app, st gkyl_msgpack_data_release(mt); } +static void +gks_src_bgk_calc_integrated_diags_disabled(gkyl_gyrokinetic_app* app, + struct gk_species *gks, struct gk_source_bgk *src, double tm) +{ + // Do nothing. +} + +static void +gks_src_bgk_calc_integrated_diags_enabled(gkyl_gyrokinetic_app* app, + struct gk_species *gks, struct gk_source_bgk *src, double tm) +{ + struct timespec wst = gkyl_wall_clock(); + + int num_mom = src->integ_mom_op.num_mom; + double avals_global[num_mom]; + + // Compute integrated moments of source term. + gk_species_moment_calc(&src->integ_moms, gks->local, app->local, src->Jrate_fmax); + app->stat.n_mom += 1; + + // Reduce (sum) over whole domain, append to diagnostics. + gkyl_array_reduce_range(src->red_integ_diag, src->integ_moms.marr, GKYL_SUM, &app->local); + gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_SUM, num_mom, + src->red_integ_diag, src->red_integ_diag_global); + if (app->use_gpu) { + gkyl_cu_memcpy(avals_global, src->red_integ_diag_global, sizeof(double[num_mom]), GKYL_CU_MEMCPY_D2H); + } + else { + memcpy(avals_global, src->red_integ_diag_global, sizeof(double[num_mom])); + } + gkyl_dynvec_append(src->integ_diag, tm, avals_global); + + app->stat.species_diag_calc_tm += gkyl_time_diff_now_sec(wst); + app->stat.n_diag += 1; +} + +static void +gks_src_bgk_write_integrated_diags_disabled(gkyl_gyrokinetic_app *app, + struct gk_species *gks, struct gk_source_bgk *src) +{ + // Empty. +} + +static void +gks_src_bgk_write_integrated_diags_enabled(gkyl_gyrokinetic_app *app, + struct gk_species *gks, struct gk_source_bgk *src) +{ + struct timespec wst = gkyl_wall_clock(); + + int rank; + gkyl_comm_get_rank(app->comm, &rank); + + if (rank == 0) { + // Write integrated diagnostic moments. + const char *fmt = "%s-%s_source_bgk_integrated_%s.gkyl"; + int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name, gkyl_distribution_moments_strs[GKYL_F_MOMENT_M0M1M2]); + char fileNm[sz+1]; // ensures no buffer overflow + snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name, gkyl_distribution_moments_strs[GKYL_F_MOMENT_M0M1M2]); + + if (src->is_first_integ_write_call) { + gkyl_dynvec_write(src->integ_diag, fileNm); + src->is_first_integ_write_call = false; + src->integ_diag_file_exists = true; + } + else { + gkyl_dynvec_awrite(src->integ_diag, fileNm); + } + } + gkyl_dynvec_clear(src->integ_diag); + app->stat.n_diag_io += 1; + + app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wst); +} + static void gk_source_bgk_write_conf_array(gkyl_gyrokinetic_app* app, struct gk_species *gks, struct gk_source_bgk *src, int frame, double stime, char* file_suffix, @@ -321,6 +352,8 @@ gk_species_source_bgk_init(struct gkyl_gyrokinetic_app *app, struct gk_species * src->write_diags_func = gk_species_source_bgk_write_diags_disabled; src->rhs_func = gk_species_source_bgk_rhs_disabled; + src->calc_integrated_diags_func = gks_src_bgk_calc_integrated_diags_disabled; + src->write_integrated_diags_func = gks_src_bgk_write_integrated_diags_disabled; if (src->source_bgk_id == GKYL_SOURCE_BGK_DEFAULT) { int vdim_phys = gks->info.vdim == 1? 1 : 3; @@ -332,6 +365,10 @@ gk_species_source_bgk_init(struct gkyl_gyrokinetic_app *app, struct gk_species * : gkyl_array_acquire(src->rate); gkyl_proj_on_basis *proj_rate = gkyl_proj_on_basis_new(&app->grid, &app->basis, app->poly_order+1, 1, gks->info.source_bgk.rate_profile, gks->info.source_bgk.rate_profile_ctx); + + // Divide rate by 3 so we can group density, momentum and energy terms. + gkyl_array_scale(rate_host, 1.0/3.0, &app->local); + gkyl_proj_on_basis_advance(proj_rate, 0.0, &app->local, rate_host); gkyl_array_copy(src->rate, rate_host); gkyl_proj_on_basis_release(proj_rate); @@ -418,10 +455,26 @@ gk_species_source_bgk_init(struct gkyl_gyrokinetic_app *app, struct gk_species * src->implicit_step = false; src->dt_implicit = 1e9; + // Integrated moments of the operator. + gk_species_moment_init(app, gks, &src->integ_mom_op, GKYL_F_MOMENT_M0M1M2, true); + + if (app->use_gpu) { + src->red_integ_diag = gkyl_cu_malloc(sizeof(double[src->integ_moms.num_mom])); + src->red_integ_diag_global = gkyl_cu_malloc(sizeof(double[src->integ_moms.num_mom])); + } + else { + src->red_integ_diag = gkyl_malloc(sizeof(double[src->integ_moms.num_mom])); + src->red_integ_diag_global = gkyl_malloc(sizeof(double[src->integ_moms.num_mom])); + } + src->integ_diag = gkyl_dynvec_new(GKYL_DOUBLE, src->integ_moms.num_mom); + src->is_first_integ_write_call = true; + // Methods chosen at runtime. src->rhs_func = gk_species_source_bgk_rhs_external_enabled; if (src->write_diagnostics) { src->write_diags_func = gk_species_source_bgk_write_diags_external_enabled; + src->calc_integrated_diags_func = gks_src_bgk_calc_integrated_diags_enabled; + src->write_integrated_diags_func = gks_src_bgk_write_integrated_diags_enabled; } } } @@ -484,6 +537,16 @@ gk_species_source_bgk_release(const struct gkyl_gyrokinetic_app *app, const stru } gkyl_bgk_collisions_release(src->bgk_op); + gk_species_moment_release(app, &src->integ_mom_op); + gkyl_dynvec_release(src->integ_diag); + if (app->use_gpu) { + gkyl_cu_free(src->red_integ_diag); + gkyl_cu_free(src->red_integ_diag_global); + } + else { + gkyl_free(src->red_integ_diag); + gkyl_free(src->red_integ_diag_global); + } } } diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 58bf6c438..c1a022847 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -852,13 +852,20 @@ struct gk_source_bgk { struct gkyl_bgk_collisions *bgk_op; // BGK operator. bool implicit_step; // Whether or not to take an implcit BGK step. double dt_implicit; // Timestep used by the implicit collisions. + struct gk_species_moment integ_mom_op; // Integrated moments. gkyl_dynvec vtsq_amp_diag; // Stores vtsq_amplitude for diagnostics. + double *red_integ_diag, *red_integ_diag_global; // Reduced integrated moments. + gkyl_dynvec integ_diag; // Integrated moments of the source. bool is_first_diag_dynvec_write_call; // Whether dynvec is being written for the first time. // Methods chosen at runtime. void (*rhs_func)(gkyl_gyrokinetic_app *app, struct gk_species *species, struct gk_source_bgk *src, const struct gkyl_array *fin, struct gkyl_array *rhs); void (*write_diags_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, struct gk_source_bgk *src, double tm, int frame); + void (*calc_integrated_diags_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, + struct gk_source_bgk *src, double tm); + void (*write_integrated_diags_func)(gkyl_gyrokinetic_app *app, struct gk_species *gks, + struct gk_source_bgk *src); }; struct gk_positivity { @@ -2967,6 +2974,27 @@ void gk_species_source_bgk_rhs(gkyl_gyrokinetic_app *app, struct gk_species *spe void gk_species_source_bgk_write_diags(gkyl_gyrokinetic_app* app, struct gk_species *gks, struct gk_source_bgk *src, double tm, int frame); +/** + * Calculate integrated diagnostics of the BGK source operator. + * + * @param app Gyrokinetic app object. + * @param gks Species object. + * @param src BGK source object. + * @param tm Current simulation time. + */ +void gk_species_source_bgk_calc_integrated_diags(gkyl_gyrokinetic_app* app, struct gk_species *gks, + struct gk_source_bgk *src, double tm); + +/** + * Write integrated diagnostics of the BGK source operator. + * + * @param app Gyrokinetic app object. + * @param gks Species object. + * @param src BGK source object. + */ +void gk_species_source_bgk_write_integrated_diags(gkyl_gyrokinetic_app *app, + struct gk_source_bgk *src, struct gk_positivity *pos); + /** * Release species source_bgk object. * From f5a1f76974dadca96737f5334bde0d178b7d2290 Mon Sep 17 00:00:00 2001 From: manauref Date: Mon, 16 Mar 2026 12:00:57 -0700 Subject: [PATCH 2/3] Fix some names and plug EIRENE coupling integrated diags to MB solver. --- gyrokinetic/apps/gk_species_source_bgk.c | 37 ++++++++++++++++-------- gyrokinetic/apps/gkyl_gyrokinetic.h | 18 ++++++++++++ gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 4 +-- gyrokinetic/apps/gyrokinetic.c | 18 +++++++++++- gyrokinetic/apps/gyrokinetic_multib.c | 21 ++++++++++++++ 5 files changed, 83 insertions(+), 15 deletions(-) diff --git a/gyrokinetic/apps/gk_species_source_bgk.c b/gyrokinetic/apps/gk_species_source_bgk.c index 2312095a4..2fc0d69ec 100644 --- a/gyrokinetic/apps/gk_species_source_bgk.c +++ b/gyrokinetic/apps/gk_species_source_bgk.c @@ -244,11 +244,11 @@ gks_src_bgk_calc_integrated_diags_enabled(gkyl_gyrokinetic_app* app, double avals_global[num_mom]; // Compute integrated moments of source term. - gk_species_moment_calc(&src->integ_moms, gks->local, app->local, src->Jrate_fmax); + gk_species_moment_calc(&src->integ_mom_op, gks->local, app->local, src->Jrate_fmax); app->stat.n_mom += 1; // Reduce (sum) over whole domain, append to diagnostics. - gkyl_array_reduce_range(src->red_integ_diag, src->integ_moms.marr, GKYL_SUM, &app->local); + gkyl_array_reduce_range(src->red_integ_diag, src->integ_mom_op.marr, GKYL_SUM, &app->local); gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_SUM, num_mom, src->red_integ_diag, src->red_integ_diag_global); if (app->use_gpu) { @@ -286,10 +286,9 @@ gks_src_bgk_write_integrated_diags_enabled(gkyl_gyrokinetic_app *app, char fileNm[sz+1]; // ensures no buffer overflow snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name, gkyl_distribution_moments_strs[GKYL_F_MOMENT_M0M1M2]); - if (src->is_first_integ_write_call) { + if (src->is_first_diag_dynvec_write_call) { gkyl_dynvec_write(src->integ_diag, fileNm); - src->is_first_integ_write_call = false; - src->integ_diag_file_exists = true; + src->is_first_diag_dynvec_write_call = false; } else { gkyl_dynvec_awrite(src->integ_diag, fileNm); @@ -367,7 +366,7 @@ gk_species_source_bgk_init(struct gkyl_gyrokinetic_app *app, struct gk_species * app->poly_order+1, 1, gks->info.source_bgk.rate_profile, gks->info.source_bgk.rate_profile_ctx); // Divide rate by 3 so we can group density, momentum and energy terms. - gkyl_array_scale(rate_host, 1.0/3.0, &app->local); + gkyl_array_scale_range(rate_host, 1.0/3.0, &app->local); gkyl_proj_on_basis_advance(proj_rate, 0.0, &app->local, rate_host); gkyl_array_copy(src->rate, rate_host); @@ -459,15 +458,15 @@ gk_species_source_bgk_init(struct gkyl_gyrokinetic_app *app, struct gk_species * gk_species_moment_init(app, gks, &src->integ_mom_op, GKYL_F_MOMENT_M0M1M2, true); if (app->use_gpu) { - src->red_integ_diag = gkyl_cu_malloc(sizeof(double[src->integ_moms.num_mom])); - src->red_integ_diag_global = gkyl_cu_malloc(sizeof(double[src->integ_moms.num_mom])); + src->red_integ_diag = gkyl_cu_malloc(sizeof(double[src->integ_mom_op.num_mom])); + src->red_integ_diag_global = gkyl_cu_malloc(sizeof(double[src->integ_mom_op.num_mom])); } else { - src->red_integ_diag = gkyl_malloc(sizeof(double[src->integ_moms.num_mom])); - src->red_integ_diag_global = gkyl_malloc(sizeof(double[src->integ_moms.num_mom])); + src->red_integ_diag = gkyl_malloc(sizeof(double[src->integ_mom_op.num_mom])); + src->red_integ_diag_global = gkyl_malloc(sizeof(double[src->integ_mom_op.num_mom])); } - src->integ_diag = gkyl_dynvec_new(GKYL_DOUBLE, src->integ_moms.num_mom); - src->is_first_integ_write_call = true; + src->integ_diag = gkyl_dynvec_new(GKYL_DOUBLE, src->integ_mom_op.num_mom); + src->is_first_diag_dynvec_write_call = true; // Methods chosen at runtime. src->rhs_func = gk_species_source_bgk_rhs_external_enabled; @@ -493,6 +492,20 @@ gk_species_source_bgk_write_diags(gkyl_gyrokinetic_app* app, struct gk_species * src->write_diags_func(app, gks, src, tm, frame); } +void +gk_species_source_bgk_calc_integrated_diags(gkyl_gyrokinetic_app* app, struct gk_species *gks, + struct gk_source_bgk *src, double tm) +{ + src->calc_integrated_diags_func(app, gks, src, tm); +} + +void +gk_species_source_bgk_write_integrated_diags(gkyl_gyrokinetic_app *app, struct gk_species *gks, + struct gk_source_bgk *src) +{ + src->write_integrated_diags_func(app, gks, src); +} + void gk_species_source_bgk_release(const struct gkyl_gyrokinetic_app *app, const struct gk_source_bgk *src) { diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h index 836b448d3..b43984f91 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic.h @@ -937,6 +937,15 @@ void gkyl_gyrokinetic_app_calc_species_source_integrated_mom(gkyl_gyrokinetic_ap */ void gkyl_gyrokinetic_app_calc_neut_species_source_integrated_mom(gkyl_gyrokinetic_app* app, int sidx, double tm); +/** + * Calculate integrated diagnostic moments for a plasma species BGK source. + * + * @param app App object. + * @param sidx Index of species whose integrated moments to write. + * @param tm Time at which integrated diagnostics are to be computed + */ +void gkyl_gyrokinetic_app_calc_species_source_bgk_integrated_diagnostics(gkyl_gyrokinetic_app* app, int sidx, double tm); + /** * Write integrated diagnostic moments for charged species source to file. Integrated * moments are appended to the same file. @@ -946,6 +955,15 @@ void gkyl_gyrokinetic_app_calc_neut_species_source_integrated_mom(gkyl_gyrokinet */ void gkyl_gyrokinetic_app_write_species_source_integrated_mom(gkyl_gyrokinetic_app *app, int sidx); +/** + * Write integrated diagnostic moments for charged species BGK source to file. Integrated + * moments are appended to the same file. + * + * @param app App object. + * @param sidx Index of species whose source integrated moments to write. + */ +void gkyl_gyrokinetic_app_write_species_source_bgk_integrated_diagnostics(gkyl_gyrokinetic_app *app, int sidx); + /** * Write integrated diagnostic moments for neutral species source to file. Integrated * moments are appended to the same file. diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index c1a022847..c44817b0f 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -2992,8 +2992,8 @@ void gk_species_source_bgk_calc_integrated_diags(gkyl_gyrokinetic_app* app, stru * @param gks Species object. * @param src BGK source object. */ -void gk_species_source_bgk_write_integrated_diags(gkyl_gyrokinetic_app *app, - struct gk_source_bgk *src, struct gk_positivity *pos); +void gk_species_source_bgk_write_integrated_diags(gkyl_gyrokinetic_app *app, struct gk_species *gks, + struct gk_source_bgk *src); /** * Release species source_bgk object. diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c index 499c75533..58dfe687a 100644 --- a/gyrokinetic/apps/gyrokinetic.c +++ b/gyrokinetic/apps/gyrokinetic.c @@ -1467,7 +1467,7 @@ gkyl_gyrokinetic_app_write_species_damping(gkyl_gyrokinetic_app* app, int sidx, } // -// ............. df/dt multiplier utputs ............... // +// ............. df/dt multiplier outputs ............... // // void gkyl_gyrokinetic_app_write_species_fdot_multiplier(gkyl_gyrokinetic_app* app, int sidx, double tm, int frame) @@ -1486,6 +1486,20 @@ gkyl_gyrokinetic_app_write_species_source_bgk_diagnostics(gkyl_gyrokinetic_app* gk_species_source_bgk_write_diags(app, gks, &gks->bgk_src, tm, frame); } +void +gkyl_gyrokinetic_app_calc_species_source_bgk_integrated_diagnostics(gkyl_gyrokinetic_app* app, int sidx, double tm) +{ + struct gk_species *gks = &app->species[sidx]; + gk_species_source_bgk_calc_integrated_diags(app, gks, &gks->bgk_src, tm); +} + +void +gkyl_gyrokinetic_app_write_species_source_bgk_integrated_diagnostics(gkyl_gyrokinetic_app *app, int sidx) +{ + struct gk_species *gks = &app->species[sidx]; + gk_species_source_bgk_write_integrated_diags(app, gks, &gks->bgk_src); +} + // // ............. EIRENE outputs ............... // // @@ -1752,6 +1766,7 @@ gkyl_gyrokinetic_app_calc_integrated_mom(gkyl_gyrokinetic_app* app, double tm) gkyl_gyrokinetic_app_calc_species_rad_integrated_mom(app, i, tm); gkyl_gyrokinetic_app_calc_species_boundary_flux_integrated_mom(app, i, tm); gkyl_gyrokinetic_app_calc_species_source_integrated_mom(app, i, tm); + gkyl_gyrokinetic_app_calc_species_source_bgk_integrated_diagnostics(app, i, tm); gkyl_gyrokinetic_app_calc_species_positivity_integrated_diagnostics(app, i, tm); } @@ -1777,6 +1792,7 @@ gkyl_gyrokinetic_app_write_integrated_mom(gkyl_gyrokinetic_app *app) for (int i=0; inum_species; ++i) { gkyl_gyrokinetic_app_write_species_integrated_mom(app, i); gkyl_gyrokinetic_app_write_species_source_integrated_mom(app, i); + gkyl_gyrokinetic_app_write_species_source_bgk_integrated_diagnostics(app, i); gkyl_gyrokinetic_app_write_species_positivity_integrated_diagnostics(app, i); gkyl_gyrokinetic_app_write_species_lte_max_corr_status(app, i); gkyl_gyrokinetic_app_write_species_rad_integrated_mom(app, i); diff --git a/gyrokinetic/apps/gyrokinetic_multib.c b/gyrokinetic/apps/gyrokinetic_multib.c index 9e722bc43..b4e9c0199 100644 --- a/gyrokinetic/apps/gyrokinetic_multib.c +++ b/gyrokinetic/apps/gyrokinetic_multib.c @@ -1394,6 +1394,25 @@ gkyl_gyrokinetic_multib_app_write_neut_species_source_integrated_mom(gkyl_gyroki } } +// +// ............. BGK Source outputs ............... // +// +void +gkyl_gyrokinetic_multib_app_calc_species_source_bgk_integrated_diagnostics(gkyl_gyrokinetic_multib_app* app, int sidx, double tm) +{ + for (int b=0; bnum_local_blocks; ++b) { + gkyl_gyrokinetic_app_calc_species_source_bgk_integrated_diagnostics(app->singleb_apps[b], sidx, tm); + } +} + +void +gkyl_gyrokinetic_multib_app_write_species_source_bgk_integrated_diagnostics(gkyl_gyrokinetic_multib_app *app, int sidx) +{ + for (int b=0; bnum_local_blocks; ++b) { + gkyl_gyrokinetic_app_write_species_source_bgk_integrated_diagnostics(app->singleb_apps[b], sidx); + } +} + // // ............. LTE outputs ............... // // @@ -1530,6 +1549,7 @@ gkyl_gyrokinetic_multib_app_calc_integrated_mom(gkyl_gyrokinetic_multib_app* app for (int i=0; inum_species; ++i) { gkyl_gyrokinetic_multib_app_calc_species_integrated_mom(app, i, tm); gkyl_gyrokinetic_multib_app_calc_species_source_integrated_mom(app, i, tm); + gkyl_gyrokinetic_multib_app_calc_species_source_bgk_integrated_diagnostics(app, i, tm); gkyl_gyrokinetic_multib_app_calc_species_rad_integrated_mom(app, i, tm); gkyl_gyrokinetic_multib_app_calc_species_boundary_flux_integrated_mom(app, i, tm); } @@ -1546,6 +1566,7 @@ gkyl_gyrokinetic_multib_app_write_integrated_mom(gkyl_gyrokinetic_multib_app *ap for (int i=0; inum_species; ++i) { gkyl_gyrokinetic_multib_app_write_species_integrated_mom(app, i); gkyl_gyrokinetic_multib_app_write_species_source_integrated_mom(app, i); + gkyl_gyrokinetic_multib_app_write_species_source_bgk_integrated_diagnostics(app, i); gkyl_gyrokinetic_multib_app_write_species_lte_max_corr_status(app, i); gkyl_gyrokinetic_multib_app_write_species_rad_integrated_mom(app, i); gkyl_gyrokinetic_multib_app_write_species_boundary_flux_integrated_mom(app, i); From ebe3752475d8e4f22803c5b3fef20a2b9583ab1c Mon Sep 17 00:00:00 2001 From: Akash Shukla Date: Fri, 20 Mar 2026 11:35:40 -0700 Subject: [PATCH 3/3] Add solps build to mkdeps --- install-deps/build-solps.sh | 196 ++++++++++++++++++++++++++++++++++++ install-deps/mkdeps.sh | 14 +++ 2 files changed, 210 insertions(+) create mode 100755 install-deps/build-solps.sh diff --git a/install-deps/build-solps.sh b/install-deps/build-solps.sh new file mode 100755 index 000000000..b5d759c3e --- /dev/null +++ b/install-deps/build-solps.sh @@ -0,0 +1,196 @@ +#!/bin/bash + +# ============================================================================== +# SOLPS-ITER Automated Installer for NERSC Perlmutter +# ============================================================================== + +# Stop on error +set -e + +echo ">>> Starting SOLPS-ITER Installation..." +cd ../../gkylsoft/ + +# Define SOLPSTOP absolute path +# We'll get the current directory and then append SOLPS-ITER to it +# This will be used in the setup script +SOLPSTOP_VAL="$(pwd)/SOLPS-ITER" + +# 1. Clone the repository +if [ ! -d "SOLPS-ITER" ]; then + echo ">>> Cloning SOLPS-ITER..." + git clone https://github.com/iterorganization/SOLPS-ITER.git +else + echo ">>> SOLPS-ITER directory already exists. Skipping clone." +fi + +cd SOLPS-ITER + +# 2. Modify .gitmodules to remove ADAS +echo ">>> Removing ADAS from .gitmodules..." +if [ -f ".gitmodules" ]; then + # Use perl to remove the ADAS submodule section completely. + # The regex is now more flexible to match different branch names or whitespace. + perl -i -0777 -pe 's/\[submodule "ADAS"\]\s+path = modules\/adas\s+url = ssh:\/\/git\@git\.iter\.org\/imex\/amns-adas\.git\s+branch = .*?\n//sg' .gitmodules +fi + +# 3. Initialize and update submodules +echo ">>> Initializing submodules..." +# If ADAS was already initialized, we might need to remove it from git config +git config --local --remove-section submodule.modules/adas 2>/dev/null || true +git config --local --remove-section submodule.ADAS 2>/dev/null || true +# Remove the directory if it exists and isn't a proper submodule anymore +if [ -d "modules/adas" ]; then + rm -rf modules/adas +fi +git submodule init || true +git submodule update + +# 4. Modify SETUP/setup.csh.NERSC.gfortran +echo ">>> Configuring SETUP/setup.csh.NERSC.gfortran..." + +setup_file="SETUP/setup.csh.NERSC.gfortran" + +# NEW_HEADER includes the corrected SOLPSTOP definition and NERSC modules +NEW_HEADER=$(cat << EOM + setenv SOLPSTOP $SOLPSTOP_VAL + +module load PrgEnv-gnu +module load cray-hdf5 +module load cray-libsci +module load cray-netcdf +module load texlive/2024 +module load cudatoolkit + +setenv GLI_HOME \${SOLPSTOP}/lib/NERSC.gfortran +setenv MSCL_ROOT \${SOLPSTOP}/lib/NERSC.gfortran +setenv GR_ROOT \${SOLPSTOP}/lib/NERSC.gfortran +setenv GR_SOURCE \${SOLPSTOP}/lib/NERSC.gfortran +setenv GKS_ROOT \${SOLPSTOP}/lib/NERSC.gfortran +setenv PATH \${PATH}:\${SOLPSTOP}/scripts +setenv LD_LIBRARY_PATH \${SOLPSTOP}/lib/NERSC.gfortran:\$LD_LIBRARY_PATH +EOM +) + +# Replace the top block up to '# Environment variables' +# Use a temporary file to safely handle the replacement with a large string +export NEW_HEADER +perl -i -0777 -pe 's/^(.*?)# Environment variables/$ENV{NEW_HEADER}\n# Environment variables/s' "$setup_file" + +# 5. Step 1: First Setup & Miniforge +echo ">>> Running first_setup and installing miniforge..." +cat << 'EOF' > step1_install.csh +#!/bin/tcsh +setenv HOST_NAME NERSC +setenv COMPILER gfortran +source SETUP/setup.csh.NERSC.gfortran +first_setup +printf "y\n" | install_miniforge +EOF + +tcsh step1_install.csh +rm step1_install.csh + +# 6. Step 2: Dependencies +echo ">>> Installing dependencies..." +# Based on instructions: +# Yes: Blas/Lapack, JSON-Fortran, NCARG, MSCL, GKS, GR +# No: NetCDF, hdf5 (these are available on Perlmutter) +# +# The expected prompt order in install_dependencies (standard for SOLPS-ITER) is: +# 1. Gnuplot (y) - but it seems it might be already installed or skipped +# 2. Blas/Lapack (y) +# 3. NetCDF (n) +# 4. JSON-Fortran (y) +# 5. NCARG (y) +# 6. MSCL (y) +# 7. GKS (y) +# 8. GR (y) +# Note: hdf5 doesn't seem to have a prompt in the script for HPC? +# Actually, I'll provide enough 'y' and 'n' to cover all bases. +INPUTS="y\ny\nn\ny\ny\ny\ny\ny\n" + +cat << EOF > step2_deps.csh +#!/bin/tcsh +setenv HOST_NAME NERSC +setenv COMPILER gfortran +source SETUP/setup.csh.NERSC.gfortran +conda activate ./lib/NERSC.gfortran/miniforge3 +setenv CC cc +printf "$INPUTS" | install_dependencies +EOF + +tcsh step2_deps.csh +rm step2_deps.csh + +# 7. Step 3: Fix JSON-Fortran Build +echo ">>> Fixing JSON-Fortran build..." +cat << 'EOF' > step3_json.csh +#!/bin/tcsh +setenv HOST_NAME NERSC +setenv COMPILER gfortran +source SETUP/setup.csh.NERSC.gfortran +conda activate solps_env + +if ( ! -d lib/NERSC.gfortran/json-fortran-9.0.2 ) then + cd lib/NERSC.gfortran + wget https://github.com/jacobwilliams/json-fortran/archive/refs/tags/9.0.2.tar.gz + tar xvf 9.0.2.tar.gz + rm 9.0.2.tar.gz + cd ../.. +endif + +cd lib/NERSC.gfortran/json-fortran-9.0.2 +if ( -d build ) then + rm -rf build +endif +mkdir build +cd build +cmake -DENABLE_IN_SOURCE_BUILDS=ON .. +make +# Also fulfill the user's manual step as it might be what they expect +cd .. +# If there is no Makefile in the root, 'make' will fail but we can ignore it +# because the build in 'build/' directory succeeded. +make || true +EOF + +tcsh step3_json.csh +rm step3_json.csh + +# 8. Step 4: Configure Carre and Compile +echo ">>> Configuring Carre and compiling SOLPS..." + +config_file="modules/Carre/config/config.NERSC.gfortran" +# Add DEFINES += -DGFORTRAN after LD_NCARG +# We'll use a more precise match that handles both ?= and = +# Check if it already exists to avoid duplicate entries +if ! grep -q "DEFINES += -DGFORTRAN" "$config_file"; then + perl -i -pe 's/(LD_NCARG\s*[\?\=]*\s*.*\n)/$1DEFINES += -DGFORTRAN\n/' "$config_file" +fi + +cat << 'EOF' > step4_compile.csh +#!/bin/tcsh +setenv HOST_NAME NERSC +setenv COMPILER gfortran +source setup.csh +conda activate solps_env + +make clean +make clean_solps + +echo ">>> Compilation starting. Piping 'q' to handle any pdflatex pauses..." +# We use 'make' instead of 'make solps' because 'make clean_solps' and 'make solps' are part of the instructions. +# 'make' in the root usually builds everything. +# Using source setup.csh instead of specific setup.csh.NERSC.gfortran because setup.csh +# handles adding scripts/ to PATH, which is where which_latex resides. +yes "q" | make solps + +echo ">>> Performing manual DivGeo build to ensure dg.exe is created..." +cd modules/DivGeo +make +EOF + +tcsh step4_compile.csh +rm step4_compile.csh + +echo ">>> Installation Complete." diff --git a/install-deps/mkdeps.sh b/install-deps/mkdeps.sh index 825f01e22..c234b77a2 100755 --- a/install-deps/mkdeps.sh +++ b/install-deps/mkdeps.sh @@ -17,6 +17,7 @@ BUILD_OPENMPI=no BUILD_LUAJIT=no BUILD_TCC=no BUILD_CUDSS=no +BUILD_SOLPS=no # by default, download as well as build packages DOWNLOAD_PKGS=yes @@ -167,6 +168,10 @@ do [ -n "$value" ] || die "Missing value in flag $key." BUILD_ADAS="$value" ;; + --build-solps) + [ -n "$value" ] || die "Missing value in flag $key." + BUILD_SOLPS="$value" + ;; *) die "Error: Unknown flag: $1" ;; @@ -261,6 +266,14 @@ build_adas() { fi } +build_solps() { + if [ "$BUILD_SOLPS" = "yes" ] + then + echo "Building SOLPS for neutral evolution" + ./build-solps.sh + fi +} + echo "Installations will be in $PREFIX" build_openmpi @@ -270,3 +283,4 @@ build_openblas build_superlu build_cudss build_adas +build_solps