Compiling code | Runs but fails tests on elastic-acoustic

Implemented Ricker source. I plan to use one of the original example from SPECFEM2D from Wave propagation near a fluid‐solid interface: A spectral‐element approach.

The tests for coupled elastic-acoustic domains are failing.
- Things to do to debug
  - Check if the Ricker source is implemented correctly. Check source time function values
  - Check source location inside the mesh (ispec value)
  - Check the compute forces routine and source interaction routine at the ispec where source is located
  - Check the values at the interface when the wave reaches the interface. Check compute coupling routine

CMakeLists.txt

@@ -160,6 +160,7 @@ target_link_libraries(
 add_library(
         source_time_function
         src/source_time_function/dirac.cpp
+        src/source_time_function/ricker.cpp
 )
 
 target_link_libraries(

examples/CMakeLists.txt

@@ -1,3 +1,4 @@
 cmake_minimum_required(VERSION 3.17.5)
 
 configure_file(homogeneous-medium-flat-topography/specfem_config.yaml.in ${CMAKE_SOURCE_DIR}/examples/homogeneous-medium-flat-topography/specfem_config.yaml)
+configure_file(fluid-solid-interface/specfem_config.yaml.in ${CMAKE_SOURCE_DIR}/examples/fluid-solid-interface/specfem_config.yaml)

examples/fluid-solid-interface/sources.yaml

@@ -0,0 +1,13 @@
+number-of-sources: 1
+sources:
+  - force:
+      x : 1575.0
+      z : 2900.0
+      source_surf: false
+      angle : 0.0
+      vx : 0.0
+      vz : 0.0
+      Ricker:
+        factor: 1e9
+        tshift: 0.0
+        f0: 10.0

examples/fluid-solid-interface/specfem_config.yaml.in

@@ -0,0 +1,50 @@
+parameters:
+
+  header:
+    ## Header information is used for logging. It is good practice to give your simulations explicit names
+    title: Heterogeneous acoustic-elastic medium with 1 acoustic-elastic interface # name for your simulation
+    # A detailed description for your simulation
+    description: |
+      Material systems : Elastic domain (1), Acoustic domain (1)
+      Interfaces : Acoustic-elastic interface (1)
+      Sources : Force source (1)
+      Boundary conditions : Neumann BCs on all edges
+
+  simulation-setup:
+    ## quadrature setup
+    quadrature:
+      alpha: 0.0
+      beta: 0.0
+      ngllx: 5
+      ngllz: 5
+
+    ## Solver setup
+    solver:
+      time-marching:
+        type-of-simulation: forward
+        time-scheme:
+          type: Newmark
+          dt: 0.85e-3
+          nstep: 800
+
+  receivers:
+    stations-file: "@CMAKE_SOURCE_DIR@/examples/fluid_solid_interface/databases/STATIONS"
+    angle: 0.0
+    seismogram-type:
+      - displacement
+      - velocity
+    nstep_between_samples: 1
+
+  ## Runtime setup
+  run-setup:
+    number-of-processors: 1
+    number-of-runs: 1
+
+  ## databases
+  databases:
+    mesh-database: "@CMAKE_SOURCE_DIR@/examples/fluid_solid_interface/databases/database.bin"
+    source-file: "@CMAKE_SOURCE_DIR@/examples/fluid_solid_interface/databases/sources.yaml"
+
+  seismogram:
+    seismogram-format: ascii
+    output-folder: "."

include/compute/coupled_interfaces.hpp

@@ -147,25 +147,16 @@ struct coupled_interfaces {
 
 namespace iterator {
 
-namespace enums {
+int npoints(const specfem::enums::coupling::edge::type &edge, const int ngllx,
+            const int ngllz);
 
-/**
- * @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);
+void self_iterator(const int &ipoint,
+                   const specfem::enums::coupling::edge::type &edge,
+                   const int ngllx, const int ngllz, int &i, int &j);
 
-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);
+void coupled_iterator(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

include/compute/interface.hpp

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

include/coupled_interface/coupled_interface.hpp

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

include/coupled_interface/coupled_interface.tpp

@@ -45,7 +45,7 @@ specfem::coupled_interface::
   }
 
   // Allocate the edges
-  edges = specfem::kokkos::DeviceView1d<specfem::coupled_interface::impl::edge<
+  edges = specfem::kokkos::DeviceView1d<specfem::coupled_interface::impl::edges::edge<
       self_domain_type, coupled_domain_type> >(
       "specfem::coupled_interfaces::coupled_interfaces::edges", num_interfaces);
 
@@ -55,7 +55,7 @@ specfem::coupled_interface::
   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,
+        specfem::coupled_interface::impl::edges::edge<self_domain_type,
                                                coupled_domain_type>(
             inum_edge, self_domain, coupled_domain, quadrature_points,
             coupled_interfaces, partial_derivatives, wxgll, wzgll, ibool);
@@ -86,7 +86,7 @@ void specfem::coupled_interface::coupled_interface<
         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(
+            specfem::compute::coupled_interfaces::iterator::npoints(
                 self_edge_type, ngllx, ngllz);
 
         // Iterate over the edges using TeamThreadRange

include/coupled_interface/impl/edge/edge.hpp

@@ -1,27 +1,55 @@
 #ifndef _COUPLED_INTERFACE_IMPL_EDGE_HPP
 #define _COUPLED_INTERFACE_IMPL_EDGE_HPP
 
+#include "compute/coupled_interfaces.hpp"
+#include "specfem_enums.hpp"
+#include <Kokkos_Core.hpp>
+
 namespace specfem {
 namespace coupled_interface {
 namespace impl {
+namespace edges {
+
+struct self_iterator {
+  specfem::enums::coupling::edge::type edge_type;
+  int ngllx;
+  int ngllz;
+
+  self_iterator() = default;
+
+  self_iterator(const specfem::enums::coupling::edge::type &edge_type,
+                const int &ngllx, const int &ngllz)
+      : edge_type(edge_type), ngllx(ngllx), ngllz(ngllz){};
+
+  KOKKOS_INLINE_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);
+    return;
+  }
+};
 
-using self_iterator_type = decltype(
-    specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges<
-        specfem::compute::coupled_interfaces::iterator::enums::edge::self>);
+struct coupled_iterator {
+  specfem::enums::coupling::edge::type edge_type;
+  int ngllx;
+  int ngllz;
 
-using coupled_iterator_type = decltype(
-    specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges<
-        specfem::compute::coupled_interfaces::iterator::enums::edge::coupled>);
+  coupled_iterator() = default;
 
-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(const specfem::enums::coupling::edge::type &edge_type,
+                   const int &ngllx, const int &ngllz)
+      : edge_type(edge_type), ngllx(ngllx), ngllz(ngllz){};
 
-coupled_iterator_type &coupled_iterator =
-    specfem::compute::coupled_interfaces::iterator::get_points_along_the_edges<
-        specfem::compute::coupled_interfaces::iterator::enums::edge::coupled>;
+  KOKKOS_INLINE_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);
+    return;
+  }
+};
 
 template <class self_domain, class coupled_domain> class edge {};
+} // namespace edges
 } // namespace impl
 } // namespace coupled_interface
 } // namespace specfem

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

@@ -11,6 +11,7 @@
 namespace specfem {
 namespace coupled_interface {
 namespace impl {
+namespace edges {
 template <typename qp_type>
 class edge<
     specfem::domain::domain<specfem::enums::element::medium::acoustic, qp_type>,
@@ -64,7 +65,11 @@ class edge<
   qp_type quadrature_points;
   specfem::kokkos::DeviceView1d<type_real> wxgll;
   specfem::kokkos::DeviceView1d<type_real> wzgll;
+
+  specfem::coupled_interface::impl::edges::self_iterator self_iterator;
+  specfem::coupled_interface::impl::edges::coupled_iterator coupled_iterator;
 };
+} // namespace edges
 } // namespace impl
 } // namespace coupled_interface
 } // namespace specfem

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

@@ -12,7 +12,7 @@
 #include <Kokkos_Core.hpp>
 
 template <typename qp_type>
-specfem::coupled_interface::impl::edge<
+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> >::
@@ -74,6 +74,12 @@ specfem::coupled_interface::impl::edge<
   coupled_edge_type =
       coupled_interfaces.elastic_acoustic.elastic_edge(inum_edge);
 
+  self_iterator = specfem::coupled_interface::impl::edges::self_iterator(
+      self_edge_type, ngllx, ngllz);
+
+  coupled_iterator = specfem::coupled_interface::impl::edges::coupled_iterator(
+      coupled_edge_type, ngllx, ngllz);
+
   self_field_dot_dot = self_domain.get_field_dot_dot();
   coupled_field = coupled_domain.get_field();
 
@@ -86,7 +92,7 @@ specfem::coupled_interface::impl::edge<
 }
 
 template <typename qp_type>
-void specfem::coupled_interface::impl::edge<
+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)
@@ -96,13 +102,13 @@ void specfem::coupled_interface::impl::edge<
   quadrature_points.get_ngll(&ngllx, &ngllz);
 
   int i, j;
-  coupled_iterator(ipoint, coupled_edge_type, ngllx, ngllz, i, j);
+  coupled_iterator(ipoint, i, j);
 
   int iglob = coupled_ibool(j, i);
   const type_real displ_x = coupled_field(iglob, 0);
   const type_real displ_z = coupled_field(iglob, 1);
 
-  self_iterator(ipoint, self_edge_type, ngllx, ngllz, i, j);
+  self_iterator(ipoint, i, j);
 
   iglob = self_ibool(j, i);
 

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

@@ -11,6 +11,7 @@
 namespace specfem {
 namespace coupled_interface {
 namespace impl {
+namespace edges {
 
 template <typename qp_type>
 class edge<
@@ -67,8 +68,11 @@ class edge<
   qp_type quadrature_points;
   specfem::kokkos::DeviceView1d<type_real> wxgll;
   specfem::kokkos::DeviceView1d<type_real> wzgll;
-};
 
+  specfem::coupled_interface::impl::edges::self_iterator self_iterator;
+  specfem::coupled_interface::impl::edges::coupled_iterator coupled_iterator;
+};
+} // namespace edges
 } // namespace impl
 } // namespace coupled_interface
 } // namespace specfem