Mondaic - Full Waveform Solutions

This tutorial assumes that you've already read the introduction and 2D meshing tutorial for this external mesh series.

First, the base geometry consisting of a cuboid and a hemisphere will be generated using

brick x 2 y 2 z 1
create sphere radius 0.5
move Volume 2 z 0.5 include_merged

In order to remove the part of the cuboid where it overlaps with the sphere, a boolean subtraction can be performed with

subtract volume 2 from volume 1 keep_tool

Next, the sphere can be webcut in order to remove the upper portion of it; this leaves behind the bottom hemisphere.

webcut volume 2  with plane from surface 9
delete volume 2

Next, overlapping surfaces are merged in order to consolidate common surface entities using

merge surface all

Before commencing with the meshing of the geometry, a series of webcuts need to be applied to the geometry. These will be used in order to guide the mesher when constructing the discretization. In particular, cutting the geometry into quadrants through the origin in the

xz

- and

yz

-planes will make the polyhedron meshing scheme particularly well-suited for this geometry. This webcutting can be achieved with

webcut volume all with plane xplane offset 0
webcut volume all with plane yplane offset 0

merge surface all

As noted previously, the polyhedron meshing scheme will be applied to all of the geometric entities with

volume all scheme polyhedron

Next, the element sizes can again be approximated using the same approach as detailed in the 2D example. Here we will consider the following material properties:

Entity	Physics	`VP`	`VS`	`RHO`
Hemisphere	Acoustic	1500	N/A	1000
Cuboid	Elastic	4500	3000	2000

We'll again discretize the domain with approximately 2.0 elements per wavelength, but for this simulation we'll use a source frequency of 20 kHz. This results in the approximate element sizes

\Delta x

\begin{aligned} \Delta x_{\text{hem}} &= \frac{1500 \; [\text{m/s}]}{20000 \; [\text{Hz}] \times 2.0} \\ &= 0.0375 \; [\text{m}] \end{aligned}

and

\begin{aligned} \Delta x_{\text{cub}} &= \frac{3000 \; [\text{m/s}]}{20000 \; [\text{Hz}] \times 2.0} \\ &= 0.075 \; [\text{m}]. \end{aligned}

These element sizes can then be applied with

volume 3 5 7 9 size 0.0375
volume 1 4 6 8  size 0.075
volume 3 5 7 9  sizing function type skeleton min_size auto max_size 0.0375 max_gradient 1 min_num_layers_3d 1 min_num_layers_2d 1 min_num_layers_1d 1
volume 1 4 6 8  sizing function type skeleton min_size auto max_size 0.075 max_gradient 1 min_num_layers_3d 1 min_num_layers_2d 1 min_num_layers_1d 1

The polyhedron scheme requires all of the regions of the domain to be meshed at once. Thus, we can simply run

mesh volume all

The uniformity of the elements can be improved using mesh smoothing through the use of

surface all smooth scheme edge length
surface smooth all

volume all smooth scheme equipotential
smooth volume all

Here's what the mesh looks like, along with a cut-away to reveal its internal structure:

Finally, each of the individual parts in the hemisphere and the cuboid can be added to a single block; this will become useful later when attempting to assign material properties to discrete blocks within the domain. This can be achieved using

block 1 add volume 3 5 7 9
block 2 add volume 1 4 6 8

As before, the mesh can be exported to Exodus using

set exodus netcdf4 on
export mesh "/path/to/file/sample_3D.e" dimension 3

Similar to the 2D example, here's a complete list of the journal commands necessary for creating the mesh:

# Create the cuboid
brick x 2 y 2 z 1

# Create the sphere and shift it upwards
create sphere radius 0.5
move Volume 2 z 0.5 include_merged

# Subtract the sphere from the cuboid
subtract volume 2 from volume 1 keep_tool

# Remove the top part of the sphere
webcut volume 2  with plane from surface 9
delete volume 2

# Combine the bowl parts of the surfaces
merge surface all

# Make webcuts so that we can use the polyhedron scheme on everything
webcut volume all with plane xplane offset 0
webcut volume all with plane yplane offset 0

# Merge the webcut surfaces
merge surface all

# Use the polyhedron scheme to mesh everything
volume all scheme polyhedron

# Set the element sizes
volume 3 5 7 9 size 0.0375
volume 1 4 6 8  size 0.075
volume 3 5 7 9  sizing function type skeleton min_size auto max_size 0.0375 max_gradient 1 min_num_layers_3d 1 min_num_layers_2d 1 min_num_layers_1d 1
volume 1 4 6 8  sizing function type skeleton min_size auto max_size 0.075 max_gradient 1 min_num_layers_3d 1 min_num_layers_2d 1 min_num_layers_1d 1

# Mesh everything
mesh volume all

# Smooth the mesh a bit
surface all smooth scheme edge length
smooth surface all
volume all smooth scheme equipotential
smooth volume all

# Assign everything as two blocks
block 1 add volume 3 5 7 9
block 2 add volume 1 4 6 8

# Export to Exodus
set exodus netcdf4 on
export mesh "/path/to/file/sample_3D.e" dimension 3

Importing meshes in 3D is extremely similar to how it is performed in 2D. In fact, we'll re-use many of the functions we defined in the previous example for the 3D mesh.

Copy

import os

import matplotlib.pyplot as plt
import numpy as np

from salvus.mesh.algorithms.unstructured_mesh.metrics import compute_time_step
import salvus.namespace as sn

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

We can again use the from_exodus constructor to directly import the Exodus mesh.

Copy

m_3d = sn.UnstructuredMesh.from_exodus(
    "data/sample_3D.e", attach_element_block_indices=True
)

# Add a side set to the mesh so that we can visualize the mesh
m_3d.find_surface()

m_3d

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

We can re-use the same material properties dictionary as we used in the previous tutorial:

Copy

# Define our properties for the mesh
block_props = {
    1: {"VP": 1500.0, "RHO": 1000.0},
    2: {
        "VP": 4500.0,
        "VS": 3000.0,
        "RHO": 2000.0,
    },
}

Let's apply the materials and perform the mesh QC.

Copy

def assign_properties_to_mesh(
    mesh: sn.UnstructuredMesh, properties: dict
) -> None:
    # Initialize the fluid flag
    mesh.attach_field("fluid", np.zeros(mesh.nelem, dtype=bool))

    # Initialize the elemental fields for the mesh
    for field in ["VP", "VS", "RHO"]:
        mesh.attach_field(field, np.zeros_like(mesh.connectivity, dtype=float))

    for block_id in properties:
        # Get the elements that belong to the given block ID
        mask = mesh.elemental_fields["element_block_index"] == block_id

        # Assign the material properties for that block
        for field, value in properties[block_id].items():
            mesh.elemental_fields[field][mask, :] = value

        # Set the fluid flag for the given block; this will tell the solver
        # whether to use the acoustic or elastic solver for that part of the
        # mesh
        if "VS" not in properties[block_id]:
            mesh.elemental_fields["fluid"][mask] = True
        else:
            mesh.elemental_fields["fluid"][mask] = False


assign_properties_to_mesh(m_3d, block_props)
m_3d.qc_test()

[SUCCESS] No issues found!

{}

Everything there looks good, so let's again proceed to looking at the resolved frequencies and the global time step that the mesh is able to obtain.

As within the 2D example, let's compute the minimum resolved frequencies and the global time step that the mesh is able to achieve.

Copy

min_velocity = np.zeros(m_3d.nelem, dtype=float)

for block_id, props in block_props.items():
    mask = m_3d.elemental_fields["element_block_index"] == block_id

    if "VS" in props:
        min_velocity[mask] = np.min(
            m_3d.elemental_fields["VS"][mask, :], axis=1
        )
    else:
        min_velocity[mask] = np.min(
            m_3d.elemental_fields["VP"][mask, :], axis=1
        )

min_frequency_3d, resolved_frequencies_3d = m_3d.estimate_resolved_frequency(
    min_velocity=min_velocity,
    elements_per_wavelength=2.0,
)

print("Minimum resolved frequency is %.2f kHz" % (min_frequency_3d / 1e3))

Minimum resolved frequency is 14.76 kHz

The minimum resolved frequency is somewhat below our desired bound of 20kHz. This implies that there are some elements in the mesh that are larger than we might like, which might cause us to accumulate more numerical dispersion as a result. Let's take a look at the resolved frequencies for the mesh:

Copy

def colored_histogram(
    data: np.ndarray,
    n_bins: int = 50,
    colormap: str = "coolwarm",
) -> None:
    counts, bins = np.histogram(data, bins=n_bins)

    fig = plt.figure(figsize=[10, 3])
    cmap = plt.get_cmap(colormap)
    colors = cmap(np.linspace(0, 1, len(bins) - 1))

    for i in range(len(bins) - 1):
        plt.bar(
            bins[i],
            counts[i],
            width=bins[i + 1] - bins[i],
            color=colors[i],
            edgecolor="black",
        )

    return fig, counts, bins


def plot_resolved_frequencies(
    resolved_frequencies: np.ndarray,
    n_bins: int = 50,
    colormap: str = "coolwarm",
) -> None:
    _, counts, bins = colored_histogram(
        resolved_frequencies / 1e3, n_bins, colormap
    )

    plt.annotate(
        "Elements\nToo Large",
        xy=(bins[0], max(counts) * 0.45),
        xytext=(bins[6], max(counts) * 0.45),
        arrowprops=dict(facecolor="black", arrowstyle="-|>", lw=1.5),
        ha="center",
        va="center",
    )
    plt.annotate(
        "Elements\nToo Small",
        xy=(bins[-1], max(counts) * 0.45),
        xytext=(bins[-7], max(counts) * 0.45),
        arrowprops=dict(facecolor="black", arrowstyle="-|>", lw=1.5),
        ha="center",
        va="center",
    )

    plt.title("Resolved Frequencies")
    plt.xlabel("Frequency [kHz]")
    plt.ylabel("Number of Elements")
    plt.show()


def plot_dt(
    dts: np.ndarray, n_bins: int = 50, colormap: str = "coolwarm_r"
) -> None:
    _, counts, bins = colored_histogram(dts, n_bins, colormap)

    plt.annotate(
        "Restricting\nTime Step",
        xy=(bins[0], max(counts) * 0.05),
        xytext=(bins[0], max(counts) * 0.3),
        arrowprops=dict(facecolor="black", arrowstyle="-|>", lw=1.5),
        ha="center",
        va="center",
    )

    plt.title("Time Steps")
    plt.xlabel(r"$\Delta t$ [$\mu$s]")
    plt.ylabel("Number of Elements")
    plt.show()


plot_resolved_frequencies(resolved_frequencies_3d)

It appears that most of the elements are, indeed, slightly below the desired bound of 20kHz. This means that we would ideally choose to mesh the domain with a slightly finer discretization in order to avoid the presence of these large elements. It is generally helpful to plot the resolved frequencies in order to identify the problematic elements in the mesh:

Copy

m_3d.attach_field("resolved_frequency", resolved_frequencies_3d)

# Change plotted field to `resolved_frequency`
m_3d

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

Note that here we are only plotting the elements on the surface of the hexahedral mesh; one could use the write_h5() method in order to visualize the internal structure of the mesh in an alternative visualization tool such as ParaView.

The tail on the right-hand-side of the resolved frequencies histogram implies that the mesh also contains some locally constricted elements, given that the upper bound of the resolved frequencies is higher than the expected maximum frequency of the source. This suggests that the time step will likely become limited by a handful of these small elements.

To verify that this is indeed the case, let's compute the time steps

\Delta t

of each element.

Copy

dt_3d, dts_3d = compute_time_step(
    mesh=m_3d,
    max_velocity=m_3d.elemental_fields["VP"].max(axis=1),
)

print(f"Global time step: {dt_3d * 1e6:.2f} us")

plot_dt(dts_3d * 1e6)

Global time step: 0.40 us

We can see in this distribution of the time steps across all of the elements in the mesh that there are a handful of elements that are restricting the time step. Let's plot the

\Delta t

on the mesh so that we can see the parts of the discretization that are problematic for the time step.

Copy

m_3d.attach_field("dt", dts_3d)

# Change plotted field to `dt`
m_3d

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

At this stage, one could return to the meshing stage in order to improve the quality of the discretization. However, for the sake of this example, we will proceed with the existing mesh that we have, but with the understanding that further improvements to the mesh could be conducted.

Let's again set up a very minimalistic forward simulation to demonstrate that this mesh actually works in a wave simulation:

Copy

p_3d = sn.Project.from_mesh(
    path="project_3D",
    mesh=m_3d,
    load_if_exists=True,
)

p_3d.add_to_project(
    sn.Event(
        event_name="event_3D",
        sources=sn.simple_config.source.cartesian.ScalarPoint3D(
            x=0.0, y=0.0, z=0.5, f=1.0
        ),
    )
)

wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=500.0e-6)
stf = sn.simple_config.stf.Ricker(center_frequency=20e3 / 2)
stf.plot()

Copy

p_3d.add_to_project(
    sn.UnstructuredMeshSimulationConfiguration(
        name="forward_simulation_3D",
        unstructured_mesh=m_3d,
        event_configuration=sn.EventConfiguration(
            waveform_simulation_configuration=wsc,
            wavelet=stf,
        ),
    )
)

Before running the simulation, we'll double check the simulation setup:

Copy

p_3d.viz.nb.simulation_setup(
    simulation_configuration="forward_simulation_3D",
    events="event_3D",
)

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

That looks good to go! A few small things to note here about this simulation:

We'll only output a very small number of time steps (i.e., with sampling_interval_in_time_steps) in order to keep the file size of the wavefield small. If you would want to create a wavefield animation for this data, decreasing the sampling_interval_in_time_steps would output more time steps. However, be warned that the files can quickly become very large!
Increasing the number of ranks in accordance with the number of CPU cores that your machine has will make the simulation run much more quickly.

Copy

p_3d.simulations.launch(
    simulation_configuration="forward_simulation_3D",
    events="event_3D",
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=RANKS_PER_JOB,
    extra_output_configuration={
        "volume_data": {
            "sampling_interval_in_time_steps": 1000,
            "fields": ["displacement", "gradient-of-phi"],
        },
    },
)

[2025-10-03 21:15:13,919] INFO: Submitting job ...
Uploading 1 files...

🚀  Submitted job_2510032115022274_07912f900d@local

Copy

p_3d.simulations.query(
    simulation_configuration="forward_simulation_3D",
    events="event_3D",
    block=True,
)

True

Finally, we can again visualize the wavefield data in a visualization software such as ParaView. As with the 2D example, the relevant state file is provided in order to plot this wavefield in ParaView. Here is a snapshot of the final time step of the wavefield: