3. Kernels and Adjoint Simulations#

In this notebook we will build upon the forward simulation training
Participants will learn about adjoint simulations and their role in seismic imaging
We will run adjoint simulations using synthetic data
Objective: Introduce the concepts of: misfit, adjoint sources, and kernels and run an adjoint simulation
Adjoint simulations are key for performing seismic imaging (Day 3) as their results guide iterative model updates during the inverse problem
These instructions should be run from inside a Docker container, using Jupyter Lab (see instructions here).

Relevant Links:

This Notebook: https://github.com/adjtomo/adjdocs/blob/main/workshops/2024-5-21_scoped_uw/3_adjoint_simulations.ipynb

Jupyter Quick Tips:

Run cells one-by-one by hitting the \(\blacktriangleright\) button at the top, or by hitting Shift + Enter
Run all cells by hitting the \(\blacktriangleright\blacktriangleright\) button at the top, or by running Run -> Run All Cells
Currently running cells that are still processing will have a [*] symbol next to them
Finished cells will have a [1] symbol next to them. The number inside the brackets represents what order this cell has been run in.
Commands that start with ! are Bash commands (i.e., commands you would run from the terminal)
Commands that start with % are Jupyter Magic commands.

0) Motivation#

In this Notebook we will be running adjoint simulations
Our aim here is to generate a simple kernel as in the following figure from Tape et al. (2007)
In this notebook, we will manually assemble all the necessary ingredients to generate our kernel

tape_2007_f5

1) Setting Up a SPECFEM2D Working Directory#

As with Day 1, before we start we want to set up a clean working directory to run SPECFEM2D
Reminder that SPECFEM only requires the bin/, DATA/ and OUTPUT_FILES/ directories to run most exectuables
This will help us preserve our cloned repository and reduce file clutter
We will be doing all our work in the directory /home/scoped/work/day_2, all the following cells assume that we are in this directory

# Python packages we will use in this notebook
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image
from scipy.integrate import simpson
from seisflows.tools.model import Model 
from pyatoa import Config, Manager, logger
from pysep.utils.io import read_sem

# Make correct dir. and move there
! mkdir -p /home/scoped/work/kernels
%cd /home/scoped/work/kernels

# Symlink the executables and copy the relevant DATA/ directory
! ln -s /home/scoped/specfem2d/bin .
! cp -r /home/scoped/specfem2d/EXAMPLES/Tape2007/DATA .
! mkdir OUTPUT_FILES

! ls

2) Experimental Setup: Tape et al. 2007 Example Problem#

As with Notebook 2, we will be working with an example problem from the Tape et al. 2007 GJI publication
We will be revisiting the homogeneous halfspace model we saw in Notebook 1.
We will also see a perturbation checkerboard model which defines a checkerboard with \(\pm\)10% velocity perturbations
The perturbation checkerboard is overlain ontop of the homogeneous halfspace model
Here we use Python (NumPy and Matplotlib) to visualize the Example problem setup
In your own research you may use their own preferred tools to visualize their experimental setup

# Incase participants get lost, each numbered section 
# contains the absolute work dir path
%cd /home/scoped/work/kernels

2a) Homogeneous Halfspace Model#

The homogeneous halfspace model example is defined in the Par_file,
It defines a region with P-wave velocity \(V_p=\)5.8km/s and S-wave velocity \(V_s\)=3.5km/s
This is the same model that we looked at in Day 1B. We plot it below for reference

def plot_homogeneous_halfspace():
    """
    Plots a representation of the SPECFEM2D 
    homogeneous halfspace model
    """
    # Sets the X and Z dimensions of our mesh
    x = np.arange(0, 480000, 4000)
    z = np.arange(0, 480000, 4000)
    
    # Reformat the 1D arrays into 2D
    xv, zv = np.meshgrid(x, z)

    # Set a homogeneous value of Vs=3.5km/s 
    vs = 3.5 * np.ones(np.shape(xv))

    # Plot the arrays as a homogeneous halfspace
    plt.tricontourf(xv.flatten(), zv.flatten(), vs.flatten(), 
                    cmap="seismic_r", vmin=3.1, vmax=4.)
    plt.colorbar(label="Vs [km/s]", format="%.1f")
    plt.xlabel("X [m]")
    plt.ylabel("Z [m]")
    plt.title("2D Homogeneous Halfspace Model\n Vs=3.5km/s")
    
# Calls the function we just defined
plot_homogeneous_halfspace()

2b) Perturbation Checkerboard Model#

The checkerboard model features smoothly varying 2D Gaussians which perturbs the homogeneous halfspace model (\(V_s\) and \(V_p\)) by roughly \(\pm\)10%
The checkerboard model is defined by an external velocity model file: DATA/model_velocity.dat_checker
SPECFEM can read external files of various formats to define the velocity model.
SPECFEM3D NOTE: In SPECFEM3D we typically use external tomography files (.xyz) which defines coordinates (X, Y, Z) and material properties (e.g., Vp, Vs, rho, Q)
The following checkerboard model file similarly defines coordinates and material properties

The columns of the checkerboard model are defined:

line_number x[m] z[m] density vp[m/s] vs[m/s]

# Look at the first few lines of the checkerboard model
! head DATA/model_velocity.dat_checker

# Grab values of X, Z, Vs and Vp for plotting
chkbd_x, chkbd_z, chkbd_vp, chkbd_vs = np.genfromtxt("DATA/model_velocity.dat_checker", dtype=float, usecols=[1,2,4,5]).T

def plot_checkerboard(x, z, c, label, alpha=1):
    """Simple re-usable model plotting function"""
    plt.tricontourf(x, z, c, levels=125, cmap="seismic_r", alpha=alpha)
    plt.xlabel("X [m]")
    plt.ylabel("Z [m]")
    plt.title("Checkerboard Model")
    plt.colorbar(label=label)
    
plot_checkerboard(chkbd_x, chkbd_z, chkbd_vp, label="Vp [m/s]")

plot_checkerboard(chkbd_x, chkbd_z, chkbd_vs, label="Vs [m/s]")

Understanding the checkerboard#

The mean (perturbed) velocity is equal to the homogeneous halfspace counterpart
red/warm colors equal relatively slow velocities
blue/cold colors equal relatively fast velocities
This is the typical color convention scheme used in seismic imaging

2c) Visualizing Source-Receiver Geometry#

It’s useful to plot SOURCES and STATIONS with respect to the checkerboard model
These figures can give us some idea of how the underlying model will affect the resulting synthetics
These data grabbing / plotting functions are the same as Day 1

# Grab coordinates from STATIONS file
sta_x, sta_z = np.genfromtxt("DATA/STATIONS_checker", dtype=float, usecols=[2, 3]).T

# Grab coordinates from SOURCE files
ev_x, ev_z = [], []
for i in range(1, 26):
    source_file = f"DATA/SOURCE_{i:0>3}"
    with open(source_file, "r") as f:
        lines = f.readlines()
    # Trying to break apart the following line
    # 'xs = 299367.72      # source location x in meters\n'
    xs = float(lines[2].split("=")[1].split("#")[0].strip())
    zs = float(lines[3].split("=")[1].split("#")[0].strip())
    
    ev_x.append(xs)
    ev_z.append(zs)
    
# Plot SOURCES and STATIONS on top of the checkerboard model
plot_checkerboard(chkbd_x, chkbd_z, chkbd_vp, label="Vp [m/s]")
plt.scatter(sta_x, sta_z, c="c", marker="v", s=40, edgecolor="k")
plt.scatter(ev_x, ev_z, c="y", marker="*", s=150, edgecolor="k")

# Plot SOURCES next to source names
plot_checkerboard(chkbd_x, chkbd_z, chkbd_vp, label="Vp [m/s]", alpha=0.1)

for i, (x, z) in enumerate(zip(ev_x, ev_z)):
    plt.scatter(x, z, c="y", marker="*", edgecolor="k", s=80)
    plt.text(x, z, f"{i+1:0>3}", c="k")  # SOURCE numbering starts at 1
plt.title(f"SOURCES; N={len(ev_x)}")

# Plot STATIONS next to source names
plot_checkerboard(chkbd_x, chkbd_z, chkbd_vp, label="Vp [m/s]", alpha=0.1)

for i, (x, z) in enumerate(zip(sta_x, sta_z)):
    plt.scatter(x, z, c="c", marker="v", s=8, edgecolor="k")
    plt.text(x, z, f"{i:0>3}", fontsize=9)
plt.title(f"STATIONS; N={len(sta_x)}")

# Plot the experimental setup we'll use for Section 2
plot_checkerboard(chkbd_x, chkbd_z, chkbd_vp, label="Vp [m/s]")

# Plot SOURCE 1, STATION 1 and a connecting line
plt.plot((sta_x[0], ev_x[0]), (sta_z[0], ev_z[0]), c="k", ls=":", lw=2, zorder=1)
plt.scatter(sta_x[0], sta_z[0], c="c", marker="v", s=100, edgecolor="k")
plt.scatter(ev_x[0], ev_z[0], c="y", marker="*", edgecolor="k", s=200)

plt.title(f"SOURCE 01; STATION 01")

Understanding Source-Receiver Figures#

Upside-down blue triangles represent 132 receivers/stations
Yellow stars represent 25 sources/events
Red colors in the checkerboard model represent relatively slow velocities
Blue colors in the checkerboard model represent relatively fast velocities

3) Setting up an Adjoint Simulation#

In this section we will learn how to set SPECFEM for an adjoint simulation
We will need to generate adjoint sources (defined below)
We will also need to make some Par_file edits

Useful Definitions#

Adjoint simulation: simulates the interaction of a forward and adjoint wavefield
Forward wavefield: the seismic wavefield propagated from the source location
Adjoint wavefield: a wavefield that propagates from receiver locations, whos time-dependent amplitude is controlled by adjoint sources
Adjoint source: time-reversed waveforms input at receiver locations. Typically they contain information about the forward wavefield (sensitivity kernels), or data-synthetic misfit (misfit kernel)
Kernel: The volumetric integration of the interaction between the forward and adjoint wavefields, highlighting regions/parameters of the model that have an effect on the wavefield or misfit

Useful Information#

SPECFEM saves the forward wavefield and propagates it backwards in time
Adjoint wavefield propagates adjoint sources from each receiver simultaneously (i.e., 1 adjoint simulation per event)
Each component of each receiver therefore is capable of generating an adjoint source
Adjoint sources can be selectively windowed to isolate certain parts of the synthetic seismogram

tape_etal_2007_f5

Tape et al. 2007 Figure 5: Sequence of interactions between the regular and adjoint wavefields during the construction of a traveltime cross‐correlation event kernel K(x). The ⋆ symbol denotes the source, and the Δ symbol denotes the receiver. Each row represents the time‐step indicated on the left… The event kernel shows the region of the current model that gives rise to the discrepancy between the data and the synthetics.

%cd /home/scoped/work/kernels

3a) Generating ‘Data’ with a Target Model#

Objective: Generate synthetic seismogram using checkerboard model
In the context of seismic imaging, adjoint sources are proportional to data-synthetic misfit
We use the checkerboard model shown above as a True or Target
Target model is used to generate synthetic waveforms that are meant to approximate real world data.
In a real seismic inversion, data are actual waveforms recorded during a seismic event
We will choose a single source and a single receiver to generate our misfit kernel

# We use SOURCE 1 as our preferred event
! cp -f DATA/SOURCE_001 DATA/SOURCE

# We take the first line of the STATIONS file as our single station
! head -1 DATA/STATIONS_checker > DATA/STATIONS

# Copy the Parameter file and make some adjustments 
! cp -f DATA/Par_file_Tape2007_132rec_checker DATA/Par_file

# Set parameters 
! seisflows sempar -P DATA/Par_file NSTEP 5000  # to match homogeneous halfspace Par_file
! seisflows sempar -P DATA/Par_file save_model binary
! seisflows sempar -P DATA/Par_file setup_with_binary_database 1

# Ensure that SPECFEM can find the checkerboard model by naming it correctly
! cp -f DATA/model_velocity.dat_checker DATA/proc000000_model_velocity.dat_input

# Run Meshfem and SPECFEM SPECFEM2D 
! mpirun -n 1 bin/xmeshfem2D > OUTPUT_FILES/output_mesher.txt
! mpirun -n 1 bin/xspecfem2D > OUTPUT_FILES/output_solver.txt

# Rename the OUTPUT_FILES so that our next run doesn't overwrite files
! mv OUTPUT_FILES OUTPUT_FILES_TRUE

# Print the log message when the simulation is finished
! tail OUTPUT_FILES_TRUE/output_solver.txt

Parameter changes explained#

NSTEP (number of time steps) is increased to 5000 to set synthetic waveform length
save_model is set to binary so that SPECFEM outputs FORTRAN binary files defining the model
setup_with_binary_database tells SPECFEM to write DATABASE files rather than discarding them

MPI NOTE: We ran this simulation in serial (NPROC==1). If you wanted to this simulation in parallel (NPROC>1) you would also need to manually partition the checkerboard model (model_velocity.dat_checker -> proc000???_model_velocity.dat_input). This is a unique caveat of SPECFEM2D’s Legacy input model type. Those using external tomography models (.xyz files) will not face this requirement.

# SPECFEM has generated one synthetic seismogram
! ls OUTPUT_FILES_TRUE/*semd

# We can look at this 'data' waveform using SeisFlows plot
! seisflows plotst OUTPUT_FILES_TRUE/AA.S000000.BXY.semd --savefig AA.S000000.BXY.png
Image("AA.S000000.BXY.png")

3b) Generating Synthetics using Initial Model#

Objective: generate synthetic waveforms using the homogeneous halfspace model.
In this example we will use the homogeneous halfspace model as our starting or initial model
We use the same source and receivers as the Target model, to ensure that we can compare the waveforms generated
We expect that the waveform generated by the Target and Initial models will differ

# The SOURCE and STATIONS files should remain the same,
# we only want to tell SPECFEM to use the homogeneous halfspace model
! cp -f DATA/Par_file_Tape2007_onerec DATA/Par_file

# Set some Par_file parameters to match the previous run
! seisflows sempar -P DATA/Par_file use_existing_stations .true.
! seisflows sempar -P DATA/Par_file save_model binary
! seisflows sempar -P DATA/Par_file setup_with_binary_database 1
! seisflows sempar -P DATA/Par_file save_forward .true.

# Run SPECFEM with the homogeneous halfspace model, defined in the Par_file
! rm -r OUTPUT_FILES
! mkdir OUTPUT_FILES
! mpirun -n 1 bin/xmeshfem2D > OUTPUT_FILES/output_mesher.txt
! mpirun -n 1 bin/xspecfem2D > OUTPUT_FILES/output_solver.txt

# Copy the OUTPUT_FILES so that our subsequent run doesn't overwrite files
# We do this because the adjoint simulation will require DATABASE files located within
! cp -r OUTPUT_FILES OUTPUT_FILES_INIT

# Print the log message when the simulation is finished
! tail OUTPUT_FILES_INIT/output_solver.txt

Parameter changes explained#

The parameters we have changed are the same as in Subsection (a)
The new save_forward parameter tells SPECFEM to save a snapshot of the last frame of the wavefield
This snapshot is used to back-construct the forward wavefield during an adjoint simulation

# Look at the resulting outputs of the `save_forward` parameter
! diff -qr OUTPUT_FILES_INIT/ OUTPUT_FILES_TRUE/ | grep 'Only'

# Again we can look at the waveform. We will compare in the next section
! ls OUTPUT_FILES_INIT/*semd

! seisflows plotst OUTPUT_FILES_INIT/AA.S000000.BXY.semd --savefig AA.S000000.BXY.png
Image("AA.S000000.BXY.png")

! head OUTPUT_FILES/AA.S000000.BXY.semd

3c) Misfit Quantification and Adjoint Sources#

We use the term misfit to define differences between data and synthetics
Misfit can be defined based on a number of factors such as: phase, amplitude, or waveform differences
Adjoint sources are waveforms proportional to the time-dependent misfit
The chosen misfit function directly influences adjoint simulations and seismic imaging studies by extension

# Read in the two-column ASCII files using NumPy where - t=time; d=data
t_init, d_init = np.loadtxt("OUTPUT_FILES_INIT/AA.S000000.BXY.semd").T
t_true, d_true = np.loadtxt("OUTPUT_FILES_TRUE/AA.S000000.BXY.semd").T

# Plot both waveforms using Matplotlib
plt.plot(t_init, d_init, c="r", label="MODEL INIT; SYNTHETIC")
plt.plot(t_true, d_true, c="k", label="MODEL TRUE; 'DATA'")
plt.xlabel("Time [s]")
plt.ylabel("Displacement [m]")
plt.legend()
plt.show()

# Zoom in on waveform figure ignoring T<50s
plt.plot(t_init, d_init, c="r", label="MODEL INIT; SYNTHETIC")
plt.plot(t_true, d_true, c="k", label="MODEL TRUE; 'DATA'")
plt.xlabel("Time [s]")
plt.ylabel("Displacement [m]")
plt.xlim((50,250))
plt.legend()
plt.show()

Comments on Default Waveform Figure#

INIT (red) synthetic shows a phase delay (negative time shift) with respect to the TRUE synthetic (black)
This reflects that longer relative path length in the red/slow potion of the checkerboard

4) Misfit Quantification#

Objective: Quantify the misfit between data and synthetic and generate adjoint source using Python and showing alternatives with Pyatoa
Misfit function will be defined by waveform misfit (also known as objective function)
Adjoint source defines the gradient of the misfit function
Acknowledging that some of this section is modified from Ridvan Orsvuran’s Jupyter Notebook

4a) Misfit Calculation with Python#

Waveform misfit function: \( \chi = \frac{1}{2} \int [s(t)-d(t)]^2 dt~\)
- d(t): time-dependent ‘data’ waveform (black trace Section 2c)
- s(t): time-dependent ‘synthetic’ waveform (red trace Section 2c)
In Python we will perform the integration step using using Simpson’s rule

# Integrate using scipy
dt = t_true[1] - t_true[0]  # dt represents the time step
misfit = 1/2 * simpson((d_true - d_init)**2, dx=dt)

print(f"misfit={misfit:.3E}")

4b) Adjoint Source with Python#

The adjoint source is a time-dependent measure of misfit
It is defined as the gradient of the misfit function
Waveform adjoint source: \(f^\dagger (t) = s(t) - d(t)\)

# Calculate the adjoint source using data and synthetic
d_adj = d_init - d_true

# Plot both waveforms using Matplotlib
plt.plot(t_init, d_init, c="r", label="MODEL INIT; SYNTHETIC")
plt.plot(t_true, d_true, c="k", label="MODEL TRUE; 'DATA'")
plt.plot(t_init, d_adj, c="g", label="ADJOINT SOURCE")
plt.xlabel("Time [s]")
plt.ylabel("Displacement [m]")
plt.xlim([50, 175])
plt.legend()
plt.show()

Understanding the waveform adjoint source figure#

We can see that the adjoint source is equal to 0 when the waveforms are the same and peaks when the waveforms are the most different
In this case we only have one component of an adjoint source (Y)
The adjoint source will be fed in at receiver locations during an adjoint simulation
The resulting adjoint wavefield interacts with the forward wavefield to illuminate parts of the model that the misfit function is sensitive to

5) Running an Adjoint Simulation#

To run an adjoint simulation in SPECFEM2D, we need to make a few changes to the parameter file
We will use our Python-derived adjoint source, which needs to be placed in a directory locateable by SPECFEM
COMPONENT NOTE: SPECFEM expects adjoint sources for ALL components, even if they do not have adjoint sources. These should just be waveforms with amplitude 0
SPECFEM3D NOTE: SPECFEM3D also requires a DATA/STATIONS_ADJOINT file, which defines stations that have corresponding adjoint sources.
TIME AXIS NOTE: SPECFEM does not expect adjoint sources to be time-reversed. It does this internally

%cd /home/scoped/work/kernels

5a) Writing Adjoint Sources#

IMPORTANT: SPECFEM expects adjoint sources in a directory called SEM/ with a specific format and filename

Adjoint sources should be defined in the same way that the synthetics have been created
In this case that means two-column ASCII files where the time axis exactly matches the the synthetic outputs
Filename must also match synthetics, but with a .adj suffix; i.e., for synthetic AA.S000000.BXY.semd, corresponding adjoint source is: AA.S000000.BXY.adj

# Make the requisite SEM/ directory
! mkdir SEM/

# Generate the two-column (time, data) format required
adjoint_source = np.vstack((t_init, d_adj)).T

# Save the .adj file as an ASCII file
np.savetxt("SEM/AA.S000000.BXY.adj", adjoint_source)
! head SEM/AA.S000000.BXY.adj

5b) Run `xspecfem2D`#

Now that we have the required setup in the Par_file, and the SEM/ directory with the .adj adjoint source files, we can run the solver xspecfem2D
We do not need to re-run the mesher as the adjoint simulations will run on the same numerical mesh that we used for our solver

Parameter file changes#

simulation_type: Simulation type 3 tells SPECFEM that we are running an adjoint simulation
save_ASCII_kernels: Tells SPECFEM to save the resulting misfit kernel in binary format (not ASCII)

# Set the parameter file for an adjoint simulation
! seisflows sempar -P DATA/Par_file simulation_type 3
! seisflows sempar -P DATA/Par_file save_ASCII_kernels .false.

! mpirun -n 1 bin/xspecfem2D > OUTPUT_FILES/output_adjsolver.txt

# The adjoint solver also produces an informative log file
! tail OUTPUT_FILES/output_adjsolver.txt

# We can check this was an adjoint simulation because SPECFEM considers adjoint sources
! echo
! head -347 OUTPUT_FILES/output_adjsolver.txt | tail -n 5

5c) Understanding Adjoint Simulation Outputs#

Adjoint simulations generate kernel files in the same format as the Databases and model files
Kernels are generated for each quantity that defines the model

Event Kernel#

The most important output of the adjoint simulation is the event kernel
These files are named something like proc*_alpha_kernel.bin (here alpha==Vp)
Kernels define the volumentric integration of the interaction between the forward and adjoint wavefields
We’ll rename these kernels to make them more intutive (alpha \(\rightarrow V_p\); beta \(\rightarrow V_s\))

# Copy over the Model files so SeisFlows has access to coordinate information
! cp -r DATA/*.bin OUTPUT_FILES

# Rename alpha and beta kernels
! cp OUTPUT_FILES/proc000000_alpha_kernel.bin OUTPUT_FILES/proc000000_vp_kernel.bin 
! cp OUTPUT_FILES/proc000000_beta_kernel.bin OUTPUT_FILES/proc000000_vs_kernel.bin 

# Use SeisFlows to visualize the kernel
m = Model(path="OUTPUT_FILES")
m.plot2d("vs_kernel", show=False)

# NOTE: The below code snippet is repeated from Section 1, make sure that has been run
# Plot SOURCE 1, STATION 1 and a connecting line
plt.plot((sta_x[0], ev_x[0]), (sta_z[0], ev_z[0]), c="k", ls=":", lw=2, zorder=1)
plt.scatter(sta_x[0], sta_z[0], c="c", marker="v", s=100, edgecolor="k")
plt.scatter(ev_x[0], ev_z[0], c="y", marker="*", edgecolor="k", s=200)

Understanding the event kernel#

Constructed via the interaction between the forward wavefield and the adjoint wavefield
Illuminates sections of the model that the objective function is sensitive to
Note that the direct ray path (dotted line) has almost zero sensitivity (white)
The kernel is ‘fat’, meaninng sensitivity extends a finite width away from the direct raypath
Generated by a single source-receiver pair here, but in larger adjoint simulations, will combine many adjoint sources

SPECFEM2D Adjoint Images#

SPECFEM2D also produces adjoint images which document the adjoint wavefield
These do not show the interaction, but rather adjoint source emanating from the receiver location (i.e. column 2 in Tape et al. Fig 5)

Image("OUTPUT_FILES/adjoint_image000000200.jpg")

Image("OUTPUT_FILES/adjoint_image000000800.jpg")

6) Smoothing/Regularizing Kernels#

The kernel above is complicated, lots of fringes seen away from the actual kernel
Regularization is a process that changes results to be simpler
There are many instances in which you might want to regularize a kernel:
- If your data is band-limited, you may know that you don’t have sensitivity to small-scale perturbations, which we can remove to get at large-scale structure
- Due to numerical discretization, areas around the source and receiver will have very large-amplitude, high-frequency kernels that are not realistic
- Smooth away these small-scale perturbations to get at large-scale structure

In SPECFEM we can smooth kernels as a form of regularization

We use the SPECFEM executable xsmooth_sem to convolve our kernel with a 2D Gaussian
We choose the X and Z half-widths of this Gaussian

The usage of xsmooth_sem is given as

USAGE:  mpirun -np NPROC bin/xsmooth_sem SIGMA_H SIGMA_V KERNEL_NAME INPUT_DIR OUTPUT_DIR GPU_MODE

We will need to choose values and directories to make this work

SIGMA_H: Horizontal smoothing length [m] representing the horizontal half-width of the Gaussian
SIGMA_Z: Vertical smoothing length [m] representing the vertical half-width of the Gaussian
KERNEL_NAME: The name of the kernel we want to smooth. Must match filename, so proc000000_vs_kernel.bin corresponds to a matching kernel name: vs_kernel
INPUT_DIR: where SPECFEM can find the kernel files
OUTPUT_DIR: Where SPECFEM should output the SMOOTHED kernels
GPU_MODE: Use GPU acceleration to speed up the smoothing operation (.true. or .false.)

# Run the smoothinig operation
! mpirun -n 1 bin/xsmooth_sem 5000 5000 vs_kernel OUTPUT_FILES/ OUTPUT_FILES/ .false.

# Smoothed kernels are renamed <KERNEL_NAMEL>_smooth
! ls OUTPUT_FILES/*kernel*.bin

# Rename kernel files so we can plot the smooth version with SeisFlows
! mv OUTPUT_FILES/proc000000_vs_kernel.bin OUTPUT_FILES/proc000000_vs_kernel_raw.bin
! mv OUTPUT_FILES/proc000000_vs_kernel_smooth.bin OUTPUT_FILES/proc000000_vs_kernel.bin

# Use SeisFlows to visualize the kernel
m = Model("OUTPUT_FILES")
m.plot2d("vs_kernel", show=False)

# Plot source and receiver on the kernel figure
plt.plot((sta_x[0], ev_x[0]), (sta_z[0], ev_z[0]), c="k", ls=":", lw=2, zorder=1)
plt.scatter(sta_x[0], sta_z[0], c="c", marker="v", s=100, edgecolor="k")
plt.scatter(ev_x[0], ev_z[0], c="y", marker="*", edgecolor="k", s=200)

Comparing RAW and SMOOTHED kernels#

The amplitudes of the smoothed kernel are smaller than the raw kernel
The smoothed kernel looks visibly smoother (of course!)
Much more of the domain in the smoothed kernel shows amplitude of 0 (small-scale fringes removed)

High Performance Seismology

3. Kernels and Adjoint Simulations

Contents

3. Kernels and Adjoint Simulations#

0) Motivation#

1) Setting Up a SPECFEM2D Working Directory#

2) Experimental Setup: Tape et al. 2007 Example Problem#

2a) Homogeneous Halfspace Model#

2b) Perturbation Checkerboard Model#

Understanding the checkerboard#

2c) Visualizing Source-Receiver Geometry#

Understanding Source-Receiver Figures#

3) Setting up an Adjoint Simulation#

Useful Definitions#

Useful Information#

3a) Generating ‘Data’ with a Target Model#

Parameter changes explained#

3b) Generating Synthetics using Initial Model#

Parameter changes explained#

3c) Misfit Quantification and Adjoint Sources#

Comments on Default Waveform Figure#

4) Misfit Quantification#

4a) Misfit Calculation with Python#

4b) Adjoint Source with Python#

Understanding the waveform adjoint source figure#

5) Running an Adjoint Simulation#

5a) Writing Adjoint Sources#

5b) Run `xspecfem2D`#

Parameter file changes#

5c) Understanding Adjoint Simulation Outputs#

Event Kernel#

Understanding the event kernel#

SPECFEM2D Adjoint Images#

6) Smoothing/Regularizing Kernels#

Comparing RAW and SMOOTHED kernels#

High Performance Seismology

3. Kernels and Adjoint Simulations

Contents

3. Kernels and Adjoint Simulations#

0) Motivation#

1) Setting Up a SPECFEM2D Working Directory#

2) Experimental Setup: Tape et al. 2007 Example Problem#

2a) Homogeneous Halfspace Model#

2b) Perturbation Checkerboard Model#

Understanding the checkerboard#

2c) Visualizing Source-Receiver Geometry#

Understanding Source-Receiver Figures#

3) Setting up an Adjoint Simulation#

Useful Definitions#

Useful Information#

3a) Generating ‘Data’ with a Target Model#

Parameter changes explained#

3b) Generating Synthetics using Initial Model#

Parameter changes explained#

3c) Misfit Quantification and Adjoint Sources#

Comments on Default Waveform Figure#

4) Misfit Quantification#

4a) Misfit Calculation with Python#

4b) Adjoint Source with Python#

Understanding the waveform adjoint source figure#

5) Running an Adjoint Simulation#

5a) Writing Adjoint Sources#

5b) Run xspecfem2D#

Parameter file changes#

5c) Understanding Adjoint Simulation Outputs#

Event Kernel#

Understanding the event kernel#

SPECFEM2D Adjoint Images#

6) Smoothing/Regularizing Kernels#

Comparing RAW and SMOOTHED kernels#

5b) Run `xspecfem2D`#