Merge pull request #1608 from PrincetonUniversity/issue-1592

Issue 1592 - Implements computation of tau_sigma for attenuation

Author: lsawade

core/specfem/CMakeLists.txt

@@ -1,6 +1,7 @@
 # Add existing subdirectories
 add_subdirectory(algorithms)
 add_subdirectory(assembly)
+add_subdirectory(attenuation)
 add_subdirectory(jacobian)
 add_subdirectory(kokkos_kernels)
 add_subdirectory(logger)

core/specfem/attenuation.hpp

@@ -0,0 +1,10 @@
+#pragma once
+#include "attenuation/compute_tau_sigma.hpp"
+#include "attenuation/compute_tau_sigma.tpp"
+
+
+/**
+ * @brief Basic functions for the computation of attenuation related parameters.
+ * 
+ */
+namespace specfem::attenuation {} // namespace specfem::attenuation

core/specfem/attenuation/CMakeLists.txt

@@ -0,0 +1,15 @@
+# Create header-only interface library for specfem::attenuation
+add_library(specfem_attenuation INTERFACE)
+
+# Create alias for consistent naming
+add_library(specfem::attenuation ALIAS specfem_attenuation)
+
+# Add include directories so the library can find its headers
+target_include_directories(specfem_attenuation
+    INTERFACE
+        ${CMAKE_SOURCE_DIR}/core
+)
+
+target_link_libraries(specfem_attenuation INTERFACE
+        Kokkos::kokkos
+)

core/specfem/attenuation/compute_tau_sigma.hpp

@@ -0,0 +1,40 @@
+#pragma once
+
+#include "constants.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+#include <cmath>
+
+namespace specfem {
+namespace attenuation {
+
+/**
+ * @brief Compute stress relaxation times \f$\tau_{\sigma}\f$ equally spaced in
+ * log10 frequency.
+ *
+ * The routine constructs a set of relaxation times corresponding to
+ * frequencies that are equally spaced in base-10 logarithmic scale. Let
+ * \f$f_1 = 1/T_{\text{max}}\f$ and \f$f_2 = 1/T_{\text{min}}\f$ be the minimum and
+ * maximum frequencies, where \f$T_{\text{min}}\f$ and \f$T_{\text{max}}\f$ are the input minimum and maximum periods.
+ * and
+ * \f[\Delta = \frac{\log_{10} f_2 - \log_{10} f_1}{N\_SLS - 1}.\f]
+ * For index \f$i\in[0, N\_SLS-1]\f$ the relaxation time is
+ * \f[\tau_{\sigma,i} = \frac{1}{2\pi\,10^{\log_{10} f_1 + i\Delta}}.\f]
+ *
+ * @tparam N_SLS Number of standard linear solids (requires \f$N\_SLS>1\f$)
+ * @param min_period Minimum period (s)
+ * @param max_period Maximum period (s)
+ * @return Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+ *
+ * @code
+ * // Example: build 4 SLS between 0.1s and 10s
+ * auto tau = specfem::attenuation::compute_tau_sigma<4>(0.1_rt, 10.0_rt);
+ * for (int i = 0; i < 4; ++i) std::cout << tau(i) << "\n";
+ * @endcode
+ */
+template <int N_SLS>
+Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+compute_tau_sigma(const type_real min_period, const type_real max_period);
+
+} // namespace attenuation
+} // namespace specfem

core/specfem/attenuation/compute_tau_sigma.tpp

@@ -0,0 +1,39 @@
+#pragma once
+#include "compute_tau_sigma.hpp"
+#include "constants.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+#include <cmath>
+
+namespace specfem {
+namespace attenuation {
+
+template <int N_SLS>
+Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+compute_tau_sigma(const type_real min_period, const type_real max_period) {
+  static_assert(N_SLS > 1, "N_SLS must be greater than 1 to avoid division by zero");
+
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace> tau_s(
+      "tau_sigma");
+
+  // min/max frequencies
+  const type_real f1 = 1.0 / max_period;
+  const type_real f2 = 1.0 / min_period;
+
+  // logarithms
+  const type_real exp1 = std::log10(f1);
+  const type_real exp2 = std::log10(f2);
+
+  // equally spaced in log10 frequency
+  const type_real dexpval =
+      (exp2 - exp1) / (static_cast<type_real>(N_SLS) - 1);
+
+  for (int i = 0; i < N_SLS; ++i) {
+    tau_s(i) = 1.0 / (pi * 2.0 * std::pow(10.0, exp1 + i * dexpval));
+  }
+
+  return tau_s;
+}
+
+} // namespace attenuation
+} // namespace specfem

docs/sections/api/specfem/attenuation/index.rst

@@ -0,0 +1,10 @@
+``specfem::attenuation``
+========================
+
+.. doxygennamespace:: specfem::attenuation
+    :desc-only:
+
+.. toctree::
+    :maxdepth: 1
+
+    tau_sigma.rst

docs/sections/api/specfem/attenuation/tau_sigma.rst

@@ -0,0 +1,4 @@
+``specfem::attenuation::compute_tau_sigma``
+===========================================
+
+.. doxygenfunction:: specfem::attenuation::compute_tau_sigma

docs/sections/api/specfem/index.rst

@@ -9,6 +9,7 @@
 
     algorithms/index
     assembly/index
+    attenuation/index
     boundary_conditions/index
     chunk_edge/index
     chunk_element/index

tests/unit-tests/CMakeLists.txt

@@ -124,6 +124,19 @@ target_link_libraries(
   Kokkos::kokkos
 )
 
+add_executable(
+  attenuation_tests
+  attenuation/runner.cpp
+  attenuation/compute_tau_sigma_tests.cpp
+)
+
+target_link_libraries(
+  attenuation_tests
+  gtest_main
+  specfem::attenuation
+  Kokkos::kokkos
+)
+
 add_executable(
   gll_tests
   gll/gll_tests.cpp
@@ -895,15 +908,14 @@ target_link_libraries(
 # )
 
 
-
-
 # Link to gtest only if MPI is disabled
 if(NOT MPI_PARALLEL)
   include(GoogleTest)
   set(TEST_TARGETS
     abort_tests
     assembly_receivers_tests
     assembly_tests
+    attenuation_tests
     # compute_acoustic_tests
     # compute_coupled_interfaces_tests
     # compute_elastic_tests

tests/unit-tests/attenuation/compute_tau_sigma_tests.cpp

@@ -0,0 +1,132 @@
+#include "specfem/attenuation.hpp"
+#include <cmath>
+#include <gtest/gtest.h>
+
+using specfem::attenuation::compute_tau_sigma;
+
+// Test that the function returns the correct number of elements
+TEST(Attenuation_ComputeTauSigma, ReturnsCorrectSize) {
+  constexpr type_real min_period = 0.01;
+  constexpr type_real max_period = 10.0;
+
+  auto tau_s_3 = compute_tau_sigma<3>(min_period, max_period);
+  EXPECT_EQ(tau_s_3.extent(0), 3);
+
+  auto tau_s_5 = compute_tau_sigma<5>(min_period, max_period);
+  EXPECT_EQ(tau_s_5.extent(0), 5);
+}
+
+// Test that tau_sigma values are positive
+TEST(Attenuation_ComputeTauSigma, ValuesArePositive) {
+  constexpr type_real min_period = 0.01;
+  constexpr type_real max_period = 10.0;
+
+  auto tau_s = compute_tau_sigma<3>(min_period, max_period);
+
+  for (int i = 0; i < 3; ++i) {
+    EXPECT_GT(tau_s(i), 0.0) << "tau_s(" << i << ") should be positive";
+  }
+}
+
+// Test that tau_sigma values are monotonically decreasing
+// (higher frequency = smaller tau)
+TEST(Attenuation_ComputeTauSigma, ValuesAreMonotonicallyDecreasing) {
+  constexpr type_real min_period = 0.01;
+  constexpr type_real max_period = 10.0;
+
+  auto tau_s = compute_tau_sigma<5>(min_period, max_period);
+
+  for (int i = 1; i < 5; ++i) {
+    EXPECT_LT(tau_s(i), tau_s(i - 1))
+        << "tau_s(" << i << ") should be less than tau_s(" << i - 1 << ")";
+  }
+}
+
+// Test that tau_sigma values are equally spaced in log10 frequency
+TEST(Attenuation_ComputeTauSigma, EquallySpacedInLog10Frequency) {
+  constexpr type_real min_period = 0.01;
+  constexpr type_real max_period = 10.0;
+  constexpr int N_SLS = 5;
+
+  auto tau_s = compute_tau_sigma<N_SLS>(min_period, max_period);
+
+  // Convert tau_s to frequencies: f = 1 / (2 * pi * tau_s)
+  // Then check that log10(f) values are equally spaced
+  std::array<type_real, N_SLS> log_freq;
+  for (int i = 0; i < N_SLS; ++i) {
+    type_real freq = 1.0 / (2.0 * pi * tau_s(i));
+    log_freq[i] = std::log10(freq);
+  }
+
+  // Check that differences between consecutive log10(freq) are equal
+  type_real expected_spacing = log_freq[1] - log_freq[0];
+  for (int i = 2; i < N_SLS; ++i) {
+    type_real actual_spacing = log_freq[i] - log_freq[i - 1];
+    EXPECT_NEAR(actual_spacing, expected_spacing, 1e-10)
+        << "Log10 frequency spacing should be constant";
+  }
+}
+
+// Test boundary values match expected frequencies
+TEST(Attenuation_ComputeTauSigma, BoundaryFrequenciesMatch) {
+  constexpr type_real min_period = 0.01;
+  constexpr type_real max_period = 10.0;
+  constexpr int N_SLS = 3;
+
+  auto tau_s = compute_tau_sigma<N_SLS>(min_period, max_period);
+
+  // Expected frequencies at boundaries
+  type_real f_min = 1.0 / max_period; // 0.1 Hz
+  type_real f_max = 1.0 / min_period; // 100 Hz
+
+  // Convert tau_s to frequency: f = 1 / (2 * pi * tau_s)
+  type_real freq_first = 1.0 / (2.0 * pi * tau_s(0));
+  type_real freq_last = 1.0 / (2.0 * pi * tau_s(N_SLS - 1));
+
+  // First tau_s corresponds to f_min, last to f_max
+  EXPECT_NEAR(freq_first, f_min, f_min * 1e-10)
+      << "First frequency should match f_min";
+  EXPECT_NEAR(freq_last, f_max, f_max * 1e-10)
+      << "Last frequency should match f_max";
+}
+
+// Test with different period ranges
+TEST(Attenuation_ComputeTauSigma, DifferentPeriodRanges) {
+  // Narrow range
+  {
+    auto tau_s = compute_tau_sigma<3>(1.0, 10.0);
+    EXPECT_GT(tau_s(0), 0.0);
+    EXPECT_GT(tau_s(1), 0.0);
+    EXPECT_GT(tau_s(2), 0.0);
+  }
+
+  // Wide range
+  {
+    auto tau_s = compute_tau_sigma<3>(0.001, 100.0);
+    EXPECT_GT(tau_s(0), 0.0);
+    EXPECT_GT(tau_s(1), 0.0);
+    EXPECT_GT(tau_s(2), 0.0);
+  }
+}
+
+// Test N_SLS = 2 edge case (minimum valid value)
+// Note: N_SLS=1 causes division by zero in dexpval calculation
+TEST(Attenuation_ComputeTauSigma, TwoSLS) {
+  constexpr type_real min_period = 0.01;
+  constexpr type_real max_period = 10.0;
+
+  auto tau_s = compute_tau_sigma<2>(min_period, max_period);
+
+  EXPECT_EQ(tau_s.extent(0), 2);
+  EXPECT_GT(tau_s(0), 0.0);
+  EXPECT_GT(tau_s(1), 0.0);
+
+  // First should correspond to f_min, second to f_max
+  type_real f_min = 1.0 / max_period;
+  type_real f_max = 1.0 / min_period;
+  type_real freq_first = 1.0 / (2.0 * pi * tau_s(0));
+  type_real freq_last = 1.0 / (2.0 * pi * tau_s(1));
+
+  EXPECT_NEAR(freq_first, f_min, f_min * 1e-10);
+  EXPECT_NEAR(freq_last, f_max, f_max * 1e-10);
+}