This documentation is not for the latest stable Salvus version.
%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")
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 = {}
return mesh
vs = 500.0
vp = 1000.0
rho = 1000.0
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
)
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
)
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
fluid
flag set to 1.# 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_2403150904860227_ed82d9ba36` running on `local` with 2 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32 -> Current Task: Time loop complete* Downloaded 529.5 KB of results to `acoustic_rhovp`. * Total run time: 1.17 seconds. * Pure simulation time: 0.64 seconds.
<matplotlib.tri._tricontour.TriContourSet at 0x7f37063fbd10>
# 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_2403150904243522_6f64c215f8` running on `local` with 2 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32 -> Current Task: Time loop complete* Downloaded 529.4 KB of results to `acoustic_linear`. * Total run time: 0.99 seconds. * Pure simulation time: 0.77 seconds.
# 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_2403150904536285_5177dfb1a3` running on `local` with 2 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32 -> Current Task: Time loop complete* Downloaded 681.1 KB of results to `elastic_vpvsrho`. * Total run time: 1.25 seconds. * Pure simulation time: 0.93 seconds.
<matplotlib.tri._tricontour.TriContourSet at 0x7f370bd1ae50>
# 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_2403150904079546_c1a22c4385` running on `local` with 2 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32 -> Current Task: Time loop complete* Downloaded 681.1 KB of results to `elastic_lambdamurho`. * Total run time: 2.34 seconds. * Pure simulation time: 2.00 seconds.
<matplotlib.tri._tricontour.TriContourSet at 0x7f373c26ff10>
# 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)
SalvusJob `job_2403150904804266_5e66dc9e41` running on `local` with 2 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32 -> Current Task: Time loop complete* Downloaded 681.2 KB of results to `elastic_kappamurho`. * Total run time: 3.79 seconds. * Pure simulation time: 1.52 seconds.
w = sn.simple_config.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
fluid
flag set to 1.# 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"][:]
SalvusJob `job_2403150904594354_0d2ee431cd` running on `local` with 2 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32
* Downloaded 66.0 MB of results to `acoustic_rhovp`. * Total run time: 7.83 seconds. * Pure simulation time: 7.30 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)
SalvusJob `job_2403150904321129_c1a101673d` running on `local` with 2 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32