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.

Layered Meshing 01: Basics

The "layered meshing" interface in Salvus is a new API that reimagines how certain model types are meshed. In this tutorial we will explore some of the basics and give an overview of the different meshing "policies" one can choose from. Let's get started!
As always, the first step is to import the Salvus Python environment. Here we'll add an extra import to the canonical Salvus namespace to bring the layered meshing API into scope.
Copy
import numpy as np
import xarray as xr

import salvus.namespace as sn
import salvus.mesh.layered_meshing as lm

Primer

The layered mesher simplifies the process of generating multi-layered meshes, with different physical material models allowed in each layer, and with each material parameterization accepting distinct n-dimensional values per parameter. Before we explore this complexity, however, the interface also allows for the generation of simple meshes in a minimal number of statements. For example, a homogeneous elastic unit square can be meshed with the following statement.
lm.mesh_from_domain(
    domain=sn.domain.dim2.BoxDomain(x0=0.0, x1=1.0, y0=0.0, y1=1.0),
    model=lm.material.elastic.Velocity.from_params(rho=1.0, vp=1.0, vs=0.5),
    mesh_resolution=sn.MeshResolution(reference_frequency=10.0),
)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f7a7e94c4d0>
Now, let's dig into the details.
In this tutorial we'll consider different idealized models defined on a unit square. First we'll make things super simple by just continuing with our definition of a homogeneous elastic unit square. To do this we'll need to select a "material" from the layered meshing toolbox, along with a parameterization.
m0 = lm.material.elastic.Velocity.from_params(rho=1.0, vp=1.0, vs=0.5)
We can explore this material a bit to see what it consists of. Important to note is that all materials are internally represented by a xarray Dataset, which can be queried by the materials .ds variable.
m0.ds
<xarray.Dataset> Size: 24B
Dimensions:  ()
Data variables:
    RHO      float64 8B 1.0
    VP       float64 8B 1.0
    VS       float64 8B 0.5
<xarray.Dataset> Size: 24B
Dimensions:  ()
Data variables:
    RHO      float64 8B 1.0
    VP       float64 8B 1.0
    VS       float64 8B 0.5
As you may expect, we simply see an xarray representation of the material parameters we specified. You might notice that no coordinates are specified here. This is no problem, as the coordinates will be "realized" when the model is eventually used in the meshing process. Realization in this context refers to ensuring that all coordinates present in a layered model are consistent with the spatial dimension and domain they are being used within. As no coordinates are present here, the realization stage will thus ensure that the model values are propagated to all locations in the host domain. In the future we will see how realization can be exploited to define models on only a subset of a domain's coordinate axes, and additionally define them relative to either the domain bounds, or to the bounds of a host layer.

Different parameterizations

The parameterization of our model in terms of velocities is not the only way we can specify an isotropic elastic model. We could have equivalently specified our parameters in the form of Lamé Parameters, which for an elastic model consists of λ\lambda, μ\mu and ρ\rho. For instance:
m0_lame = lm.material.elastic.isotropic.LameParameters.from_params(
    lam=0.5, mu=0.25, rho=1.0
)
results in the same material as we just specified above. Indeed, all materials can be transformed from one parameterization to another easily.
print(lm.material.elastic.Velocity.from_material(m0_lame))
Velocity(RHO=ConstantParameter(p=1.0), VP=ConstantParameter(p=1.0), VS=ConstantParameter(p=0.5))
This way of switching parameterizations also means that the spatially variable minumum wavelength for each parameterization can be automatically computed by Salvus. To this end, each material also contains a method named .to_wavelength_oracle() that will automatically compute the critical parameter required to determine the required mesh size. For isotropic elastic models this is simply the model's shear wave velocity, but for less symmetric materials the oracle may instead be a combination of the material's parameters.
Let's verify that the wavelength oracle for the two parameterizations we've specified is the same as expected.
print(
    f"{m0.to_wavelength_oracle(n_dim=3)}\n{m0_lame.to_wavelength_oracle(n_dim=3)}"
)
ConstantParameter(p=0.5)
ConstantParameter(p=0.5)
Indeed they are. Notice that even though we did not specify m0_lame in terms of VpV_p, VsV_s, and ρ\rho, Salvus automatically computed it for us. The presence of this oracle for each parameterization means that the tedious process of manually computing the minimum wavelength is a thing of the past -- Salvus will take care of this for you.
The model we've created above, while homogeneous, can also be seen as a layered model with a single layer. The utility of this will become apparent later in the tutorial; for now we just need to to wrap our material in a LayeredModel to tell Salvus that we are ready to proceed to the next stage of mesh generation.
lm_0 = lm.LayeredModel(m0)
The LayeredModel object has some interesting properties. First off, we can query the models it contains and see our material stored therein. Note that the models are presented as a list with one element; again, the utility of this will be apparent later on.
lm_0.models
[Velocity(RHO=ConstantParameter(p=1.0), VP=ConstantParameter(p=1.0), VS=ConstantParameter(p=0.5))]
We can also query the .interfaces property. Here, we see an empty list.
lm_0.interfaces
[]
Each layer in a layered model is defined as a material bounded by two interfaces. When we created our layered model above, however, we did not specify any interfaces. If not specified explicitly, Salvus will automatically add bounding hyperplanes to the top and bottom and bottom of the model. To see this in action, you can call the .complete() method on the layered model and inspect the interfaces added.
print("\n\n".join(str(s) for s in lm_0.complete().interfaces))
Hyperplane(da=<xarray.DataArray ()> Size: 8B
array(0.)
Attributes:
    reference_elevation:  Depth(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7f7a86c18450>, interpolation_method='linear')

Hyperplane(da=<xarray.DataArray ()> Size: 8B
array(0.)
Attributes:
    reference_elevation:  Height(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7f7a86c18450>, interpolation_method='linear')
Finally, you can inspect the complete layered model, including bounding interfaces and internal materials, by inspecting the model's strata.
print("\n\n".join(str(s) for s in lm_0.complete().strata))
Hyperplane(da=<xarray.DataArray ()> Size: 8B
array(0.)
Attributes:
    reference_elevation:  Depth(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7f7a86c18450>, interpolation_method='linear')

Velocity(RHO=ConstantParameter(p=1.0), VP=ConstantParameter(p=1.0), VS=ConstantParameter(p=0.5))

Hyperplane(da=<xarray.DataArray ()> Size: 8B
array(0.)
Attributes:
    reference_elevation:  Height(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7f7a86c18450>, interpolation_method='linear')

Reference coordinates

The absolute locations of the boundary interfaces we inspected above are not yet determined -- in fact we haven't yet specified any coordinates yet at all! The implicitly generated interfaces are therefore specified in terms of relative coordinates. Let's inspect the types of both interfaces' reference coordinates to see what is going on.
i0, i1 = lm_0.complete().interfaces

print(i0.da.reference_elevation)
print(i1.da.reference_elevation)
Depth(value=0.0)
Height(value=0.0)
Here we see that the top interface (index 0) is referenced to a depth of 0, and the bottom interface (index 1) is referenced to a height of 0. In this way we are able to defer the computation of the interfaces' absolute locations to later, for instance when a domain is defined. These relative coordinates can also become quite handy when defining generic internal boundaries with reference to the domain bounds, and we'll learn how to specify them later.
Of course we need to specify absolute coordinates at some point, and the simplest way to do this is through Salvus' Domain object. For now, let's just use a simple 2-D Box domain.
d_2d = sn.domain.dim2.BoxDomain(x0=0, x1=1, y0=0, y1=1)
d_2d.plot()