Version:

This documentation is not for the latest stable Salvus version.

This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

2-D Global Gradients without SalvusProject

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.

Copy
# 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"

Step 1: Build Mesh

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
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
{'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
    ),
)
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
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<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
                                                     
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
* 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>

Adjoint Source

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
                                  
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
* 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")
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<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.
PAGE CONTENTS