From 54b61e7ec6f8077522165095c31e568dfadb0290 Mon Sep 17 00:00:00 2001 From: Maxwell-Rosen Date: Tue, 25 Aug 2026 09:55:32 -0400 Subject: [PATCH] Add omegaH_dt diagnostic support to gyrokinetic species - Introduced omegaH_dt_diagnostic flag in gkyl_gyrokinetic_species and gkyl_gyrokinetic_multib_species structs. - Implemented logic to compute and write omegaH stable time step diagnostics in gk_species_calc_integrated_mom_dynamic and gk_species_write_integrated_mom_dynamic functions. - Updated initialization and release functions to handle omegaH_dt dynvec. - Enhanced Lua binding to support omegaH_dt_diagnostic configuration. --- gyrokinetic/apps/gk_species.c | 46 ++++++++++++++++++++++ gyrokinetic/apps/gkyl_gyrokinetic.h | 1 + gyrokinetic/apps/gkyl_gyrokinetic_multib.h | 1 + gyrokinetic/apps/gkyl_gyrokinetic_priv.h | 3 ++ gyrokinetic/apps/gyrokinetic.c | 2 + gyrokinetic/apps/gyrokinetic_lw.c | 1 + gyrokinetic/apps/gyrokinetic_multib.c | 1 + 7 files changed, 55 insertions(+) diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c index 2b93be63cf..12ebbd306c 100644 --- a/gyrokinetic/apps/gk_species.c +++ b/gyrokinetic/apps/gk_species.c @@ -494,6 +494,14 @@ gk_species_calc_integrated_mom_dynamic(gkyl_gyrokinetic_app* app, struct gk_spec } gkyl_dynvec_append(gks->integ_diag, tm, avals_global); + if (gks->info.omegaH_dt_diagnostic) { + double omegaH_dt_local = gk_species_omegaH_dt(app, gks, gks->f); + double omegaH_dt_global = DBL_MAX; + gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MIN, 1, + &omegaH_dt_local, &omegaH_dt_global); + gkyl_dynvec_append(gks->omegaH_dt, tm, &omegaH_dt_global); + } + if (gks->info.time_rate_diagnostics) { // Reduce (sum) over whole domain, append to diagnostics. gkyl_array_accumulate(gks->fdot_mom_new, -1.0, gks->fdot_mom_old); @@ -551,6 +559,37 @@ gk_species_write_integrated_mom_dynamic(gkyl_gyrokinetic_app *app, struct gk_spe } gkyl_dynvec_clear(gks->integ_diag); app->stat.n_diag_io += 1; + + if (gks->info.omegaH_dt_diagnostic) { + if (rank == 0) { + const char *fmt = "%s-%s_omegaH_dt.gkyl"; + int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name); + char fileNm[sz+1]; + snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name); + + if (gks->is_first_omegaH_dt_write_call) { + struct gkyl_msgpack_map_elem io_meta_omegaH_dt[] = { + { .key = "Description", .elem_type = GKYL_MP_STRING, + .cval = "Stable time step for the omega_H mode for this charged species." } + }; + int io_meta_len[] = {gks->io_meta_basic_len, app->gk_geom->io_meta_basic_len, 1}; + const struct gkyl_msgpack_map_elem* io_meta[] = { + gks->io_meta_basic, app->gk_geom->io_meta_basic, io_meta_omegaH_dt + }; + struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union( + sizeof(io_meta_len)/sizeof(int), io_meta_len, io_meta); + + gkyl_dynvec_write_wmeta(gks->omegaH_dt, fileNm, mt); + gks->is_first_omegaH_dt_write_call = false; + gkyl_msgpack_data_release(mt); + } + else { + gkyl_dynvec_awrite(gks->omegaH_dt, fileNm); + } + } + gkyl_dynvec_clear(gks->omegaH_dt); + app->stat.n_diag_io += 1; + } if (gks->info.time_rate_diagnostics) { if (rank == 0) { @@ -722,6 +761,8 @@ gk_species_release_dynamic(const gkyl_gyrokinetic_app* app, const struct gk_spec // Release integrated diag memory. gkyl_dynvec_release(s->integ_diag); + if (s->info.omegaH_dt_diagnostic) + gkyl_dynvec_release(s->omegaH_dt); if (app->use_gpu) { gkyl_cu_free(s->red_integ_diag); gkyl_cu_free(s->red_integ_diag_global); @@ -805,6 +846,11 @@ gk_species_init_dynamic(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app gks->integ_diag = gkyl_dynvec_new(GKYL_DOUBLE, gks->integ_moms.num_mom); gks->is_first_integ_write_call = true; + if (gks->info.omegaH_dt_diagnostic) { + gks->omegaH_dt = gkyl_dynvec_new(GKYL_DOUBLE, 1); + gks->is_first_omegaH_dt_write_call = true; + } + // Allocate dynamic-vector to store Delta f integrated moments. if (gks->info.time_rate_diagnostics) { gks->fdot_mom_old = mkarr(app->use_gpu, gks->integ_moms.marr->ncomp, gks->integ_moms.marr->size); diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h index e04e1510ba..3679ce73cf 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic.h @@ -452,6 +452,7 @@ struct gkyl_gyrokinetic_species { enum gkyl_distribution_moments integrated_diag_moments[12]; // List of integrated diagnostic moments. bool time_rate_diagnostics; // Whether to ouput df/dt diagnostics. bool write_omega_cfl; // Whether to ouput dt diagnostic for the CFL constraint. + bool omegaH_dt_diagnostic; // Whether to output the omega_H stable time step as a dynvector. struct gkyl_gyrokinetic_collisionless collisionless; // Collisionless terms. diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_multib.h b/gyrokinetic/apps/gkyl_gyrokinetic_multib.h index e47288d6f1..7adbf7156a 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_multib.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_multib.h @@ -33,6 +33,7 @@ struct gkyl_gyrokinetic_multib_species { int num_integrated_diag_moments; // Number of integrated diagnostic moments. enum gkyl_distribution_moments integrated_diag_moments[12]; // List of integrated diagnostic moments. bool time_rate_diagnostics; // Whether to ouput df/dt diagnostics. + bool omegaH_dt_diagnostic; // Whether to output the omega_H stable time step as a dynvector. struct gkyl_phase_diagnostics_inp boundary_flux_diagnostics; diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h index 5615160973..5ee193bcd4 100644 --- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h +++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h @@ -1047,6 +1047,9 @@ struct gk_species { gkyl_dynvec integ_diag; // Integrated moments reduced across grid bool is_first_integ_write_call; // Whether dynvec is being written for the first time. + gkyl_dynvec omegaH_dt; // omega_H stable time step reduced across MPI ranks. + bool is_first_omegaH_dt_write_call; // Whether omega_H dt dynvec is being written for the first time. + struct gkyl_array *fdot_mom_old, *fdot_mom_new; // Moments of f_old and f_new. gkyl_dynvec fdot_integ_diag; // Integrated moments of Delta f=f_new - f_old.. bool is_first_fdot_integ_write_call; // Whether dynvec is being written for the first time. diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c index 526aeebc67..9991850e4d 100644 --- a/gyrokinetic/apps/gyrokinetic.c +++ b/gyrokinetic/apps/gyrokinetic.c @@ -3387,6 +3387,8 @@ gkyl_gyrokinetic_app_from_frame_species(gkyl_gyrokinetic_app *app, int sidx, int // Append to existing integrated diagnostics. app->is_first_dt_write_call = false; gk_s->is_first_integ_write_call = false; + if (gk_s->info.omegaH_dt_diagnostic) + gk_s->is_first_omegaH_dt_write_call = false; gk_s->is_first_L2norm_write_call = false; for (int b=0; bbflux.num_boundaries; ++b) gk_s->bflux.is_first_intmom_write_call[b] = false; diff --git a/gyrokinetic/apps/gyrokinetic_lw.c b/gyrokinetic/apps/gyrokinetic_lw.c index 7a1a1500a7..f3f3aaf21b 100644 --- a/gyrokinetic/apps/gyrokinetic_lw.c +++ b/gyrokinetic/apps/gyrokinetic_lw.c @@ -356,6 +356,7 @@ gyrokinetic_species_lw_new(lua_State *L) gk_species.charge = glua_tbl_get_number(L, "charge", 0.0); gk_species.mass = glua_tbl_get_number(L, "mass", 1.0); gk_species.polarization_density = glua_tbl_get_number(L, "polarizationDensity", 0.0); + gk_species.omegaH_dt_diagnostic = glua_tbl_get_bool(L, "omegaHDtDiagnostic", false); with_lua_tbl_tbl(L, "cells") { vdim = glua_objlen(L); diff --git a/gyrokinetic/apps/gyrokinetic_multib.c b/gyrokinetic/apps/gyrokinetic_multib.c index de0352dee2..87b1c17bcb 100644 --- a/gyrokinetic/apps/gyrokinetic_multib.c +++ b/gyrokinetic/apps/gyrokinetic_multib.c @@ -176,6 +176,7 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid, species_inp.integrated_diag_moments[n] = sp->integrated_diag_moments[n]; } species_inp.time_rate_diagnostics = sp->time_rate_diagnostics; + species_inp.omegaH_dt_diagnostic = sp->omegaH_dt_diagnostic; species_inp.boundary_flux_diagnostics = sp->boundary_flux_diagnostics; // Choose proper block-specific species input.