Merge pull request #1615 from PrincetonUniversity/issue-1591

Issue 1591 - Maxwell Solid, Nelder-Mead-Algorithm, tau_epsilon computation

Author: lsawade

core/specfem/attenuation.hpp

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

core/specfem/attenuation/compute_tau_eps.hpp

@@ -0,0 +1,124 @@
+#pragma once
+
+#include "compute_tau_sigma.hpp"
+#include "constants.hpp"
+#include "maxwell.hpp"
+#include "specfem/optimization.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+#include <cmath>
+
+namespace specfem {
+namespace attenuation {
+
+/**
+ * @brief Number of frequencies used for evaluating attenuation objective
+ *
+ * Matching SPECFEM3D Fortran implementation which uses 100 frequencies.
+ */
+constexpr int NF_ATTENUATION = 100;
+
+/**
+ * @brief Objective function for \f$\tau_\epsilon\f$ optimization
+ *
+ * This callable struct computes the misfit between the achieved \f$Q\f$ (from
+ * the Maxwell solid model) and the target \f$Q\f$ value. It is designed to be
+ * used with the Nelder-Mead optimizer.
+ *
+ * The objective is computed as the sum of squared differences:
+ * \f[
+ *   \sum_i \left( \tan\delta(f_i) - \frac{1}{Q_\text{target}} \right)^2
+ * \f]
+ *
+ * where \f$\tan\delta = B/A\f$ from the Maxwell solid moduli.
+ *
+ * @tparam N_SLS Number of standard linear solids
+ */
+template <int N_SLS> struct AttenuationObjective {
+  type_real Q;  ///< Target quality factor \f$Q\f$
+  type_real iQ; ///< \f$1/Q\f$ (target \f$\tan\delta\f$)
+  Kokkos::View<type_real[NF_ATTENUATION], Kokkos::LayoutRight,
+               Kokkos::HostSpace>
+      f; ///< Evaluation frequencies \f$f\f$ (Hz)
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+      tau_sigma; ///< Stress relaxation times \f$\tau_\sigma\f$
+
+  /**
+   * @brief Evaluate the objective function for given \f$\tau_\epsilon\f$ values
+   *
+   * @param tau_eps Candidate strain relaxation times \f$\tau_\epsilon\f$
+   * @return Sum of squared misfit between achieved and target \f$1/Q\f$
+   */
+  type_real operator()(
+      Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+          tau_eps) const {
+
+    // Compute Maxwell moduli for all frequencies
+    auto moduli = maxwell<NF_ATTENUATION, N_SLS>(f, tau_sigma, tau_eps);
+
+    // Compute sum of squared misfit
+    type_real misfit = 0.0;
+    for (int i = 0; i < NF_ATTENUATION; ++i) {
+      // tan_delta = B / A ≈ 1/Q
+      type_real tan_delta = moduli.B(i) / moduli.A(i);
+      type_real diff = tan_delta - iQ;
+      misfit += diff * diff;
+    }
+
+    return misfit;
+  }
+};
+
+/**
+ * @brief Compute strain relaxation times via simplex optimization
+ *
+ * This function finds the strain relaxation times \f$\tau_\epsilon\f$ that
+ * achieve a target quality factor \f$Q\f$ for a generalized Maxwell solid with
+ * \f$N_\text{SLS}\f$ standard linear solids.
+ *
+ * The algorithm:
+ * 1. Computes \f$\tau_\sigma\f$ from the period range using compute_tau_sigma
+ * 2. Sets up \f$N_F=100\f$ evaluation frequencies equally spaced in
+ * \f$\log_{10}\f$
+ * 3. Uses Nelder-Mead simplex to minimize the misfit between achieved
+ *    and target \f$1/Q\f$ values over the frequency range
+ *
+ * The initial guess for \f$\tau_\epsilon\f$ is \f$\tau_\sigma\f$ (no
+ * attenuation), and the optimization typically converges within a few hundred
+ * iterations.
+ *
+ * @tparam N_SLS Number of standard linear solids (requires \f$N_\text{SLS} >
+ * 1\f$)
+ *
+ * @param Q Target quality factor \f$Q\f$ (must be positive)
+ * @param tau_sigma Pre-computed stress relaxation times \f$\tau_\sigma\f$ from
+ * compute_tau_sigma
+ * @param min_period Minimum period \f$T_\text{min}\f$ (s) for frequency range
+ * @param max_period Maximum period \f$T_\text{max}\f$ (s) for frequency range
+ *
+ * @return View containing \f$N_\text{SLS}\f$ strain relaxation times
+ * \f$\tau_\epsilon\f$
+ *
+ * @note The returned \f$\tau_\epsilon\f$ values should satisfy \f$\tau_\epsilon
+ * > \tau_\sigma\f$ for positive \f$Q\f$ (physical attenuation).
+ *
+ * @code
+ * // Example: compute tau_eps for Q=200 with 3 SLS
+ * constexpr int N_SLS = 3;
+ * type_real Q = 200.0;
+ * type_real min_period = 0.01;
+ * type_real max_period = 10.0;
+ * auto tau_sigma = compute_tau_sigma<N_SLS>(min_period, max_period);
+ * auto tau_eps = compute_tau_eps<N_SLS>(Q, tau_sigma, min_period, max_period);
+ * @endcode
+ */
+template <int N_SLS>
+Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+compute_tau_eps(
+    type_real Q,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_sigma,
+    type_real min_period, type_real max_period);
+
+} // namespace attenuation
+} // namespace specfem

core/specfem/attenuation/compute_tau_eps.tpp

@@ -0,0 +1,72 @@
+#pragma once
+#include "compute_tau_eps.hpp"
+#include "maxwell.hpp"
+#include "constants.hpp"
+#include "specfem/optimization.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_eps(
+    type_real Q,
+    Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+        tau_sigma,
+    type_real min_period, type_real max_period) {
+
+  static_assert(N_SLS > 1,
+                "N_SLS must be greater than 1 for tau_eps computation");
+
+  // Set up evaluation frequencies equally spaced in log10
+  // These span the same range as tau_sigma but with NF_ATTENUATION points
+  Kokkos::View<type_real[NF_ATTENUATION], Kokkos::LayoutRight, Kokkos::HostSpace>
+      f("frequencies");
+
+  const type_real f1 = 1.0 / max_period; // min frequency
+  const type_real f2 = 1.0 / min_period; // max frequency
+  const type_real log_f1 = std::log10(f1);
+  const type_real log_f2 = std::log10(f2);
+  const type_real d_log_f =
+      (log_f2 - log_f1) / (static_cast<type_real>(NF_ATTENUATION) - 1);
+
+  for (int i = 0; i < NF_ATTENUATION; ++i) {
+    f(i) = std::pow(10.0, log_f1 + i * d_log_f);
+  }
+
+  // Create the objective function
+  AttenuationObjective<N_SLS> objective;
+  objective.Q = Q;
+  objective.iQ = 1.0 / Q;
+  objective.f = f;
+  objective.tau_sigma = tau_sigma;
+
+  // Initial guess: tau_eps = tau_sigma * (1 + 2/Q)
+  // This matches the Fortran SPECFEM3D implementation and provides a better
+  // starting point for the optimization (derived from tan_delta ≈ 1/Q at
+  // omega*tau = 1)
+  Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace> x0(
+      "tau_eps_init");
+  for (int j = 0; j < N_SLS; ++j) {
+    x0(j) = tau_sigma(j) + (tau_sigma(j) * 2.0 / Q);
+  }
+
+  // Run Nelder-Mead optimization
+  // Use default tolerances matching Fortran SPECFEM3D
+  optimization::NelderMeadOptions<N_SLS> options;
+  options.x0 = x0;
+  options.max_iterations = -1;  // default max iterations
+  options.tol_f = 1.0e-4;
+  options.tol_x = 1.0e-4;
+
+  auto result = optimization::optimize(optimization::NelderMeadSimplex{},
+                                        objective, options);
+
+  return result();
+}
+
+} // namespace attenuation
+} // namespace specfem

core/specfem/attenuation/maxwell.hpp

@@ -0,0 +1,129 @@
+#pragma once
+
+#include "constants.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+#include <cmath>
+
+namespace specfem {
+namespace attenuation {
+
+/**
+ * @brief Result containing real (\f$A\f$) and imaginary (\f$B\f$) moduli from
+ * Maxwell solid computation
+ *
+ * For a standard linear solid (SLS), the complex modulus can be written as:
+ * \f[
+ *   M^*(\omega) = M_R \left( 1 + \sum_i \frac{(\tau_{\epsilon_i} -
+ * \tau_{\sigma_i}) \omega^2 \tau_{\sigma_i}}{1 + \omega^2 \tau_{\sigma_i}^2}
+ * \right)
+ *              + i M_R \sum_i \frac{(\tau_{\epsilon_i} - \tau_{\sigma_i})
+ * \omega}{1 + \omega^2 \tau_{\sigma_i}^2}
+ * \f]
+ *
+ * \f$A\f$ represents the real part (storage modulus ratio),
+ * \f$B\f$ represents the imaginary part (loss modulus ratio).
+ * The quality factor \f$Q = A / B = 1 / \tan\delta\f$.
+ *
+ * @tparam NF Number of frequencies
+ *
+ */
+template <int NF> struct MaxwellModuli {
+  Kokkos::View<type_real[NF], Kokkos::LayoutRight, Kokkos::HostSpace> A;
+  Kokkos::View<type_real[NF], Kokkos::LayoutRight, Kokkos::HostSpace> B;
+};
+
+/**
+ * @brief Compute Maxwell solid moduli for a series of standard linear solids
+ *
+ * This function computes the real (\f$A\f$) and imaginary (\f$B\f$) parts of
+ * the viscoelastic modulus for a generalized Maxwell solid (also known as
+ * generalized Zener model) consisting of \f$N_\text{SLS}\f$ standard linear
+ * solids in parallel.
+ *
+ * The formulas match Jeroen's Attenuation notes (43)-(44), with 1/L
+ * normalization:
+ *
+ * For angular frequency \f$\omega = 2 \pi f\f$:
+ *
+ * \f[A(\omega) = \frac{1}{L} \sum_{i=1}^{L} \frac{1 + \omega^2
+ * \tau_{\epsilon_i} \tau_{\sigma_i}}{1 + \omega^2 \tau_{\sigma_i}^2}\f]
+ *
+ * \f[B(\omega) = \frac{1}{L} \sum_{i=1}^{L} \frac{\omega (\tau_{\epsilon_i} -
+ * \tau_{\sigma_i})}{1 + \omega^2 \tau_{\sigma_i}^2}\f]
+ *
+ * where \f$L = N_\text{SLS}\f$. At low frequency, \f$A\f$ approaches \f$1\f$.
+ * The quality factor \f$Q = A/B\f$ is independent of the normalization.
+ * \f$B\f$ is independent of the normalization.
+ *
+ * The quality factor \f$Q = A / B\f$, so \f$\tan(\delta) = B / A = 1 / Q\f$
+ *
+ * @tparam NF Number of frequencies to evaluate
+ * @tparam N_SLS Number of standard linear solids
+ *
+ * @param f View containing \f$N_F\f$ frequencies \f$f\f$ (in Hz, not
+ * \f$\log_{10}\f$)
+ * @param tau_s View containing \f$N_\text{SLS}\f$ stress relaxation times
+ * \f$\tau_\sigma\f$
+ * @param tau_eps View containing \f$N_\text{SLS}\f$ strain relaxation times
+ * \f$\tau_\epsilon\f$
+ *
+ * @return MaxwellModuli containing \f$A\f$ (real) and \f$B\f$ (imaginary)
+ * moduli
+ *
+ * @code
+ * // Example: compute moduli for 3 SLS over a frequency range
+ * constexpr int NF = 100;
+ * constexpr int N_SLS = 3;
+ * Kokkos::View<type_real[NF], ...> f("f");
+ * Kokkos::View<type_real[N_SLS], ...> tau_s("tau_s");
+ * Kokkos::View<type_real[N_SLS], ...> tau_eps("tau_eps"); // ... fill in values
+ * ... auto moduli = maxwell<NF, N_SLS>(f, tau_s, tau_eps); // Q at
+ * frequency i is approximately moduli.A(i) / moduli.B(i)
+ * @endcode
+ */
+template <int NF, int N_SLS>
+MaxwellModuli<NF>
+maxwell(Kokkos::View<type_real[NF], Kokkos::LayoutRight, Kokkos::HostSpace> f,
+        Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+            tau_s,
+        Kokkos::View<type_real[N_SLS], Kokkos::LayoutRight, Kokkos::HostSpace>
+            tau_eps) {
+
+  MaxwellModuli<NF> result;
+  result.A =
+      Kokkos::View<type_real[NF], Kokkos::LayoutRight, Kokkos::HostSpace>("A");
+  result.B =
+      Kokkos::View<type_real[NF], Kokkos::LayoutRight, Kokkos::HostSpace>("B");
+
+  for (int i = 0; i < NF; ++i) {
+    // Angular frequency: w = 2 * pi * f
+    const type_real w = 2.0 * pi * f(i);
+    const type_real w2 = w * w;
+
+    type_real A_sum = 0.0;
+    type_real B_sum = 0.0;
+
+    for (int j = 0; j < N_SLS; ++j) {
+      const type_real tau_s_j = tau_s(j);
+      const type_real tau_eps_j = tau_eps(j);
+      const type_real tau_s_j_sq = tau_s_j * tau_s_j;
+      const type_real denom = 1.0 + w2 * tau_s_j_sq;
+
+      // Real part: (1 + w^2 * tau_eps * tau_s) / denom
+      A_sum += (1.0 + w2 * tau_eps_j * tau_s_j) / denom;
+
+      // Imaginary part: w * (tau_eps - tau_s) / denom
+      B_sum += w * (tau_eps_j - tau_s_j) / denom;
+    }
+
+    // Apply 1/L normalization per Jeroen Tromp's Attenuation notes
+    result.A(i) = A_sum / N_SLS;
+    result.B(i) = B_sum / N_SLS;
+  }
+
+  return result;
+}
+
+} // namespace attenuation
+} // namespace specfem