Mondaic
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.
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 lmBelow 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 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)]
)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)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f154196c290>
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 = lm.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.
lm.mesh_from_domain(domain=d, model=m_p1, mesh_resolution=mr)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f1575fb1710>
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 = lm.partition(
d,
m_p1,
model_predicate=lm.filters.filter_doubling_monotonic_top_down,
mode="continuous",
)m_p2.n_layers3
lm.mesh_from_domain(domain=d, model=m_p2, mesh_resolution=mr)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f1538e532d0>
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.
lm.mesh_from_domain(
domain=d,
model=lm.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.
lm.meshing_protocol.IntralayerVerticalRefine(
refinement_type="doubling"
),
# Enforce a constant element size.
lm.meshing_protocol.IntralayerConstant(),
],
),
mesh_resolution=mr,
)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f1538cd7690>
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.lm.mesh_from_domain(
domain=d,
model=lm.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.
lm.meshing_protocol.IntralayerVerticalRefine(
refinement_type="doubling"
),
# Enforce a constant element size.
lm.meshing_protocol.IntralayerConstant(),
],
# Enforce a constant number of elements in the horizontal direction
# across layer boundaries.
interlayer_coarsening_policy=lm.meshing_protocol.InterlayerConstant(),
),
mesh_resolution=mr,
)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f1538da1ad0>
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