Mondaic Documentation
www.mondaic.com
Mondaic
Version:
No results found.

Units in Coupled Simulations

This notebook is a compact, hands-on companion to the longer units explainer in our knowledge base. In purely acoustic simulations Salvus deliberately permits several physically consistent interpretations of the scalar potential.
Once we consider coupling to elasticity, the two interface conditions strip away that freedom: only Interpretation C on the explainer page survives because it is the unique scaling making both dimensionally consistent. This forces [ϕ]=kg/m[\phi]=\text{kg}/\text{m} and consequently pressure is recovered as p=−∂t2ϕp = -\partial_t^2 \phi without extra scaling.

In this notebook, you will:

  1. Build the simplest geometry that exercises sign conventions: an inclined interface (a horizontal one hides normal component bookkeeping).
  2. Run a coupled acoustic–elastic simulation with a single moment tensor source in the elastic half.
  3. Extract receiver time series exactly on the interface in both media.
  4. Demonstrate the two continuity conditions numerically and discuss the role of density.
With a slope the global X/Y components project non-trivially onto the interface normal, so we explicitly test that our vector/tensor projections (and sign flips) are correct—mirroring real post-processing where geometric errors easily creep in.
We will validate the two interface conditions:
  • Pressure / traction continuity: p  n^a→e=(σ⋅n^e→a)p \; \hat{n}_{a \rightarrow e} = (\boldsymbol{\sigma}\cdot\hat{n}_{e \rightarrow a}).
  • Normal displacement continuity: n^⋅ue=ρ−1∂n^ϕ\hat{n}\cdot \mathbf{u}_e = \rho^{-1} \partial_{\hat{n}} \phi, or in terms of velocity n^⋅ve=ρ−1∂n^∂tϕ\hat{n}\cdot \mathbf{v}_e = \rho^{-1} \partial_{\hat{n}} \partial_t \phi.
In these equations n^\hat{n} is a unit vector normal to the coupled interface, and the subscripts indicate which direction the normal is pointing. For the normal displacement condition, the normal is arbitrarily but consistently (same on both sides) oriented.
From Salvus, we will directly output phi_tt (∂t2ϕ\partial_t^2 \phi) and gradient-of-phi (∇ϕ\nabla \phi, which allows us to compute ∂n^ϕ\partial_{\hat{n}} \phi by ∂n^ϕ=n^⋅∇ϕ\partial_{\hat{n}} \phi = \hat{n} \cdot \nabla \phi). Reconstructing either via finite differencing (time or space) introduces phase and amplitude error—especially on coarser meshes or when sampling intervals get sparse. Saving them directly from SalvusCompute avoids numerical differentiation noise and simplifies unit reasoning.
Copy
import os

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import salvus.namespace as sn

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")

1. Overview & Geometry

We intentionally choose a modest domain and an inclined acoustic–elastic contact to force non-trivial projections. A horizontal interface would make the normal purely vertical and hide potential sign or component-mapping mistakes. Here, both X and Y contribute to the normal.
Let's create a simple 2D box domain.
domain = sn.domain.dim2.BoxDomain.from_bounds(10, 10)
Now let's define a sloping interface that separates our acoustic and elastic regions, and its related unit vectors.
We will have our acoustic region overlying our elastic region, s.t. we have a definition we can use to flip our normal vectors consistently.
# Chosen to give an irrational slope; avoids round-off coincidences.
interface_angle_deg = 20.251432
slope_coefficient = np.tan(np.radians(interface_angle_deg))

x = np.linspace(0.0, 10.0, 1000)
interface_topography = slope_coefficient * (x - 5.0) + 5.0

# We'll create the unit vectors as xarray DataArrays, which will simplify the
# later vector/tensor math.
vector_components = ["X", "Y"]

unit_vectors_global = [
    np.array([1.0, 0.0]),
    np.array([0.0, 1.0]),
]

# Normal vector pointing from acoustic to elastic on surface at incline_angle
normal_vector_ac_to_el = xr.DataArray(
    np.sin(np.radians(interface_angle_deg)) * unit_vectors_global[0]
    - np.cos(np.radians(interface_angle_deg)) * unit_vectors_global[1],
    dims=["component"],
    coords={"component": vector_components},
)

normal_vector_el_to_ac = -normal_vector_ac_to_el

normal_vector_ac_to_el  # Pointing down-ish, normal to surface
<xarray.DataArray (component: 2)> Size: 16B
array([ 0.34614051, -0.93818269])
Coordinates:
  * component  (component) <U1 8B 'X' 'Y'
<xarray.DataArray (component: 2)> Size: 16B
array([ 0.34614051, -0.93818269])
Coordinates:
  * component  (component) <U1 8B 'X' 'Y'
Let's visualize the domain, interface and normal vectors to make sure we understand the geometry.
fig, ax = plt.subplots(
    figsize=(6, 6),
)
subsampling = 100
ss_slice = slice(int(subsampling / 2), None, subsampling)

ax.fill(
    domain.bounding_box[[0, 2, 3, 1], 0],
    domain.bounding_box[[0, 2, 3, 1], 1],
    color="lightgray",
    edgecolor="black",
)


ax.scatter(
    x[ss_slice],
    interface_topography[ss_slice],
    color="black",
    s=2.5,
)

ax.quiver(
    x[ss_slice],
    interface_topography[ss_slice],
    np.ones_like(interface_topography)[ss_slice]
    * normal_vector_ac_to_el.values[0],
    np.ones_like(interface_topography)[ss_slice]
    * normal_vector_ac_to_el.values[1],
    color="#1f77b4",
    scale=15,
    width=0.005,
    label="Acoustic-to-Elastic normals",
)
ax.quiver(
    x[ss_slice],
    interface_topography[ss_slice],
    np.ones_like(interface_topography)[ss_slice]
    * normal_vector_el_to_ac.values[0],
    np.ones_like(interface_topography)[ss_slice]
    * normal_vector_el_to_ac.values[1],
    color="#d62728",
    scale=15,
    width=0.005,
    label="Elastic-to-Acoustic normals",
)

interface = plt.plot(
    x, interface_topography, "--", color="black", lw=1.5, label="Interface"
)

# Labels for acoustic and elastic regions
ax.text(1.0, 7.5, "Acoustic", fontsize=14, color="#1f77b4")
ax.text(1.0, 1.5, "Elastic", fontsize=14, color="#d62728")

ax.legend(loc="upper left", framealpha=0)

ax.set_xlabel("X [m]")
ax.set_ylabel("Y [m]")
ax.set_xlim(0, 10)
ax.set_ylim(0, 10)

ax.set_title("Domain with inclined acoustic–elastic interface")

ax.set_aspect("equal")
Now we'll use the domain and interface to create a mesh of an acoustic domain overlying an elastic domain.
We set contrasting densities to make the scaling in the displacement continuity condition visible: the acoustic normal displacement involves 1/ρacoustic1/\rho_{\text{acoustic}}.
# The layered mesher takes in DataArrays as well
interface_curve = xr.DataArray(
    data=interface_topography, dims=["x"], coords={"x": x}
)

density_acoustic = 1600.0
density_elastic = 2250.0

velocity_acoustic = 800.0
velocity_elastic_p = 2500.0
velocity_elastic_s = 1450.0

meshing_resolved_freq = 2000.0

model = [
    sn.material.acoustic.isotropic.Velocity.from_params(
        rho=density_acoustic, vp=velocity_acoustic
    ),
    sn.layered_meshing.interface.Curve(interface_curve),
    sn.material.elastic.isotropic.Velocity.from_params(
        rho=density_elastic, vp=velocity_elastic_p, vs=velocity_elastic_s
    ),
]

mr = sn.MeshResolution(
    reference_frequency=meshing_resolved_freq,
    elements_per_wavelength=2.0,
    model_order=2,
)

mesh = sn.layered_meshing.mesh_from_domain(
    domain=domain, model=model, mesh_resolution=mr
)

mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x75c9ea7a44d0>
We place a moment tensor source in the elastic half, and two receivers on the interface, one in each medium. We will output a comprehensive set of fields for exposition, but we only need a few to compute our interface conditions.
elastic_fields = [
    "displacement",
    "gradient-of-displacement",
    "velocity",
    "acceleration",
    "strain",
    "stress",
]
acoustic_fields = [
    "phi",
    "phi_t",
    "phi_tt",
    "gradient-of-phi",
]

# Receivers on the interface, one in each medium
interface_receiver_acoustic = sn.simple_config.receiver.cartesian.Point2D(
    x=5.0,
    y=5.0,
    station_code="R2A",
    fields=acoustic_fields,
)
interface_receiver_elastic = sn.simple_config.receiver.cartesian.Point2D(
    x=5.0,
    y=5.0,
    station_code="R2E",
    fields=elastic_fields,
)

source_elastic = sn.simple_config.source.cartesian.MomentTensorPoint2D(
    x=5, y=2.5, mxx=1.0, myy=1, mxy=0
)

event_1 = sn.Event(
    event_name="event_1",
    sources=source_elastic,
    receivers=[
        interface_receiver_elastic,
        interface_receiver_acoustic,
    ],
)
Next, we define our temporal parameters. We use a Ricker wavelet with half the center frequency of our meshing resolution. The end time is chosen s.t. a few reflections fit in the domain. We use a typical simulation configuration from this.
source_time_function = sn.simple_config.source.stf.Ricker(
    center_frequency=meshing_resolved_freq / 2.0
)

end_time = 1e-2

wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=end_time)

ec = sn.EventConfiguration(
    wavelet=source_time_function, waveform_simulation_configuration=wsc
)

sc = sn.UnstructuredMeshSimulationConfiguration(
    name="coupled_medium",
    unstructured_mesh=mesh,
    event_configuration=ec,
)
Finally, we create a Salvus Project, add our simulation configuration and event, and inspect it.
project_path = "project_units"

p = sn.Project.from_domain(
    path=project_path,
    domain=domain,
    load_if_exists=True,
)
p += sc
p += event_1

p
<salvus.project.project.Project at 0x75ca06683e50>
We can now launch the simulation. It should take only a few seconds on any machine.
p.simulations.launch(
    simulation_configuration=sc.name,
    events=[event_1],
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
)

p.simulations.query(block=True, ping_interval_in_seconds=1)
[2025-11-12 21:44:18,325] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2511122144411548_53c4e4c947@local

True
We now analyze both interface equalities.
# Some common setup:
# pick contrasting, readable colors
color_acoustic = "#1f77b4"  # blue
color_elastic = "#d62728"  # red

# Define colors and styles
style = {
    "acoustic": dict(color=color_acoustic, lw=1.9, ls="-"),
    "elastic": dict(color=color_elastic, lw=1.5, ls="--"),
}
# Let's get the data:
event_data = p.waveforms.get(data_name=sc.name, events=[event_1])[0]
Compute the pressure and its side of the interface condition.
pressure = -(  # Note the minus sign!
    event_data.get_waveform_data_xarray(receiver_field="phi_tt")
    .unstack()
    .sel(receiver_id="XX.R2A.", component="A")
)

pressure_interface_condition = pressure * normal_vector_ac_to_el
Now compute the elastic side of the interface condition, the normal stress vector t=σ⋅n^e→a\mathbf{t} = \boldsymbol{\sigma}\cdot \hat{n}_{e\rightarrow a}.
data_elastic = (
    event_data.get_waveform_data_xarray(receiver_field="stress").unstack()
).sel(receiver_id="XX.R2E.")

# Stress component are in the resulting data structure in the following order:
# {"0": "X", "1": "Y", "2": "XY"}

# Compute components of stress tensor \cdot normal_vector
nx = normal_vector_el_to_ac.sel(component="X")
ny = normal_vector_el_to_ac.sel(component="Y")
s_xx = data_elastic.sel(component="0")
s_yy = data_elastic.sel(component="1")
s_xy = data_elastic.sel(component="2")

# traction vector on the solid: t = σ n
t_x = s_xx * nx + s_xy * ny
t_y = s_xy * nx + s_yy * ny

normal_stress = xr.concat(
    [t_x, t_y],
    dim=xr.DataArray(["X", "Y"], dims="component"),
).assign_coords(component=["X", "Y"])
Finally, we can make a combined plot of both sides of the interface condition.
# Combine both fields into one DataArray for easier plotting
interface_conditions = xr.concat(
    [pressure_interface_condition, normal_stress],
    dim="field",
).assign_coords(field=["acoustic", "elastic"])

pressure_ifcond_label = "p \; \hat{n}_{a\\rightarrow e}"
stress_ifcond_label = "\sigma \cdot \hat{n}_{e\\rightarrow a}"

labels = {
    "acoustic": f"acoustic pressure (${pressure_ifcond_label}$)",
    "elastic": f"elastic normal traction (${stress_ifcond_label}$)",
}

# Manual loop by field (hue)
fig, axes = plt.subplots(
    nrows=len(interface_conditions.component),
    sharex=True,
    figsize=(7, 2.5 * len(interface_conditions.component)),
)

for i, comp in enumerate(interface_conditions.component.values):
    ax = axes[i] if len(interface_conditions.component) > 1 else axes
    for field_name, attrs in style.items():
        interface_conditions.sel(component=comp, field=field_name).plot(
            ax=ax, label=labels[field_name], **attrs
        )
    ax.set_title(f"component = {comp}")
    ax.set_ylabel("Pressure / Stress [Pa]")
    ax.legend(frameon=False, fontsize=9)
    ax.grid(alpha=0.3, ls=":")

axes[-1].set_xlabel("Time [s]")
fig.suptitle(
    f"Normal stress interface condition (${pressure_ifcond_label} = "
    f"{stress_ifcond_label}$)",
)
fig.tight_layout()
That looks pretty good!
Now let's verify the normal displacement continuity condition n^⋅ue=ρ−1∂nϕ\hat{n}\cdot \mathbf{u}_e = \rho^{-1} \partial_n \phi.
unit_vector_displacement = (
    # Would also work with the reverse unit_vector_el_to_ac, it's a matter of
    # sign of the direction of the displacement/velocity
    normal_vector_ac_to_el
)

# Get displacement in the global coordinate system in the elastic medium
global_displacement_elastic = (
    event_data.get_waveform_data_xarray(
        receiver_field="displacement"
    ).unstack()
).sel(receiver_id="XX.R2E.")


# Get the displacement in the global coordinate system in the acoustic medium
global_displacement_acoustic = (1.0 / density_acoustic) * (
    event_data.get_waveform_data_xarray(receiver_field="gradient-of-phi")
    .unstack()
    .sel(receiver_id="XX.R2A.")
)

# Get the normal components
normal_displacement_elastic = xr.dot(
    global_displacement_elastic, unit_vector_displacement
)
normal_displacement_acoustic = xr.dot(
    global_displacement_acoustic, unit_vector_displacement
)
Ready to plot both sides of the normal displacement:
fig, ax = plt.subplots(figsize=(7, 3))
normal_displacement_acoustic.plot(
    label="acoustic normal displacement",
    color=style["acoustic"]["color"],
    lw=style["acoustic"]["lw"],
    linestyle=style["acoustic"]["ls"],
)
normal_displacement_elastic.plot(
    linestyle=style["elastic"]["ls"],
    label="elastic normal displacement",
    color=style["elastic"]["color"],
    lw=style["elastic"]["lw"],
)
plt.legend()
plt.xlabel("Time [s]")
plt.ylabel("Normal displacement [m]")
plt.title(
    "Continuous displacement interface condition ($\\hat{n} \\cdot "
    "\\mathbf{u}_e = \\rho^{-1} \\partial_{\\hat{n}} \\phi$)"
)
Text(0.5, 1.0, 'Continuous displacement interface condition ($\\hat{n} \\cdot \\mathbf{u}_e = \\rho^{-1} \\partial_{\\hat{n}} \\phi$)')
Spot on!
In this notebook we have:
  • Verified pressure/traction continuity via p=−∂t2ϕp = - \partial_t^2 \phi.
  • Verified normal displacement continuity via density‑scaled gradient of ϕ\phi.
  • Built understanding of the units of the scalar potential in coupled simulations.
Previous
Acoustic VTI
Next
2D teleseismic simulation
We use cookies to analyze our traffic, but won't track you until explicit consent is given. Using our contact form will create relevant cookies and transmit data to third party services to keep track of your request. Read more about how we use personal information in our Privacy Policy.
PAGE CONTENTS