Salvus 2026.8.0

Add the ability to specify the number of wavefield snapshots to output for a surface or volume wavefield. For example, suppose that one would want to output 300 snapshots for a wavefield animations. The user can specify the approximate number of output wavefield snapshots using the number_of_snapshots argument:

p.simulations.launch(
    ...,
    extra_output_configuration={
        "volume_data": {
            "number_of_snapshots": 300,
            "fields": ["displacement"],
        },
    },
)

Write the starttime and sampling_rate attributes of ASDF receiver output as the scalar attributes the format requires.

They used to be written as one-dimensional attributes of length one. Recent versions of numpy refuse to convert those to a scalar, which made such files unreadable through pyasdf. Files written by earlier versions of Salvus remain readable with pyasdf>=0.8.3.

Suppress cython source message in IPython/Jupyter notebooks.

An exception raised when running Salvus in a notebook can result in a traceback that contains messages like: ‘Could not get source, probably due dynamically evaluated source code.’ (sic). This distracts from the actual error message and is not useful to the user. This change patches the IPython traceback formatter to suppress this message..

Fixed a bug that was causing 3D meshes to disappear in the mesh visualization widget when attempting to plot their side sets.

Add functions for performing Laplacian and centroidal smoothing on meshes. These mesh smoothing routines are available under salvus.mesh.algorithms.improve and offer a computationally inexpensive way of improving the mesh quality for discretizations that contain strongly deformed elements.

Added the scaled Jacobian as a mesh quality metric. This metric returns a value between -1.0 and 1.0 for each element in the mesh. In order for the mesh to be used in the solver, all elements must have a scaled Jacobian that is greater than 0.0.

Here’s an example of its usage:

import salvus.namespace as sn

m = sn.UnstructuredMesh(...)
scaled_jacobians = m.compute_mesh_quality("scaled_jacobian")

Fixed a bug that was causing side sets to be incorrectly assigned for Exodus meshes when some of the side sets in the .e file were unnamed.

Spherical volume models provided as Xarray datasets that have geospatial_vertical_units of km are now reliably rescaled to base SI units during interpolation.

Add a function to compute the minimum time step across multiple meshes. This is useful for ensuring that simulations are discretized on the same time axis. Here’s an example:

from salvus.mesh.algorithms.unstructured_mesh.metrics import (
    compute_minimum_time_step_across_meshes,
)

# A list of meshes for which to compute the minimum time step.
meshes = [
    sn.UnstructuredMesh(...),
    sn.UnstructuredMesh(...),
]

# Compute the minimum time step across the various meshes
time_step = compute_minimum_time_step_across_meshes(meshes=meshes)

Determine the number of threads when searching the side set elevation based on the user-defined execution policy.

The spherical branch of get_side_set_elevations located points with a k-d tree whose queries were handed every available core, regardless of the parallelism the caller had asked for. It now uses the thread count from the execution policy it is given, as the rest of the routine already did.

Ensure that the linear solids get (re-)attached when interpolating the model on a mesh chunk of a distributed mesh.

This fixes a bug where the attenuation parameters were not properly attached during the mesh building and subsequently caused the model validation to fail.

Fix an error that could have aborted an inversion when a trial model was rejected in any iteration but the first, reported as “Error during L-BFGS computation. Unexpected state of precomputed L-BFGS data.”.

Speed up hashing and caching for inversion workflows. By making the internal bookkeeping more efficient, this should help with querying events, misfits or gradients in an inversion workflow.

The speedups are most notable for inversions with many events.

Add a slider to mesh visualization widgets that allows for the size of sources and receivers to be modified.

Allow event-dependent event configuration with task chains.

Improve the logic for the built-in computation of the bandwidth for the linear solid approximation. If the frequency is not specified by the user, the automated bandwidth computation is now considering the reference frequency of the Q factor model and the mesh frequency to ensure both are covered.

Previously, this might have caused an error in constructing the linear solid for non-monotonic bounds in cases where the reference frequency of the linear solids significantly varies from the frequencies resolved by the mesh.

Significantly reduce the I/O cost of point (receiver) output during the time loop. The receivers HDF5 file is now kept open across periodic flushes instead of being re-opened and closed each time (closing forced a synchronous writeback of the whole file), and per-receiver interpolation weights are computed once instead of per flush. On receiver-dense simulations this removes orders of magnitude of write amplification.

Corrects the computation of spatial gradients on 2-D meshes in the case where the element geometry was irregular.

Add output.time_reduced_data — a new output group that computes per-DOF max/min reductions over the entire simulation time range and writes a single snapshot at the end of the run. Each entry names a field, a reduction, a required criterion (signed/abs, or a field-kind-specific functional such as vonmises, frobenius, l2, pressure, volumetric), and an optional record list (value/time/state). signed/abs on vector or tensor fields reduce every stored component independently (“componentwise”); all other criteria are scalar functionals defining a single extremal instant, optionally with the full field snapshotted at that instant (state). Works for both primary fields (displacement, velocity, phi, …) and derived fields (strain, stress, gradient-of-displacement, …) across coupled acoustic-elastic simulations. Results are written as multi-component HDF5 datasets (component order recorded as metadata) rather than per-component suffixed scalars. format (hdf5/hdf5-minimal/ hdf5-full) is required and follows the same convention as final_time_data: unless hdf5-minimal, an XDMF sidecar is written alongside the HDF5 file for direct ParaView consumption. GPU (CUDA) is not yet supported and will raise an error at setup.

Fix job array submissions to Flux sites. Recent Flux releases no longer accept abbreviated command line options, so all options are now spelled out in full.

Hardened the consistency checks around observed-vs-synthetic waveform traces so that mismatched data surfaces as an immediate, actionable error instead of silently corrupting or aborting a workflow late.

Improved the algorithm for band-limited interface reinterpolation for high-order (>1) meshes.

When a layered mesh is built at model order greater than one (or locally refined), its interface side sets are re-snapped to the true interface geometry. Aliasing of the interface topography by surrounding elements could result in pinched / inverted elements with invalid Jacobians. To remedy this, the target interface is now evaluated through a standalone GllPatch that samples the interface on the side-set facet GLL nodes and reconstructs it through those facets’ Lagrange basis. Several performance improvements now make this viable for general meshes:

• A new deduplication algorithm more aggressively deduplicates horizontal points, resulting in fewer queries to the high-order interpolation routines
• Horizontal queries of a deep mesh with vertical-only deformation can now use only the 1st order element edges for the reference coordinate calculation and enclosing element search, while still interpolating through the full high-order facets
• Distributed meshing workflows can now work with full-mesh surface patches directly, removing the need for a costly halo layer to surround each distributed mesh chunk

Together we expect these developments to result in more robust handling of meshes in the presence of aliased boundary topography, decreasing the occurrence of sporadic “invalid Jacobian” error messages.

Replace the dependency on ortools with built-in solutions.

The two integer optimization problems in the layered mesher that used it are now solved directly either via a closed-form solution of a shortest path problem or via dynamic programming. Both replacements are exact, and reproduce the previous results apart from a small number of cases in which several solutions are equally good and a different one of them may now be chosen.