Skip to content

Instaseis 3D Globe

Jupyter Notebook

This page is automatically generated from a Jupyter Notebook.

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

Click here for a version that is compatbile with the 0.9.X releases.

Introduction

In this notebook we will set up and run a global-scale seismic simulation, and test the simulation computed by Salvus against a semi-analytic solution computed with Instaseis. We'll attach realistic Earthquake sources and stations, and analyze point (receiver) and surface (movie) output. Before you begin this tutorial, ensure that Salvus is installed on the machine where you will run the simulations, and that Salvus Flow is properly set up to communicate with that machine.

Environment

%config Completer.use_jedi = False
%matplotlib inline

# Generic Python libraries.
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
import os

# Packages to assist with Salvus's preferred way of representing seismic data.
import obspy
import pyasdf
import instaseis
from obspy import UTCDateTime
from obspy.clients.syngine import Client

# Salvus flow objects to help set up simulations.
from salvus_flow.simple_config import source, receiver, simulation, stf

# Helper to generate a mesh.
from salvus_mesh.simple_mesh.basic_mesh import GlobalBuiltIn3D

# General Salvus flow API.
import salvus_flow.api

# General parameters.
period = 100.0
src_lat = 3.0
src_lon = 11.0
rec_lat = 7.0
rec_lon = 82.0

attenuation = bool(int(os.environ.get('ATTENUATION', True)))
anisotropic = bool(int(os.environ.get('ANISOTROPY', False)))
instaseis_db_root = str(os.environ.get('INSTASEIS_ROOT', '.'))

Mesh

# Built-in model to use for a comparison.
if anisotropic:
    model = "prem_ani_no_crust"
else:
    model = "prem_iso_no_crust"

# Generate the Salvus mesh.
mesh = GlobalBuiltIn3D(
    model=model, period=period, tensor_order=4
).create_mesh()

# Visualize the mesh in the notebook.
mesh

<salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x1c0f6b5080>

Source and receivers

# Generic source parameters.
decay_rate = 3.5
half_duration = period
time_shift = 2 * period

# Create an explosive source to attach to our mesh.
source = source.seismology.MomentTensorPoint3D(
    latitude=src_lat,
    longitude=src_lon,
    mrr=0.0,
    mtt=-1e21,
    mpp=1e21,
    mrt=0.0,
    mrp=0.0,
    mtp=0.0,
    depth_in_m=0.0,
    source_time_function=stf.GaussianRate(
        half_duration_in_seconds=half_duration,
        decay_rate=decay_rate,
        time_shift_in_seconds=time_shift,
    ),
)

# Create a receiver to attach to the mesh.
# Save velocity, displacement, and strain at the receiver location.
stations = receiver.seismology.Point3D(
    latitude=rec_lat,
    longitude=rec_lon,
    depth_in_m=0.0,
    station_code="XX",
    fields=["displacement", "velocity", "strain"],
)

Simulate

# Simulation parameters.
start_time = 0.0
end_time = 3600.0
# Initialize.
w = simulation.Waveform(mesh=mesh, sources=source, receivers=stations)

w.physics.wave_equation.time_step_in_seconds = 0.025
w.physics.wave_equation.start_time_in_seconds = start_time
w.physics.wave_equation.end_time_in_seconds = end_time
w.physics.wave_equation.attenuation = attenuation

# Make a movie.
#     w.output.surface_data.sampling_interval_in_time_steps = 5000
#     w.output.surface_data.fields = ["displacement"]
#     w.output.surface_data.filename = "output.h5"
#     w.output.surface_data.format = "hdf5"
#     w.output.surface_data.side_sets = ["r1"]
w.validate()

salvus_flow.api.run(
    ranks=1200,
    site_name="cluster",
    input_file=w,
    output_folder="output",
    get_all=True,
    overwrite=True,
    wall_time_in_seconds=1200,
)
Job `job_1904261024_43000` running on `cluster` with 1200 rank(s).
Site information:
  * Salvus version: 0.8.2-62-gf59b5c1
  * Floating point size: 32

Analyze

def gauss(
    *,
    t_array: np.ndarray,
    half_duration: float,
    t_shift: float = 0.0,
    decay: float = 3.5
):
    """
    Generate a Gaussian.
    :param t_array: Time samples in seconds at which to sample Gaussian.
    :param half_duration: Half-duration of the Gaussian.
    :param t_shift: Time shift of the Gaussian in seconds.
    :param decay: Decay rate.
    :returns: A numpy array.
    """
    return (
        decay
        / (half_duration * np.pi ** 0.5)
        * np.exp(-(decay / half_duration * (t_array - t_shift)) ** 2)
    )
def filter(*, tr: obspy.Stream, min_period: float, max_period: float):
    """
    Apply as simple filter to the recorded seismograms.
    :param min_period: Minimum period for the bandpass filter.
    :param max_period: Maxiumum period for the bandpass filter.
    :returns: A copy of a filtered obspy trace.
    """
    tr0 = tr.copy().filter(
        "lowpass", freq=1 / min_period, corners=8, zerophase=True
    )
    return tr0.filter("lowpass", freq=1 / min_period, corners=8, zerophase=True)
# Open database and get source-time function.
if attenuation and anisotropic:
    db_name = "prem_ani_light_20s_att"
elif attenuation and not anisotropic:
    db_name = "prem_iso_light_20s_att"
elif not attenuation and anisotropic:
    db_name = "prem_ani_light_20s_noatt"
else:
    db_name = "prem_iso_light_20s_noatt"

db_file = Path(instaseis_db_root, db_name)
db = instaseis.open_db(str(db_file))
t_source = np.arange(db.info.npts) * db.info.dt
stf = gauss(t_array=t_source, half_duration=period, t_shift=time_shift)

# Generate the Instaseis source object.
source = instaseis.Source(
    latitude=src_lat,
    longitude=src_lon,
    depth_in_m=0.0,
    m_rr=0.0,
    m_tt=-1e21,
    m_pp=1e21,
    m_rt=0,
    m_rp=0,
    m_tp=0,
    sliprate=stf,
    dt=db.info.dt,
)

# Generate a receiver.
receiver = instaseis.Receiver(latitude=rec_lat, longitude=rec_lon)
st_instaseis = db.get_seismograms(
    source,
    receiver,
    components="ENZ",
    remove_source_shift=False,
    reconvolve_stf=True,
    dt=db.info.dt,
    kind="velocity",
).trim(UTCDateTime(start_time), UTCDateTime(start_time) + end_time)
output_dir = Path("output")

with pyasdf.ASDFDataSet(output_dir / "receivers.h5") as fh:
    st_salvus = fh.waveforms["XX.XX"].velocity
min_period_filter = 1 * period
max_period_filter = 1000.0
f, ax = plt.subplots(nrows=3, ncols=1, figsize=(15, 10))

ax[0].plot(
    st_instaseis[0].times(),
    filter(
        tr=st_instaseis[0],
        min_period=min_period_filter,
        max_period=max_period_filter,
    ).data,
    label="Instaseis",
)
ax[0].plot(
    st_salvus[0].times(),
    filter(
        tr=st_salvus[0],
        min_period=min_period_filter,
        max_period=max_period_filter,
    ).data,
    label="Salvus",
)
ax[0].set_title("East")
ax[0].legend()

ax[1].plot(
    st_instaseis[0].times(),
    filter(
        tr=st_instaseis[2],
        min_period=min_period_filter,
        max_period=max_period_filter,
    ),
    label="Instaseis",
)
ax[1].plot(
    st_salvus[0].times(),
    filter(
        tr=st_salvus[2],
        min_period=min_period_filter,
        max_period=max_period_filter,
    ).data,
    label="Salvus",
)
ax[1].set_title("Vertical")
ax[1].legend()

ax[2].plot(
    st_instaseis[0].times(),
    filter(
        tr=st_instaseis[1],
        min_period=min_period_filter,
        max_period=max_period_filter,
    ).data,
    label="Instaseis",
)
ax[2].plot(
    st_salvus[0].times(),
    filter(
        tr=st_salvus[1],
        min_period=min_period_filter,
        max_period=max_period_filter,
    ).data,
    label="Salvus",
)
ax[2].set_title("North")
ax[2].legend()
plt.savefig('{0}.png'.format(db_name))
plt.show()

png