Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

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.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 0x79deb35cd110>
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 0x79dea4c8f690>
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-11-20 13:21:34,938] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2411201321032069_27c7eb2b00@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 0x79dea3ad3650>
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 0x79dea375df10>
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-11-20 13:21:37,432] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2411201321435915_6f230d0e05@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 0x79dee78b9e50>
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 0x79dea34323d0>
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 0x79dea3774f90>
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-11-20 13:21:38,521] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2411201321524066_805564425b@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 0x79de9d90c690>
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 0x79dea03bad90>
As show in the 2-D example, shotgathers or other point measurements can also be extracted at arbitrary locations and curves.
Next
GETTING STARTED
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