Viscoelastic wave propagator

[sframba]

This PR introduces attenuation for the isotropic Elastic Wave Propagator (SEM, 2nd order formulation).
Specifically, this PR follows the approach of [1], which is based on the satndard Fichtner approach [2] for frequency-indpendent quality factors.

This formulation is based on the standard-linear-solid (SLS) description. Compared to the referenced article, the method has been adapted to the SEM formulation with a Leap-Frog 2nd order timestepping scheme. Since the stress field $\sigma$ is not explicitly stored, this allows reducing the memory footprint of the auxiliary variables from $6 L$ to $3L + 3$, where $L$ si the number of auxiliary SLSs.

With the present formulation, the standard elastic wave equation is augmented as follows:

$\partial_t^2 u - \nabla\cdot C\varepsilon(u) + \sum\limits_{\nu=1}^L\psi^a_\nu = f$,
$\partial_t\psi^a_\nu + \omega_\nu\psi^a_\nu = \omega_\nu y_\nu\nabla\cdot C^a\varepsilon(u)$,

where $\varepsilon$ is the strain tensor, $C^a$ is the isotropic attenuative elasticity tensor (see [1]), and $(\omega_\nu, y_\nu)_{\nu=1}^L$ are a set of $L$ reference angular frequencies and adimensional anelasticity coefficients, respectively.

The attenuative matrix $C^a$ is computed from the attenuative Lamé coefficients $\lambda_a, \mu_a$, obtained from the original $\lambda, \mu$ Lamé fields using $\lambda_a+2\mu_a=(\lambda+2\mu)/Q_p$, $\mu_a=\mu/Q_s$ (see [1]). $Q_p$ and $Q_s$ are the pressure and shear elastic quality factor coefficients, which are cell-based fields required as an input and can be either defined on the mesh or given as constant fields via a FieldSpecification. A higher value of $Q$ indicates less attenuation.

Note that if $\sum\limits_{\nu=1}^L y_\nu > \mathrm{\min}(Q)$, some artifacts may appear in the solution, usually in the form of a long-term stable zero-velocity wave focused on the source point(s). If this happens, a warning is given by GEOS.

These coefficients $\omega_\nu$ and $y_\nu$ are needed to describe the behavior of $Q$ as a function of frequency, and can be computed as described in [2] given a minimal Q requirement as well as a frequency range. For example, for a frequency range $[9-110Hz]$ and $Q_{\mathrm{min}}=30$, a relatively flat $Q$ (independent of frequency) can be obtained using $\omega_1\approx 69 s^{-1}$, $\omega_2\approx 592 s^{-1}$, $y_1\approx 1.64$ and $y_2\approx 1.75$.

The SLS attenuation mechanism is activated by adding the option attenuationType="sls" to the ElasticWaveSEM solver. SLS parameters can be entered via the slsReferenceAngularFrequencies and slsAnelasticityCoefficients fields. If no SLS parameters are specified, the following default values are used:
$L=1$
$\omega_1=2\pi f$,
$y_1=\frac{2 Q_0}{Q_0-1}$,

where $f$ is the peak frequency of the source and $Q_0$ is the minimal value of $Q_p$ and $Q_s$ on the entire domain.

This mechanism works also with the VTI elastic wave propagator.

[1] P.-T. Trinh et al, "Efficient time-domain 3D elastic and viscoelastic full-waveform inversion using a spectral-element method on flexible Cartesian-based mesh", Geophysics, Vol. 84, No. 1.

[2] A. Fichtner, M. van Driel, "Models and Fréchet kernels for frequency-(in)dependent Q", Geophysical Journal International, Volume 198, Issue 3, September, 2014, Pages 1878–1889.

[codecov]

Codecov Report

Attention: Patch coverage is 97.04142% with 5 lines in your changes missing coverage. Please review.

Project coverage is 53.78%. Comparing base (19c3850) to head (021a2f5).

Files	Patch %	Lines
...sSolvers/wavePropagation/shared/WaveSolverBase.cpp	81.81%	4 Missing ⚠️
...econdOrderEqn/isotropic/ElasticWaveEquationSEM.cpp	98.27%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #3080      +/-   ##
===========================================
+ Coverage    53.70%   53.78%   +0.07%     
===========================================
  Files         1010     1011       +1     
  Lines        85639    85792     +153     
===========================================
+ Hits         45991    46139     +148     
- Misses       39648    39653       +5     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.