Version:

This documentation is not for the latest stable Salvus version.

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.

Attenuation

This micro-tutorial demonstrates how to run viscoelastic simulations with SalvusProject. It shows the concept in 2-D but it works just the same in 3-D by just adding an extra dimension.
Copy
import matplotlib.pyplot as plt
import numpy as np
import os
import pathlib
import xarray as xr

import salvus.namespace as sn

# Run simulations on this site.
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
--> Server: 'https://data.mondaic.com/license_server/licensing_server', User: '_MONDAIC_INTERNAL_', Group: '_MONDAIC_INTERNAL_'.
--> Negotiating 1 license instance(s) for 'SalvusMesh' [license version 1.0.0] for 1 seconds ...
--> Success! [Total duration: 3.13 seconds]

Parameter Model

There are different ways to define meshes and models in SalvusProject. Here we'll define the multi-parameter model as an xarray data set which will later be interpolated onto our numerical meshes.
# Grid points of the velocity model.
# NOT the grid point of the simulation!
nx, ny = 200, 200

# Create a coordinate grid.
x = np.linspace(-100.0, +100.0, nx)
y = np.linspace(-100.0, +100.0, nx)
xx, yy = np.meshgrid(x, y, indexing="ij")

# Simplistic 5 parameter model with a depth gradient and
# some sinusoidal perturbation.
# This is just an example but this would be the place to add
# your own velocity model.
vp = 1500.0 * np.ones_like(xx) - yy * 0.1 + np.sin(xx / 8 + 0.5) * 2
vs = 1000.0 * np.ones_like(xx) - yy * 0.1 + np.sin(xx / 12 - 0.3) * 2
rho = 980.0 * np.ones_like(xx) - yy * 0.1 + np.sin(xx / 15 + 1) * 2

# Assemble into a dataset without Q.
ds_without_q = xr.Dataset(
    data_vars={
        "vp": (["x", "y"], vp),
        "vs": (["x", "y"], vs),
        "rho": (["x", "y"], rho),
    },
    coords={"x": x, "y": y},
)


# Do the same again but this time with a Q model.
qkappa = 2000.0 * np.ones_like(xx) - yy * 0.1 + np.sin(xx / 6) * 2
qmu = 75.0 * np.ones_like(xx) - yy * 0.1 + np.sin(xx / 10 - 1.2) * 2

ds_with_q = xr.Dataset(
    data_vars={
        "vp": (["x", "y"], vp),
        "vs": (["x", "y"], vs),
        "rho": (["x", "y"], rho),
        "qkappa": (["x", "y"], qkappa),
        "qmu": (["x", "y"], qmu),
    },
    coords={"x": x, "y": y},
)

# Have a look at the Qmu model.
plt.figure(figsize=(10, 8))
ds_with_q.qmu.T.plot()
plt.show()
# Either load an existing project.
if pathlib.Path("project").exists():
    print("Opening existing project.")
    p = sn.Project(path="project")
# Or create a new one.
else:
    print("Creating new project.")
    p = sn.Project.from_domain(
        path="project",
        # Should match the volume model.
        domain=sn.domain.dim2.BoxDomain(x0=-100, x1=100, y0=-100, y1=100),
    )
Creating new project.
Here we'll create a single explosive source in the center with receivers in a ring around it.
src = sn.simple_config.source.cartesian.MomentTensorPoint2D(
    x=0, y=0, mxx=1.0, myy=1.0, mxy=0.0
)

recs = sn.simple_config.receiver.cartesian.collections.RingPoint2D(
    x=0, y=0, radius=75.0, count=20, fields=["displacement"]
)

p += sn.Event(event_name="event_a", sources=src, receivers=recs)

p.viz.nb.domain()
Now we'll add both velocity models.
p += sn.model.volume.cartesian.GenericModel(
    name="model_without_Q", data=ds_without_q
)

p += sn.model.volume.cartesian.GenericModel(
    name="model_with_Q", data=ds_with_q
)
# A few wave-lengths are needed to see an attenuating effect.
center_frequency = 100.0

# Set-up two simulations: One with attenuation, one without.
for attenuation in [True, False]:
    if attenuation:
        simulation_name = "with_attenuation"
        mc = sn.ModelConfiguration(
            background_model=None,
            volume_models="model_with_Q",
            # Attenuation results in slightly frequency dependent wave speed.
            # Thus one must specify at what frequency they are specified.
            # They are often given at 1 Hz.
            linear_solids=sn.LinearSolids(reference_frequency=1.0),
        )
    else:
        simulation_name = "without_attenuation"
        mc = sn.ModelConfiguration(
            background_model=None, volume_models="model_without_Q"
        )

    p += sn.SimulationConfiguration(
        name=simulation_name,
        max_frequency_in_hertz=center_frequency * 2.0,
        elements_per_wavelength=2.0,
        model_configuration=mc,
        event_configuration=sn.EventConfiguration(
            wavelet=sn.simple_config.stf.Ricker(
                center_frequency=center_frequency
            ),
            waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
                end_time_in_seconds=0.07,
                # This is to tell the simulation to actually
                # simulate with attenuation or not.
                attenuation=attenuation,
            ),
        ),
    )

p.viz.nb.simulation_setup("with_attenuation")
[2021-10-19 09:59:58,664] INFO: Creating mesh. Hang on.