Mondaic - Full Waveform Solutions

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".

Copy

import numpy as np
import xarray as xr

import salvus.namespace as sn

Below we create an input model with vp, vs, and rho varying with depth.

Copy

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.

Copy

vs.plot()

[<matplotlib.lines.Line2D at 0x7a7780c13290>]

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.

Copy

xc = np.linspace(0, 1, 101)
i0 = sn.layered_meshing.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.

Copy

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 = sn.layered_meshing.LayeredModel(
    [i0, sn.material.from_params(rho=rho, vp=vp, vs=vs)]
)

One can already throw this model into the layered mesh and see what comes out. Let's try that!

Copy

sn.layered_meshing.mesh_from_domain(domain=d, model=m, mesh_resolution=mr)

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

Hmm: while we've got a mesh that respects the minimum wavelength of the model, it doesn't look so ideal. The discontinuity near the surface is not represented well in the mesh, and we're not taking advantage of the increase in wavespeeds at depth to allow for larger elements sizes. What can we do?

The layered mesher has a set of model partitioning algorithms that can be used to partition a continuous input model into a layered one. The partition function takes as input a domain, a layered model, and at minimum a "predicate" which informs the algorithm how to split the input layered model. Below we pass a predicate that marks locations in the model where Vs is equal to 0.5; the model will be split into distinct regions for every discrete value returned by the predicate. Since we know a-priori that there is a true discontinuity that we're looking for, we specify that the partitioning should be discontinuous.

Copy

m_p1 = sn.layered_meshing.partition(
    d, m, model_predicate=lambda x: x.ds.VS == 0.5, mode="discontinuous"
)

Copy

m_p1.n_layers

The partition function spits out another layered model, which we can then immediately plug into our meshing algorithms again.

Copy

sn.layered_meshing.mesh_from_domain(domain=d, model=m_p1, mesh_resolution=mr)

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

Ah, this looks a lot better! Our discontinuity is now well represented in the mesh.

There is more we can do however. A common thing to do in geological SEM simulations is to "double" the mesh in regions where the velocity allows for larger wavelengths. To accomodate this, Salvus has several partitioning algorithms built in that can be passed as predicates. One of these is "filter_doubling_monotonic_top_down". This will identify regions in the model where doubling can safely occur. We choose a "continuous" partitioning mode here as there is no true discontinuity we're looking for, simply a logical one that will allow the mesh to add a doubling layer.

Copy

m_p2 = sn.layered_meshing.partition(
    d,
    m_p1,
    model_predicate=(
        sn.layered_meshing.filters.filter_doubling_monotonic_top_down
    ),
    mode="continuous",
)

Copy

m_p2.n_layers

Copy

sn.layered_meshing.mesh_from_domain(domain=d, model=m_p2, mesh_resolution=mr)

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

Looks good! We're now respecting the increase in velocity with depth and reducing the total element count accordingly.

You may notice that the elements in the topmost layer are squeezed between the surface and the first discontinuity. We can request that Salvus produce a more uniform element size by specifying an Intralayer Coarsening Policy. These policies dictate what types of modifications Salvus is allowed to make to the mesh within each layer. One policy is allowed per layer, and the last policy will be applied to all subsequent layers. Let's try out the policy on our mesh below, a more thorough look at the different options will be explored in a future tutorial.

Copy

sn.layered_meshing.mesh_from_domain(
    domain=d,
    model=sn.layered_meshing.MeshingProtocol(
        m_p2,
        intralayer_coarsening_policy=[
            # Allow for a variable number of elements in the vertical
            # direction. doubling: Place refinements at the top of the layer.
            sn.layered_meshing.meshing_protocol.IntralayerVerticalRefine(
                refinement_type="doubling"
            ),
            # Enforce a constant element size.
            sn.layered_meshing.meshing_protocol.IntralayerConstant(),
        ],
    ),
    mesh_resolution=mr,
)

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

Nice, things are looking even better! We've gotten rid of the squeezed elements.

In contrast with intralayer coarsening policies, in- ter -layer policies dictate how the mesh can be coarsened between layers. As we saw above, mesh doubling layers have been automatically inserted between each discontinuity as the velocities support this strategy. In fact, what happened is that the default interlayer policy is set to InterlayerDoubling(). We can change this by explicitly setting a different policy, i.e. InterlayerConstant() as below. One policy is allowed per internal interface, and the last policy will be applied to all subsequent interfaces. Let's quickly inspect the effect of changing the policy to InterlayerConstant() below, a more thorough look at the different options will be explored in a future tutorial.

Copy

sn.layered_meshing.mesh_from_domain(
    domain=d,
    model=sn.layered_meshing.MeshingProtocol(
        m_p2,
        intralayer_coarsening_policy=[
            # Allow for a variable number of elements in the vertical
            # direction. doubling: Place refinements at the top of the layer.
            sn.layered_meshing.meshing_protocol.IntralayerVerticalRefine(
                refinement_type="doubling"
            ),
            # Enforce a constant element size.
            sn.layered_meshing.meshing_protocol.IntralayerConstant(),
        ],
        # Enforce a constant number of elements in the horizontal direction
        # across layer boundaries.
        interlayer_coarsening_policy=(
            sn.layered_meshing.meshing_protocol.InterlayerConstant()
        ),
    ),
    mesh_resolution=mr,
)

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

It's always useful to understand the defaults as there are situations where they may not produce the optimal mesh. The refinements added reduce the number of elements significantly, but the cost of this is a hit to the time step due to the deformations. As such, turning off doubling may be useful in some situations.

Part 1 - Basics

Part 3 - Meshing complex geophysical domains

Site Map Contact Us

Impressum

Academic License Agreement

Credits

PAGE CONTENTS