Implemented coupled interfaces compute forces | compiled

- Implemetned coupled_interfaces class. The class is templated on self_domain_type  and coupled_domain_type  to define the domains along the interface. The class contains a single method which computes the coupling contribution and adds it to the acceleration field of self domain. The method primarily defines the kokkos parallelism and then calls an edge class which computes the coupling term for every GLL point on that interface.
- The edges are templated again on self_domain_type  and coupled_domain_type  and provides specialized implementation for various combinations of the domain types.

include/coupled_interface/coupled_interface.hpp

@@ -0,0 +1,40 @@
+#ifndef _COUPLED_INTERFACE_HPP_
+#define _COUPLED_INTERFACE_HPP_
+
+#include "compute/interface.hpp"
+#include "impl/edge/interface.hpp"
+#include "kokkos_abstractions.h"
+#include "specfem_enums.hpp"
+#include "specfem_setup.hpp"
+
+namespace specfem {
+namespace coupled_interface {
+template <class self_domain_type, class coupled_domain_type>
+class coupled_interface {
+public:
+  using self_medium = typename self_domain_type::medium_type;
+  using coupled_medium = typename coupled_domain_type::medium_type;
+  using quadrature_points_type =
+      typename self_domain_type::quadrature_points_type;
+
+  coupled_interface(
+      self_domain_type &self_domain, coupled_domain_type &coupled_domain,
+      const specfem::compute::coupled_interfaces::coupled_interfaces
+          &coupled_interfaces,
+      const quadrature_points_type &quadrature_points,
+      const specfem::compute::partial_derivatives &partial_derivatives,
+      const specfem::kokkos::DeviceView3d<int> ibool,
+      const specfem::kokkos::DeviceView1d<type_real> wxgll,
+      const specfem::kokkos::DeviceView1d<type_real> wzgll);
+  void compute_coupling();
+
+private:
+  self_domain_type self_domain;
+  coupled_domain_type coupled_domain;
+  specfem::kokkos::DeviceView1d<specfem::coupled_interface::impl::edge<
+      self_domain_type, coupled_domain_type> >
+      edges;
+};
+} // namespace coupled_interface
+} // namespace specfem
+#endif // _COUPLED_INTERFACES_HPP_

include/coupled_interface/coupled_interface.tpp

@@ -0,0 +1,95 @@
+#ifndef _COUPLED_INTERFACE_TPP
+#define _COUPLED_INTERFACE_TPP
+
+#include "compute/interface.hpp"
+#include "coupled_interface.hpp"
+#include "impl/edge/interface.hpp"
+#include "kokkos_abstractions.h"
+#include "macros.hpp"
+#include "specfem_enums.hpp"
+#include <Kokkos_Core.hpp>
+
+template <class self_domain_type, class coupled_domain_type>
+specfem::coupled_interface::
+    coupled_interface<self_domain_type, coupled_domain_type>::coupled_interface(
+        self_domain_type &self_domain, coupled_domain_type &coupled_domain,
+        const specfem::compute::coupled_interfaces::coupled_interfaces
+            &coupled_interfaces,
+        const quadrature_points_type &quadrature_points,
+        const specfem::compute::partial_derivatives &partial_derivatives,
+        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) {
+
+  static_assert(std::is_same_v<self_medium, coupled_medium> == false,
+                "Error: self_medium cannot be equal to coupled_medium");
+
+  bool constexpr elastic_acoustic_condition =
+      (std::is_same_v<self_medium, specfem::enums::element::medium::elastic> &&
+       std::is_same_v<coupled_medium,
+                      specfem::enums::element::medium::acoustic>) ||
+      (std::is_same_v<self_medium, specfem::enums::element::medium::acoustic> &&
+       std::is_same_v<coupled_medium,
+                      specfem::enums::element::medium::elastic>);
+
+  static_assert(elastic_acoustic_condition,
+                "Only acoustic-elastic coupling is supported at the moment.");
+
+  int num_interfaces;
+
+  // Iterate over all domains
+  if (elastic_acoustic_condition) {
+    // Get the number of edges
+    num_interfaces = coupled_interfaces.elastic_acoustic.num_interfaces;
+  }
+
+  // Allocate the edges
+  edges = specfem::kokkos::DeviceView1d<specfem::coupled_interface::impl::edge<
+      self_domain_type, coupled_domain_type> >(
+      "specfem::coupled_interfaces::coupled_interfaces::edges", num_interfaces);
+
+  auto h_edges = Kokkos::create_mirror_view(edges);
+
+  // Iterate over the edges
+  for (int inum_edge = 0; inum_edge < num_interfaces; ++inum_edge) {
+    // Create the edge
+    h_edges(inum_edge) =
+        specfem::coupled_interface::impl::edge<self_domain_type,
+                                               coupled_domain_type>(
+            inum_edge, self_domain, coupled_domain, quadrature_points,
+            coupled_interfaces, partial_derivatives, wxgll, wzgll, ibool);
+  }
+
+  Kokkos::deep_copy(edges, h_edges);
+
+  return;
+}
+
+template <class self_domain_type, class coupled_domain_type>
+void specfem::coupled_interface::coupled_interface<
+    self_domain_type, coupled_domain_type>::compute_coupling() {
+
+  Kokkos::parallel_for(
+      "specfem::coupled_interfaces::coupled_interfaces::compute_coupling",
+      specfem::kokkos::DeviceTeam(this->num_edges, Kokkos::AUTO, 1),
+      KOKKOS_CLASS_LAMBDA(
+          const specfem::kokkos::DeviceTeam::member_type &team_member) {
+        // Get the edge
+        auto &edge = this->edges(team_member.league_rank());
+        auto &[edge1l, edge2l] = edge.get_edges();
+        // Get the number of points along the edge
+        auto npoints =
+            specfem::compute::coupled_interfaces::iterator::get_npoints(
+                edge1l, this->ngllx, this->ngllz);
+
+        // Iterate over the edges using TeamThreadRange
+        Kokkos::parallel_for(
+            Kokkos::TeamThreadRange(team_member, npoints),
+            [=](const int &ipoint) { edge.compute_coupling(ipoint); });
+      });
+
+  return;
+}
+
+#endif // _COUPLED_INTERFACE_TPP

include/coupled_interface/impl/edge/edge.hpp

@@ -0,0 +1,12 @@
+#ifndef _COUPLED_INTERFACE_IMPL_EDGE_HPP
+#define _COUPLED_INTERFACE_IMPL_EDGE_HPP
+
+namespace specfem {
+namespace coupled_interface {
+namespace impl {
+template <class self_domain, class coupled_domain> class edge {};
+} // namespace impl
+} // namespace coupled_interface
+} // namespace specfem
+
+#endif // _COUPLED_INTERFACE_IMPL_EDGE_HPP

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

@@ -0,0 +1,63 @@
+#ifndef _COUPLED_INTERFACE_IMPL_ACOUSTIC_ELASTIC_EDGE_HPP
+#define _COUPLED_INTERFACE_IMPL_ACOUSTIC_ELASTIC_EDGE_HPP
+
+#include "coupled_interface/impl/edge/edge.hpp"
+#include "kokkos_abstractions.h"
+#include "specfem_enums.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+
+namespace specfem {
+namespace coupled_interface {
+namespace impl {
+template <typename qp_type>
+class edge<
+    specfem::domain::domain<specfem::enums::element::medium::acoustic, qp_type>,
+    specfem::domain::domain<specfem::enums::element::medium::elastic,
+                            qp_type> > {
+public:
+  using self_medium = typename specfem::domain::domain<
+      specfem::enums::element::medium::acoustic, qp_type>::medium_type;
+  using coupled_medium =
+      typename specfem::domain::domain<specfem::enums::element::medium::elastic,
+                                       qp_type>::medium_type;
+  using quadrature_points_type = qp_type;
+
+  edge() = default;
+  edge(const int &inum_edge,
+       const specfem::domain::domain<specfem::enums::element::medium::acoustic,
+                                     qp_type> &self_domain,
+       const specfem::domain::domain<specfem::enums::element::medium::elastic,
+                                     qp_type> &coupled_domain,
+       const qp_type &quadrature_points,
+       const specfem::compute::coupled_interfaces::coupled_interfaces
+           &coupled_interfaces,
+       const specfem::compute::partial_derivatives &partial_derivatives,
+       const specfem::kokkos::DeviceView1d<type_real> wxgll,
+       const specfem::kokkos::DeviceView1d<type_real> wzgll,
+       const specfem::kokkos::DeviceView3d<int> ibool);
+
+  void compute_coupling(const int &ipoint);
+
+private:
+  specfem::kokkos::DeviceView2d<int> self_ibool;
+  specfem::kokkos::DeviceView2d<int> coupled_ibool;
+  specfem::kokkos::DeviceView2d<type_real> xix;
+  specfem::kokkos::DeviceView2d<type_real> xiz;
+  specfem::kokkos::DeviceView2d<type_real> gammax;
+  specfem::kokkos::DeviceView2d<type_real> gammaz;
+  specfem::kokkos::DeviceView2d<type_real> jacobian;
+  specfem::enums::coupling::edge::type self_edge_type;
+  specfem::enums::coupling::edge::type coupled_edge_type;
+  specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>
+      self_field_dot_dot;
+  specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft> coupled_field;
+  qp_type quadrature_points;
+  specfem::kokkos::DeviceView1d<type_real> wxgll;
+  specfem::kokkos::DeviceView1d<type_real> wzgll;
+};
+} // namespace impl
+} // namespace coupled_interface
+} // namespace specfem
+
+#endif // _COUPLED_INTERFACE_IMPL_ACOUSTIC_ELASTIC_EDGE_HPP

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

@@ -0,0 +1,135 @@
+#ifndef _COUPLED_INTERFACE_IMPL_ACOUSTIC_ELASTIC_TPP
+#define _COUPLED_INTERFACE_IMPL_ACOUSTIC_ELASTIC_TPP
+
+#include "compute/interface.hpp"
+#include "coupled_interface/impl/edge/edge.hpp"
+#include "coupled_interface/impl/edge/elastic_acoustic/acoustic_elastic.hpp"
+#include "domain/interface.hpp"
+#include "kokkos_abstractions.h"
+#include "macros.hpp"
+#include "specfem_enums.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+
+template <typename qp_type>
+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> >::
+    edge(const int &inum_edge,
+         const specfem::domain::domain<
+             specfem::enums::element::medium::acoustic, qp_type> &self_domain,
+         const specfem::domain::domain<specfem::enums::element::medium::elastic,
+                                       qp_type> &coupled_domain,
+         const qp_type &quadrature_points,
+         const specfem::compute::coupled_interfaces::coupled_interfaces
+             &coupled_interfaces,
+         const specfem::compute::partial_derivatives &partial_derivatives,
+         const specfem::kokkos::DeviceView1d<type_real> wxgll,
+         const specfem::kokkos::DeviceView1d<type_real> wzgll,
+         const specfem::kokkos::DeviceView3d<int> ibool)
+    : quadrature_points(quadrature_points), wxgll(wxgll), wzgll(wzgll) {
+
+  int ngllx, ngllz;
+  quadrature_points.get_ngll(&ngllx, &ngllz);
+
+#ifndef NDEBUG
+  assert(ibool.extent(1) == ngllz);
+  assert(ibool.extent(2) == ngllx);
+  assert(partial_derivatives.xix.extent(1) == ngllz);
+  assert(partial_derivatives.xix.extent(2) == ngllx);
+  assert(partial_derivatives.xiz.extent(1) == ngllz);
+  assert(partial_derivatives.xiz.extent(2) == ngllx);
+  assert(partial_derivatives.gammax.extent(1) == ngllz);
+  assert(partial_derivatives.gammax.extent(2) == ngllx);
+  assert(partial_derivatives.gammaz.extent(1) == ngllz);
+  assert(partial_derivatives.gammaz.extent(2) == ngllx);
+  assert(partial_derivatives.jacobian.extent(1) == ngllz);
+  assert(partial_derivatives.jacobian.extent(2) == ngllx);
+  assert(wxgll.extent(0) == ngllx);
+  assert(wzgll.extent(0) == ngllz);
+#endif
+
+  const int self_ispec =
+      coupled_interfaces.elastic_acoustic.acoustic_ispec(inum_edge);
+  const int coupled_ispec =
+      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);
+
+  self_edge_type = coupled_interfaces.elastic_acoustic.acoustic_edge(inum_edge);
+  coupled_edge_type =
+      coupled_interfaces.elastic_acoustic.elastic_edge(inum_edge);
+
+  self_field_dot_dot = self_domain.get_field_dot_dot();
+  coupled_field = coupled_domain.get_field();
+
+#ifndef NDEBUG
+  assert(self_field_dot_dot.extent(1) == self_medium::components);
+  assert(coupled_field.extent(1) == coupled_medium::components);
+#endif
+
+  return;
+}
+
+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) {
+
+  const 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);
+
+  const 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);
+
+  iglob = self_ibool(j, i);
+
+  switch (self_edge_type) {
+  case specfem::enums::coupling::edge::type::LEFT:
+    Kokkos::atomic_add(&self_field_dot_dot(iglob, 0),
+                       -1.0 * wzgll(j) *
+                           (xix(j, i) * jacobian(j, i) * displ_x +
+                            xiz(j, i) * jacobian(j, i) * displ_z));
+    break;
+  case specfem::enums::coupling::edge::type::RIGHT:
+    Kokkos::atomic_add(&self_field_dot_dot(iglob, 0),
+                       wzgll(j) * (xix(j, i) * jacobian(j, i) * displ_x +
+                                   xiz(j, i) * jacobian(j, i) * displ_z));
+    break;
+  case specfem::enums::coupling::edge::type::BOTTOM:
+    Kokkos::atomic_add(&self_field_dot_dot(iglob, 0),
+                       -1.0 * wxgll(i) *
+                           (gammax(j, i) * jacobian(j, i) * displ_x +
+                            gammaz(j, i) * jacobian(j, i) * displ_z));
+
+    break;
+  case specfem::enums::coupling::edge::type::TOP:
+    Kokkos::atomic_add(&self_field_dot_dot(iglob, 0),
+                       wxgll(i) * (gammax(j, i) * jacobian(j, i) * displ_x +
+                                   gammaz(j, i) * jacobian(j, i) * displ_z));
+    break;
+  default:
+    break;
+  }
+
+  return;
+}
+
+#endif // _COUPLED_INTERFACE_IMPL_ACOUSTIC_ELASTIC_TPP