Mondaic - Full Waveform Solutions

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.

import numpy as np
import xarray as xr

import salvus.namespace as sn

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.

sn.layered_meshing.mesh_from_domain(
    domain=sn.domain.dim2.BoxDomain(x0=0.0, x1=1.0, y0=0.0, y1=1.0),
    model=sn.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 0x7ddff074d510>

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.

# We cant actually access any material much faster by letting salvus figure
# out which one we want by the arguments supplied
# m0 = sn.material.elastic.Velocity.from_params(rho=1.0, vp=1.0, vs=0.5)
m0 = sn.material.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

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.

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 = sn.material.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(sn.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{m0_lame.to_wavelength_oracle()}")

ConstantParameter(p=0.5)
ConstantParameter(p=0.5)

Indeed they are. Notice that even though we did not specify m0_lame in terms of

V_p

,

V_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 = sn.layered_meshing.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 0x7ddff85029b0>, 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 0x7ddff85029b0>, 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 0x7ddff85029b0>, 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 0x7ddff85029b0>, interpolation_method='linear')

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()

When combined with the interfaces inspected above, our material will be bounded by two planar interfaces at y = 1.0 (i.e, depth = 0) and y = 0.0 (i.e., height = 1).

The final thing we need to set before generating our first mesh is the desired resolution, conveniently encapsulated in the MeshResolution object. Here we can set the desired elements per wavelength, reference frequency, and order of model representation as follows:

mr = sn.MeshResolution(
    reference_frequency=2.0, elements_per_wavelength=1.5, model_order=2
)

where the order of model representation defaults to order 1.

Now we need to bring everything together using the layered mesher's mesh_from_domain function as follows:

mesh = sn.layered_meshing.mesh_from_domain(
    domain=d_2d, model=lm_0, mesh_resolution=mr
)
mesh

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7de0186e9090>

Note that, as we did in the "primer", wrapping the input material into a LayeredModel is optional for models with only 1 layer and no explicitly defined interfaces. Indeed, executing the above cell with the m0 model instance directly achieves the same result.

sn.layered_meshing.mesh_from_domain(domain=d_2d, model=m0, mesh_resolution=mr)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd7b8a50>

The separation of domain definitions, mesh resolution, materials, and interfaces makes it now very easy to switch between dimensions. Indeed, all that's required to create a 3-D mesh is to just switch the definition of our domain.

d_3d = sn.domain.dim3.BoxDomain(x0=0, x1=1, y0=0, y1=1, z0=0, z1=1)
sn.layered_meshing.mesh_from_domain(domain=d_3d, model=m0, mesh_resolution=mr)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd83bf10>

We can of course solve the wave equation in homogeneous elastic models analytically -- let's move on to something a bit more advanced by defining a linear velocity gradient with depth.

vp_grad = xr.DataArray(
    np.linspace(1.0, 2.0, 11), [("y", np.linspace(1.0, 0.0, 11))]
)
m1 = sn.material.from_params(rho=1.0, vp=vp_grad, vs=0.5 * vp_grad)

Note here that we are mixing variable and homogeneous parameters. This is no problem -- all parameters will be automatically expanded onto a consistent grid before any meshing takes place. Let's re-use our domain and mesh resolution definitions from above to skip directly to the mesh generation stage.

sn.layered_meshing.mesh_from_domain(domain=d_2d, model=m1, mesh_resolution=mr)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd8caf90>

To switch the above model to 3-D requires one extra step not present in the first example. Consider what happens when we just switch the domain.

sn.layered_meshing.mesh_from_domain(domain=d_3d, model=m1, mesh_resolution=mr)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd6f5cd0>

At first glance this may look like what we want, but notice that the gradient is horizontal, not in depth! This is because we explicitly requested the gradient to be in the "y" coordinate above. "y" is indeed the vertical coordinate in 2-D domains, and to get this in 3-D we simply need to rename the coordinate.

vp_grad_z = vp_grad.rename({"y": "z"})

m1_z = sn.material.from_params(rho=1.0, vp=vp_grad_z, vs=0.5 * vp_grad_z)

sn.layered_meshing.mesh_from_domain(
    domain=d_3d, model=m1_z, mesh_resolution=mr
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd596c10>

Much better :). As another note here though, you can use models with mixed dimensions. Consider the following example. What would you expect here? Does the mesh respect those expectations?

sn.layered_meshing.mesh_from_domain(
    domain=d_3d,
    model=sn.material.from_params(rho=1.0, vp=vp_grad, vs=0.5 * vp_grad_z),
    mesh_resolution=mr,
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd6eb910>

Finally, we can avoid the above issue altogether by specifying our vertical coordinate as "v" (for "vertical"), which will map between y and z depending on the domain's dimension.

vp_grad_v = xr.DataArray(
    np.linspace(1.0, 2.0, 11), [("v", np.linspace(1.0, 0.0, 11))]
)
m1_v = sn.material.from_params(rho=1.0, vp=vp_grad_v, vs=0.5 * vp_grad_v)

# 2d
sn.layered_meshing.mesh_from_domain(
    domain=d_2d, model=m1_v, mesh_resolution=mr
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd71fe90>

# 3d
sn.layered_meshing.mesh_from_domain(
    domain=d_3d, model=m1_v, mesh_resolution=mr
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd4c2c90>

Let's make things a bit more interesting by creating a mesh with two layers. First, let's define two materials.

m2_l0 = sn.material.from_params(rho=0.5, vp=0.5)
m2_l1 = m1_v

Now we need to do something new: define an interface which separates the two materials. The simplest interface is a "Hyperplane", which can be defined at a constant vertical coordinate as follows.

m2_i0 = sn.layered_meshing.interface.Hyperplane.at(0.5)

Now we have a nontrivial collection of materials. In order to tell Salvus that we want to combine these definitions into a single layered model, we need to initialize our LayeredModel object with a list of strata, defined in a top-down order.

lm_2 = sn.layered_meshing.LayeredModel([m2_l0, m2_i0, m2_l1])

Strata must be alternating instances of materials and layers, always specified in a top-down order. With this new layered model specified, we can go ahead and generate a new mesh.

sn.layered_meshing.mesh_from_domain(
    domain=d_2d, model=lm_2, mesh_resolution=mr
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd581690>

Go ahead and inspect the mesh. Are your expectations met? Why all of a sudden is the mesh irregular? Moving to 3-D is again trivial.

sn.layered_meshing.mesh_from_domain(
    domain=d_3d, model=lm_2, mesh_resolution=mr
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7de01874fc90>

For one final example we'll make things a bit more interesting by specifying more interesting interface characteristics. Let's define 2 deformed interfaces as the following.

x = np.linspace(0, 1, 51)

m3_i0 = sn.layered_meshing.interface.Curve.from_points(
    x,
    np.sin(2 * np.pi * x) * 0.1 - 0.1,
    reference_elevation=sn.layered_meshing.interface.Depth(0.0),
    axis="x",
)

m3_i1 = sn.layered_meshing.interface.Curve.from_points(
    x,
    np.sin(np.pi * x) * 0.2,
    reference_elevation=sn.layered_meshing.interface.Depth(0.5),
    axis="x",
)

Note here that we're now using a relative coordinate for our interfaces' reference elevation. We're also specifying a top bounding interface, and no longer want Salvus to add one implicitly. With this in mind, can you see why we're subtracting

0.1

for the elevation of the top interface?

Our layered model definition now looks like:

lm_3 = sn.layered_meshing.LayeredModel([m3_i0, m2_l0, m3_i1, m2_l1])

and creating a mesh proceeds almost as before.

sn.layered_meshing.mesh_from_domain(
    domain=d_2d, model=lm_3, mesh_resolution=mr
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdcf25010>

Finally, moving to 3-D:

lm_3_z = sn.layered_meshing.LayeredModel([m3_i0, m2_l0, m3_i1, m2_l1])
sn.layered_meshing.mesh_from_domain(
    domain=d_3d, model=lm_3_z, mesh_resolution=mr
)

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7ddfdd275610>

GETTING STARTED

© Mondaic AG (2024)

Site Map Contact Us

Academic License Agreement

PAGE CONTENTS