Mondaic - Full Waveform Solutions

This tutorial demonstrates the usage of the materials submodule of SalvusPy. It shows how to create and manipulate materials, transforming between compatible forms, and how to visualize them.

The tutorial will make extensive use of the materials submodule of SalvusPy, an addition from the 2024 version of Salvus that simplifies dealing with various forms of anisotropies, conversions between parameterizations and material orientations.

Please note that although all materials below are available in SalvusPy, for now the solver still needs either an elastic.isotropic.Velocity or elastic.triclinic.TensorComponents, material. This is because the solver currently only supports those forms and materials are not automatically converted. The notebook will end with a small example of how to convert to a solver-compatible form.

The materials module started out as a way to simplify the creation of meshes in the layered mesher. It has since grown into a full-fledged material library, allowing you to create and manipulate materials in a variety of ways.

We'll start out in this notebook with homogeneous materials, but we'll naturally progress to more complex, spatially varying materials. Arrays of stacked materials, i.e. layered models, are also a natural progression here, but those are fully handled by the layered mesher itself, and can be explored in the layered mesher tutorial.

# Import Salvus' material module
from salvus import material

# It is also exposed via the namespace module
import salvus.namespace as sn

sn.material

<module 'salvus.material' from '/miniconda/envs/salvus-py311/lib/python3.11/site-packages/salvus/material/__init__.py'>

Let's start by creating a simple homogeneous isotropic material. The material module allows you to create a material in a few different ways. The most common way is to use the from_params method, which takes a set of keywords and values.

material_1_elastic_velocity = material.from_params(
    vp=5000,
    vs=3500,
    rho=2500,
)

material_1_elastic_velocity

<IPython.core.display.HTML object>

salvus.material.elastic.isotropic.Velocity

One of the most important methods on materials is the to_wavelength_oracle method. This method is used internally by the mesher to create something that predicts the wavelength of the material (hence why we called it an "oracle") at 1 Hz.

For an homogeneous material in this parameterization, this is nothing too special:

material_1_elastic_velocity.to_wavelength_oracle()

Constant parameter
Spatially constant
Value: 3.5e+03

This method becomes a lot more powerful when considering materials that are not parametrized in velocities directly.

Let's create such a material, but this time we'll manually specify the material that we're creating:

material_2_elastic_Lame = (
    material.elastic.isotropic.LameParameters.from_params(
        lam=1.0e9,
        mu=30.0e9,
        rho=2500,
    )
)

material_2_elastic_Lame

<IPython.core.display.HTML object>

salvus.material.elastic.isotropic.LameParameters

This material is also isotropic and elastic, and has the same number of degrees of freedom as the previous one. However, it is parametrized in Lame parameters, which are a different set of parameters that are often used in geophysics. The Lame parameters are a set of two elastic constants that describe the elastic properties of a material.

Let's see what the wavelength of this material is at 1 Hz:

material_2_elastic_Lame.to_wavelength_oracle()

Constant parameter
Spatially constant
Value: 3.5e+03

That's exactly the same as the previous material, interesting.

Because these are compatible materials (i.e. they belong to the same equivalence class), we can convert between them. Every material in this module has a from_material method that takes another material as an argument and returns a material in that parameterization. In other words, we can make a material parametrized with velocities from a material parametrized with Lame parameters, and vice versa:

material_2_elastic_converted = (
    material.elastic.isotropic.Velocity.from_material(
        material_2_elastic_Lame,
    )
)

material_2_elastic_converted

<IPython.core.display.HTML object>

salvus.material.elastic.isotropic.Velocity

These materials where nearly the same! The P-wave velocity is just a tad lower for the second material, whereas the S-wave velocity seems to be exatly the same.

The little table that is displayed only has limited precision, so let's extract some of the parameter values:

material_2_elastic_converted.VP

Constant parameter
Spatially constant
Value: 4.9e+03

The output above is still a Salvus internal representation, but this time of a parameter. To get the underlying data, we can use the .p property of the parameter (VP in this case):

material_2_elastic_converted.VP.p, material_2_elastic_converted.VS.p

(4939.635614091388, 3464.1016151377544)

This shows that the underlying data in this case is a Python float, and the actual values of the material have high precision. The value for shear wave velocity is close to 3500 m/s, but not exactly that.

Since it can be cumbersome to always have to use the parameter name and .p property, we also implemented a way to retrieve all the data for a material as an XArray Dataset, which allows convenient numerical operations.

material_2_elastic_converted.ds

<xarray.Dataset> Size: 24B
Dimensions:  ()
Data variables:
    RHO      float64 8B 2.5e+03
    VP       float64 8B 4.94e+03
    VS       float64 8B 3.464e+03

There is quite the menagerie of parameterizations supported in Salvus. Your IDE (VSCode, Jupyter Notebook, ...) should be able to autocomplete any of the material submodules to allow easy exploration of the available materials, but we can also e.g. list all the available elastic and acoustic isotropic parameterizations:

print("Elastic isotropic parameterizations:")
for el_iso_mat in sn.material.elastic.isotropic.__all__:
    print(f"\t- {el_iso_mat}")

print("Acoustic isotropic parameterizations:")
for ac_iso_mat in sn.material.acoustic.isotropic.__all__:
    print(f"\t- {ac_iso_mat}")

Elastic isotropic parameterizations:
	- Velocity
	- LameParameters
	- BulkAndShearModulus
	- ShearModulusAndPoissonsRatio
	- TensorComponents
	- YoungsAndShearModulus
	- YoungsModulusAndPoissonsRatio
	- ShearModulusAndPoissonsRatio2D
	- YoungsAndShearModulus2D
	- YoungsModulusAndPoissonsRatio2D
Acoustic isotropic parameterizations:
	- BulkModulus
	- LinearCoefficients
	- Velocity

Every material within a single equivalence class (i.e. with the same degree of anisotropy) has the same number of degrees of freedom, so one can 'round-trip' between them:

(
    # This won't work if the material is not the same parameterization
    material_2_elastic_Lame
    == sn.material.elastic.isotropic.LameParameters.from_material(
        sn.material.elastic.isotropic.YoungsAndShearModulus.from_material(
            sn.material.elastic.isotropic.Velocity.from_material(
                material_2_elastic_Lame,
            )
        )
    )
)

False

The materials module also allows you to create materials that are not homogeneous, but rather spatially varying. It supports three different ways to introduce spatial variation:

XArray DataArrays, that carry information about the spatial variation in one or multiple dimensions.
SymPy expressions, that allow you to define a material in terms of a mathematical expression in one or multiple dimensions.
NumPy arrays, that allow you to define a material in terms of a NumPy array. This does not in itself carry spatial information, but can be used on NumPy arrays obtained from meshes and reattached, in which way unstructured spatial information can be introduced.

Let's have a look first at XArray DataArrays. The layered mesher supports a variety of coordinates, such as cartesian, depth, radius -- we will stick to the three cartesian coordinates for this example.

import xarray as xr

vp = xr.DataArray(
    data=[5000, 5500, 6000],
    dims=["z"],
    coords={"z": [0, 1000, 2000]},
)

# The materials module supports mixed parameter types, i.e. floats together
# with DataArrays.
vs = 3500
rho = 2500

material_3_elastic_velocity_spatial = sn.material.from_params(
    vp=vp,
    vs=vs,
    rho=rho,
)

material_3_elastic_velocity_spatial

<IPython.core.display.HTML object>

salvus.material.elastic.isotropic.Velocity

In the widget's output, we can see that the P-wave velocity is now a DataArray, whereas the S-wave velocity and density are still floats.

Salvus also takes care of broadcasting spatial information, so if we try to transform this material to a different parameterization, more of the parameters will be converted to DataArrays (as necessary).

sn.material.elastic.isotropic.BulkAndShearModulus.from_material(
    material_3_elastic_velocity_spatial
)

<IPython.core.display.HTML object>

salvus.material.elastic.isotropic.BulkAndShearModulus

The result of our wavelength oracle is now a DataArray as well, shouldn't it?

material_3_elastic_velocity_spatial.to_wavelength_oracle()

Constant parameter
Spatially constant
Value: 3.5e+03

No! This is due to the fact that the wavelength oracle is a function of the slowest wave speed in the material -- S-wave velocity, which we left as a float.

Funnily enough, if one first converts the material to a different parameterization (which will broadcast P-wave velocity's spatial information against all parameters that are dependent on it in the transformation), the wavelength oracle will be a DataArray as well:

sn.material.elastic.isotropic.YoungsAndShearModulus.from_material(
    material_3_elastic_velocity_spatial
).to_wavelength_oracle()

Discrete parameter
Coordinates: 'z' [0 to 2e+03]
Range: 3.5e+03 to 3.5e+03 - shape: (3,) - 24.0 bytes

The broadcasting is smart enough to broadcast the spatial information for multiple dimensions:

All parameters here still only have one dimension of spatial information. However, when converting to a parameterization that combines the original parameters, the resulting parameters will have the combined spatial information:

mu_retrieved = sn.material.elastic.isotropic.LameParameters.from_material(
    material_4_elastic_velocity_spatial
).MU  # This is only the internal Salvus representation of the parameter

mu_retrieved_dataarray = mu_retrieved.p  # This is the underlying DataArray

# Standard XArray plotting
mu_retrieved_dataarray.sel(z=400, method="nearest").plot()

<matplotlib.collections.QuadMesh at 0x70f43febed90>

domain = sn.domain.dim3.BoxDomain(x0=0, x1=2000, y0=0, y1=2000, z0=0, z1=2000)

mesh_material_4 = sn.layered_meshing.mesh_from_domain(
    domain=domain,
    model=material_4_elastic_velocity_spatial,
    mesh_resolution=sn.MeshResolution(reference_frequency=10),
)

mesh_material_4

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

sn.layered_meshing.mesh_from_domain(
    domain=domain,
    model=sn.material.elastic.isotropic.LameParameters.from_material(
        material_4_elastic_velocity_spatial
    ),
    mesh_resolution=sn.MeshResolution(reference_frequency=10),
)

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

We can also define materials in terms of SymPy expressions, using the same coordinates as above.

import sympy as sp

depth_symbol = sp.symbols("depth")

vp_sp = (
    2000 + 25 * sp.cos(depth_symbol / 50) + 100 * sp.exp(depth_symbol / 2000)
)

material_5_elastic_velocity_spatial = sn.material.from_params(
    vp=vp_sp, vs=1000, rho=2500
)
sn.layered_meshing.mesh_from_domain(
    domain=domain,
    model=material_5_elastic_velocity_spatial,
    mesh_resolution=sn.MeshResolution(reference_frequency=10),
)

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

Materials created using SymPy expressions transform just as easily:

sn.material.elastic.isotropic.YoungsModulusAndPoissonsRatio.from_material(
    material_5_elastic_velocity_spatial
)  # (Not that these expressions are particularly nice to look at)

<IPython.core.display.HTML object>

salvus.material.elastic.isotropic.YoungsModulusAndPoissonsRatio

The last way to create a material is to use NumPy arrays. This is not particularly useful in itself as a means of carrying spatial information, but it can be used to convert material parameters defined on a mesh as NumPy arrays:

vp_np, vs_np, rho_np = [
    mesh_material_4.element_nodal_fields[par] for par in ["VP", "VS", "RHO"]
]

vp_np.shape, vs_np.shape, rho_np.shape

((729, 8), (729, 8), (729, 8))

We can now use these NumPy arrays to create a material, convert it, and reattach it to the mesh. This allows us to convert a mesh from one system of parameters to another. Although Salvus currently supports only a limited number of parameterizations for the solver (isotropic VP/VS, TTI velocities, and full stiffness tensor Cij), future versions will take any of the defined materials and run the simulation directly with those materials.

material_6_converted_from_mesh = (
    sn.material.elastic.isotropic.YoungsModulusAndPoissonsRatio.from_material(
        sn.material.from_params(
            vp=vp_np,
            vs=vs_np,
            rho=rho_np,
        )
    )
)

mesh_material_6 = mesh_material_4.copy()
mesh_material_6.elemental_fields.clear()  # Remove old fields

for field_name, numpy_array in material_6_converted_from_mesh.ds.items():
    mesh_material_6.attach_field(name=field_name, data=numpy_array)

mesh_material_6

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

That's a rather roundabout way of converting a mesh. Luckily, we provide two methods that allow you to:

Convert a mesh's material to a different parameterization, and
Convert a mesh to a material.

from salvus.mesh.algorithms.unstructured_mesh import (
    material_operations as mesh_material_ops,
)

# To create a material from a mesh
material_7_from_mesh = mesh_material_ops.to_material(mesh_material_4)

# To directly convert a mesh to have a different material
mesh_material_4.transform_material(
    sn.material.elastic.transversely_isotropic.Velocity
)

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

This last material is the first anisotropic material we've seen in this notebook. Notebook 2 in this series will go into more detail about anisotropic materials.

Part 2 - Anisotropy and elliptic coefficients

© Mondaic AG (2025)

Site Map Contact Us

Academic License Agreement

PAGE CONTENTS