Merge pull request #288 from PrincetonUniversity/issue-287

Added method to compute anisotropic stacey traction
PrincetonUniversity · Dec 13, 2024 · eca3d71 · eca3d71
2 parents b07ec56 + 98af5bb
commit eca3d71
Show file tree

Hide file tree

Showing 2 changed files with 128 additions and 3 deletions.
diff --git a/include/domain/impl/boundary_conditions/stacey/stacey.tpp b/include/domain/impl/boundary_conditions/stacey/stacey.tpp
@@ -21,6 +21,11 @@ using isotropic_type =
     std::integral_constant<specfem::element::property_tag,
                            specfem::element::property_tag::isotropic>;
 
+using anisotropic_type =
+    std::integral_constant<specfem::element::property_tag,
+                           specfem::element::property_tag::anisotropic>;
+
+// Elastic Isotropic Stacey Boundary Conditions not using SIMD types
 template <
     typename PointBoundaryType, typename PointPropertyType,
     typename PointFieldType, typename ViewType,
@@ -58,6 +63,7 @@ impl_enforce_traction(const acoustic_type &, const isotropic_type &,
   return;
 }
 
+// Elastic Isotropic Stacey Boundary Conditions using SIMD types
 template <
     typename PointBoundaryType, typename PointPropertyType,
     typename PointFieldType, typename ViewType,
@@ -100,6 +106,7 @@ impl_enforce_traction(const acoustic_type &, const isotropic_type &,
   return;
 }
 
+// Elastic Isotropic Stacey Boundary Conditions not using SIMD types
 template <
     typename PointBoundaryType, typename PointPropertyType,
     typename PointFieldType, typename ViewType,
@@ -151,6 +158,7 @@ impl_enforce_traction(const elastic_type &, const isotropic_type &,
   return;
 }
 
+// Elastic Isotropic Stacey Boundary Conditions using SIMD types
 template <
     typename PointBoundaryType, typename PointPropertyType,
     typename PointFieldType, typename ViewType,
@@ -209,6 +217,119 @@ impl_enforce_traction(const elastic_type &, const isotropic_type &,
 
   return;
 }
+
+// Elastic Anisotropic stacey boundary conditions not using SIMD typess
+template <
+    typename PointBoundaryType, typename PointPropertyType,
+    typename PointFieldType, typename ViewType,
+    typename std::enable_if_t<!PointBoundaryType::simd::using_simd, int> = 0>
+KOKKOS_FUNCTION void
+impl_enforce_traction(const elastic_type &, const anisotropic_type &,
+                      const PointBoundaryType &boundary,
+                      const PointPropertyType &property,
+                      const PointFieldType &field, ViewType &traction) {
+
+  static_assert(PointBoundaryType::boundary_tag ==
+                    specfem::element::boundary_tag::stacey,
+                "Boundary tag must be stacey");
+
+  static_assert(PointPropertyType::medium_tag ==
+                    specfem::element::medium_tag::elastic,
+                "Medium tag must be elastic");
+
+  static_assert(PointPropertyType::property_tag ==
+                    specfem::element::property_tag::anisotropic,
+                "Property tag must be anisotropic");
+
+  constexpr static auto tag = PointBoundaryType::boundary_tag;
+
+  if (boundary.tag != tag)
+    return;
+
+  const auto vn =
+      specfem::algorithms::dot(field.velocity, boundary.edge_normal);
+
+  const auto &dn = boundary.edge_normal;
+
+  const auto jacobian1d = dn.l2_norm();
+
+  using datatype = std::remove_const_t<decltype(jacobian1d)>;
+
+  datatype factor[2];
+
+  for (int icomp = 0; icomp < 2; ++icomp) {
+    factor[icomp] = ((vn * dn(icomp) / (jacobian1d * jacobian1d)) *
+                     (property.rho_vp - property.rho_vs)) +
+                    field.velocity(icomp) * property.rho_vs;
+  }
+
+  traction(0) += static_cast<type_real>(-1.0) * factor[0] * jacobian1d *
+                 boundary.edge_weight;
+  traction(1) += static_cast<type_real>(-1.0) * factor[1] * jacobian1d *
+                 boundary.edge_weight;
+
+  return;
+}
+
+// Elastic Anisotropic Stacey Boundary Conditions using SIMD types
+template <
+    typename PointBoundaryType, typename PointPropertyType,
+    typename PointFieldType, typename ViewType,
+    typename std::enable_if_t<PointBoundaryType::simd::using_simd, int> = 0>
+KOKKOS_FUNCTION void
+impl_enforce_traction(const elastic_type &, const anisotropic_type &,
+                      const PointBoundaryType &boundary,
+                      const PointPropertyType &property,
+                      const PointFieldType &field, ViewType &traction) {
+
+  static_assert(PointBoundaryType::boundary_tag ==
+                    specfem::element::boundary_tag::stacey,
+                "Boundary tag must be stacey");
+
+  static_assert(PointPropertyType::medium_tag ==
+                    specfem::element::medium_tag::elastic,
+                "Medium tag must be elastic");
+
+  static_assert(PointPropertyType::property_tag ==
+                    specfem::element::property_tag::anisotropic,
+                "Property tag must be anisotropic");
+
+  constexpr int components = PointFieldType::components;
+  constexpr auto tag = PointBoundaryType::boundary_tag;
+
+  using mask_type = typename PointBoundaryType::simd::mask_type;
+
+  mask_type mask([&](std::size_t lane) { return boundary.tag[lane] == tag; });
+
+  if (Kokkos::Experimental::none_of(mask))
+    return;
+
+  const auto vn =
+      specfem::algorithms::dot(field.velocity, boundary.edge_normal);
+  const auto &dn = boundary.edge_normal;
+
+  const auto jacobian1d = dn.l2_norm();
+
+  using datatype = std::remove_const_t<decltype(jacobian1d)>;
+
+  datatype factor[2];
+
+  for (int icomp = 0; icomp < 2; ++icomp) {
+    factor[icomp] = ((vn * dn(icomp) / (jacobian1d * jacobian1d)) *
+                     (property.rho_vp - property.rho_vs)) +
+                    field.velocity(icomp) * property.rho_vs;
+  }
+
+  Kokkos::Experimental::where(mask, traction(0)) =
+      traction(0) + static_cast<type_real>(-1.0) * factor[0] * jacobian1d *
+                        boundary.edge_weight;
+
+  Kokkos::Experimental::where(mask, traction(1)) =
+      traction(1) + static_cast<type_real>(-1.0) * factor[1] * jacobian1d *
+                        boundary.edge_weight;
+
+  return;
+}
 } // namespace
 
 template <typename PointBoundaryType, typename PointPropertyType,

diff --git a/include/point/properties.hpp b/include/point/properties.hpp
@@ -132,7 +132,9 @@ struct properties<specfem::dimension::type::dim2,
   value_type c25; ///< @f$ c_{25} @f$
   ///@}
 
-  value_type rho; ///< Density @f$ \rho @f$
+  value_type rho;    ///< Density @f$ \rho @f$
+  value_type rho_vp; ///< P-wave velocity @f$ \rho v_p @f$
+  value_type rho_vs; ///< S-wave velocity @f$ \rho v_s @f$
 
 private:
   KOKKOS_FUNCTION
@@ -142,7 +144,8 @@ struct properties<specfem::dimension::type::dim2,
              const value_type &c12, const value_type &c23,
              const value_type &c25, const value_type &rho, std::false_type)
       : c11(c11), c13(c13), c15(c15), c33(c33), c35(c35), c55(c55), c12(c12),
-        c23(c23), c25(c25) {}
+        c23(c23), c25(c25), rho(rho), rho_vp(sqrt(rho * c33)),
+        rho_vs(sqrt(rho * c55)) {}
 
   KOKKOS_FUNCTION
   properties(const value_type &c11, const value_type &c13,
@@ -151,7 +154,8 @@ struct properties<specfem::dimension::type::dim2,
              const value_type &c12, const value_type &c23,
              const value_type &c25, const value_type &rho, std::true_type)
       : c11(c11), c13(c13), c15(c15), c33(c33), c35(c35), c55(c55), c12(c12),
-        c23(c23), c25(c25), rho(rho) {}
+        c23(c23), c25(c25), rho(rho), rho_vp(Kokkos::sqrt(rho * c33)),
+        rho_vs(Kokkos::sqrt(rho * c55)) {}
 
 public:
   KOKKOS_FUNCTION