Material models and physical scope¶
This page states the constitutive assumptions behind the propagators. A tensor with the right shape and finite entries is not necessarily a physical material. Deepwave validates the family conditions below when ordinary model values are readable, before propagation or planning, and never projects a model into an admissible class. Under vmap or compilation, a tensor’s values may be unavailable to these checks; validity remains the caller’s precondition while static shape/device checks still apply. A user-supplied max_vel does not validate stiffnesses or density. See Units and source conventions for sources, receivers, and their time and spatial normalisation, and Discretisation: the exact scheme for the numerical scheme.
Variables and dimensions¶
Use one consistent length, time, and mass system. In SI:
Symbol / argument | Meaning | SI unit |
|---|---|---|
\(c\), | Wave speed; elastic P and S wave speeds | \(\mathrm{m/s}\) |
\(\mathbf v\) | Particle velocity, distinct from the argument | \(\mathrm{m/s}\) |
\(\rho\), | Mass density | \(\mathrm{kg/m^3}\) |
\(b\), | Inverse density \(1/\rho\) in material; zero is a separate vacuum sentinel | \(\mathrm{m^3/kg}\) |
\(K,\lambda,\mu,c_{IJ}\) | Bulk modulus, Lamé parameters, and Voigt stiffness components | \(\mathrm{Pa}\) |
\(p,\sigma_{ij}\) | Compression-positive pressure and tension-positive stress | \(\mathrm{Pa}\) |
\(\epsilon_{ij}\) | Infinitesimal strain | Dimensionless |
\(Q,\epsilon,\delta,\gamma\) (Thomsen parameters) | Attenuation quality factor and anisotropy ratios, respectively | Dimensionless |
\(\tau_\sigma,\tau_\varepsilon\) | Stress and strain relaxation times | \(\mathrm{s}\) |
| Cyclic frequency; angular frequency is \(\omega=2\pi f\) | \(\mathrm{s}^{-1}\) (Hz) |
The scalar wavefield \(u\) has the units chosen for its interpretation; when interpreted as pressure it has units Pa and its PDE source has units \(\mathrm{Pa/m^2}\). Born perturbations have the units of their named background parameter: dv is velocity, drho density, dlamb/dmu pressure, and dbuoyancy inverse density. A Born perturbation need not itself be a positive material: it is a direction about a valid background.
Scalar and variable-density acoustic media¶
scalar solves \(\nabla^2u-u_{tt}/c^2=f\). Positive \(c\) gives the ordinary wave equation; spatial variation in \(c\) alone does not implement variable-density acoustics. acoustic instead solves
with positive density and bulk modulus in material. Here \(\mathbf f\) is force density (\(\mathrm{N/m^3}\)), and \(s\) is volume injection rate per volume (\(\mathrm{s}^{-1}\)). Zero-density/stiffness masks used for a surface are separate degeneracies; they are not ordinary material with an infinite or negative inverse density.
The numerical scalar update uses \(u_{tt}=c^2(\nabla^2u-f)\). It accepts either sign of \(c\) and the zero-speed limit. The acoustic and viscoacoustic speed inputs likewise enter through \(c^2\); their density and attenuation restrictions are separate. The CFL bound uses absolute speed. For the Born propagators, reversing both v and dv preserves the linearisation of \(c^2\).
A zero time step evaluates the discrete update at dt=0 without CFL substeps. In particular, the scalar two-level recurrence becomes \(u^{n+1}=2u^n-u^{n-1}\). Unequal initial time levels therefore still produce a changing sequence; a zero time step does not discard that initial state. With no supplied PML frequency, the zero-step limit uses zero rather than inferring a spectrum from samples with zero temporal spacing.
Isotropic elastic solids, fluids, and dimension reduction¶
The constitutive law is \(\sigma=\lambda\operatorname{tr}(\epsilon)I+2\mu\epsilon\). Writing \(\epsilon=\epsilon_{\rm dev}+\operatorname{tr}(\epsilon)I/3\) gives the 3D strain energy density
A positive-energy 3D solid therefore has \(\rho>0\), \(\mu>0\), and \(K_3>0\). The Lamé parameter \(\lambda\) may be negative. A fluid has \(\mu=0\), \(K_3>0\); its zero shear energy is deliberate. An elastic vacuum cell is represented by zero stiffnesses and zero buoyancy, not by taking \(\rho\to\infty\). The improved-vacuum treatment follows Zeng et al. (2012); its surface averaging is a distinct part of the scheme, with the restrictions described in Elastic wave equation and Boundary conditions.
In 2D, Deepwave keeps the in-plane P-SV displacement/velocity and stress components: the omitted displacement and its spatial derivatives vanish. This is a plane-strain restriction, not a plane-stress reduction of the Lamé parameters. The retained strain energy is positive for \(\mu>0\) and \(\lambda+\mu>0\). These weaker conditions suffice for the restricted 2D system; a model intended as a physical 3D material should satisfy the full 3D conditions above. In 1D only the longitudinal modulus \(M=\lambda+2\mu>0\) enters. The elastic p receiver is minus the mean of the retained normal stresses; in 2D that is not the 3D thermodynamic pressure. The source-conversion derivation in Units and source conventions keeps this distinction.
VTI material domains¶
The symmetry axis is dimension 0: y in 2D and z in 3D. The complete stiffness matrix and exact Thomsen conversion are given in Elastic VTI (anisotropic) wave equation; the parameter definitions originate in Thomsen (1986). Although the conventional parameter names come from a weak-anisotropy paper, the implemented stiffness conversion uses the exact defining algebra, not its first-order expansion.
For 2D elastic VTI, the retained normal-stiffness block is \(\left[\begin{smallmatrix}c_{11}&c_{13}\\c_{13}&c_{33}\end{smallmatrix}\right]\) and the shear stiffness is \(c_{55}\). Positive definiteness is equivalent to
In 3D, \(c_{12}=c_{11}-2c_{66}\). In an orthonormal basis containing the horizontal difference strain, the normal block splits into the scalar \(2c_{66}\) and
Thus a positive-energy 3D VTI solid has
These are positive-energy material conditions. They should not be confused with a complete classification of strong ellipticity or a fully discrete CFL condition. Fluid and vacuum limits are semidefinite degeneracies, not strictly positive-definite solids. As in isotropic elasticity, the 2D conditions only constrain the represented in-plane system.
acoustic_vti sets the shear stiffnesses to zero and keeps two stress-like fields. Its physical domain has \(b>0\), \(c_{11},c_{33}>0\), and \(c_{11}c_{33}\ge c_{13}^2\); equality includes elliptical anisotropy. The propagator checks these conditions when values can be inspected, using a scaled normal-block comparison that avoids squaring large stiffnesses. This is the pseudo-acoustic approximation of Duveneck et al. (2008), not a general anisotropic fluid: its extra pseudo-shear mode and the model-changing source_taper option are described in Pseudo-acoustic VTI wave equation.
What the interior energy argument establishes¶
The aligned velocity/stress layout follows the staggered-grid approach of Virieux (1986); a separate, inspectable implementation of this layout is the PETSc DMStag elastic-wave example. Deepwave’s relevant implementation choices are:
Normal stresses, normal strains, and their stiffness block are collocated. The supplied normal stiffnesses are used there without averaging.
A shear modulus at a staggered shear node is the harmonic mean of the four surrounding cell values, \(4/\sum_{j=1}^{4}(1/\mu_j)\). Elastic VTI uses \(c_{55}\) or \(c_{66}\) according to the shear plane. The positive-material branch has positive means; zero/vacuum masks use a separate zero branch. Values near numerical zero are guarded in the code.
Velocity-node density is the arithmetic mean of its two adjacent cell densities, so its buoyancy is \(2/(\rho_i+\rho_{i+1})\), not the arithmetic mean of their buoyancies. Elastic inputs are converted from buoyancy before this averaging.
These statements correspond to elastic.prepare_parameters, elastic_vti.prepare_parameters, and _staggered_grid.stagger_density_arithmetic_mean. Harmonic modulus and arithmetic density averaging are discussed in Moczo et al. (2002); this citation provides context, not a claim that Deepwave implements every averaging detail of that paper.
For the lossless interior on an infinite uniform grid (or with vanishing boundary terms), let \(D\) form normal strain rates and engineering shear strain rates from velocity. Its entries have units inverse length; engineering shear rate is \(\partial_i v_j+\partial_j v_i\). The paired integer/half-grid stencils in staggered_grid.h and _staggered_grid.py obey summation by parts: their reverse shifted difference is the negative transpose. This follows term by term from the antisymmetric difference weights at each supported spatial order, before applying boundary fills.
Let \(R_\rho\) be the diagonal matrix of positive staggered densities, and \(C\) the block-diagonal matrix of collocated normal stiffnesses and staggered shear stiffnesses. With sources and PML omitted, the semidiscrete equations and energy are
All terms inside the parentheses have energy-density units Pa; \(E_h\) is energy in 3D and energy per omitted-axis length in 2D. For positive-definite \(C\),
Spatially varying positive material coefficients do not break this cancellation. This is an interior, continuous-time discrete-space result. It does not prove stability of every finite timestep, PML, ghost-cell closure, or staircase vacuum interface. The inverse-stiffness expression also does not apply directly to fluid/vacuum degeneracies. Finite boundaries must supply the appropriate energy flux or a compatible closure; Boundary conditions states the accepted combinations and the tests of their reflection and energy behaviour. CPML/MPML adds memory and damping terms, requiring its own analysis. In particular this derivation does not guarantee arbitrary anisotropic PML stability or a particular interface-reflection error.
Attenuating media¶
viscoacoustic uses \(K=\rho c^2\) as the relaxed, zero-frequency bulk modulus. For one SLS, finite positive \(Q\) and positive reference frequency give positive relaxation times with \(\tau_\varepsilon\ge\tau_\sigma\). The source convention and exact sourced equations are stated in Units and source conventions; the homogeneous constitutive response and memory update are derived in Viscoacoustic (attenuating) wave equation.
With several mechanisms, \(\tau_l=w_l/Q\) is a weak-attenuation band fit. Relaxation-memory approaches to approximate constant Q are established in Blanch, Robertsson and Symes (1995). The following formula specifies Deepwave’s particular fit, rather than claiming identity with every parameter-fitting method in that literature. For \(x_l=\omega\tau_{\sigma,l}\), its normalised complex modulus is
Here the harmonic convention is \(e^{i\omega t}\); positive imaginary modulus denotes loss. The fit approximates the denominator by \(1/Q\) and treats the numerator as approximately one. It therefore does not enforce the supplied target Q exactly at low Q. Individual weights may be signed: aggregate loss, rather than the sign of one weight alone, determines whether the response is dissipative at a frequency. Large cancellation can impair numerical accuracy, particularly for many mechanisms in a narrow band. Neither the fit nor the elastic interior-energy argument proves passivity for every possible fitted configuration; the documented approximation should be checked against the attenuation required for the scientific application.
Positive infinite q represents an exactly lossless cell for both the single- and multiple-mechanism attenuation models. It may be mixed with finite positive q in the same model. The attenuation strengths are proportional to 1/q, so the existing coefficient formulas give zero strength and the limiting zero derivative with respect to q there.