Mondaic - Full Waveform Solutions

This page here has been created for the latest stable release of Salvus. You have chosen to view the documentation for another Salvus version. Please be aware there might be some small differences you have to account for.

Determining the physical units of inputs and outputs in Salvus requires understanding which interpretation of the governing equations is chosen. Salvus supports several physically consistent interpretations for the scalar acoustic potential in purely acoustic simulations, and a unique one when acoustic and elastic physics are coupled. This page lists all supported interpretations, how to translate saved fields to physically meaningful quantities (pressure, particle velocity, displacement), and the special role of the coupled case.

In the elastic wave equation, all outputs are given in their canonical SI units. This is in contrast to the acoustic case, where the units depend on the chosen interpretation of the scalar field, typically some potential. Note that in coupled acoustic-elastic simulations, only one interpretation is allowed.

The strong form we solve is written in a unified way as

\rho^{-1} c^{-2} \partial_t^2 \phi = \nabla \cdot (\rho^{-1} \nabla \phi) + f,

with density $\rho$ , wave speed $c$ , scalar potential $\phi$ , and source term $f$ . Different physical choices for the dimensions assigned to $\phi$ and $f$ lead to different practical post-processing steps. Salvus does not enforce a single convention for pure acoustic runs; you decide which table below applies.

In the first interpretation, the source term represents a volumetric injection rate with units $[f]=\text{s}^{-1}$ , one finds $[\phi]=\text{Pa}\,s$ . To get pressure and particle velocity, use:

Desired output	Expression	Field to save	Further operations
Pressure (Pa)	$\partial_t \phi$	`["phi_t"]`	None
Particle velocity ( $\text{m}\,\text{s}^{-1}$ )	$\rho^{-1} \nabla \phi$	`["gradient-of-phi"]`	Multiply by inverse density

Assume $f$ carries divergence of force density units $[f]=\text{N}\,\text{m}^{-4}$ , we find $[\phi]=\text{kg}^2\,\text{m}^{-4}\,\text{s}^{-2}$ . Then

Desired output	Expression	Field to save	Further operations
Pressure (Pa)	$\rho^{-1} \phi$	`["phi"]`	Multiply by inverse density
Particle velocity ( $\text{m}\,\text{s}^{-1}$ )	$\rho^{-2} \int \nabla \phi\, dt$	`["gradient-of-phi"]`	Time integration + multiply by $\rho^{-2}$

Motivated by interface coupling (see below) we can assign $[\phi]=\text{kg}/\text{m}$ . Plugging this into the same governing equation implies a dimensionless right hand side in the strong form (i.e., $[f]=1$ ). This means that the source term becomes a normalized monopole source strength. For more insight in this, see the section on conversions below. Practical consequences:

Desired output	Expression	Field to save	Further operations
Pressure (Pa)	$-\partial_t^2 \phi$	`["phi_tt"]`	Multiply by -1
Normal displacement (at interface) ( $\text{m}$ )	$\rho^{-1} \partial_n \phi$	`["gradient-of-phi"]`	Select normal component & multiply by $1/\rho$
Particle displacement (bulk) ( $\text{m}$ )	$\rho^{-1} \nabla \phi$	`["gradient-of-phi"]`	Multiply by $1/\rho$
Particle velocity (bulk) ( $\text{m}$ )	$\partial_t \rho^{-1} \nabla \phi$	`["gradient-of-phi"]`	Multiply by $1/\rho$ , then take the time derivative*

*Note that one could also output [phi_t] and then compute the gradient afterwards, but this is much more involved on unstructured meshes than taking the time derivative of the gradient field with e.g. a central difference (higher order) scheme.

Interpretation C is mandatory when coupling to elasticity; A and B are only valid in purely acoustic scenarios.

Choosing among A or B (pure acoustic): Pick the one aligned with how your source is specified physically. Need immediate pressure? Use A. Have a known force density history? Use B. Planning to share meshes with elastic physics? Use C from the start to avoid reinterpretation.

For full algorithmic details see our paper here. Additionally, a complete list of fields which can be output can be found in the documentation.

It is common to have a source specified in one physical description and wish to implement it under a different interpretation (or eventually migrate to a coupled simulation). The relations below let you translate between a purely acoustic physical source description and the dimensionless right-hand side f used in the coupled convention (Interpretation C). Denote the coupled-form source term by f_c(t) (dimensionless) and the time variable by t.

Given one of the following physical source descriptions in a pure acoustic setting:

Physical description	Symbol & units	Conversion to dimensionless `f_c(t)` (Interpretation C)	Notes
Volume injection rate	`q(t)` [s^-1]	`f_c(t) = - ∫ q(t) dt` (constant chosen so f_c→0 outside source)	Time integral; remove DC after integration
Pressure history (localized)	`p(t)` [Pa]	`f_c(t) = - p(t) / (ρ c^2)`	Uses local acoustic impedance (`ρ c^2`) scaling
Divergence of force density	`g(t)` [N m^-4]	`f_c(t) = - (1/ρ) ∫∫ g(t) dt^2`	Two time integrals; detrend afterwards

Each physical description corresponds to a way of injecting work or mass/volume change into the acoustic field. In the coupled Interpretation C, the governing equation has a dimensionless right-hand side, so any dimensional source must be rescaled to produce the correctly scaled potential $\phi$ with units $\text{kg}/\text{m}$ . A volume injection rate $q(t)$ is (up to constants) the time derivative of a dimensionless monopole strength—hence we integrate once to obtain its cumulative effect. A localized pressure history $p(t)$ enters as an imposed normal stress; dividing by the bulk modulus proxy $\rho c^2$ yields a dimensionless amplitude consistent with $p = -\phi_{tt}$ . Finally, a divergence of force density $g(t)$ is (in the pure-acoustic interpretations) the second time derivative of a monopole-like loading; recovering the dimensionless f_c therefore requires two time integrations (and detrending to remove integration constants). This mapping ensures that after substitution into the coupled equation, $\phi$ and its derivatives reproduce the originally intended physical pressure / traction and volume-change effects.

Conversely, once you have a dimensionless f_c(t) (coupled interpretation):

Quantity (pure acoustic meaning)	Recover from `f_c(t)`	Comment
Volume injection rate `q(t)`	`q(t) = - d f_c / dt`	Central differences recommended
Local pressure `p(t)`	`p(t) = - ρ c^2 f_c(t)`	Near-source approximation
Divergence force density `g(t)`	`g(t) = - ρ d^2 f_c / dt^2`	Central differences recommended

Implementation tips:

Always window the source time function so its value and its first derivative vanish near simulation start to avoid long-period artifacts.
Numerical differentiation/integration should use the native solver time step to avoid interpolation phase errors (or heavily supersampled signals).
After integration, remove any residual linear trend (detrending) to enforce zero mean; this prevents persistent offsets in φ.

Design guidance:

Use Interpretation A if your physical source is naturally specified as a volume injection rate and you mainly care about pressure time series.
Use Interpretation B if you have (or conceptually think in terms of) a known divergence of force density field.
Use Interpretation C if you will ever couple to elasticity, or if you prefer pressure via a second derivative and a straightforward formula for normal displacement at interfaces.
Converting previously archived A/B datasets to the coupled view can be done via the conversion tables by reconstructing f_c(t) and recomputing derived fields.

In the elastic case, Salvus solves the general elastic wave equation. For the analysis of units involved in the source term specifically, let's limit ourselves to purely elastic linear media. The strong form of this equation is:

\rho \, \partial _t^2 \mathbf{u} = \nabla \cdot (\mathbf{C} : \nabla \mathbf{u}) + \mathbf{\mathcal{f}}.

Here $\rho$ is density, $\mathbf{u}$ is displacement, and $f$ represents an external forcing term. The time and space dependency of these last three terms is taken as implicit to simplify notation. $\mathbf{C}$ is a fourth-order tensor characterizing the stiffness of the medium, and the symbol $:$ denotes a contraction over adjacent indices.

In the elastic case, it is a little less flexible to substitute units of the applied force. Let's see what units are present in the equation when one defines $\mathbf{u}$ as displacement in meters:

Symbol	Quantity (Unit)
$\rho$	Density ( $\text{kg}\cdot\text{m}^{-3}$ )
$\mathbf{u}$	Displacement ( $\text{m}$ )
$\mathbf{C}$	Elasticity ( $\text{Pa}$ )
$\partial_t$	Time derivative, variation per second ( $\text{s}^{-1}$ )
$\nabla$	Spatial derivative, variation per meter ( $\text{m}^{-1}$ )

Let's substitute these units into the equation, and figure out in what units the source term should be. In this exercise, we forego the dimensionalities of the quantities.

[\text{kg} \cdot \text{m}^{-3}] \, \partial _t^2 [\text{m}] = \nabla \cdot ([\text{Pa}] : \nabla [\text{m}]) + \mathbf{\mathcal{f}}.

Applying the space and time derivatives as $[\text{m}^{-1}]$ and $[\text{s}^{-1}]$ respectively, and dropping the vector and tensor operators as they don't influence the units:

[\text{kg} \cdot \text{m}^{-3}] \, [\text{m} \cdot \text{s}^{-2}] = [\text{m}^{-1}] \left( [\text{Pa}] [\text{m} \cdot \text{m}^{-1}]\right) + \mathbf{\mathcal{f}}.

Simplifying and using $[\text{Pa}] = [\text{N} \cdot \text{m}^{-2}]$ , where $\text{N}$ denotes force in units of Newtons gives:

[\text{kg}\cdot\text{m} \cdot \text{s}^{-2} \cdot \text{m}^{-3}] = [\text{m}^{-1}] \cdot [\text{N} \cdot \text{m}^{-2}] + \mathbf{\mathcal{f}},

So, finally using $[ \text{N} ] = [\text{kg} \cdot \text{m} \cdot \text{s}^{-2}]$ :

[\text{N} \cdot \text{m}^{-3}] = [\text{N} \cdot \text{m}^{-3}] + \mathbf{\mathcal{f}}.

Which dictates that our source term $f$ needs to be interpreted in the strong form as a source term with units $\text{N} \cdot \text{m}^{-3}$ .

One often desires point sources to be injected in a wavefield simulation. As a case study, we look at how to interpret the units of VectorPoint3D. The signature of creating this object is as follows:

Copy
class VectorPoint3D(salvus.flow.simple_config.source._BaseSource):
    def __init__(
        self,
        x: float,
        y: float,
        z: float,
        fx: float,
        fy: float,
        fz: float,
        source_time_function: Optional[
            salvus.flow.simple_config.stf._Base
        ] = None,
        rotation_on_input: Optional[Dict] = None,
    ):
        ...

To know in what units the source_time_function $\text{stf}(t)$ of these sources should be interpreted, we start with the formulation of this point source $f$ as injected into the strong form:

f(\mathbf{x}, t) = \text{stf}(t) \, \mathbf{F} \,\delta(\mathbf{x})

In this definition, the vector $\mathbf{F} = [$ fx, fy, fz $]$ allows one to orient (and scale) the source in any direction. The symbol $\delta$ indicates the (multidimensional) Dirac delta function.

The units of the Dirac delta function should be interpreted to collapse to dimensionless if integrated over its argument. In this case, we would integrate over three spatial coordinates, giving it the units $\text{m}^{-3}$ . If now we interpret the orientation vector $\mathbf{F}$ as dimensionless, the units of $\text{stf}(t)$ are found in the following way:

[\text{N} \cdot \text{m}^{-3}] = \text{stf}(t) \cdot [ 1 ] \cdot [\text{m}^{-3}]

Yielding that $\text{stf}(t)$ should be interpreted as Newtons. Conversely, if $\mathbf{F}$ is not unit length, only the combination $\text{stf}(t) \, \mathbf{F}$ should be interpreted as the exact Newtons injected in the wavefield simulation.

When an acoustic region adjoins an elastic region within the same simulation, the interface conditions uniquely determine the permissible scaling of the acoustic potential. Across any acoustic–elastic boundary with outward normal $\hat{n}$ (Nissen-Meyer et al., 2007, equation 46):

u_n = \hat{n}\cdot \mathbf{u} = \rho^{-1} \partial_n \phi, \qquad \sigma_{nn} = -p = \partial_t^2 \phi.

With elastic displacement $[\mathbf{u}] = \text{m}$ and stress $[\sigma_{nn}] = \text{Pa}$ , these two relations force $[\phi]=\text{kg}/\text{m}$ and $p = -\phi_{tt}$ . Thus only Interpretation C is consistent; A and B would break continuity. Practical recipe for coupled simulations:

Request phi_tt to obtain pressure: pressure = -phi_tt.
Request gradient-of-phi and multiply its components by 1/rho to get displacement components.
All elastic outputs retain their canonical SI units: displacement (m), velocity (m/s), stress / pressure (Pa), etc.

Why is there a minus sign in $p = -\phi_{tt}$ ? The sign follows directly from the stress balance with our positive pressure convention; in the frequency domain $\phi_{tt} \rightarrow -\omega^2 \phi$ giving the expected 180° phase relative to $\phi$ .

What are the units of the acoustic source term now? Dimensionless in the strong form; any physical forcing is mediated through traction continuity or is scaled so that resulting $\phi$ has $[\text{kg}/\text{m}]$ .