ViscoElastic¶
sweep.equations.ViscoElastic ¶
Bases: sweep.equations.base.FirstOrderEquation
First-order 2-D visco-elastic wave equation (GSLS, nearly constant Q).
Velocity-stress staggered grid as :class:~sweep.equations.elastic.Elastic
plus 3 L memory variables (r_xx_l, r_zz_l, r_xz_l, one
triple per standard linear solid). With relaxed P and shear moduli
pi_R, mu_R, strengths tp_l = tau_l(Qp), ts_l = tau_l(Qs)
and theta = dvx/dx + dvz/dz::
dsxx/dt = pi_U theta - 2 mu_U dvz/dz + sum_l r_xx_l
dszz/dt = pi_U theta - 2 mu_U dvx/dx + sum_l r_zz_l
dsxz/dt = mu_U (dvx/dz + dvz/dx) + sum_l r_xz_l
dr_xx_l/dt = -(r_xx_l + pi_R tp_l theta - 2 mu_R ts_l dvz/dz) / tau_sigma_l
(r_zz_l, r_xz_l analogous)
with unrelaxed moduli pi_U = pi_R (1 + sum_l tp_l), mu_U = mu_R (1 +
sum_l ts_l), i.e. M(w) = M_R (1 + sum_l tau_l i w tau_sigma_l / (1 + i w
tau_sigma_l)). The memory variables live at integer time levels and are
advanced with the trapezoidal rule (unconditionally stable for any
dt / tau_sigma_l); the stress picks up their time average.
Models are vp, vs, rho, Qp, Qs; vp/vs are phase velocities
at the reference frequency f_ref (SPECFEM's
READ_VELOCITIES_AT_f0 = .true.), converted exactly to the relaxed
moduli through the GSLS complex modulus at f_ref. All five models are
differentiable. Qp = Qs = inf reduces bit-exactly to Elastic.
The tau_sigma_l are fixed by (n_sls, f_band) at construction
(relaxation frequencies log-spaced over the band) and are not model
parameters; tau_l(Q) is the least-squares constant-Q fit (see
:class:GSLSFit).
Stability: the explicit stress update runs at the UNRELAXED
velocities, which exceed vp/vs (e.g. ~8% at Q = 10 for a 12x band),
so choose dt from the CFL condition on
sqrt((lam_U + 2 mu_U) / rho) (:meth:unrelaxed_velocities). Q must
exceed self.q_min (about 1.3 for the default 3 mechanisms), below which
the fit would need a negative relaxation strength; n_sls >= 5 is
refused for the same reason.
Free surface (flat, per edge): the surface-row normal strain rate is
solved so that the normal traction stays exactly zero INCLUDING the memory
variables; in the elastic limit this is Robertsson's modified coefficient
4 mu (lam + mu) / (lam + 2 mu) used by Elastic.
Backends: eager and compiled CUDA (impl='c'; forward, and the
full, chunk- and recursive-checkpoint backwards with a hand-derived
exact adjoint). Not supported: boundary saving (the attenuation is
dissipative, so reverse reconstruction is unstable; impl='c' defaults
to full), topography / APM, domain decomposition, per-shot batched
models.
References
Robertsson, J. O., Blanch, J. O. and Symes, W. W., 1994, Viscoelastic finite-difference modeling: Geophysics, 59(9), 1444-1456, doi:10.1190/1.1443701 -- the staggered-grid memory-variable scheme. Emmerich, H. and Korn, M., 1987, Incorporation of attenuation into time-domain computations of seismic wave fields: Geophysics, 52(9), 1252-1264, doi:10.1190/1.1442386 -- the constant-Q fit (linear least squares for the strengths at fixed relaxation frequencies).
Models (constructor input order)
vp(m/s): P-wave phase velocity at f_ref.vs(m/s): S-wave phase velocity at f_ref.rho(kg/m^3): Density model.Qp: P-wave quality factor (inf = no attenuation).Qs: S-wave quality factor (inf = no attenuation).
Wavefields
vx(aliases:velocity_x): Particle velocity in the x direction; default receiver.vz(aliases:velocity_z): Particle velocity in the z direction; default receiver.sxx(aliases:stress_xx): Normal stress in the x direction; default source.szz(aliases:stress_zz): Normal stress in the z direction; default source.sxz(aliases:stress_xz,shear_xz): Shear stress component.r_xx_0: GSLS memory variable xx, mechanism 0 (internal).r_zz_0: GSLS memory variable zz, mechanism 0 (internal).r_xz_0: GSLS memory variable xz, mechanism 0 (internal).r_xx_1: GSLS memory variable xx, mechanism 1 (internal).r_zz_1: GSLS memory variable zz, mechanism 1 (internal).r_xz_1: GSLS memory variable xz, mechanism 1 (internal).r_xx_2: GSLS memory variable xx, mechanism 2 (internal).r_zz_2: GSLS memory variable zz, mechanism 2 (internal).r_xz_2: GSLS memory variable xz, mechanism 2 (internal).m_vxx: CPML memory variable for dvx/dx (internal).m_vxz: CPML memory variable for dvx/dz (internal).m_vzx: CPML memory variable for dvz/dx (internal).m_vzz: CPML memory variable for dvz/dz (internal).m_txxx: CPML memory variable for dsxx/dx (internal).m_txxz: Reserved elastic auxiliary field (internal).m_tzzx: Reserved elastic auxiliary field (internal).m_tzzz: CPML memory variable for dszz/dz (internal).m_txzx: CPML memory variable for dsxz/dx (internal).m_txzz: CPML memory variable for dsxz/dz (internal).
Defaults
source_type:['sxx', 'szz']receiver_type:['vx', 'vz']pml_type:'cpmls'
Args:
spatial_order: FD order (as Elastic).
device: Device for the gradient kernels.
backend: 'torch' only.
f_ref: Reference frequency (Hz) at which vp / vs are given.
f_band: (f_min, f_max) over which Q is fitted constant. Defaults
to SPECFEM's band: a factor 12 centred (geometrically) on f_ref.
For FWI set it to the inverted band.
n_sls: Number of standard linear solids (SPECFEM default 3).
C_NAME
class-attribute
¶
str(object='') -> str str(bytes_or_buffer[, encoding[, errors]]) -> str
Create a new string object from the given object. If encoding or errors is specified, then the object must expose a data buffer that will be decoded using the given encoding and error handler. Otherwise, returns the result of object.str() (if defined) or repr(object). encoding defaults to sys.getdefaultencoding(). errors defaults to 'strict'.
prepare_models_for_c
class-attribute
¶
bool(x) -> bool
Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.
supports_boundary_saving_c
class-attribute
¶
bool(x) -> bool
Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.
supports_per_edge_free_surface
class-attribute
¶
bool(x) -> bool
Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.
supports_per_edge_free_surface_c
class-attribute
¶
bool(x) -> bool
Returns True when the argument x is true, False otherwise. The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, and cannot be subclassed.
c_eq_aux ¶
eq_aux for the CUDA kernels: one float32 CPU tensor
[a_0 .. a_{L-1}, c_0 .. c_{L-1}], evaluated exactly as the eager
step evaluates them (float32 dt).
prepare_models ¶
(vp, vs, rho, Qp, Qs) -> [lam_U, mu_U, rho, lam_U + 2 mu_U,
pi_R tp_0 .. pi_R tp_{L-1}, mu_R ts_0 .. mu_R ts_{L-1}] -- the step
coefficients shared by the eager step and the CUDA kernels (rho third,
where elastic2d's kernels keep it). The
unrelaxed Lame set is exactly Elastic._lame of the unrelaxed
velocities, so at Q = inf (tau_l = 0 -> factors exactly 1) the
elastic part is bit-identical to Elastic.
impl='c' hands (B, 1, nz, nx) repeats of one shared model (this
equation has no per-shot models): the coefficients are computed once
and repeated, and the repeat's backward sums the per-shot gradients.
relaxation_coefficients ¶
Trapezoidal memory-update weights r+ = a_l r - c_l F:
a_l = (1 - dt/2ts) / (1 + dt/2ts), c_l = (dt/ts) / (1 + dt/2ts).
dt may be the eager step's float32 0-d tensor (possibly inside a
compiled step: kept symbolic); :meth:c_eq_aux evaluates the same
expressions so both backends see identical float32 weights.
unrelaxed_velocities ¶
(vp_U, vs_U) for (vp, vs, rho, Qp, Qs): the velocities the
explicit stress update actually runs at (above vp / vs for
finite Q, equal at Q = inf) -- take dt from the CFL condition
on vp_U.
sweep.equations.visco_elastic.GSLSFit ¶
Per-mechanism strengths tau_l(Q) for fixed tau_sigma_l.
With a_l(w) = (w ts_l)^2 / (1 + (w ts_l)^2) and b_l(w) = w ts_l /
(1 + (w ts_l)^2), M(w)/M_R = 1 + sum_l tau_l (a_l + i b_l). Constant
Q means sum_l tau_l (Q b_l - a_l) = 1 on the band; the least-squares
tau solves (BB - (BA + AB) q + AA q^2) tau = q (sb - sa q) with
q = 1/Q (sampled at n_freq log-spaced frequencies). By Cramer's
rule tau_l = q N_l(q) / D(q) with polynomial D (degree 2L) and
N_l; :meth:tau evaluates that form elementwise (q = 0 gives
tau = 0 exactly), :meth:tau_np the direct solve (reference).
min_valid_q ¶
Smallest Q above which every tau_l(Q) > 0 (a physical,
dissipative relaxation spectrum). Raises if even the high-Q limit
has a non-positive strength (too many mechanisms for the band).
tau ¶
Torch: Q of any shape -> list of L tensors tau_l (same
shape/dtype as Q), differentiable in Q.
q_of_frequency ¶
Exact Q(f) realised for a target Q (numpy; checks the fit).