MondaicThis documentation is not for the latest stable Salvus version.
This tutorial demonstrates how to prepare meshes for global seismic wave simulations using the anisotropic adaptive mesh refinement rechnique (aamr), which in this context is also known as smoothiesem [1, 2]. The meshes are specifically build for a particular source location to take advantage of the lateral smoothness of the wavefield and reduce the numerical burden significantly for global and regional scale simulations.
The mesh can be refined in predefined depth, distance and azimuth regions. Additionally, the domain can be restricted to a region of interest as discussed in the separate data adaptive mesh masking tutorial.
In this tutorial we demonstrate the basic concept and most important parameters to build these meshes.
%matplotlib inline
%config Completer.use_jedi = False
from salvus.mesh.simple_mesh import SmoothieSEM
from salvus import namespace as sn
period_global = 100.0
# a quake in Turkey that we will use in this tutorial, original data from IRIS spud:
# http://service.iris.edu/fdsnws/event/1/query?eventid=2847365
# http://ds.iris.edu/spudservice/momenttensor/gcmtid/C201003241411A/quakeml#momenttensor
source = sn.simple_config.source.seismology.SideSetMomentTensorPoint3D(
latitude=38.82,
longitude=40.14,
depth_in_m=4500,
side_set_name="r1",
mrr=5.47e15,
mtt=-4.11e16,
mpp=3.56e16,
mrt=2.26e16,
mrp=-2.25e16,
mtp=1.92e16,
)First, we build some full sphere meshes to highlight the different lateral refinement options.
The simplemost smoothiesem mesh is a full sphere with no lateral refinements and the symmetry axis is aligned with the seismic source location. The tensor order controls boths the accuracy of representing the seismic velocities as well as the geometry of the domain (here: the sphere). On top, the parameter
nlat allows to vary the number of elements in the lateral direction that is used throughout the whole mesh. Compare Figure 9 in [1] to choose values appropriate for your application.sm = SmoothieSEM()
sm.basic.model = "prem_iso_one_crust"
sm.basic.min_period_in_seconds = period_global
sm.basic.elements_per_wavelength = 2.0
sm.advanced.tensor_order = 2
sm.basic.number_of_lateral_elements = 4
sm.source.latitude = source._initial_arguments["latitude"]
sm.source.longitude = source._initial_arguments["longitude"]
sm.create_mesh()
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fc27683c910>
This type of refinement is particularly useful to increase resolution of the surface towards the equator, hence to improve resolution of surface waves as well as surface topography. Here,
theta is the angular distance from the source in the range [0°, 180°]. Multiple such refinements can be combined, but refinement boundaries should not cross to ensure high quality elements: each refinement should be fully contained in all previous refinements.sm.refinement.lateral_refinements = [
{
"theta_min": 60.0,
"theta_max": 150.0,
"r_min": 6000.0,
},
{
"theta_min": 80.0,
"theta_max": 130.0,
"r_min": 6200.0,
},
]
sm.create_mesh()
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fc27512afd0>
This type of refinement may be used to increase resolution in particular areas of interest defined by a distance and azimuth range [-180°, 180°] for [
phi_min, phi_max] relative to the source azimuths measured clockwise from north. This may be a practical alternative to absorbing boundaries and a azimuthally constrained domain.sm.refinement.lateral_refinements = [
{
"theta_min": 40.0,
"theta_max": 110.0,
"r_min": 6000.0,
"phi_min": -45.0,
"phi_max": 45.0,
}
]
sm.source.azimuth = 270.0
sm.create_mesh()
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fc273d18d50>
The lateral refinements can also be contrained to a depth region not at the surface to resolve e.g. small scale structure at a certain dept of interest. This can be achieved with the
r_max parameter of ther lateral refinement object:sm.refinement.lateral_refinements = [
{
"theta_min": 40.0,
"theta_max": 110.0,
"r_min": 3000.0,
"r_max": 4000.0,
"phi_min": -45.0,
"phi_max": 45.0,
}
]
sm.create_mesh()
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fc271d5ed50>
The meshing region can be constrained in depth and distance to create regional or local meshes with higher resolution. This is provided by the following options in the simple mesh object. These parameters can freely be combined with the refinement parameters from above.
sml = SmoothieSEM()
sml.basic.model = "prem_iso_one_crust"
sml.basic.min_period_in_seconds = 10.0
sml.basic.elements_per_wavelength = 2.0
sml.spherical.min_radius = 6000.0
sml.chunk.max_colatitude = 10.0
sml.advanced.tensor_order = 2
sml.basic.number_of_lateral_elements = 4
sml.source.latitude = source._initial_arguments["latitude"]
sml.source.longitude = source._initial_arguments["longitude"]
sml.create_mesh()
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fc2750a5f10>
SmoothieSEM meshes are also capable to include smooth 3D stucture both in terms of the material properties (i.e. seismic velocities) as well as the domain shape (e.g. topography). However, special care needs to be taken to choose appropriate resolution for the mesh as well as the 3D structure and this is currently not automatized. Instead, it is the users responsibility to explicitly choose the resolution of the input data used in the meshing and verify not only it is represented without aliasing on the mesh, but also that the smoothieSEM approximation is still apropriate for the given application e.g. by a convergence test refining in azimuthal direction.
First, download the data files:
import pathlib
import requests
for fname in [
"topography_earth2014_egm2008_lmax_256.nc",
"bathymetry_earth2014_lmax_256.nc",
"moho_topography_crust_1_0_egm2008.nc",
]:
topography_file = pathlib.Path(fname)
if not topography_file.exists():
r = requests.get(
f"https://data.mondaic.com/topography-data/{fname}",
stream=True,
)
assert r.ok
with topography_file.open("wb") as f:
f.write(r.raw.read())Topography, crustal thickness, ocean loading and ellipticity can be included in the meshing using the following options:
sm = SmoothieSEM()
sm.basic.model = "prem_iso_one_crust"
sm.basic.min_period_in_seconds = 100.0
sm.basic.elements_per_wavelength = 2.0
sm.basic.number_of_lateral_elements = 6
sm.advanced.tensor_order = 2
sm.source.latitude = 38.82
sm.source.longitude = 40.14
sm.topography.topography_file = "topography_earth2014_egm2008_lmax_256.nc"
sm.topography.topography_varname = (
"topography_earth2014_egm2008_lmax_256_lmax_32"
)
sm.topography.moho_topography_file = "moho_topography_crust_1_0_egm2008.nc"
sm.topography.moho_topography_varname = (
"moho_topography_crust_1_0_egm2008_lmax_32"
)The topography files can be directly opened with paraview (using the
NetCDF Reader) to vizualize the available datasets. The variable name is composed of the filename appended with the maximum resolution in terms of the spherical harmonic degree and order lmax. For reasonable accuracy, lmax should not be larger than approximately number_of_lateral_elements * tensor_order, so the settings in the example here are very optimistic to keep it small enough.# WGS84 ellipticity value
sm.spherical.ellipticity = 0.0033528106647474805
m = sm.create_mesh()
# for vizualisation compute the radius (note that this in case includes the ellipticity)
m.attach_field("r", (m.points**2).sum(axis=1) ** 0.5)
# choose 'r' in the widget to see the 3D data
m
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fc273a01250>
[1] Driel, Martin van, C Boehm, Lion Krischer, and Michael Afanasiev. 2020. “Accelerating Numerical Wave Propagation Using Wavefield Adapted Meshes. Part I: Forward and Adjoint Modelling.” Geophysical Journal International 221 (3): 1580–90. https://doi.org/10.1093/gji/ggaa058.
[2] Thrastarson, Solvi, Martin van Driel, Lion Krischer, Christian Boehm, Michael Afanasiev, Dirk-Philip van Herwaarden, and Andreas Fichtner. 2020. “Accelerating Numerical Wave Propagation by Wavefield Adapted Meshes. Part II: Full-Waveform Inversion.” Geophysical Journal International 221 (3): 1591–1604. https://doi.org/10.1093/gji/ggaa065.
PAGE CONTENTS