Skip to content

Commit

Permalink
Reduce duplication in species initialisation
Browse files Browse the repository at this point in the history 
Closes #240

See merge request gysela-developpers/gyselalibxx!515

--------------------------------------------

Co-authored-by: Emily Bourne <[email protected]>
  • Loading branch information
@EmilyBourne
EmilyBourne committed Jun 14, 2024
1 parent 039f809 commit 8825ed4
Show file tree
Hide file tree
Showing 10 changed files with 209 additions and 256 deletions.
31 changes: 3 additions & 28 deletions simulations/geometryXVx/bump_on_tail/bumpontail_fem_uniform.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
#include "restartinitialization.hpp"
#include "singlemodeperturbinitialization.hpp"
#include "species_info.hpp"
#include "species_init.hpp"
#include "spline_interpolator.hpp"
#include "splitvlasovsolver.hpp"

Expand Down Expand Up @@ -70,8 +71,7 @@ int main(int argc, char** argv)
SplineInterpPointsVx>(conf_voicexx, "vx");
IDomainXVx const meshXVx(mesh_x, mesh_vx);

IVectSp const nb_kinspecies(PCpp_len(conf_voicexx, ".SpeciesInfo"));
IDomainSp const dom_kinsp(IndexSp(0), nb_kinspecies);
IDomainSp const dom_kinsp = init_species(conf_voicexx);

IDomainSpXVx const meshSpXVx(dom_kinsp, meshXVx);
IDomainSpVx const meshSpVx(dom_kinsp, mesh_vx);
Expand All @@ -81,48 +81,23 @@ int main(int argc, char** argv)
SplineVxBuilder const builder_vx(meshXVx);
SplineVxBuilder_1d const builder_vx_poisson(mesh_vx);

host_t<FieldSp<int>> kinetic_charges(dom_kinsp);
host_t<DFieldSp> masses(dom_kinsp);
host_t<DFieldSp> epsilon_bot(dom_kinsp);
host_t<DFieldSp> temperature_bot(dom_kinsp);
host_t<DFieldSp> mean_velocity_bot(dom_kinsp);
host_t<DFieldSp> init_perturb_amplitude(dom_kinsp);
host_t<FieldSp<int>> init_perturb_mode(dom_kinsp);
int nb_elec_adiabspecies = 1;
int nb_ion_adiabspecies = 1;

for (IndexSp const isp : dom_kinsp) {
// --> SpeciesInfo info
PC_tree_t const conf_isp = PCpp_get(conf_voicexx, ".SpeciesInfo[%d]", isp.uid());

kinetic_charges(isp) = static_cast<int>(PCpp_int(conf_isp, ".charge"));
if (kinetic_charges(isp) == -1) {
nb_elec_adiabspecies = 0;
} else {
nb_ion_adiabspecies = 0;
}

masses(isp) = PCpp_double(conf_isp, ".mass");
epsilon_bot(isp) = PCpp_double(conf_isp, ".epsilon_bot");
temperature_bot(isp) = PCpp_double(conf_isp, ".temperature_bot");
mean_velocity_bot(isp) = PCpp_double(conf_isp, ".mean_velocity_bot");
init_perturb_amplitude(isp) = PCpp_double(conf_isp, ".perturb_amplitude");
init_perturb_mode(isp) = static_cast<int>(PCpp_int(conf_isp, ".perturb_mode"));
}

// Create the domain of all species including kinetic species + adiabatic species (if existing)
IDomainSp const
dom_allsp(IndexSp(0), nb_kinspecies + nb_elec_adiabspecies + nb_ion_adiabspecies);
host_t<FieldSp<int>> charges(dom_allsp);
for (IndexSp isp : dom_kinsp) {
charges(isp) = kinetic_charges(isp);
}
if (nb_elec_adiabspecies + nb_ion_adiabspecies > 0) {
charges(dom_kinsp.back() + 1) = nb_ion_adiabspecies - nb_elec_adiabspecies;
}

// Initialization of the distribution function
ddc::init_discrete_space<IDimSp>(std::move(charges), std::move(masses));
DFieldSpVx allfequilibrium(meshSpVx);
BumpontailEquilibrium const init_fequilibrium(
std::move(epsilon_bot),
Expand Down Expand Up @@ -199,7 +174,7 @@ int main(int argc, char** argv)
expose_mesh_to_pdi("MeshX", mesh_x);
expose_mesh_to_pdi("MeshVx", mesh_vx);
ddc::expose_to_pdi("nbstep_diag", nbstep_diag);
ddc::expose_to_pdi("Nkinspecies", nb_kinspecies.value());
ddc::expose_to_pdi("Nkinspecies", dom_kinsp.size());
ddc::expose_to_pdi("fdistribu_charges", ddc::discrete_space<IDimSp>().charges()[dom_kinsp]);
ddc::expose_to_pdi("fdistribu_masses", ddc::discrete_space<IDimSp>().masses()[dom_kinsp]);
ddc::PdiEvent("initial_state").with("fdistribu_eq", allfequilibrium_host);
Expand Down
31 changes: 3 additions & 28 deletions simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
#include "restartinitialization.hpp"
#include "singlemodeperturbinitialization.hpp"
#include "species_info.hpp"
#include "species_init.hpp"
#include "spline_interpolator.hpp"
#include "splitvlasovsolver.hpp"

Expand Down Expand Up @@ -75,8 +76,7 @@ int main(int argc, char** argv)
IDomainVx mesh_vx(SplineInterpPointsVx::get_domain<IDimVx>());
IDomainXVx meshXVx(mesh_x, mesh_vx);

IVectSp const nb_kinspecies(PCpp_len(conf_voicexx, ".SpeciesInfo"));
IDomainSp const dom_kinsp(IndexSp(0), nb_kinspecies);
IDomainSp const dom_kinsp = init_species(conf_voicexx);

IDomainSpXVx const meshSpXVx(dom_kinsp, mesh_x, mesh_vx);
IDomainSpVx const meshSpVx(dom_kinsp, mesh_vx);
Expand All @@ -85,48 +85,23 @@ int main(int argc, char** argv)
SplineVxBuilder const builder_vx(meshXVx);
SplineVxBuilder_1d const builder_vx_poisson(mesh_vx);

host_t<FieldSp<int>> kinetic_charges(dom_kinsp);
host_t<DFieldSp> masses(dom_kinsp);
host_t<DFieldSp> epsilon_bot(dom_kinsp);
host_t<DFieldSp> temperature_bot(dom_kinsp);
host_t<DFieldSp> mean_velocity_bot(dom_kinsp);
host_t<DFieldSp> init_perturb_amplitude(dom_kinsp);
host_t<FieldSp<int>> init_perturb_mode(dom_kinsp);
int nb_elec_adiabspecies = 1;
int nb_ion_adiabspecies = 1;

for (IndexSp const isp : dom_kinsp) {
// --> SpeciesInfo info
PC_tree_t const conf_isp = PCpp_get(conf_voicexx, ".SpeciesInfo[%d]", isp.uid());

kinetic_charges(isp) = static_cast<int>(PCpp_int(conf_isp, ".charge"));
if (kinetic_charges(isp) == -1) {
nb_elec_adiabspecies = 0;
} else {
nb_ion_adiabspecies = 0;
}

masses(isp) = PCpp_double(conf_isp, ".mass");
epsilon_bot(isp) = PCpp_double(conf_isp, ".epsilon_bot");
temperature_bot(isp) = PCpp_double(conf_isp, ".temperature_bot");
mean_velocity_bot(isp) = PCpp_double(conf_isp, ".mean_velocity_bot");
init_perturb_amplitude(isp) = PCpp_double(conf_isp, ".perturb_amplitude");
init_perturb_mode(isp) = static_cast<int>(PCpp_int(conf_isp, ".perturb_mode"));
}

// Create the domain of all species including kinetic species + adiabatic species (if existing)
IDomainSp const
dom_allsp(IndexSp(0), nb_kinspecies + nb_elec_adiabspecies + nb_ion_adiabspecies);
host_t<FieldSp<int>> charges(dom_allsp);
for (IndexSp isp : dom_kinsp) {
charges(isp) = kinetic_charges(isp);
}
if (nb_elec_adiabspecies + nb_ion_adiabspecies > 0) {
charges(dom_kinsp.back() + 1) = nb_ion_adiabspecies - nb_elec_adiabspecies;
}

// Initialization of the distribution function
ddc::init_discrete_space<IDimSp>(std::move(charges), std::move(masses));
DFieldSpVx allfequilibrium(meshSpVx);
BumpontailEquilibrium const init_fequilibrium(
std::move(epsilon_bot),
Expand Down Expand Up @@ -198,7 +173,7 @@ int main(int argc, char** argv)
expose_mesh_to_pdi("MeshX", mesh_x);
expose_mesh_to_pdi("MeshVx", mesh_vx);
ddc::expose_to_pdi("nbstep_diag", nbstep_diag);
ddc::expose_to_pdi("Nkinspecies", nb_kinspecies.value());
ddc::expose_to_pdi("Nkinspecies", dom_kinsp.size());
ddc::expose_to_pdi("fdistribu_charges", ddc::discrete_space<IDimSp>().charges()[dom_kinsp]);
ddc::expose_to_pdi("fdistribu_masses", ddc::discrete_space<IDimSp>().masses()[dom_kinsp]);
ddc::PdiEvent("initial_state").with("fdistribu_eq", allfequilibrium_host);
Expand Down
30 changes: 3 additions & 27 deletions simulations/geometryXVx/landau/landau_fem_uniform.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@
#include "restartinitialization.hpp"
#include "singlemodeperturbinitialization.hpp"
#include "species_info.hpp"
#include "species_init.hpp"
#include "spline_interpolator.hpp"
#include "splitvlasovsolver.hpp"

Expand Down Expand Up @@ -69,8 +70,7 @@ int main(int argc, char** argv)
SplineInterpPointsVx>(conf_voicexx, "vx");
IDomainXVx const meshXVx(mesh_x, mesh_vx);

IVectSp const nb_kinspecies(PCpp_len(conf_voicexx, ".SpeciesInfo"));
IDomainSp const dom_kinsp(IndexSp(0), nb_kinspecies);
IDomainSp const dom_kinsp = init_species(conf_voicexx);

IDomainSpXVx const meshSpXVx(dom_kinsp, meshXVx);
IDomainSpVx const meshSpVx(dom_kinsp, mesh_vx);
Expand All @@ -80,48 +80,24 @@ int main(int argc, char** argv)
SplineVxBuilder const builder_vx(meshXVx);
SplineVxBuilder_1d const builder_vx_poisson(mesh_vx);

host_t<FieldSp<int>> kinetic_charges(dom_kinsp);
host_t<DFieldSp> masses(dom_kinsp);
host_t<DFieldSp> density_eq(dom_kinsp);
host_t<DFieldSp> temperature_eq(dom_kinsp);
host_t<DFieldSp> mean_velocity_eq(dom_kinsp);
host_t<DFieldSp> init_perturb_amplitude(dom_kinsp);
host_t<FieldSp<int>> init_perturb_mode(dom_kinsp);
int nb_elec_adiabspecies = 1;
int nb_ion_adiabspecies = 1;

for (IndexSp const isp : dom_kinsp) {
// --> SpeciesInfo info
PC_tree_t const conf_isp = PCpp_get(conf_voicexx, ".SpeciesInfo[%d]", isp.uid());

kinetic_charges(isp) = static_cast<int>(PCpp_int(conf_isp, ".charge"));
if (kinetic_charges(isp) == -1) {
nb_elec_adiabspecies = 0;
} else {
nb_ion_adiabspecies = 0;
}

masses(isp) = PCpp_double(conf_isp, ".mass");
density_eq(isp) = PCpp_double(conf_isp, ".density_eq");
temperature_eq(isp) = PCpp_double(conf_isp, ".temperature_eq");
mean_velocity_eq(isp) = PCpp_double(conf_isp, ".mean_velocity_eq");
init_perturb_amplitude(isp) = PCpp_double(conf_isp, ".perturb_amplitude");
init_perturb_mode(isp) = static_cast<int>(PCpp_int(conf_isp, ".perturb_mode"));
}

// Create the domain of all species including kinetic species + adiabatic species (if existing)
IDomainSp const
dom_allsp(IndexSp(0), nb_kinspecies + nb_elec_adiabspecies + nb_ion_adiabspecies);
host_t<FieldSp<int>> charges(dom_allsp);
for (IndexSp isp : dom_kinsp) {
charges(isp) = kinetic_charges(isp);
}
if (nb_elec_adiabspecies + nb_ion_adiabspecies > 0) {
charges(dom_kinsp.back() + 1) = nb_ion_adiabspecies - nb_elec_adiabspecies;
}

// Initialization of the distribution function
ddc::init_discrete_space<IDimSp>(std::move(charges), std::move(masses));
DFieldSpVx allfequilibrium(meshSpVx);
MaxwellianEquilibrium const init_fequilibrium(
std::move(density_eq),
Expand Down Expand Up @@ -198,7 +174,7 @@ int main(int argc, char** argv)
expose_mesh_to_pdi("MeshX", mesh_x);
expose_mesh_to_pdi("MeshVx", mesh_vx);
ddc::expose_to_pdi("nbstep_diag", nbstep_diag);
ddc::expose_to_pdi("Nkinspecies", nb_kinspecies.value());
ddc::expose_to_pdi("Nkinspecies", dom_kinsp.size());
ddc::expose_to_pdi("fdistribu_charges", ddc::discrete_space<IDimSp>().charges()[dom_kinsp]);
ddc::expose_to_pdi("fdistribu_masses", ddc::discrete_space<IDimSp>().masses()[dom_kinsp]);
ddc::PdiEvent("initial_state").with("fdistribu_eq", allfequilibrium_host);
Expand Down
32 changes: 3 additions & 29 deletions simulations/geometryXVx/landau/landau_fft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
#include "restartinitialization.hpp"
#include "singlemodeperturbinitialization.hpp"
#include "species_info.hpp"
#include "species_init.hpp"
#include "spline_interpolator.hpp"
#include "splitvlasovsolver.hpp"

Expand Down Expand Up @@ -67,58 +68,31 @@ int main(int argc, char** argv)
SplineInterpPointsVx>(conf_voicexx, "vx");
IDomainXVx const meshXVx(mesh_x, mesh_vx);

IVectSp const nb_kinspecies(PCpp_len(conf_voicexx, ".SpeciesInfo"));
IDomainSp const dom_kinsp(IndexSp(0), nb_kinspecies);

IDomainSp const dom_kinsp = init_species(conf_voicexx);
IDomainSpXVx const meshSpXVx(dom_kinsp, meshXVx);
IDomainSpVx const meshSpVx(dom_kinsp, mesh_vx);

SplineXBuilder const builder_x(meshXVx);
SplineVxBuilder const builder_vx(meshXVx);
SplineVxBuilder_1d const builder_vx_poisson(mesh_vx);

host_t<FieldSp<int>> kinetic_charges(dom_kinsp);
host_t<DFieldSp> masses(dom_kinsp);
host_t<DFieldSp> density_eq(dom_kinsp);
host_t<DFieldSp> temperature_eq(dom_kinsp);
host_t<DFieldSp> mean_velocity_eq(dom_kinsp);
host_t<DFieldSp> init_perturb_amplitude(dom_kinsp);
host_t<FieldSp<int>> init_perturb_mode(dom_kinsp);
int nb_elec_adiabspecies = 1;
int nb_ion_adiabspecies = 1;

for (IndexSp const isp : dom_kinsp) {
// --> SpeciesInfo info
PC_tree_t const conf_isp = PCpp_get(conf_voicexx, ".SpeciesInfo[%d]", isp.uid());

kinetic_charges(isp) = static_cast<int>(PCpp_int(conf_isp, ".charge"));
if (kinetic_charges(isp) == -1) {
nb_elec_adiabspecies = 0;
} else {
nb_ion_adiabspecies = 0;
}

masses(isp) = PCpp_double(conf_isp, ".mass");
density_eq(isp) = PCpp_double(conf_isp, ".density_eq");
temperature_eq(isp) = PCpp_double(conf_isp, ".temperature_eq");
mean_velocity_eq(isp) = PCpp_double(conf_isp, ".mean_velocity_eq");
init_perturb_amplitude(isp) = PCpp_double(conf_isp, ".perturb_amplitude");
init_perturb_mode(isp) = static_cast<int>(PCpp_int(conf_isp, ".perturb_mode"));
}

// Create the domain of all species including kinetic species + adiabatic species (if existing)
IDomainSp const
dom_allsp(IndexSp(0), nb_kinspecies + nb_elec_adiabspecies + nb_ion_adiabspecies);
host_t<FieldSp<int>> charges(dom_allsp);
for (IndexSp isp : dom_kinsp) {
charges(isp) = kinetic_charges(isp);
}
if (nb_elec_adiabspecies + nb_ion_adiabspecies > 0) {
charges(dom_kinsp.back() + 1) = nb_ion_adiabspecies - nb_elec_adiabspecies;
}

// Initialization of the distribution function
ddc::init_discrete_space<IDimSp>(std::move(charges), std::move(masses));
DFieldSpVx allfequilibrium(meshSpVx);
MaxwellianEquilibrium const init_fequilibrium(
std::move(density_eq),
Expand Down Expand Up @@ -191,7 +165,7 @@ int main(int argc, char** argv)
expose_mesh_to_pdi("MeshX", mesh_x);
expose_mesh_to_pdi("MeshVx", mesh_vx);
ddc::expose_to_pdi("nbstep_diag", nbstep_diag);
ddc::expose_to_pdi("Nkinspecies", nb_kinspecies.value());
ddc::expose_to_pdi("Nkinspecies", dom_kinsp.size());
ddc::expose_to_pdi("fdistribu_charges", ddc::discrete_space<IDimSp>().charges()[dom_kinsp]);
ddc::expose_to_pdi("fdistribu_masses", ddc::discrete_space<IDimSp>().masses()[dom_kinsp]);
ddc::PdiEvent("initial_state").with("fdistribu_eq", allfequilibrium_host);
Expand Down
Loading

0 comments on commit 8825ed4

Please sign in to comment.