Viscoelastic wave propagator

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:

$\partialt^2 u - \nabla\cdot C\varepsilon(u) + \sum\limits{\nu=1}^L\psi^a_\nu = f$, $\partialt\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.

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

GEOS-DEV / GEOS

Viscoelastic wave propagator #3080

Codecov Report