removed old style particle loops (#256)

* removed old style particle loops from tests * moved SimpleSOL neutral particles to new style particle loops * moved SimpleSOL mass recording to new style particle loops * moved H3LAPD particle system to new style particle loops * moved H3LAPD mass recorder to new style particle loops * moved Electrostatic2D3V over to new style particle loop * removed obsolete members missed from first pass of mass recording classes * moved particle boundary conditions over to new style particle loop * fixed math functions which are called on device but prefixed with std:: * formatting * fixed bad Sym(Sym(foo)) call
ExCALIBUR-NEPTUNE · Dec 10, 2024 · 618c4ec · 618c4ec
1 parent 0a13b29
commit 618c4ec
Show file tree

Hide file tree

Showing 16 changed files with 414 additions and 1,117 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -118,13 +118,12 @@ set(LIB_SRC_FILES
     ${SRC_DIR}/nektar_interface/particle_cell_mapping/map_particles_host.cpp
     ${SRC_DIR}/nektar_interface/particle_cell_mapping/nektar_graph_local_mapper.cpp
     ${SRC_DIR}/nektar_interface/utilities.cpp
-    ${SRC_DIR}/nektar_interface/solver_base/partsys_base.cpp
-    )
+    ${SRC_DIR}/nektar_interface/solver_base/partsys_base.cpp)
 
-#set(LIB_SRC_FILES_IGNORE ${SRC_DIR}/main.cpp)
+# set(LIB_SRC_FILES_IGNORE ${SRC_DIR}/main.cpp)
 check_file_list(${SRC_DIR} cpp "${LIB_SRC_FILES}" "${LIB_SRC_FILES_IGNORE}")
 
-#set(MAIN_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cpp)
+# set(MAIN_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cpp)
 set(HEADER_FILES
     ${INC_DIR}/nektar_interface/basis_evaluation.hpp
     ${INC_DIR}/nektar_interface/basis_reference.hpp
@@ -218,8 +217,7 @@ set(HEADER_FILES
     ${INC_DIR}/particle_utility/particle_initialisation_line.hpp
     ${INC_DIR}/particle_utility/position_distribution.hpp
     ${INC_DIR}/solvers/solver_callback_handler.hpp
-    ${INC_DIR}/solvers/solver_runner.hpp
-    )
+    ${INC_DIR}/solvers/solver_runner.hpp)
 
 set(HEADER_FILES_IGNORE ${INC_DIR}/revision.hpp)
 check_file_list(${INC_DIR} hpp "${HEADER_FILES}" "${HEADER_FILES_IGNORE}")
@@ -232,8 +230,7 @@ target_include_directories(
   ${NESO_LIBRARY_NAME}
   PUBLIC $<INSTALL_INTERFACE:include>
          $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>)
-target_compile_definitions(${NESO_LIBRARY_NAME} PUBLIC -D${SYCL_FLAG}
-                                                       )
+target_compile_definitions(${NESO_LIBRARY_NAME} PUBLIC -D${SYCL_FLAG})
 target_compile_options(${NESO_LIBRARY_NAME} PRIVATE ${BUILD_TYPE_COMPILE_FLAGS})
 target_link_options(${NESO_LIBRARY_NAME} PUBLIC ${BUILD_TYPE_LINK_FLAGS})
 target_link_libraries(

diff --git a/include/nektar_interface/coordinate_mapping.hpp b/include/nektar_interface/coordinate_mapping.hpp
@@ -33,8 +33,7 @@ template <typename SPECIALISATION> struct BaseCoordinateMapping3D {
   inline void loc_coord_to_loc_collapsed(const T xi0, const T xi1, const T xi2,
                                          T *eta0, T *eta1, T *eta2) {
     auto &underlying = static_cast<SPECIALISATION &>(*this);
-    underlying.loc_coord_to_loc_collapsed_v(xi0, xi1, xi2, eta0, eta1,
-                                                     eta2);
+    underlying.loc_coord_to_loc_collapsed_v(xi0, xi1, xi2, eta0, eta1, eta2);
   }
 
   /**
@@ -54,8 +53,7 @@ template <typename SPECIALISATION> struct BaseCoordinateMapping3D {
   inline void loc_collapsed_to_loc_coord(const T eta0, const T eta1,
                                          const T eta2, T *xi0, T *xi1, T *xi2) {
     auto &underlying = static_cast<SPECIALISATION &>(*this);
-    underlying.loc_collapsed_to_loc_coord_v(eta0, eta1, eta2, xi0, xi1,
-                                                     xi2);
+    underlying.loc_collapsed_to_loc_coord_v(eta0, eta1, eta2, xi0, xi1, xi2);
   }
 
   /**

diff --git a/include/nektar_interface/particle_boundary_conditions.hpp b/include/nektar_interface/particle_boundary_conditions.hpp
@@ -64,21 +64,6 @@ class NektarCartesianPeriodic {
   void execute();
 };
 
-namespace {
-struct NormalType {
-  REAL x;
-  REAL y;
-  REAL z;
-
-  inline NormalType &operator=(const int v) {
-    this->x = v;
-    this->y = v;
-    this->z = v;
-    return *this;
-  }
-};
-} // namespace
-
 /**
  * Implementation of a reflection process which truncates the particle
  * trajectory at the mesh boundary.

diff --git a/solvers/Electrostatic2D3V/Diagnostics/kinetic_energy.hpp b/solvers/Electrostatic2D3V/Diagnostics/kinetic_energy.hpp
@@ -12,8 +12,6 @@ using namespace NESO::Particles;
  */
 class KineticEnergy {
 private:
-  BufferDeviceHost<double> dh_kinetic_energy;
-
 public:
   /// ParticleGroup of interest.
   ParticleGroupSharedPtr particle_group;
@@ -34,69 +32,40 @@ class KineticEnergy {
   KineticEnergy(ParticleGroupSharedPtr particle_group,
                 const double particle_mass, MPI_Comm comm = MPI_COMM_WORLD)
       : particle_group(particle_group), particle_mass(particle_mass),
-        comm(comm), dh_kinetic_energy(particle_group->sycl_target, 1) {
+        comm(comm) {
 
     int flag;
     MPICHK(MPI_Initialized(&flag));
-    ASSERTL1(flag, "MPI is not initialised");
+    NESOASSERT(flag, "MPI is not initialised");
   }
 
   /**
    *  Compute the current kinetic energy of the ParticleGroup.
    */
   inline double compute() {
-
     auto t0 = profile_timestamp();
-    auto sycl_target = this->particle_group->sycl_target;
-    const double k_half_particle_mass = 0.5 * this->particle_mass;
-    auto k_V = (*this->particle_group)[Sym<REAL>("V")]->cell_dat.device_ptr();
-    const auto k_ndim_velocity = (*this->particle_group)[Sym<REAL>("V")]->ncomp;
-
-    const auto pl_iter_range =
-        this->particle_group->mpi_rank_dat->get_particle_loop_iter_range();
-    const auto pl_stride =
-        this->particle_group->mpi_rank_dat->get_particle_loop_cell_stride();
-    const auto pl_npart_cell =
-        this->particle_group->mpi_rank_dat->get_particle_loop_npart_cell();
-
-    this->dh_kinetic_energy.h_buffer.ptr[0] = 0.0;
-    this->dh_kinetic_energy.host_to_device();
-
-    auto k_kinetic_energy = this->dh_kinetic_energy.d_buffer.ptr;
-
-    sycl_target->queue
-        .submit([&](sycl::handler &cgh) {
-          cgh.parallel_for<>(
-              sycl::range<1>(pl_iter_range), [=](sycl::id<1> idx) {
-                NESO_PARTICLES_KERNEL_START
-                const INT cellx = NESO_PARTICLES_KERNEL_CELL;
-                const INT layerx = NESO_PARTICLES_KERNEL_LAYER;
-
-                double half_mvv = 0.0;
-                for (int vdimx = 0; vdimx < k_ndim_velocity; vdimx++) {
-                  const double V_vdimx = k_V[cellx][vdimx][layerx];
-                  half_mvv += (V_vdimx * V_vdimx);
-                }
-                half_mvv *= k_half_particle_mass;
-
-                sycl::atomic_ref<double, sycl::memory_order::relaxed,
-                                 sycl::memory_scope::device>
-                    kinetic_energy_atomic(k_kinetic_energy[0]);
-                kinetic_energy_atomic.fetch_add(half_mvv);
-
-                NESO_PARTICLES_KERNEL_END
-              });
-        })
-        .wait_and_throw();
-    sycl_target->profile_map.inc("KineticEnergy", "Execute", 1,
-                                 profile_elapsed(t0, profile_timestamp()));
-    this->dh_kinetic_energy.device_to_host();
-    const double kernel_kinetic_energy =
-        this->dh_kinetic_energy.h_buffer.ptr[0];
-
-    MPICHK(MPI_Allreduce(&kernel_kinetic_energy, &(this->energy), 1, MPI_DOUBLE,
-                         MPI_SUM, this->comm));
-
+    const REAL k_half_particle_mass = 0.5 * this->particle_mass;
+    const auto k_ndim_velocity =
+        this->particle_group->get_dat(Sym<REAL>("V"))->ncomp;
+
+    auto ga_kinetic_energy = std::make_shared<GlobalArray<REAL>>(
+        this->particle_group->sycl_target, 1, 0.0);
+
+    particle_loop(
+        "KineticEnergy::compute", this->particle_group,
+        [=](auto k_V, auto k_kinetic_energy) {
+          REAL half_mvv = 0.0;
+          for (int vdimx = 0; vdimx < k_ndim_velocity; vdimx++) {
+            const REAL V_vdimx = k_V.at(vdimx);
+            half_mvv += (V_vdimx * V_vdimx);
+          }
+          half_mvv *= k_half_particle_mass;
+          k_kinetic_energy.add(0, half_mvv);
+        },
+        Access::read(Sym<REAL>("V")), Access::add(ga_kinetic_energy))
+        ->execute();
+
+    this->energy = ga_kinetic_energy->get().at(0);
     return this->energy;
   }
 };

diff --git a/solvers/Electrostatic2D3V/Diagnostics/line_field_evaluations.hpp b/solvers/Electrostatic2D3V/Diagnostics/line_field_evaluations.hpp
@@ -205,31 +205,12 @@ template <typename T> class LineFieldEvaluations {
     this->field_evaluate->evaluate(Sym<REAL>("FIELD_EVALUATION"));
 
     if (this->mean_shift) {
-
       const double k_mean = this->field_mean->get_mean();
-
-      const auto pl_iter_range =
-          this->particle_group->mpi_rank_dat->get_particle_loop_iter_range();
-      const auto pl_stride =
-          this->particle_group->mpi_rank_dat->get_particle_loop_cell_stride();
-      const auto pl_npart_cell =
-          this->particle_group->mpi_rank_dat->get_particle_loop_npart_cell();
-
-      const auto k_EVAL = (*this->particle_group)[Sym<REAL>("FIELD_EVALUATION")]
-                              ->cell_dat.device_ptr();
-
-      this->particle_group->sycl_target->queue
-          .submit([&](sycl::handler &cgh) {
-            cgh.parallel_for<>(sycl::range<1>(pl_iter_range),
-                               [=](sycl::id<1> idx) {
-                                 NESO_PARTICLES_KERNEL_START
-                                 const INT cellx = NESO_PARTICLES_KERNEL_CELL;
-                                 const INT layerx = NESO_PARTICLES_KERNEL_LAYER;
-                                 k_EVAL[cellx][0][layerx] -= k_mean;
-                                 NESO_PARTICLES_KERNEL_END
-                               });
-          })
-          .wait_and_throw();
+      particle_loop(
+          "LineFieldEvaluations::write", this->particle_group,
+          [=](auto k_EVAL) { k_EVAL.at(0) -= k_mean; },
+          Access::write(Sym<REAL>("FIELD_EVALUATION")))
+          ->execute();
     }
 
     this->h5part->write(this->step);

diff --git a/solvers/Electrostatic2D3V/Diagnostics/potential_energy.hpp b/solvers/Electrostatic2D3V/Diagnostics/potential_energy.hpp
@@ -23,7 +23,6 @@ template <typename T> class PotentialEnergy {
   Array<OneD, NekDouble> phys_values;
   int num_quad_points;
 
-  BufferDeviceHost<double> dh_energy;
   std::shared_ptr<FieldEvaluate<T>> field_evaluate;
   std::shared_ptr<FieldMean<T>> field_mean;
 
@@ -49,8 +48,7 @@ template <typename T> class PotentialEnergy {
                   ParticleGroupSharedPtr particle_group,
                   std::shared_ptr<CellIDTranslation> cell_id_translation,
                   MPI_Comm comm = MPI_COMM_WORLD)
-      : field(field), particle_group(particle_group), comm(comm),
-        dh_energy(particle_group->sycl_target, 1) {
+      : field(field), particle_group(particle_group), comm(comm) {
 
     int flag;
     MPICHK(MPI_Initialized(&flag));
@@ -79,57 +77,24 @@ template <typename T> class PotentialEnergy {
     this->field_evaluate->evaluate(Sym<REAL>("ELEC_PIC_PE"));
 
     auto t0 = profile_timestamp();
-    auto sycl_target = this->particle_group->sycl_target;
-    const auto k_Q =
-        (*this->particle_group)[Sym<REAL>("Q")]->cell_dat.device_ptr();
-    const auto k_PHI = (*this->particle_group)[Sym<REAL>("ELEC_PIC_PE")]
-                           ->cell_dat.device_ptr();
-
-    const auto pl_iter_range =
-        this->particle_group->mpi_rank_dat->get_particle_loop_iter_range();
-    const auto pl_stride =
-        this->particle_group->mpi_rank_dat->get_particle_loop_cell_stride();
-    const auto pl_npart_cell =
-        this->particle_group->mpi_rank_dat->get_particle_loop_npart_cell();
-
-    this->dh_energy.h_buffer.ptr[0] = 0.0;
-    this->dh_energy.host_to_device();
-
-    auto k_energy = this->dh_energy.d_buffer.ptr;
+    auto ga_energy = std::make_shared<GlobalArray<REAL>>(
+        this->particle_group->sycl_target, 1, 0.0);
     const double k_potential_shift = -this->field_mean->get_mean();
 
-    sycl_target->queue
-        .submit([&](sycl::handler &cgh) {
-          cgh.parallel_for<>(
-              sycl::range<1>(pl_iter_range), [=](sycl::id<1> idx) {
-                NESO_PARTICLES_KERNEL_START
-                const INT cellx = NESO_PARTICLES_KERNEL_CELL;
-                const INT layerx = NESO_PARTICLES_KERNEL_LAYER;
-
-                const double phi = k_PHI[cellx][0][layerx];
-                const double q = k_Q[cellx][0][layerx];
-                const double tmp_contrib = q * (phi + k_potential_shift);
-
-                sycl::atomic_ref<double, sycl::memory_order::relaxed,
-                                 sycl::memory_scope::device>
-                    energy_atomic(k_energy[0]);
-                energy_atomic.fetch_add(tmp_contrib);
-
-                NESO_PARTICLES_KERNEL_END
-              });
-        })
-        .wait_and_throw();
-    sycl_target->profile_map.inc("PotentialEnergy", "Execute", 1,
-                                 profile_elapsed(t0, profile_timestamp()));
-    this->dh_energy.device_to_host();
-    const double kernel_energy = this->dh_energy.h_buffer.ptr[0];
-
-    MPICHK(MPI_Allreduce(&kernel_energy, &(this->energy), 1, MPI_DOUBLE,
-                         MPI_SUM, this->comm));
+    particle_loop(
+        "PotentialEnergy::compute", this->particle_group,
+        [=](auto k_Q, auto k_PHI, auto k_ga_energy) {
+          const REAL phi = k_PHI.at(0);
+          const REAL q = k_Q.at(0);
+          const REAL tmp_contrib = q * (phi + k_potential_shift);
+          k_ga_energy.add(0, tmp_contrib);
+        },
+        Access::read(Sym<REAL>("Q")), Access::read(Sym<REAL>("ELEC_PIC_PE")),
+        Access::add(ga_energy))
+        ->execute();
 
     // The factor of 1/2 in the electrostatic potential energy calculation.
-    this->energy *= 0.5;
-
+    this->energy = ga_energy->get().at(0) * 0.5;
     return this->energy;
   }
 };