Version:
This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

# Model Parameterizations in Salvus

In this notebook we will explore all the different options to parameterize the material properties, which are currently supported in Salvus. As always, the first step is to import the requisite Python modules.
Copy
%matplotlib inline
import os
from functools import partial
from pathlib import Path

import h5py
import matplotlib.pyplot as plt
import numpy as np

import salvus.namespace as sn
import salvus.toolbox.toolbox as st

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
First we set up a simple function to return simple meshes for testing in both 2- and 3-D dimensions. The size of the mesh is small as these examples will only run for a trivial duration, and are just meant to be used to explore the possible parameterizations.
def get_basic_mesh(dim: int, epd: int = 20) -> sn.UnstructuredMesh:
"""Get a simple mesh to outline allowed parameter types.

Parameters
----------
dim : int
Dimension of the mesh.
epd : int, optional
Elements per dimension, by default 20

Returns
-------
um.UnstructuredMesh
An unstructured mesh free of any parameters.

"""
x = 2.0
h = x / float(epd)
vp = 1000.0

if dim == 2:
mesh = sn.simple_mesh.basic_mesh.CartesianHomogeneousAcoustic2D(
x_max=x,
y_max=x,
rho=1000.0,
vp=vp,
max_frequency=0.5 * vp / h,
tensor_order=4,
).create_mesh()

elif dim == 3:
mesh = sn.simple_mesh.basic_mesh.CartesianHomogeneousAcoustic3D(
x_max=x,
y_max=x,
z_max=x,
rho=1000.0,
vp=vp,
max_frequency=0.5 * vp / h,
tensor_order=4,
).create_mesh()

# Delete the material properties as they will be added later on.
mesh.elemental_fields.clear()
return mesh
Second, we will define some basic parameter values to use (all SI units).
vs = 500.0
vp = 1000.0
rho = 1000.0
Now we will define a set of sources which can be used in all of our testing environments. We'll use a Ricker wavelet in all cases, but we'll change the spatial type of our source depending on the dimension of the problem and the physics.
stf = sn.simple_config.stf.Ricker(center_frequency=2e3)

src_scalar_2d = sn.simple_config.source.cartesian.ScalarPoint2D(
f=1, x=1, y=1, source_time_function=stf
)
src_scalar_3d = sn.simple_config.source.cartesian.ScalarPoint3D(
f=1, x=1, y=1, z=1, source_time_function=stf
)
src_vector_2d = sn.simple_config.source.cartesian.VectorPoint2D(
fx=1, fy=1, x=1, y=1, source_time_function=stf
)
src_vector_3d = sn.simple_config.source.cartesian.VectorPoint3D(
fx=1, fy=1, fz=1, x=1, y=1, z=1, source_time_function=stf
)
An finally we'll partially fill in Salvus Flow's api.run function. We'll be running the function many times with mostly the same setup, so this is a nice space saving method.
run_salvus = partial(
sn.api.run, ranks=2, get_all=True, site_name=SALVUS_FLOW_SITE_NAME
)
Now, onto the examples themselves.

## 2D domains

We'll start with all of the 2-D parameterizations. Again, to save space, we'll prepare our simulation.Waveform object with values which will be re-used.
w = sn.simple_config.simulation.Waveform()

w.domain.dimension = 2
w.output.volume_data.format = "hdf5"
w.output.volume_data.filename = "output.h5"
w.output.volume_data.sampling_interval_in_time_steps = 100

### Acoustic

Acoustic meshes can be parameterized either using p-velocity and density, or a linear parameterization using M0 and M1. Acoustic elements must have the fluid flag set to 1.

#### Velocity (sound speed) and density

# Generate the mesh.
m = get_basic_mesh(2)

# Attach parameter to the nodes of each element.
par_template = np.ones_like(m.get_element_nodes()[:, :, 0])
m.attach_field("VP", par_template * vp)
m.attach_field("RHO", par_template * rho)
m.attach_field("fluid", np.ones(m.nelem))

# Attach the mesh and set some custom output.
w.set_mesh(m)
w.output.volume_data.fields = ["phi"]
w.physics.wave_equation.point_source = [src_scalar_2d]

# Run the solver.
output_folder = Path("acoustic_rhovp")
output_file = output_folder / "output.h5"
run_salvus(input_file=w, output_folder=output_folder)

# Visualize the results.
f, ax = plt.subplots(1, 1)
ax.set_aspect("equal")
t, da0 = st.visualize_wavefield_2d(output_file, "phi")
ax.tricontourf(t, da0[-1, :])
SalvusJob job_2405151229252685_cdad9a8d24 running on local with 2 rank(s).
Site information:
* Salvus version: 2024.1.0
* Floating point size: 32
-> Current Task: Time loop complete* Downloaded 529.4 KB of results to acoustic_rhovp.
* Total run time: 1.10 seconds.
* Pure simulation time: 0.51 seconds.

<matplotlib.tri._tricontour.TriContourSet at 0x7faf455ac710>

#### (Linear) compliance parameters

# Re-parameterize as m0 and m1.
m0, m1 = (1 / rho) / vp**2, 1 / rho

# Generate the mesh.
m = get_basic_mesh(2)

# Attach parameter to the nodes of each element.
par_template = np.ones_like(m.get_element_nodes()[:, :, 0])
m.attach_field("M0", par_template * m0)
m.attach_field("M1", par_template * m1)
m.attach_field("fluid", np.ones(m.nelem))

# Attach the mesh and set some custom output.
w.set_mesh(m)
w.output.volume_data.fields = ["phi"]
w.physics.wave_equation.point_source = [src_scalar_2d]

# Run the solver.
output_folder = Path("acoustic_linear")
output_file = output_folder / "output.h5"
run_salvus(input_file=w, output_folder=output_folder)

# Visualize the results.
f, ax = plt.subplots(1, 1)
ax.set_aspect("equal")
t, da1 = st.visualize_wavefield_2d(output_file, "phi")
ax.tricontourf(t, da1[-1, :])

# All parameterizations should have produced the same output.
np.testing.assert_allclose(da0, da1, atol=1e-3)
SalvusJob job_2405151229673676_4154dba8e9 running on local with 2 rank(s).
Site information:
* Salvus version: 2024.1.0
* Floating point size: 32
-> Current Task: Time loop complete* Downloaded 529.5 KB of results to acoustic_linear.
* Total run time: 1.27 seconds.
* Pure simulation time: 0.57 seconds.


### Elastic

Isotropic elastic models can be parameterized using either:
• velocities {VP, VS, RHO}
• Lame's first parameter, shear modulus, and density {LAMBDA, MU, RHO}
• Bulk modulus, shear modulus, and density {KAPPA, MU, RHO}

#### Velocities and density

# Generate the mesh.
m = get_basic_mesh(2)

# Attach parameter to the nodes of each element.
par_template = np.ones_like(m.get_element_nodes()[:, :, 0])
m.attach_field("VP", par_template * vp)
m.attach_field("VS", par_template * vs)
m.attach_field("RHO", par_template * rho)
m.attach_field("fluid", np.zeros(m.nelem))

# Attach the mesh and set some custom output.
w.set_mesh(m)
w.output.volume_data.fields = ["displacement"]
w.physics.wave_equation.point_source = [src_vector_2d]

# Run the solver.
output_folder = Path("elastic_vpvsrho")
output_file = output_folder / "output.h5"
run_salvus(input_file=w, output_folder=output_folder)

# Visualize the results.
f, ax = plt.subplots(1, 1)
ax.set_aspect("equal")
t, de0 = st.visualize_wavefield_2d(output_file, "displacement")
ax.tricontourf(t, de0[-1, :])
SalvusJob job_2405151229430823_b90b92577f running on local with 2 rank(s).
Site information:
* Salvus version: 2024.1.0
* Floating point size: 32
-> Current Task: Time loop complete* Downloaded 681.1 KB of results to elastic_vpvsrho.
* Total run time: 1.22 seconds.
* Pure simulation time: 0.80 seconds.

<matplotlib.tri._tricontour.TriContourSet at 0x7faf45582150>

#### Lame parameters and density

# Generate the mesh.
m = get_basic_mesh(2)

# Attach parameter to the nodes of each element.
mu = rho * vs**2
lam = rho * vp**2 - 2 * mu
par_template = np.ones_like(m.get_element_nodes()[:, :, 0])
m.attach_field("LAMBDA", par_template * lam)
m.attach_field("MU", par_template * mu)
m.attach_field("RHO", par_template * rho)
m.attach_field("fluid", np.zeros(m.nelem))

# Attach the mesh and set some custom output.
w.set_mesh(m)
w.output.volume_data.fields = ["displacement"]
w.physics.wave_equation.point_source = [src_vector_2d]

# Run the solver.
output_folder = Path("elastic_lambdamurho")
output_file = output_folder / "output.h5"
run_salvus(input_file=w, output_folder=output_folder)

# Visualize the results.
f, ax = plt.subplots(1, 1)
ax.set_aspect("equal")
t, de1 = st.visualize_wavefield_2d(output_file, "displacement")
ax.tricontourf(t, de1[-1, :])
SalvusJob job_2405151229166600_1b4bfa0daa running on local with 2 rank(s).
Site information:
* Salvus version: 2024.1.0
* Floating point size: 32

* Downloaded 681.1 KB of results to elastic_lambdamurho.
* Total run time: 1.19 seconds.
* Pure simulation time: 0.66 seconds.

<matplotlib.tri._tricontour.TriContourSet at 0x7fafe27401d0>

#### Elastic moduli and density

# Generate the mesh.
m = get_basic_mesh(2)

# Attach parameter to the nodes of each element.
mu = rho * vs**2
kap = rho * (vp**2 - 4 / 3 * vs**2)
par_template = np.ones_like(m.get_element_nodes()[:, :, 0])
m.attach_field("KAPPA", par_template * kap)
m.attach_field("MU", par_template * mu)
m.attach_field("RHO", par_template * rho)
m.attach_field("fluid", np.zeros(m.nelem))

# Attach the mesh and set some custom output.
w.set_mesh(m)
w.output.volume_data.fields = ["displacement"]
w.physics.wave_equation.point_source = [src_vector_2d]

# Run the solver.
output_folder = Path("elastic_kappamurho")
output_file = output_folder / "output.h5"
run_salvus(input_file=w, output_folder=output_folder)

# Visualize the results.
f, ax = plt.subplots(1,