Mapping
class in SalvusOpt. This gives a lot of flexibility, for instance, to use different discretizations (e.g., a coarser mesh for the inversion or an event-dependent mesh for the simulation) or model parameterizations (e.g., inverting only for a subset of the physical parameters).%matplotlib inline
%config Completer.use_jedi = False
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
import matplotlib.pyplot as plt
import numpy as np
import pathlib
import time
import xarray as xr
import salvus.namespace as sn
nx, ny = 10, 10
x = np.linspace(0.0, 3000.0, nx)
y = np.linspace(-1000.0, 0.0, nx)
xx, yy = np.meshgrid(x, y, indexing="ij")
vp = 1500.0 - yy
rho = 1000.0 - yy
ds = xr.Dataset(
data_vars={
"vp": (["x", "y"], vp),
"rho": (["x", "y"], rho),
},
coords={"x": x, "y": y},
)
ds.vp.T.plot()
<matplotlib.collections.QuadMesh at 0x7fb379202610>
p = sn.Project.from_volume_model(
path="project",
volume_model=sn.model.volume.cartesian.GenericModel(name="model", data=ds),
load_if_exists=True,
)
src = sn.simple_config.source.cartesian.ScalarPoint2D(x=500.0, y=-500.0, f=1.0)
rec = sn.simple_config.receiver.cartesian.collections.ArrayPoint2D(
x=np.linspace(100.0, 2900.0, 10), y=0.0, fields=["phi"]
)
p += sn.Event(event_name="event", sources=src, receivers=rec)
p.viz.nb.domain()
ec = sn.EventConfiguration(
waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
end_time_in_seconds=2.0
),
wavelet=sn.simple_config.stf.Ricker(center_frequency=5.0),
)
p += sn.SimulationConfiguration(
name="sim_model",
elements_per_wavelength=2,
tensor_order=4,
max_frequency_in_hertz=10.0,
model_configuration=sn.ModelConfiguration(
background_model=None, volume_models="model"
),
event_configuration=ec,
absorbing_boundaries=sn.AbsorbingBoundaryParameters(
reference_velocity=2000.0,
number_of_wavelengths=3.5,
reference_frequency=5.0,
),
)
p += sn.MisfitConfiguration(
name="misfit",
observed_data=None,
misfit_function="L2_energy_no_observed_data",
receiver_field="phi",
)
gradients = {}
while not gradients:
gradients = p.actions.inversion.compute_gradients(
simulation_configuration="sim_model",
misfit_configuration="misfit",
wavefield_compression=sn.WavefieldCompression(
forward_wavefield_sampling_interval=1
),
events=p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=4,
)
time.sleep(2.0)
raw_gradient = sn.UnstructuredMesh.from_h5(gradients["event"])
[2024-07-02 14:19:48,491] INFO: Creating mesh. Hang on. [2024-07-02 14:19:48,789] INFO: Submitting job ... [2024-07-02 14:19:49,060] INFO: Launched simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:19:51,288] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2407021419292789_60c3caa599@local [2024-07-02 14:19:51,417] INFO: Launched adjoint simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:19:53,518] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:19:53,625] INFO: Some simulations are still running. Please check again to see if they are finished. [2024-07-02 14:19:55,706] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:19:55,755] INFO: Some simulations are still running. Please check again to see if they are finished. [2024-07-02 14:19:57,809] INFO: 1 events have already been submitted. They will not be submitted again.
raw_gradient
. In the widget below, we notice that the sensitivity at the source location and - to a smaller degree - at the receiver locations has significantly higher amplitudes than everywhere else in the domain. This is the result of all energy passing through these points, which is why the waveforms are clearly most sensitive to changes at those locations. However, this clearly does not look like a reasonable model update.raw_gradient
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb36f717a90>
VP
, but could do the same for RHO
.def visualize_gradient(grad, clip=None):
g = grad.model.copy()
mask = np.logical_and(
g.get_element_centroid()[:, 1] > -1000.0,
np.abs(g.get_element_centroid()[:, 0] - 1500.0) < 1500.0,
)
g = g.apply_element_mask(mask)
if clip:
scale = (
clip
* np.max(np.abs(g.elemental_fields["VP"]))
* np.ones_like(g.elemental_fields["VP"])
)
g.elemental_fields["VP"] = np.minimum(g.elemental_fields["VP"], scale)
g.elemental_fields["VP"] = np.maximum(g.elemental_fields["VP"], -scale)
display(g)
prior = p.simulations.get_mesh("sim_model")
absolute
scaling of the physical parameters, there is no difference between both discretizations. Hence the mapped gradient is the same as the raw gradient.map1 = sn.Mapping(
inversion_parameters=["VP"],
scaling="absolute",
)
grad1 = map1.adjoint(
mesh=raw_gradient.copy(),
prior=prior,
)
visualize_gradient(grad1, clip=None)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a25ac310>
VP
in locations different from the source / receivers.visualize_gradient(grad1, clip=0.1)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a25dbe50>
map2 = sn.Mapping(
inversion_parameters=["VP"],
scaling="absolute",
source_cutout_radius_in_meters=200.0,
)
grad2 = map2.adjoint(
mesh=raw_gradient.copy(),
prior=prior,
event=p.waveforms.get(data_name="sim_model", events="event")[0],
)
visualize_gradient(grad2, clip=0.1)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a2655190>
map3 = sn.Mapping(
inversion_parameters=["VP"],
scaling="absolute",
source_cutout_radius_in_meters=200.0,
receiver_cutout_radius_in_meters=100.0,
)
grad3 = map3.adjoint(
mesh=raw_gradient.copy(),
prior=prior,
event=p.waveforms.get(data_name="sim_model", events="event")[0],
)
visualize_gradient(grad3, clip=0.1)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a24bf050>
visualize_gradient(grad3, clip=0.8)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a2145110>
mesh = p.simulations.get_mesh(simulation_configuration="sim_model")
roi = np.zeros_like(mesh.connectivity)
mask = mesh.points[mesh.connectivity][:, :, 1] < -100.0
roi[mask] = 1.0
mesh.attach_field("region_of_interest", roi)
map4 = sn.Mapping(
inversion_parameters=["VP"],
scaling="absolute",
source_cutout_radius_in_meters=200.0,
region_of_interest=roi,
)
grad4 = map4.adjoint(
mesh=raw_gradient.copy(),
prior=prior,
event=p.waveforms.get(data_name="sim_model", events="event")[0],
)
visualize_gradient(grad4, clip=0.8)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a24bc1d0>
map5 = sn.Mapping(
inversion_parameters=["VP"],
scaling="relative_deviation_from_prior",
source_cutout_radius_in_meters=200.0,
region_of_interest=roi,
)
grad5 = map5.adjoint(
mesh=raw_gradient.copy(),
prior=prior,
event=p.waveforms.get(data_name="sim_model", events="event")[0],
)
visualize_gradient(grad5, clip=0.8)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a2527d50>
from salvus.opt.models import UnstructuredModel
smoothed_gradient = UnstructuredModel(
name="smoothed_gradient",
model=p.actions.inversion.smooth_model(
model=grad5.model,
smoothing_configuration=sn.ConstantSmoothing(
smoothing_lengths_in_meters={"VP": [100.0, 50.0]}
),
site_name="local",
ranks_per_job=1,
),
fields=["VP"],
)
visualize_gradient(smoothed_gradient, clip=0.8)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a25f9050>
forward_wavefield_sampling_interval
in WavefieldCompression
is thus an important tuning parameter. Depending on the application and meshing strategy a factor between 5 and 100 is typically achieved.gradients = {}
for i in [1, 5, 10, 20]:
gradient_files = {}
while not gradient_files:
gradient_files = p.actions.inversion.compute_gradients(
simulation_configuration="sim_model",
misfit_configuration="misfit",
wavefield_compression=sn.WavefieldCompression(
forward_wavefield_sampling_interval=i
),
events=p.events.list(),
site_name=SALVUS_FLOW_SITE_NAME,
ranks_per_job=4,
)
time.sleep(2.0)
gradients[i] = sn.UnstructuredMesh.from_h5(gradient_files["event"])
[2024-07-02 14:20:06,689] 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'] [2024-07-02 14:20:06,724] INFO: Submitting job ... [2024-07-02 14:20:06,829] INFO: Launched simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:20:11,174] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2407021420177930_056f2a0cd0@local [2024-07-02 14:20:11,313] INFO: Launched adjoint simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:20:13,383] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:20:13,424] INFO: Some simulations are still running. Please check again to see if they are finished. [2024-07-02 14:20:15,500] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:20:15,658] INFO: Some simulations are still running. Please check again to see if they are finished. [2024-07-02 14:20:17,800] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:20:17,925] INFO: Some simulations are still running. Please check again to see if they are finished. [2024-07-02 14:20:20,011] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:20:22,250] 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'] [2024-07-02 14:20:22,308] INFO: Submitting job ... [2024-07-02 14:20:22,407] INFO: Launched simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:20:24,555] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2407021420567643_e5076f8846@local [2024-07-02 14:20:24,824] INFO: Launched adjoint simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:20:26,887] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:20:26,924] INFO: Some simulations are still running. Please check again to see if they are finished. [2024-07-02 14:20:29,003] INFO: 1 events have already been submitted. They will not be submitted again. [2024-07-02 14:20:31,183] 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'] [2024-07-02 14:20:31,217] INFO: Submitting job ... [2024-07-02 14:20:31,305] INFO: Launched simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:20:33,529] INFO: Submitting job ... Uploading 1 files... 🚀 Submitted job_2407021420532057_ace6517413@local [2024-07-02 14:20:33,601] INFO: Launched adjoint simulations for 1 events. Please check again to see if they are finished. [2024-07-02 14:20:35,635] INFO: 1 events have already been submitted. They will not be submitted again.
for i in [1, 5, 10, 20]:
grad = map5.adjoint(
mesh=gradients[i].copy(),
prior=prior,
event=p.waveforms.get(data_name="sim_model", events="event")[0],
)
print(f"----------------------------------------------------------")
print(f"Mapped gradient with sampling interval {i}:")
visualize_gradient(grad, clip=0.8)
---------------------------------------------------------- Mapped gradient with sampling interval 1:
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb363014b90>
---------------------------------------------------------- Mapped gradient with sampling interval 5:
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb363025690>
---------------------------------------------------------- Mapped gradient with sampling interval 10:
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb3a2298fd0>
---------------------------------------------------------- Mapped gradient with sampling interval 20:
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7fb363193ed0>