Salvus 2025.1.0

SalvusCompute can now account for spatial variation of acoustic and elastic attenuation parameters (Qmu, QKappa) across elements when running on CPUs, resulting in higher fidelity modelling for attenuating media.

Receiver attachments for very large number of receivers is now significantly faster.

Added a parallel version of the get_enclosing_elements function enabling very significant speedups of a number of algorithms.

Enable coupled anisotropic simulations that mix parameterizations in tensor and velocities.

It is now possible, for instance, to have an isotropic fluid layer parameterized in VP and RHO and an anisotropic elastic medium parameterized in Cij within the same mesh.

Added a few new utility functions to help estimate spatial dispersion errors of simulations.

>>> from salvus.fem import error_estimations

# Estimate the spatial dispersion error per cycle.
>>> error_estimations.estimate_spatial_dispersion_error_per_cycle(
...     p=4,
...     elements_per_wavelength=2.0
... )
4.6692585388669274e-05

# Given an existing simulation polynomial order, and a number of elements per
# wavelength, estimate how many elements per wavelength are needed to achieve
# the same accuracy with a different polynomial order.
>>> error_estimations.estimate_epw_for_new_p(
...     original_p=4,
...     original_epw=1.5,
...     new_p=7
... )
0.6858925744742578

# Estimate the required elements per wavelength to achieve a certain max
# spatial dispersion error for a given number of cycles/propagating
# wavelengths.
# This example assumes one wants to propagate a wavefield for 100 cycles with
# a max error of 1%:
>>> error_estimations.estimate_epw_for_simulation(
...     p=4,
...     number_of_cycles=100.0,
...     acceptable_error=1e-2
... )
1.8183856099191276

The method get_data_cube has been deprecated and become more performant by internally using get_waveform_data_xarray. Passing an index to the get_data_cube is removed from the documentation, as this was never possible. Typical usage of get_data_cube should be replaced by get_waveform_data_xarray.

Added utilities to acquire JSONSchemas for the definitions of SalvusExecutors.

More stable distribution of task chains across runners on a single node.

The receiver ids are now validated to not contain duplicates upon initialization of a sn.ReceiverBlock object. This used to only raise an exception upon running a simulation - now it fails faster.

Fixing two bugs: The first prevented using certain combinations of job array sizes in Slurm while limiting the number of jobs in a job array that can run in parallel. The second causes a failure if Slurm reported non-consecutive job ids in the job array status.

Temporary files created by SalvusProject are now temporarily stored in the project file structure. This avoids issues on systems with small /tmp directories.

Improved function serialization and deserialization in Salvus.

Add DirectOrthonormalBasis and AzimuthDipRake materials to the orientation materials module, allowing unique orientation of orthotropic, monoclinic and triclinic materials.

Renamed parameters of isotropic materials in elastic and acoustic case:

• The shear modulus in all materials (except Lamé’s parameters) are now named G, to be consistent with EngineeringConstants of lower order symmetries.
• Any material that had different variations for 2D and 3D, now has 3D as the default case without any extra indications (i.e. removend the 3D appendix), while 2D materials remain special cases of those with differently named parameters, e.g. V_2D.

Added Cubic and Hexagonal engineering constants materials.

Added pretty _repr_html_() methods to the materials for Jupyter widgets.

Added pretty __repr__() methods to the materials and parameters for easy object display.

Adds some memory optimizations that improve the performance of the get_rotation_matrices function, which should be noticeable when working with large arrays of oriented materials.

Add a new module salvus.material.plugins which houses additional functionality for extending material classes. To introduce this functionality, the salvus.materials.plugins.with_approx_anisotropic_acoustic_oracle has been introduced. When passed a Material type, this plugin will overwrite the computation of anisotropic acoustic wavelength oracles with a simpler approach based on velocities along and perpendicular to the symmetry axes, as opposed to the default which solves an Eigenvalue problem on the acoustic elastic tensor. The benefit of this approximation is a significant performance gain during meshing and model partitioning, most notably for large models.

UnstructuredMesh.write_h5() can now automatically determine the numerical precision of the fields defined on the mesh and write out the mesh in the lowest possible precision instead of using the default double precision. This helps to save disk space when writing large meshes to disk and will shorten the time required to transfer meshes between local and remote sites.

To select this option, set the datatype argument to None, i.e.:

mesh.write_h5(datatype=None)

By default, double precision is used to retain existing behavior, i.e.:

import numpy as np
mesh.write_h5(datatype=np.float64)

The model precision of meshes built by and stored in a Salvus project can be configured like this:

p = sn.Project.from_domain(...)
p.config.model_precision = np.float32

This setting now defaults to single precision np.float32. Users that run the double precision Salvus binary (should only affect a small group of users) should now specifically set the configuration of their project to

p.config.model_precision = np.float64

Integer fields such as fluid and ABC_elements are now always stored at lower precision.

The write mode of meshes stored in a Salvus project can now also be configured project-wide like this:

p.config.mesh_file_mode = "minimal"

This now defaults to "minimal_fast", which should significantly reduce the file size of meshes in a project. To visualize the resulting project meshes with Paraview, please write them out like this:

p.simulations.get_mesh(...).write_h5("mesh_for_paraview.h5")

Add functions to compute volume integrals of arbitrary fields defined at the GLL points on an unstructured mesh using the Lagrangian basis.

import sympy
import salvus.namespace as sn
from salvus.mesh import layered_meshing as lm

# Polynomial order.
n = 4
domain = sn.domain.dim2.BoxDomain(x0=0, x1=1, y0=0, y1=1)

# Create a simple mesh to use as a test.
mesh = lm.mesh_from_domain(
    domain=domain,
    model=lm.material.from_params(
        # The symbol `v` is interpreted by the layered mesher as vertical,
        # agnostic of the coordinate system / dimensionality of the mesh.
        rho=sympy.symbols("v"),
        # And something funky to give the integration algorithm a run for its
        # money.
        vp=0.5 + 0.1 * sympy.sin(sympy.symbols("v") * 16),
    ),
    mesh_resolution=sn.MeshResolution(
        reference_frequency=10.0,
        elements_per_wavelength=2.0,
        model_order=n,
    ),
)

# Symbolic integral of the function f(x,y) = y, evaluated over the domain
# [0,1], [0,1]
integral_rho = sympy.integrate(
    sympy.symbols("y"),
    (sympy.symbols("x"), 0, 1),
    (sympy.symbols("y"), 0, 1),
)
integral_vp = sympy.integrate(
    0.5 + 0.1 * sympy.sin(sympy.symbols("y") * 16),
    (sympy.symbols("x"), 0, 1),
    (sympy.symbols("y"), 0, 1),
)
integrals_exact = {"RHO": float(integral_rho), "VP": float(integral_vp)}

# Now using the new methods on meshes
integrals_gll = mesh.compute_volume_integrals(fields=["RHO", "VP"])

for field in integrals_gll:
    print(
        f"Integrals of {field} using:
"
        f"	GLL points:           {integrals_gll[field]:.7e}
"
        f"	Symbolic integration: {integrals_exact[field]:.7e})"
    )

The method to extrude a side set on a 3D mesh now retains all field values and also works for higher-order meshes.

import numpy as np

import salvus.namespace as sn

domain = sn.domain.dim3.BoxDomain(
    x0=-0.5, x1=0.5, y0=-0.5, y1=0.5, z0=-0.5, z1=0.5
)

mesh = sn.layered_meshing.mesh_from_domain(
    domain=domain,
    model=sn.material.from_params(vp=1.0, rho=1.0),
    mesh_resolution=sn.MeshResolution(
        reference_frequency=5.0,
        elements_per_wavelength=2.0,
        model_order=4,
    ),
)

offsets_scalar = np.array([0, 1.0, 2.0, 3.0])
offset_direction = np.array([0, 0, -1])
offsets = offset_direction * offsets_scalar[:, None]

# Will not retain elemental nodal fields like VP
new_mesh = mesh.extrude_side_set(
    "z0", offsets=offsets
)

# Any manually creted mesh needs to be parsed for side sets manually.
new_mesh.find_side_sets(mode="cartesian")

# Show the final combination
new_mesh + mesh

To retain elemental nodal fields, one needs to average/downcast these fields to only vary per element.

Added a function that computes the normal vectors of element GLL points to a side-set bordering those elements.

import numpy as np

from salvus.mesh.algorithms.unstructured_mesh import side_sets
from salvus.mesh.algorithms._generators import spherical

mesh = spherical.cubed_sphere_chunk_3d_range(
    r_inner=25, r_outer=30, nelem=2, element_count_rad=3
).get_unstructured_mesh()

mesh.find_side_sets_enclosing(mode="spherical_chunk")

elements_r0, _ = mesh.side_sets["r0"]
mask = np.zeros((mesh.nelem,), dtype=np.bool_)

# One could set more elements in the mask to True, but since those won't border
# the side set, there won't be any normals computed.
mask[elements_r0] = True

normals = side_sets.compute_normal_vectors_towards_side_set(
    mesh, {"r0": mask}, split_per_mask=False, outward_pointing=True
)

Add an algorithm to add a side set between two elemental masks.

from salvus.mesh.algorithms.unstructured_mesh import side_sets

# Some 2- or 3-D mesh.
import salvus.namespace as sn
from salvus.mesh import layered_meshing as lm
n = 1  # polynomial order
domain = sn.domain.dim2.BoxDomain(x0=-0.5, x1=0.5, y0=-0.5, y1=0.5)

# Create a simple mesh to use as a test. Let parameters vary only in a
# single direction so that we know their gradients in advance.
mesh = lm.mesh_from_domain(
    domain=domain,
    model=lm.material.from_params(
        vp=1., rho=1.
    ),
    mesh_resolution=sn.MeshResolution(
        reference_frequency=10.0,
        elements_per_wavelength=1.0,
        model_order=n,
    ),
)

# Mask all X coordinates larger than 0.0
mask_A = mesh.get_element_centroid()[:, 0] > 0.0,
# Mask all Y coordinates smaller than 0.0
mask_B = mesh.get_element_centroid()[:, 1] < 0.0

# Create a sideset between those two masks
side_set = side_sets.create_side_set_between_masks(
    mesh,
    mask_A=mask_A,
    mask_B=mask_B,
)

# And add it to the mesh
mesh.side_sets["my_new_side_set"] = side_set

Homogeneous materials can now be directly passed to SimulationConfiguration objects:

material = sn.material.from_params(
    rho=1500, vp=1875
)

sc = sn.SimulationConfiguration(
    name="sandstone",
    model_configuration=material
    ...
)

The MeshResolution object now allows users to specify the desired elements per wavelength per dimension and per layer. For example, in a 2-D 2-layer model one could specify:

- [(1.0, 2.0), (2.0, 4.0)] to request, in the first layer, 1.0
  element per wavelength in the horizontal direction and 2.0 in the
  vertical direction, and in the second layer, 2.0 elements per
  wavelength in the horizontal direction and 4.0 in the vertical
  direction.

While in a 3-D 2-layer model one could specify:

- [(1.0, 1.0, 2.0), (2.0, 2.0, 4.0)] to request, in the first
  layer, 1.0 element per wavelength in _both_ horizontal directions
  and 2.0 in the vertical direction, and in the second layer, 2.0
  elements per wavelength in _both_ horizontal directions and 4.0
  in the vertical direction.

Broadcasting rules in 2-D, where values for not all dimensions or layers are specified, are as follows:

- 1.0 => 1.0 elements per wavelength in all dimensions for all
  layers.
- (1.0, 2.0) => 1.0 elements per wavelength in the horizontal
  direction, 2.0 in the vertical direction, for all layers.
- [(1.0, 2.0), (2.0, 4.0)] => That which is described in the 2-D
  example above, with the additional comment that the last entry
  (2.0, 4.0) will be used for all subsequent layers beyond 2.

Broadcasting rules in 3-D, where values for not all dimensions or layers are specified, are slightly different:

- 1.0 => 1.0 elements per wavelength in all dimensions for all
  layers.
- (1.0, 2.0) => 1.0 elements per wavelength in _both_ horizontal
  direction, 2.0 in the vertical direction, for all layers.
- [(1.0, 2.0), (2.0, 4.0)] => That which is described in the 3-D
  example above: note that the _horizontal_ direction is broadcast
  to become (1.0, 1.0, 2.0) and (2.0, 2.0, 4.0) respectively. The
  last entry (2.0, 4.0) will be used for all subsequent layers
  beyond 2, with the same horizontal broadcasting rules applied.

The corresponding dimensions to each of the elements per wavelength are independent of the coordinate system, and simply reference horizontal and vertical directions, respectively. Coordinate names will always be alphabetical, so:

- in 2-D: x, y
- in 3-D: x, y, z
- in 2-D spherical: latitude, radius
- in 3-D spherical: latitude, longitude, radius

Increase size of the feasible set for interlayer optimization (to from 1 to 50 times the proposed number of elements in the current layer), in response to a use case where the previous default (from 1 to 20 times the proposed number of elements in the current layer) was too small.

Adds a new experimental function interpolate_mesh_to_mesh within salvus.mesh.tools.interpolate_mesh_to_mesh which performs mesh-to-mesh interpolation between Cartesian meshes using a new optimized algorithm.

Improved the performance of higher-order mesh re-interpolation via the introduction of an approximate algorithm and other improvements. For large meshes, the new default settings can give order-of-magnitude improvements in performance for 4th-order meshes. By default, the interpolation of surface deformation to the interior of a domain is now done linearly, rather than at the model order of the mesh. This is sufficient for the majority of cases. If the resulting mesh has geometric issues (like invalid Jacobians), increasing the _interior_deformation_order argument to the MeshingProtocol to match the mesh’s model order may help create a valid mesh with smoother internal deformations. The new algorithms are only considered for Cartesian meshes.

Fixes a bug which could cause a typing error when IntralayerVerticalRefine() policies were selected.

Bugfix for the experimental periodic boundaries where a permutation of the connectivity was required when periodization occurred. This allows for the proper periodization of more complex meshes, i.e. those which are not just built from a single structured element block. Permutations for 3-D meshes are not yet implemented, and attempting to make such a mesh periodic will result in an error.

Several performance improvements were implemented during the mesh generation process which should be visible when using large model files which are mutually aligned to a common coordinate grid.

The option use_event_dependent_gradients for the trust-region method now by default uses the option “auto-detect”. When running an inversion, the appropriate setting is automatically determined depending on whether or not the forward simulations use event dependent mesh masking.

A region of interest can now optionally be provided as a callback function that takes a single mandatory input argument prior_model, corresponding to an UnstructuredMesh object and returns either an UnstructuredMesh object with a region_of_interest field or - preferred - just the numpy array of the field.

Example:

import numpy as np
import salvus.namespace as sn

def roi(prior_model):
  f = np.ones_like(prior_model.connectivity)
  f[prior_model.points[prior_model.connectivity][:, :, 0] > 0.5] = 0.0
  prior_model.attach_field("region_of_interest", f)
  return prior_model

mapping = sn.Mapping(
  scaling="absolute",
  inversion_parameters=["VP"],
  region_of_interest=roi)

Alternatively, the callback function can also directly return a numpy array of the region of interest defined as an elemental field, like this:

import numpy as np
import salvus.namespace as sn

def roi_numpy_array(prior_model):
  f = np.ones_like(prior_model.connectivity)
  f[prior_model.points[prior_model.connectivity][:, :, 0] > 0.5] = 0.0
  return f

mapping = sn.Mapping(
  scaling="absolute",
  inversion_parameters=["VP"],
  region_of_interest=roi_numpy_array)
```.

Fix a bug that prevented the use of data normalizations prior to misfit computation for misfit functions without observed data.

Ignore parameter fields outside the region of interest when scaling by relative deviations from the prior. This is relevant for inverting elastic parameters in coupled media, where some fields (e.g. VS, Cij) are non-zero only in the elastic domain.

Previously, an error would have been raised when inverting for relative changes in VS in coupled media, even if the region of interest targets only the elastic domain.

Enable gradient computations in coupled fully anisotropic elastic and isotropic acoustic media. Mixing parameterizations of Cij in the elastic domain and rho, vp in the acoustic domain was previously not supported.

Fix a bug that would lead SalvusProject to fail when attempting to retrieve jobs that could not be found in the internal database.

Custom source time functions are now automatically stored within the project. This fixes a bug where a custom source time function would sometimes not be found anymore on machines with a temporary file system because it was stored in an arbitrary directory that did not exist anymore. This is also useful when copying the project to a different location.

A project’s configuration is now no longer mutable if the project contains at least one simulation configuration.

External data stored in an xarray object can now directly be saved to a project.

import xarray as xr

ds = xr.Dataset(...)

# Adds the data + the event object to a project.
project.waveforms.add_external(
  data_name="observed_data",
  event=event.event_name,
  data=ds,
  add_event_to_project=True
)

The task chain runner workflows can now recover, even if one of the local task chain runners (e.g. the process coordinating the different runners on one node) crashes.

The MisfitConfiguration now takes an optional extra_output_callback_function argument that can be used to write out quality control information during the misfit computation.

def get_extra_output(
    *, event_misfit: sn.EventMisfit, simulation_configuration_name: str
) -> ExtraOutput | None:
    # In this example we only desire the output for "event_b"
    if event_misfit.event_name != "event_b":
        return None

    # If this event is encountered during the misfit computation
    # it will write the observed processed data to that file.
    return ExtraOutput(
        data_observed_processed=pathlib.Path("data_observed_processed.h5")
    )

mc = sn.MisfitConfiguration(
    ...
    extra_output_callback_function=get_extra_output
)

# Proceed as usual

The .from_bounds() constructor of the 2-D and 3-D BoxDomain’s has been expanded to support simplified usage:

# Domains starts at a left/bottom corner of (0/0).
d_2 = domain.dim2.BoxDomain.from_bounds(
    x=22.0, y=33.0
)
# Also possible to specify the lower bounds.
d_2d = domain.dim2.BoxDomain.from_bounds(
    x=(-2, 22.0), y=(-3.3, 33.0)
)
# In 3-D.
d_3d = domain.dim3.BoxDomain.from_bounds(
    x=22.0), y=33.0, z=44.0
)
d_3d = domain.dim3.BoxDomain.from_bounds(
    x=(-2, 22.0), y=(-3.3, 33.0), z=(-4.4, 44.0)
)

Ensure compatibility with the latest versions of scipy (1.15.0), matplotlib (3.10), and mypy (1.14).

Add support for matplotlib 3.10 and maintain backward compatibility with matplotlib >=3.5.1.

Reduce the logging output of the inversion component for reporting on ongoing tasks. Unless the ping_interval_in_seconds is explicitly set in p.inversions.iterate(), unfinished tasks will not be logged multiple times.

Add the ability to explicitly set the Courant number in the WaveformSimulationConfiguration, which will be used to compute the time step of a simulation.

Example:

wsc = WaveformSimulationConfiguration(
    end_time_in_seconds=2.0,
    courant_number=0.5,
)

Adds a Project.from_mesh to create project from mesh, so that you do not need to specify a domain first.

Add RTLD_DEEPBIND to the dynamic loader flags to protect salvus from global symbol exports done by third-party dependencies.

Preserve receiver order for the extra output and the adjoint sources.

Introduced a new SALVUS_DISABLE_NUMBA_CACHE environment variable to optionally disable the numba file cache on problematic systems. Please set it to "1" to disable it. This results in slightly longer import times of Salvus but should have no other negative effects.

Added experimental features to orient anisotropic materials in sn.materials.orientation, and methods to convert these materials to form compatible with SalvusCompute.

Added experimental support for 2D material orientations.

Interface topography and 3-D model are now handled properly within the experimental create_smoothiesem_mesh.

The latitude of longitude of a source can now be passed to the experimental create_smoothiesem_mesh. This will rotate the mesh so that it is centered at the source location.

Fix a bug where wavefield adaptive meshing would fail for model order > 1.

Adds an experimental implementation of wavefield adaptive meshes for use in spherical domains, similar to those presented in “Accelerating numerical wave propagation using wavefield adapted meshes. Part I: forward and adjoint modelling”, van Driel et. al, 2020. This experimental implementation does not yet include sources located away from the geographic poles, surface topography, or model orders >= 1.

Upgrades the experiemental implementation of periodic boundary conditions to also support 2-D meshes. Additionally, changes the required argument to periodic_axes from periodic_axes to better reflect the intended usage. Note that there are several caveats that must be considered when using this feature and it should be considered experimental. For this reason, we recommend you to reach out on the Salvus support forum if you plan to use it.

Introduces a new function mesh.utils.create_mask_from_threshold which creates a point-wise mask based on threshold bounds defined by a set of materials. For instance:

create_mask_from_threshold(
    mesh,
    sn.material.from_params(rho=1.0, vp=1.0),
    sn.material.from_params(rho=None, vp=1.5),
)

will create a masking array of shape (mesh.nelem, mesh.nodes_per_element) that is True where vp is between 1.0 and 1.5 and rho is between 1.0 and infinity. If the mesh is in an alternative parameterization to that defined by the threshold values an attempt will be made to first transform the mesh’s parameterization into that of the threshold bounds.

This functionality is useful for, among other things, selecting features in a mesh to freeze during an inversion, etc.

Isotropic elastic materials can now be parameterized in terms of:

• density, Young’s Modulus, Shear Modulus
• density, Young’s Modulus, Poisson’s Ratio
• density, Shear Modulus, Poisson’s Ratio

All of these parameterizations are available within the material.elastic.isotropic module, and can be instantiated in either 2 or 3 dimensions.

Adds a new C++ binding to compute the congruence transform of the matrix b by a, defined by aba^t. This function is handy when rotating tensor quantities, in which case a represents a rotation matrix. Support for both single and batched transforms exist, with the latter taking an optional ExecutionPolicy argument governing parallelism over the batched dimension.

This function is now used as part of the rotation routines for tilted transversely isotropic media, resulting in a significant (> 2x) speedup for large models.