Working GPU code | performance issues

CMakeLists.txt

@@ -328,10 +328,6 @@ target_link_libraries(
 
 add_library(
         domain
-        src/domain/impl/elements/acoustic/acoustic2d.cpp
-        src/domain/impl/elements/elastic/elastic2d.cpp
-        src/domain/impl/elements/element.cpp
-        src/domain/impl/elements/kernel.cpp
         src/domain/impl/boundary_conditions/none/none.cpp
         src/domain/impl/boundary_conditions/none/print.cpp
         src/domain/impl/boundary_conditions/dirichlet/dirichlet.cpp
@@ -343,6 +339,10 @@ add_library(
         src/domain/impl/boundary_conditions/composite_stacey_dirichlet/mass_matrix.cpp
         src/domain/impl/boundary_conditions/composite_stacey_dirichlet/print.cpp
         src/domain/impl/boundary_conditions/boundary_conditions.cpp
+        src/domain/impl/elements/acoustic/acoustic2d.cpp
+        src/domain/impl/elements/elastic/elastic2d.cpp
+        src/domain/impl/elements/element.cpp
+        src/domain/impl/elements/kernel.cpp
         src/domain/impl/kernels.cpp
         src/domain/domain.cpp
 )

benchmarks/elastic_domain_forward/rules/generate_mesh.smk

@@ -41,8 +41,7 @@ rule define_par_file:
 rule generate_mesh:
     input:
         par_file = os.path.join(work_directory, "mesh.par"),
-        interfaces = os.path.join(work_directory, "topography.dat"),
-        build_dir = f"{build_dir}"
+        interfaces = os.path.join(work_directory, "topography.dat")
 
     output:
         database = os.path.join(mesh_output_folder, "database.bin"),
@@ -53,5 +52,5 @@ rule generate_mesh:
 
     shell:
         '''
-            {input.build_dir}/xmeshfem2D -p {input.par_file}
+            xmeshfem2D -p {input.par_file}
         '''

benchmarks/elastic_domain_forward/rules/run_simulation.smk

@@ -11,8 +11,7 @@ rule define_specfem_config:
     params:
         mesh_database = os.path.join(mesh_output_folder, "database.bin"),
         source_file = os.path.join(snakefile_dir, "rules/templates/sources.yaml"),
-        stations_file = os.path.join(mesh_output_folder, "STATIONS"),
-        build_dir = os.path.join(work_directory, "build_elastic_domain_forward"),
+        stations_file = os.path.join(mesh_output_folder, "STATIONS")
         work_directory = f"{work_directory}"
 
     localrule:
@@ -34,8 +33,7 @@ rule run_simulation:
     input:
         specfem_config = os.path.join(work_directory, "specfem_config.yaml"),
         mesh_database = os.path.join(mesh_output_folder, "database.bin"),
-        stations = os.path.join(mesh_output_folder, "STATIONS"),
-        build_dir = f"{build_dir}"
+        stations = os.path.join(mesh_output_folder, "STATIONS")
 
     params:
         work_directory = f"{work_directory}"
@@ -56,5 +54,5 @@ rule run_simulation:
             module load boost/1.73.0
             mkdir -p {params.work_directory}/seismograms
             echo "Hostname: $(hostname)" > {output.log}
-            {input.build_dir}/specfem2d -p {input.specfem_config} >> {output.log}
+            specfem2d -p {input.specfem_config} >> {output.log}
         '''

benchmarks/elastic_domain_forward/rules/variables.smk

@@ -1,4 +1,4 @@
 
-snakefile_dir = workflow.basedir
-work_directory = os.path.join(snakefile_dir, "workdir_{nxmax}_{nzmax}_{repeat}")
-build_dir = os.path.join(snakefile_dir, "build_elastic_domain_forward")
+work_directory = os.path.join(workflow.basedir, "workdir_{nxmax}_{nzmax}_{repeat}")
+MESHFEM_BIN = "xmeshfem2D"
+SPECFEM_BIN = "specfem2d"

include/compute/coupled_interfaces/interface_container.tpp

@@ -396,8 +396,8 @@ compute_edge_factors_and_normals(
       throw std::runtime_error("Edge normals need to be opposite in direction");
     }
 
-    const std::array<type_real, 2> weights = { quadrature.gll.weights(i1),
-                                               quadrature.gll.weights(j1) };
+    const std::array<type_real, 2> weights = { quadrature.gll.h_weights(i1),
+                                               quadrature.gll.h_weights(j1) };
 
     edge_factor[ipoint] = [&]() {
       switch (edge1) {

include/compute/fields/data_access.tpp

@@ -1,7 +1,7 @@
 #pragma once
 
-#include "point/coordinates.hpp"
 #include "point/assembly_index.hpp"
+#include "point/coordinates.hpp"
 
 namespace specfem {
 namespace compute {

include/coupled_interface/coupled_interface.tpp

@@ -50,13 +50,12 @@ void specfem::coupled_interface::coupled_interface<
       this->interface_data.get_edge_type();
 
   EdgePolicyType edge_policy(self_index_mapping, coupled_index_mapping,
-                                   self_edge_type, coupled_edge_type,
-                                   this->npoints);
+                             self_edge_type, coupled_edge_type, this->npoints);
 
   Kokkos::parallel_for(
       "specfem::coupled_interfaces::compute_coupling",
       static_cast<typename EdgePolicyType::policy_type &>(edge_policy),
-      KOKKOS_LAMBDA(const typename EdgePolicyType::member_type &team_member) {
+      KOKKOS_CLASS_LAMBDA(const typename EdgePolicyType::member_type &team_member) {
         const auto iterator =
             edge_policy.league_iterator(team_member.league_rank());
 
@@ -74,15 +73,15 @@ void specfem::coupled_interface::coupled_interface<
                          edge_normal(1, iedge, ipoint));
 
               CoupledPointFieldType coupled_field;
-              specfem::compute::load_on_device(coupled_index, field,
+              specfem::compute::load_on_device(coupled_index, this->field,
                                                coupled_field);
 
               SelfPointFieldType acceleration;
               specfem::coupled_interface::impl::compute_coupling(
                   factor, normal, coupled_field, acceleration);
 
               specfem::compute::atomic_add_on_device(self_index, acceleration,
-                                                     field);
+                                                     this->field);
             });
       });
 }

include/domain/impl/boundary_conditions/boundary_conditions.hpp

@@ -5,6 +5,7 @@
 #include "enumerations/boundary.hpp"
 #include "none/none.hpp"
 #include "stacey/stacey.hpp"
+#include <Kokkos_Core.hpp>
 #include <type_traits>
 
 namespace specfem {

include/domain/impl/boundary_conditions/dirichlet/dirichlet.hpp

@@ -23,11 +23,13 @@ KOKKOS_FORCEINLINE_FUNCTION void impl_apply_boundary_conditions(
     const PointPropertyType &, const PointFieldType &,
     PointAccelerationType &acceleration) {
 
+  constexpr static auto tag = PointBoundaryType::boundary_tag;
+
   static_assert(PointBoundaryType::boundary_tag ==
                     specfem::element::boundary_tag::acoustic_free_surface,
                 "Boundary tag must be acoustic_free_surface");
 
-  if (boundary.tag != PointBoundaryType::boundary_tag)
+  if (boundary.tag != tag)
     return;
 
   constexpr int components = PointFieldType::components;

include/domain/impl/boundary_conditions/stacey/stacey.tpp

@@ -43,7 +43,9 @@ impl_enforce_traction(const acoustic_type &, const isotropic_type &,
                     specfem::element::property_tag::isotropic,
                 "Property tag must be isotropic");
 
-  if (boundary.tag != PointBoundaryType::boundary_tag)
+  constexpr static auto tag = PointBoundaryType::boundary_tag;
+
+  if (boundary.tag != tag)
     return;
 
   const auto factor = boundary.edge_weight * boundary.edge_normal.l2_norm();
@@ -120,7 +122,9 @@ impl_enforce_traction(const elastic_type &, const isotropic_type &,
                     specfem::element::property_tag::isotropic,
                 "Property tag must be isotropic");
 
-  if (boundary.tag != PointBoundaryType::boundary_tag)
+  constexpr static auto tag = PointBoundaryType::boundary_tag;
+
+  if (boundary.tag != tag)
     return;
 
   const auto vn =

include/domain/impl/elements/kernel.tpp

@@ -237,7 +237,7 @@ void specfem::domain::impl::kernels::element_kernel_base<
           specfem::algorithms::divergence(
               team, iterator, partial_derivatives, wgll,
               element_quadrature.hprime_wgll, stress_integrand.F,
-              [&](const typename ChunkPolicyType::iterator_type::index_type
+              [&, istep = istep](const typename ChunkPolicyType::iterator_type::index_type
                       &iterator_index,
                   const typename PointAccelerationType::ViewType &result) {
                 auto index = iterator_index.index;

include/enumerations/boundary.hpp

@@ -75,7 +75,7 @@ class boundary_tag_container {
         this->tag = boundary_tag::composite_stacey_dirichlet;
         break;
       default:
-        throw std::runtime_error("Invalid boundary tag");
+        Kokkos::abort("Invalid boundary tag");
         break;
       }
       break;
@@ -91,7 +91,7 @@ class boundary_tag_container {
       case boundary_tag::composite_stacey_dirichlet:
         break;
       default:
-        throw std::runtime_error("Invalid boundary tag");
+        Kokkos::abort("Invalid boundary tag");
         break;
       }
       break;
@@ -105,12 +105,12 @@ class boundary_tag_container {
       case boundary_tag::composite_stacey_dirichlet:
         break;
       default:
-        throw std::runtime_error("Invalid boundary tag");
+        Kokkos::abort("Invalid boundary tag");
         break;
       }
       break;
     default:
-      throw std::runtime_error("Invalid boundary tag");
+      Kokkos::abort("Invalid boundary tag");
       break;
     }
 
@@ -152,7 +152,7 @@ class boundary_tag_container {
       return (this->tag == boundary_tag::composite_stacey_dirichlet);
       break;
     default:
-      throw std::runtime_error("Invalid boundary tag");
+      Kokkos::abort("Invalid boundary tag");
       break;
     }
   }

include/jacobian/jacobian.hpp

@@ -59,8 +59,7 @@ compute_locations(const specfem::kokkos::HostView2d<type_real> s_coorg,
 std::tuple<type_real, type_real>
 compute_locations(const specfem::kokkos::HostTeam::member_type &teamMember,
                   const specfem::kokkos::HostScratchView2d<type_real> s_coorg,
-                  const int ngnod,
-                  const specfem::kokkos::HostView1d<type_real> shape2D);
+                  const int ngnod, const std::vector<type_real> shape2D);
 
 /**
  * @brief Compute global locations (x,z) from shape function matrix calcualted
@@ -73,8 +72,7 @@ compute_locations(const specfem::kokkos::HostTeam::member_type &teamMember,
  */
 std::tuple<type_real, type_real>
 compute_locations(const specfem::kokkos::HostView2d<type_real> s_coorg,
-                  const int ngnod,
-                  const specfem::kokkos::HostView1d<type_real> shape2D);
+                  const int ngnod, const std::vector<type_real> shape2D);
 
 /**
  * @brief Compute partial derivatives at  \f$ (\xi, \gamma) \f$