This documentation is not for the latest stable Salvus version.
import pathlib
import os
import salvus.namespace as sn
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"
# Define common filter parameters so they are consistent
# in both cases.
min_frequency_in_hertz = 1.0 / 120.0
max_frequency_in_hertz = 1.0 / 70.0
highpass_corners = 3
lowpass_corners = 3
# Project with a large 2-D box domain.
d = sn.domain.dim2.BoxDomain(x0=0, x1=1e6, y0=0, y1=1e6)
p = sn.Project.from_domain(
path=PROJECT_DIR,
domain=d,
load_if_exists=True,
)
# One simulation configuration with a Heaviside STF.
p += sn.SimulationConfiguration(
name="heaviside_stf",
tensor_order=1,
elements_per_wavelength=1.25,
min_period_in_seconds=70.0,
model_configuration=sn.ModelConfiguration(
background_model=sn.model.background.homogeneous.IsotropicElastic(
vp=2000, vs=1000.0, rho=1000
)
),
event_configuration=sn.EventConfiguration(
# Heaviside.
wavelet=sn.simple_config.stf.Heaviside(),
waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
end_time_in_seconds=600.0, time_step_in_seconds=2.0
),
),
)
# Another one with a filtered Heaviside STF. Salvus has a convenience
# function for that.
p += sn.SimulationConfiguration(
name="filtered_heaviside_stf",
tensor_order=1,
elements_per_wavelength=1.25,
min_period_in_seconds=70.0,
model_configuration=sn.ModelConfiguration(
background_model=sn.model.background.homogeneous.IsotropicElastic(
vp=2000, vs=1000.0, rho=1000
)
),
event_configuration=sn.EventConfiguration(
# Filtered Heaviside.
wavelet=sn.simple_config.stf.FilteredHeaviside(
min_frequency_in_hertz=min_frequency_in_hertz,
max_frequency_in_hertz=max_frequency_in_hertz,
end_time_in_seconds=600.0,
highpass_corners=highpass_corners,
lowpass_corners=lowpass_corners,
sampling_rate_in_hertz=0.5,
),
waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
end_time_in_seconds=600.0, time_step_in_seconds=2.0
),
),
)
# Single source + receiver.
p += sn.Event(
event_name="event",
sources=[
sn.simple_config.source.cartesian.VectorPoint2D(
x=0.5e6, y=0.5e6, fx=1e6, fy=1e6
)
],
receivers=[
sn.simple_config.receiver.cartesian.Point2D(
x=0.7e6, y=0.7e6, station_code="AA", fields=["displacement"]
)
],
)
# Launch both
for sim in ["heaviside_stf", "filtered_heaviside_stf"]:
p.simulations.launch(
sim,
p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=2,
)
p.simulations.query(block=True)
[2024-03-15 09:18:47,561] INFO: Creating mesh. Hang on. [2024-03-15 09:18:48,370] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2403150918568388_afb7b1ed54@local [2024-03-15 09:18:48,752] INFO: Submitting job ... Uploading 2 files... 🚀 Submitted job_2403150918762535_08993193d9@local
True
p.viz.nb.waveforms(["heaviside_stf"], receiver_field="displacement")
FilteredHeaviside
source time function.from salvus.project.tools.processing import (
get_bandpass_filter_processing_function,
)
p.add_to_project(
sn.processing.ProcessingConfiguration(
name="heaviside_filtered_afterwards",
data_source_name="heaviside_stf",
# Could use a custom function but this yields a fully
# compatible processing function.
processing_function=get_bandpass_filter_processing_function(
min_frequency_in_hertz=min_frequency_in_hertz,
max_frequency_in_hertz=max_frequency_in_hertz,
highpass_corners=highpass_corners,
lowpass_corners=lowpass_corners,
),
),
overwrite=True,
)
# Now this already looks much better.
p.viz.nb.waveforms(
[
"PROCESSED_DATA:heaviside_filtered_afterwards",
"SYNTHETIC_DATA:filtered_heaviside_stf",
],
receiver_field="displacement",
)
# Let's do a final numerical equality check.
import numpy as np
st1 = p.waveforms.get(
"SYNTHETIC_DATA:filtered_heaviside_stf", events=p.events.list()
)[0].get_receiver_data("XX.AA.", receiver_field="displacement")
st2 = p.waveforms.get(
"PROCESSED_DATA:heaviside_filtered_afterwards", events=p.events.list()
)[0].get_receiver_data("XX.AA.", receiver_field="displacement")
for tr1, tr2 in zip(st1, st2):
# We have to give some tolerance for numerical inaccuracies.
np.testing.assert_allclose(
tr1.data, tr2.data, rtol=1e-3, atol=1e-5 * np.abs(tr1.data).max()
)