Merge pull request #107 from PrincetonUniversity/adjoint-sources

Adjoint 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/assembly/assembly.hpp

@@ -5,10 +5,11 @@
 #include "compute/compute_mesh.hpp"
 #include "compute/compute_partial_derivatives.hpp"
 #include "compute/compute_receivers.hpp"
-#include "compute/compute_sources.hpp"
+// #include "compute/compute_sources.hpp"
 #include "compute/coupled_interfaces/coupled_interfaces.hpp"
 #include "compute/fields/fields.hpp"
 #include "compute/properties/interface.hpp"
+#include "compute/sources/sources.hpp"
 #include "enumerations/specfem_enums.hpp"
 #include "mesh/mesh.hpp"
 #include "receiver/interface.hpp"

include/compute/compute_sources.hpp

@@ -43,6 +43,9 @@ struct sources {
   specfem::kokkos::HostMirror1d<int> h_ispec_array; ///< Spectral element number
                                                     ///< where the source lies
                                                     ///< stored on host
+  specfem::kokkos::HostView1d<specfem::wavefield::type>
+      h_wavefield_type; ///< Type of wavefield that each source is associated
+                        ///< with
   /**
    * @brief Default constructor
    *

include/compute/interface.hpp

@@ -6,12 +6,14 @@
 #include "compute_mesh.hpp"
 #include "compute_partial_derivatives.hpp"
 #include "compute_receivers.hpp"
-#include "compute_sources.hpp"
+// #include "compute_sources.hpp"
 // #include "coupled_interfaces.hpp"
 #include "assembly/assembly.hpp"
 #include "coupled_interfaces/coupled_interfaces.hpp"
 #include "coupled_interfaces/interface_container.hpp"
 #include "fields/fields.hpp"
 #include "properties/interface.hpp"
+#include "sources/impl/source_medium.hpp"
+#include "sources/sources.hpp"
 
 #endif

include/compute/sources/impl/source_medium.hpp

@@ -0,0 +1,39 @@
+#ifndef _COMPUTE_SOURCES_IMPL_SOURCE_MEDIUM_HPP
+#define _COMPUTE_SOURCES_IMPL_SOURCE_MEDIUM_HPP
+
+#include "compute/compute_mesh.hpp"
+#include "source/source.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+
+namespace specfem {
+namespace compute {
+namespace impl {
+namespace sources {
+template <specfem::dimension::type Dimension,
+          specfem::element::medium_tag Medium>
+struct source_medium {
+  using medium_type = specfem::medium::medium<Dimension, Medium>;
+
+  source_medium() = default;
+
+  source_medium(
+      const std::vector<std::shared_ptr<specfem::sources::source> > &sources,
+      const specfem::compute::mesh &mesh,
+      const specfem::compute::partial_derivatives &partial_derivatives,
+      const specfem::compute::properties &properties, const type_real t0,
+      const type_real dt, const int nsteps);
+
+  specfem::kokkos::DeviceView1d<int> source_index_mapping;
+  specfem::kokkos::HostMirror1d<int> h_source_index_mapping;
+  specfem::kokkos::DeviceView3d<type_real> source_time_function;
+  specfem::kokkos::HostMirror3d<type_real> h_source_time_function;
+  specfem::kokkos::DeviceView4d<type_real> source_array;
+  specfem::kokkos::HostMirror4d<type_real> h_source_array;
+};
+} // namespace sources
+} // namespace impl
+} // namespace compute
+} // namespace specfem
+
+#endif /* _COMPUTE_SOURCES_IMPL_SOURCE_MEDIUM_HPP */

include/compute/sources/impl/source_medium.tpp

@@ -0,0 +1,52 @@
+#ifndef _COMPUTE_SOURCES_IMPL_SOURCE_MEDIUM_TPP
+#define _COMPUTE_SOURCES_IMPL_SOURCE_MEDIUM_TPP
+
+#include "algorithms/locate_point.hpp"
+#include "point/coordinates.hpp"
+#include "source_medium.hpp"
+#include <Kokkos_Core.hpp>
+
+template <specfem::dimension::type Dimension,
+          specfem::element::medium_tag Medium>
+specfem::compute::impl::sources::source_medium<Dimension, Medium>::
+    source_medium(
+        const std::vector<std::shared_ptr<specfem::sources::source> > &sources,
+        const specfem::compute::mesh &mesh,
+        const specfem::compute::partial_derivatives &partial_derivatives,
+        const specfem::compute::properties &properties, const type_real t0,
+        const type_real dt, const int nsteps)
+    : source_index_mapping("specfem::sources::source_index_mapping",
+                           sources.size()),
+      h_source_index_mapping(Kokkos::create_mirror_view(source_index_mapping)),
+      source_time_function("specfem::sources::source_time_function", nsteps,
+                           sources.size(), medium_type::components),
+      h_source_time_function(Kokkos::create_mirror_view(source_time_function)),
+      source_array("specfem::sources::source_array", sources.size(),
+                   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++) {
+    auto sv_source_array = Kokkos::subview(
+        this->h_source_array, isource, Kokkos::ALL, Kokkos::ALL, Kokkos::ALL);
+    sources[isource]->compute_source_array(mesh, partial_derivatives,
+                                           properties, sv_source_array);
+    auto sv_stf_array = Kokkos::subview(this->h_source_time_function,
+                                        Kokkos::ALL, isource, Kokkos::ALL);
+    sources[isource]->compute_source_time_function(t0, dt, nsteps,
+                                                   sv_stf_array);
+    specfem::point::gcoord2 coord(sources[isource]->get_x(),
+                                  sources[isource]->get_z());
+
+    auto lcoord = specfem::algorithms::locate_point(coord, mesh);
+    this->h_source_index_mapping(isource) = lcoord.ispec;
+  }
+
+  Kokkos::deep_copy(source_array, h_source_array);
+  Kokkos::deep_copy(source_time_function, h_source_time_function);
+  Kokkos::deep_copy(source_index_mapping, h_source_index_mapping);
+
+  return;
+}
+
+#endif /* _COMPUTE_SOURCES_IMPL_SOURCE_MEDIUM_TPP */

include/compute/sources/sources.hpp

@@ -0,0 +1,55 @@
+#ifndef _COMPUTE_SOURCES_2_HPP
+#define _COMPUTE_SOURCES_2_HPP
+
+#include "compute/compute_mesh.hpp"
+#include "compute/compute_partial_derivatives.hpp"
+#include "compute/properties/properties.hpp"
+#include "enumerations/dimension.hpp"
+#include "enumerations/wavefield.hpp"
+#include "impl/source_medium.hpp"
+#include "kokkos_abstractions.h"
+#include "source/source.hpp"
+
+namespace specfem {
+namespace compute {
+struct sources {
+
+  sources() = default;
+
+  sources(
+      const std::vector<std::shared_ptr<specfem::sources::source> > &sources,
+      const specfem::compute::mesh &mesh,
+      const specfem::compute::partial_derivatives &partial_derivatives,
+      const specfem::compute::properties &properties, const type_real t0,
+      const type_real dt, const int nsteps);
+
+  template <specfem::element::medium_tag Medium>
+  inline specfem::compute::impl::sources::source_medium<
+      specfem::dimension::type::dim2, Medium>
+  get_source_medium() const {
+    if constexpr (Medium == specfem::element::medium_tag::acoustic) {
+      return acoustic_sources;
+    } else if constexpr (Medium == specfem::element::medium_tag::elastic) {
+      return elastic_sources;
+    } else {
+      static_assert("Invalid medium type");
+    }
+  }
+
+  int nsources;
+  specfem::kokkos::HostView1d<int> source_domain_index_mapping;
+  specfem::kokkos::HostView1d<specfem::element::medium_tag>
+      source_medium_mapping;
+  specfem::kokkos::HostView1d<specfem::wavefield::type>
+      source_wavefield_mapping;
+  specfem::compute::impl::sources::source_medium<
+      specfem::dimension::type::dim2, specfem::element::medium_tag::acoustic>
+      acoustic_sources;
+  specfem::compute::impl::sources::source_medium<
+      specfem::dimension::type::dim2, specfem::element::medium_tag::elastic>
+      elastic_sources;
+};
+} // namespace compute
+} // namespace specfem
+
+#endif /* _COMPUTE_SOURCES_2_HPP */

include/domain/impl/kernels.tpp

@@ -170,41 +170,36 @@ void allocate_isotropic_sources(
                                                   medium_tag, property_tag,
                                                   qp_type> &isotropic_sources) {
 
-  const auto value = medium_tag;
-  // Create isotropic sources
-  const auto ispec_array = assembly.sources.h_ispec_array;
-
-  // Count the number of sources within this medium
-  int nsources = 0;
-  for (int isource = 0; isource < ispec_array.extent(0); isource++) {
-    const int ispec = ispec_array(isource);
-    if (assembly.properties.h_element_types(ispec) == value) {
-      nsources++;
+  const auto &sources = assembly.sources;
+  const int nsources = sources.nsources;
+
+  int nsources_in_this_domain = 0;
+  for (int isource = 0; isource < sources.nsources; isource++) {
+    if ((sources.source_medium_mapping(isource) == medium_tag) &&
+        (sources.source_wavefield_mapping(isource) == WavefieldType)) {
+      nsources_in_this_domain++;
     }
   }
 
   // Save the index for sources in this domain
-  specfem::kokkos::HostView1d<int> h_source_kernel_index_mapping(
-      "specfem::domain::domain::source_kernel_index_mapping", nsources);
-
-  specfem::kokkos::HostMirror1d<int> h_source_mapping(
-      "specfem::domain::domain::source_mapping", nsources);
+  specfem::kokkos::HostView1d<int> h_source_domain_index_mapping(
+      "specfem::domain::domain::h_source_domain_index_mapping",
+      nsources_in_this_domain);
 
   int index = 0;
-  for (int isource = 0; isource < ispec_array.extent(0); isource++) {
-    const int ispec = ispec_array(isource);
-    if (assembly.properties.h_element_types(ispec) == value) {
-      h_source_kernel_index_mapping(index) = ispec_array(isource);
-      h_source_mapping(index) = isource;
-      index++;
+  for (int isource = 0; isource < nsources; isource++) {
+    const int isources_domain = sources.source_domain_index_mapping(isource);
+    if (sources.source_medium_mapping(isource) == medium_tag &&
+        sources.source_wavefield_mapping(isource) == WavefieldType) {
+      h_source_domain_index_mapping(index) =
+          sources.source_domain_index_mapping(isource);
     }
   }
 
   // Allocate isotropic sources
   isotropic_sources = specfem::domain::impl::kernels::source_kernel<
       WavefieldType, DimensionType, medium_tag, property_tag, qp_type>(
-      assembly, h_source_kernel_index_mapping, h_source_mapping,
-      quadrature_points);
+      assembly, h_source_domain_index_mapping, quadrature_points);
 
   return;
 }

include/domain/impl/sources/acoustic/acoustic2d_isotropic.hpp

@@ -77,8 +77,10 @@ class source<
    * the source
    */
   KOKKOS_INLINE_FUNCTION void compute_interaction(
-      const type_real &stf,
-      const specfem::kokkos::array_type<type_real, 2> &lagrange_interpolant,
+      const specfem::kokkos::array_type<type_real, medium_type::components>
+          &stf,
+      const specfem::kokkos::array_type<type_real, medium_type::components>
+          &lagrange_interpolant,
       specfem::kokkos::array_type<type_real, medium_type::components>
           &acceleration) const;
 };

include/domain/impl/sources/acoustic/acoustic2d_isotropic.tpp

@@ -45,12 +45,12 @@ KOKKOS_INLINE_FUNCTION void specfem::domain::impl::sources::source<
     specfem::element::property_tag::isotropic,
     specfem::enums::element::quadrature::static_quadrature_points<NGLL> >::
     compute_interaction(
-        const type_real &stf,
-        const specfem::kokkos::array_type<type_real, 2> &lagrange_interpolant,
+        const specfem::kokkos::array_type<type_real, 1> &stf,
+        const specfem::kokkos::array_type<type_real, 1> &lagrange_interpolant,
         specfem::kokkos::array_type<type_real, medium_type::components>
             &acceleration) const {
 
-  acceleration[0] = lagrange_interpolant[0] * stf;
+  acceleration[0] = lagrange_interpolant[0] * stf[0];
 
   return;
 }

include/domain/impl/sources/elastic/elastic2d_isotropic.hpp

@@ -79,8 +79,10 @@ class source<
    * the source
    */
   KOKKOS_INLINE_FUNCTION void compute_interaction(
-      const type_real &stf,
-      const specfem::kokkos::array_type<type_real, 2> &lagrange_interpolant,
+      const specfem::kokkos::array_type<type_real, medium_type::components>
+          &stf,
+      const specfem::kokkos::array_type<type_real, medium_type::components>
+          &lagrange_interpolant,
       specfem::kokkos::array_type<type_real, medium_type::components>
           &acceleration) const;
 };

include/domain/impl/sources/elastic/elastic2d_isotropic.tpp

@@ -38,16 +38,16 @@ KOKKOS_INLINE_FUNCTION void specfem::domain::impl::sources::source<
     specfem::element::property_tag::isotropic,
     specfem::enums::element::quadrature::static_quadrature_points<NGLL> >::
     compute_interaction(
-        const type_real &stf,
+        const specfem::kokkos::array_type<type_real, 2> &stf,
         const specfem::kokkos::array_type<type_real, 2> &lagrange_interpolant,
         specfem::kokkos::array_type<type_real, medium_type::components>
             &acceleration) const {
 
   if constexpr (specfem::globals::simulation_wave == specfem::wave::p_sv) {
-    acceleration[0] = lagrange_interpolant[0] * stf;
-    acceleration[1] = lagrange_interpolant[1] * stf;
+    acceleration[0] = lagrange_interpolant[0] * stf[0];
+    acceleration[1] = lagrange_interpolant[1] * stf[1];
   } else if constexpr (specfem::globals::simulation_wave == specfem::wave::sh) {
-    acceleration[0] = lagrange_interpolant[0] * stf;
+    acceleration[0] = lagrange_interpolant[0] * stf[0];
     acceleration[1] = 0;
   }
 

include/domain/impl/sources/kernel.hpp

@@ -29,8 +29,7 @@ class source_kernel {
   source_kernel() = default;
   source_kernel(
       const specfem::compute::assembly &assembly,
-      const specfem::kokkos::HostView1d<int> h_source_kernel_index_mapping,
-      const specfem::kokkos::HostView1d<int> h_source_mapping,
+      const specfem::kokkos::HostView1d<int> h_source_domain_index_mapping,
       const quadrature_point_type quadrature_points);
 
   void compute_source_interaction(const int timestep) const;
@@ -39,15 +38,15 @@ class source_kernel {
   int nsources;
   specfem::compute::points points;
   specfem::compute::quadrature quadrature;
-  specfem::kokkos::DeviceView1d<int> source_kernel_index_mapping;
-  specfem::kokkos::HostMirror1d<int> h_source_kernel_index_mapping;
-  specfem::kokkos::DeviceView1d<int> source_mapping;
+  specfem::kokkos::DeviceView1d<int> source_domain_index_mapping;
+  specfem::kokkos::HostMirror1d<int> h_source_domain_index_mapping;
   Kokkos::View<int * [specfem::element::ntypes], Kokkos::LayoutLeft,
                specfem::kokkos::DevMemSpace>
       global_index_mapping;
   specfem::compute::properties properties;
   specfem::compute::impl::field_impl<medium_type> field;
-  specfem::compute::sources sources;
+  specfem::compute::impl::sources::source_medium<DimensionType, MediumTag>
+      sources;
   quadrature_point_type quadrature_points;
   specfem::domain::impl::sources::source<DimensionType, MediumTag, PropertyTag,
                                          quadrature_point_type>