Mondaic Documentation
www.mondaic.com
Mondaic
Version:
No results found.
W-Scan (Beta)

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].

Maximum Amplitude Distributions

When an ultrasonic probe interrogates a defect, it is often informative to know not just what the wavefield looks like at a single instant, but also how strongly each point in the specimen is insonified over the entire duration of the inspection. Mapping the peak amplitude reached at every point yields a "maximum amplitude distribution", which represents a single image that summarizes the region of the specimen that the probe most strongly excites, including any amplification or shadowing introduced by a defect.
This tutorial follows the same general workflow as the notched specimen tutorial, reusing the same thin, jagged surface-breaking Crack notch as the defect. A wedge probe is placed on the top surface of the plate, a tone-burst excitation is injected into the medium, and the resulting volumetric wavefield is post-processed into a maximum amplitude distribution.

Importing libraries and setting constants

As usual, a series of standardized libraries are imported.
Copy
import os

import build123d
import matplotlib.pyplot as plt
import numpy as np
from salvus.mesh.unstructured_mesh_utils import extract_model_to_regular_grid
import salvus.namespace as sn
import wscan
import xarray as xr
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_maximum_amplitude_tutorial"
Next, parameters that control the mesh resolution and the simulation are set.
# Maximum frequency of the source in Hz; also controls how finely the mesh is
# to be discretized
MAX_FREQUENCY = 2.5e6

# Due to the number of wavelengths that will be propagated, the
# spectral-element order is chosen as 7
SPECTRAL_ELEMENT_ORDER = 7

# 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 = 1.0

# Simulation end time in seconds
SIMULATION_END_TIME = 25.0e-6
The material properties of the various constituents can be defined using sn.material.Material objects. These can either be defined manually or with W-Scan's material library.
The specimen in this example is assumed to consist of a test piece that is composed of acrylic glass. This medium results in a relatively low impedance contrast between the wedge and the plate, which allows for the illumination of the crack to be more effectively visualized in the maximum amplitude distribution.
# The material properties of the wedge and the plate are taken from the
# material library
WEDGE_MATERIAL = wscan.library.materials.materials["Epoxy resin"].material
SPECIMEN_MATERIAL = wscan.library.materials.materials["Acrylic glass"].material
Finally, the geometric properties of the domain can be defined. The plate shares the same extents as in the notched specimen tutorial. By convention, the top of the specimen sits at y = 0 and the specimen extends downwards into the negative y half-plane; the wedge probe sits near the left edge and the crack notch breaks the bottom surface near the centre of the plate and extends upwards (the SimpleSpecimen class internally rotates the supplied notch by 180 degrees), where the obliquely incident wavefield reaches it.
# Horizontal extents of the plate in meters. The wedge probe sits near the left
# edge and the crack notch is positioned near the centre.
X_MIN = -30.0e-3
X_MAX = 10.0e-3

# Plate thickness in meters
THICKNESS = 10.0e-3

# Crack notch dimensions in meters.
NOTCH_WIDTH = 0.5e-3
NOTCH_DEPTH = 5.0e-3
The domain geometry is built programmatically using build123d. The notched specimen and the probe are defined separately before being combined into a single compound.

Crack notch

Rather than building the plate by hand, the specimen is assembled with the wscan.geometry.dim2.notches.SimpleSpecimen helper, which combines a rectangular plate (spanning X_MIN to X_MAX horizontally and extending from the surface at y = 0 down to y = -THICKNESS) with an embedded notch. The resulting face is labeled specimen, which is later used to assign material properties during the meshing process.
The defect is a thin, jagged surface-breaking Crack, constructed exactly as in the notched specimen tutorial. The Crack is fully parameterized: direction orients the crack as it propagates from the surface towards its tip, depth_in_meters and width_in_meters set its vertical extent and mouth opening, and lateral_variability controls the roughness of the crack walls.
The fixed random_seed ensures the result is reproducible. However, one could alternatively vary this within a parametric study in order to simulate a collection of different crack geometries.
specimen = wscan.geometry.dim2.notches.SimpleSpecimen(
    notch=wscan.geometry.dim2.notches.Crack(
        direction=np.array([-0.25, -1.0]),
        depth_in_meters=NOTCH_DEPTH,
        width_in_meters=NOTCH_WIDTH * 0.1,
        lateral_variability=0.05,
        random_seed=9001,
    ),
    x_left_in_meters=X_MIN,
    x_right_in_meters=X_MAX,
    thickness_in_meters=THICKNESS,
)
specimen.plot()
A WedgeProbe2D is placed on the top surface of the specimen so that it fires obliquely toward the crack. The wedge geometry, transducer layout, and damping layer follow the same dimensions as the other NDT tutorials:
  • origin places the probe near the left edge of the specimen, away from the crack region.
  • wedge_angle_in_degrees=45 produces an obliquely incident wavefront in the specimen.
  • A large number_of_elements is used to synthesise a smooth plane wave.
  • damping_depth introduces an absorbing layer on the back face of the wedge that suppresses spurious reverberations from the probe body.
probe = wscan.instruments.probe.WedgeProbe2D(
    material=wscan.library.materials.materials["Epoxy resin"].get_material(),
    origin=np.array([-0.02, 0.0]),
    transducer_position=np.array([2.4e-3, 9.4e-3]),
    front_face_distance=16.9e-3,
    back_face_distance=2.7e-3,
    height=14.2e-3,
    wedge_angle_in_degrees=45.0,
    number_of_elements=256,
    element_width=3.75e-05,
    inter_element_spacing=0.0,
    damping_depth=1.0e-3,
    receiver_fields=["displacement"],
    maximum_frequency=MAX_FREQUENCY,
)
probe.plot()
<Axes: xlabel='x [m]', ylabel='y [m]'>
Currently all of the geometric entities are disjoint from one another. In order to mesh the domain, all components need to be combined into a single build123d.Compound object:
  • specimen.geometry: the notched specimen face, with the crack void already subtracted from the plate and the specimen label preserved.
  • *probe.geometry.children: unpacks the probe's wedge and damping sub-bodies.
geometry = build123d.Compound(
    children=(*probe.geometry.children, specimen.geometry),
    label="assembly",
)

ax = wscan.geometry.dim2.utils.plot_build123d_entity(geometry, legend=True)
plt.show()
The assembly is meshed via Coreform Cubit using the notched_specimen algorithm, a meshing strategy tailored to this class of geometries. In particular, the assignment of suitable element sizes, side sets, and materials is handled automatically by this meshing strategy. The compound contains exactly three labeled faces — wedge, damping, and specimen — each of which is assigned to the corresponding material properties through the model.
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, examples of each of the meshes are provided in the supplementary data directory. However, if one were to construct these meshes themselves, one could do so using the following command:
mesh = wscan.mesh.mesh_from_geometry(
    geometry=geometry,
    model={
        "wedge": WEDGE_MATERIAL,
        "damping": WEDGE_MATERIAL,
        "specimen": SPECIMEN_MATERIAL,
    },
    mesh_resolution=sn.MeshResolution(
        reference_frequency=MAX_FREQUENCY,
        elements_per_wavelength=ELEMENTS_PER_WAVELENGTH,
        model_order=SPECTRAL_ELEMENT_ORDER,
    ),
    algorithm=wscan.mesh.algorithms.cubit.notched_specimen.Mesher(
        wscan.mesh.algorithms.cubit.notched_specimen.Options(site=SALVUS_FLOW_SITE_NAME)
    ),
)
However, for this tutorial, the meshes are instead loaded from disk.
mesh = sn.UnstructuredMesh.from_h5("./data/crack.h5")
The generated mesh can be visually inspected to confirm that the crack is properly resolved and that the wedge–specimen interface looks reasonable.
mesh
<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7d47814e9350>
A new Salvus project is initialized from the mesh.
p = sn.Project.from_mesh(PATH_TO_PROJECT, mesh)
The probe's get_event method instantiates the sources and receivers corresponding to each transducer element. The excitation_mode="longitudinal" specifies that each element radiates longitudinal (compressional) waves into the wedge.
p.add_to_project(
    probe.get_event(
        excitation_mode="longitudinal",
        source_type=sn.simple_config.source.cartesian.SideSetVectorPoint2D,
        receiver_type=sn.simple_config.receiver.cartesian.SideSetPoint2D,
    )
)
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=3 / 4 * MAX_FREQUENCY,
    num_cycles=4,
)
stf.plot()
The simulation configuration ties together the mesh, the event, the source time function, and the solver settings.
Absorbing boundary conditions are applied on the left (x0) and right (x1) plate edges, in addition to the probe damping layer (damping), in order 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="forward_simulation",
        unstructured_mesh=mesh,
        event_configuration=sn.EventConfiguration(
            wavelet=stf,
            waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
                end_time_in_seconds=SIMULATION_END_TIME,
                spectral_element_order=SPECTRAL_ELEMENT_ORDER,
                boundary_conditions=sn.simple_config.boundary.Absorbing(
                    taper_amplitude=MAX_FREQUENCY,
                    side_sets=["x0", "x1", "damping"],
                    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 on the wedge's inclined face.
ax = wscan.geometry.dim2.utils.plot_build123d_entity(
    geometry, 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()
A notable difference between this simulation from the ones performed in the aforementioned notches tutorial is that the maximum amplitude of the wavefield is to be computed for each point in space over the course of the simulation. This would normally require one to output the entire wavefield with a fine temporal resolution, which is intractable for larger-scale problems.
As an alternative to outputting the entire time-domain wavefield and computing the maxima in time as a post-processing step, one can instead compute the frequency-domain wavefield using the FFT on-the-fly. This approach allows for the Fourier-transformed wavefield to be output for a pre-specified selection of frequencies, which dramatically reduces the storage requirements for subsequently computing the maximum amplitude distributions.
The on-the-fly FFT is applied by using the frequency_domain key in the extra_output_configuration. Recall that the source time function is a tone burst with the following temporal characteristics:
stf.plot()
As can be seen above, most of the frequency content of the source is concentrated between approximately
12ωmax<ω<ωmax,\frac{1}{2} \omega_{\text{max}} < \omega < \omega_{\text{max}},
where ωmax\omega_{\text{max}} is the approximate maximum resolved frequency of the source. In the interest of keeping the computational cost and storage requirements limited for this example, 25 equispaced frequencies are selected within this range.
p.simulations.launch(
    simulation_configuration="forward_simulation",
    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={
        "frequency_domain": {
            "frequencies": np.linspace(
                MAX_FREQUENCY / 2, MAX_FREQUENCY, 25
            ).tolist(),
            "fields": ["displacement"],
            "start_time_in_seconds": 0.0,
            "end_time_in_seconds": SIMULATION_END_TIME,
        },
    },
    delete_conflicting_previous_results=True,
)
p.simulations.query(block=True, ping_interval_in_seconds=1.0)
[2026-06-16 03:58:20,213] INFO: Submitting job ...
Uploading 2 files...

🚀  Submitted job_2606160358404832_24731e4a04@local

True
The frequency-domain output written by the solver is loaded back as an UnstructuredMesh. It carries one complex-valued FFT-* elemental field per requested frequency and displacement component, alongside an auxiliary fluid field that is not needed here and is, therefore, discarded.
fft_wavefield = sn.UnstructuredMesh.from_h5(
    p.simulations.get_simulation_output_directory(
        simulation_configuration="forward_simulation", event=p.events.list()[0]
    )
    / "frequency_domain_output.h5"
)
del fft_wavefield.elemental_fields["fluid"]
A single scalar measure of how strongly each point is insonified is then assembled from the individual Fourier components. This can be achieved using Parseval's theorem, which can be applied using
umax2(x)=1N∑ω(∥u^x(x,  ω)∥2+∥u^y(x,  ω)∥2),u_{\text{max}}^2 (\mathbf{x}) = \frac{1}{N} \sum_\omega \left( \left\| \hat{u}_{x} (\mathbf{x}, \;\omega) \right\|^2 + \left\| \hat{u}_{y} (\mathbf{x}, \;\omega) \right\|^2 \right),
where umax(x)u_{\text{max}} (\mathbf{x}) is the maximum amplitude at each spatial coordinate x\mathbf{x}, NN is the number of discrete frequencies, and u^(x,y)(x,  ω)\hat{u}_{(x, y)} (\mathbf{x}, \;\omega) are the frequency-domain representations of the x- and y-components of the displacement wavefield.
The FFT on-the-fly exports four fields per frequency:
  • FFT-Re_#_X: Real part of the x-component of the displacement wavefield
  • FFT-Im_#_X: Imaginary part of the x-component of the displacement wavefield
  • FFT-Re_#_Y: Real part of the y-component of the displacement wavefield
  • FFT-Im_#_Y: Imaginary part of the y-component of the displacement wavefield
The maximum amplitude can thus be computed directly:
# Attach the field that will store the maximum amplitudes
fft_wavefield.attach_field(
    "maximum_amplitude",
    np.zeros_like(fft_wavefield.connectivity, dtype=float),
)

# Get the names of the elemental fields that contain the FFT wavefield for each
# discrete frequency
fft_fields = [
    key
    for key in fft_wavefield.elemental_fields.keys()
    if key.startswith("FFT-")
]

for field in fft_fields:
    # Add the square of the FFT field to the maximum field
    fft_wavefield.elemental_fields["maximum_amplitude"] += (
        fft_wavefield.elemental_fields[field] ** 2
    )

    # Remove the field from the mesh; this is only done to make it easier to
    # find the `maximum_amplitude` field in the mesh visualization widget
    fft_wavefield.elemental_fields.pop(field)

# Apply the normalization factor of 1 / N and take the square root
fft_wavefield.elemental_fields["maximum_amplitude"] = np.sqrt(
    fft_wavefield.elemental_fields["maximum_amplitude"] / (len(fft_fields) / 4)
)

# Normalize for ease of plotting
fft_wavefield.elemental_fields[
    "maximum_amplitude"
] /= fft_wavefield.elemental_fields["maximum_amplitude"].max()
There are two options for plotting the maximum amplitude distributions:
  1. Plot the spectral-element mesh containing the maximum_amplitude field using Salvus' mesh visualization widget
  2. Interpolate the data to a regular grid and plot this using matplotlib
Plotting the spectral-element mesh can be done by calling the mesh object directly:
fft_wavefield
<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7d47697aab10>
To visualise the result with matplotlib, the normalised maximum-amplitude field is interpolated from the spectral-element mesh onto a regular Cartesian grid spanning the full extent of the model. 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 crack void (which contains no finite elements) will remain unclaimed during the interpolation. The resulting warning from Salvus can safely be ignored.
# Define the grid onto which to interpolate the mesh
grid = xr.Dataset(
    coords={
        "x": np.linspace(
            fft_wavefield.points[:, 0].min(),
            fft_wavefield.points[:, 0].max(),
            601,
        ),
        "y": np.linspace(
            fft_wavefield.points[:, 1].min(),
            fft_wavefield.points[:, 1].max(),
            601,
        ),
    }
)

# Perform the interpolation
envelope = extract_model_to_regular_grid(
    fft_wavefield, grid, ["maximum_amplitude"]
)

# Plot the data
plt.figure(figsize=(8, 8))
envelope.maximum_amplitude.T.plot(
    cmap="plasma",
    vmin=0.0,
    vmax=0.6,
    add_colorbar=False,
)
wscan.geometry.dim2.utils.plot_build123d_entity(
    geometry,
    face_alpha=0.0,
    ax=plt.gca(),
)
plt.show()
One approach may be more suitable than the other for visualizing this type of data, depending on one's use-case.
Previous
Conical pitted specimens
Next
Part 0 - Prologue
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