Issue 1595 - Implement Modulus and scale_factor computation

core/specfem/attenuation.hpp

@@ -4,6 +4,8 @@
 #include "attenuation/maxwell.hpp"
 #include "attenuation/compute_tau_eps.hpp"
 #include "attenuation/compute_tau_eps.tpp"
+#include "attenuation/compute_factors.hpp"
+#include "attenuation/compute_factors.tpp"
 
 
 /**

core/specfem/attenuation/compute_factors.hpp

@@ -0,0 +1,80 @@
+#pragma once
+
+#include "constants.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+
+namespace specfem {
+namespace attenuation {
+
+/**
+ * @brief Attenuation property values for modulus calculations
+ *
+ * Stores coefficients for computing relaxed and unrelaxed moduli:
+ * - @f$ \beta^{\text{defect}}_i = \tau_{\epsilon_i}/\tau_{\sigma_i} - 1 @f$
+ * (modulus defect per mechanism)
+ * - @f$ \text{OneMinusSumBeta} = \sum_i
+ *   \tau_{\epsilon_i}/\tau_{\sigma_i} @f$
+ *
+ * @tparam N_SLS Number of standard linear solids
+ */
+template <int N_SLS> struct AttenuationPropertyValues {
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace> beta;
+  type_real one_minus_sum_beta;
+};
+
+/**
+ * @brief Compute attenuation property values from relaxation times
+ *
+ * Computes @f$ \beta^{\text{defect}}_i = \tau_{\epsilon_i}/\tau_{\sigma_i} - 1
+ * @f$ and
+ * @f$ \text{OneMinusSumBeta} = 1 - \sum_i \beta_i = \sum_{i=1}^{N\_SLS}
+ * \tau_{\epsilon_i}/\tau_{\sigma_i} @f$ for each standard linear solid.
+ *
+ * @tparam N_SLS Number of standard linear solids
+ * @param tau_s Stress relaxation times @f$ \tau_\sigma @f$
+ * @param tau_eps Strain relaxation times @f$ \tau_\epsilon @f$
+ * @return AttenuationPropertyValues containing @f$ \beta @f$ and
+ *         @f$ \text{one\_minus\_sum\_beta} @f$
+ */
+template <int N_SLS>
+AttenuationPropertyValues<N_SLS> get_attenuation_property_values(
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_s,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_eps);
+
+/**
+ * @brief Compute physical dispersion scaling factor for attenuation
+ *
+ * Computes @f$ \Psi @f$ to scale @f$ \mu_0 @f$ to the unrelaxed modulus:
+ * @f[
+ *   \Psi = \Psi_\mu \times \Psi_{\mu_0}
+ * @f]
+ * where:
+ * - @f$ \Psi_{\mu_0} = 1 + \frac{2}{\pi Q} \ln(f_c / f_0) @f$ corrects for
+ *   logarithmic frequency dependence (Aki & Richards 1980, eq. 5.81)
+ * - @f$ \Psi_\mu = \frac{\sum(1 + \beta_i/N_{\text{SLS}})}{\sum[1 + \beta_i/(1
+ *   + 1/(\omega\tau_{\sigma_i})^2)/N_{\text{SLS}}]} @f$ accounts for SLS
+ *   frequency dispersion
+ *
+ * @tparam N_SLS Number of standard linear solids
+ * @param f_c_source Central frequency of the source @f$ f_c @f$ (Hz)
+ * @param tau_eps Strain relaxation times @f$ \tau_\epsilon @f$
+ * @param tau_sigma Stress relaxation times @f$ \tau_\sigma @f$
+ * @param Q_val Target quality factor @f$ Q @f$
+ * @param attenuation_f0_reference Reference frequency @f$ f_0 @f$ (Hz)
+ * @return Scale factor @f$ \Psi @f$ (expected range [0.5, 1.5])
+ * @throws std::runtime_error if @f$ \Psi @f$ is outside [0.5, 1.5]
+ */
+template <int N_SLS>
+type_real get_attenuation_scale_factor(
+    type_real f_c_source,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_eps,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_sigma,
+    type_real Q_val, type_real attenuation_f0_reference);
+
+} // namespace attenuation
+} // namespace specfem

core/specfem/attenuation/compute_factors.tpp

@@ -0,0 +1,94 @@
+#pragma once
+#include "compute_factors.hpp"
+#include "constants.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+#include <cmath>
+#include <sstream>
+#include <stdexcept>
+
+namespace specfem {
+namespace attenuation {
+
+template <int N_SLS>
+AttenuationPropertyValues<N_SLS> get_attenuation_property_values(
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_s,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_eps) {
+
+  AttenuationPropertyValues<N_SLS> result;
+  result.beta =
+      Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>(
+          "beta");
+
+  // See Komatitsch & Tromp 1999, eq. (7)
+  // Coefficients beta = tau_eps / tau_s
+  for (int i = 0; i < N_SLS; ++i) {
+    result.beta(i) = tau_eps(i) / tau_s(i);
+  }
+
+  // Sum of coefficients beta, then subtract 1 from each beta
+  // to get the modulus defect
+  result.one_minus_sum_beta = 0.0;
+  for (int i = 0; i < N_SLS; ++i) {
+    result.one_minus_sum_beta += result.beta(i);
+    result.beta(i) -= 1.0;
+  }
+
+  return result;
+}
+
+template <int N_SLS>
+type_real get_attenuation_scale_factor(
+    type_real f_c_source,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_eps,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_sigma,
+    type_real Q_val, type_real attenuation_f0_reference) {
+
+  // Quantity by which to scale mu_0 to get mu
+  // See Liu et al. 1976; Aki & Richards 1980, eq. (5.81)
+  const type_real factor_scale_mu0 =
+      1.0 +
+      2.0 * std::log(f_c_source / attenuation_f0_reference) / (pi * Q_val);
+
+  // Quantity by which to scale mu to get mu_unrelaxed
+  type_real sum_unrelaxed = 1.0;
+  type_real sum_weighted = 1.0;
+
+  for (int i = 0; i < N_SLS; ++i) {
+    const type_real defect = tau_eps(i) / tau_sigma(i) - 1.0;
+    sum_unrelaxed += defect / N_SLS;
+
+    const type_real omega_tau = 2.0 * pi * f_c_source * tau_sigma(i);
+    sum_weighted +=
+        defect / (1.0 + 1.0 / (omega_tau * omega_tau)) / N_SLS;
+  }
+
+  const type_real factor_scale_mu = sum_unrelaxed / sum_weighted;
+
+  // Total factor by which to scale mu0 to get mu_unrelaxed
+  const type_real scale_factor = factor_scale_mu * factor_scale_mu0;
+
+  // Check that the correction factor is close to one
+  if (scale_factor < 0.5 || scale_factor > 1.5) {
+    std::ostringstream msg;
+    msg << "Error in get_attenuation_scale_factor(): "
+        << "scale factor = " << scale_factor
+        << " should be between 0.5 and 1.5. "
+        << "factor_scale_mu = " << factor_scale_mu
+        << ", factor_scale_mu0 = " << factor_scale_mu0
+        << ", Q = " << Q_val
+        << ", f_c_source = " << f_c_source
+        << ", attenuation_f0_reference = " << attenuation_f0_reference
+        << ". Please check your reference frequency.";
+    throw std::runtime_error(msg.str());
+  }
+
+  return scale_factor;
+}
+
+} // namespace attenuation
+} // namespace specfem

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

@@ -0,0 +1,8 @@
+``specfem::attenuation::compute_factors``
+=========================================
+
+.. doxygenfunction:: specfem::attenuation::get_attenuation_property_values
+
+.. doxygenfunction:: specfem::attenuation::get_attenuation_scale_factor
+
+.. doxygenstruct:: specfem::attenuation::AttenuationPropertyValues

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

@@ -9,5 +9,6 @@
 
     tau_sigma.rst
     compute_tau_eps.rst
+    compute_factors.rst
     maxwell.rst
-    optimization.rst
+    optimization.rst

tests/unit-tests/CMakeLists.txt

@@ -130,6 +130,7 @@ add_executable(
   attenuation/compute_tau_sigma_tests.cpp
   attenuation/maxwell_tests.cpp
   attenuation/compute_tau_eps_tests.cpp
+  attenuation/compute_factors_tests.cpp
 )
 
 target_link_libraries(