# 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)