Merge pull request #291 from PrincetonUniversity/issue-276

Added template specializations for acoustic and elastic kernels class

include/domain/impl/kernels.hpp

@@ -14,8 +14,174 @@ namespace kernels {
 template <specfem::wavefield::simulation_field WavefieldType,
           specfem::dimension::type DimensionType,
           specfem::element::medium_tag medium, typename qp_type>
-class kernels {
+class kernels;
 
+template <specfem::wavefield::simulation_field WavefieldType,
+          specfem::dimension::type DimensionType, typename qp_type>
+class kernels<WavefieldType, DimensionType,
+              specfem::element::medium_tag::elastic, qp_type> {
+public:
+  using quadrature_point_type = qp_type;
+  using dimension = specfem::dimension::dimension<DimensionType>;
+  kernels() = default;
+
+  kernels(const type_real dt, const specfem::compute::assembly &assembly,
+          const quadrature_point_type &quadrature_points);
+
+  /**
+   * @brief Compute the interaction of stiffness matrix with wavefield at a time
+   * step
+   *
+   * @param istep Time step
+   */
+  inline void compute_stiffness_interaction(const int istep) const {
+    isotropic_elements.compute_stiffness_interaction(istep);
+    anisotropic_elements.compute_stiffness_interaction(istep);
+    isotropic_elements_dirichlet.compute_stiffness_interaction(istep);
+    anisotropic_elements_dirichlet.compute_stiffness_interaction(istep);
+    isotropic_elements_stacey.compute_stiffness_interaction(istep);
+    anisotropic_elements_stacey.compute_stiffness_interaction(istep);
+    isotropic_elements_stacey_dirichlet.compute_stiffness_interaction(istep);
+    anisotropic_elements_stacey_dirichlet.compute_stiffness_interaction(istep);
+    return;
+  }
+
+  /**
+   * @brief Compute the mass matrix
+   *
+   * @param dt Time step
+   */
+  inline void compute_mass_matrix(const type_real dt) const {
+    isotropic_elements.compute_mass_matrix(dt);
+    anisotropic_elements.compute_mass_matrix(dt);
+    isotropic_elements_dirichlet.compute_mass_matrix(dt);
+    anisotropic_elements_dirichlet.compute_mass_matrix(dt);
+    isotropic_elements_stacey.compute_mass_matrix(dt);
+    anisotropic_elements_stacey.compute_mass_matrix(dt);
+    isotropic_elements_stacey_dirichlet.compute_mass_matrix(dt);
+    anisotropic_elements_stacey_dirichlet.compute_mass_matrix(dt);
+    return;
+  }
+
+  /**
+   * @brief Compute the interaction of source with wavefield at a time step
+   *
+   * @param timestep Time step
+   */
+  inline void compute_source_interaction(const int timestep) const {
+    isotropic_sources.compute_source_interaction(timestep);
+    anisotropic_sources.compute_source_interaction(timestep);
+    return;
+  }
+
+  /**
+   * @brief Compute seismograms at a time step
+   *
+   * @param isig_step Seismogram time step. Same if seismogram is computed at
+   * every time step.
+   */
+  inline void compute_seismograms(const int &isig_step) const {
+    isotropic_receivers.compute_seismograms(isig_step);
+    anisotropic_receivers.compute_seismograms(isig_step);
+    return;
+  }
+
+private:
+  constexpr static specfem::element::boundary_tag dirichlet =
+      specfem::element::boundary_tag::acoustic_free_surface;
+  constexpr static specfem::element::boundary_tag stacey =
+      specfem::element::boundary_tag::stacey;
+  constexpr static specfem::element::boundary_tag none =
+      specfem::element::boundary_tag::none;
+  constexpr static specfem::element::boundary_tag composite_stacey_dirichlet =
+      specfem::element::boundary_tag::composite_stacey_dirichlet;
+  constexpr static specfem::element::property_tag isotropic =
+      specfem::element::property_tag::isotropic;
+  constexpr static specfem::element::property_tag anisotropic =
+      specfem::element::property_tag::anisotropic;
+  constexpr static specfem::element::medium_tag elastic =
+      specfem::element::medium_tag::elastic;
+
+  constexpr static int NGLL = quadrature_point_type::NGLL;
+
+  template <specfem::dimension::type dimension,
+            specfem::element::property_tag property,
+            specfem::element::boundary_tag boundary>
+  using element_kernel = specfem::domain::impl::kernels::element_kernel<
+      WavefieldType, DimensionType, elastic, property, boundary,
+      NGLL>; ///< Underlying element kernel data structure
+
+  template <specfem::dimension::type dimension,
+            specfem::element::property_tag property>
+  using source_kernel = specfem::domain::impl::kernels::source_kernel<
+      WavefieldType, DimensionType, elastic, property,
+      quadrature_point_type>; ///< Underlying source kernel data structure
+
+  template <specfem::dimension::type dimension,
+            specfem::element::property_tag property>
+  using receiver_kernel = specfem::domain::impl::kernels::receiver_kernel<
+      WavefieldType, DimensionType, elastic, property,
+      quadrature_point_type>; ///< Underlying receiver kernel data structure
+
+  element_kernel<DimensionType, isotropic, none>
+      isotropic_elements; ///< Stiffness kernels for isotropic elements
+
+  element_kernel<DimensionType, isotropic, dirichlet>
+      isotropic_elements_dirichlet; ///< Stiffness kernels for isotropic
+                                    ///< elements with Dirichlet boundary
+                                    ///< conditions
+
+  element_kernel<DimensionType, isotropic, stacey>
+      isotropic_elements_stacey; ///< Stiffness kernels for isotropic elements
+                                 ///< with Stacey boundary conditions
+
+  element_kernel<DimensionType, isotropic, composite_stacey_dirichlet>
+      isotropic_elements_stacey_dirichlet; ///< Stiffness kernels for isotropic
+                                           ///< elements with Stacey and
+                                           ///< Dirichlet boundary conditions on
+                                           ///< the same element
+
+  element_kernel<DimensionType, anisotropic, none>
+      anisotropic_elements; ///< Stiffness kernels for anisotropic elements
+
+  element_kernel<DimensionType, anisotropic, dirichlet>
+      anisotropic_elements_dirichlet; ///< Stiffness kernels for anisotropic
+                                      ///< elements with Dirichlet boundary
+                                      ///< conditions
+
+  element_kernel<DimensionType, anisotropic, stacey>
+      anisotropic_elements_stacey; ///< Stiffness kernels for anisotropic
+  ///< elements with Stacey boundary conditions
+
+  element_kernel<DimensionType, anisotropic, composite_stacey_dirichlet>
+      anisotropic_elements_stacey_dirichlet; ///< Stiffness kernels for
+                                             ///< anisotropic elements with
+                                             ///< Stacey and Dirichlet boundary
+                                             ///< conditions on the same element
+
+  source_kernel<DimensionType, isotropic> isotropic_sources; ///< Source kernels
+                                                             ///< for isotropic
+                                                             ///< elements
+
+  source_kernel<DimensionType, anisotropic>
+      anisotropic_sources; ///< Source
+                           ///< kernels for
+                           ///< anisotropic
+                           ///< elements
+
+  receiver_kernel<DimensionType, isotropic>
+      isotropic_receivers; ///< Kernels for computing seismograms within
+                           ///< isotropic elements
+
+  receiver_kernel<DimensionType, anisotropic>
+      anisotropic_receivers; ///< Kernels for computing seismograms within
+                             ///< anisotropic elements
+};
+
+template <specfem::wavefield::simulation_field WavefieldType,
+          specfem::dimension::type DimensionType, typename qp_type>
+class kernels<WavefieldType, DimensionType,
+              specfem::element::medium_tag::acoustic, qp_type> {
 public:
   using quadrature_point_type = qp_type;
   using dimension = specfem::dimension::dimension<DimensionType>;
@@ -83,26 +249,28 @@ class kernels {
       specfem::element::boundary_tag::composite_stacey_dirichlet;
   constexpr static specfem::element::property_tag isotropic =
       specfem::element::property_tag::isotropic;
+  constexpr static specfem::element::medium_tag acoustic =
+      specfem::element::medium_tag::acoustic;
 
   constexpr static int NGLL = quadrature_point_type::NGLL;
 
   template <specfem::dimension::type dimension,
             specfem::element::property_tag property,
             specfem::element::boundary_tag boundary>
   using element_kernel = specfem::domain::impl::kernels::element_kernel<
-      WavefieldType, DimensionType, medium, property, boundary,
+      WavefieldType, DimensionType, acoustic, property, boundary,
       NGLL>; ///< Underlying element kernel data structure
 
   template <specfem::dimension::type dimension,
             specfem::element::property_tag property>
   using source_kernel = specfem::domain::impl::kernels::source_kernel<
-      WavefieldType, DimensionType, medium, property,
+      WavefieldType, DimensionType, acoustic, property,
       quadrature_point_type>; ///< Underlying source kernel data structure
 
   template <specfem::dimension::type dimension,
             specfem::element::property_tag property>
   using receiver_kernel = specfem::domain::impl::kernels::receiver_kernel<
-      WavefieldType, DimensionType, medium, property,
+      WavefieldType, DimensionType, acoustic, property,
       quadrature_point_type>; ///< Underlying receiver kernel data structure
 
   element_kernel<DimensionType, isotropic, none>

include/domain/impl/kernels.tpp

@@ -73,15 +73,16 @@ void allocate_elements(
     }
   }
 
-  if constexpr (wavefield_type == specfem::wavefield::simulation_field::forward ||
-                wavefield_type == specfem::wavefield::simulation_field::adjoint) {
+  if constexpr (wavefield_type ==
+                    specfem::wavefield::simulation_field::forward ||
+                wavefield_type ==
+                    specfem::wavefield::simulation_field::adjoint) {
 
     std::cout << "  - Element type: \n"
               << "    - dimension           : " << dimension::to_string()
               << "\n"
-              << "    - Element type        : " << specfem::element::to_string(
-                                                     medium_tag, property_tag)
-              << "\n"
+              << "    - Element type        : "
+              << specfem::element::to_string(medium_tag, property_tag) << "\n"
               << "    - Boundary Conditions : "
               << specfem::domain::impl::boundary_conditions::print_boundary_tag<
                      boundary_tag>()
@@ -190,12 +191,85 @@ void allocate_isotropic_receivers(
 } // namespace
 
 template <specfem::wavefield::simulation_field WavefieldType,
-          specfem::dimension::type DimensionType,
-          specfem::element::medium_tag medium, typename qp_type>
-specfem::domain::impl::kernels::
-    kernels<WavefieldType, DimensionType, medium, qp_type>::kernels(
-        const type_real dt, const specfem::compute::assembly &assembly,
-        const qp_type &quadrature_points) {
+          specfem::dimension::type DimensionType, typename qp_type>
+specfem::domain::impl::kernels::kernels<
+    WavefieldType, DimensionType, specfem::element::medium_tag::elastic,
+    qp_type>::kernels(const type_real dt,
+                      const specfem::compute::assembly &assembly,
+                      const qp_type &quadrature_points) {
+
+  const int nspec = assembly.mesh.nspec;
+  specfem::kokkos::HostView1d<element_tag> element_tags(
+      "specfem::domain::domain::element_tag", nspec);
+
+  // -----------------------------------------------------------
+  for (int ispec = 0; ispec < nspec; ispec++) {
+    element_tags(ispec) =
+        element_tag(assembly.properties.h_element_types(ispec),
+                    assembly.properties.h_element_property(ispec),
+                    assembly.boundaries.boundary_tags(ispec));
+  }
+
+  if constexpr (WavefieldType ==
+                    specfem::wavefield::simulation_field::forward ||
+                WavefieldType ==
+                    specfem::wavefield::simulation_field::adjoint) {
+    std::cout << " Element Statistics \n"
+              << "------------------------------\n"
+              << "- Types of elements in elastic medium :\n\n";
+  }
+
+  // -----------------------------------------------------------
+
+  // Allocate isotropic elements with dirichlet boundary conditions
+  allocate_elements(assembly, element_tags, isotropic_elements_dirichlet);
+
+  // Allocate aniostropic elements with dirichlet boundary conditions
+  allocate_elements(assembly, element_tags, anisotropic_elements_dirichlet);
+
+  // Allocate isotropic elements with stacey boundary conditions
+  allocate_elements(assembly, element_tags, isotropic_elements_stacey);
+
+  // Allocate anisotropic elements with stacey boundary conditions
+  allocate_elements(assembly, element_tags, anisotropic_elements_stacey);
+
+  // Allocate isotropic elements with stacey dirichlet boundary conditions
+  allocate_elements(assembly, element_tags,
+                    isotropic_elements_stacey_dirichlet);
+
+  // Allocate anisotropic elements with stacey dirichlet boundary conditions
+  allocate_elements(assembly, element_tags,
+                    anisotropic_elements_stacey_dirichlet);
+
+  // Allocate isotropic elements
+  allocate_elements(assembly, element_tags, isotropic_elements);
+
+  // Allocate anisotropic elements
+  allocate_elements(assembly, element_tags, anisotropic_elements);
+
+  // Allocate isotropic sources
+
+  allocate_isotropic_sources(assembly, quadrature_points, isotropic_sources);
+
+  // Allocate isotropic receivers
+
+  allocate_isotropic_receivers(assembly, quadrature_points,
+                               isotropic_receivers);
+
+  // Compute mass matrices
+
+  this->compute_mass_matrix(dt);
+
+  return;
+}
+
+template <specfem::wavefield::simulation_field WavefieldType,
+          specfem::dimension::type DimensionType, typename qp_type>
+specfem::domain::impl::kernels::kernels<
+    WavefieldType, DimensionType, specfem::element::medium_tag::acoustic,
+    qp_type>::kernels(const type_real dt,
+                      const specfem::compute::assembly &assembly,
+                      const qp_type &quadrature_points) {
 
   const int nspec = assembly.mesh.nspec;
   specfem::kokkos::HostView1d<element_tag> element_tags(
@@ -209,12 +283,13 @@ specfem::domain::impl::kernels::
                     assembly.boundaries.boundary_tags(ispec));
   }
 
-  if constexpr (WavefieldType == specfem::wavefield::simulation_field::forward ||
-                WavefieldType == specfem::wavefield::simulation_field::adjoint) {
+  if constexpr (WavefieldType ==
+                    specfem::wavefield::simulation_field::forward ||
+                WavefieldType ==
+                    specfem::wavefield::simulation_field::adjoint) {
     std::cout << " Element Statistics \n"
               << "------------------------------\n"
-              << "- Types of elements in " << specfem::element::to_string(medium)
-              << " medium :\n\n";
+              << "- Types of elements in acoustic medium :\n\n";
   }
 
   // -----------------------------------------------------------

src/domain/impl/elements/kernel.cpp

@@ -11,6 +11,7 @@ constexpr static auto elastic = specfem::element::medium_tag::elastic;
 constexpr static auto acoustic = specfem::element::medium_tag::acoustic;
 
 constexpr static auto isotropic = specfem::element::property_tag::isotropic;
+constexpr static auto anisotropic = specfem::element::property_tag::anisotropic;
 
 constexpr static auto dirichlet =
     specfem::element::boundary_tag::acoustic_free_surface;
@@ -79,8 +80,12 @@ constexpr static auto composite_stacey_dirichlet =
 
 GENERATE_KERNELS(elastic, isotropic, 5)
 
+GENERATE_KERNELS(elastic, anisotropic, 5)
+
 GENERATE_KERNELS(acoustic, isotropic, 5)
 
 GENERATE_KERNELS(elastic, isotropic, 8)
 
+GENERATE_KERNELS(elastic, anisotropic, 8)
+
 GENERATE_KERNELS(acoustic, isotropic, 8)