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.

Layered Meshing 02: Basics of partitioning and policies

The layered meshing interface introduced in the previous tutorial works well when your model is already split into distinct layers. In many cases, however, all you have to work with is a model file, usually discretized onto a rectilinear grid. In these cases, it may not be trivial to manually split out layers to feed into the layered meshing algorithm. In this tutorial we'll walk through a simple example that shows how "partitioning" models works in Salvus, as well as touching on some more advanced meshing techniques, or "policies".
import numpy as np
import xarray as xr

import salvus.namespace as sn
import salvus.mesh.layered_meshing as lm
Below we create an input model with vp, vs, and rho varying with depth.
vp = xr.DataArray(
    np.r_[np.full(11, 1.0), np.linspace(2.0, 7.5, 51)],
    [("y", np.linspace(1.0, 0.0, 62))],

vs = xr.DataArray(
    np.r_[np.full(11, 0.5), np.linspace(1.0, 3.5, 51)],
    [("y", np.linspace(1.0, 0.0, 62))],

rho = xr.DataArray(
    np.r_[np.full(11, 1.0), np.linspace(2.0, 7.5, 51)],
    [("y", np.linspace(1.0, 0.0, 62))],
Let's plot the vs array to get a feel for what we're working with.
[<matplotlib.lines.Line2D at 0x7f15577a3310>]
As we can see there is a discontinuity near the top of the model, and then a steep linear gradient with depth. Let's also define some slight deformation on the top surface to make things more interesting.
xc = np.linspace(0, 1, 101)
i0 = lm.interface.Curve.from_points(
    xc, np.sin(2 * np.pi * xc) * 0.05 - 0.05, reference_elevation=1.0, axis="x"
Finally, let's define a domain, a mesh resolution object, and bring everything together as a layered model.
d = sn.domain.dim2.BoxDomain(x0=0, x1=1, y0=0, y1=1)
mr = sn.MeshResolution(reference_frequency=10.0, elements_per_wavelength=2.0)
m = lm.LayeredModel(
    [i0, lm.material.elastic.Velocity.from_params(rho=rho, vp=vp, vs=vs)]


A naive approach

One can already throw this model into the layered mesh and see what comes out. Let's try that!
lm.mesh_from_domain(domain=d, model=m, mesh_resolution=mr)