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
from salvus.mesh import mesh_block
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
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
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 = mesh_block.generators.cartesian.rectangle_2d(
nelem_x=epd, nelem_y=epd, max_x=x, max_y=x
).get_unstructured_mesh()
elif dim == 3:
mesh = mesh_block.generators.cartesian.cube_3d(
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
vs = 500.0
vp = 1000.0
rho = 1000.0
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
)
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=SALVUS_FLOW_SITE_NAME
)
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
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, :])
Job `job_2204302001384311_8b9c92bee4` running on `local_cuda` with 2 rank(s). Site information: * Salvus version: 0.12.0 * Floating point size: 32
* Downloaded 528.9 KB of results to `acoustic_rhovp`. * Total run time: 2.14 seconds. * Pure simulation time: 1.02 seconds.
<matplotlib.tri.tricontour.TriContourSet at 0x7f678ccccad0>
# 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_2204302001649015_068cad9018` running on `local_cuda` with 2 rank(s). Site information: * Salvus version: 0.12.0 * Floating point size: 32
* Downloaded 528.8 KB of results to `acoustic_linear`. * Total run time: 1.69 seconds. * Pure simulation time: 1.35 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, :])
Job `job_2204302001592267_c3d55c3735` running on `local_cuda` with 2 rank(s). Site information: * Salvus version: 0.12.0 * Floating point size: 32
* Downloaded 679.9 KB of results to `elastic_vpvsrho`. * Total run time: 1.34 seconds. * Pure simulation time: 1.03 seconds.
<matplotlib.tri.tricontour.TriContourSet at 0x7f678ae12ed0>
# 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_2204302001186179_6080328977` running on `local_cuda` with 2 rank(s). Site information: * Salvus version: 0.12.0 * Floating point size: 32
* Downloaded 679.9 KB of results to `elastic_lambdamurho`. * Total run time: 1.64 seconds. * Pure simulation time: 1.27 seconds.
<matplotlib.tri.tricontour.TriContourSet at 0x7f678ccd9310>
# 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_2204302001189713_e0f15969ad` running on `local_cuda` with 2 rank(s). Site information: * Salvus version: 0.12.0 * Floating point size: 32
* Downloaded 679.9 KB of results to `elastic_kappamurho`. * Total run time: 2.15 seconds. * Pure simulation time: 1.71 seconds.
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
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"][:]
Job `job_2204302001087285_890c75db9f` running on `local_cuda` with 2 rank(s). Site information: * Salvus version: 0.12.0 * Floating point size: 32
* Downloaded 66.0 MB of results to `acoustic_rhovp`. * Total run time: 4.98 seconds. * Pure simulation time: 4.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_2204302001642082_e45d57f6ba` running on `local_cuda` with 2 rank(s). Site information: * Salvus version: 0.12.0 * Floating point size: 32
* Downloaded 66.0 MB of results to `acoustic_linear`. * Total run time: 6.03 seconds. * Pure simulation time: 4.69 seconds.
# 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_2204302001186122_9985d6bb1c` running on `local_cuda` with 2 rank(s).
Site in