Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

This documentation is not for the latest stable Salvus version.

Continental Scale FWI Tutorial - Part 4: Iterative inversion

It is finally time to start the inversion and compute model updates.
A quick recap of what we have done so far:
  • 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!
Copy
import os

MIN_PERIOD = 70.0
MAX_PERIOD = 120.0
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"

period_band_name = f"{int(MIN_PERIOD)}_{int(MAX_PERIOD)}_seconds"
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.
from salvus import namespace as sn
p = sn.Project(path=PROJECT_DIR)

Setting up an Inverse Problem

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_{period_band_name}",
    events=SELECTED_EVENTS[:2],
    misfit_configuration=f"tf_phase_misfit_{period_band_name}",
    wavefield_compression=sn.WavefieldCompression(
        forward_wavefield_sampling_interval=10
    ),
    prior_model=f"initial_model_{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=1.0 / MIN_PERIOD,
        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_{period_band_name}"
)
[2024-03-15 09:41:25,543] 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_{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_{period_band_name}",
)
[2024-03-15 09:41:26,309] INFO: Resuming iteration #0.
Current stage: initialize
[2024-03-15 09:41:26,309] INFO: 1 new tasks have been issued.
[2024-03-15 09:41:26,309] INFO: Processing task `misfit_and_gradient`
[2024-03-15 09:41:27,319] INFO: Submitting job array with 2 jobs ...

Uploading 1 files...
Uploading 1 files...

[2024-03-15 09:41:27,371] INFO: Launched adjoint simulations for 2 events. Please check again to see if they are finished.
[2024-03-15 09:41:27,371] INFO: Some tasks of iteration #0 are still running. Please check again later.

p.inversions.resume(
    inverse_problem_configuration=f"inv_{period_band_name}",
)
[2024-03-15 09:41:27,378] INFO: Resuming iteration #0.
Current stage: check_convergence
[2024-03-15 09:41:27,378] INFO: Processing task `misfit_and_gradient`
[2024-03-15 09:41:27,503] INFO: 2 events have already been submitted. They will not be submitted again.
[2024-03-15 09:41:27,549] INFO: Some simulations are still running. Please check again to see if they are finished.
[2024-03-15 09:41:27,550] INFO: Some tasks of iteration #0 are still running. Please check again later.
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_{period_band_name}",
)
[2024-03-15 09:41:27,556] INFO: Resuming iteration #0.
Current stage: check_convergence
[2024-03-15 09:41:27,557] INFO: Processing task `misfit_and_gradient`
[2024-03-15 09:41:27,685] INFO: 2 events have already been submitted. They will not be submitted again.
[2024-03-15 09:41:27,723] INFO: Some simulations are still running. Please check again to see if they are finished.
[2024-03-15 09:41:27,724] INFO: Some tasks of iteration #0 are still running. Please check again later.
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_{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_{period_band_name}",
    timeout_in_seconds=1800,
    ping_interval_in_seconds=100,
)
[2024-03-15 09:41:28,187] INFO: Resuming iteration #0.

[2024-03-15 09:41:28,188] INFO: Processing task `misfit_and_gradient`
[2024-03-15 09:41:28,318] INFO: 2 events have already been submitted. They will not be submitted again.
[2024-03-15 09:41:28,354] INFO: Some simulations are still running. Please check again to see if they are finished.
[2024-03-15 09:41:28,355] INFO: Some tasks of iteration #0 are still running. Please check again later.
[2024-03-15 09:43:08,357] INFO: Processing task `misfit_and_gradient`
[2024-03-15 09:43:08,507] INFO: 2 events have already been submitted. They will not be submitted again.
[2024-03-15 09:43:08,562] INFO: Some simulations are still running. Please check again to see if they are finished.
[2024-03-15 09:43:08,564] INFO: Some tasks of iteration #0 are still running. Please check again later.
[2024-03-15 09:44:48,567] INFO: Processing task `misfit_and_gradient`
[2024-03-15 09:44:48,710] INFO: 2 events have already been submitted. They will not be submitted again.
[2024-03-15 09:44:48,770] INFO: Some simulations are still running. Please check again to see if they are finished.
[2024-03-15 09:44:48,771] INFO: Some tasks of iteration #0 are still running. Please check again later.
[2024-03-15 09:46:28,773] INFO: Processing task `misfit_and_gradient`
[2024-03-15 09:46:28,902] INFO: 2 events have already been submitted. They will not be submitted again.
[2024-03-15 09:46:29,394] INFO: 
Iteration 0: Number of events: 2	 chi = 308.4717891328537	 ||g|| = 1.222967280451441e-05
pred = ---	ared = ---	norm_update = ---	tr_radius = ---
[2024-03-15 09:46:29,395] INFO: 1 new tasks have been issued.
[2024-03-15 09:46:29,396] INFO: Processing task `preconditioner`

[2024-03-15 09:46:29,667] INFO: Some tasks of iteration #0 are still running. Please check again later.
[2024-03-15 09:48:09,698] INFO: Processing task `preconditioner`
[2024-03-15 09:48:09,983] INFO: 1 new tasks have been issued.
[2024-03-15 09:48:09,984] INFO: Processing task `misfit`
[2024-03-15 09:48:10,166] INFO: Submitting job array with 2 jobs ...

[2024-03-15 09:48:10,241] INFO: Launched simulations for 2 events. Please check again to see if they are finished.
[2024-03-15 09:48:10,241] INFO: Some tasks of iteration #0 are still running. Please check again later.
[2024-03-15 09:49:50,342] INFO: Processing task `misfit`
[2024-03-15 09:49:51,087] INFO: 
old misfit control group: 308.4717891328537
new misfit control group: 220.000211222491
predicted reduction control group: -44.11213423327639
actual reduction control group: -88.4715779103627
2 out of 2 event(s) improved the misfit.
[2024-03-15 09:49:51,088] INFO: 
Model update accepted.
[2024-03-15 09:49:51,088] INFO: 1 new tasks have been issued.
[2024-03-15 09:49:51,088] INFO: Processing task `finalize_iteration`
[2024-03-15 09:49:51,258] INFO: Succesfully completed iteration #0.
[2024-03-15 09:49:51,259] 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_{period_band_name}",
)
Before continuing with more iterations, let's have a more detailed look at intermediate data of the iteration.
p.viz.nb.gradients(
    simulation_configuration=f"initial_model_{period_band_name}",
    misfit_configuration=f"tf_phase_misfit_{period_band_name}",
    wavefield_compression=sn.WavefieldCompression(
        forward_wavefield_sampling_interval=10
    ),
    events=SELECTED_EVENTS[:2],
)