Mondaic - Full Waveform Solutions

In the previous notebooks, we explored how to create and manipulate a wide variety of materials, including isotropic, anisotropic, and spatially varying materials, as well as how to orient them in space. Now, we will see how these materials can be used in actual wave propagation simulations, and what effects they have on the resulting wavefields.

This notebook demonstrates:

How to set up a Salvus simulation project using different material types.
How material symmetry, anisotropy and orientation affect wave propagation.
How to visualize the resulting wavefields using the supplied helper functions.

Note: While SalvusPy supports many material forms, the solver currently only accepts elastic.isotropic.Velocity, elastic.hexagonal.Velocity or elastic.triclinic.TensorComponents materials on the elastic side. You may need to convert your material to one of these forms before simulation.

We start by importing the necessary libraries and setting up some global parameters for the simulation. These include the domain size, frequency, and output settings. Feel free to adjust these parameters to experiment with different simulation setups.

import os
import numpy as np
import xarray as xr
from salvus import material
from salvus.material._details.material import to_solver_form
import salvus.namespace as sn
from scipy.spatial import Voronoi, cKDTree

# This is a local import for visualization helpers (provided with the notebook)
import wavefield_helpers

Here we define the main simulation parameters, such as the domain size, frequency, and output field. These will be used throughout the notebook.

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
SALVUS_FLOW_RANKS = int(os.environ.get("NUM_MPI_RANKS", 4))
OUTPUT_FIELD = "displacement"  # You can change this to e.g. 'velocity'
MAX_FREQ = 1.0  # Hz
PROJECT_DIR = f"project_{MAX_FREQ}"
DOMAIN_SIZE = 75_000  # meters
approximate_wavelength = 3500.0
RESOLUTION = int(
    7.5 * DOMAIN_SIZE / approximate_wavelength
)  # For visualization

We create a Salvus project with a cubic domain and define two events: one P-wave dominated and one S-wave dominated. This allows us to compare how different materials affect the propagation of different wave types.

# Create a cubic domain and initialize the project
box_domain = sn.domain.dim3.BoxDomain.from_bounds(
    DOMAIN_SIZE, DOMAIN_SIZE, DOMAIN_SIZE
)
p = sn.Project.from_domain(
    path=PROJECT_DIR, domain=box_domain, load_if_exists=True
)

# Define two events: P-wave and S-wave dominated
events = ["P-wave dominated", "S-wave dominated"]
location = {
    c: v for c, v in zip(["x", "y", "z"], box_domain.bounding_box.mean(axis=0))
}

# Define an isotropic moment tensor for the P-wave source
moment_P = {
    "mxx": 1.0,
    "myy": 1.0,
    "mzz": 1.0,
    "mxy": 0.0,
    "mxz": 0.0,
    "myz": 0.0,
}

source_P = sn.simple_config.source.cartesian.MomentTensorPoint3D(
    **location, **moment_P
)

# Define an antisymmetric gradient source (pure torque) for the S-wave source
source_S = sn.simple_config.source.cartesian.VectorGradientPoint3D(
    **location,
    gxx=0.0,
    gxy=1.0,
    gxz=0.0,
    gyx=-1.0,
    gyy=0.0,
    gyz=0.0,
    gzx=0.0,
    gzy=0.0,
    gzz=0.0,
)

sources = {k: v for k, v in zip(events, [source_P, source_S])}

# Add events to the project
for event in events:
    p += sn.Event(event_name=event, sources=[sources[event]])

We fix the simulation timing and output configuration to ensure results are comparable across different materials. The Ricker wavelet is used as the source time function.

stf = sn.simple_config.stf.Ricker(center_frequency=MAX_FREQ / 2.0)
number_of_snapshots = 1  # Only store the last time step
end_time = (DOMAIN_SIZE / 10_000) * (4.0 / 5.0)
start_time = stf.get_auto_start_time()
visualization_time = end_time
time_step = 2.5e-2 / MAX_FREQ
number_of_time_steps = (end_time - start_time) / time_step
elements_per_wavelength = 1.25

output_configuration = {
    "volume_data": {
        # Calculate interval as to store 1 snapshot (Salvus always stores first
        # and last).
        "sampling_interval_in_time_steps": int(
            number_of_time_steps / number_of_snapshots
        ),
        "fields": [OUTPUT_FIELD],
    },
}

wsc = sn.WaveformSimulationConfiguration(
    time_step_in_seconds=time_step,
    start_time_in_seconds=start_time,
    end_time_in_seconds=end_time,
)
ec = sn.EventConfiguration(
    wavelet=stf,
    waveform_simulation_configuration=wsc,
)
sc_common_settings = {
    "elements_per_wavelength": elements_per_wavelength,
    "event_configuration": ec,
    "max_frequency_in_hertz": MAX_FREQ,
}

Salvus supports a variety of ways to launch simulations. One would be to create meshes yourself (via e.g. the layered mesher) and add an UnstructuredMeshSimulationConfiguration object to your project. If you simply need a homogeneous simulation, you can use the simpler way of creating a ModelConfiguration object with a background model from a material, and passing this to a SimulationConfiguration object.

This function also contains an implementation of meshing a heterogeneous material, and creating an UnstructuredMeshSimulationConfiguration manually.

def create_results(project, name: str, material, mesh=None):
    """
    Collection of SalvusProject actions we need to repeat often in this
    notebook.

    Note that it uses some variables that are not passed explicitly (like
    `events`, `sc_common_settings`, `SALVUS_FLOW_SITE_NAME`, etc.). This works
    because Python creates a closure: the function “remembers” those
    outer-scope variables and can use them without you passing them in every
    time. A closure is simply a function that captures values from its defining
    scope.

    Only works for elastic materials, as it contains a material cast to
    triclinic.
    """

    if sn.material.is_oriented(material):
        # We need to convert the material to a tensor component form
        # before passing it to the background model.
        material = to_solver_form(material)

    if sn.material.is_homogeneous(material):
        # Homogeneous material, so SalvusProject will create the mesh for us.
        project.add_to_project(
            sn.SimulationConfiguration(
                name=name,
                model_configuration=sn.ModelConfiguration(
                    background_model=sn.model.background.homogeneous.FromMaterial(
                        material=sn.material.elastic.triclinic.TensorComponents.from_material(
                            material
                        )
                    )
                ),
                **sc_common_settings,
            ),
            overwrite=True,
        )
    else:
        # Heterogeneous material, so we create a mesh using the layered mesher
        # instead, and add it as an UnstructuredMeshSimulationConfiguration.
        mesh = sn.layered_meshing.mesh_from_domain(
            domain=box_domain,
            model=material,
            mesh_resolution=sn.MeshResolution(
                reference_frequency=MAX_FREQ,
                elements_per_wavelength=elements_per_wavelength,
            ),
        )

        project.add_to_project(
            sn.UnstructuredMeshSimulationConfiguration(
                name=name, unstructured_mesh=mesh, event_configuration=ec
            ),
            overwrite=True,
        )

    p.simulations.launch(
        simulation_configuration=name,
        events=events,
        site_name=SALVUS_FLOW_SITE_NAME,
        ranks_per_job=SALVUS_FLOW_RANKS,
        extra_output_configuration=output_configuration,
    )

    p.simulations.query(block=True, ping_interval_in_seconds=2.5)

We create 3 different materials:

material_1: A simple isotropic material with P and S wave velocities.
material_2: A transversely isotropic material with anisotropic P and S wave velocities.
material_3: A fully anisotropic material with a stiffness tensor defined by its stiffness coefficients. Those are chosen as to result in P and S wave velocities that are similar to the ones in the previous two materials while providing strong anisotropy -- the specific values are not meaningful.

material_1 = material.from_params(vp=5000, vs=3500, rho=2500)
material_2 = material.from_params(
    vpv=6000, vph=4000, vsv=2500, vsh=3500, rho=2500, eta=1.0
)
material_3 = material.from_params(
    **{
        # fmt: off
        "rho": 2500,

        "c11": 6.25e10,
        "c12": 1.25e9,
        "c13": 1.93e9,
        "c14": 0.7e9,
        "c15": 0.8e9,
        "c16": 0.9e9,

        "c22": 7.5625e10,
        "c23": 1.45e9,
        "c24": 0.8e9,
        "c25": 1.9e9,
        "c26": 1.0e9,
        
        "c33": 5.e10,
        "c34": 1.2e9,
        "c35": 1.9e9,
        "c36": 1.4e9,

        "c44": 2.7e10,
        "c45": 1.2e9,
        "c46": 2.2e9,
        
        "c55": 3.56e10,
        "c56": 1.6e9,

        "c66": 2.3e10,
        # fmt: on
    }
)

Running the simulations for the three materials is simplified by our helper function create_results. This function takes care of creating the simulations, launching them, and waiting for them to finish.

create_results(p, "sc_material_1", material_1)

[2025-05-21 02:43:42,162] INFO: Creating mesh. Hang on.
[2025-05-21 02:43:42,602] INFO: Submitting job array with 2 jobs ...

Because we record the volumetric vectorial displacement field, we need to visualize this field in a way that is meaningful on screens. To do this, we create 3 slices along the coordinate planes, and visualize the magnitude of the displacement field on these slices. Because our output is unstructured, we need to sample the wavefields to a regular grid before we can visualize them, on a grid defined by RESOLUTION.

_ = wavefield_helpers.visualize_magnitude(
    p,
    "sc_material_1",
    events,
    OUTPUT_FIELD,
    visualization_time,
    resolution=RESOLUTION,
)

For isotropic materials, we see that the P-wave dominated event (top row) is a pure P-wave, with the same distance travelled in all directions at a given time.

The S-wave event (bottom row) shows a clear radiation pattern, as the rotational motion in the X-Y plane can not radiate energy in the Z direction.

Now let's simulate our transversely isotropic material:

create_results(p, "sc_material_2", material_2)
_ = wavefield_helpers.visualize_magnitude(
    p,
    "sc_material_2",
    events,
    OUTPUT_FIELD,
    visualization_time,
    resolution=RESOLUTION,
)

[2025-05-21 02:45:29,660] INFO: Creating mesh. Hang on.
[2025-05-21 02:45:30,428] INFO: Submitting job array with 2 jobs ...

This material has a clear directionally dependent wave speed for both P and S waves in the X-Z and Y-Z plane. Additioanlly, the isotropic P-wave source mechanism, does generate some S-wave energy. An isotropic moment tensor source in a VTI medium does radiate shear waves, even though the moment tensor itself is purely volumetric.

Finally, we simulate our fully anisotropic material:

create_results(p, "sc_material_3", material_3)
_ = wavefield_helpers.visualize_magnitude(
    p,
    "sc_material_3",
    events,
    OUTPUT_FIELD,
    visualization_time,
    resolution=RESOLUTION,
)

[2025-05-21 02:49:21,100] INFO: Creating mesh. Hang on.
[2025-05-21 02:49:21,608] INFO: Submitting job array with 2 jobs ...

In this case, both source mechanisms generate phases of all kinds. Due to the heavily anisotropic nature of the material, phases split along many directions.

Alternatively to visualizing the magnitude of the displacement field, we can visualize the vector field itself. This is done by plotting the vector field on the slices. Note that this only plots the in-plane components of the vector field, and not the out-of-plane components.

You might notice that for the pure S-wave event, the X-Z and Y-Z slices are blank! The wave is polarized fully in the X-Y plane. Note that the Y-Z plane and X-Z still contain oscillations, but only in the X direction, and thus will not show up in in-plane the vector plots.

_ = wavefield_helpers.visualize_vector(
    p,
    "sc_material_1",
    events,
    OUTPUT_FIELD,
    visualization_time,
    resolution=RESOLUTION,
    scale=2.5e-4,
)

orientation_material_2 = sn.material.orientation.AzimuthDip.from_params(
    azimuth=90, dip=30
)

material_2_oriented = material_2.with_orientation(orientation_material_2)

# Create results will take care of casting the material to triclinic and
# "burning in" the orientation to material.
create_results(p, "sc_material_2_oriented", material_2_oriented)

[2025-05-21 02:51:59,057] INFO: Creating mesh. Hang on.
[2025-05-21 02:51:59,584] INFO: Submitting job array with 2 jobs ...

_ = wavefield_helpers.visualize_magnitude(
    p,
    "sc_material_2_oriented",
    events,
    OUTPUT_FIELD,
    visualization_time,
    resolution=RESOLUTION,
)

This visualization shows how in the Y-Z plane (i.e. azimuth = 90), the wavefields are rotated by 30 degrees. In the other two planes (X-Z and Y-Z), the wavefields are not rotated as such, but the intersection of the original wavefronts has shifted across the radiation pattern, so not exactly the same slice across the wavefront is visualized as in the unrotated case. E.g. in the left most plot on the bottom row, the two split S-waves are shown.

Lastly, because the S-wave source is anisotropic, it can be clearly noted how the source itself was not rotated, and therefor the magnitude across the radiation pattern has changed from the unrotated case.

As the grand finale of this notebook, we will create a strongly anisotropic material, which has spatially varying (heterogeneous) orientation. This is achieved by:

Discretizing the domain on an arbitrary grid.
Creating N random points in the domain, and using these as the centers of Voronoi cells.
For each Voronoi cell, we create a random local basis (orthonormal coordinate system) for the cell, and assign this to the cell.
For each grid point, we find the Voronoi cell it belongs to, and assign the local basis of that cell to the grid point.
Create DataArrays for the local basis vector components which carry the spatial information from the first discretization step.
We then create a DirectOrthonormalBasis object from the local basis vectors.
Finally, we create a new material with the same stiffness tensor as material_3, but with the new orientation.

# Create voronoi cells in 3D on a grid
# set random seed for reproducibility
np.random.seed(42)

grid_x_linear = np.linspace(0, DOMAIN_SIZE, RESOLUTION)
grid_y_linear = np.linspace(0, DOMAIN_SIZE, RESOLUTION)
grid_z_linear = np.linspace(0, DOMAIN_SIZE, RESOLUTION)
grid_x_mesh, grid_y_mesh, grid_z_mesh = np.meshgrid(
    grid_x_linear, grid_y_linear, grid_z_linear
)
grid_x = grid_x_mesh.flatten()
grid_y = grid_y_mesh.flatten()
grid_z = grid_z_mesh.flatten()

# Create random points in 3D
n_voronoi_points = 1500
points = np.random.rand(n_voronoi_points, 3) * DOMAIN_SIZE

# Get random local basis vectors for each point
v1 = np.random.randn(n_voronoi_points, 3)
v1 /= np.linalg.norm(v1, axis=1)[:, np.newaxis]
v2 = np.random.randn(n_voronoi_points, 3)
v2 /= np.linalg.norm(v2, axis=1)[:, np.newaxis]
v2 -= np.sum(v1 * v2, axis=1)[:, np.newaxis] * v1
v3 = np.cross(v1, v2)
v3 /= np.linalg.norm(v3, axis=1)[:, np.newaxis]

# Compute which grid point belongs to which voronoi cell
vor = Voronoi(points)
tree = cKDTree(points)
distances, indices = tree.query(np.c_[grid_x, grid_y, grid_z], k=1)

# Get the indices of the voronoi cells for each grid point
cell_indices = vor.point_region[indices]

# Ensure the max index is less than the number of voronoi points
cell_indices[cell_indices >= n_voronoi_points] = 0

# Create grid normal vectors

grid_v1 = np.zeros((RESOLUTION, RESOLUTION, RESOLUTION, 3))
grid_v2 = np.zeros((RESOLUTION, RESOLUTION, RESOLUTION, 3))
grid_v3 = np.zeros((RESOLUTION, RESOLUTION, RESOLUTION, 3))
for i in range(RESOLUTION):
    for j in range(RESOLUTION):
        for k in range(RESOLUTION):
            # Get the index of the cell this grid point belongs to
            index = cell_indices[
                i * RESOLUTION * RESOLUTION + j * RESOLUTION + k
            ]
            # Get the local basis vectors for this cell
            grid_v1[i, j, k] = v1[index]
            grid_v2[i, j, k] = v2[index]
            grid_v3[i, j, k] = v3[index]

Now we have the normal vectors assigned to each grid point, we need to pack these up into something which carries spatial information. For the Material module of Salvus, this is easiest done by creating DataArrays with the relevant dimensions and coordinates.

# Create xarray DataArrays for each basis vector component in a concise loop
basis_names = ["m00", "m01", "m02", "m10", "m11", "m12", "m20", "m21", "m22"]
basis_arrays = {}
for idx, (vec, name_prefix) in enumerate(
    zip([grid_v1, grid_v2, grid_v3], ["m0", "m1", "m2"])
):
    for comp in range(3):
        basis_arrays[f"{name_prefix}{comp}"] = xr.DataArray(
            vec[:, :, :, comp],
            dims=("x", "y", "z"),
            coords={
                "x": grid_x_linear,
                "y": grid_y_linear,
                "z": grid_z_linear,
            },
        )

We now have the basis vectors in a format that can be used to create a DirectOrthonormalBasis object. Behind the scenes, this class will also run Gram-Schmidt orthonormalization to ensure the basis vectors are orthonormal. One can also skip this step by setting skip_gram_schmidt=True, but this is not recommended unless you know that the basis vectors are definitely orthonormal.

# Create a DirectOrthonormalBasis object from the basis arrays
orientation_3d = sn.material.orientation.DirectOrthonormalBasis.from_params(
    **basis_arrays,
    skip_orthonormalization=False,
)

orientation_3d

<IPython.core.display.HTML object>

salvus.material.orientation.DirectOrthonormalBasis

Now we have the orientation object, we can create a new material with the same stiffness tensor as material_3, but with the new orientation. Salvus will automatically broadcast the spatial information in the orientation to the material itself, so we don't need to define a material for every grid point.

# Orient the material with the new orientation
material_3_3d_oriented = to_solver_form(
    material_3.with_orientation(orientation_3d)
)

Although create_results will take care of meshing the material, we can create a mesh for visualization purposes, to see exactly how the material is oriented in space.

# Create a dummy mesh for visualization (recreated in `create_results`)
mesh_3d_viz = sn.layered_meshing.mesh_from_domain(
    domain=box_domain,
    model=material_3_3d_oriented,
    mesh_resolution=sn.MeshResolution(reference_frequency=MAX_FREQ),
)

# Retain only C11 for visualization by removing all other elemental fields
for field in list(mesh_3d_viz.elemental_fields.keys()):
    if field != "C11":
        mesh_3d_viz.elemental_fields.pop(field)

mesh_3d_viz

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

That looks definitely like Voronoi grains!

Let's simulate and visualize the wavefield for this material.

create_results(p, "sc_material_3_oriented", material_3_3d_oriented)

[2025-05-21 02:56:12,872] INFO: Submitting job array with 2 jobs ...

wavefield_helpers.visualize_magnitude(
    p,
    "sc_material_3_oriented",
    events,
    OUTPUT_FIELD,
    visualization_time,
    resolution=RESOLUTION,
)

The resulting wavefields show strong coda from all the changes in material properties, but on a macro-level, the wavefront more resemble isotropic behavior, even when passing through an anisotropic medium. Effectively, this simulation is a demonstration what happens when propagating long wavelengths through a medium with a lot of small-scale anisotropy: in the (long-wavelength) limit, the medium behaves isotropically.

Part 3 - Oriented materials

Receiver Gathers

© Mondaic AG (2025)

Site Map Contact Us

Academic License Agreement

PAGE CONTENTS