Mondaic - Full Waveform Solutions

In Parts 1 and 2 of this series, we explored the structure of meshes within Salvus and discussed how we can work with fields in the mesh. Part 3 of this series focuses on two concepts for these types of unstructured meshes, namely:

The importance of Jacobian matrices.
Estimating the time step of the subsequent wave simulation, based on the quality of the mesh.

A particular focus of this tutorial is to analyze the influence that these factors have when considering the overall compute cost of the subsequent wave simulations.

import functools

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import numpy.typing as npt
import xarray as xr

import salvus.namespace as sn
from salvus import fem
from salvus.mesh.algorithms.unstructured_mesh import metrics
from salvus.mesh.layered_meshing import simple_post_refinement

To begin, we’ll construct a simple 4-element 2D mesh covering the square domain

[ -1, +1 ]^2

; this is extremely similar to the mesh that was used in Part 2 of this series. This basic example will serve as our starting point for exploring meshing concepts before we move on to more complex configurations.

mesh = sn.layered_meshing.mesh_from_domain(
    domain=sn.domain.dim2.BoxDomain.from_bounds((-1, +1), (-1, +1)),
    model=sn.material.from_params(vp=2.0, rho=1.0),
    mesh_resolution=sn.MeshResolution(
        reference_frequency=1.0, elements_per_wavelength=2.0
    ),
)

# Plot the mesh
mesh

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

Let's plot each of the node positions in the mesh, in addition to their indices:

fig = plt.figure(figsize=(5, 5))
ax = fig.subplots(1, 1)
ax.scatter(*mesh.points.T, c="b")

# Label each point with its id.
for i, (x, y) in enumerate(mesh.points):
    ax.annotate(
        str(i),
        (x, y),
        xytext=(5, 5),
        textcoords="offset points",
        fontsize=8,
        ha="left",
    )

ax.set_aspect("equal")
ax.set_xlabel("x [m]")
ax.set_ylabel("x [m]")
plt.show()

Recall from Part 1 that the index of each node is what we use to construct the connectivity of the mesh. What we have so far resembles a finite-difference (FD) grid: an evenly spaced set of points covering our domain of interest. What turns this point cloud into a mesh is exactly this connectivity array, which specifies how each element is built from these points. We’ll start by inspecting the node ordering for the elements in this mesh, which already reflects this connectivity.

mesh.get_element_nodes()

array([[[-1., -1.],
        [ 0., -1.],
        [-1.,  0.],
        [ 0.,  0.]],

       [[ 0., -1.],
        [ 1., -1.],
        [ 0.,  0.],
        [ 1.,  0.]],

       [[-1.,  0.],
        [ 0.,  0.],
        [-1.,  1.],
        [ 0.,  1.]],

       [[ 0.,  0.],
        [ 1.,  0.],
        [ 0.,  1.],
        [ 1.,  1.]]])

mesh.get_element_nodes() is essentially a convenience function; under the hood, it just returns mesh.points[mesh.connectivity]. This means the connectivity array is the real source of new information here, defining how points are grouped into elements.

In Salvus, the nodes of each element are ordered using a tensor product ordering. For a quadrilateral reference element spanning the coordinates

[-1, 1]^2

in the

r

(horizontal) and

s

(vertical) directions, this means nodes are arranged by stepping through one coordinate (e.g.,

r

) while holding the other (

s

) fixed, and then moving to the next value of

s

. Here's an example of what this type of ordering looks like:

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(7, 4))

axs[0].set_title("Tensor Product Ordering")
axs[1].set_title("Counterclockwise Ordering")
axs[0].set_ylabel("$s$")

corner_nodes = np.array(
    [
        [-1.0, -1.0],
        [1.0, -1.0],
        [-1.0, 1.0],
        [1.0, 1.0],
    ]
)

# Plot the edges with the order being counterclockwise (CCW)
ccw_ordering = np.array([0, 1, 3, 2])

for ax in axs:
    ax.set_aspect("equal")
    ax.fill(
        corner_nodes[ccw_ordering, 0],
        corner_nodes[ccw_ordering, 1],
        facecolor="none",
        edgecolor="k",
        linestyle="--",
    )
    ax.plot(corner_nodes[:, 0], corner_nodes[:, 1], "bo")
    ax.set_xlabel("$r$")

# Place the annotation slightly inside of the element
for i, corner_node in enumerate(corner_nodes):
    axs[0].annotate(
        i, xy=corner_node * 0.8, fontsize=20, ha="center", va="center"
    )
    axs[1].annotate(
        ccw_ordering[i],
        xy=corner_node * 0.8,
        fontsize=20,
        ha="center",
        va="center",
    )

plt.show()

This ordering is critical. Many quantities we’ll compute later, such as basis functions and interpolation operators, depend on the nodes being stored in this consistent tensor product sequence. Other FEM implementations may use different node-ordering conventions (ordering the nodes counterclockwise is common, for example). As a result, when importing meshes from external codes, it is often necessary to permute the connectivity so that the node ordering matches the expectations of Salvus. Ensuring the correct ordering is essential, as many FEM operations rely on this consistency for accurate computation. These checks are included as part of the mesh qc module.

In the following cells, we will construct a simple mesh-plotting utility using Matplotlib. While Salvus and third-party tools such as ParaView offer advanced visualization capabilities, our focus here will be on producing a minimal plotter that directly uses the mesh’s point and connectivity data. This will allow us to clearly illustrate the structure of the mesh in a controlled, step-by-step manner.

Let's start by just plotting the points of a single element in their stored (tensor product) ordering.

# Plot the points of the first element.
plt.plot(*mesh.get_element_nodes()[0].T, "--o", color="k", markerfacecolor="b")
plt.axis("equal")
plt.show()

While each vertex is plotted in the correct location, the visualization remains incomplete. What is missing is a connectivity permutation — a reordering of the points within each element to produce a more intuitive layout for plotting. A common approach is to order the nodes so that applying the right-hand rule yields a normal vector pointing out of the page. We also need to "close" the element by duplicating its origin node.

# Permute the connectivity according to the right-hand rule.
plt.plot(
    *mesh.get_element_nodes()[0, [0, 1, 3, 2, 0]].T,
    "--o",
    color="k",
    markerfacecolor="b",
)
plt.axis("equal")
plt.show()

The visualization now resembles a single element. Next, we will create a function capable of plotting multiple elements and overlaying simple fields, such as a velocity distribution, directly onto each element. The implementation details are not critical here; each line is commented and should be self-explanatory. Our focus will remain on the mechanics of the mesh data structure rather than the plotting code itself.

def plot_mesh(
    mesh: sn.UnstructuredMesh,
    values: npt.NDArray | None = None,
    annotate: bool = False,
    title: str | None = None,
    ax: plt.Axes | None = None,
):
    """
    Plot a 2D mesh with optional color mapping and value annotations.

    Creates a visualization of an unstructured mesh where each element can be
    colored according to a scalar field. Elements are filled with colors from a
    colormap, edges are outlined in black, and vertices are marked with dots.

    Args:
        mesh: The mesh to visualize. Must be a 2D quadrilateral mesh.
        values: Scalar values to map to colors for each element. If None,
            elements are colored gray. Array length must match mesh.nelem.
        annotate: If True, display the numeric value at the center of each
            element.
        title: Title of the plot. If None, no title is added.
        axes: Matplotlib axes to plot on.  If None, a new plot is created.
    """

    # Get element nodes and reorder for visualization (right-hand rule)
    nodes = mesh.get_element_nodes()[:, [0, 1, 3, 2]]

    # Close each element by appending the first node to create filled polygons
    nodes_closed = np.concatenate((nodes, nodes[:, [0]]), axis=1)

    # Create figure and axes
    if ax is None:
        f, ax = plt.subplots()
    else:
        f = plt.gcf()
    ax.set_aspect("equal")

    if values is not None:
        # Normalize values to [0,1] range for colormap mapping
        norm_values = np.divide(
            (values - values.min()).astype(float),
            (values.max() - values.min()).astype(float),
            where=~np.isclose(values.max(), values.min()),
            out=np.full_like(values, 0.5, dtype=float),
        )

        # Apply colormap to normalized values
        colors = matplotlib.colormaps["plasma_r"](norm_values)

        # Plot each element with its corresponding color
        for i, elm in enumerate(nodes_closed):
            ax.fill(*elm.T, facecolor=colors[i], edgecolor="k", alpha=0.8)
            ax.plot(*elm.T, "ko", markersize=1)

            # Optionally add text annotation at element center
            if annotate:
                # Calculate element centroid (excluding duplicate closing node)
                cx, cy = elm[:-1].mean(axis=0)

                # Display value, centered in the moddle of the element
                ax.text(cx, cy, f"{values[i]}", ha="center", va="center")

        # Create and add colorbar for value mapping
        sm = plt.cm.ScalarMappable(
            cmap="plasma_r",
            norm=plt.Normalize(vmin=values.min(), vmax=values.max()),
        )
        plt.colorbar(sm, ax=ax)
    else:
        # Default visualization: gray elements with black edges
        for elm in nodes_closed:
            ax.fill(*elm.T, facecolor="lightgray", edgecolor="k", alpha=0.5)
            ax.plot(*elm.T, "ko", markersize=6, zorder=1)

    if title is not None:
        ax.set_title(title)

    return f, ax

To test the function, we will plot the mesh with each element annotated by its ID. Here, the ID corresponds to the element’s index within the connectivity array.

plot_mesh(mesh, np.arange(mesh.nelem), annotate=True, title="Element Indices")

(<Figure size 640x480 with 2 Axes>,
 <Axes: title={'center': 'Element Indices'}>)

With the plotting tools in place, we can now move on to computing more interesting quantities on this mesh.

You may have heard the term "Jacobian" in the context of Salvus and finite or spectral elements. Here, the Jacobian refers to a matrix that describes how vectors transform between two coordinate systems:

The reference coordinate system: defined on a standard, undeformed element (a square spanning $[-1, +1]^d$ , where $d$ is the number of spatial dimensions). In 2-D, we use $r$ and $s$ to defined the two reference axes.
The physical coordinate system: defined by the actual elements in our mesh. In 2-D, we use $x$ and $y$ to define the two physical axes (which is likely not surprising).

The Jacobian matrix tells us how to convert between these two parameterizations. At each point within a 2-D element, The Jacobian matrix

J

is

J = \begin{bmatrix} \frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\ \frac{\partial x}{\partial s} & \frac{\partial y}{\partial s} \end{bmatrix}

In summary, the Jacobian matrix describes how a small movement in the reference element (in

r

or

s

) translates to a movement in the physical element (in

x

or

y

). For our simple, un-deformed mesh, the transformation is straightforward: moving 2.0 units in the reference element corresponds to moving 1.0 unit in the physical element, so each diagonal entry of the Jacobian is 0.5. This means that distances along the axes are halved when mapping from the reference to the physical element.

We can also easily see that the off-diagonal terms of the Jacobian are both 0. This reflects the fact that scaling from the reference element to the physical element only changes the length along each axis (by a factor of 0.5), without any rotation or skew. Thus, moving along one axis in the reference element maps directly to movement along the same axis in the physical element, and there is no mixing between axes.

This intuition bring us to the correct Jacobian for this mesh:

J = \begin{bmatrix} 0.5 & 0.0 \\ 0.0 & 0.5 \end{bmatrix}

from which we can trivially compute the inverse Jacobian

J^{-1}

:

J^{-1} = \begin{bmatrix} 2.0 & 0.0 \\ 0.0 & 2.0 \end{bmatrix}

The inverse Jacobian naturally describes how a small movement in the physical element (in

x

or

y

) translates to a movement in the reference element (in

r

or

s

).

The Jacobian matrices for un-deformed axis-aligned meshes are a bit too simple, as they make it hard to build an intuition for how the off-diagonal entries of the Jacobian behave. To remedy this, we now consider a mesh rotated by counterclockwise by

45\degree

.

def rotate_mesh(
    mesh: sn.UnstructuredMesh, angle: float
) -> sn.UnstructuredMesh:
    """
    Rotate a mesh by `angle` degrees.

    Args:
        mesh: The mesh. Will be copied.
        angle: The angle in degrees. Positive angles rotate clockwise.

    Returns:
        A rotated mesh.
    """

    theta = np.deg2rad(angle)

    rotation_matrix = np.array(
        [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]
    )

    _mesh = mesh.copy()
    _mesh.points[:] = _mesh.points @ rotation_matrix
    return _mesh


mesh_r = rotate_mesh(mesh, -45.0)
plot_mesh(
    mesh_r,
    values=np.arange(mesh_r.nelem),
    annotate=True,
    title="Element Indices",
)

(<Figure size 640x480 with 2 Axes>,
 <Axes: title={'center': 'Element Indices'}>)

This is now a bit more interesting. As an exercise: can we intuit what the Jacobian for each element should be here? Well, we know each edge length in the mesh is

1

, so we have

a^2 + b^2 = c^2

, and since we have an isosceles right triangle the projection of the triangle on each of the coordinate axes is

1/\sqrt{2}

. Since the reference element is twice as big as each physical element, we have an additional factor of

0.5

, which gives us

0.5 / \sqrt{2} = 0.3535...

. From this, we can estimate the Jacobian as:

J \approx \begin{bmatrix} 0.3535& 0.3535 \\ -0.3535& 0.3535 \end{bmatrix}

Why the negative sign at

J[1, 0]

? As we move along the

s

axis in the reference element, we’re actually moving backwards along the

x

axis in the physical element. More generally, this exercise shows how the Jacobian captures the relationship between the reference and physical elements, regardless of orientation. This mapping is central to the flexibility of FEM and SEM methods, enabling them to adapt seamlessly to complex geometries.

In waveform modeling and inversion, this is especially powerful: it allows us to transform computations from a deformed physical grid — for example, one shaped by topography or refined near the seafloor — to a regular reference element where operations like differentiation and integration are far simpler. This built-in geometric mapping is a key advantage of spectral and finite element methods over traditional finite-difference approaches.

We can also interpret the Jacobian geometrically. For instance, we can recognize that:

Rows of the Jacobian are vectors that point along an element's facets,
Columns of the Jacobian are vectors that point perpendicular to an element's facets

Similarly:

Columns of the inverse Jacobian are vectors that, in the reference element, are orthogonal to lines of constant $r$ or $s$ respectively, accounting for deformation,
Rows of the inverse Jacobian are vectors that, in the reference element, point along lines of constant $r$ and $s$ respectively, accounting for deformation.

These relationships hold in general, but recall that the Jacobian is a local linearization, and for elements with non-affine mappings, it must be evaluated at multiple points to capture geometric variation accurately. Even so, this perspective provides a path towards operations such as structure-aligned smoothing and shape optimization, and show how these become simpler to achieve once a mesh is aligned to features or contours of interest.

Of course "intuiting" the Jacobian or computing it analytically only works for trivially simple cases. In practice, Salvus uses the shape functions that define the geometry of each element to numerically compute the Jacobian. This functionality is exposed in Python via the salvus.fem module below. Since the transformation from reference to physical coordinates is affine here, the Jacobian will be the same no matter where in the element we compute it.

# Compute the Jacobian at a reference coordinate for each element.
j = fem.jacobian.j(np.array([0.0, 0.0]), mesh_r.get_element_nodes())
j

array([[[ 0.35355339,  0.35355339],
        [-0.35355339,  0.35355339]],

       [[ 0.35355339,  0.35355339],
        [-0.35355339,  0.35355339]],

       [[ 0.35355339,  0.35355339],
        [-0.35355339,  0.35355339]],

       [[ 0.35355339,  0.35355339],
        [-0.35355339,  0.35355339]]])

As was demonstrated in Part 2, one of the defining characteristics of spectral-element methods compared to standard finite elements is the use of a high-order basis (typically order > 2) to represent both the solution field and the geometry within each element. In Salvus, this basis is built from Lagrange polynomials defined on the Gauss–Lobatto–Legendre (GLL) points. While there are several motivations for this choice, here we simply note that in practice the basis forms a discrete, structured grid within each element, with each grid point corresponding to a single Lagrange polynomial. These points are generally not equally spaced.

This multi-layered discretization — with both the outer elements and the internal GLL points — can be confusing to those from a finite-difference (FD) background. The confusion is compounded in Salvus by our use of "elements per wavelength" as the primary resolution parameter, whereas FD methods typically use "points per wavelength". In our context, points per wavelength refers specifically to the number of GLL points per wavelength.

To avoid ambiguity, we will use the following terminology:

Embedding: The set of points that define the vertices of an element. Quadrilaterals and hexahedra are always embedded by 4 or 8 points, respectively, and this count does not change with the order of the solution field.
GLL points: The internal grid points that define the Lagrange basis functions used to approximate the solution. Each element contains $(n+1)^d$ GLL points, where $n$ is the polynomial order and $d$ is the spatial dimension.

Some GLL points coincide with embedding points — specifically, those at the corners of an element. Visualizing the location of GLL points within a physical element is often instructive. This is done by interpolating the coordinates of the embedding to the positions of the GLL points in the reference element. In practice, an additional post-processing step may be applied to snap higher-order points to geometric features, but the interpolation procedure itself is general.

Examining the toy meshes from above with varying polynomial orders highlights an important observation: for the same embedding, different approximation orders yield different points per wavelength counts. This is why meaningful comparisons between SEM and FD resolution should be made on the GLL grid rather than on the embedding alone.

def compute_gll_points(mesh: sn.UnstructuredMesh, n: int):
    """
    Compute the physical coordinates of GLL points for all elements in a mesh.

    Maps GLL points from the reference element to physical space using
    interpolation. The GLL points define where Lagrange basis functions are
    centered in spectral element methods.

    Args:
        mesh: The mesh containing element connectivity and vertex coordinates.
        n: The polynomial order of the GLL points. Number of points per
            element dimension is (n+1).
    """

    # Get coefficients to map from mesh geometry order to GLL order.
    ic = fem.jacobian.get_interpolation_coefficients_from_order(
        n_dim=mesh.ndim, from_order=mesh.shape_order, to_order=n
    )

    # Interpolate mesh vertex x-coordinates to GLL point locations.
    gll_x = ic @ mesh.get_element_nodes()[..., None, 0]

    # Interpolate mesh vertex y-coordinates to GLL point locations.
    gll_y = ic @ mesh.get_element_nodes()[..., None, 1]

    # Stack coordinates and flatten to get all GLL points as (x,y) pairs.
    return np.stack([gll_x.ravel(), gll_y.ravel()])


# Get the GLL points for a few different orders
orders = [1, 2, 4, 7, 15, 21]
_, axs = plt.subplots(2, 3, figsize=(10, 7), sharey=True, sharex=True)

for row in range(axs.shape[0]):
    for col in range(axs.shape[1]):
        # Get the i^th index for fetching the order
        i = row * axs.shape[1] + col

        # Plot the corner nodes of the elements and the GLL points
        plot_mesh(mesh, title=f"Order: {orders[i]}", ax=axs[row, col])
        axs[row, col].plot(
            *compute_gll_points(mesh, orders[i]), "ro", markersize=2, zorder=0
        )

plt.show()

In practice, using an order of 4 tends to yield excellent results for most applications; this is also the default order that is used in Salvus on the solver-side.

Salvus uses an explicit time-stepping scheme and, as such, the Courant–Friedrichs–Lewy (CFL) condition plays a large role in determining a stable time step. The CFL number for a given mesh is computed from two quantities, namely

the length interval between neighboring points, and
the fastest velocity (i.e. wave speed) between those points.

The wave speed is relatively simple to compute (i.e., it is usually

v_p

, which is often directly stored in the model, or can be otherwise computed through a simple relation), but the distance between neighboring points can be more difficult, especially in the case of deformed meshes.

To implement a robust time-step estimator in Salvus — one that remains accurate even for highly deformed meshes — we draw on a technique from crystallography involving reciprocal crystal lattices. In particular, we use the concept of Miller indices, together with reciprocal lattice vectors, to approximate the minimum distance between GLL points in a deformed mesh. This distance is then inserted directly into the CFL equation to estimate a stable time step.

At present, this estimate is global: the element with the largest ratio

v_{\max} / h_{\min}

dictates the time step for the entire simulation. This can sometimes lead to unexpected cost implications, especially if a small, highly deformed element drives the global time step.

To illustrate the approach, we will manually reproduce Salvus’ time-step estimation procedure for a few example meshes. Note that this procedure is similar to that which is computed by the compute_time_step(...) method within Salvus' mesh.algorithms.unstructured_mesh.metrics module.

First, we focus on how the minimum distance is measured, as this is often less intuitive than the other CFL ingredient (the maximum velocity per element).

As noted above, we estimate distances in a deformed mesh using reciprocal lattice theory, which effectively creates a wavenumber-domain representation of a lattice of points. In practice, Salvus constructs a set of reciprocal lattice planes from a chosen set of Miller indices, and then estimates the minimum point spacing using these planes.

Two key observations underlie this approach, one of which was originally mentioned previously:

The columns of the inverse Jacobian can be interpreted as vectors normal to lines of constant $r$ and $s$ in physical space.
This set of vectors forms a basis for a reciprocal lattice, from which wavenumbers can be computed that quantify distances in a deformed mesh.

In the cells below, we will plot reciprocal lattice vectors for some simple meshes to build geometric intuition for this concept.

def plot_reciprocal_lattice(mesh: sn.UnstructuredMesh) -> None:
    """
    Plot reciprocal lattice vectors on top of a mesh visualization.

    Visualizes the reciprocal lattice basis vectors used in time-step
    estimation. As we only compute the Jacobians at the element centers here,
    the linearization error will be apparent for non-affine elements (those
    which can not be deformed by a simple rotation and scaling).

    Args:
        mesh: The mesh to analyze. Jacobians are computed at element centers.
    """

    j = fem.jacobian.j(np.array([0.0, 0.0]), mesh.get_element_nodes())

    # Plot vectors at the center of each element, oriented along the columns of
    # j
    for (x, y), j_elm in zip(mesh.get_element_nodes().mean(axis=1), j):
        # The columns of inverse Jacobian can be interpreted as vectors normal
        # to constant [r, s] isolines in physical space. This acts as a basis
        # for the reciprocal lattice.
        a0, a1 = np.linalg.inv(j_elm).T

        # The Miller indices to construct a diagonal lattice in a 2-D crystal
        # are [+1, +1] and [+1, -1]. Construct the two lattice vectors here.
        r0 = 1 * a0 + 1 * a1
        r1 = 1 * a0 - 1 * a1

        # Get the direction of the normal planes to r0 and r1 for plotting.
        l0 = np.array([-r0[1], r0[0]]) / np.linalg.det(np.linalg.inv(j_elm.T))
        l1 = np.array([-r1[1], r1[0]]) / np.linalg.det(np.linalg.inv(j_elm.T))

        # Add quiver plot for l0 and l1 vectors at element centers. These
        # should point along the lattice plane directions.
        plt.quiver(x, y, *l0)
        plt.quiver(x, y, *l1)

        # Plot an estimate of the lattice itself. Since we're only computing
        # the Jacobian at a single point (the element center), this is only
        # accurate for affine element mappings (i.e. constant Jacobians). It
        # can nevertheless be an interesting exercise to see how the lattice is
        # distorted in non-affine elements.
        plt.plot([x - l0[0], x, x + l0[0]], [y - l0[1], y, y + l0[1]], "r-.")
        plt.plot([x - l1[0], x, x + l1[0]], [y - l1[1], y, y + l1[1]], "m-.")


plot_mesh(mesh)
plot_reciprocal_lattice(mesh)
plt.gca().set_aspect("equal")

With the lattice vectors constructed, we can now finish the exercise of computing the minimum distance between points in a deformed mesh. Specifically how this is done for 2-D elements can be seen in the compute_minimum_distance function below, which each line annotated in order to give additional context. This is the exact same algorithm that executes in the solver upon simulation startup, and which is a key factor in determining the ever-so-important time step of a simulation.

def compute_minimum_distance(mesh: sn.UnstructuredMesh, n: int):
    """
    Compute minimum distance between GLL points using a reciprocal lattice.

    Estimates the smallest physical spacing between neighboring GLL points in a
    deformed mesh using crystallographic techniques.

    Args:
        mesh: The mesh containing element connectivity and geometry.
        n: The polynomial order.
    """
    # Get the positions of the GLL points in reference coordinates.
    gll_1d = fem.utils.gll_points_from_order(n)
    gll_nd = fem.utils.tensorized_quadrature_points(mesh.ndim, n)

    # Compute the Jacobian at each GLL point.
    j = fem.jacobian.j(gll_nd, mesh.get_element_nodes())

    # The columns of the transposed inverse jacobian are the reciprocal lattice
    # vectors.
    r_lat = np.swapaxes(np.linalg.inv(j), -1, -2)

    # Multiply the reciprocal lattice vectors by the diagonal Miller indices on
    # a quad to get the reciprocal wavevectors that encode the spacing between
    # adjacent lattice planes.
    r_wv = r_lat @ np.stack(([+1, +1], [+1, -1])).T

    # Compute the norm of each lattice vector. Get the maximum value for all
    # GLL points and lattice directions per element. This represents the
    # largest reciprocal wavevector, which in turn represents the smallest
    # physical distance between two lattice planes.
    r_wv_max_norm = np.max(np.linalg.norm(r_wv, axis=-2), axis=(1, 2))

    # Scale this wavevector by the minimum GLL spacing, and multiply by sqrt(2)
    # to adjust the scaling factor to account for edge-length.
    return np.min(np.diff(gll_1d)) * np.sqrt(2) / r_wv_max_norm

With the minimum distance calculation handled, the next step is to estimate the maximum velocities per element, which is the second half of the CFL equation. This is fortunately a lot simpler that the distance computation, and the result is often just the

v_p

velocity itself. For more complex anisotropic media, especially those which are not axis aligned, we do take a more involved approach that involves estimating the largest eigenvalue of the elastic tensor. However, this is not covered here in the interest of brevity.

For the rest of this tutorial we'll keep things simple and just assume a constant maximum velocity of 1.0 m/s. In real-life applications the heterogeneity of the model will, of course, come into play. What is important to note though is that this computation is always reduced to the maximum value per element.

With the two primary ingredients of the time-step estimation covered, the next and final step is to compute the time step itself. This needs one additional piece of information, the CFL number for the SEM, which has been found heuristically. Note that this varies per dimension. We can try to compute the time step for some meshes that we've created above, and will revisit this later in the context of more complicated domains.

def compute_time_step(
    n: int, mesh: sn.UnstructuredMesh, max_v: float | npt.NDArray = 1.0
):
    """
    Compute stable time step for explicit time integration using CFL condition.

    Args:
        n: The polynomial order.
        mesh: The mesh.
        max_v: Maximum velocity per element. Can be scalar (constant) or array
            with length matching mesh.nelem for spatially varying velocities.
    """

    cfl = {2: 0.6, 3: 0.47}
    min_dist = compute_minimum_distance(mesh, n)

    return np.asarray(cfl[mesh.ndim] * min_dist / max_v)


n = 1
plot_mesh(
    mesh,
    compute_time_step(n, mesh),
    annotate=True,
    title="Time step [s]",
)

plt.plot(*compute_gll_points(mesh, n), "ro", markersize=2)

[<matplotlib.lines.Line2D at 0x71931d16fe10>]

We can also check that the version of the time-step computation that is built into salvus gives the same result.

smallest_time_step, time_step_per_elm = metrics.compute_time_step(
    mesh,
    1.0,
    simulation_order=n,
)

plot_mesh(
    mesh,
    time_step_per_elm,
    annotate=True,
    title="Time step [s] (from salvus)",
)

plt.plot(*compute_gll_points(mesh, n), "ro", markersize=2)

[<matplotlib.lines.Line2D at 0x71931886a690>]

In this section, we move beyond the simple, manually generated meshes used so far. Our focus shifts to a two-layer mesh designed to mimic an ocean-bottom setting, introducing a more realistic geometry and material configuration.

All of the tools and procedures developed earlier — mesh plotting, GLL point visualization, Jacobian interpretation, and time-step estimation — will be applied directly to these more complex examples.

First, we define a function that generates a two-layer ocean-bottom mesh, assigning distinct velocities to each layer.

def generate_nontrivial_mesh(
    interlayer_coarsening_policy: (
        sn.layered_meshing.meshing_protocol.coarsening_policy.InterlayerCoarseningPolicy
        | None
    ) = None,
    intralayer_coarsening_policy: (
        sn.layered_meshing.meshing_protocol.coarsening_policy.IntralayerCoarseningPolicy
        | None
    ) = None,
    local_refinement_policy: (
        sn.layered_meshing.LocalRefinementPolicy | None
    ) = None,
    add_low_vel: bool = False,
):
    """
    Generate a complex two-layer ocean-bottom mesh with optional velocity
    heterogeneity.

    Creates a simple 2D mesh representing an ocean-sediment system with a
    water layer above and a sediment layer below. Includes configurable mesh
    refinement policies and an optional Gaussian velocity anomaly to
    demonstrate adaptive meshing strategies and their impact on computational
    cost.

    Args:
        interlayer_coarsening_policy: Policy for mesh transitions between
            material layers. Controls element size changes at material
            boundaries.
        intralayer_coarsening_policy: Policy for mesh refinement within
            individual layers. Can be single policy or list for layer-specific
            control.
        local_refinement_policy: Policy for targeted mesh refinement in
            specific regions. Typically used around velocity anomalies or
            interfaces.
        add_low_vel: If True, adds a Gaussian low-velocity anomaly in the
            sediment layer to create challenging time-step conditions.
    """

    # Define spatial grid for velocity field interpolation x-axis spans 5 km
    # horizontally, y-axis spans 2.5 km vertically.
    x, y = np.linspace(0, 5_000, 101), np.linspace(-2500.0, 0.0, 101)

    # Create 2D coordinate meshes for field evaluation Use 'ij' indexing so
    # xx[i,j] corresponds to x[i], yy[i,j] corresponds to y[j].
    xx, yy = np.meshgrid(x, y, indexing="ij")

    # Define the center of our anomaly.
    cx, cy = 2500.0, -1750.0

    # Create velocity perturbation that results in a challenging low-velocity
    # zone that will reduce time steps.
    anom = np.where(
        (np.abs(np.hypot((xx - cx), yy - cy) - 500.0) <= np.max(np.diff(x))),
        -2000.0,
        0.0,
    )

    # Define computational domain.
    domain = sn.domain.dim2.BoxDomain.from_bounds((0, 5_000), (-2_500, 0))

    # Construct layered velocity model with water-sediment system.
    model = sn.layered_meshing.LayeredModel(
        [
            # Upper layer: seawater with typical oceanic properties
            # Density 1050 kg/m³, P-wave velocity 1500 m/s.
            sn.material.from_params(rho=1050.0, vp=1500.0),
            # Interface between water and sediment layers. Creates linear-ramp
            # seafloor with 500m elevation variation.
            sn.layered_meshing.interface.Curve.from_points(
                [0.0, 5_000.0],
                [-250.0, +250.0],
                axis="x",
                reference_elevation=-500.0,
            ),
            # Lower layer: sediment with higher velocity and optional anomaly
            # density 2000 kg/m³, base P-wave velocity 3000 m/s.
            sn.material.from_params(
                rho=2000.0,
                # Velocity field: base velocity + optional perturbation.
                vp=xr.DataArray(
                    3000 + (anom if add_low_vel else 0.0),
                    [("x", x), ("v", y)],
                ),
            ),
        ]
    )

    # Configure local refinement if policy is provided; Local refinement
    # targets specific regions (e.g., around velocity anomalies).
    if local_refinement_policy is not None:
        # Extract oracle filter and refinement policy from user-provided
        # function. Oracle filter masks the minimum velocity in the meshing
        # process; local refinement policy specifies where to refine.
        of, lrp = local_refinement_policy(domain=domain, layered_model=model)
    else:
        # Use default behavior: no local refinement
        of, lrp = sn.layered_meshing.filters.no_filter, None

    # Generate the final mesh using Salvus layered meshing system.
    return sn.layered_meshing.mesh_from_domain(
        domain=domain,
        # Meshing protocol encapsulates all mesh generation parameters.
        model=sn.layered_meshing.MeshingProtocol(
            # The layered material model defined above.
            lm=model,
            # Control mesh transitions between different material layers.
            interlayer_coarsening_policy=interlayer_coarsening_policy,
            # Control mesh refinement within individual layers.
            intralayer_coarsening_policy=intralayer_coarsening_policy,
            # Apply local refinement in targeted regions.
            local_refinement_policy=lrp,
            # Oracle filter determines refinement locations.
            oracle_filter=of,
        ),
        mesh_resolution=sn.MeshResolution(
            reference_frequency=10.0,
            elements_per_wavelength=1.0,
        ),
    )

We'll also create a simple function that allows us to compute a "cost" factor to make it easier to judge how well different policies interact. This is just a simple division of the number of elements by the time step, so a higher number means that the mesh will be more expensive to simulate. We take the natural logarithm of this to make the numbers a bit more manageable.

def compute_cost_factor(n_elm: int, time_step: float) -> float:
    """
    Compute computational cost factor for mesh and time-step combination.

    Estimates the relative computational cost by combining number of elements
    with time-stepping constraints. Uses logarithmic scaling to handle the wide
    range of values encountered in practice.

    Args:
        n_elm: Number of elements in the mesh.
        time_step: Minimum stable time step for the mesh.
    """

    return np.log(n_elm / time_step)

Finally, we'll also create a small wrapper function to plot the meshes with the time step overlain, as well as print the cost factor and relevant numbers.

def plot_and_report_cost(
    mesh: sn.UnstructuredMesh,
    n: int = 1,
    show_gll: bool = True,
    title: str | None = None,
):
    """
    Visualize mesh with the time-step and report computational cost metrics.

    Args:
        mesh: The mesh to analyze and visualize.
        n: Polynomial order for GLL point computation and time-step
            calculation.
        show_gll: If True, overlay red dots showing GLL point locations.
        title: Title of the plot.  If None, no title is added.
    """

    max_vp = np.max(mesh.element_nodal_fields["VP"], axis=-1)
    _, time_step = metrics.compute_time_step(
        mesh=mesh,
        simulation_order=n,
        max_velocity=max_vp,
    )

    plot_mesh(mesh, time_step, title=title)

    if show_gll:
        plt.plot(*compute_gll_points(mesh, n), "ro", markersize=1)

    print(f"NUM ELM: {mesh.nelem}")
    print(f"MIN TIME STEP: {np.min(time_step)}")
    print(f"COST FACTOR: {compute_cost_factor(mesh.nelem, np.min(time_step))}")

Here we generate a mesh with the default settings, and just plot the velocity on top of it to get a sense for the domain.

nontrivial_mesh = generate_nontrivial_mesh()
plot_mesh(
    nontrivial_mesh,
    np.mean(nontrivial_mesh.element_nodal_fields["VP"], axis=-1),
    title="$v_p$ [m/s]",
)

(<Figure size 640x480 with 2 Axes>, <Axes: title={'center': '$v_p$ [m/s]'}>)

Now, let's look at the cost factor of this mesh to get a baseline.

plot_and_report_cost(
    nontrivial_mesh, n=1, show_gll=True, title="Time step [s]"
)

NUM ELM: 340
MIN TIME STEP: 0.014636505628749122
COST FACTOR: 10.0531821031102

Here we turn off doubling at the sea bottom and enforce a constant horizontal element size throughout the mesh. Is the cost factor surprising?

plot_and_report_cost(
    generate_nontrivial_mesh(
        interlayer_coarsening_policy=sn.layered_meshing.InterlayerConstant()
    ),
    title="Time step [s]",
)

NUM ELM: 442
MIN TIME STEP: 0.024566243134638042
COST FACTOR: 9.79769189079098

Now we add vertical refinements to accommodate the thinning of the layers.

plot_and_report_cost(
    generate_nontrivial_mesh(
        intralayer_coarsening_policy=sn.layered_meshing.IntralayerVerticalRefine()
    ),
    title="Time step [s]",
)

NUM ELM: 340
MIN TIME STEP: 0.011171871522389187
COST FACTOR: 10.323301748540702

We think here that we were too aggressive in our prescription of refinements, so remove them in the subsurface and only keep them in the ocean.

plot_and_report_cost(
    generate_nontrivial_mesh(
        intralayer_coarsening_policy=[
            sn.layered_meshing.IntralayerVerticalRefine(),
            sn.layered_meshing.IntralayerConstant(),
        ]
    ),
    n=1,
    show_gll=True,
    title="Time step [s]",
)

NUM ELM: 324
MIN TIME STEP: 0.011171871522389187
COST FACTOR: 10.275099646722825

We disable vertical refinements in the subsurface and enable a special type of vertical refinements in the ocean.

plot_and_report_cost(
    generate_nontrivial_mesh(
        intralayer_coarsening_policy=[
            sn.layered_meshing.IntralayerVerticalRefine(
                refinement_type="doubling"
            ),
            sn.layered_meshing.IntralayerConstant(),
        ]
    ),
    n=1,
    show_gll=True,
    title="Time step [s]",
)

NUM ELM: 304
MIN TIME STEP: 0.014636505628749122
COST FACTOR: 9.941264186906217

Here we examine the effect of a variety of unstructured local refinements on the time step and cost factor. For this we need to produce a somewhat more interesting model, and we insert a low-velocity anomaly in the subsurface to achieve that.

mesh_anom = generate_nontrivial_mesh(
    intralayer_coarsening_policy=[
        sn.layered_meshing.IntralayerVerticalRefine(
            refinement_type="doubling"
        ),
        sn.layered_meshing.IntralayerConstant(),
    ],
    local_refinement_policy=functools.partial(
        simple_post_refinement,
        fac=[2.0],
        restrict=[1],
        refinement_style="unstable",
    ),
    add_low_vel=True,
)

plot_mesh(mesh_anom, np.mean(mesh_anom.element_nodal_fields["VP"], axis=-1))
plot_and_report_cost(mesh_anom, n=1, show_gll=False, title="Time step [s]")

NUM ELM: 950
MIN TIME STEP: 0.006302712208122705
COST FACTOR: 11.923237213595918

As can be seen within this mesh, the time step is dictated by the constricted elements that make up the circular inclusion. Thus, in order to improve the global time step of the simulation, specialized external meshing tools may be necessary in order to further optimize this type of discretization; see the external meshes tutorial series for more information.

Here is a summary of the key takeaways of this tutorial:

Meshes are relatively simple data structures that are composed of points and connectivity.
GLL points are interpolated onto a mesh, but do not contribute to its primary structure. Nevertheless, it is the GLL points that ultimately determine the resolution and time step of a simulation.
Fully automatic quad / hex meshing of arbitrary domains and velocities is a notorious unsolved problem in computational geometry.
Some simplifications in Salvus make the problem more tractable, but it is still not a silver bullet.
Experimentation and checking the "cost factor" for a given mesh is often worthwhile.
In addition to their use in wavefield simulations, meshes can also be used for geometric operations, such as for determining the local alignment of geological substructures.

Part 2 - Averaging and integrating mesh parameters

Part 1 - Basics

© Mondaic AG (2025)

Site Map Contact Us

Academic License Agreement

PAGE CONTENTS