ViscoAcoustic

sweep.equations.ViscoAcoustic

ViscoAcoustic(
    spatial_order=4,
    device="cpu",
    backend="torch",
    dim=2,
    phase_shift=True,
    amplitude_damping=True,
)

Bases: sweep.equations.base.SecondOrderEquation

Second-order 2-D nearly constant-Q visco-acoustic wave equation.

An attenuating counterpart to :class:~sweep.equations.acoustic.Acoustic: the pressure-like field h1 is driven by vp**2 · Laplace(u) and absorbed by the same CPML formulation (cpmlr by default), with attenuation added on top as two decoupled terms:

• phase shift -- velocity dispersion, Zhu & Harris (2014) eq. 10's fractional Laplacian -c^2 eta_hat (-lap)^(gamma+1). Non-dissipative: the phase velocity follows Kjartansson's power law c_p = c0 (w/w0)^gamma (higher frequencies genuinely travel faster), so it moves the wavefront without removing energy. Implemented as the CPML FD Laplacian at the paper's velocity c = vp*cos(pi*gamma/2) plus a spectral remainder that vanishes identically as gamma -> 0.
• amplitude damping -- attenuation, the paper's tau d/dt (-lap)^(gamma+1/2) term. The dissipative term: a k^(2*gbar+1) wavenumber filter on du/dt removes energy, absorbing high frequencies more strongly (the medium acts as a distance-dependent low-pass). It weakens the wavefront without moving it.

The two switch independently via phase_shift / amplitude_damping (both default on). Both off reduces bit-exactly to Acoustic; both on gives the full visco-acoustic response. Heterogeneous media follow the paper's freezing-unfreezing treatment: every coefficient map varies in space, only the fractional EXPONENT is frozen at the average gbar (whose Q-derivative is also frozen, on every backend).

Because the decoupling is what makes the switches meaningful, note that the two effects are not physically independent -- causality ties dispersion and attenuation together through the Kramers-Kronig relations. The decoupled form is an approximation that buys the ability to treat them separately; the fully coupled constant-Q equation is more faithful but harder to solve.

Parameter order: vp, Q, omega. Wavefields: (h1, h2) plus the CPML auxiliaries. omega is the reference angular frequency that anchors the power law -- a real model parameter with a genuine gradient (it enters through (vp/omega)^(2*gamma)).

Backends: eager (CPU/GPU, torch + jax) and compiled CUDA (impl='c'). The CUDA path reuses the acoustic2d CPML kernels with the prepared vp_step and applies the amplitude damping per step via ATen/cuFFT; gradient memory strategies are full and ckpt (chunk or recursive). Boundary saving is not supported — the damping term is dissipative and global, so the reverse-time reconstruction that boundary saving relies on does not exist; the impl='c' default silently falls back to full and an explicit boundary request raises. Per-edge free surfaces work on both backends; topography and domain decomposition are eager-only.

References

Zhu, T. and Harris, J. M., 2014, Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians: Geophysics, 79(3), T105-T116, doi:10.1190/geo2013-0245.1 -- implements the paper's eq. 10 with coefficients eq. 11 (dispersion validated against the Kjartansson power law to 0.4% over the band; attenuation to 2%).

Models (constructor input order)

• vp (m/s): Visco-acoustic wave velocity model.
• Q: Quality factor controlling attenuation.
• omega (rad/s): Angular reference frequency for visco-acoustic attenuation.

Wavefields

• h1 (aliases: pressure, p): Primary visco-acoustic pressure-like wavefield; default source and receiver.
• h2 (aliases: pressure_prev): Previous-step pressure-like wavefield (internal).
• psix: CPML memory variable for the x-derivative term (internal).
• psiz: CPML memory variable for the z-derivative term (internal).
• zetax: CPML auxiliary wavefield for the x-direction update (internal).
• zetaz: CPML auxiliary wavefield for the z-direction update (internal).

Defaults

• source_type: ['h1']
• receiver_type: ['h1']
• pml_type: 'cpmlr'

Visco-acoustic wave equation solver.

Parameters:

• spatial_order (int, default: 4 ) –

  The order of the taylor expansion(Must be even). Defaults to 4.

• phase_shift (bool, default: True ) –

  Apply the (non-dissipative) dispersion term that shifts the wavefront. Defaults to True.

• amplitude_damping (bool, default: True ) –

  Apply the dissipative attenuation term that decays the wavefront. Defaults to True. Both off = acoustic; both on = full visco-acoustic.

C_NAME class-attribute

C_NAME = 'visco_acoustic2d'

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

prepare_models_for_c = True

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_batched_models class-attribute

supports_batched_models = True

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

supports_boundary_saving_c = False

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_image_topography class-attribute

supports_image_topography = True

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

supports_per_edge_free_surface = True

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

supports_per_edge_free_surface_c = True

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

c_eq_aux(prop)

ForwardInput.eq_aux for the CUDA kernels: the spectral filter grids on the runtime grid. The composition encodes the active terms (see visco_acoustic2d_spectral_grids): (D_loss,) damping only, (D_k2, D_frac) dispersion only, all three for both, () for none. |k| is raised to the frozen average exponent gbar from the latest :meth:prepare_models call (the paper's freezing-unfreezing: on impl='c' the exponent is DATA, matching the eager side where the exponent's gradient path is detached too). Built with the same :func:~sweep.equations.base.init_wavenumbers rule as the eager path and cached per (shape, h, device, gbar, switches).

init_abc

init_abc(type='cpml', **kwargs)

Set up the PML profiles, then the FFT wavenumber grid.

shape here is the runtime grid the step actually sees (PML pad plus stencil halo) and grid_spacing is a plain Python value, so the wavenumber grid is built entirely outside any jit trace.

prepare_models

prepare_models(models)

Map the user models (vp, Q, omega) onto the step coefficients (vp_step, B1, B2, A) for Zhu & Harris (2014) eq. 10 with coefficients eq. 11 (gamma = arctan(1/Q)/pi):

• vp_step = vp*cos(pi*gamma/2) -- the paper's wave-equation velocity c; the CPML FD Laplacian carries -c^2 k^2.
• B1 = c^2 and B2 = c^2 * eta_hat with eta_hat = cos(pi*gamma) * (vp/omega)^(2*gamma) -- the spectral remainder B1*k^2 - B2*k^(2*gbar+2) upgrades the FD term to the paper's dispersion operator -c^2 eta_hat k^(2*gamma+2).
• A = c^2 * tau_hat with tau_hat = (vp/omega)^(2*gamma) * sin(pi*gamma) / vp -- the loss coefficient for the k^(2*gbar+1) filter.

The spatially varying gamma enters every coefficient locally; only the fractional EXPONENT is frozen at the average gbar (stashed on self for :meth:_spectral_grids), which is the paper's freezing-unfreezing treatment of heterogeneous media. All ops are differentiable, so autograd routes the returned gradients back to vp / Q / omega.