Working with Salvus’ wavefield output files

Salvus contains a powerful set of tools that enable users to output simulated dynamic quantities, such as velocity, displacement, or strain, within volumes and along surfaces. These outputs are, in general, saved on unstructured meshes, making them somewhat awkward to work with using standard tools such as numpy or xarray. To remedy this, Salvus contains a lightweight WavefieldOutput class that enables the simple extraction of stored quantities to regular grids. Let’s dive in and take a look.

import os

import numpy as np
import pathlib
import salvus.namespace as sn
from salvus.mesh.unstructured_mesh_utils import read_model_from_h5
from salvus.toolbox.helpers import wavefield_output

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

Setting up a simple simulation

We’ll demonstrate the basic concepts with a 2-D isotropic elastic mesh created using Salvus’ layered_meshing API. The use of this API, however, is not required to use the functionality presented herein, and any volumetric simulation output from a Salvus simulation will suffice.

domain_2d = sn.domain.dim2.BoxDomain(x0=0, x1=1, y0=0, y1=1)

material = sn.layered_meshing.material.elastic.Velocity.from_params(
    rho=1.0, vp=1.0, vs=np.sqrt(1 / 3)
)

mesh_2d = sn.layered_meshing.mesh_from_domain(
    domain=domain_2d,
    model=material,
    mesh_resolution=sn.MeshResolution(
        reference_frequency=10.0, elements_per_wavelength=2.0
    ),
)

mesh_2d

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f7d284b0bd0>

We’ll then create a project …

p_2d = sn.Project.from_domain(
    path="project_2d", domain=domain_2d, load_if_exists=True
)

… with source at the top of the domain …

src_2d = sn.simple_config.source.cartesian.VectorPoint2D(
    x=0.5,
    y=1.0,
    fx=0.0,
    fy=-1.0,
)

p_2d.add_to_project(sn.Event(event_name="event", sources=[src_2d]))

… and set up a simple simulation.

p_2d.add_to_project(
    sn.UnstructuredMeshSimulationConfiguration(
        name="sim_2d",
        unstructured_mesh=mesh_2d,
        event_configuration=sn.EventConfiguration(
            wavelet=sn.simple_config.stf.Ricker(center_frequency=2.5),
            waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
                end_time_in_seconds=2.0
            ),
        ),
    ),
)
p_2d.viz.nb.simulation_setup("sim_2d", events=["event"])

<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7f7c153e0c90>

Finally, let’s run the simulation with additional outputs. Let’s save velocity, and sample the output every 10 time steps.

p_2d.simulations.launch(
    simulation_configuration="sim_2d",
    events=p_2d.events.list(),
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
    extra_output_configuration={
        "volume_data": {
            "sampling_interval_in_time_steps": 10,
            "fields": ["velocity"],
        },
    },
)
p_2d.simulations.query(block=True)

[2026-10-12 20:35:45,147] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2610122035238303_424713ec83@local

True

Analyzing the wavefield output

We can use Salvus’ wavefield_output module to work with the volumetric data that we just output.

out_2d = p_2d.simulations.get_simulation_output_directory("sim_2d", "event")

vel_wo = wavefield_output.WavefieldOutput.from_file(
    pathlib.Path(
        out_2d,
        "volume_data_output.h5",
    ),
    "velocity",
    "volume",
)

The simplest thing we might want to do is resample the wavefield to a regular grid in order to visualize the wavefield directly in the notebook. To do this we can use the wavefield_output_to_xarray function which, as one might expect, returns the wavefield interpolated onto an xarray.DataArray. Below we pass a list of length 2 containing points we want to extract to in each dimension. Passing a list in this manner will result in an extraction to a rectilinear grid spanned by the cartesian product of the input arrays, similar to the behavior of np.meshgrid.

vel_2d = wavefield_output.wavefield_output_to_xarray(
    vel_wo,
    points=[np.linspace(0, 1, 101), np.linspace(0, 1, 101)],
)

That’s all there is to it. We can now slice vel along arbitrary axes in space and time using xarray’s powerful indexing tools, and plot the result directly in the notebook. Below we plot the y- (second) component of velocity at the 35th output slice.

vel_2d.isel(c=1, t=35).T.plot(shading="gouraud", infer_intervals=False)

<matplotlib.collections.QuadMesh at 0x7f7c15220e50>

In addition to outputting values on rectilinear grids, we can also pass a n-dimensional array of points to get “shotgathers” with complex geometries. For instance, below well get a shotgather that emulates receiver locations following a figure eight.

a, b = 0.5, 0.5
u = np.linspace(0, 2 * np.pi, 1001)
x = a * np.sin(u) / (1 + np.cos(u) ** 2) + 0.5
y = b * np.sin(u) * np.cos(u) / (1 + np.cos(u) ** 2) + 0.5

vel_sg8 = wavefield_output.wavefield_output_to_xarray(
    vel_wo, points=np.column_stack((x, y))
)

vel_sg8.isel(c=1).plot(
    shading="gouraud", infer_intervals=False, yincrease=False
)

<matplotlib.collections.QuadMesh at 0x7f7c15141190>

Output data at final time

Saving the last snapshot of a simulation can be interesting for several use cases such as checkpointing or restarting a simulation with inhomogeneous initial values.

Salvus supports outputting the final time values on an unstructured mesh. The output format is independent of the partition and ranks used in the simulation.

To demonstrate its use, we add another simulation and reduce the end time yo 0.5 seconds.

p_2d.add_to_project(
    sn.UnstructuredMeshSimulationConfiguration(
        name="sim_2d_short",
        unstructured_mesh=mesh_2d,
        event_configuration=sn.EventConfiguration(
            wavelet=sn.simple_config.stf.Ricker(center_frequency=2.5),
            waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
                end_time_in_seconds=0.5
            ),
        ),
    ),
)

Similar to volumetric output, final time data can be selected through the extra_output_configuration.

p_2d.simulations.launch(
    simulation_configuration="sim_2d_short",
    events=p_2d.events.list(),
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
    extra_output_configuration={
        "final_time_data": {"fields": ["velocity"]},
    },
)
p_2d.simulations.query(block=True)

[2026-10-12 20:35:47,866] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2610122035869737_f78156b7ab@local

True

To visualize the results, we can simply load the data and attach them to the mesh. Note that the output uses the polynomial degree of the wavefield (4) and we thus have to make sure the mesh has the same tensor order if we want to visualize the fields.

mesh_out = mesh_2d.copy()
mesh_out.change_tensor_order(4)

out_2d = p_2d.simulations.get_simulation_output_directory(
    "sim_2d_short", "event"
)
final_output = read_model_from_h5(out_2d / "final_time_data_output.h5")

mesh_out.elemental_fields.clear()
mesh_out.attach_field("velocity_Y", final_output["velocity_Y"])

mesh_out

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f7c14fb7010>

3-D simulations

In 3-D, the concepts are the same. Let’s again create a mesh using the layered_meshing API.

domain_3d = sn.domain.dim3.BoxDomain(x0=0, x1=1, y0=0, y1=1, z0=0, z1=1)
mesh_3d = sn.layered_meshing.mesh_from_domain(
    domain=domain_3d,
    model=material,
    mesh_resolution=sn.MeshResolution(
        reference_frequency=5.0, elements_per_wavelength=2.0
    ),
)

mesh_3d

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f7c15200050>

We’ll again create a porject with a source at the top of the domain …

p_3d = sn.Project.from_domain(
    path="project_3d", domain=domain_3d, load_if_exists=True
)

src_3d = sn.simple_config.source.cartesian.VectorPoint3D(
    x=0.5,
    y=0.5,
    z=1.0,
    fx=0.0,
    fy=0.0,
    fz=-1.0,
)


p_3d.add_to_project(sn.Event(event_name="event", sources=[src_3d]))

… and set up a simple simulation.

p_3d.add_to_project(
    sn.UnstructuredMeshSimulationConfiguration(
        name="sim_3d",
        unstructured_mesh=mesh_3d,
        event_configuration=sn.EventConfiguration(
            wavelet=sn.simple_config.stf.Ricker(center_frequency=2.5),
            waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
                end_time_in_seconds=1.0
            ),
        ),
    ),
)
p_3d.viz.nb.simulation_setup("sim_3d", events=["event"])

<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7f7c14e58750>

Let’s again add some volumetric output, but this time also save the wavefield along the top mesh surface.

p_3d.simulations.launch(
    simulation_configuration="sim_3d",
    events=p_3d.events.list(),
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=1,
    extra_output_configuration={
        "volume_data": {
            "sampling_interval_in_time_steps": 10,
            "fields": ["velocity"],
        },
        "surface_data": {
            "sampling_interval_in_time_steps": 10,
            "fields": ["velocity"],
            "side_sets": ["z1"],
        },
    },
)
p_3d.simulations.query(block=True)

[2026-10-12 20:35:48,810] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2610122035813210_2a0e16de4f@local

True

Analyzing wavefield output

As before we can use the WavefieldOutput class to help us manipulate the data we just created. Here we create 2 objects: one for volume and one for surface data.

output_dir = p_3d.simulations.get_simulation_output_directory(
    "sim_3d", "event"
)

vel_3d_vol_wo = wavefield_output.WavefieldOutput.from_file(
    pathlib.Path(output_dir, "volume_data_output.h5"), "velocity", "volume"
)

vel_3d_srf_wo = wavefield_output.WavefieldOutput.from_file(
    pathlib.Path(output_dir, "surface_data_output.h5"), "velocity", "surface"
)

Let’s take a look at the surface output first. Since the surface is a true 3-D surface it contains 3 spatial coordinate values, rather than 2 as you might expect. This is important as the surface we are analyzing may not be axis aligned, and in fact may not be planar at all. To emulate a truly 2-d surface we can use the drop_coordinate method to drop a coordinate from the output. Since we know we stored the surface output along the "z1", or top, side set, we can drop the z, or last, dimension to pull out a surface.

vel_3d_srf_wo_xy = vel_3d_srf_wo.drop_dimension(2)

We can then extract the data values to a rectilinear grid and plot using xarray.

vel_3d_srf = wavefield_output.wavefield_output_to_xarray(
    vel_3d_srf_wo_xy, [np.linspace(0, 1, 101)] * 2
)

vel_3d_srf.isel(c=2, t=25).plot(shading="gouraud", infer_intervals=False)

<matplotlib.collections.QuadMesh at 0x7f7c101f2950>

We can of alternatively the same result by just extracting the 3-D output and along the top surface of the mesh.

vel_3d_vol = wavefield_output.wavefield_output_to_xarray(
    vel_3d_vol_wo, [np.linspace(0, 1, 101)] * 2 + [1.0]
)

As we expect, the values are the same.

vel_3d_vol.isel(c=2, t=25).sel(z=1.0).plot(
    shading="gouraud", infer_intervals=False
)

<matplotlib.collections.QuadMesh at 0x7f7c112c9790>

As show in the 2-D example, shotgathers or other point measurements can also be extracted at arbitrary locations and curves.