Skip to content

Download 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.

Benchmark Salvus vs Instaseis

This notebook creates a benchmark of a global simulation using Salvus with a reference solution extracted from Instaseis.

# Python libraries.
import os
import copy
import json
import toml
import shutil
import numpy as np
import pandas as panda

from pathlib import Path
import matplotlib.pyplot as plt

%matplotlib inline
# Seismology toolbox.
import pyasdf
import instaseis
# Salvus modules
import salvus_seismo

from salvus_mesh.mesh.run_mesher import run_mesher
from salvus_mesh.interface.main import invoke_subcommand
from salvus_seismo.source_time_function import Ricker, GaussianRate

Input parameters

instaseisDB = 'prem_light_200s_new_repacked'
model = 'prem_iso_no_crust'
period = 200
salvus_model_order = 2

# choose gaussian or ricker here
salvus_stf = 'gaussian'   

# Parameters for our source.
src_depth_in_m = 10e3
src_lat, src_lon = 0., 0.
m_rr, m_tt, m_pp, m_rt, m_rp, m_tp = 1e20, 0.0, 0.0, 0.0, 0.0, 0.0

# source time function
src_decay_rate = 3.5
src_time_shift = 2 * period
src_half_duration = period

# receiver locations
num_rec = 60
rec_lon = np.linspace(0., 180., num_rec, endpoint=True)
rec_lat = np.zeros(rec_lon.shape)


output_folder = 'out_' + salvus_stf

Mesh generation

# Get default settings.
mesh_inputs = invoke_subcommand('Globe3D', [])[0]

# Customize for our comparison.
mesh_inputs['basic']['model'] = model
mesh_inputs['basic']['period'] = period
mesh_inputs['advanced']['tensor_order'] = salvus_model_order

mesh_filename = "$mesh_type$_$model_name$_$period$.h5"

# Create mesh
m, sk_info, mesh_filename = run_mesher(
    mesh_inputs, output_filename=mesh_filename, verbose=True, overwrite_file=True)

m.write_h5_tensorized_model(mesh_filename)

Prepare Salvus run

def generate_src_rec(mesh_file):
    """
    Generate source and receiver dictionaries.
    :param mesh_file: Mesh file to attach sources and receivers to.
    """

    # Switch sources to new Salvus format.
    if salvus_stf == "gaussian":
        stf = GaussianRate(
            decay_rate=src_decay_rate, 
            half_duration_in_seconds=src_half_duration,
            time_shift_in_seconds=src_time_shift)
    elif salvus_stf == "ricker":
        stf = Ricker(
            time_shift_in_seconds=src_time_shift,
            center_frequency=src_decay_rate/(np.pi * period))        

    # Generate our source.
    sources = [salvus_seismo.Source(source_time_function=stf,
        latitude=src_lat, longitude=src_lon, depth_in_m=src_depth_in_m, 
        m_rr=m_rr, m_tt=m_tt, m_pp=m_pp, m_rt=m_rt, m_rp=m_rp, m_tp=m_tp)]

    # Generate receivers.
    receivers = []
    for i in range(0,len(rec_lon)):
        receivers.append(salvus_seismo.Receiver(
            latitude=rec_lat[i], longitude=rec_lon[i], depth_in_m=0., station='%03d' % i))

    # Attach sources and receivers to the mesh.
    src_seismo, rec_seismo = salvus_seismo.generate_sources_and_receivers(
        sources=sources, receivers=receivers, mesh=mesh_file, surface_side_set="r1")

    return src_seismo, rec_seismo
def generate_salvus_input(mesh_file, src, rec):
    """
    Generate a dictionary which represents a Salvus input file.
    :param mesh_file: Mesh file to run simulation on.
    :param src: Sources to attach.
    :param rec: Receivers to attach.
    """

    # Basic salvus input groups.
    salvus_input = {

        # Compute architecture.
        "hardware": {
            "gpu": False
        },

        # Simulation domain.
        "domain": {
            "dimension": 3,
            "polynomial-order": 4,
            "mesh": {
                "filename": mesh_file,
                "format": "hdf5"
            },            
            "model": {
                "filename": mesh_file,
                "format": "hdf5"
            },
            "geometry": {
                "filename": mesh_file,
                "format": "hdf5"
            }
        },

        # Wave propagation physics.
        "physics": {
            "wave-equation": {
                "start-time-in-seconds": 0.,
                "end-time-in-seconds": 1000.,
                "time-stepping-scheme": "newmark",
                "attenuation": False
            }
        },

        # Output.
        "output": {
            "point-data": {
                "format": "asdf",
                "sampling-interval-in-time-steps": 1,
                "filename": os.path.splitext(mesh_file)[0] + "_receiver.h5"
            }
        }
    }

    # Add sources and receivers.
    salvus_input["physics"]["wave-equation"]["point-source"] = src["point-source"]
    salvus_input["output"]["point-data"].update(rec["output"]["point-data"])

    return salvus_input
# Create the input files, and attach sources and receivers to the mesh.
salvus_inputs = generate_salvus_input(mesh_filename, *generate_src_rec(mesh_filename))
with open('run.toml', "wt") as fh:
    toml.dump(salvus_inputs, fh)

Salvus simulation

If Salvus_Flow_ is configured on your system, you can launch a simulation on your local machine or on a remote cluster using the command below. Alternatively, manually upload the mesh file and the toml file to your target machine and launch a job. Afterwards download the output data (meta.json and the receiver file) into a subdirectory out on your local machine.

import salvus_flow.api
salvus_flow.api.run(
    site_name="local", output_folder=output_folder, wall_time_in_seconds=1000,
    input_file=salvus_inputs, ranks=4, overwrite=True)

Compare to Instaseis solution

# Open the pre-computed Instaseis database.
db = instaseis.open_db(instaseisDB)
# source time function in time domain
def gauss(t, t0, tshift=0., decay=3.5):
    return decay / (t0 * np.pi ** 0.5) * np.exp(-(decay / t0 * (t - tshift)) ** 2)
# Ensure that we will use the same source in instaseis.
assert(len(salvus_inputs['physics']['wave-equation']['point-source']) == 1)
src = salvus_inputs['physics']['wave-equation']['point-source'][0]
assert(src['spatial-type'] == "moment_tensor")
t_source = np.arange(db.info.npts) * db.info.dt
stf = gauss(t_source, src_half_duration, tshift=src_time_shift, decay=src_decay_rate)

plt.plot(t_source, stf)
plt.show()

# Generate the Instaseis source object.
source = instaseis.Source(
    latitude=src_lat, longitude=src_lon, depth_in_m=src_depth_in_m,
    m_rr=m_rr, m_tt=m_tt, m_pp=m_pp, m_rt=m_rt, m_rp=m_rp, m_tp=m_tp,
    sliprate=stf, dt=db.info.dt)
def get_traces(station_idx, output_folder, fmin, fmax):
    station = f"XX_{station_idx:03d}"
    rec_file = Path(output_folder) / str(salvus_inputs['output']['point-data']['filename'])
    with pyasdf.ASDFDataSet(rec_file, mode="r") as ds:
        sd = ds.waveforms[station].StationXML[0][0]
        st_salvus = ds.waveforms[station].displacement
        dt = ds.waveforms[station].displacement[0].stats.delta        

    # Query instaseis database.
    station_lat, station_lon = sd.latitude, sd.longitude
    receiver = instaseis.Receiver(latitude=station_lat, longitude=station_lon)
    st_instaseis = db.get_seismograms(
        source, receiver, components='ZNE', remove_source_shift=False, reconvolve_stf=True, 
        dt=dt, kind='velocity')

    # Filter traces.
    st_instaseis.filter('lowpass', freq=fmax, corners=8, zerophase=True)
    st_instaseis.filter('highpass', freq=fmin, corners=8, zerophase=True)

    st_salvus.filter('lowpass', freq=fmax, corners=8, zerophase=True)
    st_salvus.filter('highpass', freq=fmin, corners=8, zerophase=True)

    if (salvus_stf == 'ricker'):
        # if the salvus stf is a ricker, we need to integrate the 
        # output twice to obtain velocity
        # Note that the salvus stf currently specifies slip, not slip rate
        st_instaseis.filter('highpass', freq=fmin, corners=8, zerophase=True)
        st_salvus.integrate()
        st_salvus.integrate()
        st_salvus.filter('highpass', freq=fmin, corners=8, zerophase=True)

        for tr in st_salvus:
            # correct amplitude for ricker source
            alpha = src_decay_rate / period
            amp_gauss = alpha / (np.pi ** 0.5)
            alpha_salvus = alpha
            fac = - 2 * alpha_salvus ** 2 * amp_gauss
            tr.data *= fac

    return dt, st_salvus, st_instaseis

def plot_station(station_idx, output_folder, fmin, fmax):

    dt, st_salvus, st_instaseis = get_traces(station_idx, output_folder, fmin, fmax)
    # Plot the data.
    f, axes = plt.subplots(3, 1, sharex=True, sharey=True, figsize=(15, 10))
    tab = []
    for c, a in zip(["Z", "N", "E"], axes):
        tr_salvus = st_salvus.select(component=c)[0]
        tr_instaseis = st_instaseis.select(component=c)[0].trim(endtime=tr_salvus.stats.endtime)
        tr_diff = tr_salvus.data - tr_instaseis.data        
        a.plot(tr_salvus.times(), tr_salvus.data, label='Salvus')
        a.plot(tr_instaseis.times(), tr_instaseis.data, label='Instaseis')
        a.set_ylabel("Velocity (m/s)")
        a.legend()
        l2_error_abs = np.linalg.norm(dt * tr_diff)
        l2_error_rel = l2_error_abs / np.linalg.norm(dt * tr_instaseis.data)        
        linf_error = np.max(tr_diff)
        tab.append({'component': c, 
                    'L2 error (abs)': l2_error_abs, 
                    'L2 error (rel)': l2_error_rel, 
                    'Linf error' : linf_error})              
#         assert(linf_error < 1e-8)        
#         assert(l2_error_rel < 0.02 or l2_error_abs < 1e-10)

    axes[-1].set_xlabel("Time (s)")    
    return f, axes


def compute_error(station_idx, output_folder, fmin, fmax):

    dt, st_salvus, st_instaseis = get_traces(station_idx, output_folder, fmin, fmax)
    # Plot the data.
    tab = []
    for c, a in zip(["Z", "N", "E"], axes):
        tr_salvus = st_salvus.select(component=c)[0]
        tr_instaseis = st_instaseis.select(component=c)[0].trim(endtime=tr_salvus.stats.endtime)
        tr_diff = tr_salvus.data - tr_instaseis.data        
        l2_error_abs = np.linalg.norm(dt * tr_diff)
        l2_error_rel = l2_error_abs / np.linalg.norm(dt * tr_instaseis.data)        
        linf_error = np.max(tr_diff)
        tab.append({'component': c, 
                    'L2 error (abs)': l2_error_abs, 
                    'L2 error (rel)': l2_error_rel, 
                    'Linf error' : linf_error})              
        assert(linf_error < 1e-8)        
        assert(l2_error_rel < 0.02 or l2_error_abs < 1e-10)
    return panda.DataFrame(tab)        
# Open the output file from the solver so we can read the simulation parameters.
with open(Path(output_folder) / "meta.json", mode="r") as f:
    metaInfo = json.load(f)

f_max = 1 / period
f_min = 0.2 * f_max
f, axes = plot_station(0, output_folder, f_min, f_max)
plt.show()
# Open the output file from the solver so we can read the simulation parameters.
with open(Path(output_folder) / "meta.json", mode="r") as f:
    metaInfo = json.load(f)

f_max = 1 / period
f_min = 0.2 * f_max
for s in range(10,26):
    df = compute_error(s, output_folder, f_min, f_max)