From 7f28234edfc5d3d68b2df6c70f1d483e41c7a0d5 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 22 Jul 2026 11:48:51 +0100 Subject: [PATCH 01/10] Move particle pushing to transform Had to move a lot of things into private members of the class. Also tidied up paths and some types. Particle advection works but is now called per RHS evaluation which means they will get called twice in the tests causing them to fail. --- include/vantage.hxx | 83 ++++++++---- src/vantage.cxx | 320 +++++++++++++++++++++++--------------------- 2 files changed, 227 insertions(+), 176 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index e4586649e..6d6694d36 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,64 @@ private: Options& units; }; +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: + + bool test_mass_conservation; + BoutReal particle_time; + 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 + 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 + + + + PetscLib petsc_lib; // Ensures PETSc is initialized for the lifetime of this component + + 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::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 + + // Needed for Vantage::apply_boundary_conditions + std::shared_ptrreflection; // 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 + 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..8803a9b2d 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -85,10 +85,10 @@ void calculate_neutral_density_in_place( // extrapolate -> Neumann } -double calculate_total_mass(Field2D& density, +BoutReal calculate_total_mass(Field2D& density, std::shared_ptr& neso_mesh) { - double local_mass = 0.0; - double total_mass = 0.0; + 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++) { @@ -107,7 +107,7 @@ initialise_diagnostics(Options& alloptions, Mesh* bout_mesh, Field2D& neutral_density, Field2D& ion_density, std::shared_ptr& neso_mesh, - std::string particle_data_filename) { + 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"]); @@ -209,14 +209,14 @@ initialise_diagnostics(Options& alloptions, {"long_name", "Gyro-radius length normalisation"} }); - bout::OptionsIO::create(particle_data_filename)->write(bout_output_data); + 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, std::shared_ptr& neso_mesh, - Options& bout_output_data, std::string particle_data_filename, + Options& bout_output_data, std::string vantage_dump_filepath, BoutReal particle_time) { // update density in Options object and write bout_output_data["neutral_density"] = neutral_density; @@ -232,7 +232,7 @@ 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}}) + bout::OptionsIO::create({{"file", vantage_dump_filepath}, {"append", true}}) ->write(bout_output_data); } @@ -375,9 +375,9 @@ void check_cell_centres(Options& alloptions, std::shared_ptr(units["meters"]); BoutReal seconds = get(units["seconds"]); - BoutReal N_w = options["N_w"] + 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.") @@ -491,18 +491,26 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) Options::root()["units"]["N_w"] = N_w; Options::root()["units"]["N_w"].setConditionallyUsed(); - Mesh* bout_mesh = bout::globals::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 @@ -529,11 +537,14 @@ 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( @@ -579,17 +590,16 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // Other settings const int ndim = 2; - const REAL dt = options["dt"] + dt = options["dt"] .doc("Timestep to use for VANTAGE kinetic neutrals (normalised units)") .withDefault(0.01); - const int nsteps = options["nsteps"] + 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 +610,7 @@ 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(); // 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)) { @@ -638,11 +648,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; @@ -689,7 +699,7 @@ 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( + dg0 = std::make_shared( neso_mesh, sycl_target, "DG", 0); // RNG kernel @@ -710,7 +720,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 = + marker_group = std::make_shared(domain, particle_spec, sycl_target); // Give particle group initial kinetic values (positions and velocities) @@ -776,7 +786,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // 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,9 +820,9 @@ 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, + this->source_manager->add_source("Siz", "ION_SOURCE_DENSITY", accumulator_transform_iz, A_particle_group, ion_source_density_zeroer); // Recombination transforms and controller @@ -826,10 +836,10 @@ 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, + this->source_manager->add_source("Srec", "ION_SOURCE_DENSITY", accumulator_transform_rec, marker_group, ion_source_density_zeroer); // Ionisation reaction @@ -847,7 +857,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 { @@ -874,7 +884,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,7 +928,7 @@ 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( + this->recombination_controller->add_reaction( std::make_shared(rec_reaction)); } @@ -955,7 +965,7 @@ 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( + this->recombination_controller->add_reaction( std::make_shared(rec_reaction)); } @@ -971,139 +981,147 @@ 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); - 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); - } + }; - // Advection & rest of code + // Initialisation // ------------------------------------------------------------------------------ - - 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); - lambda_pre_advection(aa); - lambda_apply_advection_step(aa); - lambda_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)); + // 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); + 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); + 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 - 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); + bout_output_data = + initialise_diagnostics(alloptions, bout_mesh, neutral_density, ion_density, + neso_mesh, vantage_dump_filepath); // 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); + total_density = neutral_density + ion_density; + total_mass_initial = calculate_total_mass(total_density, neso_mesh); + + // Initialise particle time + particle_time = 0.0; + } - h5part->close(); - // 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); + +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); + } +} + +void Vantage::transform_impl(GuardedOptions& UNUSED(state)) { + // 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(); + }; + 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); + 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); + this->apply_boundary_conditions(aa); + aa = lambda_find_partial_moves(aa); } + }; + + + + // Initialise h5part just before writing - earlier leads to a NESO assert error that + // can hide other bugs + auto h5part = std::make_shared( + particle_data_filepath, 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); + 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; + + // diagnose timestep stepx + update_diagnostics(neutral_density, ion_density, Siz, Srec, neso_mesh, + bout_output_data, vantage_dump_filepath, particle_time); } -} + h5part->close(); -void Vantage::transform_impl(GuardedOptions& UNUSED(state)) {} + // 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); + } +} void Vantage::finally(const Options& UNUSED(state)) {} From 43eb1a1853619000568606b0ade83d2e3398e65e Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 23 Jul 2026 14:53:33 +0100 Subject: [PATCH 02/10] Use monitor to schedule VANTAGE calls Now the VANTAGE push is triggered every RHS evaluation and no longer runs in the constructor. Updated tests - now we run 1 real timestep. Had to change BCs to Neumann to prevent the fluid side from crashing. --- include/vantage.hxx | 58 ++-- src/vantage.cxx | 295 +++++++++--------- .../dmplex-vertex-coordinates/data/BOUT.inp | 3 +- tests/integrated/particle-pusher/runtest | 4 +- .../vantage-iz-rec-balance/data/BOUT.inp | 7 +- 5 files changed, 194 insertions(+), 173 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index 6d6694d36..2d8080d0c 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -61,6 +61,22 @@ private: Options& units; }; +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); @@ -69,56 +85,58 @@ struct Vantage : public Component { void transform_impl(GuardedOptions& state) override; void outputVars(Options& state) override; -private: + // Function with the kinetic loop. + // time is the normalised VANTAGE monitor frequency. + int advance_vantage(BoutReal time); +private: bool test_mass_conservation; BoutReal particle_time; BoutReal N_w; REAL dt; int nsteps; - int num_cells_owned; // Number of VANTAGE cells owned per rank + int num_cells_owned; // Number of VANTAGE cells owned per rank Options bout_output_data; // Options object to hold output data for VANTAGE diagnostics - int mpi_rank; // Current rank ID - Mesh* bout_mesh; // Pointer to the BOUT++ mesh object + 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 + 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 - - + std::string dmplex_filepath, vantage_dump_filepath, + particle_data_filepath; // Path for output files PetscLib petsc_lib; // Ensures PETSc is initialized for the lifetime of this component - + 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::vector h_project1; // Buffer for scalar projection/evaluation of NESO-Particles properties + std::shared_ptr + b2d; // Boundary interaction object + std::shared_ptr + dg0; // DMPlex projection 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 // Needed for Vantage::apply_boundary_conditions - std::shared_ptrreflection; // Boundary reflection object + 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 + 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 8803a9b2d..f47a20e9d 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,8 +85,9 @@ void calculate_neutral_density_in_place( // extrapolate -> Neumann } -BoutReal calculate_total_mass(Field2D& density, - std::shared_ptr& neso_mesh) { +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(); @@ -103,16 +104,15 @@ BoutReal 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 vantage_dump_filepath) { +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,49 +170,43 @@ 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"} - }); + 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 vantage_dump_filepath, BoutReal particle_time) { @@ -241,8 +233,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 +241,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++; } @@ -316,7 +306,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 @@ -454,7 +445,7 @@ void VantageSourceManager::update_all_sources(double dt) { } } -Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) +Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) : Component({readOnly("species:d+:density", Regions::Interior), readWrite("species:d+:density")}) { @@ -472,9 +463,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"]); @@ -483,10 +473,10 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) BoutReal seconds = get(units["seconds"]); 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); + .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(); @@ -504,13 +494,14 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) 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); + 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, dmplex_filepath); + + 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 @@ -518,7 +509,7 @@ 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"; @@ -537,20 +528,20 @@ 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); + test_mass_conservation = + options["test_mass_conservation"] + .doc("Check that mass is conserved at runtime. Default = true") + .withDefault(true); // Normalisations // 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") @@ -585,17 +576,16 @@ 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 const int ndim = 2; dt = options["dt"] - .doc("Timestep to use for VANTAGE kinetic neutrals (normalised units)") - .withDefault(0.01); + .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); + .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); @@ -618,8 +608,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)); } @@ -652,7 +641,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) std::make_shared(domain, particle_spec, sycl_target); // Create some particle data - + std::mt19937 rng_pos(static_cast(52234234 + mpi_rank)); std::mt19937 rng_vel(static_cast(52234231 + mpi_rank)); std::vector> positions; @@ -684,8 +673,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++) { @@ -699,8 +690,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 - 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 @@ -709,8 +700,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 @@ -720,8 +709,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // Options and constants // Make new particle group just for the markers - 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 @@ -786,7 +774,8 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) // Wrappers & controllers // ------------------------------------------------------------------------------ - this->source_manager = std::make_unique(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); @@ -820,10 +809,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}; - this->reaction_controller = std::make_unique(parent_transforms_iz, child_transforms); + this->reaction_controller = + std::make_unique(parent_transforms_iz, child_transforms); - this->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 @@ -839,13 +830,14 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) this->recombination_controller = std::make_unique(parent_transforms_rec, child_transforms); - this->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 // --------------------------- @@ -866,7 +858,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(); @@ -929,10 +921,9 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) rec_data, rec_reaction_kernel, rec_data_calc_obj); this->recombination_controller->add_reaction( - std::make_shared(rec_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); @@ -966,11 +957,9 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) rec_data, rec_reaction_kernel, rec_data_calc_obj); this->recombination_controller->add_reaction( - std::make_shared(rec_reaction)); + std::make_shared(rec_reaction)); } - - // Boundary handling // ------------------------------------------------------------------------------ @@ -982,47 +971,59 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* UNUSED(solver)) b2d = std::make_shared(sycl_target, neso_mesh, boundary_groups); 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); + // 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); + // mass for conservation check + total_density = neutral_density + ion_density; + total_mass_initial = calculate_total_mass(total_density, neso_mesh); - // diagnose the initial condition - bout_output_data = - initialise_diagnostics(alloptions, bout_mesh, neutral_density, ion_density, - neso_mesh, vantage_dump_filepath); - // 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; - // Initialise particle time - particle_time = 0.0; + // Register VANTAGE timestep scheduler. + // By default, it runs before the dump is written. (FRONT mode) + solver->addMonitor(&this->monitor, Solver::FRONT); +} - } +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); + } +} +int VantageMonitor::call(Solver* UNUSED(solver), BoutReal time, int iter, + int UNUSED(nout)) { -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); - } + // 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); } -void Vantage::transform_impl(GuardedOptions& UNUSED(state)) { +// Function called by the Monitor to advance kinetic neutrals for some number of VANTAGE timesteps +int Vantage::advance_vantage(BoutReal UNUSED(time)) { // Advection & rest of code // ------------------------------------------------------------------------------ @@ -1077,13 +1078,10 @@ void Vantage::transform_impl(GuardedOptions& UNUSED(state)) { } }; - - // Initialise h5part just before writing - earlier leads to a NESO assert error that // can hide other bugs - auto h5part = std::make_shared( - particle_data_filepath, A_particle_group, - Sym("POSITION"), Sym("VELOCITY")); + auto h5part = std::make_shared(particle_data_filepath, A_particle_group, + Sym("POSITION"), Sym("VELOCITY")); // begin timestepping for (int stepx = 0; stepx < nsteps; stepx++) { @@ -1121,8 +1119,11 @@ void Vantage::transform_impl(GuardedOptions& UNUSED(state)) { 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::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..b238b427d 100644 --- a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp +++ b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp @@ -1,4 +1,4 @@ -nout = 0 +nout = 1 timestep = 1 [mesh] @@ -25,6 +25,7 @@ charge = 1 [Nd+] function = 1 +bndry_all = neumann [e] type = quasineutral diff --git a/tests/integrated/particle-pusher/runtest b/tests/integrated/particle-pusher/runtest index 18682d787..a754e44fb 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,7 +29,7 @@ def particle_push_input( ): input_file_string = f""" - nout = 0 + nout = 1 timestep = 1 [mesh] @@ -73,6 +72,7 @@ def particle_push_input( [Nd+] function = 1 + bndry_all = neumann [e] type = quasineutral diff --git a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp index 286176cc3..4e77ae287 100644 --- a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp +++ b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp @@ -1,4 +1,4 @@ -nout = 0 +nout = 1 timestep = 1 [mesh] @@ -41,6 +41,7 @@ charge = 1 [Nd+] function = 1 +bndry_all = neumann [e] type = quasineutral @@ -49,7 +50,7 @@ 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 @@ -69,4 +70,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 From 6cc38296bc62ccbdeed65affc58c0e936b05d013 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 23 Jul 2026 15:38:22 +0100 Subject: [PATCH 03/10] Prevent particle output being overwritten each step Previously h5part would be recreated every RHS. --- include/vantage.hxx | 1 + src/vantage.cxx | 8 +++++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index 2d8080d0c..3e1996c5b 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -115,6 +115,7 @@ private: 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 diff --git a/src/vantage.cxx b/src/vantage.cxx index f47a20e9d..2ba54b254 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -1079,9 +1079,11 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { }; // Initialise h5part just before writing - earlier leads to a NESO assert error that - // can hide other bugs - auto h5part = std::make_shared(particle_data_filepath, A_particle_group, - Sym("POSITION"), Sym("VELOCITY")); + // can hide other bugs if it crashes with the file open. + if (!h5part) { + h5part = std::make_shared(particle_data_filepath, A_particle_group, + Sym("POSITION"), Sym("VELOCITY")); + } // begin timestepping for (int stepx = 0; stepx < nsteps; stepx++) { From 57bb928915f08adfbd191ce47ca9b99afc9c9c83 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 23 Jul 2026 16:30:43 +0100 Subject: [PATCH 04/10] Tidy VANTAGE console output Previously the IO object would be created every step and spammed the console. Now we create it once then write and flush as needed. Also made a slightly nicer timestep printout. --- include/vantage.hxx | 9 ++++++--- src/vantage.cxx | 18 +++++++++++------- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index 3e1996c5b..847ad7c1f 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -98,8 +98,11 @@ private: int num_cells_owned; // Number of VANTAGE cells owned per rank Options bout_output_data; // Options object to hold output data for VANTAGE diagnostics - int mpi_rank; // Current rank ID - Mesh* bout_mesh; // Pointer to the BOUT++ mesh object + 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; @@ -114,7 +117,7 @@ private: std::shared_ptr b2d; // Boundary interaction object std::shared_ptr - dg0; // DMPlex projection object + 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 diff --git a/src/vantage.cxx b/src/vantage.cxx index 2ba54b254..31a3ec08e 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -208,7 +208,7 @@ initialise_diagnostics(Options& alloptions, Mesh* bout_mesh, Field2D& neutral_de void update_diagnostics(Field2D& neutral_density, Field2D& ion_density, Field2D& Siz, Field2D& Srec, std::shared_ptr& neso_mesh, - Options& bout_output_data, std::string vantage_dump_filepath, + 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; @@ -224,8 +224,9 @@ void update_diagnostics(Field2D& neutral_density, Field2D& ion_density, Field2D& bout_output_data["t_array"] = particle_time; // bout_output_data["t_array"] = 0.0; // Append data to file - bout::OptionsIO::create({{"file", vantage_dump_filepath}, {"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( @@ -512,7 +513,6 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) PETSCCHK(VecScale(coords, 1 / meters)); } - output << "Begin particle push \n"; /* * @@ -991,6 +991,9 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) bout_output_data = initialise_diagnostics(alloptions, bout_mesh, neutral_density, ion_density, neso_mesh, vantage_dump_filepath); + + vantage_dump_writer = bout::OptionsIO::create({{"file", vantage_dump_filepath}, {"append", true}}); + // mass for conservation check total_density = neutral_density + ion_density; total_mass_initial = calculate_total_mass(total_density, neso_mesh); @@ -1086,9 +1089,10 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { } // begin timestepping + output << "\nBegin VANTAGE iterations \n"; for (int stepx = 0; stepx < nsteps; stepx++) { - // nprint("step:", stepx); - output << "step:" << std::to_string(stepx) << std::endl; + + output << "Particle time: " << std::to_string(particle_time) << std::endl; particle_time += dt; A_particle_group->hybrid_move(); A_particle_group->cell_move(); @@ -1111,7 +1115,7 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { // diagnose timestep stepx update_diagnostics(neutral_density, ion_density, Siz, Srec, neso_mesh, - bout_output_data, vantage_dump_filepath, particle_time); + bout_output_data, *vantage_dump_writer, particle_time); } h5part->close(); From 16c8f08c14e4ea1f3b46e904b3ab085a326a68cb Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 23 Jul 2026 16:32:02 +0100 Subject: [PATCH 05/10] Fix a thing deprecated in C++20 Lambda functions are no longer allowed to capture private members with "=". I added a lot of members to move the push out of the constructor, they now need to be explicitly made local. This follows NESO-Particles convention. --- src/vantage.cxx | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/vantage.cxx b/src/vantage.cxx index 31a3ec08e..05f2f712a 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -513,7 +513,6 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) PETSCCHK(VecScale(coords, 1 / meters)); } - /* * * @@ -754,6 +753,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) } const INT num_cells = marker_group->domain->mesh->get_cell_count(); + const BoutReal N_w_local = this->N_w; // Calculate weight for each marker particle // based on FLUID_DENSITY and N_CELL properties contained in same particle. @@ -761,9 +761,9 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) 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) { + [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 * n_cell); + 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")), @@ -992,7 +992,8 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) initialise_diagnostics(alloptions, bout_mesh, neutral_density, ion_density, neso_mesh, vantage_dump_filepath); - vantage_dump_writer = bout::OptionsIO::create({{"file", vantage_dump_filepath}, {"append", true}}); + vantage_dump_writer = + bout::OptionsIO::create({{"file", vantage_dump_filepath}, {"append", true}}); // mass for conservation check total_density = neutral_density + ion_density; @@ -1040,16 +1041,19 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { 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, - [=](auto VELOCITY, auto POSITION, auto TSP) { - const REAL dt_left = dt - TSP.at(0); + [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; + TSP.at(0) = dt_local; TSP.at(1) = dt_left; } }, @@ -1060,7 +1064,8 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { 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"))); + 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)); From 9d38a5e8f7c828e091a448492fa0dd979b65f248 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 19 Aug 2026 12:47:59 +0100 Subject: [PATCH 06/10] Address review comments - Extra comments for context - h5part now created and closed straight away in constructor, turns out it is automatically opened on-demand. - Unit test fixed by not adding a monitor when there is no solver present. --- include/vantage.hxx | 4 ++++ src/vantage.cxx | 31 ++++++++++++++++++++----------- 2 files changed, 24 insertions(+), 11 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index 847ad7c1f..dc9f08751 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -61,6 +61,10 @@ 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. diff --git a/src/vantage.cxx b/src/vantage.cxx index 05f2f712a..47210f450 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -992,9 +992,17 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) 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); @@ -1003,8 +1011,12 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) 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) - solver->addMonitor(&this->monitor, Solver::FRONT); + // If statement prevents segfault in the unit test where there is no solver. + if (solver != nullptr) { + solver->addMonitor(&this->monitor, Solver::FRONT); + } } void Vantage::apply_boundary_conditions(ParticleSubGroupSharedPtr aa) { @@ -1086,13 +1098,6 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { } }; - // Initialise h5part just before writing - earlier leads to a NESO assert error that - // can hide other bugs if it crashes with the file open. - if (!h5part) { - h5part = std::make_shared(particle_data_filepath, A_particle_group, - Sym("POSITION"), Sym("VELOCITY")); - } - // begin timestepping output << "\nBegin VANTAGE iterations \n"; for (int stepx = 0; stepx < nsteps; stepx++) { @@ -1105,8 +1110,6 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { // 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); @@ -1118,10 +1121,16 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { // Sources are in normalised m^-3 s^-1, so need to multiply by dt ion_density += (Siz + Srec) * dt; - // diagnose timestep stepx + // 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 From 5d613ef95c2c4989b0672fe843b87d2750e17a6d Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 3 Sep 2026 10:47:50 +0100 Subject: [PATCH 07/10] Enforce 2D domain VANTAGE supports only 2D domains at the moment. Our integrated tests were running in 3D by accident. Now an exception prevent this and all tests were updated to 2D. --- src/vantage.cxx | 14 +++++++++++--- .../dmplex-vertex-coordinates/data/BOUT.inp | 1 + tests/integrated/particle-pusher/runtest | 5 +++-- .../vantage-iz-rec-balance/data/BOUT.inp | 1 + tests/unit/test_vantage.cxx | 9 ++++++--- 5 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/vantage.cxx b/src/vantage.cxx index 47210f450..219bbcea0 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -450,6 +450,17 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) : Component({readOnly("species:d+:density", Regions::Interior), readWrite("species:d+:density")}) { + //Get BOUT++ mesh and comm, throw if mesh not 2D + bout_mesh = bout::globals::mesh; + sycl_target = std::make_shared(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) { @@ -482,9 +493,6 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) Options::root()["units"]["N_w"] = N_w; Options::root()["units"]["N_w"].setConditionallyUsed(); - 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"] diff --git a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp index b238b427d..39566ba60 100644 --- a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp +++ b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp @@ -1,5 +1,6 @@ nout = 1 timestep = 1 +MZ = 1 [mesh] file = example_usn.grd.nc diff --git a/tests/integrated/particle-pusher/runtest b/tests/integrated/particle-pusher/runtest index a754e44fb..1233a0e1e 100755 --- a/tests/integrated/particle-pusher/runtest +++ b/tests/integrated/particle-pusher/runtest @@ -31,12 +31,13 @@ def particle_push_input( 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 diff --git a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp index 4e77ae287..4d7b2aadc 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 = 1 timestep = 1 +MZ = 1 [mesh] J = 1 diff --git a/tests/unit/test_vantage.cxx b/tests/unit/test_vantage.cxx index fdb25eb09..9a28e818d 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: @@ -110,4 +113,4 @@ TEST_F(VantageTest, CreateComponent) { Options alloptions = MakeOptions(); Vantage const component("vantage", alloptions, nullptr); -} \ No newline at end of file +} From f68ce1161071f1680c832519809bea1d9c0ecee0 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 3 Sep 2026 10:57:40 +0100 Subject: [PATCH 08/10] Add facility to enable or disable plasma coupling The user must specify an ion species for VANTAGE plasma coupling to be enabled, and the status of this is printed to console at construction. We now also require the neutral species to be explicitly specified. Tests have been updated. --- include/vantage.hxx | 9 +++++-- src/vantage.cxx | 27 ++++++++++++++++++- .../dmplex-vertex-coordinates/data/BOUT.inp | 2 ++ tests/integrated/particle-pusher/runtest | 1 + .../vantage-iz-rec-balance/data/BOUT.inp | 2 ++ tests/unit/test_vantage.cxx | 4 ++- 6 files changed, 41 insertions(+), 4 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index dc9f08751..082ae2252 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -95,7 +95,6 @@ struct Vantage : public Component { private: bool test_mass_conservation; - BoutReal particle_time; BoutReal N_w; REAL dt; int nsteps; @@ -113,8 +112,14 @@ private: std::string dmplex_filepath, vantage_dump_filepath, particle_data_filepath; // Path for output files - PetscLib petsc_lib; // Ensures PETSc is initialized for the lifetime of this component + // 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 + // DMPLex/VANTAGE stuff DM dm; std::shared_ptr neso_mesh; std::shared_ptr sycl_target; diff --git a/src/vantage.cxx b/src/vantage.cxx index 219bbcea0..32f03cf9c 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -586,6 +586,25 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) 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; dt = options["dt"] .doc("Timestep to use for VANTAGE kinetic neutrals (normalised units)") @@ -1152,7 +1171,13 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { 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 39566ba60..fdb64693f 100644 --- a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp +++ b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp @@ -37,6 +37,8 @@ charge = -1 dt = 0.005 nsteps = 3 +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 1233a0e1e..0b0063bc9 100755 --- a/tests/integrated/particle-pusher/runtest +++ b/tests/integrated/particle-pusher/runtest @@ -83,6 +83,7 @@ def particle_push_input( [vantage] dt = {dt} nsteps = {nsteps} + neutral_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 4d7b2aadc..1b706ac1e 100644 --- a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp +++ b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp @@ -56,6 +56,8 @@ charge = -1 dt = 700 nsteps = 3 +neutral_species = d + # This is the final expected steady state density initial_neutral_density = 8e18 diff --git a/tests/unit/test_vantage.cxx b/tests/unit/test_vantage.cxx index 9a28e818d..d6853b65a 100644 --- a/tests/unit/test_vantage.cxx +++ b/tests/unit/test_vantage.cxx @@ -103,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; } @@ -112,5 +114,5 @@ Options MakeOptions() { TEST_F(VantageTest, CreateComponent) { Options alloptions = MakeOptions(); - Vantage const component("vantage", alloptions, nullptr); + const Vantage component("vantage", alloptions, nullptr); } From 8ac98cca079cd62ee57b848d8f4d99f32a7571ff Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Fri, 4 Sep 2026 14:52:43 +0100 Subject: [PATCH 09/10] Encapsulate recombination marker calculation Recombination markers must be made from up-to-date plasma properties. This is now a function which runs at the start of each VANTAGE call. --- include/vantage.hxx | 2 + src/vantage.cxx | 132 +++++++++++++++++++++++++------------------- 2 files changed, 77 insertions(+), 57 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index 082ae2252..74427e04b 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -118,6 +118,8 @@ private: 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; diff --git a/src/vantage.cxx b/src/vantage.cxx index 32f03cf9c..45c022f2a 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -258,6 +258,73 @@ void set_initial_particle_weights( dg0->evaluate(A_particle_group, Sym("WEIGHT")); } +// 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& neso_mesh, BoutReal N_w, + BoutReal background_ion_density, BoutReal background_ion_temperature, + 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 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 < nvel; 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(); + 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); + } +} + void check_cell_volumes(std::shared_ptr& neso_mesh, Mesh*& bout_mesh, Options& alloptions) { Coordinates* coord = bout_mesh->getCoordinates(); @@ -555,12 +622,12 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) .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) @@ -571,7 +638,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* 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 = @@ -744,60 +811,6 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* 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(); - const BoutReal N_w_local = this->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); - } - // Wrappers & controllers // ------------------------------------------------------------------------------ @@ -1067,6 +1080,11 @@ int VantageMonitor::call(Solver* UNUSED(solver), BoutReal time, int iter, // 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 recombination markers + update_recombination_markers(marker_group, neso_mesh, N_w, background_ion_density, + background_ion_temperature, V_background); + // Advection & rest of code // ------------------------------------------------------------------------------ From 2e095b19179dae11bbc374da62da302b7a0370e0 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 7 Sep 2026 10:15:03 +0100 Subject: [PATCH 10/10] Ability to send plasma data to device, update particle background send_plasma_data creates a vector of Hermes-3 plasma fields which gets sent to the SYCL target. update_recombination_markers and update_particle_background read this on-device and update the particle background. Markers can be done per VANTAGE call, but particles must be updated per timestep as they move around. For now, the values for the plasma are still taken from scalars. --- include/vantage.hxx | 2 + src/vantage.cxx | 84 +++++++++++++++++-- .../dmplex-vertex-coordinates/data/BOUT.inp | 1 + tests/integrated/particle-pusher/runtest | 1 + .../vantage-iz-rec-balance/data/BOUT.inp | 1 + 5 files changed, 80 insertions(+), 9 deletions(-) diff --git a/include/vantage.hxx b/include/vantage.hxx index 74427e04b..57963ab60 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -134,6 +134,8 @@ private: 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 diff --git a/src/vantage.cxx b/src/vantage.cxx index 45c022f2a..8787d45df 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -257,14 +257,43 @@ 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, - BoutReal background_ion_density, BoutReal background_ion_temperature, const std::vector& V_background) { // number of velocity dimensions @@ -275,17 +304,18 @@ void update_recombination_markers( // 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; + [=](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::write(Sym("FLUID_DENSITY")), + Access::read(plasma_data), Access::write(Sym("FLUID_DENSITY")), Access::write(Sym("FLUID_TEMPERATURE")), Access::write(Sym("ELECTRON_DENSITY")), Access::write(Sym("ELECTRON_TEMPERATURE")), @@ -325,6 +355,25 @@ void update_recombination_markers( } } +// 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) { Coordinates* coord = bout_mesh->getCoordinates(); @@ -694,6 +743,12 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) auto domain = std::make_shared(neso_mesh, mapper); // Get the number of cells in the mesh owned on this process 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)) { @@ -1081,9 +1136,14 @@ int VantageMonitor::call(Solver* UNUSED(solver), BoutReal time, int iter, // 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, neso_mesh, N_w, background_ion_density, - background_ion_temperature, V_background); + update_recombination_markers(marker_group, plasma_data, neso_mesh, N_w, V_background); // Advection & rest of code // ------------------------------------------------------------------------------ @@ -1148,10 +1208,16 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { for (int stepx = 0; stepx < nsteps; stepx++) { 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); diff --git a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp index fdb64693f..474d91051 100644 --- a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp +++ b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp @@ -37,6 +37,7 @@ charge = -1 dt = 0.005 nsteps = 3 +ion_species = d+ neutral_species = d test_mass_conservation = true diff --git a/tests/integrated/particle-pusher/runtest b/tests/integrated/particle-pusher/runtest index 0b0063bc9..6129460ab 100755 --- a/tests/integrated/particle-pusher/runtest +++ b/tests/integrated/particle-pusher/runtest @@ -84,6 +84,7 @@ def particle_push_input( 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 1b706ac1e..bf1a25c63 100644 --- a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp +++ b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp @@ -56,6 +56,7 @@ charge = -1 dt = 700 nsteps = 3 +ion_species = d+ neutral_species = d # This is the final expected steady state density