diff --git a/include/vantage.hxx b/include/vantage.hxx index e4586649e..57963ab60 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -9,31 +9,6 @@ using namespace NESO::Particles; using namespace VANTAGE::Reactions; -struct Vantage : public Component { - Vantage(std::string name, Options& options, Solver* solver); - - ~Vantage(); // Destructor for VANTAGE related cleanup - void finally(const Options& state) override; - void transform_impl(GuardedOptions& state) override; - void outputVars(Options& state) override; - -private: - PetscLib petsc_lib; // Ensures PETSc is initialized for the lifetime of this component - std::string name; // Component name - DM dm; - std::shared_ptr neso_mesh; - std::shared_ptr sycl_target; - std::shared_ptr b2d; - - Field2D ion_density; - Field2D neutral_density; - BoutReal particle_time; -}; - -namespace { -RegisterComponent registercomponentvantage("vantage"); -} - /// @brief Data struct to hold information about a reaction source. /// @param reaction_name Name of the reaction, e.g. "ionistaion" /// @param source_name Name of the source, e.g. Siz (ion density source due to @@ -86,6 +61,99 @@ private: Options& units; }; +// Need to declare empty struct because the Monitor needs it and it must be +// before component construction as it has the monitor as a member. +// See https://bout-dev.readthedocs.io/en/latest/user_docs/time_integration.html#monitoring-the-simulation-output +// for more information on monitors. +struct Vantage; + +/// @brief Monitor to schedule kinetic iterations. +/// @param solver pointer to the solver. +/// @param time normalised sim time provided by BOUT++ +/// @param iter current iteration number provided by BOUT++ +/// @param nout number of outputs provided by BOUT++ +class VantageMonitor : public Monitor { +public: + explicit VantageMonitor(Vantage* vantage) : vantage(vantage) {} + int call(Solver* solver, BoutReal time, int iter, int nout) override; + +private: + Vantage* vantage; +}; + +struct Vantage : public Component { + Vantage(std::string name, Options& options, Solver* solver); + + ~Vantage(); // Destructor for VANTAGE related cleanup + void finally(const Options& state) override; + void transform_impl(GuardedOptions& state) override; + void outputVars(Options& state) override; + + // Function with the kinetic loop. + // time is the normalised VANTAGE monitor frequency. + int advance_vantage(BoutReal time); + +private: + bool test_mass_conservation; + BoutReal N_w; + REAL dt; + int nsteps; + int num_cells_owned; // Number of VANTAGE cells owned per rank + + Options bout_output_data; // Options object to hold output data for VANTAGE diagnostics + std::unique_ptr + vantage_dump_writer; // OptionsIO object to write VANTAGE diagnostics + + int mpi_rank; // Current rank ID + Mesh* bout_mesh; // Pointer to the BOUT++ mesh object + Field2D ion_density, neutral_density, total_density; + Field2D initial_neutral_density; // Initial VANTAGE kinetic neutral density + BoutReal total_mass_initial, total_mass; + std::string dmplex_filepath, vantage_dump_filepath, + particle_data_filepath; // Path for output files + + // Physics + BoutReal particle_time; + std::string neutral_species; // Neutral species to simulate with VANTAGE + std::string ion_species; + bool plasma_coupling; // Whether to read plasma fields from the state or not + PetscLib petsc_lib; // Ensures PETSc is initialized for the lifetime of this component + BoutReal background_ion_temperature, background_ion_density; + std::vector V_background; + + // DMPLex/VANTAGE stuff + DM dm; + std::shared_ptr neso_mesh; + std::shared_ptr sycl_target; + std::shared_ptr + b2d; // Boundary interaction object + std::shared_ptr + dg0; // DMPlex projection object + std::shared_ptr h5part; // HDF5 particle output object + std::vector + h_project1; // Buffer for scalar projection/evaluation of NESO-Particles properties + std::shared_ptr A_particle_group; // Particle group for main neutrals + std::shared_ptr marker_group; // Particle group for rec markers + std::shared_ptr> + plasma_data; // Data structure for plasma in VANTAGE grid + + // Needed for Vantage::apply_boundary_conditions + std::shared_ptr reflection; // Boundary reflection object + void apply_boundary_conditions(ParticleSubGroupSharedPtr aa); + + // These classes don't have a default constructor so need to be initialised as a unique_ptr + std::unique_ptr + source_manager; // Manager for VANTAGE reaction sources + VantageMonitor monitor{this}; // Output monitor to schedule VANTAGE iterations + + std::unique_ptr reaction_controller; + std::unique_ptr recombination_controller; +}; + +namespace { +RegisterComponent registercomponentvantage("vantage"); +} + /** * @brief Function to calculate cell volumes. * diff --git a/src/vantage.cxx b/src/vantage.cxx index def6d6ddc..8787d45df 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -22,9 +22,9 @@ #include #include // for reactions integration +#include "../include/amjuel_data.hxx" #include "../include/vantage.hxx" #include "../include/vantage_dmplex.hxx" -#include "../include/amjuel_data.hxx" #include #ifndef NESO_PARTICLES_PETSC @@ -42,8 +42,8 @@ inline void ASSERT_EQ(T t, U u) { // Helper function to convert AMJUEL rate from Hermes-3 to Reactions format // Hermes-3: vector of vectors of BoutReal // Reactions: array of REAL -std::array, 9> convert_amjuel_format( - const std::vector> & coeffs) { +std::array, 9> +convert_amjuel_format(const std::vector>& coeffs) { std::array, 9> out{}; for (std::size_t i = 0; i < 9; ++i) { @@ -85,10 +85,11 @@ void calculate_neutral_density_in_place( // extrapolate -> Neumann } -double calculate_total_mass(Field2D& density, - std::shared_ptr& neso_mesh) { - double local_mass = 0.0; - double total_mass = 0.0; +BoutReal +calculate_total_mass(Field2D& density, + std::shared_ptr& neso_mesh) { + BoutReal local_mass = 0.0; + BoutReal total_mass = 0.0; Mesh* bout_mesh = density.getMesh(); PetscInt ic = 0; for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { @@ -103,16 +104,15 @@ double calculate_total_mass(Field2D& density, } Options -initialise_diagnostics(Options& alloptions, - Mesh* bout_mesh, - Field2D& neutral_density, Field2D& ion_density, - std::shared_ptr& neso_mesh, - std::string particle_data_filename) { +initialise_diagnostics(Options& alloptions, Mesh* bout_mesh, Field2D& neutral_density, + Field2D& ion_density, + std::shared_ptr& neso_mesh, + std::string vantage_dump_filepath) { // Options object to use to write out diagnostic data of fluid quantities auto Nnorm = get(alloptions["units"]["inv_meters_cubed"]); auto Tnorm = get(alloptions["units"]["eV"]); - auto Omega_ci = 1/get(alloptions["units"]["seconds"]); + auto Omega_ci = 1 / get(alloptions["units"]["seconds"]); auto rho_s0 = get(alloptions["units"]["meters"]); auto Bnorm = get(alloptions["units"]["Tesla"]); auto Cs0 = get(alloptions["units"]["meters"]) @@ -122,8 +122,7 @@ initialise_diagnostics(Options& alloptions, set_with_attrs(bout_output_data["neutral_density"], neutral_density, {{"time_dimension", "t"}}); - set_with_attrs(bout_output_data["ion_density"], ion_density, - {{"time_dimension", "t"}}); + set_with_attrs(bout_output_data["ion_density"], ion_density, {{"time_dimension", "t"}}); set_with_attrs(bout_output_data["Nn"], neutral_density, {{"time_dimension", "t"}, @@ -163,8 +162,7 @@ initialise_diagnostics(Options& alloptions, {{"time_dimension", "t"}}); set_with_attrs(bout_output_data["total_ion_mass"], - calculate_total_mass(ion_density, neso_mesh), - {{"time_dimension", "t"}}); + calculate_total_mass(ion_density, neso_mesh), {{"time_dimension", "t"}}); set_with_attrs(bout_output_data["t_array"], 0.0, {{"time_dimension", "t"}}); @@ -172,51 +170,45 @@ initialise_diagnostics(Options& alloptions, bout_mesh->outputVars(bout_output_data); // Add metadata with normalisation factors - set_with_attrs(bout_output_data["Tnorm"], Tnorm, { - {"units", "eV"}, - {"conversion", 1}, // Already in SI units - {"standard_name", "temperature normalisation"}, - {"long_name", "temperature normalisation"} - }); - set_with_attrs(bout_output_data["Nnorm"], Nnorm, { - {"units", "m^-3"}, - {"conversion", 1}, - {"standard_name", "density normalisation"}, - {"long_name", "Number density normalisation"} - }); - set_with_attrs(bout_output_data["Bnorm"], Bnorm, { - {"units", "T"}, - {"conversion", 1}, - {"standard_name", "magnetic field normalisation"}, - {"long_name", "Magnetic field normalisation"} - }); - set_with_attrs(bout_output_data["Cs0"], Cs0, { - {"units", "m/s"}, - {"conversion", 1}, - {"standard_name", "velocity normalisation"}, - {"long_name", "Sound speed normalisation"} - }); - set_with_attrs(bout_output_data["Omega_ci"], Omega_ci, { - {"units", "s^-1"}, - {"conversion", 1}, - {"standard_name", "frequency normalisation"}, - {"long_name", "Cyclotron frequency normalisation"} - }); - set_with_attrs(bout_output_data["rho_s0"], rho_s0, { - {"units", "m"}, - {"conversion", 1}, - {"standard_name", "length normalisation"}, - {"long_name", "Gyro-radius length normalisation"} - }); - - bout::OptionsIO::create(particle_data_filename)->write(bout_output_data); + set_with_attrs(bout_output_data["Tnorm"], Tnorm, + {{"units", "eV"}, + {"conversion", 1}, // Already in SI units + {"standard_name", "temperature normalisation"}, + {"long_name", "temperature normalisation"}}); + set_with_attrs(bout_output_data["Nnorm"], Nnorm, + {{"units", "m^-3"}, + {"conversion", 1}, + {"standard_name", "density normalisation"}, + {"long_name", "Number density normalisation"}}); + set_with_attrs(bout_output_data["Bnorm"], Bnorm, + {{"units", "T"}, + {"conversion", 1}, + {"standard_name", "magnetic field normalisation"}, + {"long_name", "Magnetic field normalisation"}}); + set_with_attrs(bout_output_data["Cs0"], Cs0, + {{"units", "m/s"}, + {"conversion", 1}, + {"standard_name", "velocity normalisation"}, + {"long_name", "Sound speed normalisation"}}); + set_with_attrs(bout_output_data["Omega_ci"], Omega_ci, + {{"units", "s^-1"}, + {"conversion", 1}, + {"standard_name", "frequency normalisation"}, + {"long_name", "Cyclotron frequency normalisation"}}); + set_with_attrs(bout_output_data["rho_s0"], rho_s0, + {{"units", "m"}, + {"conversion", 1}, + {"standard_name", "length normalisation"}, + {"long_name", "Gyro-radius length normalisation"}}); + + bout::OptionsIO::create(vantage_dump_filepath)->write(bout_output_data); return bout_output_data; } -void update_diagnostics(Field2D& neutral_density, Field2D& ion_density, - Field2D& Siz, Field2D& Srec, +void update_diagnostics(Field2D& neutral_density, Field2D& ion_density, Field2D& Siz, + Field2D& Srec, std::shared_ptr& neso_mesh, - Options& bout_output_data, std::string particle_data_filename, + Options& bout_output_data, bout::OptionsIO& vantage_dump_writer, BoutReal particle_time) { // update density in Options object and write bout_output_data["neutral_density"] = neutral_density; @@ -232,8 +224,9 @@ void update_diagnostics(Field2D& neutral_density, Field2D& ion_density, bout_output_data["t_array"] = particle_time; // bout_output_data["t_array"] = 0.0; // Append data to file - bout::OptionsIO::create({{"file", particle_data_filename}, {"append", true}}) - ->write(bout_output_data); + vantage_dump_writer.write(bout_output_data); + vantage_dump_writer + .flush(); // Ensure buffer is written to disk to avoid crash data loss } void set_initial_particle_weights( @@ -241,8 +234,7 @@ void set_initial_particle_weights( std::shared_ptr& dg0, std::shared_ptr& A_particle_group, std::shared_ptr& neso_mesh, - std::vector& h_project1, - BoutReal N_w) { + std::vector& h_project1, BoutReal N_w) { Mesh* bout_mesh = initial_neutral_density.getMesh(); PetscInt ixy = 0; for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { @@ -250,13 +242,12 @@ void set_initial_particle_weights( // particle_weights are copied to all particles in this cell. // we multiply the initial density by the volume to get particle number, // then divide by markers per cell to divide them between the requested markers, - // then divide by N_w to get the weight of each marker. + // then divide by N_w to get the weight of each marker. const REAL cell_volume = neso_mesh->dmh->get_cell_volume(static_cast(ixy)); const INT nmarkers_per_cell = A_particle_group->get_npart_cell(static_cast(ixy)); const REAL particle_weights = initial_neutral_density(ix, iy) * cell_volume - / static_cast(nmarkers_per_cell) - / N_w; + / static_cast(nmarkers_per_cell) / N_w; h_project1.at(static_cast(ixy)) = particle_weights; ixy++; } @@ -266,6 +257,122 @@ void set_initial_particle_weights( // set the data from internal variables into the weights dg0->evaluate(A_particle_group, Sym("WEIGHT")); } +// Enum storing row IDs for the plasma data +enum PlasmaRow { row_id_ne = 0, row_id_te, row_id_ni, row_id_ti, plasma_nrows }; + +// Project plasma data onto VANTAGE grid +// Copy the plasma data into the VANTAGE grid, done once per VANTAGE call. +// Sends the data to SYCL target. +void send_plasma_data(std::shared_ptr>& plasma_data, + std::shared_ptr& sycl_target, Mesh* bout_mesh, + int num_cells_owned, BoutReal background_ion_density, + BoutReal background_ion_temperature) { + + // Assemble vector of plasma data + std::vector> cell_data; + cell_data.reserve(static_cast(num_cells_owned)); + + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + auto cell = std::make_shared>(sycl_target, plasma_nrows, 1); + cell->at(row_id_ne, 0) = background_ion_density; + cell->at(row_id_te, 0) = background_ion_temperature; + cell->at(row_id_ni, 0) = background_ion_density; + cell->at(row_id_ti, 0) = background_ion_temperature; + cell_data.push_back(cell); + } + } + + // Upload cell data to device + plasma_data->set_all_cells(cell_data); +} + +// Calculate recombination marker properties +// Set distribution based on plasma properties +// Also send plasma properties themselves +void update_recombination_markers( + std::shared_ptr& marker_group, + std::shared_ptr>& plasma_data, + std::shared_ptr& neso_mesh, BoutReal N_w, + const std::vector& V_background) { + + // number of velocity dimensions + const int nvel = static_cast(V_background.size()); + + // Give particle group initial fluid values: markers will contain background + // plasma properties + // From demo app "set_init_fluid_values" + particle_loop( + "set init fluid values", marker_group, + [=](auto plasma_data, auto n, auto T, auto ne, auto Te, auto speed) { + n.at(0) = plasma_data.at(row_id_ni, 0); + ne.at(0) = plasma_data.at(row_id_ne, 0); + T.at(0) = plasma_data.at(row_id_ti, 0); + Te.at(0) = plasma_data.at(row_id_te, 0); + + // Keep background velocity as scalar for now + for (int i = 0; i < nvel; i++) { + speed.at(i) = V_background[static_cast(i)]; + } + }, + Access::read(plasma_data), Access::write(Sym("FLUID_DENSITY")), + Access::write(Sym("FLUID_TEMPERATURE")), + Access::write(Sym("ELECTRON_DENSITY")), + Access::write(Sym("ELECTRON_TEMPERATURE")), + Access::write(Sym("FLUID_FLOW_SPEED"))) + ->execute(); + + // Calculate marker weights + + // Add particle property: number of particles in the local cell + // From demo app: "distribute_n_part_cell" + for (int ic = 0; ic < marker_group->domain->mesh->get_cell_count(); ic++) { + INT n_part_cell = marker_group->get_npart_cell(ic); + particle_loop( + "Update N_CELL prop", marker_group, + [=](auto n_cell_prop) { n_cell_prop.at(0) = n_part_cell; }, + Access::write(Sym("N_CELL"))) + ->execute(ic); + } + + const INT num_cells = marker_group->domain->mesh->get_cell_count(); + const BoutReal N_w_local = N_w; + + // Calculate weight for each marker particle + // based on FLUID_DENSITY and N_CELL properties contained in same particle. + for (int ic = 0; ic < num_cells; ic++) { + REAL V_cell = neso_mesh->dmh->get_cell_volume(ic); + particle_loop( + "Update weight of ions", marker_group, + [N_w_local, V_cell](auto n_cell_prop, auto ion_dens_prop, auto weight_prop) { + const BoutReal n_cell = static_cast(n_cell_prop.at(0)); + auto updated_weight = (ion_dens_prop.at(0) * V_cell) / (N_w_local * n_cell); + weight_prop.at(0) = updated_weight; + }, + Access::read(Sym("N_CELL")), Access::read(Sym("FLUID_DENSITY")), + Access::write(Sym("WEIGHT"))) + ->execute(ic); + } +} + +// Update plasma data projection +// Plasma data stays on device and particles read it every timestep. +// This is necessary as particles move between cells. +void update_particle_background(std::shared_ptr& A_particle_group, + std::shared_ptr>& plasma_data) { + + particle_loop( + "update neutral plasma properties", A_particle_group, + [=](auto plasma_data, auto ni, auto ne, auto te) { + ni.at(0) = plasma_data.at(row_id_ni, 0); + ne.at(0) = plasma_data.at(row_id_ne, 0); + te.at(0) = plasma_data.at(row_id_te, 0); + }, + Access::read(plasma_data), Access::write(Sym("ION_DENSITY")), + Access::write(Sym("ELECTRON_DENSITY")), + Access::write(Sym("ELECTRON_TEMPERATURE"))) + ->execute(); +} void check_cell_volumes(std::shared_ptr& neso_mesh, Mesh*& bout_mesh, Options& alloptions) { @@ -316,7 +423,8 @@ REAL cell_length(std::vector>& cell_vertices, std::size_t iv1, return length; } -void check_cell_centres(Options& alloptions, std::shared_ptr& neso_mesh, +void check_cell_centres(Options& alloptions, + std::shared_ptr& neso_mesh, Mesh*& bout_mesh, BoutReal absolute_tolerance, BoutReal relative_tolerance) { // get (R,Z) of cell centres in Hypnotoad grid @@ -375,9 +483,9 @@ void check_cell_centres(Options& alloptions, std::shared_ptr(0, BoutComm::get()); + + if (bout_mesh->GlobalNz != 1) { + throw BoutException( + "VANTAGE currently only supports 2D grids, but the provided grid has " + "GlobalNz = {}. Set MZ = 1 on the top of the input file (root level).", + bout_mesh->GlobalNz); + } + // TODO: Put proper permissions in // int main(int argc, char** argv) { @@ -472,9 +591,8 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // Mesh* bout_mesh = Mesh::create(&Options::root()["mesh"]); // TODO: tidy up the above - Options& mesh_options = alloptions["dmplex"]; // [mesh] - Options& options = alloptions[name]; // [vantage] + Options& options = alloptions[name]; // [vantage] Options& units = alloptions["units"]; BoutReal inv_meters_cubed = get(units["inv_meters_cubed"]); @@ -482,27 +600,33 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) BoutReal meters = get(units["meters"]); BoutReal seconds = get(units["seconds"]); - BoutReal N_w = options["N_w"] - .doc("Normalisation parameter: number of particles (normalised) per " - "unit weight. Default = 1.1 as a value close but different to unity" - "to make sure an incorrect implementation would show up in tests.") - .withDefault(1.1); + N_w = options["N_w"] + .doc("Normalisation parameter: number of particles (normalised) per " + "unit weight. Default = 1.1 as a value close but different to unity" + "to make sure an incorrect implementation would show up in tests.") + .withDefault(1.1); Options::root()["units"]["N_w"] = N_w; Options::root()["units"]["N_w"].setConditionallyUsed(); - Mesh* bout_mesh = bout::globals::mesh; - sycl_target = std::make_shared(0, BoutComm::get()); // keep dmplex_h5_filename in vantage.cxx to retain access to make_output_path() // which should presumably not need to exist within the hermes-3 library std::string dmplex_h5_filename = mesh_options["dmplex_h5_filename"] .doc("Filename to use for saving the DMPlex mesh") .withDefault("hypnotoad_dmplex_mesh_output.h5"); + // Create paths + dmplex_filepath = make_output_path(dmplex_h5_filename, alloptions); + mpi_rank = sycl_target->comm_pair.rank_parent; + vantage_dump_filepath = + make_output_path(fmt::format("BOUT.dmp.vantage.{}.nc", mpi_rank), alloptions); + particle_data_filepath = make_output_path("particle_trajectories.h5part", alloptions); + // Create and save DMPlex // This is in SI units. - dm = create_dmplex_from_Bout_mesh(bout_mesh, mesh_options, sycl_target, - make_output_path(dmplex_h5_filename, alloptions)); + + dm = + create_dmplex_from_Bout_mesh(bout_mesh, mesh_options, sycl_target, dmplex_filepath); // Normalise DMPlex after creation // Get local coords object (i.e. per rank) and scale it - this scales entire mesh @@ -510,11 +634,9 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) Vec coords = nullptr; PETSCCHK(DMGetCoordinatesLocal(dm, &coords)); if (coords != nullptr) { - PETSCCHK(VecScale(coords, 1/meters)); + PETSCCHK(VecScale(coords, 1 / meters)); } - output << "Begin particle push \n"; - /* * * @@ -531,27 +653,30 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) */ { + test_mass_conservation = + options["test_mass_conservation"] + .doc("Check that mass is conserved at runtime. Default = true") + .withDefault(true); + // Normalisations - // Initial neutral parameters - Field2D initial_neutral_density{bout_mesh}; + // Initial neutral parameters initial_neutral_density = options["initial_neutral_density"] - .doc( - "Initial neutral density for VANTAGE kinetic neutrals [m^-3]") + .doc("Initial neutral density for VANTAGE kinetic neutrals [m^-3]") .as() - / inv_meters_cubed; + / inv_meters_cubed; const int npart_per_cell = options["npart_per_cell"] .doc("Number of VANTAGE kinetic neutral particles per " "cell during initialisation") .withDefault(1); // Plasma parameters - const BoutReal background_ion_temperature = + background_ion_temperature = options["background_ion_temperature"] .doc("Background ion temp override [eV], default = 10") .withDefault(10) / eV; - const BoutReal background_ion_density = + background_ion_density = options["background_ion_density"] .doc("Background density override [m^-3], default = 1.0e19") .withDefault(1.0e19) @@ -562,7 +687,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) const BoutReal background_ion_Vx = options["background_ion_Vx"].withDefault(0.0); const BoutReal background_ion_Vy = options["background_ion_Vy"].withDefault(0.0); - const std::vector V_background = {background_ion_Vx, background_ion_Vy}; + V_background = {background_ion_Vx, background_ion_Vy}; // Reaction settings const REAL iz_rate_override = @@ -574,22 +699,39 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) options["rec_rate_override"] .doc("Recombination rate override (weight s^-1, normalised units).") .withDefault(-1.0); - const int rec_markers_per_cell = - options["rec_markers_per_cell"].withDefault(1000); + const int rec_markers_per_cell = options["rec_markers_per_cell"].withDefault(1000); // Other settings + neutral_species = + options["neutral_species"].doc("Name of the neutral species to use in VANTAGE"); + + // Note we're using Options& to avoid copying the Options + Options& ion_species_opt = + options["ion_species"].doc("Name of the ion correspondng to the VANTAGE neutral " + "species. If unset, plasma coupling is disabled"); + + // Disable plasma coupling if ions unset. + // ion_species is the option content from ion_species_opt. + plasma_coupling = ion_species_opt.isSet(); + if (plasma_coupling) { + ion_species = ion_species_opt.as(); + output_info.write("\tVANTAGE: plasma coupling with ion species '{:s}' enabled!\n", + ion_species); + } else { + output_info.write("\tVANTAGE: no ion_species set, plasma coupling disabled!\n"); + } + const int ndim = 2; - const REAL dt = options["dt"] - .doc("Timestep to use for VANTAGE kinetic neutrals (normalised units)") - .withDefault(0.01); - const int nsteps = options["nsteps"] - .doc("Number of timesteps to use for VANTAGE kinetic neutrals") - .withDefault(10); + dt = options["dt"] + .doc("Timestep to use for VANTAGE kinetic neutrals (normalised units)") + .withDefault(0.01); + nsteps = options["nsteps"] + .doc("Number of timesteps to use for VANTAGE kinetic neutrals") + .withDefault(10); const int rng_samples = options["rng_samples"] .doc("Number of RNG samples to prepare per-particle") .withDefault(40); - BoutReal particle_time = 0.0; ion_density = Field2D(background_ion_density, bout_mesh); neutral_density = Field2D(0.0, bout_mesh); // Create a mesh interface from the DM @@ -600,7 +742,13 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // Create a domain from the neso_mesh and the mapper. auto domain = std::make_shared(neso_mesh, mapper); // Get the number of cells in the mesh owned on this process - int num_cells_owned = neso_mesh->get_cell_count(); + num_cells_owned = neso_mesh->get_cell_count(); + + // Data structure for plasma data on VANTAGE grid. One row per variable. + // Save ID of each row to a variable. + plasma_data = std::make_shared>(sycl_target, num_cells_owned, + plasma_nrows, 1); + // if requested, check that neso_mesh cell volumes are identical // to bout_mesh cell volumes, otherwise, exit. if (mesh_options["test_dmplex_cell_volumes"].withDefault(true)) { @@ -608,8 +756,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) } if (mesh_options["test_dmplex_cell_centres"].withDefault(true)) { check_cell_centres( - alloptions, - neso_mesh, bout_mesh, + alloptions, neso_mesh, bout_mesh, mesh_options["dmplex_cell_centre_absolute_tolerance"].withDefault(1.0e-12), mesh_options["dmplex_cell_centre_relative_tolerance"].withDefault(0.0)); } @@ -638,11 +785,11 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) ParticleSpec particle_spec = particle_spec_builder.get_particle_spec(); // Create a Particle group with our specied particle properties. - auto A_particle_group = + A_particle_group = std::make_shared(domain, particle_spec, sycl_target); // Create some particle data - const int mpi_rank = sycl_target->comm_pair.rank_parent; + std::mt19937 rng_pos(static_cast(52234234 + mpi_rank)); std::mt19937 rng_vel(static_cast(52234231 + mpi_rank)); std::vector> positions; @@ -674,8 +821,10 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) initial_distribution[Sym("ION_DENSITY")][px][0] = background_ion_density; initial_distribution[Sym("ION_SOURCE_DENSITY")][px][0] = 0.0; initial_distribution[Sym("ION_SOURCE_ENERGY")][px][0] = 0.0; - initial_distribution[Sym("ELECTRON_DENSITY")][px][0] = background_electron_density; - initial_distribution[Sym("ELECTRON_TEMPERATURE")][px][0] = background_electron_temperature; + initial_distribution[Sym("ELECTRON_DENSITY")][px][0] = + background_electron_density; + initial_distribution[Sym("ELECTRON_TEMPERATURE")][px][0] = + background_electron_temperature; initial_distribution[Sym("ELECTRON_SOURCE_DENSITY")][px][0] = 0.0; initial_distribution[Sym("ELECTRON_SOURCE_ENERGY")][px][0] = 0.0; for (int dimx = 0; dimx < ndim; dimx++) { @@ -689,8 +838,8 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // Add the new particles to the particle group A_particle_group->add_particles_local(initial_distribution); // make pointer to projection object - auto dg0 = std::make_shared( - neso_mesh, sycl_target, "DG", 0); + dg0 = std::make_shared(neso_mesh, + sycl_target, "DG", 0); // RNG kernel // Used for sampling from velocity distribution for REC/CX @@ -699,8 +848,6 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) auto rng_kernel = get_uniform_rng_kernel(sycl_target, static_cast(rng_samples)); - - // Recombination reaction // ------------------------------------------------------------------------------ // Create Maxwellian distribution of recombination markers according to fluid @@ -710,8 +857,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // Options and constants // Make new particle group just for the markers - auto marker_group = - std::make_shared(domain, particle_spec, sycl_target); + marker_group = std::make_shared(domain, particle_spec, sycl_target); // Give particle group initial kinetic values (positions and velocities) // Numerical settings: weight, stdev, species ID @@ -720,63 +866,11 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) marker_group->add_particles_local(maxwellian_markers); - // Give particle group initial fluid values: markers will contain background - // plasma properties - // From demo app "set_init_fluid_values" - particle_loop( - "set init fluid values", marker_group, - [=](auto n, auto T, auto ne, auto Te, auto speed) { - n.at(0) = background_ion_density; - ne.at(0) = background_ion_density; - T.at(0) = background_ion_temperature; - Te.at(0) = background_ion_temperature; - - for (int i = 0; i < ndim; i++) { - speed.at(i) = V_background[static_cast(i)]; - } - }, - Access::write(Sym("FLUID_DENSITY")), - Access::write(Sym("FLUID_TEMPERATURE")), - Access::write(Sym("ELECTRON_DENSITY")), - Access::write(Sym("ELECTRON_TEMPERATURE")), - Access::write(Sym("FLUID_FLOW_SPEED"))) - ->execute(); - - // Calculate marker weights - - // Add particle property: number of particles in the local cell - // From demo app: "distribute_n_part_cell" - for (int ic = 0; ic < marker_group->domain->mesh->get_cell_count(); ic++) { - INT n_part_cell = marker_group->get_npart_cell(ic); - particle_loop( - "Update N_CELL prop", marker_group, - [=](auto n_cell_prop) { n_cell_prop.at(0) = n_part_cell; }, - Access::write(Sym("N_CELL"))) - ->execute(ic); - } - - const INT num_cells = marker_group->domain->mesh->get_cell_count(); - - // Calculate weight for each marker particle - // based on FLUID_DENSITY and N_CELL properties contained in same particle. - for (int ic = 0; ic < num_cells; ic++) { - REAL V_cell = neso_mesh->dmh->get_cell_volume(ic); - particle_loop( - "Update weight of ions", marker_group, - [=](auto n_cell_prop, auto ion_dens_prop, auto weight_prop) { - const BoutReal n_cell = static_cast(n_cell_prop.at(0)); - auto updated_weight = (ion_dens_prop.at(0) * V_cell) / (N_w * n_cell); - weight_prop.at(0) = updated_weight; - }, - Access::read(Sym("N_CELL")), Access::read(Sym("FLUID_DENSITY")), - Access::write(Sym("WEIGHT"))) - ->execute(ic); - } - // Wrappers & controllers // ------------------------------------------------------------------------------ - VantageSourceManager source_manager(neso_mesh, bout_mesh, units); + this->source_manager = + std::make_unique(neso_mesh, bout_mesh, units); const REAL remove_threshold = options["remove_threshold"].withDefault(1.0e-10); const REAL merge_threshold = options["merge_threshold"].withDefault(1.0e-2); @@ -810,10 +904,12 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) std::vector> parent_transforms_iz = std::vector{iz_accumulator_real_transform_wrapper, remove_wrapper, merge_wrapper}; - auto reaction_controller = ReactionController(parent_transforms_iz, child_transforms); + this->reaction_controller = + std::make_unique(parent_transforms_iz, child_transforms); - source_manager.add_source("Siz", "ION_SOURCE_DENSITY", accumulator_transform_iz, - A_particle_group, ion_source_density_zeroer); + this->source_manager->add_source("Siz", "ION_SOURCE_DENSITY", + accumulator_transform_iz, A_particle_group, + ion_source_density_zeroer); // Recombination transforms and controller @@ -826,16 +922,17 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) std::vector> parent_transforms_rec = std::vector{recomb_accumulator_transform_wrapper, remove_wrapper, merge_wrapper}; - auto recombination_controller = - ReactionController(parent_transforms_rec, child_transforms); + this->recombination_controller = + std::make_unique(parent_transforms_rec, child_transforms); - source_manager.add_source("Srec", "ION_SOURCE_DENSITY", accumulator_transform_rec, - marker_group, ion_source_density_zeroer); + this->source_manager->add_source("Srec", "ION_SOURCE_DENSITY", + accumulator_transform_rec, marker_group, + ion_source_density_zeroer); // Ionisation reaction // ------------------------------------------------------------------------------ main_species.set_id(0); - + // Reaction rates // --------------------------- @@ -847,7 +944,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) A_particle_group->sycl_target, iz_rate_data, iz_rate_data, main_species, electron_species); - reaction_controller.add_reaction( + this->reaction_controller->add_reaction( std::make_shared(ionisation_reaction)); } else { @@ -856,7 +953,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) const hermes::AmjuelData iz_energy_rate_amjuel("H.10_2.1.5", alloptions); auto iz_rate_coeffs = convert_amjuel_format(iz_rate_amjuel.get_coeffs()); auto iz_energy_rate_coeffs = - convert_amjuel_format(iz_energy_rate_amjuel.get_coeffs()); + convert_amjuel_format(iz_energy_rate_amjuel.get_coeffs()); // Remap names: ionisation reaction expects "FLUID_TEMPERATURE" etc. auto ionisation_rate_map = get_default_map(); @@ -874,7 +971,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) A_particle_group->sycl_target, iz_rate_data, iz_energy_rate_data, main_species, electron_species, ionisation_rate_map); - reaction_controller.add_reaction( + this->reaction_controller->add_reaction( std::make_shared(ionisation_reaction)); } @@ -918,11 +1015,10 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) sycl_target, static_cast(rec_marker_species.get_id()), rec_out_states, rec_data, rec_reaction_kernel, rec_data_calc_obj); - recombination_controller.add_reaction( - std::make_shared(rec_reaction)); + this->recombination_controller->add_reaction( + std::make_shared(rec_reaction)); - } - else { + } else { // AMJUEL derived rate and energy rate const hermes::AmjuelData rec_rate_amjuel("H.4_2.1.8", alloptions); @@ -955,12 +1051,10 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) sycl_target, static_cast(rec_marker_species.get_id()), rec_out_states, rec_data, rec_reaction_kernel, rec_data_calc_obj); - recombination_controller.add_reaction( - std::make_shared(rec_reaction)); + this->recombination_controller->add_reaction( + std::make_shared(rec_reaction)); } - - // Boundary handling // ------------------------------------------------------------------------------ @@ -971,141 +1065,203 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) b2d = std::make_shared(sycl_target, neso_mesh, boundary_groups); - auto reflection = std::make_shared(ndim, 1.0e-10); + reflection = std::make_shared(ndim, 1.0e-10); + }; + + // Initialisation + // ------------------------------------------------------------------------------ + // allocate buffer vector for scalar projection/evaluation of NESO-Particles + // properties. Resize from default construction. + h_project1.resize(static_cast(num_cells_owned)); + // set weights from a Field2D from BOUT + set_initial_particle_weights(initial_neutral_density, dg0, A_particle_group, neso_mesh, + h_project1, N_w); + + // Calculate neutral density and sources for initial condition + calculate_neutral_density_in_place(neutral_density, dg0, A_particle_group, h_project1, + N_w); + this->source_manager->update_all_sources(dt); + + // diagnose the initial condition + bout_output_data = + initialise_diagnostics(alloptions, bout_mesh, neutral_density, ion_density, + neso_mesh, vantage_dump_filepath); + + // Object for VANTAGE dump files + vantage_dump_writer = + bout::OptionsIO::create({{"file", vantage_dump_filepath}, {"append", true}}); + + // Object for particle_trajectories.h5part file. + // Close it straight away to ensure that it's closed in event of a crash. + // h5part->write re-opens it when needed. + h5part = std::make_shared(particle_data_filepath, A_particle_group, + Sym("POSITION"), Sym("VELOCITY")); + h5part->close(); + + // mass for conservation check + total_density = neutral_density + ion_density; + total_mass_initial = calculate_total_mass(total_density, neso_mesh); + + // Initialise particle time + particle_time = 0.0; + + // Register VANTAGE timestep scheduler. + // https://bout-dev.readthedocs.io/en/latest/user_docs/time_integration.html#monitoring-the-simulation-output + // By default, it runs before the dump is written. (FRONT mode) + // If statement prevents segfault in the unit test where there is no solver. + if (solver != nullptr) { + solver->addMonitor(&this->monitor, Solver::FRONT); + } +} - auto lambda_apply_boundary_conditions = [&](auto aa) { - auto sub_groups = b2d->post_integration(aa); - for (auto& gx : sub_groups) { - reflection->execute(gx.second, Sym("POSITION"), Sym("VELOCITY"), - Sym("TSP"), b2d->previous_position_sym); - } - }; +void Vantage::apply_boundary_conditions(ParticleSubGroupSharedPtr aa) { + auto sub_groups = b2d->post_integration(aa); + for (auto& gx : sub_groups) { + reflection->execute(gx.second, Sym("POSITION"), Sym("VELOCITY"), + Sym("TSP"), b2d->previous_position_sym); + } +} - // Advection & rest of code - // ------------------------------------------------------------------------------ +int VantageMonitor::call(Solver* UNUSED(solver), BoutReal time, int iter, + int UNUSED(nout)) { + + // Prevent monitor from firing at the 0th timestep (before first timestep). + // Initialisation is ran in the constructor already. + if (iter == 0) { + return 0; + } + return vantage->advance_vantage(time); +} + +// Function called by the Monitor to advance kinetic neutrals for some number of VANTAGE timesteps +int Vantage::advance_vantage(BoutReal UNUSED(time)) { + + // Send plasma data to VANTAGE grid + if (plasma_coupling) { + send_plasma_data(plasma_data, sycl_target, bout_mesh, num_cells_owned, + background_ion_density, background_ion_temperature); + } + + // Send plasma data to recombination markers + update_recombination_markers(marker_group, plasma_data, neso_mesh, N_w, V_background); - auto lambda_apply_timestep_reset = [&](auto aa) { - particle_loop( - aa, - [=](auto TSP) { - TSP.at(0) = 0.0; - TSP.at(1) = 0.0; - }, - Access::write(Sym("TSP"))) - ->execute(); - }; - auto lambda_apply_advection_step = - [=](ParticleSubGroupSharedPtr iteration_set) -> void { - particle_loop( - "euler_advection", iteration_set, - [=](auto VELOCITY, auto POSITION, auto TSP) { - const REAL dt_left = dt - TSP.at(0); - if (dt_left > 0.0) { - POSITION.at(0) += dt_left * VELOCITY.at(0); - POSITION.at(1) += dt_left * VELOCITY.at(1); - TSP.at(0) = dt; - TSP.at(1) = dt_left; - } - }, - Access::read(Sym("VELOCITY")), Access::write(Sym("POSITION")), - Access::write(Sym("TSP"))) - ->execute(); - }; - auto lambda_pre_advection = [&](auto aa) { b2d->pre_integration(aa); }; - auto lambda_find_partial_moves = [&](auto aa) { - return static_particle_sub_group( - aa, [=](auto TSP) { return TSP.at(0) < dt; }, Access::read(Sym("TSP"))); - }; - auto lambda_partial_moves_remaining = [&](auto aa) -> bool { - const int size = static_cast(get_npart_global(aa)); - ; - return size > 0; - }; - auto lambda_apply_timestep = [&](auto aa) { - lambda_apply_timestep_reset(aa); + // Advection & rest of code + // ------------------------------------------------------------------------------ + + auto lambda_apply_timestep_reset = [&](auto aa) { + particle_loop( + aa, + [=](auto TSP) { + TSP.at(0) = 0.0; + TSP.at(1) = 0.0; + }, + Access::write(Sym("TSP"))) + ->execute(); + }; + + const BoutReal dt_local = this->dt; + + auto lambda_apply_advection_step = + [=](ParticleSubGroupSharedPtr iteration_set) -> void { + particle_loop( + "euler_advection", iteration_set, + [dt_local](auto VELOCITY, auto POSITION, auto TSP) { + const REAL dt_left = dt_local - TSP.at(0); + if (dt_left > 0.0) { + POSITION.at(0) += dt_left * VELOCITY.at(0); + POSITION.at(1) += dt_left * VELOCITY.at(1); + TSP.at(0) = dt_local; + TSP.at(1) = dt_left; + } + }, + Access::read(Sym("VELOCITY")), Access::write(Sym("POSITION")), + Access::write(Sym("TSP"))) + ->execute(); + }; + auto lambda_pre_advection = [&](auto aa) { b2d->pre_integration(aa); }; + auto lambda_find_partial_moves = [&](auto aa) { + return static_particle_sub_group( + aa, [dt_local](auto TSP) { return TSP.at(0) < dt_local; }, + Access::read(Sym("TSP"))); + }; + auto lambda_partial_moves_remaining = [&](auto aa) -> bool { + const int size = static_cast(get_npart_global(aa)); + ; + return size > 0; + }; + auto lambda_apply_timestep = [&](auto aa) { + lambda_apply_timestep_reset(aa); + lambda_pre_advection(aa); + lambda_apply_advection_step(aa); + this->apply_boundary_conditions(aa); + aa = lambda_find_partial_moves(aa); + while (lambda_partial_moves_remaining(aa)) { lambda_pre_advection(aa); lambda_apply_advection_step(aa); - lambda_apply_boundary_conditions(aa); + this->apply_boundary_conditions(aa); aa = lambda_find_partial_moves(aa); - while (lambda_partial_moves_remaining(aa)) { - lambda_pre_advection(aa); - lambda_apply_advection_step(aa); - lambda_apply_boundary_conditions(aa); - aa = lambda_find_partial_moves(aa); - } - }; - - // allocate buffer vector for scalar projection/evaluation of NESO-Particles - // properties - std::vector h_project1(static_cast(num_cells_owned)); - // set weights from a Field2D from BOUT - set_initial_particle_weights(initial_neutral_density, dg0, A_particle_group, - neso_mesh, h_project1, N_w); - - // Calculate neutral density and sources for initial condition - calculate_neutral_density_in_place(neutral_density, dg0, A_particle_group, - h_project1, N_w); - source_manager.update_all_sources(dt); - - // diagnose the initial condition - std::string particle_data_filename = - make_output_path(fmt::format("BOUT.dmp.vantage.{}.nc", mpi_rank), alloptions); - Options bout_output_data = initialise_diagnostics( - alloptions, bout_mesh, neutral_density, ion_density, neso_mesh, particle_data_filename); - // mass for conservation check - Field2D total_density = neutral_density + ion_density; - double total_mass_initial = calculate_total_mass(total_density, neso_mesh); - - // Initialise h5part just before writing - earlier leads to a NESO assert error that - // can hide other bugs - auto h5part = std::make_shared( - make_output_path("particle_trajectories.h5part", alloptions), A_particle_group, - Sym("POSITION"), Sym("VELOCITY")); - - // begin timestepping - for (int stepx = 0; stepx < nsteps; stepx++) { - // nprint("step:", stepx); - output << "step:" << std::to_string(stepx) << std::endl; - particle_time += dt; - A_particle_group->hybrid_move(); - A_particle_group->cell_move(); - lambda_apply_timestep(static_particle_sub_group(A_particle_group)); - // apply reactions - reaction_controller.apply(A_particle_group, dt, ControllerMode::standard_mode); - recombination_controller.apply(marker_group, dt, A_particle_group); - // uncomment to write a trajectory - h5part->write(); - - calculate_neutral_density_in_place(neutral_density, dg0, A_particle_group, - h_project1, N_w); - source_manager.update_all_sources(dt); - Field2D Siz = source_manager.get_data("Siz"); - Field2D Srec = source_manager.get_data("Srec"); - - // "Solve" density - // Sources are in normalised m^-3 s^-1, so need to multiply by dt - ion_density += (Siz + Srec) * dt; - - // diagnose timestep stepx - update_diagnostics(neutral_density, ion_density, - Siz, Srec, - neso_mesh, bout_output_data, - particle_data_filename, particle_time); } - h5part->close(); + }; + + // begin timestepping + output << "\nBegin VANTAGE iterations \n"; + for (int stepx = 0; stepx < nsteps; stepx++) { - // mass for conservation check - total_density = neutral_density + ion_density; - double total_mass_final = calculate_total_mass(total_density, neso_mesh); - if (options["test_mass_conservation"].withDefault(true)) { - check_mass_conservation(total_mass_final, total_mass_initial); + output << "Particle time: " << std::to_string(particle_time) << std::endl; + + particle_time += dt; + A_particle_group->hybrid_move(); + A_particle_group->cell_move(); + lambda_apply_timestep(static_particle_sub_group(A_particle_group)); + + if (plasma_coupling) { + update_particle_background(A_particle_group, plasma_data); } + + // apply reactions + reaction_controller->apply(A_particle_group, dt, ControllerMode::standard_mode); + recombination_controller->apply(marker_group, dt, A_particle_group); + + calculate_neutral_density_in_place(neutral_density, dg0, A_particle_group, h_project1, + N_w); + this->source_manager->update_all_sources(dt); + Field2D Siz = this->source_manager->get_data("Siz"); + Field2D Srec = this->source_manager->get_data("Srec"); + + // "Solve" density + // Sources are in normalised m^-3 s^-1, so need to multiply by dt + ion_density += (Siz + Srec) * dt; + + // Write to VANTAGE dump files + update_diagnostics(neutral_density, ion_density, Siz, Srec, neso_mesh, + bout_output_data, *vantage_dump_writer, particle_time); + + // Write to particle_trajectories file + h5part->write(); + } + // Warning: if h5part gets destroyed while open due to crash, you will get a NESO_ASSERT + // warning which will mask the actual backtrace. In this event + // you could move this into the loop, but it will have an IO cost. + h5part->close(); + + // mass for conservation check + total_density = neutral_density + ion_density; + BoutReal total_mass_final = calculate_total_mass(total_density, neso_mesh); + if (test_mass_conservation) { + check_mass_conservation(total_mass_final, total_mass_initial); } + return 0; } void Vantage::transform_impl(GuardedOptions& UNUSED(state)) {} -void Vantage::finally(const Options& UNUSED(state)) {} +void Vantage::finally(const Options& state) { + + // Do not read from state if VANTAGE running standalone + if (!plasma_coupling) { + return; + } +} void Vantage::outputVars(Options& UNUSED(state)) {} diff --git a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp index a18f2af4f..474d91051 100644 --- a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp +++ b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp @@ -1,5 +1,6 @@ -nout = 0 +nout = 1 timestep = 1 +MZ = 1 [mesh] file = example_usn.grd.nc @@ -25,6 +26,7 @@ charge = 1 [Nd+] function = 1 +bndry_all = neumann [e] type = quasineutral @@ -35,6 +37,9 @@ charge = -1 dt = 0.005 nsteps = 3 +ion_species = d+ +neutral_species = d + test_mass_conservation = true initial_neutral_density = 20 diff --git a/tests/integrated/particle-pusher/runtest b/tests/integrated/particle-pusher/runtest index 18682d787..6129460ab 100755 --- a/tests/integrated/particle-pusher/runtest +++ b/tests/integrated/particle-pusher/runtest @@ -6,7 +6,6 @@ from boututils.run_wrapper import shell, launch_safe import h5py import numpy as np -from petsc4py import PETSc import netCDF4 as nc verbose = False @@ -30,14 +29,15 @@ def particle_push_input( ): input_file_string = f""" - nout = 0 + nout = 1 timestep = 1 + MZ = 1 [mesh] J = 1 - nx = 40 # X grid size - ny = 36 # Y grid size + nx = {nx} # X grid size + ny = {ny} # Y grid size dx = 1.0/(nx-4) # X mesh spacing dy = 2*pi/ny # Y mesh spacing @@ -73,6 +73,7 @@ def particle_push_input( [Nd+] function = 1 + bndry_all = neumann [e] type = quasineutral @@ -82,6 +83,8 @@ def particle_push_input( [vantage] dt = {dt} nsteps = {nsteps} + neutral_species = d + ion_species = d+ test_mass_conservation = true initial_neutral_density = 1e19 npart_per_cell = 20 diff --git a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp index 286176cc3..bf1a25c63 100644 --- a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp +++ b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp @@ -1,5 +1,6 @@ -nout = 0 +nout = 1 timestep = 1 +MZ = 1 [mesh] J = 1 @@ -41,6 +42,7 @@ charge = 1 [Nd+] function = 1 +bndry_all = neumann [e] type = quasineutral @@ -49,11 +51,14 @@ charge = -1 [vantage] # Trying to get as close to equilibrium in only 3 steps. -# This is not straightforward because the step can't be too big +# This is not straightforward because the step can't be too big # so that it consumes all of the neutral particles in one go. dt = 700 nsteps = 3 +ion_species = d+ +neutral_species = d + # This is the final expected steady state density initial_neutral_density = 8e18 @@ -69,4 +74,4 @@ background_ion_Vy = 0 remove_threshold = 1e-10 merge_threshold = 0 iz_rate_override = -1 -rec_rate_override = -1 \ No newline at end of file +rec_rate_override = -1 diff --git a/tests/unit/test_vantage.cxx b/tests/unit/test_vantage.cxx index fdb25eb09..d6853b65a 100644 --- a/tests/unit/test_vantage.cxx +++ b/tests/unit/test_vantage.cxx @@ -20,11 +20,14 @@ extern Mesh* mesh; // The unit tests use the global mesh using namespace bout::globals; +// Create a FakeMeshFixture in 2D only +using AxisymmetricMeshFixture = FakeMeshFixture_tmpl<3, 5, 1>; + // Constructor for VantageTest which inherits from FakeMeshFixture. // Construct a blank options and add fields to the mesh. -class VantageTest : public FakeMeshFixture { +class VantageTest : public AxisymmetricMeshFixture { public: - VantageTest() : FakeMeshFixture() {} + VantageTest() : AxisymmetricMeshFixture() {} Options alloptions; protected: @@ -100,6 +103,8 @@ Options MakeOptions() { alloptions["vantage"]["npart_per_cell"] = 1; alloptions["vantage"]["nsteps"] = 0; alloptions["vantage"]["dt"] = 0.001; + + alloptions["vantage"]["neutral_species"] = "d"; return alloptions; } @@ -109,5 +114,5 @@ Options MakeOptions() { TEST_F(VantageTest, CreateComponent) { Options alloptions = MakeOptions(); - Vantage const component("vantage", alloptions, nullptr); -} \ No newline at end of file + const Vantage component("vantage", alloptions, nullptr); +}