Merge pull request #77 from ExCALIBUR-NEPTUNE/feature/dmplex-interface

Feature/dmplex interface

CMakeLists.txt

@@ -17,6 +17,7 @@ option(
   ON)
 option(NESO_PARTICLES_LEGACY_DEVICE_SELECTORS
        "Option to use SYCL 1.2 selectors if 2020 ones are not supported." OFF)
+option(NESO_PARTICLES_ENABLE_PETSC "Look for PETSc." OFF)
 option(
   NESO_PARTICLES_DEVICE_AWARE_MPI
   "Enable device aware MPI by default, i.e. pass device pointers to MPI rather than copying to the host."
@@ -45,6 +46,10 @@ set(HEADER_FILES
     ${INCLUDE_DIR_NESO_PARTICLES}/cell_dat_move.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/cell_dat_move_impl.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/communication.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/communication/comm_pair.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/communication/communication_edges_counter.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/communication/communication_typedefs.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/communication/communication_utility.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/compute_target.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/blocked_binary_tree.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/cell_dat.hpp
@@ -59,15 +64,43 @@ set(HEADER_FILES
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/descendant_products.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/global_array.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/local_array.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/containers/lookup_table.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/nd_local_array.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/product_matrix.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/sym_vector.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/sym_vector_impl.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/containers/tuple.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/departing_particle_identification.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/departing_particle_identification_impl.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/device_functions.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/domain.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/error_propagate.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/common/bounding_box.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/common/common.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/common/coordinate_mapping.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/common/local_claim.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/common/overlay_cartesian_mesh.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/common/dof_mapper_dg.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/common/quadrature_point_mapper.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/boundary_interaction/boundary_reflection.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/boundary_interaction/boundary_interaction.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/boundary_interaction/boundary_interaction_2d.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/boundary_interaction/boundary_interaction_common.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/dmplex_cell_serialise.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/dmplex_helper.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/dmplex_interface.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/dmplex_local_mapper.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/dmplex_project_evaluate.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/dmplex_utility.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/particle_cell_mapping/dmplex_2d_mapper.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/particle_cell_mapping/dmplex_host_mapper.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/petsc_common.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/petsc_utility.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/petsc_interface.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/project_evaluate/dmplex_project_evaluate_base.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/project_evaluate/dmplex_project_evaluate_barycentric.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/petsc/project_evaluate/dmplex_project_evaluate_dg.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/external_interfaces/vtk/vtk.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/global_mapping.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/global_mapping_impl.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/global_move.hpp
@@ -81,6 +114,10 @@ set(HEADER_FILES
     ${INCLUDE_DIR_NESO_PARTICLES}/loop/particle_loop_index.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/loop/pli_particle_dat.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/mesh_hierarchy.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/mesh_hierarchy_data/mesh_hierarchy_container.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/mesh_hierarchy_data/mesh_hierarchy_data.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/mesh_hierarchy_data/serial_container.hpp
+    ${INCLUDE_DIR_NESO_PARTICLES}/mesh_hierarchy_data/serial_interface.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/mesh_interface.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/mesh_interface_local_decomp.hpp
     ${INCLUDE_DIR_NESO_PARTICLES}/packing_unpacking.hpp
@@ -168,6 +205,34 @@ if(NESO_PARTICLES_ENABLE_HDF5)
   endif()
 endif()
 
+# set the restrict keyword
+target_compile_definitions(NESO-Particles INTERFACE RESTRICT=${RESTRICT})
+
+# now look for PETSc
+if(NESO_PARTICLES_ENABLE_PETSC)
+  find_package(PkgConfig)
+  if(PKG_CONFIG_FOUND)
+    pkg_check_modules(PETSC PETSc)
+  endif()
+
+  if(PETSC_FOUND)
+    message(STATUS "PETSc found")
+    # TODO CHECK AND CLEANUP
+    target_compile_definitions(NESO-Particles INTERFACE NESO_PARTICLES_PETSC
+                                                        ${PETSC_DEFINITIONS})
+    add_definitions(-DNESO_PARTICLES_PETSC ${PETSC_DEFINITIONS})
+    target_link_libraries(NESO-Particles INTERFACE ${PETSC_LINK_LIBRARIES})
+    target_include_directories(NESO-Particles INTERFACE ${PETSC_INCLUDE_DIRS})
+    message(STATUS PETSC_LINK_LIBRARIES ${PETSC_LINK_LIBRARIES})
+    message(STATUS PETSC_LDFLAGS ${PETSC_LDFLAGS})
+    message(STATUS PETSC_INCLUDE_DIRS ${PETSC_INCLUDE_DIRS})
+    message(STATUS PETSC_LIBRARY_DIRS ${PETSC_LIBRARY_DIRS})
+    message(STATUS PETSC_CFLAGS ${PETSC_CFLAGS})
+  else()
+    message(STATUS "PETSc NOT found")
+  endif()
+endif()
+
 # Find SYCL
 include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/SYCL.cmake)
 if(NESO_PARTICLES_ENABLE_FIND_SYCL)

cmake/restrict-keyword.cmake

@@ -37,4 +37,3 @@ if(NOT DEFINED RESTRICT)
 endif()
 
 message(STATUS "Using restrict keyword: " ${RESTRICT})
-add_compile_definitions(RESTRICT=${RESTRICT})

include/neso_particles.hpp

@@ -13,20 +13,26 @@
 #include "neso_particles/containers/descendant_products.hpp"
 #include "neso_particles/containers/global_array.hpp"
 #include "neso_particles/containers/local_array.hpp"
+#include "neso_particles/containers/lookup_table.hpp"
 #include "neso_particles/containers/nd_index.hpp"
 #include "neso_particles/containers/nd_local_array.hpp"
 #include "neso_particles/containers/product_matrix.hpp"
 #include "neso_particles/containers/rng/rng.hpp"
 #include "neso_particles/containers/sym_vector.hpp"
 #include "neso_particles/containers/sym_vector_impl.hpp"
 #include "neso_particles/containers/tuple.hpp"
+#include "neso_particles/device_functions.hpp"
 #include "neso_particles/domain.hpp"
 #include "neso_particles/error_propagate.hpp"
+#include "neso_particles/external_interfaces/common/common.hpp"
+#include "neso_particles/external_interfaces/petsc/petsc_interface.hpp"
+#include "neso_particles/external_interfaces/vtk/vtk.hpp"
 #include "neso_particles/global_mapping.hpp"
 #include "neso_particles/local_mapping.hpp"
 #include "neso_particles/local_move.hpp"
 #include "neso_particles/loop/particle_loop.hpp"
 #include "neso_particles/mesh_hierarchy.hpp"
+#include "neso_particles/mesh_hierarchy_data/mesh_hierarchy_data.hpp"
 #include "neso_particles/mesh_interface.hpp"
 #include "neso_particles/mesh_interface_local_decomp.hpp"
 #include "neso_particles/parallel_initialisation.hpp"

include/neso_particles/cartesian_mesh.hpp

@@ -4,6 +4,7 @@
 #include <memory>
 #include <mpi.h>
 
+#include "cell_binning.hpp"
 #include "compute_target.hpp"
 #include "local_mapping.hpp"
 #include "loop/particle_loop.hpp"
@@ -37,6 +38,8 @@ class CartesianHMeshLocalMapperT : public LocalMapper {
   std::unique_ptr<BufferDeviceHost<int>> dh_dims;
   // Cell counts of the underlying mesh.
   std::unique_ptr<BufferDeviceHost<int>> dh_cell_counts;
+  // Cell binning for cell move
+  std::unique_ptr<CartesianCellBin> cart_cell_bin;
 
 public:
   /// Disable (implicit) copies.
@@ -224,6 +227,8 @@ class CartesianHMeshLocalMapperT : public LocalMapper {
         }
       }
     }
+    this->cart_cell_bin =
+        std::make_unique<CartesianCellBin>(this->sycl_target, this->mesh);
   };
 
   /**
@@ -232,7 +237,8 @@ class CartesianHMeshLocalMapperT : public LocalMapper {
    *
    *  @param particle_group ParticleGroup to use.
    */
-  inline void map(ParticleGroup &particle_group, const int map_cell = -1) {
+  inline void map(ParticleGroup &particle_group,
+                  const int map_cell = -1) override {
     auto r = ProfileRegion("CartesianHMeshLocalMapper", "map");
 
     ParticleDatSharedPtr<REAL> &position_dat = particle_group.position_dat;
@@ -284,12 +290,24 @@ class CartesianHMeshLocalMapperT : public LocalMapper {
     this->sycl_target->profile_map.add_region(r);
   };
 
+  /**
+   *  Map positions to owning cells. Positions should be within the domain
+   *  prior to calling map, i.e. particles should be within the domain extents.
+   *
+   *  @param particle_group ParticleGroup to use.
+   *  @param map_cell Cell to map.
+   */
+  inline void map_cells(ParticleGroup &particle_group,
+                        const int map_cell = -1) override {
+    this->cart_cell_bin->map_cells(particle_group, map_cell);
+  }
+
   /**
    *  No-op implementation of callback.
    *
    *  @param particle_group ParticleGroup.
    */
-  inline void particle_group_callback(ParticleGroup &particle_group){};
+  inline void particle_group_callback(ParticleGroup &particle_group) override{};
 };
 
 inline std::shared_ptr<CartesianHMeshLocalMapperT>

include/neso_particles/cell_binning.hpp

@@ -16,7 +16,7 @@ namespace NESO::Particles {
  * Bin particle positions into the cells of a CartesianHMesh.
  */
 class CartesianCellBin {
-private:
+protected:
   BufferDevice<int> d_cell_counts;
   BufferDevice<int> d_cell_starts;
   BufferDevice<int> d_cell_ends;
@@ -25,7 +25,43 @@ class CartesianCellBin {
   CartesianHMeshSharedPtr mesh;
   ParticleDatSharedPtr<REAL> position_dat;
   ParticleDatSharedPtr<INT> cell_id_dat;
-  ParticleLoopSharedPtr loop;
+
+  inline ParticleLoopSharedPtr get_loop(ParticleDatSharedPtr<REAL> position_dat,
+                                        ParticleDatSharedPtr<INT> cell_id_dat) {
+    const int k_ndim = this->mesh->ndim;
+    NESOASSERT(((k_ndim > 0) && (k_ndim < 4)), "Bad number of dimensions");
+    auto k_inverse_cell_width_fine = this->mesh->inverse_cell_width_fine;
+    auto k_cell_width_fine = this->mesh->cell_width_fine;
+
+    auto k_cell_counts = this->d_cell_counts.ptr;
+    auto k_cell_starts = this->d_cell_starts.ptr;
+    auto k_cell_ends = this->d_cell_ends.ptr;
+
+    return particle_loop(
+        "CartesianCellBin", position_dat,
+        [=](auto positions, auto cell_id) {
+          int cell_tmps[3];
+
+          for (int dimx = 0; dimx < k_ndim; dimx++) {
+            const REAL pos =
+                positions[dimx] - k_cell_starts[dimx] * k_cell_width_fine;
+            int cell_tmp = ((REAL)pos * k_inverse_cell_width_fine);
+            cell_tmp = (cell_tmp < 0) ? 0 : cell_tmp;
+            cell_tmp = (cell_tmp >= k_cell_ends[dimx]) ? k_cell_ends[dimx] - 1
+                                                       : cell_tmp;
+            cell_tmps[dimx] = cell_tmp;
+          }
+
+          // convert to linear index
+          int linear_index = cell_tmps[k_ndim - 1];
+          for (int dimx = k_ndim - 2; dimx >= 0; dimx--) {
+            linear_index *= k_cell_counts[dimx];
+            linear_index += cell_tmps[dimx];
+          }
+          cell_id[0] = linear_index;
+        },
+        Access::read(position_dat), Access::write(cell_id_dat));
+  }
 
 public:
   /// Disable (implicit) copies.
@@ -45,11 +81,9 @@ class CartesianCellBin {
    * @param cell_id_dat ParticleDat to write particle cell ids to.
    */
   CartesianCellBin(SYCLTargetSharedPtr sycl_target,
-                   CartesianHMeshSharedPtr mesh,
-                   ParticleDatSharedPtr<REAL> position_dat,
-                   ParticleDatSharedPtr<INT> cell_id_dat)
-      : sycl_target(sycl_target), mesh(mesh), position_dat(position_dat),
-        cell_id_dat(cell_id_dat), d_cell_counts(sycl_target, 3),
+                   CartesianHMeshSharedPtr mesh)
+      : sycl_target(sycl_target), mesh(mesh), position_dat(nullptr),
+        cell_id_dat(nullptr), d_cell_counts(sycl_target, 3),
         d_cell_starts(sycl_target, 3), d_cell_ends(sycl_target, 3) {
 
     NESOASSERT(mesh->ndim <= 3, "bad mesh ndim");
@@ -69,53 +103,54 @@ class CartesianCellBin {
         .memcpy(this->d_cell_ends.ptr, mesh->cell_ends,
                 mesh->ndim * sizeof(int))
         .wait_and_throw();
+  }
 
-    const int k_ndim = this->mesh->ndim;
-
-    NESOASSERT(((k_ndim > 0) && (k_ndim < 4)), "Bad number of dimensions");
-    auto k_inverse_cell_width_fine = this->mesh->inverse_cell_width_fine;
-    auto k_cell_width_fine = this->mesh->cell_width_fine;
-
-    auto k_cell_counts = this->d_cell_counts.ptr;
-    auto k_cell_starts = this->d_cell_starts.ptr;
-    auto k_cell_ends = this->d_cell_ends.ptr;
-
-    this->loop = particle_loop(
-        "CartesianCellBin", position_dat,
-        [=](auto positions, auto cell_id) {
-          int cell_tmps[3];
-
-          for (int dimx = 0; dimx < k_ndim; dimx++) {
-            const REAL pos =
-                positions[dimx] - k_cell_starts[dimx] * k_cell_width_fine;
-            int cell_tmp = ((REAL)pos * k_inverse_cell_width_fine);
-            cell_tmp = (cell_tmp < 0) ? 0 : cell_tmp;
-            cell_tmp = (cell_tmp >= k_cell_ends[dimx]) ? k_cell_ends[dimx] - 1
-                                                       : cell_tmp;
-            cell_tmps[dimx] = cell_tmp;
-          }
+  /**
+   * Create instance to bin particles into cells of a CartesianHMesh.
+   *
+   * @param sycl_target SYCLTargetSharedPtr to use as compute device.
+   * @param mesh CartesianHMeshSharedPtr to containing the particles.
+   * @param position_dat ParticleDat with components equal to the mesh dimension
+   * containing particle positions.
+   * @param cell_id_dat ParticleDat to write particle cell ids to.
+   */
+  CartesianCellBin(SYCLTargetSharedPtr sycl_target,
+                   CartesianHMeshSharedPtr mesh,
+                   ParticleDatSharedPtr<REAL> position_dat,
+                   ParticleDatSharedPtr<INT> cell_id_dat)
+      : CartesianCellBin(sycl_target, mesh) {
 
-          // convert to linear index
-          int linear_index = cell_tmps[k_ndim - 1];
-          for (int dimx = k_ndim - 2; dimx >= 0; dimx--) {
-            linear_index *= k_cell_counts[dimx];
-            linear_index += cell_tmps[dimx];
-          }
-          cell_id[0] = linear_index;
-        },
-        Access::read(position_dat), Access::write(cell_id_dat));
+    this->position_dat = position_dat;
+    this->cell_id_dat = cell_id_dat;
   };
 
   /**
    *  Apply the cell binning kernel to each particle stored on this MPI rank.
    *  Particles must be within the domain region owned by this MPI rank.
    */
   inline void execute() {
+    NESOASSERT(this->cell_id_dat != nullptr,
+               "Cell ID Dat not set. Maybe wrong constructor?");
+    NESOASSERT(this->position_dat != nullptr,
+               "Position Dat not set. Maybe wrong constructor?");
     auto r = ProfileRegion("CartesianCellBin", "execute");
-    this->loop->execute();
+    auto loop = this->get_loop(position_dat, cell_id_dat);
+    loop->execute();
     r.end();
     this->sycl_target->profile_map.add_region(r);
   }
+
+  /**
+   * Map call for LocalMapper.
+   *
+   *  @param particle_group ParticleGroup to use.
+   *  @param map_cell Cell to map.
+   */
+  inline void map_cells(ParticleGroup &particle_group,
+                        const int map_cell = -1) {
+    this->get_loop(particle_group.position_dat, particle_group.cell_id_dat)
+        ->execute();
+  }
 };
 
 } // namespace NESO::Particles

include/neso_particles/cell_dat_move.hpp

@@ -17,6 +17,7 @@
 #include "profiling.hpp"
 #include "sycl_typedefs.hpp"
 #include "typedefs.hpp"
+#include <iomanip>
 
 namespace NESO::Particles {
 
@@ -120,6 +121,27 @@ class CellMove {
     event_stack.wait();
   }
 
+  inline void print_particle(const int cell, const int layer) {
+    nprint("Particle info, cell:", cell, "layer:", layer);
+    std::cout << std::setprecision(18);
+    auto lambda_print_dat = [&](auto sym, auto dat) {
+      std::cout << "\t" << sym.name << ": ";
+      auto data = dat->cell_dat.get_cell(cell);
+      auto ncomp = dat->ncomp;
+      for (int cx = 0; cx < ncomp; cx++) {
+        std::cout << data->at(layer, cx) << " ";
+      }
+      std::cout << std::endl;
+    };
+
+    for (auto d : this->particle_dats_int) {
+      lambda_print_dat(d.first, d.second);
+    }
+    for (auto d : this->particle_dats_real) {
+      lambda_print_dat(d.first, d.second);
+    }
+  }
+
 public:
   /// Disable (implicit) copies.
   CellMove(const CellMove &st) = delete;