diff --git a/CMakeLists.txt b/CMakeLists.txt index 8a3eab18b..8f3d53b7c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -189,15 +189,24 @@ set(HERMES_SOURCES include/zero_current.hxx include/transform.hxx include/fixed_fraction_radiation.hxx - include/simple_pump.hxx -) + include/simple_pump.hxx) if(HERMES_USE_VANTAGE) # include Hermes-VANTAGE include and source files - list(APPEND HERMES_SOURCES - src/vantage_dmplex.cxx - include/vantage_dmplex.hxx - src/vantage.cxx + list( + APPEND + HERMES_SOURCES + src/vantage_dmplex.cxx + include/vantage_dmplex.hxx + src/vantage_helperfunctions.cxx + include/vantage_helperfunctions.hxx + src/vantage_datatransfer.cxx + include/vantage_datatransfer.hxx + src/vantage_sources.cxx + include/vantage_sources.hxx + src/vantage_diagnostics.cxx + include/vantage_diagnostics.hxx + src/vantage.cxx include/vantage.hxx) endif() @@ -406,11 +415,12 @@ if(HERMES_TESTS) option(HERMES_UNIT_TESTS "Build the unit tests" ON) option(SPACK_GTEST "Use spack-installed googletest for unit tests" ON) - # "gtest" and "GTest::gtest" are not the same, and the latter is preferred as being - # more complete. So, instead of just linking to "gtest", we will try to find the package - # using CMake if SPACK_GTEST is enabled. If not, then if gtest isn't already a target, then - # it will get it from the submodule. The GTest target is saved as its own variable - # so that it can be passed to both hermes_unit_tests. + # "gtest" and "GTest::gtest" are not the same, and the latter is preferred as + # being more complete. So, instead of just linking to "gtest", we will try to + # find the package using CMake if SPACK_GTEST is enabled. If not, then if + # gtest isn't already a target, then it will get it from the submodule. The + # GTest target is saved as its own variable so that it can be passed to both + # hermes_unit_tests. set(HERMES_GTEST_TARGET "") if(HERMES_UNIT_TESTS) if(SPACK_GTEST) diff --git a/include/vantage.hxx b/include/vantage.hxx index dc9f08751..1b91c7e90 100644 --- a/include/vantage.hxx +++ b/include/vantage.hxx @@ -2,65 +2,20 @@ #include "../include/component.hxx" #include "bout/bout.hxx" #include "bout/petsclib.hxx" +#include +#include +#include #include #include #include +#include +#include "../include/vantage_diagnostics.hxx" +#include "../include/vantage_sources.hxx" +#include "vantage_datatransfer.hxx" using namespace NESO::Particles; using namespace VANTAGE::Reactions; -/// @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 -/// ionisation). -/// @param accumulator CellwiseAccumulator to use to accumulate the source term for this -/// reaction. -/// @param particle_group ParticleGroup to which this source applies. -/// @param zeroer TransformationStrategy to use to zero the source term dat after -/// accumulation. -struct VantageSource { - std::string hermes_source_name; - std::string vantage_source_name; - std::shared_ptr> accumulator; - std::shared_ptr particle_group; - std::shared_ptr zeroer; - Field2D source_data; -}; - -/// @brief Class to manage reaction channel sources from VANTAGE. -/// Source terms from VANTAGE are extracted from the accumulator. -/// These are then converted to actual sources, e.g. units of m^-3 s^-1 for a -/// density source. -class VantageSourceManager { -public: - VantageSourceManager(std::shared_ptr& neso_mesh, - Mesh* bout_mesh, Options& units); - - Mesh* bout_mesh; - - // Register new source - void add_source(const std::string& hermes_source_name, - const std::string& vantage_source_name, - std::shared_ptr> accumulator, - std::shared_ptr particle_group, - std::shared_ptr zeroer); - - // Update the Hermes-3 source field using the accumulated data from corresponding - // VANTAGE source - void update_source(const std::string& hermes_source_name, double dt); - - // Call update_source on all sources - void update_all_sources(double dt); - - // Return data for a given Hermes-3 source name - Field2D get_data(const std::string& hermes_source_name); - -private: - std::map sources; - std::shared_ptr neso_mesh; - 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 @@ -98,33 +53,45 @@ private: BoutReal particle_time; BoutReal N_w; REAL dt; + BoutReal charge; + BoutReal AA; // mass + BoutReal initial_neutral_density; 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 + // Diagnostic variables on the kinetic mesh, for testing + std::vector neutral_density, total_density; + std::vector ion_density_kmsh; + // Field2D for storing plasma data coming from the plasma grid, that will be evaluated + // on to the kinetic mesh, and then on to the particles themselves. + Field2D electron_density, electron_temperature; + Field2D ion_density, ion_temperature; + std::vector ion_velocity; + // a threshold density, for reactions between neutrals and the plasma + REAL electron_density_threshold; BoutReal total_mass_initial, total_mass; std::string dmplex_filepath, vantage_dump_filepath, particle_data_filepath; // Path for output files + // volumes of neso_mesh cells (from neso_mesh->dmh->get_cell_volume()) + // in a vector of size of Nx*Ny, where Nx and Ny are the number of local + // x and y cells in the BOUT++ mesh (excluding guards) + std::vector neso_mesh_cell_volumes_on_plasma_grid; + 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::shared_ptr b2d; + std::shared_ptr project_eval_dg0; + std::shared_ptr mesh_coupler_dg0; + std::vector dof_kinetic_mesh_scalar; + std::vector dof_bout_mesh_scalar; + std::vector kinetic_mesh_map; // variable for recording the map from serial to parallelised DMPlex cells in terms of a vector of integers 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 @@ -132,9 +99,12 @@ private: 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 + std::shared_ptr data_transfer; // Manager for VANTAGE data transfer + std::unique_ptr + diagnostics_manager; // Manager for VANTAGE diagnostics + std::shared_ptr source_manager; // Manager for VANTAGE reaction sources + // These classes don't have a default constructor so need to be initialised as a unique_ptr VantageMonitor monitor{this}; // Output monitor to schedule VANTAGE iterations std::unique_ptr reaction_controller; diff --git a/include/vantage_datatransfer.hxx b/include/vantage_datatransfer.hxx new file mode 100644 index 000000000..40318e850 --- /dev/null +++ b/include/vantage_datatransfer.hxx @@ -0,0 +1,73 @@ +#pragma once +#include "bout/bout.hxx" +#include +#include +#include + +using namespace NESO::Particles; + +/// @brief Class to handle transfer of data between kinetic mesh and plasma grid in VANTAGE. +class VantageDataTransfer { +public: + VantageDataTransfer(std::shared_ptr& neso_mesh, + std::shared_ptr& project_eval_dg0, + std::shared_ptr& mesh_coupler, + Mesh* bout_mesh, size_t ndim_vector); + + // transfer functions for physics scalars + void transfer_scalar_to_plasma_grid( + std::vector& scalar_kinetic_mesh, + Field2D& scalar_plasma_grid); + void transfer_scalar_to_kinetic_mesh( + Field2D& scalar_plasma_grid, + std::vector& scalar_kinetic_mesh); + void transfer_scalar_to_particle_property( + std::vector& scalar_kinetic_mesh, + std::shared_ptr& A_particle_group, + std::string particle_property); + void transfer_scalar_to_particle_property( + Field2D& scalar_plasma_grid, + std::shared_ptr& A_particle_group, + std::string particle_property); + void transfer_particle_property_to_scalar( + std::shared_ptr& A_particle_group, + std::string particle_property, + std::vector& scalar_kinetic_mesh); + void transfer_particle_property_to_vector( + std::shared_ptr& A_particle_group, + std::string particle_property, + std::vector& vector_kinetic_mesh); + + // transfer functions for physics vectors + void transfer_vector_to_particle_property( + std::vector& vector_plasma_grid, + std::shared_ptr& A_particle_group, + std::string particle_property); + void transfer_vector_to_particle_property( + std::vector& vector_kinetic_mesh, + std::shared_ptr& A_particle_group, + std::string particle_property); + void transfer_vector_to_kinetic_mesh( + std::vector& vector_plasma_grid, + std::vector& vector_kinetic_mesh); + + +private: + // internal variables needed for data transfer + std::shared_ptr neso_mesh; + std::shared_ptr project_eval_dg0; + std::shared_ptr mesh_coupler; + + // the bout_mesh variable needed for Hermes-3/BOUT++ diagnostics + Mesh* bout_mesh; + + // variable used to transfer dofs from kinetic + // to BOUT++ meshes + size_t num_cells_owned_bout_grid; + std::vector dof_bout_grid_scalar; + std::vector dof_kinetic_mesh_scalar; + // number of physics vector components + size_t ndim_vector; + std::vector dof_bout_grid_vector; + std::vector dof_kinetic_mesh_vector; +}; diff --git a/include/vantage_diagnostics.hxx b/include/vantage_diagnostics.hxx new file mode 100644 index 000000000..55299c82e --- /dev/null +++ b/include/vantage_diagnostics.hxx @@ -0,0 +1,87 @@ +#pragma once +#include "bout/bout.hxx" +#include +#include +#include +#include +#include +#include "../include/vantage_datatransfer.hxx" +#include "../include/vantage_sources.hxx" + +using namespace NESO::Particles; + +REAL calculate_total_mass(std::vector& density, + std::shared_ptr& neso_mesh); +REAL calculate_total_mass(Field2D& density, + std::vector& neso_cell_volume_on_bout_mesh); + +/// @brief Class to manage diagnostics from VANTAGE. +class VantageDiagnosticsManager { +public: + VantageDiagnosticsManager(std::string vtkhdf_filename, + std::shared_ptr& neso_mesh, + // neso_mesh cell volumes on the BOUT++ mesh + std::vector& neso_cell_volumes, + std::shared_ptr& A_particle_group, + std::shared_ptr& data_transfer, + std::shared_ptr& source_manager, + BoutReal N_w, BoutReal mass, + Mesh* bout_mesh, Options& units, + std::string vantage_dump_filepath); + + // compute the kinetic velocity moments and + // store in private variables + void update_kinetic_velocity_moments(); + // write kinetic diagnostics to a vtkhdf file + void write_kinetic_velocity_moment_diagnostics(int istep, std::vector& ion_density); + // transfer kinetic moments to BOUT++ grid + void transfer_moments_to_plasma_grid(); + // write BOUT++ style diagnostics on the BOUT++ grid + void write_bout_diagnostics( + Field2D& ion_density, + BoutReal particle_time); + std::vector get_density_kinetic_mesh(); + // Field2D transfer_scalar_to_plasma_grid(std::vector& scalar_field); + +private: + // internal variables needed for diagnostics + std::string vtkhdf_filename; + std::shared_ptr neso_mesh; + std::shared_ptr A_particle_group; + std::shared_ptr data_transfer; + // pointer to vantage source_manager for diagnostics + std::shared_ptr source_manager; + BoutReal N_w; + BoutReal mass; + + // the bout_mesh variable needed for Hermes-3/BOUT++ diagnostics + Mesh* bout_mesh; + Options& units; + std::string vantage_dump_filepath; + 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 + + // variables used to store the moments of + // the neutral distribution function, on + // the kinetic mesh + const size_t ndimv = 2; // number of velocity dimensions + std::vector density; + std::vector energy; + std::vector gamma; + std::vector uvector; + std::vector pressure; + std::vector temperature; + + // variable used to transfer dofs from kinetic + // to BOUT++ meshes + std::vector dof_bout_mesh_scalar; + // variables used to store the moments of the + // neutral distribution function projected on + // to the BOUT++ grid + Field2D density_plasma_grid; + Field2D energy_plasma_grid; + Field2D pressure_plasma_grid; + Field2D temperature_plasma_grid; + // n.b. only treat scalar variables for now +}; diff --git a/include/vantage_dmplex.hxx b/include/vantage_dmplex.hxx index 88584ebd1..e64830689 100644 --- a/include/vantage_dmplex.hxx +++ b/include/vantage_dmplex.hxx @@ -50,10 +50,21 @@ std::vector cells_definition_from_RZ_ivertex( Field2D& ivertex_lower_right_corners, Field2D& ivertex_upper_right_corners, Field2D& ivertex_upper_left_corners); -DM create_dmplex_from_Bout_mesh(Mesh* bout_mesh, Options& mesh_options, - std::shared_ptr sycl_target, - std::string dmplex_h5_filename); +void create_dmplex_from_Bout_mesh(DM* dm, Mesh* bout_mesh, Options& mesh_options, + std::shared_ptr sycl_target); -#endif +void create_dmplex_from_GMSH_msh(DM* dm, std::string msh_file); + +void write_dmplex_to_file(DM dm, std::string dmplex_name, std::string dmplex_h5_filename); + +BoutReal get_triangle_area(size_t itriangle, + const std::vector& vertices, + const std::vector& tri_cell_vertices); + +std::vector get_triangle_vertices(); + +std::vector get_triangle_cell_definition(); + +#endif #endif // VANTAGE_DMPLEX_H diff --git a/include/vantage_helperfunctions.hxx b/include/vantage_helperfunctions.hxx new file mode 100644 index 000000000..9c8cf924f --- /dev/null +++ b/include/vantage_helperfunctions.hxx @@ -0,0 +1,33 @@ +#pragma once +#include "bout/bout.hxx" +#include + +using namespace NESO::Particles; + +// get the number of cells owned by the local +// BOUT++ mesh, excluding guard cells +size_t get_num_cells_owned_bout_mesh(Mesh*& bout_mesh); + +// get the cell volumes from the NESO-Particles mesh +// onto the same degrees of freedom owned by the local BOUT++ mesh +std::vector get_cell_volumes_on_plasma_grid( + DM& dm, std::vector& kinetic_mesh_map, + std::shared_ptr& neso_mesh, Mesh*& bout_mesh); + +// functions used for checks of the DMPlex + +// Check that the x, y cell volumes of BOUT++ match the +// inferred quad cell volumes from the NESO-Particles mesh +void check_cell_volumes(std::vector neso_cell_volumes_bmsh, Mesh*& bout_mesh, + Options& alloptions); + +// Check that the Rxy, Zxy cell centres of BOUT++ match the +// inferred quad cell centres from the NESO-Particles mesh +void check_cell_centres(Options& alloptions, DM& dm, + std::vector& kinetic_mesh_map, + std::shared_ptr& neso_mesh, + Mesh*& bout_mesh, BoutReal absolute_tolerance, + BoutReal relative_tolerance); + +// Function to check mass conservation at the end of particle pushing +void check_mass_conservation(REAL total_mass_final, REAL total_mass_initial); diff --git a/include/vantage_sources.hxx b/include/vantage_sources.hxx new file mode 100644 index 000000000..8b2552903 --- /dev/null +++ b/include/vantage_sources.hxx @@ -0,0 +1,95 @@ +#pragma once +#include "bout/bout.hxx" +#include +#include +#include +#include +#include +#include +#include +#include +#include "../include/vantage_datatransfer.hxx" + +using namespace NESO::Particles; +using namespace VANTAGE::Reactions; + +/// @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 +/// ionisation). +/// @param accumulator CellwiseAccumulator to use to accumulate the source term for this +/// reaction. +/// @param particle_group ParticleGroup to which this source applies. +/// @param zeroer TransformationStrategy to use to zero the source term dat after +/// accumulation. +struct VantageSource { + std::string hermes_source_name; + std::string vantage_source_name; + std::string bout_diagnostic_units; + BoutReal bout_diagnostic_conversion; + std::string bout_diagnostic_standard_name; + std::string bout_diagnostic_long_name; + std::shared_ptr> accumulator; + std::shared_ptr particle_group; + std::shared_ptr zeroer; + Field2D source_data_plasma_grid; + std::vector source_data_kinetic_mesh; +}; + +/// @brief Class to manage reaction channel sources from VANTAGE. +/// Source terms from VANTAGE are extracted from the accumulator. +/// These are then converted to actual sources, e.g. units of m^-3 s^-1 for a +/// density source. +class VantageSourceManager { +public: + VantageSourceManager(std::shared_ptr& neso_mesh, + std::shared_ptr& data_transfer, + Mesh* bout_mesh, Options& units); + + Mesh* bout_mesh; + + // Register new source + void add_source(const std::string& hermes_source_name, + const std::string& vantage_source_name, + const std::string& bout_diagnostic_units, + BoutReal bout_diagnostic_conversion, + const std::string& bout_diagnostic_standard_name, + const std::string& bout_diagnostic_long_name, + std::shared_ptr> accumulator, + std::shared_ptr particle_group, + std::shared_ptr zeroer); + + // Update the Hermes-3 source field using the accumulated data from corresponding + // VANTAGE source + void update_source(const std::string& hermes_source_name, double dt); + + // Call update_source on all sources + void update_all_sources(double dt); + + // Get vector of source names + std::vector get_source_names(); + + // Return diagnostic units data + std::string get_units(const std::string& hermes_source_name); + + // Return diagnostic units data + BoutReal get_conversion(const std::string& hermes_source_name); + + // Return diagnostic units data + std::string get_standard_name(const std::string& hermes_source_name); + + // Return diagnostic units data + std::string get_long_name(const std::string& hermes_source_name); + + // Return data for a given Hermes-3 source name on the plasma grid + Field2D get_plasma_grid_data(const std::string& hermes_source_name); + + // Return data for a given Hermes-3 source name on the kinetic mesh + std::vector get_kinetic_mesh_data(const std::string& hermes_source_name); + +private: + std::map sources; + std::shared_ptr neso_mesh; + std::shared_ptr data_transfer; + Options& units; +}; \ No newline at end of file diff --git a/scripts/dmplex_tools_for_particle_pusher/particle_animator.py b/scripts/dmplex_tools_for_particle_pusher/particle_animator.py index 44ccc1633..22dd2c869 100644 --- a/scripts/dmplex_tools_for_particle_pusher/particle_animator.py +++ b/scripts/dmplex_tools_for_particle_pusher/particle_animator.py @@ -5,6 +5,7 @@ from matplotlib.animation import FuncAnimation from petsc4py import PETSc import argparse +import xhermes parser = argparse.ArgumentParser( description="Animate particles moving on a DMPlex mesh." @@ -19,6 +20,21 @@ type=str, help="The path to the HDF5 file representing the particle data", ) +parser.add_argument( + "BOUT_file_path", + type=str, + help="The path to the BOUT.dmp.0.nc file associated with the particle data", +) +parser.add_argument( + "--equal-aspect", + action="store_true", + help="Use equal aspect ratio R, Z axes", +) +parser.add_argument( + "--set-xlim-zero", + action="store_true", + help="Use 0 as the minimum R on the axes", +) args = parser.parse_args() print( @@ -26,6 +42,10 @@ ) +def get_length_normalisation(BOUT_file_path): + ds = xhermes.open(BOUT_file_path) + return ds.attrs["metadata"]["rho_s0"] + def load_dmplex(file_path): dm = PETSc.DMPlex().create() viewer = PETSc.Viewer().createHDF5(file_path, "r") @@ -53,6 +73,8 @@ def get_mesh_edges(dm): edges.append((x0, x1)) return edges +# normalisation for particle data +meters = get_length_normalisation(args.BOUT_file_path) dm = load_dmplex(args.dmplex_h5_file_path) # dm = load_dmplex('dmplex/expected_nonorthogonal.grd.nc.mesh.h5') @@ -78,7 +100,7 @@ def load_particle_data(file_path): pdata = np.zeros((nparticles, 2)) pdata[:, 0] = P_0 pdata[:, 1] = P_1 - particle_positions.append(pdata) + particle_positions.append(np.multiply(pdata,meters)) except KeyError as error: print(f"No particles at time step {it}: {error}") # assign empty particle data @@ -107,6 +129,10 @@ def update_plot(i, data, scat): ax.set_title("Particle Positions") ax.set_xlabel("R") ax.set_ylabel("Z") +if args.equal_aspect: + ax.set_aspect("equal",adjustable="box") +if args.set_xlim_zero: + ax.set_xlim(0.0,None) def update(frame): @@ -118,6 +144,6 @@ def update(frame): ani = FuncAnimation(fig, update, frames=nstep, interval=50, blit=True) output_path = args.particle_trajectory_h5_file_path + ".animation.gif" -ani.save(output_path) +ani.save(output_path, dpi=400) print(f"Saving animation of particle paths to {output_path}") # plt.show() diff --git a/src/vantage.cxx b/src/vantage.cxx index 47210f450..8fbd81776 100644 --- a/src/vantage.cxx +++ b/src/vantage.cxx @@ -3,16 +3,21 @@ #include "bout/field2d.hxx" #include "bout/output.hxx" #include "bout/petsclib.hxx" +#include +#include #include #include #include +#include #include #include +#include #include #include #include #include #include +#include #include #include #include @@ -23,8 +28,12 @@ #include // for reactions integration #include "../include/amjuel_data.hxx" +#include "../include/component.hxx" #include "../include/vantage.hxx" +#include "../include/vantage_datatransfer.hxx" +#include "../include/vantage_diagnostics.hxx" #include "../include/vantage_dmplex.hxx" +#include "../include/vantage_helperfunctions.hxx" #include #ifndef NESO_PARTICLES_PETSC @@ -62,388 +71,54 @@ std::string make_output_path(const std::string& filename, Options& alloptions) { return fmt::format("{}/{}", output_dir, filename); } -void calculate_neutral_density_in_place( - Field2D& density, std::shared_ptr& dg0, - std::shared_ptr& A_particle_group, std::vector& h_project1, - BoutReal N_w) { - Mesh* bout_mesh = density.getMesh(); - // get a density by projecting the particle property WEIGHT to the bout_mesh - dg0->project(A_particle_group, Sym("WEIGHT")); - // std::vector h_project1; - dg0->get_dofs(1, h_project1); - std::size_t ic = 0; - for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { - for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { - density(ix, iy) = h_project1.at(ic) * N_w; - ic++; - } - } - // this fills internal guards - bout_mesh->communicate(density); - // apply boundary conditions to fill external guards - // density.applyBoundary(); - // extrapolate -> Neumann -} - -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++) { - for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { - local_mass += density(ix, iy) * neso_mesh->dmh->get_cell_volume(ic); - ic++; - } - } - MPICHK( - MPI_Allreduce(&local_mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, BoutComm::get())); - return total_mass; -} - -Options -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 rho_s0 = get(alloptions["units"]["meters"]); - auto Bnorm = get(alloptions["units"]["Tesla"]); - auto Cs0 = get(alloptions["units"]["meters"]) - / get(alloptions["units"]["seconds"]); - - Options bout_output_data; - 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["Nn"], neutral_density, - {{"time_dimension", "t"}, - {"units", "m^-3"}, - {"conversion", Nnorm}, - {"standard_name", "Density"}, - {"long_name", "Kinetic neutral density"}, - {"species", "kinetic neutrals"}, - {"source", "vantage"}}); - - set_with_attrs(bout_output_data["Siz"], Field2D{0.0, bout_mesh}, - {{"time_dimension", "t"}, - {"units", "m^-3 s^-1"}, - {"conversion", Nnorm * Omega_ci}, - {"standard_name", "Density source"}, - {"long_name", "Ionisation density source"}, - {"species", "kinetic neutrals"}, - {"source", "vantage"}}); - - set_with_attrs(bout_output_data["Srec"], Field2D{0.0, bout_mesh}, - {{"time_dimension", "t"}, - {"units", "m^-3 s^-1"}, - {"conversion", Nnorm * Omega_ci}, - {"standard_name", "Density source"}, - {"long_name", "Recombination density source"}, - {"species", "kinetic neutrals"}, - {"source", "vantage"}}); - - // Integrals - Field2D total_density = ion_density + neutral_density; - set_with_attrs(bout_output_data["total_mass"], - calculate_total_mass(total_density, neso_mesh), - {{"time_dimension", "t"}}); - - set_with_attrs(bout_output_data["total_neutral_mass"], - calculate_total_mass(neutral_density, neso_mesh), - {{"time_dimension", "t"}}); - - set_with_attrs(bout_output_data["total_ion_mass"], - calculate_total_mass(ion_density, neso_mesh), {{"time_dimension", "t"}}); - - set_with_attrs(bout_output_data["t_array"], 0.0, {{"time_dimension", "t"}}); - - // Add metadata from mesh, e.g. branch cuts - 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(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, bout::OptionsIO& vantage_dump_writer, - BoutReal particle_time) { - // update density in Options object and write - bout_output_data["neutral_density"] = neutral_density; - bout_output_data["Nn"] = neutral_density; - bout_output_data["Siz"] = Siz; - bout_output_data["Srec"] = Srec; - bout_output_data["ion_density"] = ion_density; - Field2D total_density = ion_density + neutral_density; - bout_output_data["total_mass"] = calculate_total_mass(total_density, neso_mesh); - bout_output_data["total_neutral_mass"] = - calculate_total_mass(neutral_density, neso_mesh); - bout_output_data["total_ion_mass"] = calculate_total_mass(ion_density, neso_mesh); - bout_output_data["t_array"] = particle_time; - // bout_output_data["t_array"] = 0.0; - // Append data to file - 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( - Field2D& initial_neutral_density, - std::shared_ptr& dg0, - std::shared_ptr& A_particle_group, + BoutReal& initial_neutral_density, std::shared_ptr& A_particle_group, std::shared_ptr& neso_mesh, - 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++) { - for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { - // 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. - 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; - h_project1.at(static_cast(ixy)) = particle_weights; - ixy++; - } + std::vector& dof_kinetic_mesh_scalar, + std::shared_ptr& data_transfer, BoutReal N_w) { + // set a constant density across the entire kinetic mesh + const size_t ncell = dof_kinetic_mesh_scalar.size(); + for (size_t ic = 0; ic < ncell; ic++) { + // 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. + const REAL cell_volume = neso_mesh->dmh->get_cell_volume(static_cast(ic)); + const INT nmarkers_per_cell = A_particle_group->get_npart_cell(static_cast(ic)); + const REAL particle_weights = initial_neutral_density * cell_volume + / static_cast(nmarkers_per_cell) / N_w; + dof_kinetic_mesh_scalar.at(ic) = particle_weights; } // now copy the data to internal variables - dg0->set_dofs(1, h_project1); - // set the data from internal variables into the weights - dg0->evaluate(A_particle_group, Sym("WEIGHT")); -} - -void check_cell_volumes(std::shared_ptr& neso_mesh, - Mesh*& bout_mesh, Options& alloptions) { - Coordinates* coord = bout_mesh->getCoordinates(); - PetscInt ixy = 0; - const REAL tolerance = 1.0e-12; - for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { - for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { - - const BoutReal meters = get(alloptions["units"]["meters"]); - const BoutReal meters_squared = meters * meters; - const BoutReal meters_cubed = meters * meters * meters; - - // Convert to SI: dx is m^2 T, J is m/T, dy is unitless, skip dz - // so J * dx * dy = m^3, technically per radian toroidal angle due to missing dz - const BoutReal bout_cell_area = - coord->J(ix, iy) * coord->dx(ix, iy) * coord->dy(ix, iy) * meters_cubed; - - // Straight up 2D grid, needs m^2 - const REAL neso_cell_area = neso_mesh->dmh->get_cell_volume(ixy) * meters_squared; - - const bool volumes_match = (abs(bout_cell_area - neso_cell_area) < tolerance); - - // exit if we fail to find a match - NESOASSERT(volumes_match, - fmt::format("BOUT++ mesh volume {} does not match NESO-Particles mesh " - "volume {} for ix = {} iy = {} \n Ignore this message by " - "setting [neso_particles] test_cell_volumes = false", - bout_cell_area, neso_cell_area, ix, iy)); - ixy++; - } - } + data_transfer->transfer_scalar_to_particle_property(dof_kinetic_mesh_scalar, + A_particle_group, "WEIGHT"); } -REAL cell_length(std::vector>& cell_vertices, std::size_t iv1, - std::size_t iv2, std::size_t iv3, std::size_t iv4) { - const REAL Rlength2 = - std::pow(0.5 - * (cell_vertices.at(iv1).at(0) + cell_vertices.at(iv2).at(0) - - cell_vertices.at(iv3).at(0) - cell_vertices.at(iv4).at(0)), - 2.0); - const REAL Zlength2 = - std::pow(0.5 - * (cell_vertices.at(iv1).at(1) + cell_vertices.at(iv2).at(1) - - cell_vertices.at(iv3).at(1) - cell_vertices.at(iv4).at(1)), - 2.0); - REAL length = std::pow(Zlength2 + Rlength2, 0.5); - return length; +void update_particle_properties_from_plasma( + std::shared_ptr& data_transfer, + std::shared_ptr& A_particle_group, Field2D& ion_density, + Field2D& ion_temperature, std::vector& ion_velocity, + Field2D& electron_density, Field2D& electron_temperature) { + data_transfer->transfer_scalar_to_particle_property(ion_density, A_particle_group, + "FLUID_DENSITY"); + data_transfer->transfer_scalar_to_particle_property(ion_temperature, A_particle_group, + "FLUID_TEMPERATURE"); + data_transfer->transfer_vector_to_particle_property(ion_velocity, A_particle_group, + "FLUID_FLOW_SPEED"); + data_transfer->transfer_scalar_to_particle_property(electron_density, A_particle_group, + "ELECTRON_DENSITY"); + data_transfer->transfer_scalar_to_particle_property( + electron_temperature, A_particle_group, "ELECTRON_TEMPERATURE"); } -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 - Field2D Rxy; - Field2D Zxy; - bout_mesh->get(Rxy, "Rxy"); - bout_mesh->get(Zxy, "Zxy"); - - BoutReal meters = get(alloptions["units"]["meters"]); - - // compare to cell centres calculated from cell corners - std::vector> cell_vertices; - PetscInt ixy = 0; - for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { - for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { - const REAL bout_Rxy = Rxy(ix, iy); - const REAL bout_Zxy = Zxy(ix, iy); - neso_mesh->dmh->get_cell_vertices(ixy, cell_vertices); - - REAL neso_Rxy = 0.0; - REAL neso_Zxy = 0.0; - for (std::size_t iv = 0; iv < 4; iv++) { - // DMPlex is stored in normalised units, need conversion to [m] - neso_Rxy += cell_vertices.at(iv).at(0) * meters; - neso_Zxy += cell_vertices.at(iv).at(1) * meters; - } - neso_Rxy /= 4.0; - neso_Zxy /= 4.0; - // get lengths of cell across the two dimensions - const REAL cell_length_a = cell_length(cell_vertices, 0, 1, 2, 3) * meters; - const REAL cell_length_b = cell_length(cell_vertices, 0, 3, 2, 1) * meters; - const REAL min_cell_length = std::min(cell_length_a, cell_length_b); - // we compare the difference in cell centres to the absolute tolerance and - // the relative tolerance formed by comparing to the smallest length across the cell - const REAL tolerance = absolute_tolerance + min_cell_length * relative_tolerance; - const bool centres_match = (abs(neso_Rxy - bout_Rxy) < tolerance) - && (abs(neso_Zxy - bout_Zxy) < tolerance); - // exit if we fail to find a match - NESOASSERT( - centres_match, - fmt::format("Hypnotoad/BOUT++ cell centre (R, Z) ({}, {}) does not match " - "NESO-Particles mesh " - "cell centre ({}, {}) for ix = {} iy = {} \n" - "The cell height and width are {} {} \n" - "The displacements in R and Z are {} {} \n" - "Ignore this message by " - "setting [neso_particles] test_cell_centres = false\n Relax the " - "tolerance used in this check by increasing\n" - "[neso_particles] cell_centre_absolute_tolerance = {}\n" - "[neso_particles] cell_centre_relative_tolerance = {}", - bout_Rxy, bout_Zxy, neso_Rxy, neso_Zxy, ix, iy, cell_length_a, - cell_length_b, abs(neso_Rxy - bout_Rxy), abs(neso_Zxy - bout_Zxy), - absolute_tolerance, relative_tolerance)); - ixy++; - } - } -} - -void check_mass_conservation(BoutReal total_mass_final, BoutReal total_mass_initial) { - BoutReal rtol = 1.0e-13; - BoutReal mass_conserved = - (abs(total_mass_final - total_mass_initial) < rtol * total_mass_initial); - // exit if we fail to find conservation - NESOASSERT(mass_conserved, - fmt::format("Initial total mass {} does not match " - "final total mass {} \n Ignore this message by " - "setting [neso_particles] test_mass_conservation = false", - total_mass_initial, total_mass_final)); -} - -// VANTAGE source manager implementation -// ------------------------------------------------------------------------------ -VantageSourceManager::VantageSourceManager( - std::shared_ptr& neso_mesh, Mesh* bout_mesh, - Options& units) - : bout_mesh(bout_mesh), neso_mesh(neso_mesh), units(units) {} - -// Register new source with the manager and initialise its data -void VantageSourceManager::add_source( - const std::string& hermes_source_name, const std::string& vantage_source_name, - std::shared_ptr> accumulator, - std::shared_ptr particle_group, - std::shared_ptr zeroer) { - - Field2D source_data{bout_mesh}; - source_data = 0.0; - - VantageSource source{ - hermes_source_name, vantage_source_name, accumulator, particle_group, zeroer, - source_data}; - - this->sources[hermes_source_name] = source; -} - -// Return source data -Field2D VantageSourceManager::get_data(const std::string& hermes_source_name) { - return this->sources[hermes_source_name].source_data; -} - -// Update the source from VANTAGE and reset the VANTAGE data/accumulator -void VantageSourceManager::update_source(const std::string& hermes_source_name, - double dt) { - - VantageSource& source = this->sources[hermes_source_name]; - BoutReal N_w = get(units["N_w"]); - - std::vector> accumulated_1d = - source.accumulator->get_cell_data(source.vantage_source_name); - - std::size_t ic = 0; - for (int ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { - for (int iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { - source.source_data(ix, iy) = - accumulated_1d[ic]->at(0, 0) // Total weight - * N_w // Total particles - / neso_mesh->dmh->get_cell_volume(static_cast(ic)) // Total density - / dt; // Density source - - ic++; - } - } - - // Fill internal guards - bout_mesh->communicate(source.source_data); - // Reset the accumulator object - source.accumulator->zero_buffer(source.vantage_source_name); - // Reset the accumulated source data on the particle - source.zeroer->transform(std::make_shared(source.particle_group)); -} - -// Update all sources -void VantageSourceManager::update_all_sources(double dt) { - for (auto& [hermes_source_name, source] : this->sources) { - update_source(hermes_source_name, dt); - } +// Create a ParticleSubGroup from particles that are in a cell with nonzero electron_density. +ParticleSubGroupSharedPtr create_particle_sub_group_in_plasma_volume( + std::shared_ptr& A_particle_group, + const REAL electron_density_threshold) { + ParticleSubGroupSharedPtr particle_group_in_plasma = particle_sub_group( + A_particle_group, [=](auto ne) { return (ne[0] > electron_density_threshold); }, + Access::read(Sym("ELECTRON_DENSITY"))); + return particle_group_in_plasma; } Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) @@ -470,6 +145,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) Options& units = alloptions["units"]; BoutReal inv_meters_cubed = get(units["inv_meters_cubed"]); BoutReal eV = get(units["eV"]); + BoutReal pascal = SI::qe * eV * inv_meters_cubed; BoutReal meters = get(units["meters"]); BoutReal seconds = get(units["seconds"]); @@ -478,6 +154,12 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) "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); + electron_density_threshold = + options["electron_density_reaction_threshold"] + .doc("Parameter controlling the minimum (normalised) electron density " + "at which the reactions between neutrals and charged plasma species are " + "applied.") + .withDefault(1.0e-12); Options::root()["units"]["N_w"] = N_w; Options::root()["units"]["N_w"].setConditionallyUsed(); @@ -487,24 +169,47 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) // 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_name = mesh_options["dmplex_name"] + .doc("DMPlex object name.") + .withDefault("hypnotoad_dmplex_mesh"); 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 + bool use_external_msh = + mesh_options["use_external_msh"] + .doc("Use an externally generated .msh file for the kinetic mesh. " + "Not default and recommendation is false.") + .withDefault(false); + // Create and save DMPlex + // DM dm; // pointer to DMPlex, initialised below + // This DM is created in SI units without boundary labels + if (use_external_msh) { + std::string msh_file = + mesh_options["msh_file"] + .doc("Path to an externally generated .msh file for the kinetic mesh. ") + .withDefault("kinetic.msh"); + // create a DMPlex in serial + create_dmplex_from_GMSH_msh(&dm, msh_file); + PetscSF + sf_kinetic_mesh; // Petsc variable that records map of vertices from original vector to distributed vector indices + PetscInterface::generic_distribute(&dm, BoutComm::get(), 1, &sf_kinetic_mesh); + kinetic_mesh_map = + PetscInterface::get_global_distributed_points_map(dm, sf_kinetic_mesh); + } else { + create_dmplex_from_Bout_mesh(&dm, bout_mesh, mesh_options, sycl_target); + } + // label DMPlex boundaries + PetscInterface::label_all_dmplex_boundaries(dm, PetscInterface::face_sets_label, 100); + // diagnose the DMPlex by writing to file dmplex_filepath = make_output_path(dmplex_h5_filename, alloptions); + write_dmplex_to_file(dm, dmplex_name, dmplex_filepath); + // Create paths for other diagnostics 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, dmplex_filepath); - - // Normalise DMPlex after creation + // Normalise DMPlex after creation to go from SI to normalised units // Get local coords object (i.e. per rank) and scale it - this scales entire mesh // All following interactions with the DMPlex will be in normalised units. Vec coords = nullptr; @@ -535,12 +240,32 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) .withDefault(true); // Normalisations + // Charge for ionised species in IZ reaction and mass of ion and neutral + charge = options["charge"].doc("Particle charge. electrons = -1").withDefault(1.0); + AA = options["AA"].doc("Particle atomic mass. Proton = 1").withDefault(1.0); + // check mass positive + ASSERT1(AA > 0.0); // Initial neutral parameters - initial_neutral_density = - options["initial_neutral_density"] - .doc("Initial neutral density for VANTAGE kinetic neutrals [m^-3]") - .as() - / inv_meters_cubed; + const BoutReal initial_neutral_pressure = + options["initial_neutral_pressure"] + .doc( + "Initial neutral pressure for VANTAGE kinetic neutrals [Pa], default = 1") + .withDefault(1.0) + / pascal; + const BoutReal initial_neutral_temperature = + options["initial_neutral_temperature"] + .doc("Initial neutral temperature for VANTAGE kinetic neutrals [eV], default " + "= 1") + .withDefault(1.0) + / eV; + // check initial neutral pressure is greater than or equal to zero + ASSERT1(initial_neutral_pressure >= 0.0); + // checking initial temperature greater than zero before division + ASSERT1(initial_neutral_temperature > 0.0); + initial_neutral_density = initial_neutral_pressure / initial_neutral_temperature; + // standard deviation (thermal speed) from initial condition + const BoutReal initial_neutral_thermal_speed = + std::sqrt(initial_neutral_temperature / AA); const int npart_per_cell = options["npart_per_cell"] .doc("Number of VANTAGE kinetic neutral particles per " "cell during initialisation") @@ -589,8 +314,6 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) .doc("Number of RNG samples to prepare per-particle") .withDefault(40); - ion_density = Field2D(background_ion_density, bout_mesh); - neutral_density = Field2D(0.0, bout_mesh); // Create a mesh interface from the DM neso_mesh = std::make_shared(dm, 0, BoutComm::get()); // Create a mapper for mapping particles into cells. @@ -598,16 +321,17 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) std::make_shared(sycl_target, neso_mesh); // 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 - num_cells_owned = neso_mesh->get_cell_count(); + // get the cell volumes from neso_mesh on the plasma grid, in the compound index + neso_mesh_cell_volumes_on_plasma_grid = + get_cell_volumes_on_plasma_grid(dm, kinetic_mesh_map, neso_mesh, bout_mesh); // 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)) { - check_cell_volumes(neso_mesh, bout_mesh, alloptions); + check_cell_volumes(neso_mesh_cell_volumes_on_plasma_grid, bout_mesh, alloptions); } if (mesh_options["test_dmplex_cell_centres"].withDefault(true)) { check_cell_centres( - alloptions, neso_mesh, bout_mesh, + alloptions, dm, kinetic_mesh_map, 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)); } @@ -615,7 +339,7 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) // create a Reactions particle spec auto particle_spec_builder = ParticleSpecBuilder(ndim); auto electron_species = Species("ELECTRON"); - auto main_species = Species("ION", 1.0, 0.0, 0); + auto main_species = Species("ION", AA, charge, 0); std::vector fluid_species = {electron_species, main_species}; particle_spec_builder.add_particle_prop(Properties( fluid_species, @@ -626,11 +350,15 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) Properties(fluid_species, std::vector{default_properties.source_momentum}), ndim); - ParticleSpec additional_props{ParticleProp(Sym("TSP"), 2), - ParticleProp(Sym("FLUID_DENSITY"), 1), - ParticleProp(Sym("FLUID_FLOW_SPEED"), ndim), - ParticleProp(Sym("FLUID_TEMPERATURE"), 1), - ParticleProp(Sym("N_CELL"), 1)}; + ParticleSpec additional_props{ + ParticleProp(Sym("TSP"), 2), + ParticleProp(Sym("FLUID_DENSITY"), 1), + ParticleProp(Sym("FLUID_FLOW_SPEED"), ndim), + ParticleProp(Sym("FLUID_TEMPERATURE"), 1), + ParticleProp(Sym("N_CELL"), 1), + ParticleProp(Sym("WEIGHT_V2"), 1), + ParticleProp(Sym("WEIGHT_V"), ndim), + }; particle_spec_builder.add_particle_spec(additional_props); ParticleSpec particle_spec = particle_spec_builder.get_particle_spec(); @@ -650,8 +378,9 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) &rng_pos); const int N_actual = static_cast(particle_cell_ids.size()); - auto velocities = - NESO::Particles::normal_distribution(N_actual, 2, 0.0, 1.0, rng_vel); + // use the 3D definition of sigma here, but note ndim = 2 for now + auto velocities = NESO::Particles::normal_distribution( + N_actual, 2, 0.0, initial_neutral_thermal_speed, rng_vel); int id_offset = 0; MPICHK(MPI_Exscan(&N_actual, &id_offset, 1, MPI_INT, MPI_SUM, @@ -685,13 +414,91 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) initial_distribution[Sym("CELL_ID")][px][0] = particle_cell_ids.at(pxu); initial_distribution[Sym("ID")][px][0] = px + id_offset; initial_distribution[Sym("WEIGHT")][px][0] = 1.0; + // these diagnostic properties are updated by write_kinetic_velocity_moment_diagnostics() + initial_distribution[Sym("WEIGHT_V2")][px][0] = 0.0; + for (int dimx = 0; dimx < ndim; dimx++) { + initial_distribution[Sym("WEIGHT_V")][px][dimx] = 0.0; + } } // Add the new particles to the particle group A_particle_group->add_particles_local(initial_distribution); + const size_t num_cells_owned_bout_mesh = get_num_cells_owned_bout_mesh(bout_mesh); + // Get the number of cells in the kinetic (neutral) mesh owned on this process + const int num_cells_owned_kinetic_mesh = neso_mesh->get_cell_count(); + // allocate buffer vector for scalar projection/evaluation of NESO-Particles + // properties on to the kinetic mesh + dof_kinetic_mesh_scalar = + std::vector(static_cast(num_cells_owned_kinetic_mesh)); + // allocate buffer vector for scalar projection/evaluation of NESO-Particles + // properties on to the bout mesh + dof_bout_mesh_scalar = + std::vector(static_cast(num_cells_owned_bout_mesh)); // make pointer to projection object - dg0 = std::make_shared(neso_mesh, - sycl_target, "DG", 0); - + if (use_external_msh) { + // create the dg0 variable using a constructor that + // respects the kinetic mesh external definition + std::vector> coupler_map( + static_cast(num_cells_owned_bout_mesh)); + Field2D map_RZ_to_itriangle_0; + Field2D map_RZ_to_itriangle_1; + bout_mesh->get(map_RZ_to_itriangle_0, "map_RZ_to_itriangle_0"); + bout_mesh->get(map_RZ_to_itriangle_1, "map_RZ_to_itriangle_1"); + // get data that defines triangular cells + const std::vector vertices = get_triangle_vertices(); + const std::vector tri_cell_vertices = get_triangle_cell_definition(); + int icell = 0; + for (int ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (int iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + // get triangle areas, and total area for ratio in the backward weights + const int itri_0 = static_cast(map_RZ_to_itriangle_0(ix, iy)); + const REAL area_0 = + get_triangle_area(static_cast(itri_0), vertices, tri_cell_vertices); + const int itri_1 = static_cast(map_RZ_to_itriangle_1(ix, iy)); + const REAL area_1 = + get_triangle_area(static_cast(itri_1), vertices, tri_cell_vertices); + const REAL total_area = area_0 + area_1; + // std::cout << "total area: " << total_area << " area_0: " << area_0 << " area_1: " << area_1 << " area_0/total_area: " << area_0/total_area << " area_1/total_area: " << area_1/total_area <<'\n'; + ASSERT1(total_area > 0.0); + // lower triangle + coupler_map.at(static_cast(icell)) + .push_back({kinetic_mesh_map.at(static_cast(itri_0)), 1.0, + area_0 / total_area}); + // upper triangle + coupler_map.at(static_cast(icell)) + .push_back({kinetic_mesh_map.at(static_cast(itri_1)), 1.0, + area_1 / total_area}); + icell += 1; + } + } + // object for transferring data between kinetic and bout mesh degree-of-freedom vectors + mesh_coupler_dg0 = + std::make_shared(dm, coupler_map); + } + // if (mesh_coupler_dg0 == nullptr){ + // output << "mesh_coupler_dg0 is a nullptr" << std::endl; + // } + // object for evaluating/projecting particle properties + // between the kinetic mesh degree-of-freedom vector and particles + project_eval_dg0 = std::make_shared( + neso_mesh, sycl_target, "DG", 0); + + // Object for transferring data between BOUT++ and NESO-Particles data formats + this->data_transfer = std::make_shared( + neso_mesh, project_eval_dg0, mesh_coupler_dg0, bout_mesh, ndim); + + // vectors for storing an ion density on the kinetic mesh + ion_density_kmsh = std::vector( + static_cast(num_cells_owned_kinetic_mesh), background_ion_density); + total_density = + std::vector(static_cast(num_cells_owned_kinetic_mesh), 0.0); + // Field2D for storing plasma data coming from the plasma grid + // that will be evaluated on to the particle properties + ion_density = Field2D{background_ion_density, bout_mesh}; + electron_density = Field2D{background_electron_density, bout_mesh}; + ion_temperature = Field2D{background_ion_temperature, bout_mesh}; + electron_temperature = Field2D{background_electron_temperature, bout_mesh}; + ion_velocity = std::vector{Field2D{background_ion_Vx, bout_mesh}, + Field2D{background_ion_Vy, bout_mesh}}; // RNG kernel // Used for sampling from velocity distribution for REC/CX // ------------------------------------------------------------------------------ @@ -712,32 +519,21 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) // Give particle group initial kinetic values (positions and velocities) // Numerical settings: weight, stdev, species ID + // use the same standard deviation for markers as in the initial distribution of velocities + // we should consider if marker distribution should evolve with time to track the neutral/ion temperature + // we should consider if marker distrubution should be initialised using Field2D information + const REAL initial_ion_thermal_speed = std::sqrt(background_ion_temperature / AA); ParticleSet maxwellian_markers = uniform_cellwise_maxwellian( - sycl_target, neso_mesh, particle_spec, rec_markers_per_cell, 1.0, 0.5, -1); + sycl_target, neso_mesh, particle_spec, rec_markers_per_cell, 1.0, + initial_ion_thermal_speed, -1); 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(); + // Give particle group initial fluid values: + // markers will contain background plasma properties + update_particle_properties_from_plasma(data_transfer, marker_group, ion_density, + ion_temperature, ion_velocity, + electron_density, electron_temperature); // Calculate marker weights @@ -774,8 +570,11 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) // Wrappers & controllers // ------------------------------------------------------------------------------ - this->source_manager = - std::make_unique(neso_mesh, bout_mesh, units); + this->source_manager = std::make_shared( + neso_mesh, data_transfer, bout_mesh, units); + // extract the units for the sources below + const BoutReal Nnorm = get(units["inv_meters_cubed"]); + const BoutReal Omega_ci = 1 / get(units["seconds"]); const REAL remove_threshold = options["remove_threshold"].withDefault(1.0e-10); const REAL merge_threshold = options["merge_threshold"].withDefault(1.0e-2); @@ -812,9 +611,10 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) 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", "m^-3 s^-1", Nnorm * Omega_ci, "Density source", + "Ionisation density source", accumulator_transform_iz, A_particle_group, + ion_source_density_zeroer); // Recombination transforms and controller @@ -830,9 +630,10 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* 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", "m^-3 s^-1", Nnorm * Omega_ci, "Density source", + "Recombination density source", accumulator_transform_rec, marker_group, + ion_source_density_zeroer); // Ionisation reaction // ------------------------------------------------------------------------------ @@ -975,27 +776,23 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) // 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}}); + // set weights from a constant initial density + set_initial_particle_weights(initial_neutral_density, A_particle_group, neso_mesh, + dof_kinetic_mesh_scalar, data_transfer, N_w); + // update particle properties from the plasma + update_particle_properties_from_plasma(data_transfer, A_particle_group, ion_density, + ion_temperature, ion_velocity, electron_density, + electron_temperature); + // write velocity moment diagnostics + diagnostics_manager = std::make_unique( + make_output_path("BOUT.dmp.vantage.particle.moments", alloptions), neso_mesh, + neso_mesh_cell_volumes_on_plasma_grid, A_particle_group, data_transfer, + this->source_manager, N_w, AA, bout_mesh, units, vantage_dump_filepath); + diagnostics_manager->update_kinetic_velocity_moments(); + diagnostics_manager->write_kinetic_velocity_moment_diagnostics(0, ion_density_kmsh); + diagnostics_manager->transfer_moments_to_plasma_grid(); + this->source_manager->update_all_sources(dt); // 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. @@ -1004,11 +801,17 @@ Vantage::Vantage(std::string name, Options& alloptions, Solver* solver) h5part->close(); // mass for conservation check - total_density = neutral_density + ion_density; - total_mass_initial = calculate_total_mass(total_density, neso_mesh); + neutral_density = diagnostics_manager->get_density_kinetic_mesh(); + // for (size_t ic=0; ic< static_cast(neso_mesh->get_cell_count());ic++){ + // total_density.at(ic) = neutral_density.at(ic) + ion_density_kmsh.at(ic); + // } + total_mass_initial = calculate_total_mass(neutral_density, neso_mesh); + total_mass_initial += + calculate_total_mass(ion_density, this->neso_mesh_cell_volumes_on_plasma_grid); // Initialise particle time particle_time = 0.0; + diagnostics_manager->write_bout_diagnostics(ion_density, particle_time); // Register VANTAGE timestep scheduler. // https://bout-dev.readthedocs.io/en/latest/user_docs/time_integration.html#monitoring-the-simulation-output @@ -1097,6 +900,15 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { aa = lambda_find_partial_moves(aa); } }; + // Create a ParticleSubGroup from particles that are in a cell with nonzero electron_density. + // This makes sure reactions are only applied where the neutrals are within the plasma volume + const REAL electron_density_threshold = this->electron_density_threshold; + ParticleSubGroupSharedPtr marker_group_in_plasma = + create_particle_sub_group_in_plasma_volume(marker_group, + electron_density_threshold); + ParticleSubGroupSharedPtr A_particle_group_in_plasma = + create_particle_sub_group_in_plasma_volume(A_particle_group, + electron_density_threshold); // begin timestepping output << "\nBegin VANTAGE iterations \n"; @@ -1107,24 +919,40 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { 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); + // update plasma properties on particles based on their new locations + update_particle_properties_from_plasma(data_transfer, A_particle_group, ion_density, + ion_temperature, ion_velocity, + electron_density, electron_temperature); + update_particle_properties_from_plasma(data_transfer, marker_group, ion_density, + ion_temperature, ion_velocity, + electron_density, electron_temperature); + // apply reactions to particles with a non-zero electron density property (those neutrals in the plasma) + reaction_controller->apply(A_particle_group_in_plasma, dt, + ControllerMode::standard_mode); + recombination_controller->apply(marker_group_in_plasma, 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"); + Field2D Siz = this->source_manager->get_plasma_grid_data("Siz"); + Field2D Srec = this->source_manager->get_plasma_grid_data("Srec"); + const std::vector Siz_kmsh = this->source_manager->get_kinetic_mesh_data("Siz"); + const std::vector Srec_kmsh = + this->source_manager->get_kinetic_mesh_data("Srec"); // "Solve" density // Sources are in normalised m^-3 s^-1, so need to multiply by dt ion_density += (Siz + Srec) * dt; - + data_transfer->transfer_scalar_to_kinetic_mesh(ion_density, ion_density_kmsh); + // for (size_t ic=0; ic < static_cast(neso_mesh->get_cell_count());ic++){ + // ion_density_kmsh.at(ic) += (Siz_kmsh.at(ic) + Srec_kmsh.at(ic)) * dt; + // } + + diagnostics_manager->update_kinetic_velocity_moments(); + diagnostics_manager->write_kinetic_velocity_moment_diagnostics(stepx + 1, + ion_density_kmsh); + diagnostics_manager->transfer_moments_to_plasma_grid(); // Write to VANTAGE dump files - update_diagnostics(neutral_density, ion_density, Siz, Srec, neso_mesh, - bout_output_data, *vantage_dump_writer, particle_time); - + // data_transfer->transfer_scalar_to_plasma_grid(ion_density_kmsh, ion_density); + diagnostics_manager->write_bout_diagnostics(ion_density, particle_time); // Write to particle_trajectories file h5part->write(); } @@ -1134,8 +962,10 @@ int Vantage::advance_vantage(BoutReal UNUSED(time)) { h5part->close(); // mass for conservation check - total_density = neutral_density + ion_density; - BoutReal total_mass_final = calculate_total_mass(total_density, neso_mesh); + neutral_density = diagnostics_manager->get_density_kinetic_mesh(); + REAL total_mass_final = calculate_total_mass(neutral_density, neso_mesh); + total_mass_final += + calculate_total_mass(ion_density, this->neso_mesh_cell_volumes_on_plasma_grid); if (test_mass_conservation) { check_mass_conservation(total_mass_final, total_mass_initial); } diff --git a/src/vantage_datatransfer.cxx b/src/vantage_datatransfer.cxx new file mode 100644 index 000000000..e0dbf066b --- /dev/null +++ b/src/vantage_datatransfer.cxx @@ -0,0 +1,179 @@ +#include "bout/bout.hxx" +#include +#include +#include +#include +#include "../include/vantage_datatransfer.hxx" + +using namespace NESO::Particles; + +// VANTAGE data transfer implementation +// ------------------------------------------------------------------------------ +VantageDataTransfer::VantageDataTransfer( + std::shared_ptr& neso_mesh, + std::shared_ptr& project_eval_dg0, + std::shared_ptr& mesh_coupler, + Mesh* bout_mesh, size_t ndim_vector) + : neso_mesh(neso_mesh), + project_eval_dg0(project_eval_dg0), + mesh_coupler(mesh_coupler), + bout_mesh(bout_mesh), + dof_kinetic_mesh_scalar(std::vector(static_cast(neso_mesh->get_cell_count()))), + ndim_vector(ndim_vector), + dof_kinetic_mesh_vector(std::vector(ndim_vector*static_cast(neso_mesh->get_cell_count()))) { + // local number of BOUT++ x cells, excluding guards + const int Nx = bout_mesh->xend - bout_mesh->xstart + 1; + // local number of BOUT++ y cells, excluding guards + const int Ny = bout_mesh->yend - bout_mesh->ystart + 1; + // Get the number of cells in the bout (plasma) mesh owned on this process, excluding guard cells + num_cells_owned_bout_grid = static_cast(Nx*Ny); + // a vector used to receive scalar BOUT++ data from the kinetic mesh + dof_bout_grid_scalar = std::vector(num_cells_owned_bout_grid); + dof_bout_grid_vector = std::vector(num_cells_owned_bout_grid*ndim_vector); + } + +void VantageDataTransfer::transfer_scalar_to_plasma_grid( + std::vector& scalar_kinetic_mesh, + Field2D& scalar_plasma_grid) { + // some ASSERT required here to check bout_mesh the same + if (this->mesh_coupler != nullptr){ + ASSERT1(scalar_kinetic_mesh.size() == static_cast(this->neso_mesh->get_cell_count())); + ASSERT1(scalar_kinetic_mesh.size() > this->dof_bout_grid_scalar.size()); + // we need to port data from the kinetic mesh dofs to the dofs expected by BOUT++ in the loop below + this->mesh_coupler->backward_transfer(scalar_kinetic_mesh, 1, dof_bout_grid_scalar); + } else { + ASSERT1(scalar_kinetic_mesh.size() == this->dof_bout_grid_scalar.size()); + this->dof_bout_grid_scalar = scalar_kinetic_mesh; + } + std::size_t ic = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + scalar_plasma_grid(ix, iy) = this->dof_bout_grid_scalar.at(ic); + ic++; + } + } + // this fills internal guards + this->bout_mesh->communicate(scalar_plasma_grid); + // apply boundary conditions to fill external guards + // scalar_field_plasma_grid.applyBoundary(); + // extrapolate -> Neumann +} + +void VantageDataTransfer::transfer_scalar_to_kinetic_mesh( + Field2D& scalar_plasma_grid, + std::vector& scalar_kinetic_mesh) { + // get scalar from plasma grid into the dummy vector + std::size_t ic = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + this->dof_bout_grid_scalar.at(ic) = scalar_plasma_grid(ix, iy); + ic++; + } + } + // transfer data from the dummy vector into the kinetic mesh dofs + if (this->mesh_coupler != nullptr){ + ASSERT1(scalar_kinetic_mesh.size() == static_cast(this->neso_mesh->get_cell_count())); + ASSERT1(scalar_kinetic_mesh.size() > this->dof_bout_grid_scalar.size()); + // we need to port data from the kinetic mesh dofs to the dofs expected by BOUT++ in the loop below + this->mesh_coupler->forward_transfer(dof_bout_grid_scalar, 1, scalar_kinetic_mesh); + } else { + ASSERT1(scalar_kinetic_mesh.size() == this->dof_bout_grid_scalar.size()); + scalar_kinetic_mesh = this->dof_bout_grid_scalar; + } +} + +void VantageDataTransfer::transfer_scalar_to_particle_property( + std::vector& scalar_kinetic_mesh, + std::shared_ptr& A_particle_group, + std::string particle_property){ + // check dimensions + ASSERT1(scalar_kinetic_mesh.size() == static_cast(this->neso_mesh->get_cell_count())) + // set the kinetic mesh property to NESO-Particles internal variables + this->project_eval_dg0->set_dofs(1, scalar_kinetic_mesh); + // set the data from internal variables into the weights + this->project_eval_dg0->evaluate(A_particle_group, Sym(particle_property)); +} + +void VantageDataTransfer::transfer_scalar_to_particle_property( + Field2D& scalar_plasma_grid, + std::shared_ptr& A_particle_group, + std::string particle_property){ + // check dimensions + // need some ASSERT to check bout_mesh is the same variables + this->transfer_scalar_to_kinetic_mesh(scalar_plasma_grid,this->dof_kinetic_mesh_scalar); + this->transfer_scalar_to_particle_property(this->dof_kinetic_mesh_scalar, A_particle_group, particle_property); +} + +void VantageDataTransfer::transfer_particle_property_to_scalar( + std::shared_ptr& A_particle_group, + std::string particle_property, + std::vector& scalar_kinetic_mesh){ + // check dimensions + ASSERT1(scalar_kinetic_mesh.size() == static_cast(this->neso_mesh->get_cell_count())) + // set the particle property to NESO-Particles internal variables + // some ASSERT to check particle property corresponds to a scalar? + this->project_eval_dg0->project(A_particle_group, Sym(particle_property)); + // project to the kinetic dof vector + project_eval_dg0->get_dofs(1, scalar_kinetic_mesh); +} + +void VantageDataTransfer::transfer_particle_property_to_vector( + std::shared_ptr& A_particle_group, + std::string particle_property, + std::vector& vector_kinetic_mesh){ + // check dimensions + ASSERT1(vector_kinetic_mesh.size() == this->ndim_vector*static_cast(this->neso_mesh->get_cell_count())) + // set the particle property to NESO-Particles internal variables + // some ASSERT to check particle property corresponds to a vector? + this->project_eval_dg0->project(A_particle_group, Sym(particle_property)); + // project to the kinetic dof vector + project_eval_dg0->get_dofs(static_cast(ndim_vector), vector_kinetic_mesh); +} + +void VantageDataTransfer::transfer_vector_to_particle_property( + std::vector& vector_kinetic_mesh, + std::shared_ptr& A_particle_group, + std::string particle_property){ + // check dimensions + ASSERT1(vector_kinetic_mesh.size() == this->ndim_vector*static_cast(this->neso_mesh->get_cell_count())) + // set the kinetic mesh property to NESO-Particles internal variables + this->project_eval_dg0->set_dofs(static_cast(this->ndim_vector), vector_kinetic_mesh); + // set the data from internal variables into the weights + this->project_eval_dg0->evaluate(A_particle_group, Sym(particle_property)); +} + +void VantageDataTransfer::transfer_vector_to_kinetic_mesh( + std::vector& vector_plasma_grid, std::vector& vector_kinetic_mesh){ + // check dimensions + ASSERT1(vector_plasma_grid.size() == this->ndim_vector); + std::size_t ic = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + for (size_t idim=0; idim < this->ndim_vector; idim++){ + const size_t jc = ic*this->ndim_vector + idim; + this->dof_bout_grid_vector.at(jc) = (vector_plasma_grid.at(idim))(ix, iy); + } + ic++; + } + } + // transfer data from the dummy vector into the kinetic mesh dofs + if (this->mesh_coupler != nullptr){ + ASSERT1(dof_kinetic_mesh_vector.size() == ndim_vector*static_cast(this->neso_mesh->get_cell_count())); + ASSERT1(dof_kinetic_mesh_vector.size() > this->dof_bout_grid_vector.size()); + // we need to port data from the kinetic mesh dofs to the dofs expected by BOUT++ in the loop below + this->mesh_coupler->forward_transfer( + this->dof_bout_grid_vector, + static_cast(this->ndim_vector), + this->dof_kinetic_mesh_vector); + } else { + ASSERT1(vector_kinetic_mesh.size() == this->dof_bout_grid_vector.size()); + vector_kinetic_mesh = this->dof_bout_grid_vector; + } +} +void VantageDataTransfer::transfer_vector_to_particle_property( + std::vector& vector_plasma_grid, + std::shared_ptr& A_particle_group, + std::string particle_property){ + this->transfer_vector_to_kinetic_mesh(vector_plasma_grid, this->dof_kinetic_mesh_vector); + this->transfer_vector_to_particle_property(this->dof_kinetic_mesh_vector, A_particle_group, particle_property); +} \ No newline at end of file diff --git a/src/vantage_diagnostics.cxx b/src/vantage_diagnostics.cxx new file mode 100644 index 000000000..dc0f14117 --- /dev/null +++ b/src/vantage_diagnostics.cxx @@ -0,0 +1,468 @@ +#include "bout/bout.hxx" +#include "bout/bout_types.hxx" +#include +#include +#include "../include/component.hxx" + +#include +#include +#include +#include +#include "../include/vantage_datatransfer.hxx" +#include "../include/vantage_diagnostics.hxx" +#include "../include/vantage_sources.hxx" + +using namespace NESO::Particles; + +// helper functions for diagnostics + +REAL calculate_total_mass(std::vector& density, + std::shared_ptr& neso_mesh) { + ASSERT1(density.size() == static_cast(neso_mesh->get_cell_count())); + REAL local_mass = 0.0; + REAL total_mass = 0.0; + for (size_t ic = 0;ic < static_cast(neso_mesh->get_cell_count()); ic++) { + local_mass += density.at(ic) * neso_mesh->dmh->get_cell_volume(static_cast(ic)); + } + MPICHK( + MPI_Allreduce(&local_mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, BoutComm::get())); + return total_mass; +} + +REAL calculate_total_mass(Field2D& density, + std::vector& neso_cell_volume_on_bout_mesh) { + Mesh* bout_mesh = density.getMesh(); + // local number of BOUT++ x cells, excluding guards + const int Nx = bout_mesh->xend - bout_mesh->xstart + 1; + // local number of BOUT++ y cells, excluding guards + const int Ny = bout_mesh->yend - bout_mesh->ystart + 1; + // Get the number of cells in the bout (plasma) mesh owned on this process, excluding guard cells + const size_t num_cells_owned_bout_mesh = static_cast(Nx*Ny); + // Get the number of cells in the kinetic (neutral) mesh owned on this process + ASSERT1(neso_cell_volume_on_bout_mesh.size() == num_cells_owned_bout_mesh); + // sum over the density on the BOUT++ grid, using NESO-Particles cell volumes + REAL local_mass = 0.0; + REAL total_mass = 0.0; + // sum over the ion density on the plasma mesh + size_t ixy=0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + local_mass += density(ix, iy) * neso_cell_volume_on_bout_mesh.at(ixy); + ixy++; + } + } + // sum contributions from different MPI ranks + MPICHK( + MPI_Allreduce(&local_mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, BoutComm::get())); + return total_mass; +} + +// helper function to initialise the plasma grid (BOUT++ mesh) diagnostics +Options initialise_plasma_grid_diagnostics(Options& units, Mesh* bout_mesh, + std::vector& neso_cell_volumes, + std::string vantage_dump_filepath) { + + // local number of BOUT++ x cells, excluding guards + const int Nx = bout_mesh->xend - bout_mesh->xstart + 1; + // local number of BOUT++ y cells, excluding guards + const int Ny = bout_mesh->yend - bout_mesh->ystart + 1; + // Get the number of cells in the bout (plasma) mesh owned on this process, excluding guard cells + const size_t num_cells_owned_bout_mesh = static_cast(Nx*Ny); + // Get the number of cells in the kinetic (neutral) mesh owned on this process + ASSERT1(neso_cell_volumes.size() == num_cells_owned_bout_mesh); + + // Options object to use to write out diagnostic data of fluid quantities + + const BoutReal Nnorm = get(units["inv_meters_cubed"]); + const BoutReal Tnorm = get(units["eV"]); + const BoutReal Omega_ci = 1 / get(units["seconds"]); + const BoutReal rho_s0 = get(units["meters"]); + const BoutReal Bnorm = get(units["Tesla"]); + const BoutReal Cs0 = get(units["meters"]) + / get(units["seconds"]); + + // save the area of each 2D cell where it aligns with the BOUT++ mesh + Field2D neso_cell_areas{0.0, bout_mesh}; + size_t ixy=0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + neso_cell_areas(ix,iy) = neso_cell_volumes.at(ixy); + ixy++; + } + } + + Options bout_output_data; + // Add metadata from mesh, e.g. branch cuts + bout_mesh->outputVars(bout_output_data); + // Add Rxy, Zxy coordinate data + Field2D Rxy; + Field2D Rxy_corners; + Field2D Rxy_lower_right_corners; + Field2D Rxy_upper_right_corners; + Field2D Rxy_upper_left_corners; + Field2D Zxy; + Field2D Zxy_corners; + Field2D Zxy_lower_right_corners; + Field2D Zxy_upper_right_corners; + Field2D Zxy_upper_left_corners; + // mesh->get(ivertex, "ivertex_lower_left_corners"); + bout_mesh->get(Rxy, "Rxy"); + bout_mesh->get(Rxy_corners, "Rxy_corners"); + bout_mesh->get(Rxy_lower_right_corners, "Rxy_lower_right_corners"); + bout_mesh->get(Rxy_upper_right_corners, "Rxy_upper_right_corners"); + bout_mesh->get(Rxy_upper_left_corners, "Rxy_upper_left_corners"); + bout_mesh->get(Zxy, "Zxy"); + bout_mesh->get(Zxy_corners, "Zxy_corners"); + bout_mesh->get(Zxy_lower_right_corners, "Zxy_lower_right_corners"); + bout_mesh->get(Zxy_upper_right_corners, "Zxy_upper_right_corners"); + bout_mesh->get(Zxy_upper_left_corners, "Zxy_upper_left_corners"); + set_with_attrs(bout_output_data["Rxy"], Rxy, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Rxy_corners"], Rxy_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Rxy_lower_right_corners"], Rxy_lower_right_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Rxy_upper_right_corners"], Rxy_upper_right_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Rxy_upper_left_corners"], Rxy_upper_left_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Zxy"], Zxy, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Zxy_corners"], Zxy_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Zxy_lower_right_corners"], Zxy_lower_right_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Zxy_upper_right_corners"], Zxy_upper_right_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["Zxy_upper_left_corners"], Zxy_upper_left_corners, { + {"units", "m"}, + {"conversion", 1}, // Already in SI units + }); + set_with_attrs(bout_output_data["neso_cell_areas"], neso_cell_areas, { + {"units", "m^2"}, + {"conversion", rho_s0*rho_s0}, // Already in SI units + }); + set_with_attrs(bout_output_data["y_boundary_guards"], 2, { + {"source", "vantage -- should be provided by BOUT++"} + }); + + // 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(vantage_dump_filepath)->write(bout_output_data); + return bout_output_data; +} + +// VANTAGE diagnostics manager implementation +// ------------------------------------------------------------------------------ +VantageDiagnosticsManager::VantageDiagnosticsManager( + std::string vtkhdf_filename, + std::shared_ptr& neso_mesh, + std::vector& neso_cell_volumes, + std::shared_ptr& A_particle_group, + std::shared_ptr& data_transfer, + std::shared_ptr& source_manager, + BoutReal N_w, BoutReal mass, Mesh* bout_mesh, + Options& units, std::string vantage_dump_filepath) + : vtkhdf_filename(vtkhdf_filename), + neso_mesh(neso_mesh), + A_particle_group(A_particle_group), + data_transfer(data_transfer), + source_manager(source_manager), + N_w(N_w), + mass(mass), + bout_mesh(bout_mesh), + units(units), + vantage_dump_filepath(vantage_dump_filepath) { + // initialise vectors for storing the moments on the kinetic mesh + const size_t ndimv = this->ndimv; + const std::size_t num_cells_owned_kinetic_mesh = static_cast(neso_mesh->get_cell_count()); + density = std::vector(num_cells_owned_kinetic_mesh); + energy = std::vector(num_cells_owned_kinetic_mesh); + gamma = std::vector(ndimv*num_cells_owned_kinetic_mesh); + uvector = std::vector(ndimv*num_cells_owned_kinetic_mesh); + pressure = std::vector(num_cells_owned_kinetic_mesh); + temperature = std::vector(num_cells_owned_kinetic_mesh); + // initialise BOUT++ diagnostic on BOUT++ mesh + // Object for VANTAGE dump files + vantage_dump_writer = + bout::OptionsIO::create({{"file", vantage_dump_filepath}, {"append", true}}); + bout_output_data = + initialise_plasma_grid_diagnostics(units, bout_mesh, neso_cell_volumes, vantage_dump_filepath); + // initialise plasma grid variables + density_plasma_grid = Field2D(0.0, bout_mesh); + energy_plasma_grid = Field2D(0.0, bout_mesh); + pressure_plasma_grid = Field2D(0.0, bout_mesh); + temperature_plasma_grid = Field2D(0.0, bout_mesh); + // local number of BOUT++ x cells, excluding guards + const int Nx = bout_mesh->xend - bout_mesh->xstart + 1; + // local number of BOUT++ y cells, excluding guards + const int Ny = bout_mesh->yend - bout_mesh->ystart + 1; + // Get the number of cells in the bout (plasma) mesh owned on this process, excluding guard cells + const size_t num_cells_owned_bout_mesh = static_cast(Nx*Ny); + // a vector used to receive scalar BOUT++ data from the kinetic mesh + dof_bout_mesh_scalar = std::vector(num_cells_owned_bout_mesh); + } + +// Functions for diagnostics on the kinetic mesh +void VantageDiagnosticsManager::update_kinetic_velocity_moments(){ + // get the necessary inputs from the class + std::shared_ptr neso_mesh = this->neso_mesh; + std::shared_ptr A_particle_group = this->A_particle_group; + BoutReal N_w = this->N_w; + BoutReal mass = this->mass; + // get the necessary private diagnostic variables + // vectors to hold the diagnosed moments + // use references here since we want to update the members + std::vector& density = this->density; + std::vector& energy = this->energy; + std::vector& gamma = this->gamma; + std::vector& uvector = this->uvector; + std::vector& pressure = this->pressure; + std::vector& temperature = this->temperature; + + // update the necessary particle properties for the moments + // define the lambda updating the moments + const size_t ndimv = this->ndimv; // number of velocity dimensions + const std::size_t num_cells_owned_kinetic_mesh = static_cast(neso_mesh->get_cell_count()); + auto lambda_update_moment_kernels = + [=](ParticleSubGroupSharedPtr aa) -> void { + particle_loop( + "update_moment_kernels", aa, + [=](auto VELOCITY, auto WEIGHT, auto WEIGHT_V2, auto WEIGHT_V) { + WEIGHT_V2.at(0) = WEIGHT.at(0) * (VELOCITY.at(0) * VELOCITY.at(0) + VELOCITY.at(1) * VELOCITY.at(1)); + for (int dim = 0; dim < static_cast(ndimv); dim++){ + WEIGHT_V.at(dim) = WEIGHT.at(0) * VELOCITY.at(dim); + } + }, + Access::read(Sym("VELOCITY")), + Access::read(Sym("WEIGHT")), + Access::write(Sym("WEIGHT_V2")), + Access::write(Sym("WEIGHT_V"))) + ->execute(); + }; + // call the particle loop + lambda_update_moment_kernels(static_particle_sub_group(A_particle_group)); + // extract density + this->data_transfer->transfer_particle_property_to_scalar( + A_particle_group, "WEIGHT", density); + // energy + this->data_transfer->transfer_particle_property_to_scalar( + A_particle_group, "WEIGHT_V2", energy); + // mean flow Gamma = nu + this->data_transfer->transfer_particle_property_to_vector( + A_particle_group, "WEIGHT_V", gamma); + // scalar variables + for (size_t ic=0; ic < num_cells_owned_kinetic_mesh; ic++){ + // multiply by any factors not handled in the project step + density.at(ic) *= N_w; + // (weight factor * mass / 2) + energy.at(ic) *= 0.5*N_w*mass; + } + // vector variables + for (size_t ic=0; ic < num_cells_owned_kinetic_mesh; ic++){ + for (size_t dim=0; dim < ndimv; dim++){ + const size_t jc = ic*ndimv + dim; // compound index covering all cells and dimensions + // multiply by any factors not handled in the project step + gamma.at(jc) *= N_w; + // obtain the derived quantity uvector + uvector.at(jc) = gamma.at(jc) / density.at(ic); + } + } + // derived scalar variables + for (size_t ic=0; ic < num_cells_owned_kinetic_mesh; ic++){ + // calculate pressure = (2E - m n u^2 ) / ndimv + pressure.at(ic) = 2.0*energy.at(ic); + for (size_t dim=0; dim < ndimv; dim++){ + const size_t jc = (ic*ndimv) + dim; + pressure.at(ic) -= mass*density.at(ic)*uvector.at(jc)*uvector.at(jc); + } + pressure.at(ic) /= static_cast(ndimv); + temperature.at(ic) = pressure.at(ic)/density.at(ic); + } +} + +// Function to save a VTKHDF file, writing the private member +// velocity moments and the mesh in VTK compatible format +// we pass in the ion_density here to enable testing of mass conservation +void VantageDiagnosticsManager::write_kinetic_velocity_moment_diagnostics(int istep, + std::vector& ion_density){ + // get the necessary inputs from the class + const std::string vtkhdf_filename = fmt::format("{}.istep.{}.vtkhdf",this->vtkhdf_filename,istep); + std::shared_ptr neso_mesh = this->neso_mesh; + + // write the data + VTK::VTKHDF vtk_writer(vtkhdf_filename, neso_mesh->get_comm()); + const std::size_t num_cells_owned_kinetic_mesh = static_cast(neso_mesh->get_cell_count()); + // vectors to hold the diagnosed moments + std::vector density = this->density; + std::vector energy = this->energy; + std::vector gamma = this->gamma; + std::vector uvector = this->uvector; + std::vector pressure = this->pressure; + std::vector temperature = this->temperature; + // mesh data only CellData not yet filled on each cell + std::vector dvtk0 = neso_mesh->dmh->get_vtk_cell_data(); + std::vector> cell_data(num_cells_owned_kinetic_mesh); + // scalar variables + for (size_t ic=0; ic < num_cells_owned_kinetic_mesh; ic++){ + // insert map entries at this ic + cell_data.at(ic).insert({"density", density.at(ic)}); + cell_data.at(ic).insert({"energy", energy.at(ic)}); + cell_data.at(ic).insert({"pressure", pressure.at(ic)}); + cell_data.at(ic).insert({"temperature", temperature.at(ic)}); + // write cell volume for convenience in later post-processing analysis + cell_data.at(ic).insert({"cellvolume", neso_mesh->dmh->get_cell_volume(static_cast(ic))}); + // write the "ion density" for testing purposes only + cell_data.at(ic).insert({"ion_density", ion_density.at(ic)}); + } + // vector variables + for (size_t ic=0; ic < num_cells_owned_kinetic_mesh; ic++){ + for (size_t dim=0; dim < ndimv; dim++){ + const size_t jc = ic*ndimv + dim; // compound index covering all cells and dimensions + // insert a map entry at this ic + cell_data.at(ic).insert({fmt::format("gamma_{}",dim), gamma.at(jc)}); + cell_data.at(ic).insert({fmt::format("uvector_{}",dim), uvector.at(jc)}); + } + } + // source variables (stored as scalars -> vector sources would require a refactor) + const std::vector source_names = this->source_manager->get_source_names(); + for (size_t is=0; is < source_names.size(); is++){ + const std::vector source_kinetic_mesh = this->source_manager->get_kinetic_mesh_data( + source_names.at(is)); + for (size_t ic=0; ic < num_cells_owned_kinetic_mesh; ic++){ + // insert map entries at this ic + cell_data.at(ic).insert({source_names.at(is), source_kinetic_mesh.at(ic)}); + } + } + for (size_t ic=0; ic < static_cast(num_cells_owned_kinetic_mesh); ic++){ + // fill the VTK::UnstructuredCell value appropriately + dvtk0.at(ic).cell_data = cell_data.at(ic); + } + vtk_writer.write(dvtk0); + vtk_writer.close(); +} + +void VantageDiagnosticsManager::transfer_moments_to_plasma_grid(){ + this->data_transfer->transfer_scalar_to_plasma_grid( + this->density, + this->density_plasma_grid); + this->data_transfer->transfer_scalar_to_plasma_grid( + this->energy, + this->energy_plasma_grid); + this->data_transfer->transfer_scalar_to_plasma_grid( + this->pressure, + this->pressure_plasma_grid); + this->data_transfer->transfer_scalar_to_plasma_grid( + this->temperature, + this->temperature_plasma_grid); +} + +void VantageDiagnosticsManager::write_bout_diagnostics(Field2D& ion_density, + // Field2D& Siz, Field2D& Srec, + BoutReal particle_time) { + // extract the units + const BoutReal Nnorm = get(this->units["inv_meters_cubed"]); + // const BoutReal Tnorm = get(units["eV"]); + // const BoutReal Omega_ci = 1 / get(this->units["seconds"]); + // const BoutReal rho_s0 = get(units["meters"]); + // const BoutReal Bnorm = get(units["Tesla"]); + // const BoutReal Cs0 = get(units["meters"]) + // / get(units["seconds"]); + Field2D neutral_density = this->density_plasma_grid; + set_with_attrs(this->bout_output_data["neutral_density"], neutral_density, + {{"time_dimension", "t"}}); + + set_with_attrs(this->bout_output_data["ion_density"], ion_density, {{"time_dimension", "t"}}); + + set_with_attrs(this->bout_output_data["Nn"], neutral_density, + {{"time_dimension", "t"}, + {"units", "m^-3"}, + {"conversion", Nnorm}, + {"standard_name", "Density"}, + {"long_name", "Kinetic neutral density"}, + {"species", "kinetic neutrals"}, + {"source", "vantage"}}); + // diagnose sources (scalars only -> vectors require a refactor) + const std::vector source_names = this->source_manager->get_source_names(); + for (size_t is=0; is < source_names.size(); is++){ + const Field2D source = this->source_manager->get_plasma_grid_data( + source_names.at(is)); + const std::string units_description = this->source_manager->get_units( + source_names.at(is)); + const BoutReal conversion = this->source_manager->get_conversion( + source_names.at(is)); + const std::string long_name = this->source_manager->get_long_name( + source_names.at(is)); + const std::string standard_name = this->source_manager->get_long_name( + source_names.at(is)); + set_with_attrs(this->bout_output_data[source_names.at(is)], source, + {{"time_dimension", "t"}, + {"units", units_description}, + {"conversion", conversion}, + {"standard_name", standard_name}, + {"long_name", long_name}, + {"species", "kinetic neutrals"}, + {"source", "vantage"}}); + } + + set_with_attrs(this->bout_output_data["t_array"], particle_time, {{"time_dimension", "t"}}); + + // Append data to file + this->vantage_dump_writer->write(this->bout_output_data); + // Ensure buffer is written to disk to avoid crash data loss + this->vantage_dump_writer->flush(); +} + +std::vector VantageDiagnosticsManager::get_density_kinetic_mesh(){ + return this->density; +} + + diff --git a/src/vantage_dmplex.cxx b/src/vantage_dmplex.cxx index 582d51061..476947798 100644 --- a/src/vantage_dmplex.cxx +++ b/src/vantage_dmplex.cxx @@ -6,7 +6,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -15,6 +17,7 @@ #include #include #include +#include #include "../include/vantage_dmplex.hxx" #ifndef NESO_PARTICLES_PETSC @@ -192,27 +195,45 @@ std::vector cells_definition_from_RZ_ivertex( return cells; } -DM create_dmplex_from_Bout_mesh(Mesh* bout_mesh, Options& mesh_options, - std::shared_ptr sycl_target, - std::string dmplex_h5_filename) { +void write_dmplex_to_file(DM dm, std::string dmplex_name, std::string dmplex_h5_filename){ + // save a HDF5 file containing the DM for diagnostics + PetscViewer viewer; + // Set a name for the DMPlex object (important for HDF5) + PetscObjectSetName(reinterpret_cast(dm), dmplex_name.c_str()); + // Create an HDF5 viewer + PetscViewerHDF5Open(BoutComm::get(), dmplex_h5_filename.c_str(), FILE_MODE_WRITE, + &viewer); + // Set viewer format to PETSC_VIEWER_HDF5_PETSC for compatibility + PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_PETSC); + // Save the DMPlex to the HDF5 file + DMView(dm, viewer); + // Clean up + PetscViewerDestroy(&viewer); + output << "Finished DMPlex diagnostic \n"; +} + +void create_dmplex_from_GMSH_msh(DM* dm, std::string msh_file){ + PETSCCHK(DMPlexCreateGmshFromFile(BoutComm::get(), msh_file.c_str(), + static_cast(1), dm)); +} + +void create_dmplex_from_Bout_mesh(DM* dm, Mesh* bout_mesh, Options& mesh_options, + std::shared_ptr sycl_target) { bool use_cxx_ivertex = mesh_options["use_cxx_ivertex"] - .doc("Use C++ based DMPlex creation routine instead of " + .doc("Use C++ based DMPlex creation routine instead of " "loading an external DMPlex? " "Default and recommendation is true.") - .withDefault(true); - std::string dmplex_name = mesh_options["dmplex_name"] - .doc("DMPlex object name.") - .withDefault("hypnotoad_dmplex_mesh"); + .withDefault(true); // DMPlex vertex distance tolerance for duplicate Hypnotoad vertices const BoutReal dmplex_vertex_tolerance = mesh_options["dmplex_vertex_tolerance"] .doc("Tolerance for determining duplicate vertices when creating DMPlex from " - "BOUT++ mesh.") + "BOUT++ mesh.") .withDefault(1.0e-8); output << fmt::format("Using option use_cxx_ivertex = {}", use_cxx_ivertex) - << std::endl; + << std::endl; Field2D Rxy_lower_left_corners; Field2D Rxy_lower_right_corners; Field2D Rxy_upper_right_corners; @@ -371,17 +392,17 @@ DM create_dmplex_from_Bout_mesh(Mesh* bout_mesh, Options& mesh_options, Field2D ivertex_upper_left_corners_cxx{-1, bout_mesh}; // now fill ivertex_corners arrays RZ_to_ivertex_vector(ivertex_lower_left_corners_cxx, global_Z_vertices, - global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_lower_left_corners, - Zxy_lower_left_corners); + global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_lower_left_corners, + Zxy_lower_left_corners); RZ_to_ivertex_vector(ivertex_lower_right_corners_cxx, global_Z_vertices, - global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_lower_right_corners, - Zxy_lower_right_corners); + global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_lower_right_corners, + Zxy_lower_right_corners); RZ_to_ivertex_vector(ivertex_upper_right_corners_cxx, global_Z_vertices, - global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_upper_right_corners, - Zxy_upper_right_corners); + global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_upper_right_corners, + Zxy_upper_right_corners); RZ_to_ivertex_vector(ivertex_upper_left_corners_cxx, global_Z_vertices, - global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_upper_left_corners, - Zxy_upper_left_corners); + global_R_vertices, dmplex_vertex_tolerance, bout_mesh, Rxy_upper_left_corners, + Zxy_upper_left_corners); // First we setup the topology of the mesh. PetscInt num_cells_owned = Nx * Ny; @@ -427,18 +448,18 @@ DM create_dmplex_from_Bout_mesh(Mesh* bout_mesh, Options& mesh_options, // int ivertex_minimum = mpi_rank * nvertex_per_process; // int ivertex_maximum = mpi_rank * nvertex_per_process + nvertex_this_process - 1; /* - * Each rank owns a contiguous block of global indices. We label our indices - * lexicographically (row-wise). Sorting out the global vertex indexing is - * probably one of the more tedious parts. - */ + * Each rank owns a contiguous block of global indices. We label our indices + * lexicographically (row-wise). Sorting out the global vertex indexing is + * probably one of the more tedious parts. + */ PetscInt num_vertices_owned = nvertex_this_process; /* - * Create the coordinates for the block of vertices we pass to petsc. For an - * existing mesh in memory this step will probably involve some MPI - * communication to gather the blocks of coordinates on the ranks which pass - * them to PETSc. - */ + * Create the coordinates for the block of vertices we pass to petsc. For an + * existing mesh in memory this step will probably involve some MPI + * communication to gather the blocks of coordinates on the ranks which pass + * them to PETSc. + */ std::vector vertex_coords(static_cast(num_vertices_owned * 2)); // shift due to differing rank size_t ishift; @@ -447,16 +468,13 @@ DM create_dmplex_from_Bout_mesh(Mesh* bout_mesh, Options& mesh_options, vertex_coords.at(iv * 2 + 0) = global_vertex_list_R[iv + ishift]; vertex_coords.at(iv * 2 + 1) = global_vertex_list_Z[iv + ishift]; } - // This DM will contain the DMPlex after we call the creation routine. - DM dm; // Create the DMPlex from the cells and coordinates. PETSCCHK(DMPlexCreateFromCellListParallelPetsc( BoutComm::get(), 2, num_cells_owned, num_vertices_owned, PETSC_DECIDE, 4, - PETSC_TRUE, cells.data(), 2, vertex_coords.data(), NULL, NULL, &dm)); - + PETSC_TRUE, cells.data(), 2, vertex_coords.data(), NULL, NULL, dm)); // Label all of the boundary faces with 100 in the "Face Sets" label by using // the helper function label_all_dmplex_boundaries. - PetscInterface::label_all_dmplex_boundaries(dm, PetscInterface::face_sets_label, 100); + // PetscInterface::label_all_dmplex_boundaries(dm, PetscInterface::face_sets_label, 100); // // Label subsections of the boundary by specifing pairs of vertices and using // // the label_dmplex_edges helper function. @@ -475,22 +493,101 @@ DM create_dmplex_from_Bout_mesh(Mesh* bout_mesh, Options& mesh_options, // PetscInterface::label_dmplex_edges(dm, PetscInterface::face_sets_label, // vertex_starts, vertex_ends, edge_labels); +} - // save a HDF5 file containing the DM for diagnostics - PetscViewer viewer; - // Set a name for the DMPlex object (important for HDF5) - PetscObjectSetName(reinterpret_cast(dm), dmplex_name.c_str()); - // Create an HDF5 viewer - PetscViewerHDF5Open(BoutComm::get(), dmplex_h5_filename.c_str(), FILE_MODE_WRITE, - &viewer); - // Set viewer format to PETSC_VIEWER_HDF5_PETSC for compatibility - PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_PETSC); - // Save the DMPlex to the HDF5 file - DMView(dm, viewer); - // Clean up - PetscViewerDestroy(&viewer); - output << "Finished DMPlex creation and diagnostic \n"; - return dm; +std::vector get_triangle_vertices(){ + // read data from netcdf for global vertices in mesh + // Open the NetCDF file in read-only mode + const std::string filename = Options::root()["mesh"]["file"]; + netCDF::NcFile dataFile(filename, netCDF::NcFile::read); + + // Get the vertices variable + // vertices is a global list of vertex coordinates + // std::string varName = "vertices"; + netCDF::NcVar dataVar_vertices = dataFile.getVar("vertices"); + NESOASSERT(!dataVar_vertices.isNull(), "vertices not found in file."); + std::vector dims_vertices = dataVar_vertices.getDims(); + size_t nvertices = dims_vertices[0].getSize(); + size_t ncomp = dims_vertices[1].getSize(); + // Read the data into a vector + std::vector vertices(nvertices*ncomp); + dataVar_vertices.getVar(vertices.data()); + + // close the netcdf file + dataFile.close(); + + return vertices; +} + +std::vector get_triangle_cell_definition(){ + // read data from netcdf for global vertices in mesh + // Open the NetCDF file in read-only mode + const std::string filename = Options::root()["mesh"]["file"]; + netCDF::NcFile dataFile(filename, netCDF::NcFile::read); + + // Get the tri_cell_vertices variable + // a list of integers defining each triangular cell + // in terms of indices that + // index the "vertices" list loaded above + // std::string varName_tri_cell = "tri_cell_vertices"; + netCDF::NcVar dataVar_tri_cell = dataFile.getVar("tri_cell_vertices"); + NESOASSERT(!dataVar_tri_cell.isNull(), "tri_cell_vertices not found in file."); + std::vector dims_tri_cell = dataVar_tri_cell.getDims(); + size_t ntriangle = dims_tri_cell[0].getSize(); + size_t ntricorners = dims_tri_cell[1].getSize(); + // Read the data into a vector + std::vector tri_cell_vertices(ntriangle*ntricorners); + dataVar_tri_cell.getVar(tri_cell_vertices.data()); + // close the netcdf file + dataFile.close(); + // std::cout << "tri_cell_verticies" << "\n"; + // for (size_t it=0; it < ntriangle; it++){ + // std::cout << fmt::format("local_cell.at({}): ",it); + // for (size_t iv=0; iv < 3; iv++){ + // std::cout << " " << tri_cell_vertices.at((it*3) + iv) << ", "; + // } + // std::cout << "\n "; + // } + return tri_cell_vertices; +} + +REAL get_triangle_area(size_t itriangle, + const std::vector& vertices, + const std::vector& tri_cell_vertices){ + // compute the area for this triangle + // use result of vector product for area + // A = 1/2 | u x v | + // where u and v are vectors defining two sides of the triangle + + // three vertices per triangle + const size_t ntri = 3; + std::vector local_cell(ntri); + // obtain the global vertex integers which define the local triangular cell + for (size_t iv=0; iv < local_cell.size(); iv++){ + local_cell.at(iv) = tri_cell_vertices.at((ntri*itriangle) + iv); + // std::cout << fmt::format("local_cell.at({}): ",iv) << local_cell.at(iv) << '\n'; + } + // std::cout << "local_cell: " << local_cell.data() << '\n'; + // expect two vector components per vertex, mesh is 2D + const size_t ncomp = 2; + std::vector local_vertices(ntri*ncomp); + for (size_t iv=0; iv < local_cell.size(); iv++){ + for (size_t ic=0; ic < ncomp; ic++){ + const size_t jc = (iv*ncomp) + ic; + local_vertices.at(jc) = vertices.at((static_cast(local_cell.at(iv))*ncomp) + ic); + // std::cout << fmt::format("local_vertices.at({}): ",jc) << local_vertices.at(jc) << '\n'; + } + } + const size_t iv0 = 0; + const size_t iv1 = 1; + const size_t iv2 = 2; + const REAL ux = local_vertices.at(iv1*ncomp) - local_vertices.at(iv0); + const REAL uy = local_vertices.at((iv1*ncomp) + 1) - local_vertices.at(iv0 + 1); + const REAL vx = local_vertices.at(iv2*ncomp) - local_vertices.at(iv0); + const REAL vy = local_vertices.at((iv2*ncomp) + 1) - local_vertices.at(iv0 + 1); + const REAL area = 0.5*std::abs((ux*vy) - (uy*vx)); + // std::cout << "area: " << area << '\n'; + return area; } #endif diff --git a/src/vantage_helperfunctions.cxx b/src/vantage_helperfunctions.cxx new file mode 100644 index 000000000..c03949882 --- /dev/null +++ b/src/vantage_helperfunctions.cxx @@ -0,0 +1,351 @@ +#include "../include/vantage_helperfunctions.hxx" +#include "../include/component.hxx" +#include "bout/bout.hxx" +#include "bout/bout_types.hxx" +#include + +using namespace NESO::Particles; + +size_t get_num_cells_owned_bout_mesh(Mesh*& bout_mesh) { + // local number of BOUT++ x cells, excluding guards + const int Nx = bout_mesh->xend - bout_mesh->xstart + 1; + // local number of BOUT++ y cells, excluding guards + const int Ny = bout_mesh->yend - bout_mesh->ystart + 1; + // Get the number of cells in the bout (plasma) mesh owned on this process, excluding guard cells + const size_t num_cells_owned_bout_mesh = static_cast(Nx * Ny); + return num_cells_owned_bout_mesh; +} + +std::shared_ptr +get_mesh_coupler_constant_weights(DM& dm, std::vector& kinetic_mesh_map, + Mesh*& bout_mesh, REAL backward_weight_0, + REAL backward_weight_1) { + const size_t num_cells_owned_bout_mesh = get_num_cells_owned_bout_mesh(bout_mesh); + std::vector> coupler_map_0( + static_cast(num_cells_owned_bout_mesh)); + Field2D map_RZ_to_itriangle_0; + Field2D map_RZ_to_itriangle_1; + bout_mesh->get(map_RZ_to_itriangle_0, "map_RZ_to_itriangle_0"); + bout_mesh->get(map_RZ_to_itriangle_1, "map_RZ_to_itriangle_1"); + int icell = 0; + for (int ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (int iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + // lower triangle + coupler_map_0.at(static_cast(icell)) + .push_back( + {kinetic_mesh_map.at(static_cast(map_RZ_to_itriangle_0(ix, iy))), + 1.0, backward_weight_0}); + // upper triangle + coupler_map_0.at(static_cast(icell)) + .push_back( + {kinetic_mesh_map.at(static_cast(map_RZ_to_itriangle_1(ix, iy))), + 1.0, backward_weight_1}); + icell += 1; + } + } + // object for transferring data between kinetic and bout mesh degree-of-freedom vectors + std::shared_ptr mesh_coupler = + std::make_shared(dm, coupler_map_0); + return mesh_coupler; +} + +std::vector get_cell_volumes_on_plasma_grid( + DM& dm, std::vector& kinetic_mesh_map, + std::shared_ptr& neso_mesh, Mesh*& bout_mesh) { + const size_t num_cells_owned_bout_mesh = get_num_cells_owned_bout_mesh(bout_mesh); + // Get the number of cells in the kinetic (neutral) mesh owned on this process + const size_t num_cells_owned_kinetic_mesh = + static_cast(neso_mesh->get_cell_count()); + // neso_mesh cell volumes on BOUT++ mesh indices + std::vector neso_cell_volumes_bmsh(num_cells_owned_bout_mesh); + // the checks + if (num_cells_owned_kinetic_mesh == num_cells_owned_bout_mesh) { + // zero the compound index + size_t ixy = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + neso_cell_volumes_bmsh.at(ixy) = + neso_mesh->dmh->get_cell_volume(static_cast(ixy)); + ixy++; + } + } + } else if (num_cells_owned_kinetic_mesh > num_cells_owned_bout_mesh) { + // assume that this corresponds to the case where the BOUT++ mesh is decomposed + // to triangles and there are also cells representing the region beyond the simulated plasma + // ------------------------------------------- + // first, make a mesh_coupler_dg0 object with unit weights + const std::shared_ptr mesh_coupler_unit_weight = + get_mesh_coupler_constant_weights(dm, kinetic_mesh_map, bout_mesh, 1.0, 1.0); + // obtain a list of kinetic mesh cell volumes + std::vector neso_cell_volumes_kmsh(num_cells_owned_kinetic_mesh); + for (size_t ic = 0; ic < num_cells_owned_kinetic_mesh; ic++) { + neso_cell_volumes_kmsh.at(ic) = + neso_mesh->dmh->get_cell_volume(static_cast(ic)); + } + // move these cell volumes to the bout mesh + mesh_coupler_unit_weight->backward_transfer(neso_cell_volumes_kmsh, 1, + neso_cell_volumes_bmsh); + } + return neso_cell_volumes_bmsh; +} + +std::vector get_cell_vertices_on_plasma_grid( + DM& dm, std::vector& kinetic_mesh_map, + std::shared_ptr& neso_mesh, Mesh*& bout_mesh) { + const size_t num_cells_owned_bout_mesh = get_num_cells_owned_bout_mesh(bout_mesh); + // Get the number of cells in the kinetic (neutral) mesh owned on this process + const size_t num_cells_owned_kinetic_mesh = + static_cast(neso_mesh->get_cell_count()); + // neso_mesh cell volumes on BOUT++ mesh indices + const size_t nquad_vertices = 4; + const size_t ntri_vertices = 3; + const size_t ndim = 2; // number of position coordinates expected + std::vector quad_cell_vertices_bmsh(nquad_vertices * ndim + * num_cells_owned_bout_mesh); + // get the cell vertices in flattened vectors, + // without attempting to respect anti-clockwise vertex ordering + if (num_cells_owned_kinetic_mesh == num_cells_owned_bout_mesh) { + std::vector> cell_vertices; + // zero the compound index + size_t ixy = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + neso_mesh->dmh->get_cell_vertices(static_cast(ixy), cell_vertices); + for (size_t iv = 0; iv < nquad_vertices; iv++) { + for (size_t idim = 0; idim < ndim; idim++) { + const size_t jc = (ndim * ((nquad_vertices * ixy) + iv)) + idim; + quad_cell_vertices_bmsh.at(jc) = cell_vertices.at(iv).at(idim); + } + } + ixy++; + } + } + } else if (num_cells_owned_kinetic_mesh > num_cells_owned_bout_mesh) { + // assume that this corresponds to the case where the BOUT++ mesh is decomposed + // to triangles and there are also cells representing the region beyond the simulated plasma + // ------------------------------------------- + // we need to get the triangular cell coordinates from each upper and lower triangle + // on to the local BOUT++ grid, then resolve which coordinates are unique to form + // the coordinates for the quadrilateral cell which the pair of triangles represent + // ------------------------------------------- + // first, make a mesh_coupler_dg0 object with unit weights from the lower triangle, and zero weight + // for the upper triangle + const std::shared_ptr mesh_coupler_0 = + get_mesh_coupler_constant_weights(dm, kinetic_mesh_map, bout_mesh, 1.0, 0.0); + // second, make a mesh_coupler_dg0 object with unit weights from the upper triangle, and zero weight + // for the lower triangle + const std::shared_ptr mesh_coupler_1 = + get_mesh_coupler_constant_weights(dm, kinetic_mesh_map, bout_mesh, 0.0, 1.0); + // obtain the cell coordinates for lower and upper triangles on the kinetic mesh + std::vector> cell_vertices; + std::vector tri_cell_vertices_kmsh(ntri_vertices * ndim + * num_cells_owned_kinetic_mesh); + for (size_t ic = 0; ic < num_cells_owned_kinetic_mesh; ic++) { + neso_mesh->dmh->get_cell_vertices(static_cast(ic), cell_vertices); + // fill in results to flattened vector + for (size_t iv = 0; iv < ntri_vertices; iv++) { + for (size_t idim = 0; idim < ndim; idim++) { + const size_t jc = (ndim * ((ntri_vertices * ic) + iv)) + idim; + tri_cell_vertices_kmsh.at(jc) = cell_vertices.at(iv).at(idim); + } + } + } + // transfer these results to vectors for the lower and upper triangles + std::vector tri_cell_vertices_0_bmsh(ntri_vertices * ndim + * num_cells_owned_bout_mesh); + std::vector tri_cell_vertices_1_bmsh(ntri_vertices * ndim + * num_cells_owned_bout_mesh); + mesh_coupler_0->backward_transfer(tri_cell_vertices_kmsh, ntri_vertices * ndim, + tri_cell_vertices_0_bmsh); + mesh_coupler_1->backward_transfer(tri_cell_vertices_kmsh, ntri_vertices * ndim, + tri_cell_vertices_1_bmsh); + // fill in data for quad cell vertices + // no requirement for the cell centre check to list in anti-clockwise order + size_t ixy = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + // first three vertices from lower triangle are definitely unqiue vertices for the quad + // (though perhaps in an incorrect order) + for (size_t iv = 0; iv < ntri_vertices; iv++) { + for (size_t idim = 0; idim < ndim; idim++) { + const size_t jc_quad = (ndim * ((nquad_vertices * ixy) + iv)) + idim; + const size_t jc_tri = (ndim * ((ntri_vertices * ixy) + iv)) + idim; + quad_cell_vertices_bmsh.at(jc_quad) = tri_cell_vertices_0_bmsh.at(jc_tri); + } + } + // the final unique coordinate must be determined by checking for uniqueness + const size_t ivquad = 3; + const REAL atol = 1.0e-12; + std::vector unique(ntri_vertices); + for (size_t ivp = 0; ivp < ntri_vertices; ivp++) { + const size_t jcp_tri = (ndim * ((ntri_vertices * ixy) + ivp)); + // initially presume that this index is unique + unique.at(ivp) = true; + for (size_t iv = 0; iv < ntri_vertices; iv++) { + const size_t jc_tri = (ndim * ((ntri_vertices * ixy) + iv)); + REAL sumsqr = 0.0; + // sum the squared lengths measuring the distance of this vertex from another + for (size_t idim = 0; idim < ndim; idim++) { + sumsqr += std::pow(tri_cell_vertices_0_bmsh.at(jc_tri + idim) + - tri_cell_vertices_1_bmsh.at(jcp_tri + idim), + 2); + } + const REAL l2norm = std::sqrt(sumsqr); + if (l2norm < atol) { + unique.at(ivp) = false; + } + } + if (unique.at(ivp)) { + // this vertex has proved to be unique by not matching any other vertex + for (size_t idim = 0; idim < ndim; idim++) { + const size_t jc_quad = (ndim * ((nquad_vertices * ixy) + ivquad)) + idim; + quad_cell_vertices_bmsh.at(jc_quad) = + tri_cell_vertices_1_bmsh.at(jcp_tri + idim); + } + // only one vertex can be unique + break; + } + } + ixy++; + } + } + } + return quad_cell_vertices_bmsh; +} + +void check_cell_volumes(std::vector neso_cell_volumes_bmsh, Mesh*& bout_mesh, + Options& alloptions) { + Coordinates* coord = bout_mesh->getCoordinates(); + size_t ixy = 0; + const REAL tolerance = 1.0e-12; + const size_t num_cells_owned_bout_mesh = get_num_cells_owned_bout_mesh(bout_mesh); + // Get the number of cells in the kinetic (neutral) mesh owned on this process + ASSERT1(neso_cell_volumes_bmsh.size() == num_cells_owned_bout_mesh); + // dimensional units + const BoutReal meters = get(alloptions["units"]["meters"]); + const BoutReal meters_squared = meters * meters; + const BoutReal meters_cubed = meters * meters * meters; + // the checks of cell volumes + // zero the compound index + ixy = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + // Convert to SI: dx is m^2 T, J is m/T, dy is unitless, skip dz + // so J * dx * dy = m^3, technically per radian toroidal angle due to missing dz + const BoutReal bout_cell_area = + coord->J(ix, iy) * coord->dx(ix, iy) * coord->dy(ix, iy) * meters_cubed; + // neso_mesh is a 2D grid, needs m^2 + const REAL neso_cell_area = neso_cell_volumes_bmsh.at(ixy) * meters_squared; + const bool volumes_match = (abs(bout_cell_area - neso_cell_area) < tolerance); + // exit if we fail to find a match + NESOASSERT(volumes_match, + fmt::format("BOUT++ mesh volume {} does not match NESO-Particles mesh " + "volume {} for ix = {} iy = {} \n Ignore this message by " + "setting [dmplex] test_dmplex_cell_volumes = false", + bout_cell_area, neso_cell_area, ix, iy)); + ixy++; + } + } +} + +REAL cell_length(std::vector>& cell_vertices, std::size_t iv1, + std::size_t iv2, std::size_t iv3, std::size_t iv4) { + const REAL Rlength2 = + std::pow(0.5 + * (cell_vertices.at(iv1).at(0) + cell_vertices.at(iv2).at(0) + - cell_vertices.at(iv3).at(0) - cell_vertices.at(iv4).at(0)), + 2.0); + const REAL Zlength2 = + std::pow(0.5 + * (cell_vertices.at(iv1).at(1) + cell_vertices.at(iv2).at(1) + - cell_vertices.at(iv3).at(1) - cell_vertices.at(iv4).at(1)), + 2.0); + REAL length = std::pow(Zlength2 + Rlength2, 0.5); + return length; +} + +void check_cell_centres(Options& alloptions, DM& dm, + std::vector& kinetic_mesh_map, + std::shared_ptr& neso_mesh, + Mesh*& bout_mesh, BoutReal absolute_tolerance, + BoutReal relative_tolerance) { + // get (R,Z) of cell centres in Hypnotoad grid + Field2D Rxy; + Field2D Zxy; + bout_mesh->get(Rxy, "Rxy"); + bout_mesh->get(Zxy, "Zxy"); + + BoutReal meters = get(alloptions["units"]["meters"]); + std::vector neso_cell_vertices_plasma_grid = + get_cell_vertices_on_plasma_grid(dm, kinetic_mesh_map, neso_mesh, bout_mesh); + // number of vertices per quad + const size_t nquad_vertices = 4; + // expected dimensionality + const size_t ndim = 2; + // compare to cell centres calculated from cell corners + std::vector> cell_vertices(nquad_vertices, std::vector(ndim)); + + PetscInt ixy = 0; + for (PetscInt ix = bout_mesh->xstart; ix <= bout_mesh->xend; ix++) { + for (PetscInt iy = bout_mesh->ystart; iy <= bout_mesh->yend; iy++) { + const REAL bout_Rxy = Rxy(ix, iy); + const REAL bout_Zxy = Zxy(ix, iy); + + // fill in the vertices from the flattened vector + for (std::size_t iv = 0; iv < nquad_vertices; iv++) { + for (size_t idim = 0; idim < ndim; idim++) { + const size_t jc_quad = + (ndim * ((nquad_vertices * static_cast(ixy)) + iv)) + idim; + cell_vertices.at(iv).at(idim) = neso_cell_vertices_plasma_grid.at(jc_quad); + } + } + REAL neso_Rxy = 0.0; + REAL neso_Zxy = 0.0; + for (std::size_t iv = 0; iv < nquad_vertices; iv++) { + // DMPlex is stored in normalised units, need conversion to [m] + neso_Rxy += cell_vertices.at(iv).at(0) * meters; + neso_Zxy += cell_vertices.at(iv).at(1) * meters; + } + neso_Rxy /= 4.0; + neso_Zxy /= 4.0; + // get lengths of cell across the two dimensions + const REAL cell_length_a = cell_length(cell_vertices, 0, 1, 2, 3) * meters; + const REAL cell_length_b = cell_length(cell_vertices, 0, 3, 2, 1) * meters; + const REAL min_cell_length = std::min(cell_length_a, cell_length_b); + // we compare the difference in cell centres to the absolute tolerance and + // the relative tolerance formed by comparing to the smallest length across the cell + const REAL tolerance = absolute_tolerance + min_cell_length * relative_tolerance; + const bool centres_match = (abs(neso_Rxy - bout_Rxy) < tolerance) + && (abs(neso_Zxy - bout_Zxy) < tolerance); + // exit if we fail to find a match + NESOASSERT( + centres_match, + fmt::format("Hypnotoad/BOUT++ cell centre (R, Z) ({}, {}) does not match " + "NESO-Particles mesh inferred quad " + "cell centre ({}, {}) for ix = {} iy = {} \n" + "The cell height and width are {} {} \n" + "The displacements in R and Z are {} {} \n" + "Ignore this message by " + "setting [dmplex] test_dmplex_cell_centres = false\n Relax the " + "tolerance used in this check by increasing\n" + "[dmplex] dmplex_cell_centre_absolute_tolerance = {}\n" + "[dmplex] dmplex_cell_centre_relative_tolerance = {}", + bout_Rxy, bout_Zxy, neso_Rxy, neso_Zxy, ix, iy, cell_length_a, + cell_length_b, abs(neso_Rxy - bout_Rxy), abs(neso_Zxy - bout_Zxy), + absolute_tolerance, relative_tolerance)); + ixy++; + } + } +} + +void check_mass_conservation(REAL total_mass_final, REAL total_mass_initial) { + REAL rtol = 1.0e-13; + REAL mass_conserved = + (abs(total_mass_final - total_mass_initial) < rtol * total_mass_initial); + // exit if we fail to find conservation + NESOASSERT(mass_conserved, + fmt::format("Initial total mass {} does not match " + "final total mass {} \n Ignore this message by " + "setting [vantage] test_mass_conservation = false", + total_mass_initial, total_mass_final)); +} diff --git a/src/vantage_sources.cxx b/src/vantage_sources.cxx new file mode 100644 index 000000000..0daf71686 --- /dev/null +++ b/src/vantage_sources.cxx @@ -0,0 +1,133 @@ +#include "bout/bout.hxx" +#include "bout/bout_types.hxx" +#include +#include "../include/component.hxx" +#include +#include +#include +#include "../include/vantage_sources.hxx" +#include +#include +#include +#include "../include/vantage_datatransfer.hxx" + + +using namespace NESO::Particles; +using namespace VANTAGE::Reactions; + +// VANTAGE source manager implementation +// ------------------------------------------------------------------------------ +VantageSourceManager::VantageSourceManager( + std::shared_ptr& neso_mesh, + std::shared_ptr& data_transfer, + Mesh* bout_mesh, + Options& units) + : bout_mesh(bout_mesh), + neso_mesh(neso_mesh), + data_transfer(data_transfer), + units(units) {} + +// Register new source with the manager and initialise its data +void VantageSourceManager::add_source( + const std::string& hermes_source_name, const std::string& vantage_source_name, + const std::string& bout_diagnostic_units, + const BoutReal bout_diagnostic_conversion, + const std::string& bout_diagnostic_standard_name, + const std::string& bout_diagnostic_long_name, + std::shared_ptr> accumulator, + std::shared_ptr particle_group, + std::shared_ptr zeroer) { + + Field2D source_data_plasma_grid{bout_mesh}; + source_data_plasma_grid = 0.0; + const int num_cells_owned_kinetic_mesh = neso_mesh->get_cell_count(); + std::vector source_data_kinetic_mesh(static_cast(num_cells_owned_kinetic_mesh)); + + VantageSource source{ + hermes_source_name, + vantage_source_name, + bout_diagnostic_units, + bout_diagnostic_conversion, + bout_diagnostic_standard_name, + bout_diagnostic_long_name, + accumulator, particle_group, zeroer, + source_data_plasma_grid, + source_data_kinetic_mesh}; + + this->sources[hermes_source_name] = source; +} + +// Return source data on plasma grid +Field2D VantageSourceManager::get_plasma_grid_data(const std::string& hermes_source_name) { + return this->sources[hermes_source_name].source_data_plasma_grid; +} + +// Return source data on kinetic mesh +std::vector VantageSourceManager::get_kinetic_mesh_data(const std::string& hermes_source_name) { + return this->sources[hermes_source_name].source_data_kinetic_mesh; +} + +// Return diagnostic units data +std::string VantageSourceManager::get_units(const std::string& hermes_source_name) { + return this->sources[hermes_source_name].bout_diagnostic_units; +} + +// Return diagnostic units data +BoutReal VantageSourceManager::get_conversion(const std::string& hermes_source_name) { + return this->sources[hermes_source_name].bout_diagnostic_conversion; +} + +// Return diagnostic units data +std::string VantageSourceManager::get_standard_name(const std::string& hermes_source_name) { + return this->sources[hermes_source_name].bout_diagnostic_standard_name; +} + +// Return diagnostic units data +std::string VantageSourceManager::get_long_name(const std::string& hermes_source_name) { + return this->sources[hermes_source_name].bout_diagnostic_long_name; +} + +// Update the source from VANTAGE and reset the VANTAGE data/accumulator +void VantageSourceManager::update_source(const std::string& hermes_source_name, + double dt) { + + VantageSource& source = this->sources[hermes_source_name]; + BoutReal N_w = get(units["N_w"]); + + std::vector> accumulated_1d = + source.accumulator->get_cell_data(source.vantage_source_name); + size_t naccumulated = accumulated_1d.size(); + ASSERT1(naccumulated == source.source_data_kinetic_mesh.size()); + for (size_t ic = 0; ic < naccumulated; ic++){ + source.source_data_kinetic_mesh.at(ic) = accumulated_1d[ic]->at(0, 0) // Total weight + * N_w // Total particles + / neso_mesh->dmh->get_cell_volume(static_cast(ic)) // Total density + / dt; // Density source; + } + // copy accumulated data into the BOUT++ Field2D variable + this->data_transfer->transfer_scalar_to_plasma_grid( + source.source_data_kinetic_mesh, source.source_data_plasma_grid); + // Reset the accumulator object + source.accumulator->zero_buffer(source.vantage_source_name); + // Reset the accumulated source data on the particle + source.zeroer->transform(std::make_shared(source.particle_group)); +} + +// Update all sources +void VantageSourceManager::update_all_sources(double dt) { + for (auto& [hermes_source_name, source] : this->sources) { + update_source(hermes_source_name, dt); + } +} + +// get list of (hermes-3) source names +std::vector VantageSourceManager::get_source_names() { + const size_t nsources=this->sources.size(); + std::vector source_names(nsources); + size_t is=0; + for (auto& [hermes_source_name, source] : this->sources) { + source_names.at(is) = hermes_source_name; + is += 1; + } + return source_names; +} diff --git a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp index b238b427d..c94bab4ce 100644 --- a/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp +++ b/tests/integrated/dmplex-vertex-coordinates/data/BOUT.inp @@ -37,7 +37,6 @@ dt = 0.005 nsteps = 3 test_mass_conservation = true -initial_neutral_density = 20 npart_per_cell = 20 remove_threshold = 0.0 diff --git a/tests/integrated/particle-pusher/extended_kinetic_mesh_data/BOUT.dmp.vantage.0.kinetic.msh b/tests/integrated/particle-pusher/extended_kinetic_mesh_data/BOUT.dmp.vantage.0.kinetic.msh new file mode 100644 index 000000000..f87d62453 --- /dev/null +++ b/tests/integrated/particle-pusher/extended_kinetic_mesh_data/BOUT.dmp.vantage.0.kinetic.msh @@ -0,0 +1,319 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +0 0 1 0 +1 -1 0 0 2 6.283185307179586 0 0 0 +$EndEntities +$Nodes +1 87 1 87 +2 1 0 87 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +32 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +0 0 0 +0 1.570796326794897 0 +0 3.141592653589793 0 +0 4.71238898038469 0 +0.25 0 0 +0.25 1.570796326794897 0 +0.25 3.141592653589793 0 +0.25 4.71238898038469 0 +0.5 0 0 +0.5 1.570796326794897 0 +0.5 3.141592653589793 0 +0.5 4.71238898038469 0 +0.75 0 0 +0.75 1.570796326794897 0 +0.75 3.141592653589793 0 +0.75 4.71238898038469 0 +1 0 0 +1 1.570796326794897 0 +1 3.141592653589793 0 +1 4.71238898038469 0 +0.25 6.283185307179586 0 +0.5 6.283185307179586 0 +0.75 6.283185307179586 0 +1 6.283185307179586 0 +0 6.283185307179586 0 +-1 6.283185307179586 0 +-1 0 0 +2 0 0 +2 6.283185307179586 0 +-0.4999999999986942 6.283185307179586 0 +-1 5.799863360474505 0 +-1 5.316541413769544 0 +-1 4.833219467065129 0 +-1 4.349897520360367 0 +-1 3.866575573654959 0 +-1 3.383253626950145 0 +-1 2.899931680244037 0 +-1 2.41660973353672 0 +-1 1.933287786829305 0 +-1 1.449965840122101 0 +-1 0.9666438934146182 0 +-1 0.4833219467074823 0 +-0.500000000002059 0 0 +1.5 0 0 +2 0.4833219467051673 0 +2 0.9666438934102256 0 +2 1.44996584011498 0 +2 1.933287786819744 0 +2 2.416609733524625 0 +2 2.899931680229439 0 +2 3.383253626935547 0 +2 3.866575573642864 0 +2 4.349897520350281 0 +2 4.833219467057484 0 +2 5.316541413764967 0 +2 5.799863360472104 0 +1.5 6.283185307179586 0 +-0.6071522086797074 0.784761011285081 0 +-0.607152208680398 5.498424295897261 0 +-0.5271436386483882 2.235364003519466 0 +-0.5271436386490111 4.047821303667518 0 +-0.6377166249249197 1.29275877701848 0 +-0.6377166249255449 4.990426530165874 0 +-0.581430915947327 4.591558493712746 0 +-0.6217149109190194 2.767018144897421 0 +-0.5814309159448909 1.691626813475703 0 +-0.621714910919144 3.516167162291901 0 +-0.6486859643676327 3.14159265359466 0 +-0.6960731294749422 5.95330761966039 0 +-0.3642360307921816 5.790607309349078 0 +-0.696073129475413 0.3298776875203251 0 +-0.3642360307925395 0.4925779978312037 0 +1.607152208679471 5.498424295894448 0 +1.60715220868047 0.7847610112824089 0 +1.527143638648388 4.047821303660118 0 +1.527143638649029 2.235364003512242 0 +1.637716624924872 4.990426530161094 0 +1.63771662492552 1.292758777013975 0 +1.581430915947132 1.691626813467362 0 +1.621714910919019 3.516167162282163 0 +1.581430915944891 4.59155849370388 0 +1.621714910919148 2.767018144887718 0 +1.648685964367633 3.141592653584932 0 +1.696073129475332 0.3298776875191727 0 +1.364236030792596 0.4925779978304361 0 +1.696073129474763 5.953307619659157 0 +1.364236030791846 5.790607309348226 0 +$EndNodes +$Elements +1 130 1 130 +2 1 2 130 +1 1 5 6 +2 1 6 2 +3 2 6 7 +4 2 7 3 +5 3 7 8 +6 3 8 4 +7 4 8 21 +8 4 21 25 +9 5 9 10 +10 5 10 6 +11 6 10 11 +12 6 11 7 +13 7 11 12 +14 7 12 8 +15 8 12 22 +16 8 22 21 +17 9 13 14 +18 9 14 10 +19 10 14 15 +20 10 15 11 +21 11 15 16 +22 11 16 12 +23 12 16 23 +24 12 23 22 +25 13 17 18 +26 13 18 14 +27 14 18 19 +28 14 19 15 +29 15 19 20 +30 15 20 16 +31 16 20 24 +32 16 24 23 +33 2 60 66 +34 61 4 64 +35 41 42 58 +36 31 32 59 +37 40 41 62 +38 32 33 63 +39 30 26 69 +40 27 43 71 +41 39 40 66 +42 33 34 64 +43 38 39 60 +44 34 35 61 +45 37 38 65 +46 36 37 68 +47 35 36 67 +48 42 27 71 +49 26 31 69 +50 58 42 71 +51 31 59 69 +52 41 58 62 +53 59 32 63 +54 3 67 68 +55 65 3 68 +56 40 62 66 +57 63 33 64 +58 60 39 66 +59 34 61 64 +60 4 63 64 +61 62 2 66 +62 38 60 65 +63 61 35 67 +64 37 65 68 +65 67 36 68 +66 58 71 72 +67 69 59 70 +68 43 1 72 +69 25 30 70 +70 71 43 72 +71 30 69 70 +72 1 2 72 +73 2 3 60 +74 2 58 72 +75 60 3 65 +76 3 61 67 +77 59 4 70 +78 58 2 62 +79 3 4 61 +80 4 25 70 +81 4 59 63 +82 20 75 81 +83 76 18 79 +84 55 56 73 +85 45 46 74 +86 54 55 77 +87 46 47 78 +88 44 28 84 +89 29 57 86 +90 53 54 81 +91 47 48 79 +92 52 53 75 +93 48 49 76 +94 51 52 80 +95 50 51 83 +96 49 50 82 +97 56 29 86 +98 28 45 84 +99 73 56 86 +100 45 74 84 +101 55 73 77 +102 74 46 78 +103 19 82 83 +104 80 19 83 +105 54 77 81 +106 78 47 79 +107 75 53 81 +108 48 76 79 +109 18 78 79 +110 77 20 81 +111 52 75 80 +112 76 49 82 +113 51 80 83 +114 82 50 83 +115 73 86 87 +116 84 74 85 +117 17 44 85 +118 57 24 87 +119 86 57 87 +120 44 84 85 +121 24 20 87 +122 20 19 75 +123 20 73 87 +124 75 19 80 +125 19 76 82 +126 74 18 85 +127 73 20 77 +128 19 18 76 +129 18 17 85 +130 18 74 78 +$EndElements diff --git a/tests/integrated/particle-pusher/extended_kinetic_mesh_tools.py b/tests/integrated/particle-pusher/extended_kinetic_mesh_tools.py new file mode 100644 index 000000000..a13b6fc21 --- /dev/null +++ b/tests/integrated/particle-pusher/extended_kinetic_mesh_tools.py @@ -0,0 +1,446 @@ +from boututils.run_wrapper import shell, launch_safe +from netCDF4 import Dataset +import numpy as np + + +def kinetic_mesh_path(): + return "extended_kinetic_mesh_data/BOUT.dmp.vantage.0.kinetic.msh" + + +def extended_kinetic_mesh_test_input( + bout_grid_file, + msh_file, + dt, + nsteps, + iz_rate, + rec_rate, + remove_threshold, + merge_threshold, +): + + input_file_string = f""" + nout = 1 + timestep = 1 + + [mesh] + file="{bout_grid_file}" + extrapolate_y=false + extrapolate_x=false + + [dmplex] + test_dmplex_cell_volumes = true + test_dmplex_cell_centres = true + use_external_msh = true + msh_file = "{msh_file}" + + [solver] + type = pvode + + [hermes] + components = (d+, e, vantage) + Nnorm = 1e19 + normalise_metric=false + + [d+] + type = evolve_density + AA = 1 + charge = 1 + + [Nd+] + function = 1 + + [e] + type = quasineutral + AA = 1/1836 + charge = -1 + + [vantage] + dt = {dt} + nsteps = {nsteps} + test_mass_conservation = true + initial_neutral_pressure = 1 + initial_neutral_temperature = 1 + npart_per_cell = 20 + background_ion_density = 1e19 + background_ion_temperature = 50 + background_ion_Vy = 1 + remove_threshold = {remove_threshold} + merge_threshold = {merge_threshold} + iz_rate_override = {iz_rate} + rec_rate_override = {rec_rate} + """ + return input_file_string + + +def generate_BOUT_grid_data(base_grid_dir, kinetic_nc_file_path, verbose=True): + # make the directory for the basic hermes-3 run which + # makes a slab grid + cmd = f"mkdir {base_grid_dir}" + shell(cmd) + # this input file grid resolutions here + # cannot be modified without changing + # the mesh variables + # map_RZ_to_itriangle_0 + # map_RZ_to_itriangle_1 + # vertices + # tri_cell_vertices + input_file_string = """ + nout = 0 + timestep = 1 + + [mesh] + J = 1 + + nx = 8 # X grid size + ny = 4 # Y grid size + + dx = 1.0/(nx-4) # X mesh spacing + dy = 2*pi/ny # Y mesh spacing + dz = 1 # Unity in toroidal direction + + Rxy = x + Rxy_corners = x - 0.5*dx + Rxy_lower_right_corners = x + 0.5*dx + Rxy_upper_left_corners = x - 0.5*dx + Rxy_upper_right_corners = x + 0.5*dx + Zxy = y + Zxy_corners = y - 0.5*dy + Zxy_lower_right_corners = y - 0.5*dy + Zxy_upper_left_corners = y + 0.5*dy + Zxy_upper_right_corners = y + 0.5*dy + + [dmplex] + use_cxx_ivertex=true + test_dmplex_cell_volumes = true + test_dmplex_cell_centres = true + + [solver] + type = pvode + + [hermes] + components = (d+, e, vantage) + Nnorm = 1e19 + + [d+] + type = evolve_density + AA = 1 + charge = 1 + + [Nd+] + function = 1 + + [e] + type = quasineutral + AA = 1/1836 + charge = -1 + + [vantage] + nsteps = 0 + test_mass_conservation = true + npart_per_cell = 1 + """ + file = f"{base_grid_dir}/BOUT.inp" + with open(file, "w") as file: + file.write(input_file_string) + # Command to run + cmd = f"./hermes-3 -d {base_grid_dir}" + nproc = 1 + # Launch using MPI, with OMP_NUM_THREADS=1 + if verbose: + print(f"execute: {cmd}") + s, out = launch_safe(cmd, nproc=nproc, mthread=1, pipe=True) + + # copy the file to be used as a grid file, and insert the necessary + # mesh variables that should be computed in preprocessing by the gridding/meshing workflow + cmd = f"cp {base_grid_dir}/BOUT.dmp.vantage.0.nc {kinetic_nc_file_path}" + if verbose: + print(f"execute: {cmd}") + shell(cmd) + + map_RZ_to_itriangle_0 = np.array( + [ + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + [4.0, 6.0, 0.0, 2.0, 4.0, 6.0, 0.0, 2.0], + [12.0, 14.0, 8.0, 10.0, 12.0, 14.0, 8.0, 10.0], + [20.0, 22.0, 16.0, 18.0, 20.0, 22.0, 16.0, 18.0], + [28.0, 30.0, 24.0, 26.0, 28.0, 30.0, 24.0, 26.0], + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + ] + ) + map_RZ_to_itriangle_1 = np.array( + [ + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + [5.0, 7.0, 1.0, 3.0, 5.0, 7.0, 1.0, 3.0], + [13.0, 15.0, 9.0, 11.0, 13.0, 15.0, 9.0, 11.0], + [21.0, 23.0, 17.0, 19.0, 21.0, 23.0, 17.0, 19.0], + [29.0, 31.0, 25.0, 27.0, 29.0, 31.0, 25.0, 27.0], + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + [-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0], + ] + ) + vertices = np.array( + [ + [0.0, 0.0], + [0.0, 1.5707963267948966], + [0.0, 3.141592653589793], + [0.0, 4.71238898038469], + [0.25, 0.0], + [0.25, 1.5707963267948966], + [0.25, 3.141592653589793], + [0.25, 4.71238898038469], + [0.5, 0.0], + [0.5, 1.5707963267948966], + [0.5, 3.141592653589793], + [0.5, 4.71238898038469], + [0.75, 0.0], + [0.75, 1.5707963267948966], + [0.75, 3.141592653589793], + [0.75, 4.71238898038469], + [1.0, 0.0], + [1.0, 1.5707963267948966], + [1.0, 3.141592653589793], + [1.0, 4.71238898038469], + [0.25, 6.283185307179586], + [0.5, 6.283185307179586], + [0.75, 6.283185307179586], + [1.0, 6.283185307179586], + [0.0, 6.283185307179586], + [-1.0, 6.283185307179586], + [-1.0, 0.0], + [2.0, 0.0], + [2.0, 6.283185307179586], + [-0.4999999999986942, 6.283185307179586], + [-1.0, 5.799863360474505], + [-1.0, 5.316541413769544], + [-1.0, 4.833219467065129], + [-1.0, 4.349897520360367], + [-1.0, 3.866575573654959], + [-1.0, 3.383253626950145], + [-1.0, 2.899931680244037], + [-1.0, 2.41660973353672], + [-1.0, 1.933287786829305], + [-1.0, 1.449965840122101], + [-1.0, 0.9666438934146182], + [-1.0, 0.4833219467074823], + [-0.500000000002059, 0.0], + [1.5, 0.0], + [2.0, 0.4833219467051673], + [2.0, 0.9666438934102256], + [2.0, 1.44996584011498], + [2.0, 1.933287786819744], + [2.0, 2.416609733524625], + [2.0, 2.899931680229439], + [2.0, 3.383253626935547], + [2.0, 3.866575573642864], + [2.0, 4.349897520350281], + [2.0, 4.833219467057484], + [2.0, 5.316541413764967], + [2.0, 5.799863360472104], + [1.5, 6.283185307179586], + [-0.6071522086797074, 0.784761011285081], + [-0.607152208680398, 5.498424295897261], + [-0.5271436386483882, 2.235364003519466], + [-0.5271436386490111, 4.047821303667518], + [-0.6377166249249197, 1.29275877701848], + [-0.6377166249255449, 4.990426530165874], + [-0.581430915947327, 4.591558493712746], + [-0.6217149109190194, 2.767018144897421], + [-0.5814309159448909, 1.691626813475703], + [-0.621714910919144, 3.516167162291901], + [-0.6486859643676327, 3.14159265359466], + [-0.6960731294749422, 5.95330761966039], + [-0.3642360307921816, 5.790607309349078], + [-0.696073129475413, 0.3298776875203251], + [-0.3642360307925395, 0.4925779978312037], + [1.607152208679471, 5.498424295894448], + [1.60715220868047, 0.7847610112824089], + [1.527143638648388, 4.047821303660118], + [1.527143638649029, 2.235364003512242], + [1.637716624924872, 4.990426530161094], + [1.63771662492552, 1.292758777013975], + [1.581430915947132, 1.691626813467362], + [1.621714910919019, 3.516167162282163], + [1.581430915944891, 4.59155849370388], + [1.621714910919148, 2.767018144887718], + [1.648685964367633, 3.141592653584932], + [1.696073129475332, 0.3298776875191727], + [1.364236030792596, 0.4925779978304361], + [1.696073129474763, 5.953307619659157], + [1.364236030791846, 5.790607309348226], + ] + ) + tri_cell_vertices = np.array( + [ + [0, 4, 5], + [0, 5, 1], + [1, 5, 6], + [1, 6, 2], + [2, 6, 7], + [2, 7, 3], + [3, 7, 20], + [3, 20, 24], + [4, 8, 9], + [4, 9, 5], + [5, 9, 10], + [5, 10, 6], + [6, 10, 11], + [6, 11, 7], + [7, 11, 21], + [7, 21, 20], + [8, 12, 13], + [8, 13, 9], + [9, 13, 14], + [9, 14, 10], + [10, 14, 15], + [10, 15, 11], + [11, 15, 22], + [11, 22, 21], + [12, 16, 17], + [12, 17, 13], + [13, 17, 18], + [13, 18, 14], + [14, 18, 19], + [14, 19, 15], + [15, 19, 23], + [15, 23, 22], + [1, 59, 65], + [60, 3, 63], + [40, 41, 57], + [30, 31, 58], + [39, 40, 61], + [31, 32, 62], + [29, 25, 68], + [26, 42, 70], + [38, 39, 65], + [32, 33, 63], + [37, 38, 59], + [33, 34, 60], + [36, 37, 64], + [35, 36, 67], + [34, 35, 66], + [41, 26, 70], + [25, 30, 68], + [57, 41, 70], + [30, 58, 68], + [40, 57, 61], + [58, 31, 62], + [2, 66, 67], + [64, 2, 67], + [39, 61, 65], + [62, 32, 63], + [59, 38, 65], + [33, 60, 63], + [3, 62, 63], + [61, 1, 65], + [37, 59, 64], + [60, 34, 66], + [36, 64, 67], + [66, 35, 67], + [57, 70, 71], + [68, 58, 69], + [42, 0, 71], + [24, 29, 69], + [70, 42, 71], + [29, 68, 69], + [0, 1, 71], + [1, 2, 59], + [1, 57, 71], + [59, 2, 64], + [2, 60, 66], + [58, 3, 69], + [57, 1, 61], + [2, 3, 60], + [3, 24, 69], + [3, 58, 62], + [19, 74, 80], + [75, 17, 78], + [54, 55, 72], + [44, 45, 73], + [53, 54, 76], + [45, 46, 77], + [43, 27, 83], + [28, 56, 85], + [52, 53, 80], + [46, 47, 78], + [51, 52, 74], + [47, 48, 75], + [50, 51, 79], + [49, 50, 82], + [48, 49, 81], + [55, 28, 85], + [27, 44, 83], + [72, 55, 85], + [44, 73, 83], + [54, 72, 76], + [73, 45, 77], + [18, 81, 82], + [79, 18, 82], + [53, 76, 80], + [77, 46, 78], + [74, 52, 80], + [47, 75, 78], + [17, 77, 78], + [76, 19, 80], + [51, 74, 79], + [75, 48, 81], + [50, 79, 82], + [81, 49, 82], + [72, 85, 86], + [83, 73, 84], + [16, 43, 84], + [56, 23, 86], + [85, 56, 86], + [43, 83, 84], + [23, 19, 86], + [19, 18, 74], + [19, 72, 86], + [74, 18, 79], + [18, 75, 81], + [73, 17, 84], + [72, 19, 76], + [18, 17, 75], + [17, 16, 84], + [17, 73, 77], + ] + ) + + # open the kinetic nc file and append mesh data + with Dataset(kinetic_nc_file_path, mode="a") as ncdataset: + if verbose: + print(f"append: map_RZ_to_itriangle_0 to {kinetic_nc_file_path}") + ptr_map_RZ_to_itriangle_0 = ncdataset.createVariable( + "map_RZ_to_itriangle_0", "f8", ("x", "y") + ) + ptr_map_RZ_to_itriangle_0[:] = map_RZ_to_itriangle_0 + + if verbose: + print(f"append: map_RZ_to_itriangle_1 to {kinetic_nc_file_path}") + ptr_map_RZ_to_itriangle_1 = ncdataset.createVariable( + "map_RZ_to_itriangle_1", "f8", ("x", "y") + ) + ptr_map_RZ_to_itriangle_1[:] = map_RZ_to_itriangle_1 + + if verbose: + print(f"append: vertices to {kinetic_nc_file_path}") + nvertices, vertexdim = np.shape(vertices) + ncdataset.createDimension("nvertices", nvertices) + ncdataset.createDimension("vertexdim", vertexdim) + ptr_vertices = ncdataset.createVariable( + "vertices", "f8", ("nvertices", "vertexdim") + ) + ptr_vertices[:] = vertices + + if verbose: + print(f"append: tri_cell_vertices to {kinetic_nc_file_path}") + ntriangle, tricorners = np.shape(tri_cell_vertices) + ncdataset.createDimension("ntriangle", ntriangle) + ncdataset.createDimension("tricorners", tricorners) + ptr_tri_cell_vertices = ncdataset.createVariable( + "tri_cell_vertices", "i4", ("ntriangle", "tricorners") + ) + ptr_tri_cell_vertices[:] = tri_cell_vertices + + return None diff --git a/tests/integrated/particle-pusher/runtest b/tests/integrated/particle-pusher/runtest index a754e44fb..e8dadcf11 100755 --- a/tests/integrated/particle-pusher/runtest +++ b/tests/integrated/particle-pusher/runtest @@ -3,23 +3,69 @@ # Python script to run and analyse particle pushing test from boututils.run_wrapper import shell, launch_safe - -import h5py import numpy as np +import h5py import netCDF4 as nc +# tools for extended kinetic mesh test +from extended_kinetic_mesh_tools import ( + extended_kinetic_mesh_test_input, + generate_BOUT_grid_data, + kinetic_mesh_path, +) + verbose = False # Link to the executable shell("ln -s ../../../hermes-3 hermes-3") -# make the run directory -shell("mkdir particle-push-slab-test") + + +def integrate_cellvalue(cellvalues, volumes): + integrand = cellvalues * volumes + integral = sum(integrand) + return integral + + +def get_neutral_mass(vtkhdf_file_path): + # assume that a single .vtkhdf file corresponds to a single time slice + with h5py.File(vtkhdf_file_path, "r") as dataset: + density_cellvalues = np.copy(dataset["/VTKHDF/CellData/density"][:]) + # ion_density_cellvalues = np.copy(dataset[f"/VTKHDF/CellData/ion_density"][:]) + volumes = np.copy(dataset["/VTKHDF/CellData/cellvolume"][:]) + # get the masses from this time + total_neutral_mass = integrate_cellvalue(density_cellvalues, volumes) + return total_neutral_mass + + +# consider how to parallelise this when the test is run in parallel +def get_ion_mass(ncfile_path): + with nc.Dataset(ncfile_path, mode="r") as ncdataset: + ion_density = np.copy(ncdataset.variables["ion_density"][:]) + cellareas = np.copy(ncdataset.variables["neso_cell_areas"][:]) + ion_mass_0 = np.sum(ion_density[0, :, :] * cellareas[:, :]) + ion_mass_end = np.sum(ion_density[-1, :, :] * cellareas[:, :]) + return ion_mass_0, ion_mass_end + + +def get_ionised_mass(ncfile_path): + with nc.Dataset(ncfile_path, mode="r") as ncdataset: + time = np.copy(ncdataset.variables["t_array"][:]) + Siz = np.copy(ncdataset.variables["Siz"][:]) + Srec = np.copy(ncdataset.variables["Srec"][:]) + cellareas = np.copy(ncdataset.variables["neso_cell_areas"][:]) + (ntime,) = time.shape + ionised_mass = 0.0 + for it in range(1, ntime): + dt = time[it] - time[it - 1] + ionised_mass += dt * np.sum(Siz[it, :, :] * cellareas[:, :]) + ionised_mass += dt * np.sum(Srec[it, :, :] * cellareas[:, :]) + return 0, ionised_mass # A function to define the BOUT.inp file contents def particle_push_input( - nx=40, - ny=36, + nx=8, + ny=4, dt=0.005, nsteps=1, iz_rate=0.0, @@ -35,8 +81,8 @@ def particle_push_input( [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 @@ -83,7 +129,8 @@ def particle_push_input( dt = {dt} nsteps = {nsteps} test_mass_conservation = true - initial_neutral_density = 1e19 + initial_neutral_pressure = 1 + initial_neutral_temperature = 1 npart_per_cell = 20 background_ion_density = 1e19 background_ion_temperature = 50 @@ -123,13 +170,12 @@ def load_particle_data(file_path): def test_particle_push( - file, + simdir, mode, - nsteps=10, - remove_threshold=0.0, - merge_threshold=0.0, - iz_rate=0.0, - rec_rate=0.0, + input_file_string, + nsteps, + iz_rate, + rec_rate, ): """ Set up particle push test and pass settings to the input file. @@ -139,30 +185,40 @@ def test_particle_push( - ionisation: ionisation rate > 0, recombination rate = 0. Particle count should increase - recombination: ionisation rate = 0, recombination rate > 0. Particle count should decrease """ + # check correct usage of test function + if mode == "advection": + np.testing.assert_equal( + iz_rate, 0.0, err_msg=f"iz_rate /= 0.0 in test mode={mode}" + ) + np.testing.assert_equal( + rec_rate, 0.0, err_msg=f"rec_rate /= 0.0 in test mode={mode}" + ) + elif mode == "ionisation": + np.testing.assert_equal( + rec_rate, 0.0, err_msg=f"rec_rate /= 0.0 in test mode={mode}" + ) + np.testing.assert_array_less( + 0.0, iz_rate, err_msg=f"iz_rate > 0.0 not satisfied in test mode={mode}" + ) + elif mode == "recombination": + np.testing.assert_equal( + iz_rate, 0.0, err_msg=f"iz_rate /= 0.0 in test mode={mode}" + ) + np.testing.assert_array_less( + 0.0, rec_rate, err_msg=f"rec_rate > 0.0 not satisfied in test mode={mode}" + ) + file = f"{simdir}/BOUT.inp" with open(file, "w") as file: - file.write( - particle_push_input( - nx=40, - ny=36, - dt=0.005, - nsteps=nsteps, - iz_rate=iz_rate, - rec_rate=rec_rate, - remove_threshold=remove_threshold, - merge_threshold=merge_threshold, - ) - ) + file.write(input_file_string) # Command to run - cmd = "./hermes-3 -d particle-push-slab-test" + cmd = f"./hermes-3 -d {simdir}" nproc = 1 # Launch using MPI s, out = launch_safe(cmd, nproc=nproc, mthread=1, pipe=True) - particle_positions_file_path = ( - r"particle-push-slab-test/particle_trajectories.h5part" - ) + particle_positions_file_path = rf"{simdir}/particle_trajectories.h5part" particle_positions = load_particle_data(particle_positions_file_path) nsteps_out = len(particle_positions) @@ -184,52 +240,136 @@ def test_particle_push( ) # check the mass diagnostic - ncfile_path = r"particle-push-slab-test/BOUT.dmp.vantage.0.nc" + vtkhdf_file_path_time_0 = ( + rf"{simdir}/BOUT.dmp.vantage.particle.moments.istep.0.vtkhdf" + ) + vtkhdf_file_path_time_end = ( + rf"{simdir}/BOUT.dmp.vantage.particle.moments.istep.{nsteps}.vtkhdf" + ) + ncfile_path = rf"{simdir}/BOUT.dmp.vantage.0.nc" - with nc.Dataset(ncfile_path, mode="r") as ncdataset: - total_mass = np.copy(ncdataset.variables["total_mass"][:]) - total_ion_mass = np.copy(ncdataset.variables["total_ion_mass"][:]) - total_neutral_mass = np.copy(ncdataset.variables["total_neutral_mass"][:]) + total_neutral_mass_0 = get_neutral_mass(vtkhdf_file_path_time_0) + total_neutral_mass_end = get_neutral_mass(vtkhdf_file_path_time_end) + total_ion_mass_0, total_ion_mass_end = get_ion_mass(ncfile_path) + total_ionised_mass_0, total_ionised_mass_end = get_ionised_mass(ncfile_path) + + total_mass_0 = total_neutral_mass_0 + total_ion_mass_0 + total_mass_end = total_neutral_mass_end + total_ion_mass_end + + neutral_plus_ionised_mass_0 = total_neutral_mass_0 + total_ionised_mass_0 + neutral_plus_ionised_mass_end = total_neutral_mass_end + total_ionised_mass_end # total mass should always be conserved np.testing.assert_allclose( - total_mass[0], - total_mass[-1], + total_mass_0, + total_mass_end, atol=1.0e-11, err_msg=f"Total mass not conserved in {mode} test", ) + np.testing.assert_allclose( + neutral_plus_ionised_mass_0, + neutral_plus_ionised_mass_end, + atol=1.0e-11, + err_msg=f"Neutral + ionised mass not conserved in {mode} test", + ) # if we ionise particles there should be mass exchange between neutral and ion if mode == "ionisation": np.testing.assert_array_less( - total_ion_mass[0], - total_ion_mass[-1], + total_ion_mass_0, + total_ion_mass_end, err_msg="Ion mass did not increase in ionisation test", ) np.testing.assert_array_less( - total_neutral_mass[-1], - total_neutral_mass[0], + total_ionised_mass_0, + total_ionised_mass_end, + err_msg="Ionised mass did not increase in ionisation test", + ) + np.testing.assert_array_less( + total_neutral_mass_end, + total_neutral_mass_0, err_msg="Neutral mass did not decrease in ionisation test", ) # Recombination is the opposite if mode == "recombination": np.testing.assert_array_less( - total_ion_mass[-1], - total_ion_mass[0], + total_ion_mass_end, + total_ion_mass_0, err_msg="Ion mass did not decrease in recombination test", ) np.testing.assert_array_less( - total_neutral_mass[0], - total_neutral_mass[-1], + total_ionised_mass_end, + total_ionised_mass_0, + err_msg="Ionised mass did not decrease in recombination test", + ) + np.testing.assert_array_less( + total_neutral_mass_0, + total_neutral_mass_end, err_msg="Neutral mass did not increase in recombination test", ) + return None + + +# function to test the mesh with triangular cells, aligned with BOUT++ grid, +# and with a finite volume outside the BOUT++ domain +def test_particle_push_extended_kinetic_mesh( + simdir, + mode, + bout_grid_file, + msh_file, + nsteps=10, + remove_threshold=0.0, + merge_threshold=0.0, + iz_rate=0.0, + rec_rate=0.0, +): + input_file_string = extended_kinetic_mesh_test_input( + bout_grid_file, + msh_file, + 0.005, + nsteps, + iz_rate, + rec_rate, + remove_threshold, + merge_threshold, + ) + return test_particle_push( + simdir, mode, input_file_string, nsteps, iz_rate, rec_rate + ) -file = "particle-push-slab-test/BOUT.inp" +# function to test the mesh with quad cells, aligned with BOUT++ grid +def test_particle_push_bout_quad_kinetic_mesh( + simdir, + mode, + nsteps=10, + remove_threshold=0.0, + merge_threshold=0.0, + iz_rate=0.0, + rec_rate=0.0, +): + input_file_string = particle_push_input( + nx=8, + ny=4, + dt=0.005, + nsteps=nsteps, + iz_rate=iz_rate, + rec_rate=rec_rate, + remove_threshold=remove_threshold, + merge_threshold=merge_threshold, + ) + return test_particle_push( + simdir, mode, input_file_string, nsteps, iz_rate, rec_rate + ) -test_particle_push( - file, + +# make the run directory for the test on the BOUT++/quadrilaterals grid +simdir = "particle-push-slab-test" +shell(f"mkdir {simdir}") + +test_particle_push_bout_quad_kinetic_mesh( + simdir, "advection", nsteps=3, remove_threshold=1.0e-10, @@ -240,8 +380,8 @@ test_particle_push( if verbose: print("Advection test passed") -test_particle_push( - file, +test_particle_push_bout_quad_kinetic_mesh( + simdir, "ionisation", nsteps=3, remove_threshold=1.0e-10, @@ -252,8 +392,8 @@ test_particle_push( if verbose: print("Ionisation test passed") -test_particle_push( - file, +test_particle_push_bout_quad_kinetic_mesh( + simdir, "recombination", nsteps=3, remove_threshold=1.0e-10, @@ -264,4 +404,58 @@ test_particle_push( if verbose: print("Recombination test passed") +# now test on a kinetic mesh of triangles, where a subset of the mesh +# overlaps with the BOUT++ "quadrilaterals" +base_grid_dir = "basic_slab" +# path of modified .nc file with extra kinetic mesh data +kinetic_nc_file_path = f"{base_grid_dir}/BOUT.dmp.vantage.0.with_kinetic_mesh.nc" +generate_BOUT_grid_data(base_grid_dir, kinetic_nc_file_path, verbose=verbose) +# path to the pre-computed GMSH .msh file for the kinetic mesh +kinetic_mesh_file_path = kinetic_mesh_path() +# use the generated files to run the test +simdir = "particle-push-extended-kinetic-mesh-slab-test" +shell(f"mkdir {simdir}") +test_particle_push_extended_kinetic_mesh( + simdir, + "advection", + kinetic_nc_file_path, + kinetic_mesh_file_path, + nsteps=3, + remove_threshold=1.0e-10, + merge_threshold=0, + iz_rate=0.00, + rec_rate=0.00, +) +if verbose: + print("Advection test passed -- extended kinetic mesh") + +test_particle_push_extended_kinetic_mesh( + simdir, + "ionisation", + kinetic_nc_file_path, + kinetic_mesh_file_path, + nsteps=3, + remove_threshold=1.0e-10, + merge_threshold=0, + iz_rate=0.1, + rec_rate=0.00, +) +if verbose: + print("Ionisation test passed -- extended kinetic mesh") + +test_particle_push_extended_kinetic_mesh( + simdir, + "recombination", + kinetic_nc_file_path, + kinetic_mesh_file_path, + nsteps=3, + remove_threshold=1.0e-10, + merge_threshold=0, + iz_rate=0.00, + rec_rate=0.1, +) +if verbose: + print("Recombination test passed -- extended kinetic mesh") + + print(" => Test passed") diff --git a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp index 4e77ae287..5a01f69a6 100644 --- a/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp +++ b/tests/integrated/vantage-iz-rec-balance/data/BOUT.inp @@ -56,8 +56,9 @@ dt = 700 nsteps = 3 # This is the final expected steady state density -initial_neutral_density = 8e18 - +# initial_neutral_density = 8e18 +initial_neutral_pressure = 1.2817413072 +initial_neutral_temperature = 1 test_mass_conservation = true