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)