%matplotlib inline
# This line helps with tab-completion of the simple_config objects.
# The IPython/Jupyter project default to a differnet inference based
# tab completion engine which unfortunately does not yet fully work
# with SalvusFlow. This is completely optional and a convenience
# option.
%config Completer.use_jedi = False
# Standard Python packages
import os
import toml
import numpy as np
import matplotlib.pyplot as plt
# Workflow management.
from salvus.flow import api
# Specific objects to aid in setting up simulations.
from pyasdf import ASDFDataSet
from salvus.mesh.simple_mesh import basic_mesh
import salvus.flow.simple_config as sc
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
basic_mesh
package from Salvus Mesh. This allows us to quickly generate simple layered and homogeneous meshes in 2- and 3-D. For a detailed overview of the features available in the mesher, please check out the relevant tutorials online. We'll now set some basic parameters which will describe our simulation. As described above, we've preselected some elastic parameters which will suit our purposes well.# Domain setup (m).
max_x = 2000.0 # Distance in meters.
max_y = 1000.0 # Distance in meters.
max_frequency = 25.0 # Frequency in Hz.
# Material properties.
vp = 3200.0 # P-wave velocity in m/s.
vs = 1847.5 # S-wave velocity in m/s.
rho = 2200.0 # Density in kg / m^3.
CartesianIsotropicElastic2D
to auotmatically generate and return our mesh object. Basic in-notebook visualization is available for all UnstructuredMesh
objects returned from Salvus Mesh -- an example is given in the cell below. Of course nothing really beats the fully-featured 3-D visualization options offered in a package like Paraview. To visualize Salvus meshes in this and related programs, you can write the mesh out to a file on disk. The function write_h5_tensorized_model
will write the mesh data to an HDF5 file of your choice, along with a XDMF file which provides a high-level description of the HDF5 data. It is this XDMF file which should be opened in Paraview.# Generate a mesh using the "simplified" meshing interface.
mesh = basic_mesh.CartesianHomogeneousIsotropicElastic2D(
vp=vp,
vs=vs,
rho=rho,
x_max=max_x,
y_max=max_y,
max_frequency=max_frequency,
).create_mesh()
# Visualize the mesh in the notebook.
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x72071bd977d0>
# Write the mesh to a file which can be visualized in Paraview.
# Open the associated lambs_problem_mesh.xdmf to visualize.
mesh.write_h5("lambs_problem_mesh.h5")
# Sources.
fx, fy = 0.0, -1e10 # Source components (Nm)
sx, sy = 1000.0, 500.0 # Source position (m)
stf = sc.source.stf.Ricker(center_frequency=14.5) # Center frequency in Hz.
source = sc.source.cartesian.VectorPoint2D(
x=sx, y=sy, fx=fx, fy=fy, source_time_function=stf
)
# Receivers.
nr = 5 # Number of receivers.
ry0 = 800.0 # Receiver y-value.
rx0 = 1010.0 # x-value of first receiver.
rx1 = 1410.0 # x-value of last receiver.
receivers = [
sc.receiver.cartesian.Point2D(
x=x, y=ry0, station_code=f"{_i:03d}", fields=["displacement"]
)
for _i, x in enumerate(np.linspace(rx0, rx1, nr))
]
.validate()
method which will do some quick checks to ensure that the input file will be accepted by Salvus Compute.# Set up a waveform simulation object.
sim = sc.simulation.Waveform(mesh=mesh, sources=source, receivers=receivers)
# Modify the start-time, end-time, and time-step of the simulation.
sim.physics.wave_equation.end_time_in_seconds = 0.52
sim.physics.wave_equation.time_step_in_seconds = 1e-3
sim.physics.wave_equation.start_time_in_seconds = -0.08
# sim.output.volume_data.filename = "wavefield.h5"
# sim.output.volume_data.format = "hdf5"
# sim.output.volume_data.fields = ["displacement"]
# sim.output.volume_data.sampling_interval_in_time_steps = 10
# Make sure that Salvus _Compute_ will accept the current options.
sim.validate()
api.run
command takes a few arguments worth describing:site_name
: This is an identifier which tells Flow whether you're running on your local machine, some remote cluster, or perhaps the old chess computer in your grandfather's basement. As long as Salvus has been set up correctly on the specified site all data transfers to / from the local or remote machine will happen automatically. Additionally, if a job management system is present on the remote site Flow will monitor the job queue. Here we're using the site name "local". If you gave your site a different name, please update the site_name
parameter below.input_file
: This expects the WaveformSimluation object which we generated above.ranks
: This is the number of MPI ranks the job will run on, and can range from 1 to whatever your license will allow.output_folder
: It is here where all the output files of the simulation will be copied.get_all
: Salvus can produce a lot of output! The get_all flag tells Flow that we want to fetch everthing which the simlulation creates.api.run(
input_file=sim,
site_name=SALVUS_FLOW_SITE_NAME,
ranks=4,
wall_time_in_seconds=60,
output_folder="output",
get_all=True,
)
SalvusJob `job_2411062239828273_8e518719f8` running on `local` with 4 rank(s). Site information: * Salvus version: 2024.1.2 * Floating point size: 32 -> Current Task: Time loop complete* Downloaded 51.6 KB of results to `output`. * Total run time: 0.66 seconds. * Pure simulation time: 0.47 seconds.
<salvus.flow.executors.salvus_job.SalvusJob at 0x720701d375d0>
def plot_asdf_seismograms(asdf_file, axes):
"""
A simple utility function to plot traces side-by side with a
semi-analytic solution to Lamb's problem.
:param asdf_file: ASDF datafile containing receivers.
:param axes: Axes on which to plot.
"""
for i, waveform in enumerate(asdf_file.waveforms):
for j, cmp in enumerate(["X", "Y"]):
trace = waveform.displacement.select(component=cmp)[0]
axes[j, i].plot(trace.data)
for a in axes[:, 0]:
a.set_ylabel("Displacement (m)")
for a in axes[1, :]:
a.set_xlabel("Time sample")
fig, axes = plt.subplots(2, 5, sharex="col", sharey="row", figsize=(20, 5))
with ASDFDataSet("data/force_verticale.h5", mode="r") as ex2ddir, ASDFDataSet(
"output/receivers.h5", mode="r"
) as simdata:
plot_asdf_seismograms(simdata, axes)
plot_asdf_seismograms(ex2ddir, axes)