Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

This documentation is not for the latest stable Salvus 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.

Mesh Quality

Introduction

In this notebook, you build a piecewise structured mesh using a 1D model read from a file and automatic placement of refinements. Play with the input parameters to find out:
  • Which way do you get the best mesh?
  • Does this change for 3D?
  • Does it depend on the frequency?
The automatic placement considers a number of criteria. If any of them is not met, the refinement is pushed further downwards. This is based on the assumption, that velocities increase with depth in most models (which we enforce by making the size function monotonous before calculating element sizes). The criteria are:
  • mimimum resolution in horizontal direction
  • no refinement directly at the surface or the bottom
  • no multiple refinements at the same depth
  • no refinement in very thin layers
Copy
# set up the notebook
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["figure.figsize"] = (10, 8)
The model file, edit as you like with colums in units of km, km/s and kg/m^3.
%%writefile three_layer.bm
NAME         three_layer
UNITS        km
COLUMNS      depth rho vp vs
    0.0   2.6     1.7     1.0
   10.0   2.6     1.7     1.0
   10.0   3.0     3.5     2.2
   15.0   3.0     3.5     2.2
   15.0   3.5     3.8     2.6
  100.0   3.5     3.8     2.6
Writing three_layer.bm
Read the model file and plot the seismic velocities.
from salvus.mesh.models_1D import model

model = model.read("three_layer.bm")
model.plot_vp_vs_profile(depth=True)
The model provides the discontinuities and functionality to compute element sizes according to the resolution criterion. Internally, we work with normalized coordinates, hence the need to scale.
print(
    "discontinuities:",
    ["{:.1f}".format(i) for i in model.discontinuities * model.scale],
)
print(
    "element size:   ",
    [
        "{:.1f}".format(i)
        for i in model.get_edgelengths(
            dominant_period=1.0, elements_per_wavelength=2
        )
        * model.scale
    ],
)
discontinuities: ['0.0', '85000.0', '90000.0', '100000.0']
element size:    ['1300.0', '1100.0', '500.0']
Exercise: Vary hmax_refinement, refinement_style and refinement_top_down to make the best mesh.
Note: Top down approach means minimizing number of elements at the surface at the cost of more elements at the bottom (default). If False, bottom up approach is used, that is minimizing number of elements at the bottom at the cost of more elements at the surface. Top down leads to fewer refinement. Which one is more efficient depends on the velocity model and refinement style.
frequency = 0.1  # maximum frequency in Hz
max_x = 200000.0  # Domain size in horizontal direction in m
hmax_refinement = (
    1.5  # critertion to avoid refinements in thin layers, need to be > 1.0,
)
# default is 1.5, smaller value = more aggressive
refinement_style = "doubling"  # 'doubling' or 'tripling'
refinement_top_down = True  # True or False
ndim = 2  # 2 or 3
from salvus.mesh import mesh_block

if ndim == 2:
    horizontal_boundaries = (np.array([0]), np.array([max_x / model.scale]))
elif ndim == 3:
    horizontal_boundaries = (
        np.array([0, 0]),
        np.array([max_x / model.scale, max_x / model.scale]),
    )

sk = mesh_block.generators.cartesian.mesh_block_collection(
    model.discontinuities,
    model.get_edgelengths(1.0 / frequency),
    hmax_refinement=hmax_refinement,
    horizontal_boundaries=horizontal_boundaries,
    refinement_top_down=refinement_top_down,
    refinement_style=refinement_style,
    ndim=ndim,
)
m = sk.get_unstructured_mesh()
m.find_side_sets(mode="cartesian")
m
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f15f561f190>

a) Equiangular Skewness

A popular quality measure in the community is the equiangular skewness, which is defined as
\begin{align} \text{skew} = \max \left( \frac{\theta_{\max} - \theta_{e}}{180 - \theta_{e}}, \frac{\theta_{e} - \theta_{\min}}{\theta_{e}} \right). \end{align}
Quality meshes must not have skewed elements (skewness <~ 0.75), a single bad element can cause instability in the time extrapolation.
m.plot_quality("equiangular_skewness")
Locate skewed elements visually:
m.attach_field(
    "equiangular_skewness", m.compute_mesh_quality("equiangular_skewness")
)
m