Extending the kinetic mesh beyond the BOUT++ plasma domain

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)

include/vantage.hxx

@@ -2,65 +2,20 @@
 #include "../include/component.hxx"
 #include "bout/bout.hxx"
 #include "bout/petsclib.hxx"
+#include <bout/bout_types.hxx>
+#include <bout/field2d.hxx>
+#include <memory>
 #include <neso_particles.hpp>
 #include <neso_rng_toolkit.hpp>
 #include <reactions/reactions.hpp>
+#include <vector>
+#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<CellwiseAccumulator<REAL>> accumulator;
-  std::shared_ptr<ParticleGroup> particle_group;
-  std::shared_ptr<TransformationStrategy> 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<PetscInterface::DMPlexInterface>& 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<CellwiseAccumulator<REAL>> accumulator,
-                  std::shared_ptr<ParticleGroup> particle_group,
-                  std::shared_ptr<TransformationStrategy> 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<std::string, VantageSource> sources;
-  std::shared_ptr<PetscInterface::DMPlexInterface> 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,43 +53,58 @@ 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<bout::OptionsIO>
-      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<REAL> neutral_density, total_density;
+  std::vector<REAL> 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<Field2D> 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<REAL> 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<PetscInterface::DMPlexInterface> neso_mesh;
   std::shared_ptr<SYCLTarget> sycl_target;
-  std::shared_ptr<PetscInterface::BoundaryInteraction2D>
-      b2d; // Boundary interaction object
-  std::shared_ptr<PetscInterface::DMPlexProjectEvaluateDG>
-      dg0;                        // DMPlex projection object
+  std::shared_ptr<PetscInterface::BoundaryInteraction2D> b2d;
+  std::shared_ptr<PetscInterface::DMPlexProjectEvaluateDG> project_eval_dg0;
+  std::shared_ptr<PetscInterface::DMPlexMeshCouplerDG0> mesh_coupler_dg0;
+  std::vector<REAL> dof_kinetic_mesh_scalar;
+  std::vector<REAL> dof_bout_mesh_scalar;
+  std::vector<PetscInt> 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> h5part; // HDF5 particle output object
-  std::vector<REAL>
-      h_project1; // Buffer for scalar projection/evaluation of NESO-Particles properties
   std::shared_ptr<ParticleGroup> A_particle_group; // Particle group for main neutrals
   std::shared_ptr<ParticleGroup> marker_group;     // Particle group for rec markers
 
   // Needed for Vantage::apply_boundary_conditions
   std::shared_ptr<BoundaryReflection> 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<VantageSourceManager>
+  std::shared_ptr<VantageDataTransfer> data_transfer; // Manager for VANTAGE data transfer
+  std::unique_ptr<VantageDiagnosticsManager>
+      diagnostics_manager;           // Manager for VANTAGE diagnostics
+  std::shared_ptr<VantageSourceManager>
       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<ReactionController> reaction_controller;

include/vantage_datatransfer.hxx

@@ -0,0 +1,73 @@
+#pragma once
+#include "bout/bout.hxx"
+#include <bout/bout_types.hxx>
+#include <neso_particles.hpp>
+#include <vector>
+
+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<PetscInterface::DMPlexInterface>& neso_mesh,
+    std::shared_ptr<PetscInterface::DMPlexProjectEvaluateDG>& project_eval_dg0,
+    std::shared_ptr<PetscInterface::DMPlexMeshCouplerDG0>& mesh_coupler,
+    Mesh* bout_mesh, size_t ndim_vector);
+
+  // transfer functions for physics scalars
+  void transfer_scalar_to_plasma_grid(
+    std::vector<REAL>& scalar_kinetic_mesh,
+    Field2D& scalar_plasma_grid);
+  void transfer_scalar_to_kinetic_mesh(
+    Field2D& scalar_plasma_grid,
+    std::vector<REAL>& scalar_kinetic_mesh);
+  void transfer_scalar_to_particle_property(
+    std::vector<REAL>& scalar_kinetic_mesh,
+    std::shared_ptr<ParticleGroup>& A_particle_group,
+    std::string particle_property);
+  void transfer_scalar_to_particle_property(
+    Field2D& scalar_plasma_grid,
+    std::shared_ptr<ParticleGroup>& A_particle_group,
+    std::string particle_property);
+  void transfer_particle_property_to_scalar(
+    std::shared_ptr<ParticleGroup>& A_particle_group,
+    std::string particle_property,
+    std::vector<REAL>& scalar_kinetic_mesh);
+  void transfer_particle_property_to_vector(
+    std::shared_ptr<ParticleGroup>& A_particle_group,
+    std::string particle_property,
+    std::vector<REAL>& vector_kinetic_mesh);
+
+  // transfer functions for physics vectors
+  void transfer_vector_to_particle_property(
+    std::vector<Field2D>& vector_plasma_grid,
+    std::shared_ptr<ParticleGroup>& A_particle_group,
+    std::string particle_property);
+  void transfer_vector_to_particle_property(
+    std::vector<REAL>& vector_kinetic_mesh,
+    std::shared_ptr<ParticleGroup>& A_particle_group,
+    std::string particle_property);
+  void transfer_vector_to_kinetic_mesh(
+    std::vector<Field2D>& vector_plasma_grid,
+    std::vector<REAL>& vector_kinetic_mesh);
+
+
+private:
+  // internal variables needed for data transfer
+  std::shared_ptr<PetscInterface::DMPlexInterface> neso_mesh;
+  std::shared_ptr<PetscInterface::DMPlexProjectEvaluateDG> project_eval_dg0;
+  std::shared_ptr<PetscInterface::DMPlexMeshCouplerDG0> 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<REAL> dof_bout_grid_scalar;
+  std::vector<REAL> dof_kinetic_mesh_scalar;
+  // number of physics vector components
+  size_t ndim_vector;
+  std::vector<REAL> dof_bout_grid_vector;
+  std::vector<REAL> dof_kinetic_mesh_vector;
+};

include/vantage_diagnostics.hxx

@@ -0,0 +1,87 @@
+#pragma once
+#include "bout/bout.hxx"
+#include <bout/bout_types.hxx>
+#include <memory>
+#include <neso_particles.hpp>
+#include <string>
+#include <vector>
+#include "../include/vantage_datatransfer.hxx"
+#include "../include/vantage_sources.hxx"
+
+using namespace NESO::Particles;
+
+REAL calculate_total_mass(std::vector<REAL>& density,
+          std::shared_ptr<PetscInterface::DMPlexInterface>& neso_mesh);
+REAL calculate_total_mass(Field2D& density,
+          std::vector<REAL>& neso_cell_volume_on_bout_mesh);
+
+/// @brief  Class to manage diagnostics from VANTAGE.
+class VantageDiagnosticsManager {
+public:
+  VantageDiagnosticsManager(std::string vtkhdf_filename,
+    std::shared_ptr<PetscInterface::DMPlexInterface>& neso_mesh,
+    // neso_mesh cell volumes on the BOUT++ mesh
+    std::vector<REAL>& neso_cell_volumes,
+    std::shared_ptr<ParticleGroup>& A_particle_group,
+    std::shared_ptr<VantageDataTransfer>& data_transfer,
+    std::shared_ptr<VantageSourceManager>& 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<REAL>& 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<REAL> get_density_kinetic_mesh();
+  // Field2D transfer_scalar_to_plasma_grid(std::vector<REAL>& scalar_field);
+
+private:
+  // internal variables needed for diagnostics
+  std::string vtkhdf_filename;
+  std::shared_ptr<PetscInterface::DMPlexInterface> neso_mesh;
+  std::shared_ptr<ParticleGroup> A_particle_group;
+  std::shared_ptr<VantageDataTransfer> data_transfer;
+  // pointer to vantage source_manager for diagnostics
+  std::shared_ptr<VantageSourceManager> 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<bout::OptionsIO>
+      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<REAL> density;
+  std::vector<REAL> energy;
+  std::vector<REAL> gamma;
+  std::vector<REAL> uvector;
+  std::vector<REAL> pressure;
+  std::vector<REAL> temperature;
+
+  // variable used to transfer dofs from kinetic
+  // to BOUT++ meshes
+  std::vector<REAL> 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
+};

include/vantage_dmplex.hxx

@@ -50,10 +50,21 @@ std::vector<PetscInt> 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<SYCLTarget> sycl_target,
-                                std::string dmplex_h5_filename);
+void create_dmplex_from_Bout_mesh(DM* dm, Mesh* bout_mesh, Options& mesh_options,
+                                std::shared_ptr<SYCLTarget> 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<double>& vertices,
+    const std::vector<int>& tri_cell_vertices);
+
+std::vector<double> get_triangle_vertices();
+
+std::vector<int> get_triangle_cell_definition();
+
+#endif
 
 #endif // VANTAGE_DMPLEX_H