3 Vertical Refinement

Download Jupyter Notebook

This page is automatically generated from a Jupyter Notebook. You can view and download it here and here.

Note that everything here is fully automatically tested against - and thus guaranteed to work - only for the latest repository versions of all the Salvus packages. So please update if anything works differently on your machine.

Beyond piewise structured grids: vertical refinement in layers with strong thickness variation

This notebooks shows how problems with small time steps in layers with strongly varying thickness can be avoided. One prominent example is the crustal thickness variations.

# initialize notebook
%matplotlib inline
from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (8, 6)
from pymesher.structured_grid_2D import StructuredGrid2D
from pymesher.skeleton import Skeleton

# general input
max_x = 10000.          # domain width in m
fmax = 2.               # maximum frequency in Hz
elements_per_wavelength = 2.
vs = 1000.              # s-wave velocity in m / s
vp = 1700.              # p-wave velocity in m / s

refine_x = 1.2          # extra security factor to make elements smaller before
refine_y = 1.2          # deformation

# compute edgelengths
hmax = vs / fmax / elements_per_wavelength

# some artificial topography model - h(x)
nelem_x = int(np.ceil(max_x / hmax * refine_x))
h0 = .6
h1 = .5
h2 = -.1
h3 = .15

x = np.linspace(0., max_x, nelem_x + 1)
norm_x = x / max_x * 2 * np.pi
h = h0 - h1 * np.cos(norm_x) - h2 * np.cos(2 * norm_x) - h3 * np.sin(3 * norm_x)
h = h * max_x / 2 / np.pi

# number of vertical elements needed for each element in horizontal direction
nelem_y = np.ceil(h / hmax * refine_y).astype('int')

# create box mesh with refinements
sg1 = StructuredGrid2D.rectangle_vertical_refine_doubling(
    nelem_x, nelem_y, min_x=0, max_x=max_x, min_y=0., max_y=np.max(h))

# create box mesh without refinements
sg2 = StructuredGrid2D.rectangle(
    nelem_x, max(nelem_y), min_x=0, max_x=max_x, min_y=0., max_y=np.max(h))

# make unstructured meshes
m1 = sg1.get_unstructured_mesh()
m2 = sg2.get_unstructured_mesh()

# deform the boxmesh according to the topography
for m in [m1, m2]:
    m.add_dem_2D(x, h - np.max(h))
    m.apply_dem()

# compute dt
dt1, dt1_elem = m1.compute_dt(vp)
dt2, dt2_elem = m2.compute_dt(vp)

# plot dt over the mesh
m1.plot(data= dt1_elem, show=False)
m2.plot(data= dt2_elem, show=False)

plt.show()
# compute the cost
cost1 = m1.nelem / dt1
cost2 = m2.nelem / dt2

print('speedup: ', cost2 / cost1)