Issue 1594 part 3 - Implements Assembly::Attenuation and Assembly::FieldDerivativeStorage

core/specfem/assembly/assembly.hpp

@@ -16,7 +16,7 @@
  * specfem::assembly::mesh, @ref specfem::assembly::jacobian_matrix, etc) that
  * store data computed at all GLL points
  * - Data access functions : Each container provides functions to load/store
- * data on device/host (e.g., @c load_on_device , @c store_on_device )
+ * data on device/host (e.g., @c lod_on_device , @c store_on_device )
  *
  */
 namespace specfem::assembly {

core/specfem/assembly/assembly/dim2/assembly.cpp

@@ -29,6 +29,13 @@ specfem::assembly::assembly<specfem::element::dimension_tag::dim2>::assembly(
                                   this->mesh.element_grid.ngllz, this->mesh,
                                   this->element_types, flux_scheme_config };
   this->jacobian_matrix = { this->mesh };
+
+  this->field_derivative_storage = { this->element_types, this->mesh.nspec,
+                                     this->mesh.element_grid.ngllz,
+                                     this->mesh.element_grid.ngllx };
+
+  // this->attenuation = { reference};
+
   this->properties = { this->mesh.nspec,
                        this->mesh.element_grid.ngllz,
                        this->mesh.element_grid.ngllx,

core/specfem/assembly/assembly/dim2/assembly.hpp

@@ -1,10 +1,12 @@
 #pragma once
 
+#include "specfem/assembly/attenuation.hpp"
 #include "specfem/assembly/boundaries.hpp"
 #include "specfem/assembly/boundary_values.hpp"
 #include "specfem/assembly/conforming_interfaces.hpp"
 #include "specfem/assembly/element_intersections.hpp"
 #include "specfem/assembly/element_types.hpp"
+#include "specfem/assembly/field_derivative_storage.hpp"
 #include "specfem/assembly/fields.hpp"
 #include "specfem/assembly/info.hpp"
 #include "specfem/assembly/jacobian_matrix.hpp"
@@ -78,6 +80,17 @@ template <> struct assembly<specfem::element::dimension_tag::dim2> {
    */
   specfem::assembly::jacobian_matrix<dimension_tag> jacobian_matrix;
 
+  /**
+   * @brief Storage for field derivatives at every quadrature point for elements
+   * that require it (e.g., attenuating elements).
+   */
+  specfem::assembly::FieldDerivativeStorage<dimension_tag>
+      field_derivative_storage;
+  /**
+   * @brief Attenuation properties for the mesh at every quadrature point
+   */
+  //   specfem::assembly::attenuation<dimension_tag> attenuation;
+
   /**
    * @brief Material properties for the mesh at every quadrature point
    *

core/specfem/assembly/assembly/dim3/assembly.cpp

@@ -41,6 +41,11 @@ specfem::assembly::assembly<specfem::element::dimension_tag::dim3>::assembly(
 
   this->jacobian_matrix = { this->mesh };
 
+  this->field_derivative_storage = { this->element_types, this->mesh.nspec,
+                                     this->mesh.element_grid.ngllz,
+                                     this->mesh.element_grid.nglly,
+                                     this->mesh.element_grid.ngllx };
+
   this->properties = { nspec, ngllz,          nglly,
                        ngllx, mesh.materials, this->element_types };
 

core/specfem/assembly/assembly/dim3/assembly.hpp

@@ -1,9 +1,11 @@
 #pragma once
 
+#include "specfem/assembly/attenuation.hpp"
 #include "specfem/assembly/boundaries.hpp"
 #include "specfem/assembly/boundary_values.hpp"
 #include "specfem/assembly/conforming_interfaces.hpp"
 #include "specfem/assembly/element_intersections.hpp"
+#include "specfem/assembly/field_derivative_storage.hpp"
 #include "specfem/assembly/fields.hpp"
 #include "specfem/assembly/info.hpp"
 #include "specfem/assembly/jacobian_matrix.hpp"
@@ -71,6 +73,18 @@ template <> struct assembly<specfem::element::dimension_tag::dim3> {
    */
   specfem::assembly::jacobian_matrix<dimension_tag> jacobian_matrix;
 
+  /**
+   * @brief Storage for field derivatives at every quadrature point for elements
+   * that require it (e.g., attenuating elements).
+   */
+  specfem::assembly::FieldDerivativeStorage<dimension_tag>
+      field_derivative_storage;
+
+  /**
+   * @brief Attenuation properties for the mesh at every quadrature point
+   */
+  //   specfem::assembly::attenuation<dimension_tag> attenuation;
+
   /**
    * @brief Material properties for the mesh at every quadrature point
    *

core/specfem/assembly/attenuation.hpp

@@ -0,0 +1,14 @@
+#pragma once
+
+#include "specfem/enums.hpp"
+
+namespace specfem::assembly {
+
+template <specfem::element::dimension_tag DimensionTag> struct Attenuation;
+
+} // namespace specfem::assembly
+
+#include "specfem/assembly/attenuation/dim2/attenuation.hpp"
+#include "specfem/assembly/attenuation/dim3/attenuation.hpp"
+#include "specfem/assembly/attenuation/load_on_device.hpp"
+#include "specfem/assembly/attenuation/store_on_device.hpp"

core/specfem/assembly/attenuation/dim2/attenuation.cpp

@@ -0,0 +1,87 @@
+#include "specfem/assembly/attenuation/dim2/attenuation.hpp"
+#include "specfem/attenuation.hpp"
+#include "specfem/setup.hpp"
+#include "specfem/utilities/band.hpp"
+
+#include <type_traits>
+
+void specfem::assembly::Attenuation<specfem::element::dimension_tag::dim2>::
+    init_memory_variables(
+        const specfem::assembly::element_types<
+            specfem::element::dimension_tag::dim2> &element_types,
+        const specfem::assembly::mesh<specfem::element::dimension_tag::dim2>
+            &mesh,
+        const specfem::mesh::materials<specfem::element::dimension_tag::dim2>
+            &materials,
+        const specfem::units::Hertz fc, const specfem::units::Hertz f0,
+        const specfem::utilities::Band<specfem::units::Hertz> &band,
+        const Kokkos::View<type_real[N_SLS], Kokkos::DefaultHostExecutionSpace>
+            &tau_sigma){
+
+      FOR_EACH_IN_PRODUCT(
+          (DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV),
+           PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(CONSTANT_ISOTROPIC)),
+          CAPTURE(attn_medium) {
+            auto elements = element_types.get_elements_on_host(
+                _medium_tag_, _property_tag_, _attenuation_tag_);
+            _attn_medium_ = { elements, mesh, materials, ngllz,    ngllx,
+                              fc,       f0,   band,      tau_sigma };
+          })
+    }
+
+specfem::assembly::Attenuation<specfem::element::dimension_tag::dim2>::
+    Attenuation(
+        const specfem::units::Hertz reference_frequency,
+        const specfem::utilities::Band<specfem::units::Hertz> band_in,
+        const bool auto_compute_attenuation_band, const type_real deltat,
+        const specfem::assembly::mesh<specfem::element::dimension_tag::dim2>
+            &mesh,
+        const specfem::assembly::element_types<
+            specfem::element::dimension_tag::dim2> &element_types,
+        const specfem::assembly::Info<specfem::element::dimension_tag::dim2>
+            &info,
+        const specfem::mesh::materials<specfem::element::dimension_tag::dim2>
+            &materials)
+    : ngllz(mesh.element_grid.ngllz), ngllx(mesh.element_grid.ngllx),
+      nspec(mesh.nspec), f0(reference_frequency),
+      auto_compute_attenuation_band(auto_compute_attenuation_band) {
+
+  specfem::utilities::Band<specfem::units::Hertz> band =
+      auto_compute_attenuation_band
+          ? attenuation::compute_band<N_SLS>(
+                specfem::units::Seconds(info.largest_minimum_period))
+          : band_in;
+
+  using specfem::units::unit_symbols::Hz;
+
+  // Compute the band center frequency for attenuation scaling (logarithmic
+  // center)
+  auto fc =
+      specfem::utilities::logarithmic_center(band.min.raw(), band.max.raw()) *
+      Hz;
+
+  // Compute tau_sigma once (shared across all elements)
+  auto tau_sigma = specfem::attenuation::compute_tau_sigma<N_SLS>(band);
+
+  // Compute Runge-Kutta memory-variable update coefficients
+  auto rk = specfem::attenuation::compute_integration_factors<N_SLS>(tau_sigma,
+                                                                     deltat);
+
+  this->alpha_rk = rk.alpha;
+  this->beta_rk = rk.beta;
+  this->gamma_rk = rk.gamma;
+
+  // Copy RK coefficients to device for use in device kernels
+  this->d_alpha_rk = Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+                                  Kokkos::DefaultExecutionSpace>("d_alpha_rk");
+  Kokkos::deep_copy(this->d_alpha_rk, this->alpha_rk);
+  this->d_beta_rk = Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+                                 Kokkos::DefaultExecutionSpace>("d_beta_rk");
+  Kokkos::deep_copy(this->d_beta_rk, this->beta_rk);
+  this->d_gamma_rk = Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+                                  Kokkos::DefaultExecutionSpace>("d_gamma_rk");
+  Kokkos::deep_copy(this->d_gamma_rk, this->gamma_rk);
+
+  init_memory_variables(element_types, mesh, materials, fc, f0, band,
+                        tau_sigma);
+}

core/specfem/assembly/attenuation/dim2/attenuation.hpp

@@ -0,0 +1,171 @@
+#pragma once
+
+#include "specfem/assembly/attenuation.hpp"
+#include "specfem/assembly/attenuation/impl/attenuation_medium.hpp"
+#include "specfem/assembly/element_types.hpp"
+#include "specfem/assembly/info.hpp"
+#include "specfem/assembly/mesh.hpp"
+#include "specfem/constants.hpp"
+#include "specfem/data_access/container.hpp"
+#include "specfem/enums.hpp"
+#include "specfem/macros.hpp"
+#include "specfem/mesh/dim2/materials/materials.hpp"
+#include "specfem/setup.hpp"
+#include <Kokkos_Core.hpp>
+
+namespace specfem::assembly {
+
+template <>
+struct Attenuation<specfem::element::dimension_tag::dim2>
+    : public specfem::data_access::Container<
+          specfem::data_access::ContainerType::domain,
+          specfem::data_access::DataClassType::attenuation,
+          specfem::element::dimension_tag::dim2> {
+
+  /**
+   * @name Type Definitions
+   *
+   */
+  ///@{
+
+  /**
+   * @brief Base container type providing data access infrastructure
+   *
+   * @see specfem::data_access::Container
+   */
+  using base_type = specfem::data_access::Container<
+      specfem::data_access::ContainerType::domain,
+      specfem::data_access::DataClassType::attenuation,
+      specfem::element::dimension_tag::dim2>;
+
+  /**
+   * @brief Kokkos view type for per-element index mapping arrays
+   */
+  using IndexViewType = Kokkos::View<int *, Kokkos::DefaultExecutionSpace>;
+  ///@}
+
+  /**
+   * @name Compile Time Definitions
+   *
+   */
+  ///@{
+
+  constexpr static auto dimension_tag =
+      specfem::element::dimension_tag::dim2; ///< Dimension tag for 2D
+  constexpr static int N_SLS =
+      specfem::constants::N_SLS; ///< Number of standard linear solids
+
+  ///@}
+
+  // Number of GLL points and elements
+  int ngllz; ///< Number of GLL points in the z-direction
+  int ngllx; ///< Number of GLL points in the x-direction
+  int nspec; ///< Total number of spectral elements
+
+  // Attenuation parameters
+  specfem::units::Hertz f0;           ///< Reference frequency
+  bool auto_compute_attenuation_band; ///< Whether to auto-compute the
+                                      ///< attenuation band
+
+  // Runge-Kutta attenuation factors (one coefficient per SLS mechanism)
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+               Kokkos::DefaultHostExecutionSpace>
+      alpha_rk;
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+               Kokkos::DefaultHostExecutionSpace>
+      beta_rk;
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+               Kokkos::DefaultHostExecutionSpace>
+      gamma_rk;
+
+  // Device-accessible Runge-Kutta factors for use in device kernels
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+               Kokkos::DefaultExecutionSpace>
+      d_alpha_rk;
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+               Kokkos::DefaultExecutionSpace>
+      d_beta_rk;
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight,
+               Kokkos::DefaultExecutionSpace>
+      d_gamma_rk;
+
+  // One attenuation_medium member per (medium, property) combination
+  FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV),
+                       PROPERTY_TAG(ISOTROPIC),
+                       ATTENUATION_TAG(CONSTANT_ISOTROPIC)),
+                      DECLARE(((specfem::assembly::impl::attenuation_medium,
+                                (_DIMENSION_TAG_, _MEDIUM_TAG_, _PROPERTY_TAG_,
+                                 _ATTENUATION_TAG_)),
+                               attn_medium)))
+
+  Attenuation() = default;
+
+  Attenuation(
+      const specfem::units::Hertz reference_frequency,
+      const specfem::utilities::Band<specfem::units::Hertz> band_in,
+      const bool auto_compute_attenuation_band, const type_real deltat,
+      const specfem::assembly::mesh<specfem::element::dimension_tag::dim2>
+          &mesh,
+      const specfem::assembly::element_types<
+          specfem::element::dimension_tag::dim2> &element_types,
+      const specfem::assembly::Info<specfem::element::dimension_tag::dim2>
+          &info,
+      const specfem::mesh::materials<specfem::element::dimension_tag::dim2>
+          &materials);
+
+  void init_memory_variables(
+      const specfem::assembly::element_types<
+          specfem::element::dimension_tag::dim2> &element_types,
+      const specfem::assembly::mesh<specfem::element::dimension_tag::dim2>
+          &mesh,
+      const specfem::mesh::materials<specfem::element::dimension_tag::dim2>
+          &materials,
+      const specfem::units::Hertz fc, const specfem::units::Hertz f0,
+      const specfem::utilities::Band<specfem::units::Hertz> &band,
+      const Kokkos::View<type_real[N_SLS], Kokkos::DefaultHostExecutionSpace>
+          &tau_sigma);
+
+  /**
+   * @brief Access the attenuation_medium for a given medium and property.
+   */
+  template <specfem::element::medium_tag MediumTag,
+            specfem::element::property_tag PropertyTag>
+  KOKKOS_INLINE_FUNCTION constexpr specfem::assembly::impl::attenuation_medium<
+      specfem::element::dimension_tag::dim2, MediumTag, PropertyTag,
+      specfem::element::attenuation_tag::constant_isotropic> const &
+  get_container() const {
+    FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV),
+                         PROPERTY_TAG(ISOTROPIC),
+                         ATTENUATION_TAG(CONSTANT_ISOTROPIC)),
+                        CAPTURE(attn_medium) {
+                          if constexpr (_medium_tag_ == MediumTag &&
+                                        _property_tag_ == PropertyTag) {
+                            return _attn_medium_;
+                          }
+                        })
+    Kokkos::abort("Invalid medium type detected in attenuation");
+    SUPPRESS_UNREACHABLE(return {};)
+  }
+
+  /**
+   * @brief Copy attenuation data to host mirrors.
+   */
+  void copy_to_host() {
+    FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV),
+                         PROPERTY_TAG(ISOTROPIC),
+                         ATTENUATION_TAG(CONSTANT_ISOTROPIC)),
+                        CAPTURE(attn_medium) { _attn_medium_.copy_to_host(); })
+  }
+
+  /**
+   * @brief Copy attenuation data to device views.
+   */
+  void copy_to_device() {
+    FOR_EACH_IN_PRODUCT(
+        (DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV), PROPERTY_TAG(ISOTROPIC),
+         ATTENUATION_TAG(CONSTANT_ISOTROPIC)),
+        CAPTURE(attn_medium) { _attn_medium_.copy_to_device(); })
+  }
+};
+
+} // namespace specfem::assembly