diff --git a/core/specfem/attenuation.hpp b/core/specfem/attenuation.hpp index cfd70db8f..67dacb49f 100644 --- a/core/specfem/attenuation.hpp +++ b/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" /** diff --git a/core/specfem/attenuation/compute_factors.hpp b/core/specfem/attenuation/compute_factors.hpp new file mode 100644 index 000000000..24e89f468 --- /dev/null +++ b/core/specfem/attenuation/compute_factors.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include "constants.hpp" +#include "specfem_setup.hpp" +#include + +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 struct AttenuationPropertyValues { + Kokkos::View 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 +AttenuationPropertyValues get_attenuation_property_values( + Kokkos::View + tau_s, + Kokkos::View + 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 +type_real get_attenuation_scale_factor( + type_real f_c_source, + Kokkos::View + tau_eps, + Kokkos::View + tau_sigma, + type_real Q_val, type_real attenuation_f0_reference); + +} // namespace attenuation +} // namespace specfem diff --git a/core/specfem/attenuation/compute_factors.tpp b/core/specfem/attenuation/compute_factors.tpp new file mode 100644 index 000000000..e9b55845a --- /dev/null +++ b/core/specfem/attenuation/compute_factors.tpp @@ -0,0 +1,94 @@ +#pragma once +#include "compute_factors.hpp" +#include "constants.hpp" +#include "specfem_setup.hpp" +#include +#include +#include +#include + +namespace specfem { +namespace attenuation { + +template +AttenuationPropertyValues get_attenuation_property_values( + Kokkos::View + tau_s, + Kokkos::View + tau_eps) { + + AttenuationPropertyValues result; + result.beta = + Kokkos::View( + "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 +type_real get_attenuation_scale_factor( + type_real f_c_source, + Kokkos::View + tau_eps, + Kokkos::View + 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 diff --git a/docs/sections/api/specfem/attenuation/compute_factors.rst b/docs/sections/api/specfem/attenuation/compute_factors.rst new file mode 100644 index 000000000..6fc9f14ed --- /dev/null +++ b/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 diff --git a/docs/sections/api/specfem/attenuation/index.rst b/docs/sections/api/specfem/attenuation/index.rst index 684987efe..223f7eb93 100644 --- a/docs/sections/api/specfem/attenuation/index.rst +++ b/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 \ No newline at end of file + optimization.rst diff --git a/tests/unit-tests/CMakeLists.txt b/tests/unit-tests/CMakeLists.txt index d1afb4c5c..18dd6c8e9 100644 --- a/tests/unit-tests/CMakeLists.txt +++ b/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( diff --git a/tests/unit-tests/attenuation/compute_factors_tests.cpp b/tests/unit-tests/attenuation/compute_factors_tests.cpp new file mode 100644 index 000000000..f9a1d189b --- /dev/null +++ b/tests/unit-tests/attenuation/compute_factors_tests.cpp @@ -0,0 +1,381 @@ +#include "specfem/attenuation.hpp" +#include "specfem/utilities/is_close.hpp" +#include "test_macros.hpp" +#include +#include + +using specfem::attenuation::get_attenuation_property_values; +using specfem::attenuation::get_attenuation_scale_factor; +using specfem::utilities::is_close; + +// ============================================================================= +// get_attenuation_property_values tests +// +// All tests use hand-picked values with analytically known results. +// ============================================================================= + +// Simple integer ratios: +// tau_s = [1, 2, 5], tau_eps = [2, 3, 10] +// ratio = [2, 1.5, 2] +// beta = [1, 0.5, 1] (ratio - 1) +// one_minus_sum_beta = 2 + 1.5 + 2 = 5.5 +TEST(Attenuation_PropertyValues, SimpleIntegerRatios) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_s("tau_s"); + Kokkos::View + tau_eps("tau_eps"); + + tau_s(0) = 1.0; + tau_s(1) = 2.0; + tau_s(2) = 5.0; + tau_eps(0) = 2.0; + tau_eps(1) = 3.0; + tau_eps(2) = 10.0; + + auto result = get_attenuation_property_values(tau_s, tau_eps); + + EXPECT_EQ(result.beta.extent(0), N_SLS); + EXPECT_TRUE(is_close(result.beta(0), static_cast(1.0))) + << expected_got(static_cast(1.0), result.beta(0)); + EXPECT_TRUE(is_close(result.beta(1), static_cast(0.5))) + << expected_got(static_cast(0.5), result.beta(1)); + EXPECT_TRUE(is_close(result.beta(2), static_cast(1.0))) + << expected_got(static_cast(1.0), result.beta(2)); + EXPECT_TRUE(is_close(result.one_minus_sum_beta, + static_cast(5.5))) + << expected_got(static_cast(5.5), result.one_minus_sum_beta); +} + +// Uniform ratio: +// tau_s = [2, 4, 5], tau_eps = [4, 8, 10] +// ratio = [2, 2, 2] +// beta = [1, 1, 1] +// one_minus_sum_beta = 6 +TEST(Attenuation_PropertyValues, UniformRatio) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_s("tau_s"); + Kokkos::View + tau_eps("tau_eps"); + + tau_s(0) = 2.0; + tau_s(1) = 4.0; + tau_s(2) = 5.0; + tau_eps(0) = 4.0; + tau_eps(1) = 8.0; + tau_eps(2) = 10.0; + + auto result = get_attenuation_property_values(tau_s, tau_eps); + + EXPECT_TRUE(is_close(result.beta(0), static_cast(1.0))) + << expected_got(static_cast(1.0), result.beta(0)); + EXPECT_TRUE(is_close(result.beta(1), static_cast(1.0))) + << expected_got(static_cast(1.0), result.beta(1)); + EXPECT_TRUE(is_close(result.beta(2), static_cast(1.0))) + << expected_got(static_cast(1.0), result.beta(2)); + EXPECT_TRUE(is_close(result.one_minus_sum_beta, + static_cast(6.0))) + << expected_got(static_cast(6.0), result.one_minus_sum_beta); +} + +// No attenuation (tau_eps = tau_s): +// ratio = [1, 1, 1] +// beta = [0, 0, 0] +// one_minus_sum_beta = 3 +TEST(Attenuation_PropertyValues, NoAttenuation) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_s("tau_s"); + Kokkos::View + tau_eps("tau_eps"); + + tau_s(0) = 1.0; + tau_s(1) = 2.0; + tau_s(2) = 5.0; + tau_eps(0) = 1.0; + tau_eps(1) = 2.0; + tau_eps(2) = 5.0; + + auto result = get_attenuation_property_values(tau_s, tau_eps); + + for (int i = 0; i < N_SLS; ++i) { + EXPECT_TRUE(is_close(result.beta(i), static_cast(0.0))) + << "beta(" << i << ") should be 0 with no attenuation: " + << expected_got(static_cast(0.0), result.beta(i)); + } + EXPECT_TRUE(is_close(result.one_minus_sum_beta, + static_cast(3.0))) + << expected_got(static_cast(3.0), result.one_minus_sum_beta); +} + +// N_SLS = 2: +// tau_s = [2, 5], tau_eps = [6, 10] +// ratio = [3, 2] +// beta = [2, 1] +// one_minus_sum_beta = 5 +TEST(Attenuation_PropertyValues, TwoSLS) { + constexpr int N_SLS = 2; + + Kokkos::View + tau_s("tau_s"); + Kokkos::View + tau_eps("tau_eps"); + + tau_s(0) = 2.0; + tau_s(1) = 5.0; + tau_eps(0) = 6.0; + tau_eps(1) = 10.0; + + auto result = get_attenuation_property_values(tau_s, tau_eps); + + EXPECT_EQ(result.beta.extent(0), 2); + EXPECT_TRUE(is_close(result.beta(0), static_cast(2.0))) + << expected_got(static_cast(2.0), result.beta(0)); + EXPECT_TRUE(is_close(result.beta(1), static_cast(1.0))) + << expected_got(static_cast(1.0), result.beta(1)); + EXPECT_TRUE(is_close(result.one_minus_sum_beta, + static_cast(5.0))) + << expected_got(static_cast(5.0), result.one_minus_sum_beta); +} + +// N_SLS = 5: +// tau_s = [1, 1, 1, 1, 1], tau_eps = [1.5, 1.5, 1.5, 1.5, 1.5] +// ratio = [1.5, 1.5, 1.5, 1.5, 1.5] +// beta = [0.5, 0.5, 0.5, 0.5, 0.5] +// one_minus_sum_beta = 7.5 +TEST(Attenuation_PropertyValues, FiveSLS) { + constexpr int N_SLS = 5; + + Kokkos::View + tau_s("tau_s"); + Kokkos::View + tau_eps("tau_eps"); + + for (int j = 0; j < N_SLS; ++j) { + tau_s(j) = 1.0; + tau_eps(j) = 1.5; + } + + auto result = get_attenuation_property_values(tau_s, tau_eps); + + EXPECT_EQ(result.beta.extent(0), 5); + for (int i = 0; i < N_SLS; ++i) { + EXPECT_TRUE(is_close(result.beta(i), static_cast(0.5))) + << expected_got(static_cast(0.5), result.beta(i)); + } + EXPECT_TRUE(is_close(result.one_minus_sum_beta, + static_cast(7.5))) + << expected_got(static_cast(7.5), result.one_minus_sum_beta); +} + +// ============================================================================= +// get_attenuation_scale_factor tests +// +// All tests use hand-picked values with analytically known results. +// +// The scale factor is: +// scale_factor = factor_scale_mu * factor_scale_mu0 +// where: +// factor_scale_mu0 = 1 + 2*ln(f_c/f_0) / (pi*Q) +// factor_scale_mu = sum_unrelaxed / sum_weighted +// sum_unrelaxed = 1 + sum_i(defect_i / N_SLS) +// sum_weighted = 1 + sum_i(defect_i / (1 + 1/(w*tau_i)^2) / N_SLS) +// defect_i = tau_eps_i/tau_sigma_i - 1 +// w = 2*pi*f_c +// ============================================================================= + +// No attenuation + f_c = f_0: +// tau_eps = tau_sigma => defect = 0 => factor_scale_mu = 1 +// f_c = f_0 => ln(1) = 0 => factor_scale_mu0 = 1 +// scale_factor = 1.0 exactly +TEST(Attenuation_ScaleFactor, NoAttenuationSameFreq) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_eps("tau_eps"); + Kokkos::View + tau_sigma("tau_sigma"); + + tau_sigma(0) = 1.0; + tau_sigma(1) = 2.0; + tau_sigma(2) = 5.0; + tau_eps(0) = 1.0; + tau_eps(1) = 2.0; + tau_eps(2) = 5.0; + + type_real scale_factor = get_attenuation_scale_factor( + 1.0, tau_eps, tau_sigma, 100.0, 1.0); + + EXPECT_TRUE(is_close(scale_factor, static_cast(1.0))) + << expected_got(static_cast(1.0), scale_factor); +} + +// Pure log correction (no defect): +// tau_eps = tau_sigma => factor_scale_mu = 1 +// f_c = e, f_0 = 1, Q = 200 +// factor_scale_mu0 = 1 + 2*ln(e) / (pi*200) = 1 + 2/(200*pi) = 1 + 1/(100*pi) +// scale_factor = 1 + 1/(100*pi) +TEST(Attenuation_ScaleFactor, PureLogCorrection) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_eps("tau_eps"); + Kokkos::View + tau_sigma("tau_sigma"); + + tau_sigma(0) = 1.0; + tau_sigma(1) = 2.0; + tau_sigma(2) = 5.0; + tau_eps(0) = 1.0; + tau_eps(1) = 2.0; + tau_eps(2) = 5.0; + + const type_real e = std::exp(1.0); + const type_real expected = 1.0 + 1.0 / (100.0 * pi); + + type_real scale_factor = get_attenuation_scale_factor( + e, tau_eps, tau_sigma, 200.0, 1.0); + + EXPECT_TRUE(is_close(scale_factor, expected)) + << expected_got(expected, scale_factor); +} + +// Uniform defect with omega*tau = 1: +// N_SLS = 3, all tau_sigma = 1/(2*pi), all tau_eps = 4/(3*2*pi) +// ratio = 4/3, defect = 1/3 +// f_c = 1, f_0 = 1, Q = 100 +// +// omega_tau = 2*pi * 1 * 1/(2*pi) = 1 +// sum_unrelaxed = 1 + (1/3)/3 * 3 = 1 + 1/3 = 4/3 +// weight = 1/(1 + 1/1) = 1/2 +// sum_weighted = 1 + (1/3)*(1/2)/3 * 3 = 1 + 1/6 = 7/6 +// factor_scale_mu = (4/3) / (7/6) = 8/7 +// factor_scale_mu0 = 1 (since f_c = f_0) +// scale_factor = 8/7 +TEST(Attenuation_ScaleFactor, UniformDefectOmegaTauOne) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_eps("tau_eps"); + Kokkos::View + tau_sigma("tau_sigma"); + + const type_real tau_s_val = 1.0 / (2.0 * pi); + const type_real tau_e_val = 4.0 / (3.0 * 2.0 * pi); + + for (int j = 0; j < N_SLS; ++j) { + tau_sigma(j) = tau_s_val; + tau_eps(j) = tau_e_val; + } + + type_real scale_factor = get_attenuation_scale_factor( + 1.0, tau_eps, tau_sigma, 100.0, 1.0); + + const type_real expected = static_cast(8.0 / 7.0); + EXPECT_TRUE(is_close(scale_factor, expected)) + << expected_got(expected, scale_factor); +} + +// Combined log correction and defect: +// Same uniform defect as above (scale_mu = 8/7) +// but now f_c = e, f_0 = 1, Q = 200 +// factor_scale_mu0 = 1 + 1/(100*pi) +// scale_factor = (8/7) * (1 + 1/(100*pi)) +TEST(Attenuation_ScaleFactor, CombinedLogAndDefect) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_eps("tau_eps"); + Kokkos::View + tau_sigma("tau_sigma"); + + const type_real tau_s_val = 1.0 / (2.0 * pi); + const type_real tau_e_val = 4.0 / (3.0 * 2.0 * pi); + + for (int j = 0; j < N_SLS; ++j) { + tau_sigma(j) = tau_s_val; + tau_eps(j) = tau_e_val; + } + + const type_real e = std::exp(1.0); + // factor_scale_mu0 uses the same tau_sigma but now f_c = e + // omega_tau = 2*pi*e*1/(2*pi) = e, so omega_tau^2 = e^2 + // Recompute sum_weighted with omega_tau = e: + // weight_denom = 1 + 1/e^2 + // sum_weighted = 1 + (1/3) / (1 + 1/e^2) = 1 + (1/3)*e^2/(e^2+1) + // sum_unrelaxed stays the same: 4/3 + const type_real e2 = e * e; + const type_real expected_sum_weighted = 1.0 + (1.0 / 3.0) * e2 / (e2 + 1.0); + const type_real expected_factor_mu = (4.0 / 3.0) / expected_sum_weighted; + const type_real expected_factor_mu0 = 1.0 + 1.0 / (100.0 * pi); + const type_real expected = expected_factor_mu * expected_factor_mu0; + + type_real scale_factor = get_attenuation_scale_factor( + e, tau_eps, tau_sigma, 200.0, 1.0); + + EXPECT_TRUE(is_close(scale_factor, expected)) + << expected_got(expected, scale_factor); +} + +// N_SLS = 2 with distinct tau_sigma values: +// tau_sigma = [1, 2], tau_eps = [1.5, 2.5] +// defect = [0.5, 0.25] +// f_c = 1/(2*pi), f_0 = 1/(2*pi), Q = 100 +// +// omega = 2*pi * 1/(2*pi) = 1 +// omega_tau_0 = 1*1 = 1, omega_tau_1 = 1*2 = 2 +// sum_unrelaxed = 1 + (0.5 + 0.25)/2 = 1 + 0.375 = 1.375 +// w0 = 1/(1 + 1/1) = 0.5, w1 = 1/(1 + 1/4) = 4/5 +// sum_weighted = 1 + (0.5*0.5 + 0.25*0.8)/2 = 1 + (0.25 + 0.2)/2 = 1 + 0.225 = 1.225 +// factor_scale_mu = 1.375 / 1.225 = 55/49 +// factor_scale_mu0 = 1 (f_c = f_0) +// scale_factor = 55/49 +TEST(Attenuation_ScaleFactor, TwoSLSDistinctTau) { + constexpr int N_SLS = 2; + + Kokkos::View + tau_eps("tau_eps"); + Kokkos::View + tau_sigma("tau_sigma"); + + tau_sigma(0) = 1.0; + tau_sigma(1) = 2.0; + tau_eps(0) = 1.5; + tau_eps(1) = 2.5; + + const type_real f_c = 1.0 / (2.0 * pi); + + type_real scale_factor = get_attenuation_scale_factor( + f_c, tau_eps, tau_sigma, 100.0, f_c); + + const type_real expected = static_cast(55.0 / 49.0); + EXPECT_TRUE(is_close(scale_factor, expected)) + << expected_got(expected, scale_factor); +} + +// Throws when scale_factor > 1.5: +// tau_eps = tau_sigma (factor_scale_mu = 1) +// f_c = 10, f_0 = 1, Q = 2 +// factor_scale_mu0 = 1 + 2*ln(10)/(pi*2) ≈ 1.733 > 1.5 +TEST(Attenuation_ScaleFactor, ThrowsWhenOutOfRange) { + constexpr int N_SLS = 3; + + Kokkos::View + tau_eps("tau_eps"); + Kokkos::View + tau_sigma("tau_sigma"); + + for (int j = 0; j < N_SLS; ++j) { + tau_sigma(j) = 1.0; + tau_eps(j) = 1.0; + } + + // factor_scale_mu0 = 1 + 2*ln(10)/(pi*2) ≈ 1.733 > 1.5 + EXPECT_THROW( + get_attenuation_scale_factor(10.0, tau_eps, tau_sigma, 2.0, 1.0), + std::runtime_error); +}