Mondaic DocumentationIt is finally time to start the inversion and compute model updates.
A quick recap of what we have done so far:
- Part 0 defined the preferred frequency band for the inversion.
- Part 1 collected data for several earthquakes and stations all across Europe.
- Part 2 set up the domain and physics we need for creating simulated data.
- Part 3 prepared a dataset for the inversion and configured how to compare simulated seismograms to the observed data
Full-waveform inversion is a non-linear problem that can only be solved iteratively.
This entails several tasks of orchestrating the recurring simulations to create synthetic seismograms, compute misfits and sensitivities with respect to changes of the model, preconditioning and updating the model.
The scientific question is how we can formulate this as an optimization problem such that we converge toward the "true Earth" with as little computational cost as possible.
Let's get started!
import os
import os, pathlib
from frequency_band import FrequencyBand
from salvus import namespace as sn
SIMULATION_TIME_IN_SECONDS = 1200.0
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
RANKS_PER_JOB = 4
PROJECT_DIR = "project_dir_central_europe"
fband_file = pathlib.Path("./frequency_band_70_120.pkl")
fband = FrequencyBand.load(fband_file)
fbandFrequency band: 8.33 mHz - 14.3 mHz, 70.0 s - 120.0 s
FrequencyBand(min_frequency_in_hertz=0.008333333333333333, max_frequency_in_hertz=0.014285714285714285)
To safe some time, we will only iterate on a subset of events.
Make sure to use the same events for which we have selected windows in the
previous notebook.
SELECTED_EVENTS = [
"event_CRETE_GREECE_Mag_6.15_2015-04-16-18-07",
"event_ICELAND_REGION_Mag_5.75_2020-06-20-19-26",
"event_STRAIT_OF_GIBRALTAR_Mag_6.38_2016-01-25-04-22",
]# Load the project with the data from the previous parts.
p = sn.Project(path=PROJECT_DIR)
Inverse problems are "ill-posed", which makes them hard to solve.
Specifically, inverse modelling requires to carefully set up the interplay of several components in order to retrieve a reasonable model as the result of the inversion. This entails several design choices:
Data
- Which events to consider? (
events) - How to compare synthetic to observed data? (
misfit_configuration)
Model
- What is the prior model? (
prior_model) - Which parameters to invert for? (
mapping)
Optimization
- How to find a better model? (
method) - Which additional information can be used to constrain the model? (
preconditioner)
Workflow
- Where to run the jobs?
(
job_submission)
The
InverseProblemConfiguration collects these parameters to setup the inverse problem.
Often enough, it is hard to choose the best setting a-prior, which is why the inverse problem is constructed as a tree that allows to branch off an ongoing inversion with changed parameters, or to run several scenarios concurrently.Let's consider only two events at first. We will demonstrate how to add the third event to an ongoing inversion later on.
p += sn.InverseProblemConfiguration(
name=f"inv_{fband.period_band_name}",
events=SELECTED_EVENTS[:2],
misfit_configuration=f"tf_phase_misfit_{fband.period_band_name}",
wavefield_compression=sn.WavefieldCompression(
forward_wavefield_sampling_interval=10
),
prior_model=f"initial_model_{fband.period_band_name}",
mapping=sn.Mapping(
scaling="relative_deviation_from_prior",
inversion_parameters=["VSH", "VSV", "VP"],
map_to_physical_parameters={
"VPV": "VP",
"VPH": "VP",
},
source_cutout_radius_in_meters=300000,
receiver_cutout_radius_in_meters=50000,
),
method=sn.TrustRegion(initial_trust_region_linf=0.02),
preconditioner=sn.ModelDependentSmoothing(
smoothing_lengths_in_wavelengths={
"VSH": [0.2, 0.5, 0.5],
"VSV": [0.2, 0.5, 0.5],
"VP": [0.2, 0.5, 0.5],
},
reference_frequency_in_hertz=fband.max_frequency_in_hertz,
reference_model="prior",
reference_velocities={"VSH": "VSH", "VSV": "VSH", "VP": "VSH"},
),
job_submission={
"forward": sn.SiteConfig(
site_name=SALVUS_FLOW_SITE_NAME, ranks_per_job=RANKS_PER_JOB
),
"adjoint": sn.SiteConfig(
site_name=SALVUS_FLOW_SITE_NAME, ranks_per_job=RANKS_PER_JOB
),
"preconditioner": sn.SiteConfig(
site_name="local", ranks_per_job=RANKS_PER_JOB
),
},
)Let's start by adding the first iteration.
p.inversions.add_iteration(
inverse_problem_configuration=f"inv_{fband.period_band_name}"
)[2025-05-21 03:51:59,348] INFO: Adding new iteration #0.
True
The progress of the inversion can be visualized in a widget. At the moment, there is not a lot to show, obviously...
p.viz.nb.inversion(
inverse_problem_configuration=f"inv_{fband.period_band_name}"
)
Many steps within a single iteration involve expensive simulations, e.g., for computing misfits or adjoint simulations to compute gradients. This is why
Salvus operates a task-based workflow.In order to be able to closely monitor the progress,
SalvusOpt steps through an iteration, and automatically dispatches simulations whenever necessary. The function resume will return whenever SalvusOpt is waiting for other tasks to finish first. Calling it several time, will step through the iteration in sequence. The log messages inform about the current status and tasks.p.inversions.resume(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
)[2025-05-21 03:52:00,235] INFO: Resuming iteration #0. Current stage: initialize [2025-05-21 03:52:00,240] INFO: 1 new tasks have been issued. [2025-05-21 03:52:00,240] INFO: Processing task `misfit_and_gradient`
[2025-05-21 03:52:04,687] INFO: Submitting job array with 2 jobs ... [2025-05-21 03:52:04,770] INFO: Launched adjoint simulations for 2 events. Please check again to see if they are finished. [2025-05-21 03:52:04,771] INFO: Some tasks of iteration #0 are still running. Please check again later.
True
p.inversions.resume(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
)[2025-05-21 03:52:04,780] INFO: Resuming iteration #0. Current stage: check_convergence [2025-05-21 03:52:04,781] INFO: Processing task `misfit_and_gradient` [2025-05-21 03:52:04,985] INFO: Some tasks of iteration #0 are still running. Please check again later.
False
While the simulations are running, we can take a closer look at tasks within an iteration, and how they can efficiently be dealt with.
p.inversions.resume(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
)[2025-05-21 03:52:04,992] INFO: Resuming iteration #0. Current stage: check_convergence [2025-05-21 03:52:04,992] INFO: Processing task `misfit_and_gradient` [2025-05-21 03:52:05,187] INFO: Some tasks of iteration #0 are still running. Please check again later.
False
You can always visualize the current state of an iteration in a widget, where you find a lot of information about the model, misfits and gradients. The contents will be updated as you step through the iteration. The widget will be updated once gradients are available.
p.viz.nb.inversion(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
)If you feel confident and don't want to be bothered with every single task, you can also tell
Salvus to run an entire iteration at once. Note the parameter timeout_in_seconds, which will force the cell to return even if the iteration has not been completed yet, and there might still be a few simulations running in the back.
Again, you can execute the cell several times or mix it with calls to the previous one until the iteration is complete.p.inversions.iterate(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
timeout_in_seconds=1800,
ping_interval_in_seconds=100,
)[2025-05-21 03:52:05,825] INFO: Resuming iteration #0. [2025-05-21 03:52:05,825] INFO: Processing task `misfit_and_gradient` [2025-05-21 03:52:06,056] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-05-21 03:53:46,058] INFO: Processing task `misfit_and_gradient` [2025-05-21 03:53:46,229] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-05-21 03:55:26,231] INFO: Processing task `misfit_and_gradient` [2025-05-21 03:55:26,618] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-05-21 03:57:06,623] INFO: Processing task `misfit_and_gradient` [2025-05-21 03:57:07,671] INFO: Iteration 0: Number of events: 2 chi = 310.66383072917364 ||g|| = 1.2199151493577768e-05 pred = --- ared = --- norm_update = --- tr_radius = --- [2025-05-21 03:57:07,676] INFO: 1 new tasks have been issued. [2025-05-21 03:57:07,677] INFO: Processing task `preconditioner` [2025-05-21 03:57:08,116] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-05-21 03:58:48,192] INFO: Processing task `preconditioner` [2025-05-21 03:58:48,675] INFO: 1 new tasks have been issued. [2025-05-21 03:58:48,676] INFO: Processing task `misfit` [2025-05-21 03:58:48,918] INFO: Submitting job array with 2 jobs ... [2025-05-21 03:58:49,030] INFO: Launched simulations for 2 events. Please check again to see if they are finished. [2025-05-21 03:58:49,031] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-05-21 04:00:29,194] INFO: Processing task `misfit` [2025-05-21 04:00:29,365] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-05-21 04:02:09,691] INFO: Processing task `misfit`
[2025-05-21 04:02:12,845] INFO: old misfit control group: 310.66383072917364 new misfit control group: 224.70255574734372 predicted reduction control group: -43.37518858909607 actual reduction control group: -85.96127498182992 2 out of 2 event(s) improved the misfit. [2025-05-21 04:02:12,845] INFO: Model update accepted. [2025-05-21 04:02:12,845] INFO: 1 new tasks have been issued. [2025-05-21 04:02:12,846] INFO: Processing task `finalize_iteration` [2025-05-21 04:02:12,991] INFO: Successfully completed iteration #0. [2025-05-21 04:02:12,992] INFO: Adding new iteration #1.
Here are some more details on the preconditioner, and how the trust region method determines model updates.
At the end of an iteration, a new iteration with the previous model update is automatically added.
We can now visualize the entire inversion in a widget to get an overview of the evolution of models and misfits during the iterations.
p.viz.nb.inversion(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
)
Before continuing with more iterations, let's have a more detailed look at intermediate data of the iteration.
The widget above provides only a preview of the model, and there are a lot more details to discover when visualizing models and gradients in Paraview. The cell below prints the model directory of the first iteration. Copy the path and use
scp to transfer the data to your local machine.print(
p.inversions.get_iteration_directory(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
iteration_id=0,
).absolute()
)/tmp/tmpwxbjt51t/project_dir_central_europe/INVERSIONS/inv_70.0_120.0_seconds/00000
Before we explore the model in Paraview, let's continue to run the inversion in the background.
Instead of continuing with the previous set of events, we can now add the
third event. We can simply do so by taking all from the previous iteration
and adding the desired events.
Here we first add the event to the inversion at large, and subsequently, we
create a new iteration that explicitly also includes new events. By default,
the created new iteration will inherit the events selected by its parent.
p.inversions.add_events(
f"inv_{fband.period_band_name}",
SELECTED_EVENTS[-1],
)
previous_events = p.inversions.get_iteration(
f"inv_{fband.period_band_name}", 0
).events
new_events = [
SELECTED_EVENTS[-1],
]
p.inversions.add_iteration(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
parent_id=0,
events=previous_events + new_events,
)[2025-05-21 04:02:13,878] INFO: Adding new iteration #2.
True
Let's have a look at the inversion tree. You will notice that we now have two iterations that inherit from the first iteration. While iteration #1 was automatically created previously, we have now added iteration #2 which operates on the same model, but contains the three events. You can hover over the green nodes in the graph to see this info. By default,
p.inversion.iterate will always continue on the last added iteration marked by the green edge in the inversion tree.The misfit for the third event has not yet been computed, but
SalvusOpt will automatically trigger the simulation once we continue with the inversion.p.viz.nb.inversion(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
)
The first iteration was only one step in model space to reduce the misfit, but the nonlinear problem will take many iterations to converge, of course. We can loop around the
iterate function and create a few more model updates.for i in range(2):
p.inversions.iterate(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
timeout_in_seconds=3600,
ping_interval_in_seconds=120,
)[2025-05-21 04:02:14,652] INFO: Resuming iteration #2. [2025-05-21 04:02:14,655] INFO: 2 new tasks have been issued. [2025-05-21 04:02:14,655] INFO: Processing task `misfit_and_gradient` [2025-05-21 04:02:14,724] INFO: Submitting job ... [2025-05-21 04:02:14,789] INFO: Launched simulations for 1 events. Please check again to see if they are finished. [2025-05-21 04:02:14,790] INFO: Processing task `gradient` [2025-05-21 04:02:14,942] INFO: Submitting job array with 2 jobs ... [2025-05-21 04:02:14,986] INFO: Launched adjoint simulations for 2 events. Please check again to see if they are finished. [2025-05-21 04:02:14,987] INFO: Some tasks of iteration #2 are still running. Please check again later. [2025-05-21 04:04:14,989] INFO: Processing task `misfit_and_gradient`
[2025-05-21 04:04:16,915] INFO: Submitting job ... [2025-05-21 04:04:16,991] INFO: Launched adjoint simulations for 1 events. Please check again to see if they are finished. [2025-05-21 04:04:16,993] INFO: Processing task `gradient` [2025-05-21 04:04:17,130] INFO: Some tasks of iteration #2 are still running. Please check again later. [2025-05-21 04:06:17,134] INFO: Processing task `misfit_and_gradient` [2025-05-21 04:06:17,489] INFO: Processing task `gradient` [2025-05-21 04:06:17,875] INFO: Some tasks of iteration #2 are still running. Please check again later. [2025-05-21 04:08:17,882] INFO: Processing task `misfit_and_gradient` [2025-05-21 04:08:18,193] INFO: Processing task `gradient` [2025-05-21 04:08:18,602] INFO: Some tasks of iteration #2 are still running. Please check again later. [2025-05-21 04:10:18,607] INFO: Processing task `misfit_and_gradient` [2025-05-21 04:10:18,912] INFO: Processing task `gradient` [2025-05-21 04:10:19,460] INFO: Iteration 2: Number of events: 3 chi = 266.5135110524476 ||g|| = 1.2226538931377635e-05 pred = -43.37518858909607 ared = -85.96127498182992 norm_update = 9488488.336141221 tr_radius = 9488488.336141221 [2025-05-21 04:10:19,587] INFO: 1 new tasks have been issued. [2025-05-21 04:10:19,588] INFO: Processing task `preconditioner` [2025-05-21 04:10:19,836] INFO: Some tasks of iteration #2 are still running. Please check again later. [2025-05-21 04:12:19,907] INFO: Processing task `preconditioner` [2025-05-21 04:12:20,199] INFO: 1 new tasks have been issued. [2025-05-21 04:12:20,200] INFO: Processing task `misfit` [2025-05-21 04:12:20,358] INFO: Submitting job array with 3 jobs ... [2025-05-21 04:12:20,426] INFO: Launched simulations for 3 events. Please check again to see if they are finished. [2025-05-21 04:12:20,427] INFO: Some tasks of iteration #2 are still running. Please check again later. [2025-05-21 04:14:20,545] INFO: Processing task `misfit` [2025-05-21 04:14:20,636] INFO: Some tasks of iteration #2 are still running. Please check again later. [2025-05-21 04:16:20,929] INFO: Processing task `misfit`
[2025-05-21 04:16:25,995] INFO: old misfit control group: 266.51351105244765 new misfit control group: 175.1552295132188 predicted reduction control group: -64.42640113830566 actual reduction control group: -91.35828153922884 3 out of 3 event(s) improved the misfit. [2025-05-21 04:16:25,995] INFO: Model update accepted. [2025-05-21 04:16:25,996] INFO: 1 new tasks have been issued. [2025-05-21 04:16:25,996] INFO: Processing task `finalize_iteration` [2025-05-21 04:16:26,312] INFO: Successfully completed iteration #2. [2025-05-21 04:16:26,315] INFO: Adding new iteration #3. [2025-05-21 04:16:26,331] INFO: Resuming iteration #3. [2025-05-21 04:16:26,336] INFO: 1 new tasks have been issued. [2025-05-21 04:16:26,336] INFO: Processing task `gradient` [2025-05-21 04:16:26,508] INFO: Submitting job array with 3 jobs ... [2025-05-21 04:16:26,559] INFO: Launched adjoint simulations for 3 events. Please check again to see if they are finished. [2025-05-21 04:16:26,560] INFO: Some tasks of iteration #3 are still running. Please check again later. [2025-05-21 04:18:26,562] INFO: Processing task `gradient` [2025-05-21 04:18:26,722] INFO: Some tasks of iteration #3 are still running. Please check again later. [2025-05-21 04:20:26,726] INFO: Processing task `gradient` [2025-05-21 04:20:26,857] INFO: Some tasks of iteration #3 are still running. Please check again later. [2025-05-21 04:22:26,860] INFO: Processing task `gradient` [2025-05-21 04:22:27,596] INFO: Iteration 3: Number of events: 3 chi = 175.15522951321879 ||g|| = 9.515498283449166e-06 pred = -64.42640113830566 ared = -91.35828153922884 norm_update = 18976976.368379027 tr_radius = 18976976.672282442 [2025-05-21 04:22:27,775] INFO: 1 new tasks have been issued. [2025-05-21 04:22:27,775] INFO: Processing task `preconditioner` [2025-05-21 04:22:28,050] INFO: Some tasks of iteration #3 are still running. Please check again later. [2025-05-21 04:24:28,153] INFO: Processing task `preconditioner` [2025-05-21 04:24:28,465] INFO: 1 new tasks have been issued. [2025-05-21 04:24:28,466] INFO: Processing task `misfit` [2025-05-21 04:24:28,607] INFO: Submitting job array with 3 jobs ... [2025-05-21 04:24:28,677] INFO: Launched simulations for 3 events. Please check again to see if they are finished. [2025-05-21 04:24:28,677] INFO: Some tasks of iteration #3 are still running. Please check again later. [2025-05-21 04:26:28,818] INFO: Processing task `misfit`
[2025-05-21 04:26:33,396] INFO: old misfit control group: 175.1552295132188 new misfit control group: 145.53113358556269 predicted reduction control group: -20.57927644252777 actual reduction control group: -29.62409592765613 3 out of 3 event(s) improved the misfit. [2025-05-21 04:26:33,396] INFO: Model update accepted. [2025-05-21 04:26:33,397] INFO: 1 new tasks have been issued. [2025-05-21 04:26:33,397] INFO: Processing task `finalize_iteration` [2025-05-21 04:26:33,675] INFO: Successfully completed iteration #3. [2025-05-21 04:26:33,678] INFO: Adding new iteration #4.
p.viz.nb.inversion(
inverse_problem_configuration=f"inv_{fband.period_band_name}"
)
The trust-region method ensures that the cumulative misfit of all selected events goes down. This is of course only an incomplete picture of what happened. SalvusProject offers a few ways to analyze the data in detail.
One option is to have a look at the misfit reduction from the starting model to the final iteration at each receiver on a map. Keep in mind that the misfit at each receiver is still summed over all components and windows.
p.viz.nb.misfit_comparison(
reference_data=f"initial_model_{fband.period_band_name}",
other_data=[
p.inversions.get_simulation_name(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
iteration_id=4,
)
],
misfit_configuration=f"tf_phase_misfit_{fband.period_band_name}",
event=SELECTED_EVENTS[-1],
)
| ini... (reference) | inv... | Reduction inv... | |
|---|---|---|---|
| II.BORG.00 | 5.7352e+01 | 4.6660e+01 | 1.0692e+01 |
| II.ESK.00 | 6.1329e+01 | 2.6602e+01 | 3.4727e+01 |
| II.LVZ.00 | 3.6222e+01 | 2.1644e+01 | 1.4578e+01 |
| II.OBN.00 | 1.7580e+01 | 8.3885e+00 | 9.1913e+00 |
| IU.GRFO. | 5.9097e+01 | 2.3879e+01 | 3.5217e+01 |
| IU.KEV.00 | 5.9052e+01 | 1.6879e+01 | 4.2173e+01 |
| IU.KONO.00 | 3.5901e+01 | 2.4002e+01 | 1.1899e+01 |
| IU.MACI. | 5.6439e+01 | 4.1371e+01 | 1.5068e+01 |
| IU.PAB.00 | 2.5167e+00 | 1.9639e+00 | 5.5285e-01 |
<pandas.io.formats.style.Styler at 0x7951838b5890>
Another option is to look at misfit statistics in the form of histograms.
p.viz.misfit_histogram(
simulation_configuration_a=f"initial_model_{fband.period_band_name}",
simulation_configuration_b=p.inversions.get_simulation_name(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
iteration_id=4,
),
misfit_configuration=f"tf_phase_misfit_{fband.period_band_name}",
events=SELECTED_EVENTS[-1],
merge_all_components=False,
)
Last but not least the most important thing is that the waveform got better so let's now also take a look at the waveforms after the three iterations. The misfit has been reduced for all events.
And although the changes are small, we do observe a better phase match on average and on some stations the difference is significant. Check out the event
STRAIT_OF_GIBRALTAR and station GB.GAL1 for a good example.
Note that only the parts of the waveforms within the shaded windows informed the inverse problem about the model space.
Interestingly, however, as we approach better models some wiggles that were not selected also match the data better now then they did initially; see station RO.BUR31 for an example.
We could make use of this, re-pick windows and then continue with more iterations and a modified MisfitConfiguration.p.viz.nb.waveforms(
[
f"PROCESSED_DATA:{fband.period_band_name}",
f"initial_model_{fband.period_band_name}",
p.inversions.get_simulation_name(
inverse_problem_configuration=f"inv_{fband.period_band_name}",
iteration_id=4,
),
],
receiver_field="displacement",
data_selection_configuration=f"initial_selection_{fband.period_band_name}",
)
That's the end of the 4-part tutorial on continental waveform inversion.
Although we restricted the dataset to only three events, we have seen all the ingredients to automate the entire toolchain for full-waveform inversion from composing a dataset to computing model updates.
In fact, this is not the end, but rather the starting point to dive into more complicated domains and settings.
Keep in mind that there are tons of options to play around with that will help to find the best setup for the inversion.
PAGE CONTENTS