# Magic
%matplotlib inline
%config Completer.use_jedi = False
SalvusToolbox
to generate a parameterized model of a building.import os
import matplotlib.pyplot as plt
import numpy as np
import salvus.flow.simple_config as config
import salvus.toolbox.toolbox as st
import salvus.toolbox.toolbox_geotech as st_geo
from salvus.flow import api
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
get_simple_building
function we've imported above. The function takes a variety of values which parameterize how the mesh is generated. The values below produce a 50 "story" building with a basement foundation. As with the rest of Salvus, all units are SI.n_stories = 50
wall_width = 1.0
story_height = 3.0
ceiling_height = 0.3
building_width = 20.0
basement_depth = 20.0
basement_width = 50.0
f_max = 500.0
vs_min = 1500.0
nelem_per_wavelength = 2.0
(mesh, bnd)
. The first entry here is indeed our mesh -- we'll plot it in a few cells -- and the bnd
value represents the minimum distance to which the mesh was extruded to attach absorbing boundary layers. This value will be helpful later when we're setting up the simulation.mesh, bnd = st_geo.get_simple_building(
f_max=f_max,
vs_min=vs_min,
n_stories=n_stories,
wall_width=wall_width,
story_height=story_height,
basement_width=basement_width,
basement_depth=basement_depth,
building_width=building_width,
ceiling_height=ceiling_height,
nelem_per_wavelength=nelem_per_wavelength,
)
vs_min
value we defined above for our shear wave velocity, and use a simple scaling relation such that the p-wave velocity is twice the s-wave velocity. These parameters of course can be changed, and can also be made heterogeneous, but one must always keep in mind the resolution criterion (here: 2 elements per wavelength given a minimum velocity of m/s) in mind.nodes = mesh.get_element_nodes()[:, :, 0]
mesh.attach_field("RHO", np.ones_like(nodes) * 1000)
mesh.attach_field("VS", np.ones_like(nodes) * vs_min)
mesh.attach_field("VP", np.ones_like(nodes) * vs_min * 2)
mesh
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x74b4eed2b290>
stf = config.stf.Ricker(center_frequency=f_max / 2)
source = config.source.cartesian.VectorPoint2D(
x=-20, y=0.0, fx=1.0, fy=1.0, source_time_function=stf
)
n_rec = 100
x_rec = np.ones(n_rec) * -building_width / 2
y_rec = np.linspace(0, n_stories * story_height, n_rec)
receivers = [
config.receiver.cartesian.Point2D(
x=x,
y=y,
station_code=f"{_i:03d}",
fields=["velocity", "gradient-of-displacement"],
)
for _i, (x, y) in enumerate(zip(x_rec, y_rec))
]
# Initialize with sources and receivers.
w = config.simulation.Waveform(mesh=mesh, sources=source, receivers=receivers)
# Set end time, start time is automatically determined.
w.physics.wave_equation.end_time_in_seconds = 0.1
# Define and attach our absorbing boundaries.
ab = config.boundary.Absorbing(
side_sets=["x0", "x1", "y0"], width_in_meters=bnd, taper_amplitude=f_max
)
w.physics.wave_equation.boundaries = [ab]
# Use simplified HDF5 source output.
w.output.point_data.format = "hdf5"
# Save the volumetric wavefield for visualization purposes.
# w.output.volume_data.format = "hdf5"
# w.output.volume_data.filename = "output.h5"
# w.output.volume_data.fields = ["displacement"]
# w.output.volume_data.sampling_interval_in_time_steps = 10
# Ensure that Salvus will accept our parameters.
w.validate()
# Plot the mesh with the sources and receivers.
w
<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x74b4ece761d0>
simulation.Waveform()
object gives us a visual representation of where our sources and sensors were placed. As expected, we've instrumented the left wall with a series of velocity / gradient sensors, and placed a source just outside the building.api.run
command as below. This simulation is rather small, so running on 2 ranks locally should allow for the simulation to complete within 10 - 30 seconds. Of course for bigger simulations (or simulations in 3D) we could install Salvus on a remote supercomputing site. In this case, none of the above workflow would change, and all we would do is specify the remote site for our site_name
below. Salvus Flow would take care of shuttling data to and from the remote machine as required.api.run(
ranks=2,
get_all=True,
input_file=w,
overwrite=False,
site_name=SALVUS_FLOW_SITE_NAME,
output_folder="output",
)
SalvusJob `job_2411201339318279_9907be683a` running on `local` with 2 rank(s). Site information: * Salvus version: 2024.1.2 * Floating point size: 32
* Downloaded 17.4 MB of results to `output`. * Total run time: 8.99 seconds. * Pure simulation time: 8.06 seconds.
<salvus.flow.executors.salvus_job.SalvusJob at 0x74b4e4955510>
f, ax = plt.subplots(2, 2, figsize=(15, 10), sharex=True, sharey=True)
# Read in data.
rec_file = "output/receivers.h5"
field = "gradient-of-displacement"
uxx, dt, extent = st.get_shotgather(rec_file, field=field, cmp=0, axis=1)
uxy, _, _ = st.get_shotgather(rec_file, field=field, cmp=1, axis=1)
uyx, _, _ = st.get_shotgather(rec_file, field=field, cmp=2, axis=1)
uyy, _, _ = st.get_shotgather(rec_file, field=field, cmp=3, axis=1)
# Normalize and plot the shotgather.
uxx_min, uxx_max = 0.1 * uxx.min(), 0.1 * uxx.max()
uxy_min, uxy_max = 0.1 * uxy.min(), 0.1 * uxy.max()
uyx_min, uyx_max = 0.1 * uyx.min(), 0.1 * uyx.max()
uyy_min, uyy_max = 0.1 * uyy.min(), 0.1 * uyy.max()
# Plot the different fields.
ax[0, 0].imshow(uxx, vmin=uxx_min, vmax=uxx_max, extent=extent, aspect="auto")
ax[0, 1].imshow(uxy, vmin=uxx_min, vmax=uxx_max, extent=extent, aspect="auto")
ax[1, 0].imshow(uyx, vmin=uxx_min, vmax=uxx_max, extent=extent, aspect="auto")
ax[1, 1].imshow(uyy, vmin=uxx_min, vmax=uxx_max, extent=extent, aspect="auto")
# Label plots
ax[0, 0].set_title(r"$\partial u_x / \partial x$")
ax[0, 1].set_title(r"$\partial u_x / \partial y$")
ax[1, 0].set_title(r"$\partial u_y / \partial x$")
ax[1, 1].set_title(r"$\partial u_y / \partial y$")
# Label axes.
ax[0, 0].set_ylabel("Time (s)")
ax[1, 0].set_ylabel("Time (s)")
ax[1, 0].set_xlabel("y (m)")
ax[1, 1].set_xlabel("y (m)")
plt.show()