Working code on GPU | tests passing

include/compute/coupled_interfaces.hpp

@@ -147,13 +147,16 @@ struct coupled_interfaces {
 
 namespace iterator {
 
+KOKKOS_FUNCTION
 int npoints(const specfem::enums::coupling::edge::type &edge, const int ngllx,
             const int ngllz);
 
+KOKKOS_FUNCTION
 void self_iterator(const int &ipoint,
                    const specfem::enums::coupling::edge::type &edge,
                    const int ngllx, const int ngllz, int &i, int &j);
 
+KOKKOS_FUNCTION
 void coupled_iterator(const int &ipoint,
                       const specfem::enums::coupling::edge::type &edge,
                       const int ngllx, const int ngllz, int &i, int &j);

include/coupled_interface/coupled_interface.tpp

@@ -70,16 +70,19 @@ template <class self_domain_type, class coupled_domain_type>
 void specfem::coupled_interface::coupled_interface<
     self_domain_type, coupled_domain_type>::compute_coupling() {
 
+  const int nedges = this->edges.extent(0);
+
   Kokkos::parallel_for(
       "specfem::coupled_interfaces::coupled_interfaces::compute_coupling",
-      specfem::kokkos::DeviceTeam(this->edges.extent(0), Kokkos::AUTO, 1),
+      specfem::kokkos::DeviceTeam(nedges, 5, 1),
       KOKKOS_CLASS_LAMBDA(
           const specfem::kokkos::DeviceTeam::member_type &team_member) {
         // Get number of quadrature points
         int ngllx, ngllz;
         quadrature_points.get_ngll(&ngllx, &ngllz);
+        int iedge = team_member.league_rank();
         // Get the edge
-        auto &edge = this->edges(team_member.league_rank());
+        auto edge = this->edges(iedge);
         // Get the edge types
         specfem::enums::coupling::edge::type self_edge_type;
         specfem::enums::coupling::edge::type coupled_edge_type;
@@ -92,9 +95,11 @@ void specfem::coupled_interface::coupled_interface<
         // Iterate over the edges using TeamThreadRange
         Kokkos::parallel_for(
             Kokkos::TeamThreadRange(team_member, npoints),
-            [=](const int &ipoint) { edge.compute_coupling(ipoint); });
+            [=](const int ipoint) { edge.compute_coupling(ipoint); });
       });
 
+  Kokkos::fence();
+
   return;
 }
 

include/coupled_interface/impl/edge/edge.hpp

@@ -21,7 +21,7 @@ struct self_iterator {
                 const int &ngllx, const int &ngllz)
       : edge_type(edge_type), ngllx(ngllx), ngllz(ngllz){};
 
-  KOKKOS_INLINE_FUNCTION
+  KOKKOS_FUNCTION
   void operator()(const int &ipoint, int &i, int &j) const {
     specfem::compute::coupled_interfaces::iterator::self_iterator(
         ipoint, this->edge_type, this->ngllx, this->ngllz, i, j);
@@ -40,7 +40,7 @@ struct coupled_iterator {
                    const int &ngllx, const int &ngllz)
       : edge_type(edge_type), ngllx(ngllx), ngllz(ngllz){};
 
-  KOKKOS_INLINE_FUNCTION
+  KOKKOS_FUNCTION
   void operator()(const int &ipoint, int &i, int &j) const {
     specfem::compute::coupled_interfaces::iterator::coupled_iterator(
         ipoint, this->edge_type, this->ngllx, this->ngllz, i, j);

include/coupled_interface/impl/edge/elastic_acoustic/acoustic_elastic.hpp

@@ -26,6 +26,7 @@ class edge<
   using quadrature_points_type = qp_type;
 
   edge(){};
+
   edge(const int &inum_edge,
        const specfem::domain::domain<specfem::enums::element::medium::acoustic,
                                      qp_type> &self_domain,
@@ -39,9 +40,10 @@ class edge<
        const specfem::kokkos::DeviceView1d<type_real> wzgll,
        const specfem::kokkos::DeviceView3d<int> ibool);
 
+  KOKKOS_FUNCTION
   void compute_coupling(const int &ipoint) const;
 
-  void
+  KOKKOS_FUNCTION void
   get_edges(specfem::enums::coupling::edge::type &self_edge_type,
             specfem::enums::coupling::edge::type &coupled_edge_type) const {
     self_edge_type = this->acoustic_edge;

include/coupled_interface/impl/edge/elastic_acoustic/acoustic_elastic.tpp

@@ -51,9 +51,9 @@ specfem::coupled_interface::impl::edges::edge<
 #endif
 
   this->acoustic_ispec =
-      coupled_interfaces.elastic_acoustic.acoustic_ispec(inum_edge);
+      coupled_interfaces.elastic_acoustic.h_acoustic_ispec(inum_edge);
   this->elastic_ispec =
-      coupled_interfaces.elastic_acoustic.elastic_ispec(inum_edge);
+      coupled_interfaces.elastic_acoustic.h_elastic_ispec(inum_edge);
 
   self_ibool = Kokkos::subview(ibool, acoustic_ispec, Kokkos::ALL, Kokkos::ALL);
   coupled_ibool =
@@ -70,8 +70,8 @@ specfem::coupled_interface::impl::edges::edge<
   jacobian = Kokkos::subview(partial_derivatives.jacobian, acoustic_ispec,
                              Kokkos::ALL, Kokkos::ALL);
 
-  acoustic_edge = coupled_interfaces.elastic_acoustic.acoustic_edge(inum_edge);
-  elastic_edge = coupled_interfaces.elastic_acoustic.elastic_edge(inum_edge);
+  acoustic_edge = coupled_interfaces.elastic_acoustic.h_acoustic_edge(inum_edge);
+  elastic_edge = coupled_interfaces.elastic_acoustic.h_elastic_edge(inum_edge);
 
   self_iterator = specfem::coupled_interface::impl::edges::self_iterator(
       acoustic_edge, ngllx, ngllz);
@@ -91,7 +91,7 @@ specfem::coupled_interface::impl::edges::edge<
 }
 
 template <typename qp_type>
-void specfem::coupled_interface::impl::edges::edge<
+KOKKOS_FUNCTION void specfem::coupled_interface::impl::edges::edge<
     specfem::domain::domain<specfem::enums::element::medium::acoustic, qp_type>,
     specfem::domain::domain<specfem::enums::element::medium::elastic,
                             qp_type> >::compute_coupling(const int &ipoint)

include/coupled_interface/impl/edge/elastic_acoustic/elastic_acoustic.hpp

@@ -42,9 +42,10 @@ class edge<
        const specfem::kokkos::DeviceView1d<type_real> wzgll,
        const specfem::kokkos::DeviceView3d<int> ibool);
 
+  KOKKOS_FUNCTION
   void compute_coupling(const int &ipoint) const;
 
-  void
+  KOKKOS_FUNCTION void
   get_edges(specfem::enums::coupling::edge::type &self_edge_type,
             specfem::enums::coupling::edge::type &coupled_edge_type) const {
     self_edge_type = this->elastic_edge;

include/coupled_interface/impl/edge/elastic_acoustic/elastic_acoustic.tpp

@@ -51,8 +51,8 @@ specfem::coupled_interface::impl::edges::edge<
   assert(wzgll.extent(0) == ngllz);
 #endif
 
-  ispec_acoustic = coupled_interfaces.elastic_acoustic.acoustic_ispec(inum_edge);
-  ispec_elastic = coupled_interfaces.elastic_acoustic.elastic_ispec(inum_edge);
+  ispec_acoustic = coupled_interfaces.elastic_acoustic.h_acoustic_ispec(inum_edge);
+  ispec_elastic = coupled_interfaces.elastic_acoustic.h_elastic_ispec(inum_edge);
 
   self_ibool = Kokkos::subview(ibool, ispec_elastic, Kokkos::ALL, Kokkos::ALL);
   coupled_ibool =
@@ -69,9 +69,9 @@ specfem::coupled_interface::impl::edges::edge<
   jacobian = Kokkos::subview(partial_derivatives.jacobian, ispec_acoustic,
                              Kokkos::ALL, Kokkos::ALL);
 
-  elastic_edge = coupled_interfaces.elastic_acoustic.elastic_edge(inum_edge);
+  elastic_edge = coupled_interfaces.elastic_acoustic.h_elastic_edge(inum_edge);
   acoustic_edge =
-      coupled_interfaces.elastic_acoustic.acoustic_edge(inum_edge);
+      coupled_interfaces.elastic_acoustic.h_acoustic_edge(inum_edge);
 
   self_iterator = specfem::coupled_interface::impl::edges::self_iterator(
       elastic_edge, ngllx, ngllz);
@@ -90,7 +90,7 @@ specfem::coupled_interface::impl::edges::edge<
 }
 
 template <typename qp_type>
-void specfem::coupled_interface::impl::edges::edge<
+KOKKOS_FUNCTION void specfem::coupled_interface::impl::edges::edge<
     specfem::domain::domain<specfem::enums::element::medium::elastic, qp_type>,
     specfem::domain::domain<specfem::enums::element::medium::acoustic,
                             qp_type> >::compute_coupling(const int &ipoint)

src/compute/coupled_interfaces.cpp

@@ -64,6 +64,7 @@ void get_edge_range(const specfem::enums::coupling::edge::type &edge,
 
 // Given an edge, return the number of points along the edge
 // This ends up being important when ngllx != ngllz
+KOKKOS_FUNCTION
 int specfem::compute::coupled_interfaces::iterator::npoints(
     const specfem::enums::coupling::edge::type &edge, const int ngllx,
     const int ngllz) {
@@ -78,10 +79,11 @@ int specfem::compute::coupled_interfaces::iterator::npoints(
     return ngllz;
     break;
   default:
-    throw std::runtime_error("Invalid edge type");
+    assert(false && "Invalid edge type");
   }
 }
 
+KOKKOS_FUNCTION
 void specfem::compute::coupled_interfaces::iterator::self_iterator(
     const int &ipoint, const specfem::enums::coupling::edge::type &edge,
     const int ngllx, const int ngllz, int &i, int &j) {
@@ -104,10 +106,11 @@ void specfem::compute::coupled_interfaces::iterator::self_iterator(
     j = ngllz - 1 - ipoint;
     break;
   default:
-    throw std::runtime_error("Invalid edge type");
+    assert(false && "Invalid edge type");
   }
 }
 
+KOKKOS_FUNCTION
 void specfem::compute::coupled_interfaces::iterator::coupled_iterator(
     const int &ipoint, const specfem::enums::coupling::edge::type &edge,
     const int ngllx, const int ngllz, int &i, int &j) {
@@ -130,7 +133,7 @@ void specfem::compute::coupled_interfaces::iterator::coupled_iterator(
     j = ipoint;
     break;
   default:
-    throw std::runtime_error("Invalid edge type");
+    assert(false && "Invalid edge type");
   }
 }
 

src/source/utilities.cpp

@@ -37,7 +37,7 @@ specfem::forcing_function::stf *assign_stf(std::string forcing_type,
 
     Kokkos::fence();
   } else {
-    throw std::runtime_error("Unknown forcing function type.");
+    assert(false && "Unknown forcing type");
   }
 
   return forcing_function;

src/specfem2d.cpp

@@ -198,11 +198,11 @@ void execute(const std::string &parameter_file, const std::string &default_file,
   // Instantiate coupled interfaces
   specfem::coupled_interface::coupled_interface acoustic_elastic_interface(
       acoustic_domain_static, elastic_domain_static, coupled_interfaces, qp5,
-      partial_derivatives, compute.h_ibool, gllx->get_w(), gllz->get_w());
+      partial_derivatives, compute.ibool, gllx->get_w(), gllz->get_w());
 
   specfem::coupled_interface::coupled_interface elastic_acoustic_interface(
       elastic_domain_static, acoustic_domain_static, coupled_interfaces, qp5,
-      partial_derivatives, compute.h_ibool, gllx->get_w(), gllz->get_w());
+      partial_derivatives, compute.ibool, gllx->get_w(), gllz->get_w());
 
   // Instantiate the writer
   auto writer =

tests/unit-tests/displacement_tests/Newmark/newmark_tests.cpp

@@ -189,12 +189,12 @@ TEST(DISPLACEMENT_TESTS, newmark_scheme_tests) {
       // Instantiate coupled interfaces
       specfem::coupled_interface::coupled_interface acoustic_elastic_interface(
           acoustic_domain_static, elastic_domain_static, coupled_interfaces,
-          qp5, partial_derivatives, compute.h_ibool, gllx->get_w(),
+          qp5, partial_derivatives, compute.ibool, gllx->get_w(),
           gllz->get_w());
 
       specfem::coupled_interface::coupled_interface elastic_acoustic_interface(
           elastic_domain_static, acoustic_domain_static, coupled_interfaces,
-          qp5, partial_derivatives, compute.h_ibool, gllx->get_w(),
+          qp5, partial_derivatives, compute.ibool, gllx->get_w(),
           gllz->get_w());
 
       specfem::solver::solver *solver = new specfem::solver::time_marching<