This documentation is not for the latest stable Salvus version.
This tutorial demonstrates how to compute gradients in a 2-D global Earth model only using SalvusFlow and SalvusMesh. This is useful in cases where you want full control and individual runs are small enough that there is no advantage to use the persistent on-disc storage SalvusProject provides.
# Global site name for all runs
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
%matplotlib inline
from salvus.mesh.simple_mesh import basic_mesh
from salvus.flow import simple_config
from salvus.flow import api
import salvus.namespace as sn
import json
import h5py
import matplotlib.pyplot as plt
import numpy as np
import obspy
import typing
# Print all available 1-D Earth models. You can of course
# provide your own.
from salvus.mesh.models_1D import model
model.get_builtin_models()
['iasp91', 'prem_iso_no_crust', 'ak135', 'prem_crust20_global', 'prem_ani_ocean', 'ak135f', 'prem_crust20_cont', 'prem_ani_one_crust', 'csem', 'VPREMOON_noLVL', 'Moon_Weber2011', 'VPREMOON', 'prem_ani_no_crust', 'prem_ani', 'prem_iso_one_crust', 'mars_sohl', 'prem_iso', 'prem_crust20_ocean']
# The most important settings are collected here.
PERIOD = 100.0
# A good idea as the gradients are always computed
# relative to the input parameterization. This will
# just result in a much smoother looking gradient.
# Planet will also be actually round.
TENSOR_ORDER = 4
# Probably up this a bit for any final results but works
# just fine for testing purposes.
ELEMENTS_PER_WAVELENGTH = 1.0
# No latitude in 2-D. Receiver always at the surface.
SRC_LONGITUDE = 0.0
SRC_DEPTH_KM = 100.0
REC_LONGITUDE = 135.0
# It will by default select a window around this phase.
PHASE = "Pdiff"
m = basic_mesh.Circular2D()
m.basic.model = "prem_iso_no_crust"
m.basic.min_period_in_seconds = PERIOD
m.basic.elements_per_wavelength = ELEMENTS_PER_WAVELENGTH
m.advanced.tensor_order = TENSOR_ORDER
m
{'absorbing_boundaries': {'number_of_wavelengths': 0.0, 'reference_frequency': inf, 'reference_velocity': 0.0, 'side_sets': ['none']}, 'advanced': {'compression': True, 'courant_number': 0.6, 'model_spline_order': 3, 'tensor_order': 4, 'velocity_model_representation': 'element_nodes'}, 'attenuation': {'auto_band': True, 'number_of_linear_solids': 5, 'power_law_alpha': 0.0, 'power_law_reference_frequency': 1.0}, 'basic': {'elements_per_wavelength': 1.0, 'min_period_in_seconds': 100.0, 'model': 'prem_iso_no_crust'}, 'chunk': {'max_colatitude': 180.0}, 'mesh_type': 'Circular2D', 'refinement': {'hmax_refinement': 1.5, 'horizontal_max_slowness_layer_index': 0, 'refinement_bottom_up': False, 'refinement_style': 'doubling'}, 'spherical': {'min_radius': 0.0}}
Available sources: https://mondaic.com/docs/references/python_apis/salvus_flow/simple_config/source/seismology
Definitely choose a SideSet...
one as this will guarantee they are exactly relative to the local mesh surface.
Available source time functions: https://mondaic.com/docs/references/python_apis/salvus_flow/simple_config/stf
src = simple_config.source.seismology.SideSetMomentTensorPoint2D(
longitude=SRC_LONGITUDE,
depth_in_m=SRC_DEPTH_KM * 1000.0,
mrr=0.0,
mpp=0.0,
mrp=1e20,
side_set_name="r1",
source_time_function=simple_config.stf.GaussianRate(
half_duration_in_seconds=PERIOD / 2.0
),
)
Same logic with the receiver: https://mondaic.com/docs/references/python_apis/salvus_flow/simple_config/receiver/seismology
rec = simple_config.receiver.seismology.SideSetPoint2D(
longitude=REC_LONGITUDE,
station_code="AA",
side_set_name="r1",
fields=["velocity"],
)
# Create the simulation object.
w = simple_config.simulation.Waveform(
mesh=m.create_mesh(),
sources=src,
receivers=rec,
end_time_in_seconds=2000.0,
# Necessary to be able to later run an adjoint
# simulation.
store_adjoint_checkpoints=True,
)
w.validate()
w
<salvus.flow.simple_config.simulation.Waveform object at 0x7fe1c163d990>
!rm -rf output
job_forward = api.run(
input_file=w,
site_name=SALVUS_FLOW_SITE_NAME,
output_folder="output",
ranks=4,
delete_remote_files=False,
)
Job `job_2010191442075609_0ce2aa6a25` running on `local` with 4 rank(s). Site information: * Salvus version: 0.11.19 * Floating point size: 32
* Downloaded 44.9 KB of results to `output`. * Total run time: 3.87 seconds. * Pure simulation time: 3.56 seconds.
event = job_forward.get_as_event()
event.get_receiver_data(
receiver_name="XX.AA.", receiver_field="velocity"
).plot()
<Figure size 800x500 with 2 Axes>
Use ObsPy to automatically pick a window around a certain phase - but of course this can also be done manually.
from obspy.taup import TauPyModel
from obspy.geodetics import locations2degrees
# Same model as in the Salvus meshing constructor.
m = TauPyModel(model="prem")
distance = locations2degrees(0.0, SRC_LONGITUDE, 0.0, REC_LONGITUDE)
arrivals = m.get_ray_paths(
distance_in_degree=distance,
source_depth_in_km=SRC_DEPTH_KM,
phase_list=[PHASE],
)
arrivals.plot_rays()
window_center = arrivals[0].time
print(f"Distance: {distance} degrees")
print(f"Window center: {window_center:.2f} seconds")
Distance: 135.0 degrees Window center: 967.30 seconds
In the most general case Salvus supports applying temporal weights with a custom function. If processing is done this way, Salvus can later assure to do the correct thing when actually computing the adjoint source.
from salvus.project.tools.data_selection import compute_window
def temporal_weights_fct(
st: obspy.Stream,
receiver: simple_config.receiver._Base,
sources: typing.List[simple_config.source._Base],
) -> typing.Dict[str, typing.List[typing.Dict]]:
# Absolute time array.
t = st[0].times() + st[0].stats.starttime.timestamp
weights = {}
for tr in st:
component = tr.stats.channel[-1]
# Only keep a window on the vertical component.
if component != "Z":
continue
weights[component] = [
{
"values": compute_window(
t=t,
window_center=window_center,
window_width=2.0 * PERIOD,
ramp_width=PERIOD,
)
}
]
return weights
event.register_temporal_weights_function(temporal_weights_fct)
event.get_receiver_data(
receiver_name="XX.AA.", receiver_field="velocity"
).plot()
<Figure size 800x250 with 1 Axes>
Salvus comes with a couple of built-in adjoint source. You can also define your own. Please refer to Salvus' documentation for details. Here we'll write a function to compute a simple no-data cross correlation adjoint source.
Please not that you do not have to invert it in time but rather specify it in a forward time axis.
def cc_misfit(data_synthetic: np.ndarray, sampling_rate_in_hertz: float):
dt = 1.0 / sampling_rate_in_hertz
s_dt = np.gradient(data_synthetic, dt)
s_dt_2 = np.gradient(s_dt, dt)
N = (data_synthetic * s_dt_2).sum() * dt
adj_src = -1.0 * (1.0 / N * s_dt)
# The first value is the misfit which does not really have a meaning here.
return 1.0, adj_src
# Finally construct the misfit object.
misfit = sn.EventMisfit(
# This particular misfit function does
# not require data.
observed_event=None,
synthetic_event=event,
# Passed as a function here for demonstration purposes.
# One could also pass a string to use Salvus' built-in
# misfit functions.
misfit_function=cc_misfit,
receiver_field="velocity",
)
# Generate the adjoint simuluation object from the
# adjoint source object.
w_adjoint = misfit.generate_adjoint_simulation_object(
gradient_parameterization="rho-vp-vs"
)
# Finally run it to produce the gradient.
job_adjoint = api.run(
input_file=w_adjoint,
site_name=SALVUS_FLOW_SITE_NAME,
output_folder="output_adjoint",
overwrite=True,
ranks=4,
)
Job `job_2010191442989042_0c92b38048` running on `local` with 4 rank(s). Site information: * Salvus version: 0.11.19 * Floating point size: 32
* Downloaded 1.8 MB of results to `output_adjoint`. * Total run time: 6.02 seconds. * Pure simulation time: 5.82 seconds.
# Visualize the gradients. Make sure to play with the color
# scale to get a nicer result.
from salvus.mesh.unstructured_mesh import UnstructuredMesh
UnstructuredMesh.from_h5("output_adjoint/gradient.h5")
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fe1bee925d0>
# The forward job can now be deleted.
job_forward.delete()
🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/tmp/job_2010191442075609_0ce2aa6a25/volume_acoustic_volume.xdmf ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/tmp/job_2010191442075609_0ce2aa6a25/volume_elastic_volume.xdmf ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/tmp/job_2010191442075609_0ce2aa6a25/volume.h5 ... 🗑 Deleting folder /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/tmp/job_2010191442075609_0ce2aa6a25 ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/input/mesh.h5 ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/input/input.toml ... 🗑 Deleting folder /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/input ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/PID.txt ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/stdout ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/stderr ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/run_binary.sh ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/output/receivers.h5 ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/output/progress.json ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/output/meta.json ... 🗑 Deleting folder /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/output ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/SUCCESS ... 🗑 Deleting file /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25/double_fork.py ... 🗑 Deleting folder /builds/Mondaic/Tutorials/TutorialsAndIntegrationTests/SALVUS_INSTALL/SalvusFlow/run/job_2010191442075609_0ce2aa6a25 ... The job has been added to the `archive` job group.