Compiled code | coupling not tested

include/compute/coupled_interfaces.hpp

@@ -2,6 +2,7 @@
 #define _COMPUTE_COUPLED_INTERFACES_HPP
 
 #include "kokkos_abstractions.h"
+#include "macros.hpp"
 #include "mesh/coupled_interfaces/coupled_interfaces.hpp"
 #include "specfem_enums.hpp"
 
@@ -146,12 +147,26 @@ struct coupled_interfaces {
 
 namespace iterator {
 
+namespace enums {
+
+/**
+ * @brief Tags for edges of the elements on either side of the interface
+ *
+ */
+enum class edge {
+  self,   ///< The edge of the element that is coupled
+  coupled ///< The edge of the element that is coupled to
+};
+} // namespace enums
+
 int get_npoints(const specfem::enums::coupling::edge::type &edge,
                 const int ngllx, const int ngllz);
 
+template <class edge_interface_type>
 void get_points_along_the_edges(
     const int &ipoint, const specfem::enums::coupling::edge::type &edge,
     const int &ngllx, const int &ngllz, int &i, int &j);
+
 } // namespace iterator
 } // namespace coupled_interfaces
 } // namespace compute

include/compute/coupled_interfaces.tpp

@@ -0,0 +1,70 @@
+#ifndef _COMPUTE_COUPLED_INTERFACES_TPP
+#define _COMPUTE_COUPLED_INTERFACES_TPP
+
+#include "compute/coupled_interfaces.hpp"
+
+template <specfem::compute::coupled_interfaces::iterator::enums::edge
+              edge_interface_type>
+void specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges(
+    const int &ipoint, const specfem::enums::coupling::edge::type &edge,
+    const int &ngllx, const int &ngllz, int &i, int &j) {
+
+  bool constexpr self_edge =
+      (edge_interface_type ==
+       specfem::compute::coupled_interfaces::iterator::enums::edge::self);
+
+  bool constexpr coupled_edge =
+      (edge_interface_type ==
+       specfem::compute::coupled_interfaces::iterator::enums::edge::coupled);
+
+  static_assert(!(self_edge && coupled_edge),
+                "Invalid edge type defined for iterator");
+
+  if constexpr (self_edge) {
+    switch (edge) {
+    case specfem::enums::coupling::edge::type::BOTTOM:
+      i = ipoint;
+      j = 0;
+      break;
+    case specfem::enums::coupling::edge::type::TOP:
+      i = ngllx - 1 - ipoint;
+      j = ngllz - 1;
+      break;
+    case specfem::enums::coupling::edge::type::LEFT:
+      i = 0;
+      j = ipoint;
+      break;
+    case specfem::enums::coupling::edge::type::RIGHT:
+      i = ngllx - 1;
+      j = ngllz - 1 - ipoint;
+      break;
+    default:
+      throw std::runtime_error("Invalid edge type");
+    }
+  } else if constexpr (coupled_edge) {
+    switch (edge) {
+    case specfem::enums::coupling::edge::type::BOTTOM:
+      i = ngllx - 1 - ipoint;
+      j = 0;
+      break;
+    case specfem::enums::coupling::edge::type::TOP:
+      i = ipoint;
+      j = ngllz - 1;
+      break;
+    case specfem::enums::coupling::edge::type::LEFT:
+      i = ngllx - 1;
+      j = ngllz - 1 - ipoint;
+      break;
+    case specfem::enums::coupling::edge::type::RIGHT:
+      i = 0;
+      j = ipoint;
+      break;
+    default:
+      throw std::runtime_error("Invalid edge type");
+    }
+  }
+
+  return;
+}
+
+#endif // _COMPUTE_COUPLED_INTERFACES_TPP

include/compute/interface.hpp

@@ -7,5 +7,6 @@
 #include "compute_receivers.hpp"
 #include "compute_sources.hpp"
 #include "coupled_interfaces.hpp"
+#include "coupled_interfaces.tpp"
 
 #endif

include/coupled_interface/coupled_interface.hpp

@@ -31,6 +31,7 @@ class coupled_interface {
 private:
   self_domain_type self_domain;
   coupled_domain_type coupled_domain;
+  quadrature_points_type quadrature_points;
   specfem::kokkos::DeviceView1d<specfem::coupled_interface::impl::edge<
       self_domain_type, coupled_domain_type> >
       edges;

include/coupled_interface/coupled_interface.tpp

@@ -20,7 +20,7 @@ specfem::coupled_interface::
         const specfem::kokkos::DeviceView3d<int> ibool,
         const specfem::kokkos::DeviceView1d<type_real> wxgll,
         const specfem::kokkos::DeviceView1d<type_real> wzgll)
-    : self_domain(self_domain), coupled_domain(coupled_domain) {
+    : self_domain(self_domain), coupled_domain(coupled_domain), quadrature_points(quadrature_points) {
 
   static_assert(std::is_same_v<self_medium, coupled_medium> == false,
                 "Error: self_medium cannot be equal to coupled_medium");
@@ -72,16 +72,22 @@ void specfem::coupled_interface::coupled_interface<
 
   Kokkos::parallel_for(
       "specfem::coupled_interfaces::coupled_interfaces::compute_coupling",
-      specfem::kokkos::DeviceTeam(this->num_edges, Kokkos::AUTO, 1),
+      specfem::kokkos::DeviceTeam(this->edges.extent(0), Kokkos::AUTO, 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);
         // Get the edge
         auto &edge = this->edges(team_member.league_rank());
-        auto &[edge1l, edge2l] = edge.get_edges();
+        // Get the edge types
+        specfem::enums::coupling::edge::type self_edge_type;
+        specfem::enums::coupling::edge::type coupled_edge_type;
+        edge.get_edges(self_edge_type, coupled_edge_type);
         // Get the number of points along the edge
         auto npoints =
             specfem::compute::coupled_interfaces::iterator::get_npoints(
-                edge1l, this->ngllx, this->ngllz);
+                self_edge_type, ngllx, ngllz);
 
         // Iterate over the edges using TeamThreadRange
         Kokkos::parallel_for(

include/coupled_interface/impl/edge/edge.hpp

@@ -4,6 +4,23 @@
 namespace specfem {
 namespace coupled_interface {
 namespace impl {
+
+using self_iterator_type = decltype(
+    specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges<
+        specfem::compute::coupled_interfaces::iterator::enums::edge::self>);
+
+using coupled_iterator_type = decltype(
+    specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges<
+        specfem::compute::coupled_interfaces::iterator::enums::edge::coupled>);
+
+self_iterator_type &self_iterator =
+    specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges<
+        specfem::compute::coupled_interfaces::iterator::enums::edge::self>;
+
+coupled_iterator_type &coupled_iterator =
+    specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges<
+        specfem::compute::coupled_interfaces::iterator::enums::edge::coupled>;
+
 template <class self_domain, class coupled_domain> class edge {};
 } // namespace impl
 } // namespace coupled_interface

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

@@ -2,6 +2,7 @@
 #define _COUPLED_INTERFACE_IMPL_ACOUSTIC_ELASTIC_EDGE_HPP
 
 #include "coupled_interface/impl/edge/edge.hpp"
+#include "domain/interface.hpp"
 #include "kokkos_abstractions.h"
 #include "specfem_enums.hpp"
 #include "specfem_setup.hpp"
@@ -23,7 +24,7 @@ class edge<
                                        qp_type>::medium_type;
   using quadrature_points_type = qp_type;
 
-  edge() = default;
+  edge(){};
   edge(const int &inum_edge,
        const specfem::domain::domain<specfem::enums::element::medium::acoustic,
                                      qp_type> &self_domain,
@@ -37,7 +38,15 @@ class edge<
        const specfem::kokkos::DeviceView1d<type_real> wzgll,
        const specfem::kokkos::DeviceView3d<int> ibool);
 
-  void compute_coupling(const int &ipoint);
+  void compute_coupling(const int &ipoint) const;
+
+  void
+  get_edges(specfem::enums::coupling::edge::type &self_edge_type,
+            specfem::enums::coupling::edge::type &coupled_edge_type) const {
+    self_edge_type = this->self_edge_type;
+    coupled_edge_type = this->coupled_edge_type;
+    return;
+  }
 
 private:
   specfem::kokkos::DeviceView2d<int> self_ibool;

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

@@ -56,13 +56,19 @@ specfem::coupled_interface::impl::edge<
       coupled_interfaces.elastic_acoustic.elastic_ispec(inum_edge);
 
   self_ibool = Kokkos::subview(ibool, self_ispec, Kokkos::ALL, Kokkos::ALL);
-  coupled_ibool = Kokkos::subview(ibool, coupled_ispec, Kokkos::ALL, Kokkos::ALL);
-
-  xix = Kokkos::subview(partial_derivatives.xix, self_ispec, Kokkos::ALL, Kokkos::ALL);
-  xiz = Kokkos::subview(partial_derivatives.xiz, self_ispec, Kokkos::ALL, Kokkos::ALL);
-  gammax = Kokkos::subview(partial_derivatives.gammax, self_ispec, Kokkos::ALL, Kokkos::ALL);
-  gammaz = Kokkos::subview(partial_derivatives.gammaz, self_ispec, Kokkos::ALL, Kokkos::ALL);
-  jacobian = Kokkos::subview(partial_derivatives.jacobian, self_ispec, Kokkos::ALL, Kokkos::ALL);
+  coupled_ibool =
+      Kokkos::subview(ibool, coupled_ispec, Kokkos::ALL, Kokkos::ALL);
+
+  xix = Kokkos::subview(partial_derivatives.xix, self_ispec, Kokkos::ALL,
+                        Kokkos::ALL);
+  xiz = Kokkos::subview(partial_derivatives.xiz, self_ispec, Kokkos::ALL,
+                        Kokkos::ALL);
+  gammax = Kokkos::subview(partial_derivatives.gammax, self_ispec, Kokkos::ALL,
+                           Kokkos::ALL);
+  gammaz = Kokkos::subview(partial_derivatives.gammaz, self_ispec, Kokkos::ALL,
+                           Kokkos::ALL);
+  jacobian = Kokkos::subview(partial_derivatives.jacobian, self_ispec,
+                             Kokkos::ALL, Kokkos::ALL);
 
   self_edge_type = coupled_interfaces.elastic_acoustic.acoustic_edge(inum_edge);
   coupled_edge_type =
@@ -83,21 +89,20 @@ template <typename qp_type>
 void specfem::coupled_interface::impl::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) {
+                            qp_type> >::compute_coupling(const int &ipoint)
+    const {
 
-  const int ngllx, ngllz;
+  int ngllx, ngllz;
   quadrature_points.get_ngll(&ngllx, &ngllz);
 
   int i, j;
-  specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges(
-      ipoint, coupled_edge_type, ngllx, ngllz, i, j);
+  coupled_iterator(ipoint, coupled_edge_type, ngllx, ngllz, i, j);
 
-  const int iglob = coupled_ibool(j, i);
+  int iglob = coupled_ibool(j, i);
   const type_real displ_x = coupled_field(iglob, 0);
   const type_real displ_z = coupled_field(iglob, 1);
 
-  specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges(
-      ipoint, self_edge_type, ngllx, ngllz, i, j);
+  self_iterator(ipoint, self_edge_type, ngllx, ngllz, i, j);
 
   iglob = self_ibool(j, i);
 

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

@@ -26,7 +26,7 @@ class edge<
       specfem::enums::element::medium::acoustic, qp_type>::medium_type;
   using quadrature_points_type = qp_type;
 
-  edge() = default;
+  edge(){};
 
   edge(const int &inum_edge,
        const specfem::domain::domain<specfem::enums::element::medium::elastic,
@@ -41,12 +41,14 @@ class edge<
        const specfem::kokkos::DeviceView1d<type_real> wzgll,
        const specfem::kokkos::DeviceView3d<int> ibool);
 
-  void compute_coupling(const int &ipoint);
+  void compute_coupling(const int &ipoint) const;
 
-  std::pair<specfem::enums::coupling::edge::type,
-            specfem::enums::coupling::edge::type>
-  get_edges() const {
-    return std::make_pair(self_edge_type, coupled_edge_type);
+  void
+  get_edges(specfem::enums::coupling::edge::type &self_edge_type,
+            specfem::enums::coupling::edge::type &coupled_edge_type) const {
+    self_edge_type = this->self_edge_type;
+    coupled_edge_type = this->coupled_edge_type;
+    return;
   }
 
 private: