Download Jupyter Notebook

This page is automatically generated from a Jupyter Notebook. You can view and download it here and here.

Note that everything here is fully automatically tested against - and thus guaranteed to work - only for the latest repository versions of all the Salvus packages. So please update if anything works differently on your machine.

%matplotlib inline

# Import modules and setup paths.
import os
import h5py
import matplotlib.pyplot as plt
import pyasdf

# Paths are grabbed from environment variables.
PARAVIEW_BIN = os.environ["PARAVIEW_BIN"]
SALVUS_BIN = os.environ["SALVUS_BIN"]

# This is currently built together with Salvus.
MISFIT_BIN = os.path.join(os.path.dirname(SALVUS_BIN), 'compute_misfit')
MESH2D = 'Quad_IsotropicElastic2D_Nodal_44x44.e'

In this example we will use time reversal to simulate adjoint equations and to compute sensitivity kernels.

Sensitivity kernels

The sensitivity kernel, i.e., the first derivative of the misfit functional with respect to the material parameters, indicate in which parts of the domain the misfit is sensitive to changes of the material parameters. In other words, the kernel tells us where we need to increase or decrease the model values to reduce the misfit.

At this stage we don't worry about designing a mesh from scratch, and will use a mesh that we have already used in the previous example together with the same source-receiver pair. First, we need to recreate the toml file to specify the source and receiver.

with open("salvus.toml", "w") as text_file:
    recLocX = 500
    recLoxY = 1000
    print("\n[[receiver]]", file=text_file)
    print("network = \"AB\"", file=text_file)
    print("station = \"CED\"", file=text_file)
    print("medium = \"solid\"", file=text_file)
    print("salvus_coordinates = [{0}, {1}]".format(recLocX, recLoxY), file=text_file)

    name = "source1"
    srcLocX = -500.0
    srcLocY = 1000.0
    spatial_type = "vector"
    temporal_type = "ricker"
    scaleX = 1e9
    scaleY = 1e9
    center_frequency = 14.5
    print("\n[[source]]", file=text_file)
    print("name = \"{0}\"".format(name), file=text_file)
    print("location = [{0}, {1}]".format(srcLocX, srcLocY), file=text_file)
    print("spatial_type = \"{0}\"".format(spatial_type), file=text_file)
    print("temporal_type = \"{0}\"".format(temporal_type), file=text_file)
    print("center_frequency = {0}".format(center_frequency), file=text_file)
    print("scale = [{0}, {1}]".format(scaleX,scaleY), file=text_file)

Adjoint methods are a powerful tool to compute sensitivities with respect to all kinds of parameters in PDE-constrained optimization problems.

In this example we consider the isotropic elastic wave equation

\[ \rho\frac{\partial^2}{\partial t^2} u - \nabla \cdot \left(2\cdot \varepsilon(u) + \lambda(\nabla \cdot u) I\right) = f. \]

So we have three parameters here, \(\rho, \lambda\) and \(\mu\), and the kernels are given by the following formulas:

\[ \begin{equation}\begin{aligned}K_{\rho} (\mathbf{x}) &= \int_0^T \frac{\partial^2}{\partial t^2} u(\mathbf{x},t) \cdot u^\dagger(\mathbf{x},t) \,dt,\\ K_{\lambda} (\mathbf{x}) &= \int_0^T (\nabla \cdot u(\mathbf{x},t)) \cdot(\nabla \cdot u^\dagger(\mathbf{x},t)) \,dt,\\ K_{\mu} (\mathbf{x}) &= \int_0^T \varepsilon(u)(\mathbf{x},t) : \varepsilon(u^\dagger)(\mathbf{x},t) \,dt.\end{aligned}\end{equation} \]

Let's recall the steps we need to follow to compute a sensitivity kernel: 1. Run a forward simulation and store the wavefield 2. Compute the adjoint source(s) 3. Run an adjoint simulation and compute the kernel

Which fields do we need to store during the forward simulation?

It can be tedious to keep track of that and Salvus takes care of saving the right fields automatically if you tell it to --save-fields kernel.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 \
 --start-time -0.07 --end-time 0.91 --time-step 0.0007 \
 --source-toml salvus.toml \
 --receiver-toml salvus.toml --receiver-file-name receiver.h5 \
 --receiver-fields u_ELASTIC \
 --save-fields kernel --save-wavefield-file "wavefield.h5" \
 --absorbing-boundaries x0,x1,y0,y1
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Sub-optimal start time detected. This may cause high-frequency artefacts in the simulation. Suggested start time would is -0.107544.
Timestep set by user as: 0.0007 s.
Automated timestep would be set to: 0.000735822 s.
Sub-optimal start time detected. This may cause high-frequency artefacts in the simulation. Suggested start time would is -0.107544.
Begin time loop.
Time loop progress [100%].
Time loop completed in 16.4245 seconds.
Begin post processing.
Problem complete.

Note that we need to set --start-time, --end-time and the --time-step explicitly to make sure we can use the same discretization for the adjoint run again. Check the contents of the hdf5 file to find out which fields Salvus saved.

!h5ls -vr wavefield.h5
Opened "wavefield.h5" with sec2 driver.
/                        Group
    Location:  1:96
    Links:     1
/ELASTIC                 Group
    Location:  1:800
    Links:     1
/ELASTIC/coordinates     Dataset {1936/1936, 25/25, 2/2}
    Location:  1:1832
    Links:     1
    Storage:   774400 logical bytes, 774400 allocated bytes, 100.00% utilization
    Type:      native double
/ELASTIC/data            Dataset {1401/1401, 1936/1936, 5/5, 25/25}
    Attribute: DIMENSION_LABELS {4}
        Type:      variable-length null-terminated ASCII string
        Data:  "time", "element",
               "[ a_ELASTIC_CMP_0 | a_ELASTIC_CMP_1 | strain_11 | strain_22 | strain_12 ]",
               "point"
    Location:  1:776872
    Links:     1
    Storage:   2712336000 logical bytes, 2712336000 allocated bytes, 100.00% utilization
    Type:      native double
/ELASTIC/time            Dataset {1401/1401}
    Location:  1:777216
    Links:     1
    Storage:   11208 logical bytes, 11208 allocated bytes, 100.00% utilization
    Type:      native double

Next, we need to compute the adjoint sources. In this simple example we do not consider observed data to compare with our synthetics, but just try to minimize the following misfit functional

\[ \frac{1}{2}\int_0^T u(\mathbf{x}_{r},t)^2 \,dt, \]

which is some form of ''energy norm''. The sensitivity with respect to the material parameters will show which parts of the domain influence the amount energy that is propagated to the receiver location. The adjoint source is the negative Frechet derivative of the misfit functional with respect to the displacement, which - in this simple case - is just the displacement with a flipped sign, i.e., \(-u(\mathbf{x}_{r},t)\).

Execute the following line to compute the misfit and the adjoint source for this example:

!$MISFIT_BIN salvus.toml receiver.h5 adjoint
Misfit = 1.87591e-07
Adjoint source written to 'adjointSource.h5' and 'adjointSource.toml'

Let's check the contents of the two files that have been created.

! h5ls -vr adjointSource.h5
Opened "adjointSource.h5" with sec2 driver.
/                        Group
    Location:  1:96
    Links:     1
/AuxiliaryData           Group
    Location:  1:800
    Links:     1
/AuxiliaryData/AdjointSources Group
    Location:  1:1832
    Links:     1
/AuxiliaryData/AdjointSources/AB.CED Dataset {1401/1401, 2/2}
    Attribute: dt scalar
        Type:      native double
        Data:  0.0007
    Attribute: location {2}
        Type:      native double
        Data:  500, 1000
    Attribute: spatial-type scalar
        Type:      7-byte null-terminated ASCII string
        Data:  "vector"
    Attribute: starttime scalar
        Type:      native double
        Data:  -0.07
    Location:  1:2864
    Links:     1
    Storage:   22416 logical bytes, 22416 allocated bytes, 100.00% utilization
    Type:      native double

There is a hdf5 file that contains the adjoint sources and a toml file that we can use as input for the adjoint simulation. The toml file specifies the source_input_file where the adjoint traces are stored as well as the location (groups and datasets) in this file. Location, spatial type and time are stored directly in the hdf5 file.

!cat adjointSource.toml
source_input_file = "adjointSource.h5"

[[source]]
name = "AB.CED"
dataset_name = "/AuxiliaryData/AdjointSources/AB.CED"

Next, we should take a look at the adjoint sources. Plot them with pyasdf, but remember that they are not stored as Waveforms, but in the group AuxiliaryData.

with pyasdf.ASDFDataSet("adjointSource.h5") as ds_adjSrc, \
        pyasdf.ASDFDataSet("receiver.h5") as ds_rec:
    f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, sharey=False, figsize=(20,10))
    ax1.plot(ds_adjSrc.auxiliary_data.AdjointSources["AB.CED"].data[:, 0], color="red")
    ax1.set_title('adjoint sources, component 1')
    ax2.plot(ds_adjSrc.auxiliary_data.AdjointSources["AB.CED"].data[:, 1], color="green")
    ax2.set_title('adjoint sources, component 2')
    ax3.plot(ds_rec.waveforms.AB_CED.displacement.select(component="Z")[0].data, color="red")
    ax3.set_title('receiver, component 1')
    ax4.plot(ds_rec.waveforms.AB_CED.displacement.select(component="E")[0].data, color="green")
    ax4.set_title('receiver, component 2')

png

Did you notice that the amplitudes are different? Do you know why this is the case?

After discretizing the problem, the discrete misfit looks as follows:

\[ \frac{1}{2}\int_0^T u(\mathbf{x}_{r},t)^2 \,dt \approx \frac{\Delta t}{2}\sum_{i=0}^{N_t} u(\mathbf{x}_{r},t_i)^2, \]

so the partial derivative with respect to \(u(\mathbf{x}_{r},t_i)\) is \(\Delta t \cdot u(\mathbf{x}_{r},t_i)\).

Ok, almost there. Now we only need to run the adjoint simulation. First we need to tell Salvus, that we want to run a simulation in adjoint mode by setting --adjoint true. Next, we specify the adjoint source with the toml file we created. And finally, we need to tell Salvus to load the forward wavefield and the name of the file to which it should write the kernel using the following options:

--load-wavefield-file wavefield.h5 --load-fields kernel --kernel-file kernel.e

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 \
 --start-time -0.07 --end-time 0.91 --time-step 0.0007 \
 --source-toml adjointSource.toml \
 --adjoint true \
 --load-wavefield-file wavefield.h5 --load-fields kernel --kernel-file kernel.e \
 --absorbing-boundaries x0,x1,y0,y1
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Timestep set by user as: 0.0007 s.
Automated timestep would be set to: 0.000735822 s.
Begin time loop.
Time loop progress [100%].
Time loop completed in 9.48623 seconds.
Begin post processing.
Problem complete.

Are you excited how the kernel will look like?

!$PARAVIEW_BIN kernel.e

The kernel lives in the same vector space as the model so it is stored in the same format. By default, Salvus computes the derivatives with respect to the parameters that enter linearly into the wave equation. But maybe you are interested in the \(v_p\)- and \(v_s\)-kernel instead? OK, we can do that to. Just use the option --kernel-fields RHO,VP,VS to compute the kernels using this parameterization.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 \
 --start-time -0.07 --end-time 0.91 --time-step 0.0007 \
 --source-toml adjointSource.toml \
 --adjoint true \
 --load-wavefield-file wavefield.h5 --load-fields kernel --kernel-file kernel.e \
 --kernel-fields RHO,VP,VS \
 --absorbing-boundaries x0,x1,y0,y1
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Timestep set by user as: 0.0007 s.
Automated timestep would be set to: 0.000735822 s.
Begin time loop.
Time loop progress [100%].
Time loop completed in 10.1612 seconds.
Begin post processing.
Problem complete.

Now take a look at the new kernel file.

!$PARAVIEW_BIN kernel.e

Are you wondering how the kernel evolves over time? You can visualize it easily if you tell Salvus to save the wavefield during the adjoint simulation. Run the adjoint simulation again and store the displacement and the time-dependent kernel which can be accessed with instaKernel. Just make sure to write to a different file, wavefield.h5 is already taken.

# XXX: Currently some issue with the instaKernel option.
# !$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
# --polynomial-order 4 --dimension 2 \
# --start-time -0.07 --end-time 0.91 --time-step 0.0007 \
# --source-toml adjointSource.toml \
# --adjoint true \
# --load-wavefield-file wavefield.h5 --load-fields kernel --kernel-file kernel.e \
# --save-wavefield-file wavefield_adj.h5 --save-fields u_ELASTIC,instaKernel \
# --kernel-fields RHO,VP,VS \
# --absorbing-boundaries x0,x1,y0,y1

Open the xdmf file to plot how the adjoint wavefield and the gradient evolves over time.

# !$PARAVIEW_BIN wavefield_adj_ELASTIC.xdmf

Now redo the example but this time with a free surface (i.e., no absorbing condition at y1).

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 \
 --start-time -0.07 --end-time 0.91 --time-step 0.0007 \
 --source-toml salvus.toml \
 --receiver-toml salvus.toml --receiver-file-name receiver.h5 --receiver-fields u_ELASTIC \
 --save-fields kernel --save-wavefield-file "wavefield.h5" \
 --absorbing-boundaries x0,x1,y0

!$MISFIT_BIN salvus.toml receiver.h5 adjoint

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 \
 --start-time -0.07 --end-time 0.91 --time-step 0.0007 \
 --source-toml adjointSource.toml \
 --adjoint true \
 --load-wavefield-file wavefield.h5 --load-fields kernel --kernel-file kernel.e \
 --save-wavefield-file wavefield_adj.h5 --save-fields u_ELASTIC \
 --absorbing-boundaries x0,x1,y0
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Sub-optimal start time detected. This may cause high-frequency artefacts in the simulation. Suggested start time would is -0.107544.
Timestep set by user as: 0.0007 s.
Automated timestep would be set to: 0.000735822 s.
Sub-optimal start time detected. This may cause high-frequency artefacts in the simulation. Suggested start time would is -0.107544.
Begin time loop.
Time loop progress [100%].
Time loop completed in 15.6224 seconds.
Begin post processing.
Problem complete.

Misfit = 2.00186e-07
Adjoint source written to 'adjointSource.h5' and 'adjointSource.toml'

====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Timestep set by user as: 0.0007 s.
Automated timestep would be set to: 0.000735822 s.
Begin time loop.
Time loop progress [100%].
Time loop completed in 11.4266 seconds.
Begin post processing.
Problem complete.
!$PARAVIEW_BIN kernel.e

Not only the kernel looks different but also the misfit increased, because due to the reflecting surface, more energy is now propagated to the receiver. Have a look at the seismograms and adjoint sources.

with pyasdf.ASDFDataSet("adjointSource.h5") as ds_adjSrc, \
        pyasdf.ASDFDataSet("receiver.h5") as ds_rec:
    f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, sharey=False, figsize=(20,10))
    ax1.plot(ds_adjSrc.auxiliary_data.AdjointSources["AB.CED"].data[:, 0], color="red")
    ax1.set_title('adjoint sources, component 1')
    ax2.plot(ds_adjSrc.auxiliary_data.AdjointSources["AB.CED"].data[:, 1], color="green")
    ax2.set_title('adjoint sources, component 2')
    ax3.plot(ds_rec.waveforms.AB_CED.displacement.select(component="Z")[0].data, color="red")
    ax3.set_title('receiver, component 1')
    ax4.plot(ds_rec.waveforms.AB_CED.displacement.select(component="E")[0].data, color="green")
    ax4.set_title('receiver, component 2')

png

You can see the arrival of the second wave, which explains why the kernel looks different. Now try to restrict the misfit only to the first arrival. You can call the misfit tool with the following arguments

<receiver toml> <receiver hdf5> <misfit start time> <misfit end time> adjoint

!$MISFIT_BIN salvus.toml receiver.h5 0.0 0.4 adjoint
Misfit = 3.02883e-08
Adjoint source written to 'adjointSource.h5' and 'adjointSource.toml'

The misfit is significantly smaller and the adjoint sources changed, too.

with pyasdf.ASDFDataSet("adjointSource.h5") as ds_adjSrc, \
        pyasdf.ASDFDataSet("receiver.h5") as ds_rec:
    f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, sharey=False, figsize=(20,10))
    ax1.plot(ds_adjSrc.auxiliary_data.AdjointSources["AB.CED"].data[:, 0], color="red")
    ax1.set_title('adjoint sources, component 1')
    ax2.plot(ds_adjSrc.auxiliary_data.AdjointSources["AB.CED"].data[:, 1], color="green")
    ax2.set_title('adjoint sources, component 2')
    ax3.plot(ds_rec.waveforms.AB_CED.displacement.select(component="Z")[0].data, color="red")
    ax3.set_title('receiver, component 1')
    ax4.plot(ds_rec.waveforms.AB_CED.displacement.select(component="E")[0].data, color="green")
    ax4.set_title('receiver, component 2')

png

Next, we need to run the adjoint simulation again and check the new kernel.

!$SALVUS_BIN --mesh-file $MESH2D --model-file $MESH2D \
 --polynomial-order 4 --dimension 2 \
 --start-time -0.07 --end-time 0.91 --time-step 0.0007 \
 --source-toml adjointSource.toml \
 --adjoint true \
 --load-wavefield-file wavefield.h5 --load-fields kernel --kernel-file kernel.e \
 --save-wavefield-file wavefield_adj.h5 --save-fields u_ELASTIC \
 --kernel-fields RHO,VP,VS \
 --absorbing-boundaries x0,x1,y0
====================================================
Salvus version 0.0.1-686-g644f5f44-dirty
Floating point size: 64
Compiled for GLL orders: 4
Running on 1 rank(s).
====================================================

Initializing problem.
Timestep set by user as: 0.0007 s.
Automated timestep would be set to: 0.000735822 s.
Begin time loop.
Time loop progress [100%].
Time loop completed in 10.5535 seconds.
Begin post processing.
Problem complete.
!$PARAVIEW_BIN kernel.e

Ente gut, alles gut.