Skip to content

Add loading a source from a .gkyl file #1087

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Maxwell-Rosen wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add loading a source from a .gkyl file #1087

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion core/zero/gkyl_eqn_type.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,8 @@ enum gkyl_source_id {
GKYL_NO_SOURCE = 0, // No source. This is default.
GKYL_FUNC_SOURCE, // Function source.
GKYL_PROJ_SOURCE, // Source given by projection object determined by gkyl_projection_id.
GKYL_BFLUX_SOURCE // Source which scales to boundary fluxes.
GKYL_BFLUX_SOURCE, // Source which scales to boundary fluxes.
GKYL_SOURCE_FROMFILE, // Source shape imported directly from a file.
};

// Identifiers for specific bgk source object types.
Expand Down
27 changes: 17 additions & 10 deletions gyrokinetic/apps/gk_neut_species_kinetic.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -375,10 +375,10 @@ gk_neut_species_kinetic_init_static(struct gkyl_gk *gk, struct gkyl_gyrokinetic_
}

void
gk_neut_species_kinetic_file_import_init(struct gkyl_gyrokinetic_app *app, struct gk_neut_species *s,
struct gkyl_gyrokinetic_ic_import inp)
gk_neut_species_kinetic_read_distf_from_file(struct gkyl_gyrokinetic_app *app, struct gk_neut_species *s,
struct gkyl_gyrokinetic_ic_import inp, struct gkyl_array *fout)
{
// Import initial condition from a file. Intended options include importing:
// Read a distribution function from a file onto fout. Intended options include importing:
// 1) ICs with same grid.
// 2) ICs one dimensionality lower (e.g. 2x2v for 3x2v sim).
// 3) ICs with same grid extents but different resolution (NYI).
Expand Down Expand Up @@ -538,47 +538,54 @@ gk_neut_species_kinetic_file_import_init(struct gkyl_gyrokinetic_app *app, struc
if (pdim_do == pdim-1) {
struct gkyl_translate_dim* transdim = gkyl_translate_dim_new(cdim_do,
basis_do, cdim, s->basis, -1, GKYL_NO_EDGE, app->use_gpu);
gkyl_translate_dim_advance(transdim, &local_do, &s->local, fdo, 1, s->f);
gkyl_translate_dim_advance(transdim, &local_do, &s->local, fdo, 1, fout);
gkyl_translate_dim_release(transdim);
}
else {
if (same_res) {
gkyl_array_copy(s->f, fdo);
gkyl_array_copy(fout, fdo);
}
else {
// Interpolate the donor distribution to the target grid.
struct gkyl_dg_interpolate *interp = gkyl_dg_interpolate_new(app->cdim, &s->basis,
&grid_do, &grid, &local_do, &s->local, ghost_do, app->use_gpu);
gkyl_dg_interpolate_advance(interp, fdo, s->f);
gkyl_dg_interpolate_advance(interp, fdo, fout);
gkyl_dg_interpolate_release(interp);
}
}

if (inp.type == GKYL_IC_IMPORT_AF) {
// Scale f by a conf-space factor.
// Scale fout by a conf-space factor.
gkyl_proj_on_basis *proj_conf_scale = gkyl_proj_on_basis_new(&app->grid, &app->basis,
poly_order+1, 1, inp.conf_scale, inp.conf_scale_ctx);
struct gkyl_array *xfac = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
struct gkyl_array *xfac_ho = app->use_gpu? mkarr(false, app->basis.num_basis, app->local_ext.volume)
: gkyl_array_acquire(xfac_ho);
gkyl_proj_on_basis_advance(proj_conf_scale, 0.0, &app->local, xfac_ho);
gkyl_array_copy(xfac, xfac_ho);
gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, s->f, xfac, s->f, &app->local, &s->local);
gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, fout, xfac, fout, &app->local, &s->local);
gkyl_proj_on_basis_release(proj_conf_scale);
gkyl_array_release(xfac_ho);
gkyl_array_release(xfac);
}

// Multiply f by the Jacobian.
// Multiply fout by the Jacobian.
if (scale_by_jacobgeo)
gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, s->f, app->gk_geom->geo_int.jacobgeo, s->f, &app->local, &s->local);
gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, fout, app->gk_geom->geo_int.jacobgeo, fout, &app->local, &s->local);

gkyl_rect_decomp_release(decomp_do);
gkyl_comm_release(comm_do);
gkyl_array_release(fdo);
gkyl_array_release(fdo_host);
}

void
gk_neut_species_kinetic_file_import_init(struct gkyl_gyrokinetic_app *app, struct gk_neut_species *s,
struct gkyl_gyrokinetic_ic_import inp)
{
gk_neut_species_kinetic_read_distf_from_file(app, s, inp, s->f);
}

static void
gkyl_array_move_comp(struct gkyl_array *out, int cout, struct gkyl_array *in, int cin, struct gkyl_array *tmp)
{
Expand Down
29 changes: 23 additions & 6 deletions gyrokinetic/apps/gk_neut_species_source.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,22 @@ gk_neut_species_source_init(struct gkyl_gyrokinetic_app *app, struct gk_neut_spe
src->evolve = s->info.source.evolve; // Whether the source is time dependent.

src->num_sources = s->info.source.num_sources;
for (int k=0; k<s->info.source.num_sources; k++)
gk_neut_species_projection_init(app, s, s->info.source.projection[k], &src->proj_source[k]);
if (src->source_id == GKYL_SOURCE_FROMFILE) {
// A file-imported source shape is static: no time dependence.
assert(!src->evolve);
struct gkyl_array *fdo = mkarr(app->use_gpu, s->basis.num_basis, s->local_ext.volume);
gkyl_array_clear(src->source, 0.0);
for (int k=0; k<s->info.source.num_sources; k++) {
assert(!s->info.source.source_import[k].enforce_positivity); // Not supported for sources yet.
gk_neut_species_kinetic_read_distf_from_file(app, s, s->info.source.source_import[k], fdo);
gkyl_array_accumulate(src->source, 1., fdo);
}
gkyl_array_release(fdo);
}
else {
for (int k=0; k<s->info.source.num_sources; k++)
gk_neut_species_projection_init(app, s, s->info.source.projection[k], &src->proj_source[k]);
}

// Allocate data and updaters for diagnostic moments.
src->num_diag_mom = s->info.num_diag_moments;
Expand All @@ -45,10 +59,11 @@ gk_neut_species_source_init(struct gkyl_gyrokinetic_app *app, struct gk_neut_spe
}

void
gk_neut_species_source_calc(gkyl_gyrokinetic_app *app, struct gk_neut_species *s,
gk_neut_species_source_calc(gkyl_gyrokinetic_app *app, struct gk_neut_species *s,
struct gk_source *src, struct gkyl_array *f_buffer, double tm)
{
if (src->source_id) {
// A file-imported source's shape was already read once in gk_neut_species_source_init.
if (src->source_id && src->source_id != GKYL_SOURCE_FROMFILE) {
for (int k=0; k<s->info.source.num_sources; k++) {
gk_neut_species_projection_calc(app, s, &src->proj_source[k], f_buffer, tm);
gkyl_array_accumulate(src->source, 1., f_buffer);
Expand Down Expand Up @@ -225,8 +240,10 @@ gk_neut_species_source_release(const struct gkyl_gyrokinetic_app *app, const str
if (app->use_gpu) {
gkyl_array_release(src->source_host);
}
for (int k=0; k<src->num_sources; k++) {
gk_neut_species_projection_release(app, &src->proj_source[k]);
if (src->source_id != GKYL_SOURCE_FROMFILE) {
for (int k=0; k<src->num_sources; k++) {
gk_neut_species_projection_release(app, &src->proj_source[k]);
}
}

// Release moment data.
Expand Down
41 changes: 24 additions & 17 deletions gyrokinetic/apps/gk_species.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1088,10 +1088,10 @@ gk_species_init_static(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *
}

void
gk_species_file_import_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
struct gkyl_gyrokinetic_ic_import inp)
gk_species_read_distf_from_file(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
struct gkyl_gyrokinetic_ic_import inp, struct gkyl_array *fout)
{
// Import initial condition from a file. Intended options include importing:
// Read a distribution function from a file onto fout. Intended options include importing:
// 1) ICs with same grid.
// 2) ICs one dimensionality lower (e.g. 2x2v for 3x2v sim).
// 3) ICs with same grid extents but different resolution (NYI).
Expand Down Expand Up @@ -1267,44 +1267,56 @@ gk_species_file_import_init(struct gkyl_gyrokinetic_app *app, struct gk_species
if (pdim_do == pdim-1) {
struct gkyl_translate_dim* transdim = gkyl_translate_dim_new(cdim_do,
basis_do, cdim, gks->basis, -1, GKYL_NO_EDGE, app->use_gpu);
gkyl_translate_dim_advance(transdim, &local_do, &gks->local, fdo, 1, gks->f);
gkyl_translate_dim_advance(transdim, &local_do, &gks->local, fdo, 1, fout);
gkyl_translate_dim_release(transdim);
}
else {
if (same_res) {
gkyl_array_copy(gks->f, fdo);
gkyl_array_copy(fout, fdo);
}
else {
// Interpolate the donor distribution to the target grid.
struct gkyl_dg_interpolate *interp = gkyl_dg_interpolate_new(app->cdim, &gks->basis,
&grid_do, &grid, &local_do, &gks->local, ghost_do, app->use_gpu);
gkyl_dg_interpolate_advance(interp, fdo, gks->f);
gkyl_dg_interpolate_advance(interp, fdo, fout);
gkyl_dg_interpolate_release(interp);
}
}

if (inp.type == GKYL_IC_IMPORT_AF) {
// Scale f by a conf-space factor.
// Scale fout by a conf-space factor.
gkyl_proj_on_basis *proj_conf_scale = gkyl_proj_on_basis_new(&app->grid, &app->basis,
poly_order+1, 1, inp.conf_scale, inp.conf_scale_ctx);
struct gkyl_array *xfac = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
struct gkyl_array *xfac_ho = app->use_gpu? mkarr(false, app->basis.num_basis, app->local_ext.volume)
: gkyl_array_acquire(xfac);
gkyl_proj_on_basis_advance(proj_conf_scale, 0.0, &app->local, xfac_ho);
gkyl_array_copy(xfac, xfac_ho);
gkyl_dg_mul_conf_phase_op_range(&app->basis, &gks->basis, gks->f, xfac, gks->f, &app->local, &gks->local);
gkyl_dg_mul_conf_phase_op_range(&app->basis, &gks->basis, fout, xfac, fout, &app->local, &gks->local);
gkyl_proj_on_basis_release(proj_conf_scale);
gkyl_array_release(xfac_ho);
gkyl_array_release(xfac);
}

// Multiply f by the Jacobian.
// Multiply fout by the Jacobian.
if (scale_by_jacobtot)
gkyl_dg_mul_conf_phase_op_range(&app->basis, &gks->basis, gks->f, app->gk_geom->geo_int.jacobtot, gks->f, &app->local, &gks->local);
gkyl_dg_mul_conf_phase_op_range(&app->basis, &gks->basis, fout, app->gk_geom->geo_int.jacobtot, fout, &app->local, &gks->local);

// Multiply f by the velocity space Jacobian.
// Multiply fout by the velocity space Jacobian.
if (scale_by_jacobvel)
gkyl_array_scale_by_cell(gks->f, gks->vel_map->jacobvel);
gkyl_array_scale_by_cell(fout, gks->vel_map->jacobvel);

gkyl_rect_decomp_release(decomp_do);
gkyl_comm_release(comm_do);
gkyl_array_release(fdo);
gkyl_array_release(fdo_host);
}

void
gk_species_file_import_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
struct gkyl_gyrokinetic_ic_import inp)
{
gk_species_read_distf_from_file(app, gks, inp, gks->f);

if (inp.enforce_positivity) {
// Positivity enforcing by shifting f (ps=positivity shift).
Expand All @@ -1316,11 +1328,6 @@ gk_species_file_import_init(struct gkyl_gyrokinetic_app *app, struct gk_species

gkyl_positivity_shift_gyrokinetic_release(pos_shift_op);
}

gkyl_rect_decomp_release(decomp_do);
gkyl_comm_release(comm_do);
gkyl_array_release(fdo);
gkyl_array_release(fdo_host);
}

static bool
Expand Down
29 changes: 23 additions & 6 deletions gyrokinetic/apps/gk_species_source.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -264,10 +264,11 @@ gk_species_source_write_integrated_mom_enabled(gkyl_gyrokinetic_app* app, struct


void
gk_species_source_calc(gkyl_gyrokinetic_app *app, struct gk_species *s,
gk_species_source_calc(gkyl_gyrokinetic_app *app, struct gk_species *s,
struct gk_source *src, struct gkyl_array *f_buffer, double tm)
{
if (src->source_id) {
// A file-imported source's shape was already read once in gk_species_source_init.
if (src->source_id && src->source_id != GKYL_SOURCE_FROMFILE) {
gkyl_array_clear(src->source, 0.0);
for (int k=0; k<s->info.source.num_sources; k++) {
gk_species_projection_calc(app, s, &src->proj_source[k], f_buffer, tm);
Expand Down Expand Up @@ -411,8 +412,22 @@ gk_species_source_init(struct gkyl_gyrokinetic_app *app, struct gk_species *s,
src->evolve = s->info.source.evolve || s->info.source.num_adapt_sources > 0; // Whether the source is time dependent.

src->num_sources = s->info.source.num_sources;
for (int k=0; k<s->info.source.num_sources; k++)
gk_species_projection_init(app, s, s->info.source.projection[k], &src->proj_source[k]);
if (src->source_id == GKYL_SOURCE_FROMFILE) {
// A file-imported source shape is static: no time dependence or adaptive scaling.
assert(!src->evolve);
struct gkyl_array *fdo = mkarr(app->use_gpu, s->basis.num_basis, s->local_ext.volume);
gkyl_array_clear(src->source, 0.0);
for (int k=0; k<s->info.source.num_sources; k++) {
assert(!s->info.source.source_import[k].enforce_positivity); // Not supported for sources yet.
gk_species_read_distf_from_file(app, s, s->info.source.source_import[k], fdo);
gkyl_array_accumulate(src->source, 1., fdo);
}
gkyl_array_release(fdo);
}
else {
for (int k=0; k<s->info.source.num_sources; k++)
gk_species_projection_init(app, s, s->info.source.projection[k], &src->proj_source[k]);
}

// Allocate data and updaters for diagnostic moments.
src->num_diag_mom = s->info.source.diagnostics.num_diag_moments;
Expand Down Expand Up @@ -609,8 +624,10 @@ gk_species_source_release(const struct gkyl_gyrokinetic_app *app, const struct g
if (app->use_gpu) {
gkyl_array_release(src->source_host);
}
for (int k=0; k<src->num_sources; k++) {
gk_species_projection_release(app, &src->proj_source[k]);
if (src->source_id != GKYL_SOURCE_FROMFILE) {
for (int k=0; k<src->num_sources; k++) {
gk_species_projection_release(app, &src->proj_source[k]);
}
}

// Release moment data.
Expand Down
3 changes: 3 additions & 0 deletions gyrokinetic/apps/gkyl_gyrokinetic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,9 @@ struct gkyl_gyrokinetic_source {
struct gkyl_gyrokinetic_adapt_source adapt[GKYL_MAX_SOURCES]; // Adaptive source parameters
// sources using projection routine
struct gkyl_gyrokinetic_projection projection[GKYL_MAX_SOURCES];
// Sources whose shape is imported directly from a file (used only when
// source_id == GKYL_SOURCE_FROMFILE).
struct gkyl_gyrokinetic_ic_import source_import[GKYL_MAX_SOURCES];

struct gkyl_phase_diagnostics_inp diagnostics;
};
Expand Down
25 changes: 25 additions & 0 deletions gyrokinetic/apps/gkyl_gyrokinetic_priv.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2805,6 +2805,19 @@ void gk_species_projection_calc(gkyl_gyrokinetic_app *app, struct gk_species *sp
*/
void gk_species_projection_release(const struct gkyl_gyrokinetic_app *app, const struct gk_proj *proj);

/**
* Read a distribution function from a file (optionally rescaling by
* 1/Jacobian and the velocity-space Jacobian, and translating/interpolating
* to the target grid) into aout.
*
* @param app Gyrokinetic app object.
* @param gks Species object (defines the target grid/basis/geometry).
* @param inp Input struct describing the file(s) to import.
* @param aout Output array to write the imported distribution into.
*/
void gk_species_read_distf_from_file(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
struct gkyl_gyrokinetic_ic_import inp, struct gkyl_array *aout);

/** gk_species_source API */

/**
Expand Down Expand Up @@ -3876,6 +3889,18 @@ void gk_neut_species_projection_calc(gkyl_gyrokinetic_app *app, struct gk_neut_s
*/
void gk_neut_species_projection_release(const struct gkyl_gyrokinetic_app *app, const struct gk_proj *proj);

/**
* Read a distribution function from a file (optionally rescaling by
* 1/Jacobian and translating/interpolating to the target grid) into aout.
*
* @param app Gyrokinetic app object.
* @param s Neutral species object (defines the target grid/basis/geometry).
* @param inp Input struct describing the file(s) to import.
* @param aout Output array to write the imported distribution into.
*/
void gk_neut_species_kinetic_read_distf_from_file(struct gkyl_gyrokinetic_app *app, struct gk_neut_species *s,
struct gkyl_gyrokinetic_ic_import inp, struct gkyl_array *aout);

/** gk_neut_species_source API */

/**
Expand Down
Loading
Loading