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

from functools import partial
from pathlib import Path

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

import salvus_mesh.structured_grid_2D as sg2d
import salvus_mesh.structured_grid_3D as sg3d
import salvus_mesh.unstructured_mesh as um
import salvus_toolbox.toolbox as st
from salvus_flow import api
from salvus_flow import simple_config as sc

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) -> um.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 10

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

    """
    x = 2.0
    if dim == 2:

        mesh = sg2d.StructuredGrid2D.rectangle(
            nelem_x=epd, nelem_y=epd, max_x=x, max_y=x
        ).get_unstructured_mesh()

    elif dim == 3:

        mesh = sg3d.StructuredGrid3D.cube(
            nelem_x=epd, nelem_y=epd, nelem_z=epd, max_x=x, max_y=x, max_z=x
        ).get_unstructured_mesh()

    mesh.find_side_sets()
    mesh.change_tensor_order(4)
    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 = sc.stf.Ricker(center_frequency=2e3)

src_scalar_2d = sc.source.cartesian.ScalarPoint2D(
    f=1, x=1, y=1, source_time_function=stf
)
src_scalar_3d = sc.source.cartesian.ScalarPoint3D(
    f=1, x=1, y=1, z=1, source_time_function=stf
)
src_vector_2d = sc.source.cartesian.VectorPoint2D(
    fx=1, fy=1, x=1, y=1, source_time_function=stf
)
src_vector_3d = sc.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(api.run, ranks=2, get_all=True, site_name="local")

Now, onto the examples themselves.

2D

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 = sc.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.

VpRho // 2D

# 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, :])
Job `job_1910230022296643_2fc865e4cf` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 517.9 KB of results to `acoustic_rhovp`.
* Total run time: 0.90 seconds.
* Pure simulation time: 0.41 seconds.
<matplotlib.tri.tricontour.TriContourSet at 0x7f6a35d57780>
# 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)
Job `job_1910230022575223_2f26e8227a` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 517.9 KB of results to `acoustic_linear`.
* Total run time: 0.94 seconds.
* Pure simulation time: 0.45 seconds.

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}

VpVsRho // 2D

# 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, :])
Job `job_1910230022878501_a607a6328e` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 669.5 KB of results to `elastic_vpvsrho`.
* Total run time: 1.29 seconds.
* Pure simulation time: 0.44 seconds.
<matplotlib.tri.tricontour.TriContourSet at 0x7f6a35a82080>
# 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, :])
Job `job_1910230022503950_6852fc2430` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 669.5 KB of results to `elastic_lambdamurho`.
* Total run time: 0.90 seconds.
* Pure simulation time: 0.46 seconds.
<matplotlib.tri.tricontour.TriContourSet at 0x7f6a35898f98>
# 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, 1)
ax.set_aspect("equal")
t, de2 = st.visualize_wavefield_2d(output_file, "displacement")
ax.tricontourf(t, de2[-1, :])

# All parameterizations should have produced the same output.
np.testing.assert_allclose(de0, de1, atol=1e-7)
np.testing.assert_allclose(de1, de2, atol=1e-7)
Job `job_1910230022740231_7de2a6ffe2` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 669.5 KB of results to `elastic_kappamurho`.
* Total run time: 0.91 seconds.
* Pure simulation time: 0.44 seconds.

We'll generate a new partially-filled waveform simulation object for 3-D.

w = sc.simulation.Waveform()

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

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.

VpRho // 3D

# Generate the mesh.
m = get_basic_mesh(3)

# 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.
w.set_mesh(m)
w.output.volume_data.fields = ["phi"]
w.physics.wave_equation.point_source = [src_scalar_3d]

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

# Read the results.
with h5py.File(output_file, "r") as fh:
    da0 = fh["/volume/phi"][:]
Job `job_1910230022645876_0cc4eeb8a0` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
                                  
* Downloaded 65.9 MB of results to `acoustic_rhovp`.
* Total run time: 8.41 seconds.
* Pure simulation time: 7.57 seconds.
# Re-parameterize as m0 and m1.
m0, m1 = (1 / rho) / vp ** 2, 1 / rho

# Generate the mesh.
m = get_basic_mesh(3)

# 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_3d]

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

# Visualize the results.
with h5py.File(output_file, "r") as fh:
    da1 = fh["/volume/phi"][:]

# All parameterizations should have produced the same output.
da0_l2 = np.linalg.norm(da0)
da1_l2 = np.linalg.norm(da1)
np.testing.assert_allclose(da0_l2, da1_l2, atol=da0_l2 * 1e-6)
Job `job_1910230022639887_f7e9d7927a` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
                                  
* Downloaded 65.9 MB of results to `acoustic_linear`.
* Total run time: 8.11 seconds.
* Pure simulation time: 7.14 seconds.

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}

VpVsRho // 3D

# Generate the mesh.
m = get_basic_mesh(3)

# 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_3d]

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

# Visualize the results.
with h5py.File(output_file, "r") as fh:
    de0 = fh["/volume/displacement"][:]
Job `job_1910230022317171_fd931b4d8f` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
                                  
* Downloaded 89.3 MB of results to `elastic_vpvsrho`.
* Total run time: 12.67 seconds.
* Pure simulation time: 11.51 seconds.
# Generate the mesh.
m = get_basic_mesh(3)

# 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_3d]

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

# # Visualize the results.
with h5py.File(output_file, "r") as fh:
    de1 = fh["/volume/displacement"][:]
Job `job_1910230022601054_60c25053f3` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
                                  
* Downloaded 89.3 MB of results to `elastic_lambdamurho`.
* Total run time: 12.26 seconds.
* Pure simulation time: 11.20 seconds.
# Generate the mesh.
m = get_basic_mesh(3)

# 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_3d]

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

# Visualize the results.
with h5py.File(output_file, "r") as fh:
    de2 = fh["/volume/displacement"][:]

# All parameterizations should have produced the same output.
np.testing.assert_allclose(de0, de1, atol=1e-7)
np.testing.assert_allclose(de1, de2, atol=1e-7)
Job `job_1910230023520651_a1ed08c702` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.10.6
  * Floating point size: 32
                                  
* Downloaded 89.3 MB of results to `elastic_kappamurho`.
* Total run time: 12.49 seconds.
* Pure simulation time: 11.35 seconds.
PAGE CONTENTS