2. Forward Simulations#

  • Here we build upon material learned in Notebook 1

  • This notebook allows Users to play around with their own SPECFEM2D homogeneous halfspace example in an exercise

  • Objective: Familiarize Users with setting SOURCE and STATION attributes, adjusting velocity model parameters, and assessing simulation results.

  • These instructions should be run from inside a Docker container, using Jupyter Lab (see instructions here).

Relevant Links:

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.

  • To time a task, put a %time before the command (e.g., %time ! ls)

1) Set Up a SPECFEM Working Directory#

  • It is often desireable to run SPECFEM outside of the cloned repository, in order to keep files and outputs manageable

  • The trick here is that SPECFEM only requires 3 compenents for a sucessful simulation: bin/, DATA/, and OUTPUT_FILES/

# Python packages required for this notebook
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image

# Create the correct dir. and move there
! mkdir -p /home/scoped/work/specfem2d_workdir
%cd /home/scoped/work/specfem2d_workdir

# Symlink the executables, copy example DATA/, create empty OUTPUT_FILES
! 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#

! ls DATA/

DATA/ Directory#

  • Par_file for a homogeneous halfspace model in Par_file_Tape2007_onerec

  • Par_file for a checkerboard model in Par_file_Tape2007_132rec_checker

  • Mesh files in: interfaces_Tape2007.dat and the Par_file_*

  • Model files in: model_velocity.dat_checker

  • Source files in: the 25 SOURCE_??? files

  • Stations: in the STATIONS_checker file

2a) The Homogeneous Halfspace Model#

  • The homogeneous halfspace model in this example is defined internally, in the Par_file

  • A homogeneous halfspace defines single set of properties for the entire domain

  • We can also use the utility seisflows sempar velocity_model command to look at model values

# Bash commands to look at the Par_file
! head -287 DATA/Par_file_Tape2007_onerec | tail -16
! echo "..."

# SeisFlows utility function to look at the Par_file
! seisflows sempar -P DATA/Par_file_Tape2007_onerec velocity_model

Understanding the Velocity Model#

According to the Par_file comments, the model parameter values represent the following:

model_number 1 rho Vp Vs 0 0 QKappa Qmu  0 0 0 0 0 0
1 1 2600.d0 5800.d0 3500.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0

  • The homogeneous halfspace model defines a region with P-wave velocity \(V_p=5.8\)km/s and S-wave velocity \(V_s=3.5\)km/s.

  • The halfspace is also defined by density and attenuation

  • We can understand the structure of the mesh by looking at the Par_file and the interfaces_Tape2007.dat file.

# Look at Mesh parameters to view the size of the domain
! head -306 DATA/Par_file_Tape2007_onerec | tail -n 15

# Look at interface parameters 
! cat DATA/interfaces_Tape2007.dat

Understanding the Mesh parameters#

  • From the files above, we can see that the X and Z dimensions of our mesh range from 0 to 480000.0m

  • Each dimension is separated into 40 elements (defined by nxmin, nxmax etc. in the Par_file and defined by the layer numbers in the interfaces file)

  • That means each spectral element in our 2D mesh spans: 480000m / 40 elements = 12000m / element (or 12km / element)

  • Also note that above we previously learned that the \(V_s\) model has a homogeneous value of 3.5 km/s

Visualizing the Model#

We can make a simple plot using Matplotlib to illustrate what our mesh might look like

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, 480001, 4000)
    z = np.arange(0, 480001, 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.title("2D Homogeneous Halfspace Model\n Vs=3.5km/s")
    
# Calls the function we just defined
plot_homogeneous_halfspace()

# Plot grid lines representing each of the spectral elements
for i in range(12000, 480000, 12000):
    plt.axvline(i, c="k", lw=0.5)
    plt.axhline(i, c="k", lw=0.5)

2b) Visualizing Source-Receiver Geometry#

  • We can similarly plot the SOURCES and STATIONS available to see what the experiemental setup looks like

  • This is the same Python-based approach we took in the Day 1A notebook

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

# Grab coordinates from each SOURCE file
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 together. Annotate names
plot_homogeneous_halfspace()
plt.scatter(ev_x, ev_z, c="y", marker="*", s=100, edgecolor="k")
plt.scatter(sta_x, sta_z, c="c", marker="v", s=20, edgecolor="k")
plt.title("SOURCE-RECEIVER GEOMETRY")

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

# Plot STATIONS with their names
# Because STATIONS are so close, numbers will be jumbled.
plot_homogeneous_halfspace()
for i, (x, z) in enumerate(zip(sta_x, sta_z)):
    plt.scatter(x, z, c="c", marker="v", s=12, edgecolor="k")
    plt.text(x, z, f"{i:0>3}", fontsize=9)
plt.title(f"STATIONS; N={len(sta_x)}")

  • Upside-down blue triangles represent the 132 receivers in this example

  • The 25 yellow stars are the sources.

3) Running SPECFEM2D#

  • Before we run the example, we need to do some organizational bookkeeping

  • We will choose ONE of our source files to act as our source

  • We will choose ONE stations to act as our station file

3a) Choosing a Source file#

  • SPECFEM2D will look for a file named SOURCE in the DATA/ directory to define its source

# Choose SOURCE_001 as our SOURCE File
! cp -f DATA/SOURCE_001 DATA/SOURCE

# > Make sure that the SOURCE name printed below matches choice
! head -1 DATA/SOURCE

3b) Choosing Stations#

  • SPECFEM2D will look for a file named STATIONS in the DATA/ directory to define its stations

  • The STATIONS_checker file defines 132 different station locations, we only want one

  • Remember: Individual synthetic seismograms simply extract the simulated wavefield at a location, i.e., computational expense is not tied to the number of stations.

# Write out a NEW stations file by choosing station numbers
! head -1 DATA/STATIONS_checker > DATA/STATIONS
! cat DATA/STATIONS

3c) Setting up the Par_file#

  • SPECFEM2D will look for a file called Par_file in the DATA/ directory to set its parameters

  • We will copy over the Par_file_Tape2007_onerec to define our parameter file

  • We need to change a few key parameters in the Par_file to run SPECFEM2D with desired behavior

  • We use the seisflows sempar command to make the changes but this can be done with a text editor, Bash etc.

# Copy in the Example parameter file
! cp -f DATA/Par_file_Tape2007_onerec DATA/Par_file

# Set some necessary parameters for later in the Par_file
! seisflows sempar -P DATA/Par_file nproc 4
! seisflows sempar -P DATA/Par_file use_existing_stations .true.

Understanding Parameter Changes#

NPROC: Sets the number of MPI processors to partition the mesh and run the simulation with. This must match the value following -n in the MPI
use_existing_STATIONS: Use the STATIONS file we created, as opposed to the Par_file definition of stations

3d) Run SPECFEM#

  • Now that we have set the Par_file, the SOURCE and STATIONS file, we are able to run xmeshfem2D and xspecfem2D to run our forward simulation.

  • We use 4 MPI processes to run this homogeneous halfspace simulation

  • We expect only one synthetic seismogram to be output from this simulation

# Ensures we're running with a clean OUTPUT directory
! rm -rf OUTPUT_FILES
! mkdir OUTPUT_FILES

! mpirun -n 4 bin/xmeshfem2D > OUTPUT_FILES/output_meshfem.txt
! mpirun -n 4 bin/xspecfem2D > OUTPUT_FILES/output_solver.txt

! tail OUTPUT_FILES/output_solver.txt

3e) Examine Output Files#

  • Let’s confirm that we have created one displacement seismogram

  • Then we’ll look at the forward simulation figures to see if things make sense

! ls OUTPUT_FILES/
! echo
! ls OUTPUT_FILES/*.semd

# We can use SeisFlows to plot our waveform
! seisflows plotst OUTPUT_FILES/AA.S000000.BXY.semd --savefig AA.S000000.BXY.png
Image("AA.S000000.BXY.png")

# We can also look at the wavefield snapshots
Image("OUTPUT_FILES/forward_image000000800.jpg")

# We can also look at the wavefield snapshots
Image("OUTPUT_FILES/forward_image000001200.jpg")

# We can also look at the wavefield snapshots
Image("OUTPUT_FILES/forward_image000002200.jpg")

4) Forward Simulation Exercise#

  • Participants will now be asked to edit simulation parameters to run their own simulation

  • Some things that you are asked to try include:

    1. Change the parameters of the homogeneous halfspace model defined in the Par_file

    2. Define a STATIONS file with a uniform grid of stations to record synthetics throughout the domain

    3. Choose a different source, or increase the energy released by the source (using the moment tensor)

    4. Re-run the mesher and solver to get new synthetics

    5. Analyze the new results in comparison to the old results

  • First we set up a working directory for you

! rm -rf /home/scoped/work/exercise_1
! mkdir -p /home/scoped/work/exercise_1
%cd /home/scoped/work/exercise_1

# Symlink the executables, copy example DATA/, create empty OUTPUT_FILES
! ln -s /home/scoped/specfem2d/bin .
! cp -r /home/scoped/specfem2d/EXAMPLES/Tape2007/DATA .
! mkdir OUTPUT_FILES

# Set the Par_file
! cp DATA/Par_file_Tape2007_onerec DATA/Par_file

! ls

Task 1: Edit the Velocity Model#

  • Change the velocity model parameters in the homogeneous halfspace model

  • Remember, the velocity model is defined in the Par_file

  • Try increasing seismic velocity (Vp and Vs) by 10%

  • You can use Python, Bash, seisflows sempar or a Text Editor to do this

Task 2: Create a New STATIONS File#

  • Define a STATIONS file that covers the entire domain with a uniform grid spacing of:

    • dx = 80km

    • dz = 80km

    • x_start = 0km

    • z_start = 0km

  • Or Create your own station configuration. Some examples: spiral, concentric rings, dense linear array (like DAS)

  • You can find the X and Z dimensions of the mesh in the Par_file and the interfaces file, respectively

  • Use Python/NumPy to loop values, or simply write out a text file manually with the text editor

  • Look at DATA/STATIONS_checker for an example of how the file should look

  • NOTE: The last two columns (burial, elevation) can be set to 0

Task 3: Choose and edit a SOURCE file#

  • Use one of the original sources as a template for your new source

  • Set the location of your source in the exact middle of your domain (or a location of your choice!)

  • Set the moment tensor (Mxx, Mzz, Mxz) of your event to make this an explosive source (or a mechanism of your choice!)

  • Don’t change the scaling on the moment tensor

Task 4: Run the Solver and Analyze Outputs#

  • Run the mesher and solver with your new experimental setup and 4 MPI processes

  • Remember to tell SPECFEM to use your STATIONS file and not its internal representation of stations

  • Remember to tell SPECFEM that we want to run this with 4 processors

  • Look at the source images to see if your explosion makes sense

  • Plot waveforms output from your gridded stations