MondaicThis micro-tutorial demonstrates how to set up a first full waveform inversion with Salvus.
The simulations are very light and do not require any serious resources.
Some familiary with Python as well as knowledge of how to set up forward simulations with Salvus is assumed.
For an in-depth tutorial of running full waveform inversions with Salvus, please have a look at one of our other tutorials.
import os
import numpy as np
import xarray as xr
import salvus.namespace as sn
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")We'll generate a ring of transducers, five of which act as sources. This mimics an acquisition setup common in medical devices.
# Generate a ring with 5 scalar sources.
sources = sn.simple_config.source.cartesian.collections.ScalarPoint2DRing(
x=0, y=0, radius=0.09, count=5, f=1.0
)
events = []
# Use the same receiver ring for each but mask out the receivers close to each
# source.
for i, s in enumerate(sources._sources):
all_receivers = (
sn.simple_config.receiver.cartesian.collections.RingPoint2D(
x=0, y=0, radius=0.09, count=100, fields=["phi"]
)
)
receivers = [
r
for r in all_receivers._receivers
if np.sqrt(
(s.location[0] - r.location[0]) ** 2
+ (s.location[1] - r.location[1]) ** 2
)
> 0.03
]
events.append(
sn.Event(event_name=f"event_{i}", sources=s, receivers=receivers)
)For FWI we'll need something to invert for. For most applications this will be recorded data from some real world setting. Here we'll define a synthetic inversion target and generate our own data.
def get_spherical_inclusion() -> xr.Dataset:
nx, ny = 200, 200
x = np.linspace(-0.1, +0.1, nx)
y = np.linspace(-0.1, +0.1, ny)
xx, yy = np.meshgrid(x, y, indexing="ij")
# Add 3 spherical inclusions
vp = 1500.0 * np.ones_like(xx)
rho = 980.0 * np.ones_like(xx)
mask = np.sqrt(xx**2 + yy**2) < 0.05
vp[mask] = 1480.0
rho[mask] = 1000.0
mask = np.sqrt(xx**2 + (yy - 0.025) ** 2) < 0.015
vp[mask] = 1550.0
rho[mask] = 1040.0
mask = np.sqrt(xx**2 + (yy + 0.025) ** 2) < 0.015
vp[mask] = 1460.0
rho[mask] = 1010.0
ds = xr.Dataset(
data_vars={
"vp": (["x", "y"], vp),
"rho": (["x", "y"], rho),
},
coords={"x": x, "y": y},
)
return ds
true_model = get_spherical_inclusion()
ax = true_model.vp.T.plot(figsize=(10, 6))
ax.axes.set_aspect("equal")
# Use the volume model to set up the project.
volume_model = sn.model.volume.cartesian.GenericModel(
name="true_model", data=true_model
)
p = sn.Project.from_volume_model(
path="project_full_waveform_inversion",
volume_model=volume_model,
load_if_exists=True,
)
# Add the events.
for event in events:
p.add_to_project(event)
p.visualizations.notebook.domain()
# Generate the data for the target model.
wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00015)
wsc.physics.wave_equation.time_step_in_seconds = 4.0e-7
ec = sn.EventConfiguration(
waveform_simulation_configuration=wsc,
wavelet=sn.simple_config.stf.Ricker(center_frequency=50000.0),
)
p += sn.SimulationConfiguration(
name="true_model_100kHz",
elements_per_wavelength=2,
tensor_order=4,
max_frequency_in_hertz=100000.0,
model_configuration=sn.ModelConfiguration(
background_model=None, volume_models="true_model"
),
event_configuration=ec,
)
p.simulations.launch(
simulation_configuration="true_model_100kHz",
events=p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=1,
)
p.simulations.query(block=True)[2025-01-09 21:39:25,789] INFO: Creating mesh. Hang on. [2025-01-09 21:39:25,883] INFO: Submitting job array with 5 jobs ...
True
Now we do exactly the same thing, but we'll use a homogeneous initial model.
bm = sn.model.background.homogeneous.IsotropicAcoustic(vp=1500.0, rho=980.0)
mc = sn.ModelConfiguration(background_model=bm)
wsc = sn.WaveformSimulationConfiguration(end_time_in_seconds=0.00015)
p += sn.SimulationConfiguration(
name="initial_model",
elements_per_wavelength=2,
tensor_order=2,
max_frequency_in_hertz=100000.0,
model_configuration=mc,
event_configuration=sn.EventConfiguration(
waveform_simulation_configuration=wsc,
wavelet=sn.simple_config.stf.Ricker(center_frequency=50000.0),
),
)
p.simulations.launch(
simulation_configuration="initial_model",
events=p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=1,
)
p.simulations.query(block=True)[2025-01-09 21:39:30,590] INFO: Creating mesh. Hang on. [2025-01-09 21:39:30,633] INFO: Submitting job array with 5 jobs ...
True
Salvus allows for a lot of control and flexibility in how exactly inverse problems are tackled. Things like data processing, window selection, custom weighting schemes, data normalization, bound constraints, wavelength dependent smoothing, and many more are supported but we don't show them in this micro-tutorial. Please refer to some advanced tutorials for more information.
In all cases, an inverse problem requires the definition of a misfit, e.g. the optimization target, as well the actual configuration of the inverse problem.
# The misfit configuration defines how synthetics are compared to observed data.
p += sn.MisfitConfiguration(
name="L2",
# Could be observed data. Here we compare to the synthetic target.
observed_data="true_model_100kHz",
# Salvus comes with a variety of misfit functions. You can
# also define your own.
misfit_function="L2",
# This is an acoustic simulation so we'll use recordings of phi.
receiver_field="phi",
)
# Now we define the actual inverse problem
p += sn.InverseProblemConfiguration(
name="my_inversion",
# Starting model.
prior_model="initial_model",
# The events to use.
events=p.events.list(),
# What parameters to invert for.
mapping=sn.Mapping(scaling="absolute", inversion_parameters=["VP", "RHO"]),
# The smoothing.
preconditioner=sn.ConstantSmoothing({"VP": 0.01, "RHO": 0.01}),
# The inversion method.
method=sn.TrustRegion(initial_trust_region_linf=10.0),
# The misfit configuration we defined above.
misfit_configuration="L2",
# Compress the forward wavefield by subsampling in time.
wavefield_compression=sn.WavefieldCompression(
forward_wavefield_sampling_interval=10
),
# Job submission settings.
job_submission=sn.SiteConfig(
site_name=SALVUS_FLOW_SITE_NAME, ranks_per_job=1
),
)# Lastly we perform two iterations, and have a look at the results.
for i in range(2):
p.inversions.iterate(
inverse_problem_configuration="my_inversion",
timeout_in_seconds=360,
ping_interval_in_seconds=2,
delete_disposable_files="all",
)[2025-01-09 21:39:35,080] INFO: Adding new iteration #0. [2025-01-09 21:39:35,086] INFO: Resuming iteration #0. [2025-01-09 21:39:35,087] INFO: 1 new tasks have been issued. [2025-01-09 21:39:35,087] INFO: Processing task `misfit_and_gradient` [2025-01-09 21:39:35,150] INFO: The following events have been simulated before, but checkpoints are not available for this combination of `site_name` and `ranks_per_job` and wavefield compression settings. They will be run again: ['event_0', 'event_1', 'event_2', 'event_3', 'event_4'] [2025-01-09 21:39:35,170] INFO: Submitting job array with 5 jobs ... [2025-01-09 21:39:35,258] INFO: Launched simulations for 5 events. Please check again to see if they are finished. [2025-01-09 21:39:35,259] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-01-09 21:39:37,260] INFO: Processing task `misfit_and_gradient` [2025-01-09 21:39:37,346] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-01-09 21:39:39,348] INFO: Processing task `misfit_and_gradient` [2025-01-09 21:39:41,070] INFO: Submitting job array with 5 jobs ...
Uploading 1 files... Uploading 1 files... Uploading 1 files... Uploading 1 files... Uploading 1 files... [2025-01-09 21:39:41,122] INFO: Launched adjoint simulations for 5 events. Please check again to see if they are finished. [2025-01-09 21:39:41,123] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-01-09 21:39:43,124] INFO: Processing task `misfit_and_gradient` [2025-01-09 21:39:43,355] INFO: 5 events have already been submitted. They will not be submitted again. [2025-01-09 21:39:43,417] INFO: Some simulations are still running. Please check again to see if they are finished. [2025-01-09 21:39:43,418] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-01-09 21:39:45,419] INFO: Processing task `misfit_and_gradient` [2025-01-09 21:39:45,641] INFO: 5 events have already been submitted. They will not be submitted again. [2025-01-09 21:39:46,228] INFO: Iteration 0: Number of events: 5 chi = 0.01769007842184961 ||g|| = 0.022565013801719613 pred = --- ared = --- norm_update = --- tr_radius = --- [2025-01-09 21:39:46,229] INFO: 1 new tasks have been issued. [2025-01-09 21:39:46,229] INFO: Processing task `preconditioner` [2025-01-09 21:39:46,295] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-01-09 21:39:48,305] INFO: Processing task `preconditioner` [2025-01-09 21:39:48,422] INFO: 1 new tasks have been issued. [2025-01-09 21:39:48,422] INFO: Processing task `misfit` [2025-01-09 21:39:48,464] INFO: Submitting job array with 5 jobs ... [2025-01-09 21:39:48,551] INFO: Launched simulations for 5 events. Please check again to see if they are finished. [2025-01-09 21:39:48,552] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-01-09 21:39:50,574] INFO: Processing task `misfit` [2025-01-09 21:39:50,621] INFO: Some tasks of iteration #0 are still running. Please check again later. [2025-01-09 21:39:52,643] INFO: Processing task `misfit` [2025-01-09 21:39:54,109] INFO: old misfit control group: 0.01769007842184961 new misfit control group: 0.006215538423932317 predicted reduction control group: -0.009679727771306901 actual reduction control group: -0.011474539997917293 5 out of 5 event(s) improved the misfit. [2025-01-09 21:39:54,110] INFO: Model update accepted. [2025-01-09 21:39:54,110] INFO: 1 new tasks have been issued. [2025-01-09 21:39:54,110] INFO: Processing task `finalize_iteration` [2025-01-09 21:39:54,152] INFO: ... searching for obsolete files in project_full_waveform_inversion/INVERSIONS/my_inversion/00000 [2025-01-09 21:39:54,162] INFO: Freed up 1.2 MB of space. [2025-01-09 21:39:54,162] INFO: Succesfully completed iteration #0. [2025-01-09 21:39:54,164] INFO: Adding new iteration #1. [2025-01-09 21:39:54,168] INFO: Resuming iteration #1. [2025-01-09 21:39:54,169] INFO: 1 new tasks have been issued. [2025-01-09 21:39:54,169] INFO: Processing task `gradient` [2025-01-09 21:39:54,350] INFO: Submitting job array with 5 jobs ...
Uploading 1 files... Uploading 1 files... Uploading 1 files... Uploading 1 files... Uploading 1 files... [2025-01-09 21:39:54,407] INFO: Launched adjoint simulations for 5 events. Please check again to see if they are finished. [2025-01-09 21:39:54,408] INFO: Some tasks of iteration #1 are still running. Please check again later. [2025-01-09 21:39:56,410] INFO: Processing task `gradient` [2025-01-09 21:39:56,535] INFO: 5 events have already been submitted. They will not be submitted again. [2025-01-09 21:39:56,598] INFO: Some simulations are still running. Please check again to see if they are finished. [2025-01-09 21:39:56,600] INFO: Some tasks of iteration #1 are still running. Please check again later. [2025-01-09 21:39:58,601] INFO: Processing task `gradient` [2025-01-09 21:39:58,728] INFO: 5 events have already been submitted. They will not be submitted again. [2025-01-09 21:39:59,119] INFO: Iteration 1: Number of events: 5 chi = 0.006215538423932317 ||g|| = 0.009661533661238208 pred = -0.009679727771306901 ared = -0.011474539997917293 norm_update = 0.9534161427441019 tr_radius = 0.9534163618503646 [2025-01-09 21:39:59,129] INFO: 1 new tasks have been issued. [2025-01-09 21:39:59,129] INFO: Processing task `preconditioner` [2025-01-09 21:39:59,204] INFO: Some tasks of iteration #1 are still running. Please check again later. [2025-01-09 21:40:01,224] INFO: Processing task `preconditioner` [2025-01-09 21:40:01,528] INFO: 1 new tasks have been issued. [2025-01-09 21:40:01,528] INFO: Processing task `misfit` [2025-01-09 21:40:01,567] INFO: Submitting job array with 5 jobs ... [2025-01-09 21:40:01,657] INFO: Launched simulations for 5 events. Please check again to see if they are finished. [2025-01-09 21:40:01,658] INFO: Some tasks of iteration #1 are still running. Please check again later. [2025-01-09 21:40:03,691] INFO: Processing task `misfit` [2025-01-09 21:40:03,738] INFO: Some tasks of iteration #1 are still running. Please check again later. [2025-01-09 21:40:05,773] INFO: Processing task `misfit` [2025-01-09 21:40:07,243] INFO: old misfit control group: 0.006215538423932317 new misfit control group: 0.004515146195660988 predicted reduction control group: -0.0008915795911773164 actual reduction control group: -0.0017003922282713291 5 out of 5 event(s) improved the misfit. [2025-01-09 21:40:07,243] INFO: Model update accepted. [2025-01-09 21:40:07,243] INFO: 1 new tasks have been issued. [2025-01-09 21:40:07,244] INFO: Processing task `finalize_iteration` [2025-01-09 21:40:07,305] INFO: ... searching for obsolete files in project_full_waveform_inversion/INVERSIONS/my_inversion/00001 [2025-01-09 21:40:07,416] INFO: Freed up 4.8 MB of space. [2025-01-09 21:40:07,416] INFO: Succesfully completed iteration #1. [2025-01-09 21:40:07,418] INFO: Adding new iteration #2.
p.viz.nb.inversion(inverse_problem_configuration="my_inversion")
# + {"tags": ["remove"]}
# Will be removed for the website.
p.inversions.set_stopping_criteria(
inverse_problem_configuration="my_inversion",
criteria={"max_iterations_global": 4},
)
# -
PAGE CONTENTS