Mondaic - Full Waveform Solutions

This tutorial uses the W-Scan module, which is a separate Python package from Salvus. It contains a suite of next-generation tools that are tailored to use-cases in non-destructive testing (NDT).

W-Scan is currently available as a closed beta; APIs may change without notice. To become an early user of W-Scan, please contact [email protected].

Pitting corrosion is a localized form of corrosion that produces small, surface-breaking depressions in metallic specimens. In 2D, a pit is well approximated by a triangular cross-section: a wide mouth at the surface that tapers linearly to a single point at depth. This tutorial demonstrates how to set up an ultrasonic NDT simulation in which a contact probe illuminates a conical pit on the bottom surface of a steel plate.

As usual, a series of standardized libraries are imported.

import os

import build123d
import matplotlib.pyplot as plt
import numpy as np
import salvus.namespace as sn
import wscan

For convenience, a number of constants are defined below. Constant parameters are denoted with variable names in all capital letters.

First, the simulation site and the location to which the project should be saved are defined.

# Salvus site on which simulations are to be run
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")

# Path where the project should be saved
PATH_TO_PROJECT = "project_cone_pitting_tutorial"

Next, parameters that control the mesh resolution and the source time function are set.

# Maximum frequency of the source in Hz; also controls how finely the mesh
# is to be discretized
MAX_FREQUENCY = 10.0e6
CENTER_FREQUENCY = 3 / 4 * MAX_FREQUENCY

# Number of elements per wavelength; one of the key parameters that controls
# how fine the mesh is.  Refer to
# https://docs.mondaic.com/knowledge_base/spectral_element_modelling/hp_refinement
# for more details.
ELEMENTS_PER_WAVELENGTH = 2.0

# Simulation end time in seconds
SIMULATION_END_TIME = 5.0e-6

The material properties of the steel plate are taken from W-Scan's material library.

SPECIMEN_MATERIAL = wscan.library.materials.materials["steel"].material

Finally, the geometric properties of the domain are defined. The plate's top surface sits at y = 0 and the specimen extends downwards into the negative y half-plane. The pit is introduced on the bottom surface, centred on x = 0, and opens downward into the void below the plate.

# Horizontal extents of the plate in meters.
X_MIN = -10.0e-3
X_MAX = 10.0e-3

# Plate thickness in meters
THICKNESS = 8.0e-3

# Diameter of the pit mouth in meters
PIT_DIAMETER = 1.0e-3

# Maximum depth of the pit in meters
PIT_DEPTH = 1.0e-3

The conical pit is constructed using wscan.geometry.dim2.defects.ConePit. The origin is placed at the centre of the bottom edge of the plate, and the normal is the outward-pointing surface normal. Similar to the various cracks implementations in the wscan.geometry.dim2.defects submodule, ConePit entities have their surface roughness randomly vary.

pit = wscan.geometry.dim2.defects.ConePit(
    diameter=PIT_DIAMETER,
    depth=PIT_DEPTH,
    origin=np.array([0.0, -THICKNESS]),
    normal=np.array([0.0, -1.0]),
    n_samples=9,
    lateral_variability=0.05,
    random_seed=42,
)
pit.plot()

The steel plate is constructed as a rectangular build123d face, and the pit's polygonal cross-section is subtracted from it to carve out the void at the bottom surface.

# Create the plate geometry; get the build123d face as a `Face` object for the
# subsequent Boolean subtraction
plate = build123d.Polygon(
    [
        (X_MIN, 0.0),
        (X_MIN, -THICKNESS),
        (X_MAX, -THICKNESS),
        (X_MAX, 0.0),
    ]
).face()

# Add the pitting to the plate.  The pit's `geometry` property returns a
# `build123d.Face` object that we can use to perform a Boolean subtraction with
# the plate.
specimen = plate - pit.geometry

# Plot the entity
specimen.label = "specimen"
wscan.geometry.dim2.utils.plot_build123d_entity(specimen)
plt.show()

The assembly is meshed via Coreform Cubit using the pave algorithm. As with the other tutorials, the resolution is governed by reference_frequency and elements_per_wavelength, which together determine the element size required to resolve waves at MAX_FREQUENCY.

In order for this tutorial to be run by users that do not have Coreform Cubit installed, an example mesh is provided in the supplementary data directory. However, if one were to construct the mesh themselves, one could do so using the following command:

mesh = wscan.mesh.mesh_from_geometry(
    geometry=specimen,
    model=SPECIMEN_MATERIAL,
    mesh_resolution=sn.MeshResolution(
        reference_frequency=MAX_FREQUENCY,
        elements_per_wavelength=ELEMENTS_PER_WAVELENGTH,
    ),
    algorithm=wscan.mesh.algorithms.cubit.pave.Mesher(
        wscan.mesh.algorithms.cubit.pave.Options(
            site=SALVUS_FLOW_SITE_NAME,
            composite_curves=False,
        )
    ),
)

For this tutorial, however, the mesh is loaded from disk.

mesh = sn.UnstructuredMesh.from_h5("data/plate_with_pitting.h5")

The generated mesh can be visually inspected to confirm that the pit geometry is correctly resolved.

mesh.find_side_sets()
mesh

<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7368e8aa4850>

A new Salvus project is initialized from the mesh.

p = sn.Project.from_mesh(PATH_TO_PROJECT, mesh)

A ContactTransducer2D probe is placed on the top surface of the specimen so that it fires vertically toward the pit on the bottom. This type of probe injects a plane wave into the specimen, which propagates to the cone pit at the base of the specimen.

probe = wscan.instruments.transducer.ContactTransducer2D(
    center=np.array([0.0, 0.0]),
    normal=np.array([0.0, 1.0]),
    width=10.0e-3,
    maximum_frequency=MAX_FREQUENCY,
    reference_velocity=SPECIMEN_MATERIAL.VS.p,
    sampling_factor=5.0,
)

The probe's get_sources method instantiates the sources that correspond to each transducer element. The excitation_mode="longitudinal" specifies that each element radiates longitudinal (compressional) waves into the specimen.

p.add_to_project(
    sn.Event(
        sources=probe.get_sources(
            excitation_mode="longitudinal",
            source_type=sn.simple_config.source.cartesian.SideSetVectorPoint2D,
            side_set_name="y1",
        ),
    )
)

A ToneBurst source time function is used to produce a narrowband excitation centred at approximately three-quarters of the maximum frequency.

stf = sn.simple_config.stf.ToneBurst(
    center_frequency=CENTER_FREQUENCY,
    num_cycles=4,
)
stf.plot()

An UnstructuredMeshSimulationConfiguration is added for the pit geometry, with absorbing boundary conditions on the left (x0) and right (x1) plate edges to suppress spurious reflections from those boundaries. For the sake of simplicity, sponge layers are disabled by setting width_in_meters=0.0; only first-order absorbing conditions are applied directly on the boundary faces. Any exterior boundaries that are not explicitly assigned a boundary condition are implicitly treated as free surfaces.

p.add_to_project(
    sn.UnstructuredMeshSimulationConfiguration(
        name="pit",
        unstructured_mesh=mesh,
        event_configuration=sn.EventConfiguration(
            wavelet=stf,
            waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
                end_time_in_seconds=SIMULATION_END_TIME,
                boundary_conditions=sn.simple_config.boundary.Absorbing(
                    taper_amplitude=MAX_FREQUENCY,
                    side_sets=["x0", "x1"],
                    width_in_meters=0.0,
                ),
            ),
        ),
    )
)

As a quick sanity check, the geometry is plotted together with the transducer source positions to confirm that the elements lie over the cone pit.

ax = wscan.geometry.dim2.utils.plot_build123d_entity(
    specimen, legend=True, figsize=(12, 12)
)

for source in p.events.get(p.events.list()[0]).sources:
    plt.plot(source.point[0], source.point[1], "k.")

plt.show()

In order to keep the volumetric wavefield file size manageable, the volumetric output is downsampled in time via sampling_interval_in_time_steps, which is chosen here so that roughly 300 snapshots are produced across the full simulation window. The time step is computed for the mesh under the Courant-Friedrichs-Lewy (CFL) condition to be able to approximate the appropriate sampling interval.

time_step, _ = mesh.compute_time_step(
    max_velocity=np.max(mesh.elemental_fields["VP"], axis=1)
)

print(f"Time step: {time_step} s")

Time step: 7.983279681674821e-10 s

The simulation is launched and the query(block=True) call blocks the cell until the simulation has finished.

p.simulations.launch(
    simulation_configuration="pit",
    events=p.events.list(),
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=sn.api.get_site(SALVUS_FLOW_SITE_NAME).config[
        "default_ranks"
    ],
    extra_output_configuration={
        "volume_data": {
            "sampling_interval_in_time_steps": int(
                (SIMULATION_END_TIME - stf.get_auto_start_time())
                / time_step
                / 300
            ),
            "fields": ["displacement"],
        },
    },
)
p.simulations.query(block=True)

[2026-06-09 18:26:24,486] INFO: Submitting job ...
Uploading 2 files...

🚀  Submitted job_2606091826624243_55857f602f@local

True

Once the simulation has finished, the volumetric wavefield is interpolated onto a regular Cartesian grid covering the region around the pit. The interpolation is carried out in the spectral-element basis and, therefore, inherits the full accuracy of the spectral-element solution.

Note: Points in the grid that fall inside the void of the pit, which contains no finite elements, will remain unclaimed during the interpolation. The resulting warning from Salvus can safely be ignored in this case.

event_data = p.waveforms.get("pit", p.events.list()[0])[0]

wavefield = sn.wavefield_output_to_xarray(
    wo=event_data.get_wavefield_output(
        output_type="volume",
        field="displacement",
    ),
    points=[
        np.linspace(-probe.width / 2, probe.width / 2, 301),
        np.linspace(-THICKNESS, 0.0, 301),
    ],
)

[2026-06-09 18:26:33,576] WARNING - salvus.mesh.algorithms.unstructured_mesh.utils: 577 points were not claimed by enclosing elements. Depending on your use case, this may not be an issue.

Several wavefield snapshots are plotted in order to be able to better understand the interactions of the plane wave as it propagates from one side of the plate to the other. In particular, the waves that reflect off of the cone pit are of interest.

# Approximate time at which the wave arrives at the pit
t_plot = np.linspace(0.5e-6, 3.0e-6, 6)

# Fixed color range for the displacement field
vrange = 1.0e-11

# Create the plot
_, axs = plt.subplots(
    nrows=2, ncols=3, figsize=(10, 6), sharex=True, sharey=True
)

for i, ax in enumerate(axs.flatten()):
    # Plot the outline of the geometry
    ax = wscan.geometry.dim2.utils.plot_build123d_entity(
        specimen,
        face_alpha=0.0,
        ax=ax,
    )

    # Get the time index that is closest to `t_plot`
    t_ind = np.argmin(np.abs(wavefield.t.to_numpy() - t_plot[i]))

    # Plot the wavefield; `c=1` is the y-component of displacement, `c=0` is
    # the x-component
    wavefield.sel(c=1, t=wavefield.t[t_ind]).T.plot(
        shading="gouraud",
        infer_intervals=False,
        ax=ax,
        cmap="coolwarm",
        zorder=-1,
        add_colorbar=False,
        vmin=-vrange,
        vmax=vrange,
    )

    # Formatting
    ax.set_aspect("equal")
    ax.set_xlabel(None)
    ax.set_ylabel(None)

# Add x- and y-labels
for ax in axs[:, 0]:
    ax.set_ylabel("y [m]")
for ax in axs[-1, :]:
    ax.set_xlabel("x [m]")

plt.tight_layout()
plt.show()

In order to visualize the entire wavefield, one can alternatively view the raw simulation output in a third-party visualization tool, such as ParaView. The path to the simulation output directory is listed below:

print(
    p.simulations.get_simulation_output_directory(
        simulation_configuration="pit", event=p.events.list()[0]
    )
)

project_cone_pitting_tutorial/EVENTS/432251f5-375a-431a-bd3c-e2a016427c5d/WAVEFORM_DATA/INTERNAL/73/bb/28ac7b96e0ad

For more details on visualizing wavefields with ParaView, see here.

Notched specimens

Part 0 - Prologue

© Mondaic AG (2025)

Site Map Contact Us

Academic License Agreement

PAGE CONTENTS