Mondaic DocumentationThe 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 snBelow 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.
vs.plot()[<matplotlib.lines.Line2D at 0x72579f3174d0>]
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 = 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.
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!
sn.layered_meshing.mesh_from_domain(domain=d, model=m, mesh_resolution=mr)
<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x72579e0e9650>
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?
Partitioning
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.
m_p1 = sn.layered_meshing.partition(
d, m, model_predicate=lambda x: x.ds.VS == 0.5, mode="discontinuous"
)m_p1.n_layers2
The partition function spits out another layered model, which we can then
immediately plug into our meshing algorithms again.
sn.layered_meshing.mesh_from_domain(domain=d, model=m_p1, mesh_resolution=mr)
<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7257c2186c90>
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.m_p2 = sn.layered_meshing.partition(
d,
m_p1,
model_predicate=(
sn.layered_meshing.filters.filter_doubling_monotonic_top_down
),
mode="continuous",
)m_p2.n_layers3
sn.layered_meshing.mesh_from_domain(domain=d, model=m_p2, mesh_resolution=mr)
<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7257c2004810>
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.
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,
)
<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7257c2394050>
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.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,
)
<unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7257c2130d90>
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.
PAGE CONTENTS