Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

This documentation is not for the latest stable Salvus version.

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.
Copy
import os

import numpy as np
import pathlib
import salvus.namespace as sn
from salvus.mesh import layered_meshing as lm
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 = lm.material.elastic.Velocity.from_params(
    rho=1.0, vp=1.0, vs=np.sqrt(1 / 3)
)

mesh_2d = lm.mesh_from_domain(
    domain_2d,
    material,
    sn.MeshResolution(reference_frequency=10.0, elements_per_wavelength=2.0),
)

mesh_2d
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fb6588214d0>
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 0x7fb6903186d0>
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)
[2024-03-15 09:23:39,189] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150923295463_a6b2b7ad8d@local

True
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 0x7fb64d09d850>
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 0x7fb64d0f3cd0>
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)
[2024-03-15 09:23:41,901] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150923904585_e6ad6d03f2@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 = {}
mesh_out.attach_field("velocity_Y", final_output["velocity_Y"])

mesh_out
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fb64ccfd950>
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 = lm.mesh_from_domain(
    domain_3d,
    material,
    sn.MeshResolution(reference_frequency=5.0, elements_per_wavelength=2.0),
)

mesh_3d
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fb64d2c02d0>
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 0x7fb64cbfff50>
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)
[2024-03-15 09:23:43,151] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2403150923154805_15a9154071@local

True
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 0x7fb60342e510>
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 0x7fb603b41290>
As show in the 2-D example, shotgathers or other point measurements can also be extracted at arbitrary locations and curves.
Previous
Receiver Gathers
Next
Part 1 - Setup
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