Added adjoint sources

- Methods to read adjoint source file
- User Defined sources

CMakeLists.txt

@@ -233,6 +233,7 @@ add_library(
         src/source/force_source.cpp
         src/source/moment_tensor_source.cpp
         src/source/adjoint_source.cpp
+        src/source/external.cpp
         src/source/read_sources.cpp
 )
 

include/compute/sources/impl/source_medium.tpp

@@ -22,8 +22,8 @@ specfem::compute::impl::sources::source_medium<Dimension, Medium>::
                            sources.size(), medium_type::components),
       h_source_time_function(Kokkos::create_mirror_view(source_time_function)),
       source_array("specfem::sources::source_array", sources.size(),
-                   mesh.quadratures.gll.N, mesh.quadratures.gll.N,
-                   medium_type::components),
+                   medium_type::components, mesh.quadratures.gll.N,
+                   mesh.quadratures.gll.N),
       h_source_array(Kokkos::create_mirror_view(source_array)) {
 
   for (int isource = 0; isource < sources.size(); isource++) {

include/domain/impl/sources/kernel.tpp

@@ -98,7 +98,7 @@ void specfem::domain::impl::kernels::source_kernel<
               const specfem::kokkos::array_type<type_real,
                                                 medium_type::components>
                   lagrange_interpolant(Kokkos::subview(
-                      sources.source_array, isource_l, iz, ix, Kokkos::ALL));
+                      sources.source_array, isource_l, Kokkos::ALL, iz, ix));
 
               // Source time function
               // For acoustic medium, forward simulation, divide by kappa

include/kokkos_abstractions.h

@@ -572,6 +572,9 @@ template <typename T, int N> struct array_type {
   template <typename MemorySpace, typename Layout, typename MemoryTraits>
   KOKKOS_INLINE_FUNCTION
   array_type(const Kokkos::View<T *, Layout, MemorySpace, MemoryTraits> view) {
+#ifndef NDEBUG
+    assert(view.extent(0) == N);
+#endif
     for (int i = 0; i < N; ++i) {
       data[i] = view(i);
     }

include/source/adjoint_source.hpp

@@ -28,6 +28,8 @@ class adjoint_source : public source {
     return specfem::wavefield::type::adjoint;
   }
 
+  std::string print() const override;
+
 private:
   std::string station_name;
   std::string network_name;

include/source/external.hpp

@@ -0,0 +1,39 @@
+#ifndef _SPECFEM_SOURCES_EXTERNAL_HPP1_
+#define _SPECFEM_SOURCES_EXTERNAL_HPP1_
+
+#include "compute/compute_mesh.hpp"
+#include "compute/compute_partial_derivatives.hpp"
+#include "compute/properties/properties.hpp"
+#include "source.hpp"
+#include "yaml-cpp/yaml.h"
+
+namespace specfem {
+namespace sources {
+class external : public source {
+public:
+  external(){};
+
+  external(YAML::Node &Node, const int nsteps, const type_real dt,
+           const specfem::wavefield::type wavefield_type)
+      : wavefield_type(wavefield_type), specfem::sources::source(Node, nsteps,
+                                                                 dt){};
+
+  void compute_source_array(
+      const specfem::compute::mesh &mesh,
+      const specfem::compute::partial_derivatives &partial_derivatives,
+      const specfem::compute::properties &properties,
+      specfem::kokkos::HostView3d<type_real> source_array) override;
+
+  specfem::wavefield::type get_wavefield_type() const override {
+    return wavefield_type;
+  }
+
+  std::string print() const override;
+
+private:
+  specfem::wavefield::type wavefield_type;
+};
+} // namespace sources
+} // namespace specfem
+
+#endif /* _SPECFEM_SOURCES_EXTERNAL_HPP_ */

include/source/interface.hpp

@@ -2,6 +2,7 @@
 #define _SOURCE_INTERFACE_HPP
 
 #include "adjoint_source.hpp"
+#include "external.hpp"
 #include "force_source.hpp"
 #include "moment_tensor_source.hpp"
 #include "read_sources.hpp"

include/source_time_function/external.hpp

@@ -34,6 +34,8 @@ class external : public stf {
     return ss.str();
   }
 
+  type_real get_t0() const override { return this->__t0; }
+
 private:
   int __nsteps;
   type_real __t0;

src/reader/seismogram.cpp

@@ -20,14 +20,15 @@ void specfem::reader::seismogram::read() {
 
   while (std::getline(file, line)) {
     std::istringstream iss(line);
-    type_real time, value;
+    double time, value;
     if (!(iss >> time >> value)) {
       throw std::runtime_error("Seismogram file " + filename +
                                " is not formatted correctly");
     }
 
     source_time_function(nsteps, 0) = time;
     source_time_function(nsteps, 1) = value;
+    nsteps++;
   }
 
   file.close();

src/source/adjoint_source.cpp

@@ -16,6 +16,7 @@ void specfem::sources::adjoint_source::compute_source_array(
   const auto N = mesh.quadratures.gll.N;
 
   const auto el_type = properties.h_element_types(lcoord.ispec);
+  const int ncomponents = source_array.extent(0);
 
   // Compute lagrange interpolants at the source location
   auto [hxi_source, hpxi_source] =
@@ -31,17 +32,38 @@ void specfem::sources::adjoint_source::compute_source_array(
     for (int ix = 0; ix < N; ++ix) {
       hlagrange = hxi_source(ix) * hgamma_source(iz);
 
-      if (el_type == specfem::element::medium_tag::acoustic ||
-          (el_type == specfem::element::medium_tag::elastic &&
-           specfem::globals::simulation_wave == specfem::wave::sh)) {
+      if (el_type == specfem::element::medium_tag::acoustic) {
+        if (ncomponents != 1) {
+          throw std::runtime_error(
+              "Force source requires 1 component for acoustic medium");
+        }
         source_array(0, iz, ix) = hlagrange;
-        source_array(1, iz, ix) = 0;
-      } else if ((el_type == specfem::element::medium_tag::elastic &&
-                  specfem::globals::simulation_wave == specfem::wave::p_sv) ||
-                 el_type == specfem::element::medium_tag::poroelastic) {
-        source_array(0, iz, ix) = hlagrange;
-        source_array(1, iz, ix) = hlagrange;
+      } else if ((el_type == specfem::element::medium_tag::elastic) ||
+                 (el_type == specfem::element::medium_tag::poroelastic)) {
+        if (ncomponents != 2) {
+          throw std::runtime_error(
+              "Force source requires 2 components for elastic medium");
+        }
+        if (specfem::globals::simulation_wave == specfem::wave::sh) {
+          source_array(0, iz, ix) = hlagrange;
+          source_array(1, iz, ix) = 0;
+        } else {
+          source_array(0, iz, ix) = hlagrange;
+          source_array(1, iz, ix) = hlagrange;
+        }
       }
     }
   }
 };
+
+std::string specfem::sources::adjoint_source::print() const {
+
+  std::ostringstream message;
+  message << "- Adjoint Source: \n"
+          << "    Source Location: \n"
+          << "      x = " << type_real(this->x) << "\n"
+          << "      z = " << type_real(this->z) << "\n"
+          << "    Source Time Function: \n"
+          << this->forcing_function->print() << "\n";
+  return message.str();
+}

src/source/external.cpp

@@ -0,0 +1,83 @@
+#include "algorithms/locate_point.hpp"
+#include "compute/compute_mesh.hpp"
+#include "compute/compute_partial_derivatives.hpp"
+#include "compute/properties/properties.hpp"
+#include "enumerations/specfem_enums.hpp"
+#include "globals.h"
+#include "kokkos_abstractions.h"
+#include "point/coordinates.hpp"
+#include "quadrature/interface.hpp"
+#include "source/interface.hpp"
+#include "source_time_function/interface.hpp"
+#include "specfem_mpi/interface.hpp"
+#include "specfem_setup.hpp"
+// #include "utilities.cpp"
+#include "yaml-cpp/yaml.h"
+#include <cmath>
+
+void specfem::sources::external::compute_source_array(
+    const specfem::compute::mesh &mesh,
+    const specfem::compute::partial_derivatives &partial_derivatives,
+    const specfem::compute::properties &properties,
+    specfem::kokkos::HostView3d<type_real> source_array) {
+
+  specfem::point::gcoord2 coord = specfem::point::gcoord2(this->x, this->z);
+  auto lcoord = specfem::algorithms::locate_point(coord, mesh);
+
+  const auto xi = mesh.quadratures.gll.h_xi;
+  const auto gamma = mesh.quadratures.gll.h_xi;
+  const auto N = mesh.quadratures.gll.N;
+
+  const auto el_type = properties.h_element_types(lcoord.ispec);
+  const int ncomponents = source_array.extent(0);
+
+  // Compute lagrange interpolants at the source location
+  auto [hxi_source, hpxi_source] =
+      specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
+          lcoord.xi, N, xi);
+  auto [hgamma_source, hpgamma_source] =
+      specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
+          lcoord.gamma, N, gamma);
+
+  type_real hlagrange;
+
+  for (int iz = 0; iz < N; ++iz) {
+    for (int ix = 0; ix < N; ++ix) {
+      hlagrange = hxi_source(ix) * hgamma_source(iz);
+
+      if (el_type == specfem::element::medium_tag::acoustic) {
+        if (ncomponents != 1) {
+          throw std::runtime_error(
+              "Force source requires 1 component for acoustic medium");
+        }
+        source_array(0, iz, ix) = hlagrange;
+      } else if ((el_type == specfem::element::medium_tag::elastic) ||
+                 (el_type == specfem::element::medium_tag::poroelastic)) {
+        if (ncomponents != 2) {
+          throw std::runtime_error(
+              "Force source requires 2 components for elastic medium");
+        }
+        if (specfem::globals::simulation_wave == specfem::wave::sh) {
+          source_array(0, iz, ix) = hlagrange;
+          source_array(1, iz, ix) = 0;
+        } else {
+          source_array(0, iz, ix) = hlagrange;
+          source_array(1, iz, ix) = hlagrange;
+        }
+      }
+    }
+  }
+};
+
+std::string specfem::sources::external::print() const {
+
+  std::ostringstream message;
+  message << "- External Source: \n"
+          << "    Source Location: \n"
+          << "      x = " << type_real(this->x) << "\n"
+          << "      z = " << type_real(this->z) << "\n"
+          << "    Source Time Function: \n"
+          << this->forcing_function->print() << "\n";
+
+  return message.str();
+}

src/source/force_source.cpp

@@ -29,7 +29,7 @@ void specfem::sources::force::compute_source_array(
   const auto N = mesh.quadratures.gll.N;
 
   const auto el_type = properties.h_element_types(lcoord.ispec);
-  const int ncomponents = source_array.extent(2);
+  const int ncomponents = source_array.extent(0);
 
   // Compute lagrange interpolants at the source location
   auto [hxi_source, hpxi_source] =

src/source/moment_tensor_source.cpp

@@ -33,7 +33,7 @@ void specfem::sources::moment_tensor::compute_source_array(
         "Moment tensor source not implemented for acoustic medium");
   }
 
-  const int ncomponents = source_array.extent(2);
+  const int ncomponents = source_array.extent(0);
   if (ncomponents != 2) {
     throw std::runtime_error(
         "Moment tensor source requires 2 components for elastic medium");

src/source/read_sources.cpp

@@ -1,3 +1,5 @@
+#include "source/external.hpp"
+#include "source/force_source.hpp"
 #include "source/interface.hpp"
 #include "source/moment_tensor_source.hpp"
 #include "specfem_setup.hpp"
@@ -39,6 +41,9 @@ specfem::sources::read_sources(
     } else if (YAML::Node moment_tensor_source = N["moment-tensor"]) {
       sources.push_back(std::make_shared<specfem::sources::moment_tensor>(
           moment_tensor_source, nsteps, dt, source_wavefield_type));
+    } else if (YAML::Node external = N["user-defined"]) {
+      sources.push_back(std::make_shared<specfem::sources::external>(
+          external, nsteps, dt, source_wavefield_type));
     } else if (YAML::Node adjoint_source = N["adjoint-source"]) {
       sources.push_back(std::make_shared<specfem::sources::adjoint_source>(
           adjoint_source, nsteps, dt));
@@ -59,6 +64,7 @@ specfem::sources::read_sources(
   type_real t0 = std::numeric_limits<type_real>::max();
   for (auto &source : sources) {
     type_real cur_t0 = source->get_t0();
+    std::cout << cur_t0 << std::endl;
     if (cur_t0 < t0) {
       t0 = cur_t0;
     }

src/source_time_function/external.cpp

@@ -10,7 +10,7 @@ specfem::forcing_function::external::external(const YAML::Node &external,
                                               const int nsteps,
                                               const type_real dt)
     : __nsteps(nsteps), __dt(dt),
-      __filename(external["filename"].as<std::string>()) {
+      __filename(external["stf-file"].as<std::string>()) {
 
   if ((external["format"].as<std::string>() == "ascii") ||
       (external["format"].as<std::string>() == "ASCII")) {
@@ -19,7 +19,9 @@ specfem::forcing_function::external::external(const YAML::Node &external,
     throw std::runtime_error("Only ASCII format is supported");
   }
 
-  std::vector<std::string> extentions = { "BXX.semd", "BXY.semd" };
+  std::vector<std::string> extentions = { ".BXX.semd", ".BXY.semd" };
+
+  int file_components = 0;
 
   // get t0 from file
   for (int icomp = 0; icomp < 2; ++icomp) {
@@ -35,10 +37,17 @@ specfem::forcing_function::external::external(const YAML::Node &external,
                                  " is not formatted correctly");
       }
       this->__t0 = time;
+      file_components++;
     }
     file.close();
   }
 
+  if (file_components != 1) {
+    throw std::runtime_error(
+        "Wrong number for files found for external source: " +
+        this->__filename);
+  }
+
   return;
 }
 
@@ -50,9 +59,9 @@ void specfem::forcing_function::external::compute_source_time_function(
 
   const auto extention = [&ncomponents]() -> std::vector<std::string> {
     if (ncomponents == 2) {
-      return { "BXX.semd", "BXZ.semd" };
+      return { ".BXX.semd", ".BXZ.semd" };
     } else if (ncomponents == 1) {
-      return { "BXY.semd" };
+      return { ".BXY.semd" };
     } else {
       throw std::runtime_error("Invalid number of components");
     }
@@ -69,7 +78,7 @@ void specfem::forcing_function::external::compute_source_time_function(
   }
 
   for (int icomp = 0; icomp < ncomponents; ++icomp) {
-    specfem::kokkos::HostView2d<type_real> data("external", nsteps);
+    specfem::kokkos::HostView2d<type_real> data("external", nsteps, 2);
     specfem::reader::seismogram reader(
         this->__filename + extention[icomp],
         specfem::enums::seismogram::format::ascii, data);