Mondaic - Full Waveform Solutions

This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

Download Jupyter Notebook [124 kB]

This integration test compares seismograms in a 2D elastic domain with an ocean layer on top for three different scenarios:

The ocean layer is completely ignored and replaced by a free surface
The ocean layer is approximated by a load applied to the surface of the elastic domain
The ocean is modelled as an acoustic medium in a coupled solid/fluid simulation.

Depending on the height (weight) of the water column, we expect the solutions of [1] and [2] to be similar for (very) thin ocean layers. Furthermore, the solutions to [2] and [3] are similar, if the height of the water column corresponds to a fraction of a wavelength (here: 10%).

The ocean load is modelled as a boundary condition, which is readily available from simple_config. There are two required steps to run a simulation with ocean loading:

Attach a scalar field OCEANLOAD to the mesh, which contains the weight of the water column, i.e., the ocean height times its density. The field is only required (and read) at the ocean surface. However, we currently don't have the concept of material properties at interface, so it still needs to have the same size as all other material parameters.
Add ocean loading as boundary condition to the waveform simulation object:

import salvus.flow.simple_config as sc
sim = sc.simulation.Waveform(mesh=mesh)
oceanload = sc.boundary.OceanLoading(side_sets=[...])
sim.add_boundary_conditions(oceanload)

Copy

%matplotlib inline
%config Completer.use_jedi = False
import os
import toml
import numpy as np
import matplotlib.pyplot as plt
import salvus.flow.api
from pyasdf import ASDFDataSet
from salvus.mesh.simple_mesh import basic_mesh
import salvus.flow.simple_config as sc

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
RANKS_PER_JOB = 2

Copy

# Domain setup (m).
max_x = 2500.0  # Distance in meters.
max_y = 700.0  # Distance in meters.
max_frequency = 40.0  # Frequency in Hz.

vp_fluid = 1500.0
rho_fluid = 1200.0

Copy

def generate_meshes(ocean_height):
    header = [
        "NAME         ocean_benchmark\n",
        "UNITS        km\n",
        "COLUMNS      radius rho vp vs\n",
    ]

    solid_model = [
        "    0.7   2.6     1.7     1.0\n",
        "    0.0   7.6     5.1     3.0\n",
    ]

    ocean_layer = [
        f"    {1 + ocean_height/1000} {rho_fluid/1000}     {vp_fluid/1000}     0.0\n"
        f"    0.7   {rho_fluid/1000}     {vp_fluid/1000}     0.0\n"
    ]

    with open("solid.bm", "w") as f:
        f.writelines(header)
        f.writelines(solid_model)

    with open("coupled.bm", "w") as f:
        f.writelines(header)
        f.writelines(ocean_layer)
        f.writelines(solid_model)

    mesh = basic_mesh.CartesianFromBm2D(
        bm_file="solid.bm",
        x_max=max_x,
        y_max=max_y,
        max_frequency=max_frequency,
    ).create_mesh()
    mesh_coupled = basic_mesh.CartesianFromBm2D(
        bm_file="coupled.bm",
        x_max=max_x,
        y_max=max_y + ocean_height,
        max_frequency=max_frequency,
    ).create_mesh()

    mesh.attach_field(
        "OCEANLOAD",
        ocean_height
        * rho_fluid
        * np.ones_like(mesh.element_nodal_fields["VP"]),
    )

    return mesh, mesh_coupled

Copy

# Sources.
fx, fy = 1e9, -1e10
sx, sy = 500.0, 500.0
stf = sc.source.stf.Ricker(center_frequency=14.5)
source = sc.source.cartesian.VectorPoint2D(
    x=sx, y=sy, fx=fx, fy=fy, source_time_function=stf
)

# Receivers.
nr = 10  # Number of receivers.
ry0 = 690.0  # Receiver y-value.
rx0 = 1500.0  # x-value of first receiver.
rx1 = 2000.0  # x-value of last receiver.
receivers = [
    sc.receiver.cartesian.Point2D(
        x=x, y=ry0, station_code=f"{_i:03d}", fields=["displacement"]
    )
    for _i, x in enumerate(np.linspace(rx0, rx1, nr))
]

Copy

def get_simulation_template(mesh, output=False):
    sim = sc.simulation.Waveform(
        mesh=mesh, sources=source, receivers=receivers
    )

    sim.physics.wave_equation.end_time_in_seconds = 0.75

    if output:
        sim.output.volume_data.filename = "wavefield.h5"
        sim.output.volume_data.format = "hdf5"
        sim.output.volume_data.fields = ["displacement"]
        sim.output.volume_data.sampling_interval_in_time_steps = 50

    boundaries = sc.boundary.Absorbing(
        side_sets=["x0", "x1", "y0"], taper_amplitude=0.0, width_in_meters=0.0
    )
    sim.add_boundary_conditions(boundaries)
    sim.validate()

    return sim

Copy

def compare_seismograms(asdf_file, asdf_reference, freqmax, freqmin):
    with ASDFDataSet(asdf_file, mode="r") as data, ASDFDataSet(
        asdf_reference, mode="r"
    ) as ref:
        total_misfit = 0
        for u, uref in zip(data.waveforms, ref.waveforms):
            tr_misfit = []
            for j, cmp in enumerate(["X", "Y"]):
                dt = 1e-4
                tr = u.displacement.select(component=cmp)[0]
                tr.interpolate(method="lanczos", a=12, sampling_rate=1 / dt)
                tr1 = tr.copy().filter(
                    "bandpass",
                    freqmax=freqmax,
                    freqmin=freqmin,
                    corners=8,
                    zerophase=True,
                )
                tr = uref.displacement.select(component=cmp)[0]
                tr.interpolate(method="lanczos", a=12, sampling_rate=1 / dt)
                tr2 = tr.copy().filter(
                    "bandpass",
                    freqmax=freqmax,
                    freqmin=freqmin,
                    corners=8,
                    zerophase=True,
                )

                diff = (tr1.data - tr2.data) ** 2
                tr_misfit.append(np.sqrt(diff.sum() / (tr2.data**2).sum()))
            total_misfit += (tr_misfit[0] + tr_misfit[1]) / 2
        avg_misfit = total_misfit / len(data.waveforms)

    return avg_misfit

Copy

ocean_height = 3.75

mesh, mesh_coupled = generate_meshes(ocean_height)
sim = get_simulation_template(mesh)

sim_ocean = sim.copy()
oceanload = sc.boundary.OceanLoading(side_sets=["y1"])
sim_ocean.add_boundary_conditions(oceanload)
sim_ocean.validate()

sim_coupled = sim.copy()
sim_coupled.set_mesh(mesh_coupled)
boundaries = sc.boundary.HomogeneousDirichlet(side_sets=["y1"])
sim_coupled.add_boundary_conditions(boundaries)
sim_coupled.validate()

salvus.flow.api.run(
    site_name=SALVUS_FLOW_SITE_NAME,
    input_file=sim_ocean,
    ranks=RANKS_PER_JOB,
    output_folder="ocean",
    get_all=True,
    overwrite=True,
)
salvus.flow.api.run(
    site_name=SALVUS_FLOW_SITE_NAME,
    input_file=sim_coupled,
    ranks=RANKS_PER_JOB,
    output_folder="coupled",
    get_all=True,
    overwrite=True,
)

chi = compare_seismograms(
    "coupled/receivers.h5", "ocean/receivers.h5", 40.0, 0.1
)
print("average misfit between ocean load and coupled simulation: ", chi)
assert chi < 0.01

SalvusJob `job_2403150921939263_aeaa66fc70` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

* Downloaded 310.7 KB of results to `ocean`.
* Total run time: 23.61 seconds.
* Pure simulation time: 23.18 seconds.
SalvusJob `job_2403150921513106_4e850c372d` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

* Downloaded 310.7 KB of results to `coupled`.
* Total run time: 19.89 seconds.
* Pure simulation time: 19.67 seconds.
average misfit between ocean load and coupled simulation:  0.002808999064986973

Note that the average misfit here is on the same order of magnitude as in the previous simulation despite the fact that the ocean height is only 0.01 instead of 3.75

Copy

ocean_height = 0.01

mesh, mesh_coupled = generate_meshes(ocean_height)
sim = get_simulation_template(mesh)

sim_ocean = sim.copy()
oceanload = sc.boundary.OceanLoading(side_sets=["y1"])
sim_ocean.add_boundary_conditions(oceanload)
sim_ocean.validate()

salvus.flow.api.run(
    site_name=SALVUS_FLOW_SITE_NAME,
    input_file=sim,
    ranks=RANKS_PER_JOB,
    output_folder="free_surface",
    get_all=True,
    overwrite=True,
)
salvus.flow.api.run(
    site_name=SALVUS_FLOW_SITE_NAME,
    input_file=sim_ocean,
    ranks=RANKS_PER_JOB,
    output_folder="ocean",
    get_all=True,
    overwrite=True,
)
chi = compare_seismograms(
    "free_surface/receivers.h5", "ocean/receivers.h5", 40.0, 0.1
)
print("average misfit between free surface and coupled simulation: ", chi)
assert chi < 0.01

SalvusJob `job_2403150922182505_9000e5da03` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

* Downloaded 310.5 KB of results to `free_surface`.
* Total run time: 19.16 seconds.
* Pure simulation time: 18.82 seconds.
SalvusJob `job_2403150922345866_6377aef77d` running on `local` with 2 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

* Downloaded 310.7 KB of results to `ocean`.
* Total run time: 16.55 seconds.
* Pure simulation time: 16.01 seconds.
average misfit between free surface and coupled simulation:  0.0004279385239773756

2D teleseismic simulation

Receiver Gathers

Site Map Contact Us

Impressum

Academic License Agreement

Credits

PAGE CONTENTS