# Topography and refinements

Download Jupyter Notebook

This page is automatically generated from a Jupyter Notebook.

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

## Beyond piecewise 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.

Salvus provides structured grid classes in 2D and 3D, that allow to vary the number of elements in one dimension, creating anisotropic element sizes. When combined with a corresponding topography model, the elements end up with isotropic sizes again and avoid the small time step that would result from a regular grid.

The workflow is as follows:

- create a topography model
- compute the number of elements needed in vertical direction according to the resolution criterion
- create a structured grid with anisotropic refinement
- convert to an unstructured mesh and deform the elements according to the topography

# initialize notebook %matplotlib inline import matplotlib.pyplot as plt import numpy as np plt.style.use('ggplot') plt.rcParams['figure.figsize'] = (12, 8)

from salvus_mesh.structured_grid_2D import StructuredGrid2D from salvus_mesh.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 # create 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)) # make unstructured mesh m1 = sg1.get_unstructured_mesh() m1.find_side_sets('cartesian') m1.plot(linewidths=1.)

At this point we have a somewhat unusal looking mesh on a cartesian domain and varying number of elements in vertical direction. Once we deform it according to the topography, it looks much more reasonable

m1.add_dem_2D(x, h - np.max(h)) m1.apply_dem() m1.plot(linewidths=1.)

For reference we also create a 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)) m2 = sg2.get_unstructured_mesh() m2.find_side_sets('cartesian') m2.plot(linewidths=1.)

which has very thin elements after applying topography:

m2.add_dem_2D(x, h - np.max(h)) m2.apply_dem() m2.plot(linewidths=1.)

These result in small time steps due to the Courant criterion discussed in exercise 2. We again use the accurate algorithm to estimate \Delta t, which accounts for the shape of the elements:

dt1, dt1_elem = m1.compute_dt(vp, fast=False) dt2, dt2_elem = m2.compute_dt(vp, fast=False) print('dt1: {:10.4f}'.format(dt1)) print('dt2: {:10.4f}'.format(dt2)) print('dt1 / dt2: {:10.4f}'.format(dt1 / dt2)) # plot dt over the mesh m1.attach_field('dt', dt1_elem) m2.attach_field('dt', dt2_elem)

dt1: 0.0199 dt2: 0.0061 dt1 / dt2: 3.2646

If we include the number of elements, which is also reduced with the anistropic refinement, we can estimate the resulting speed up of the simulation:

cost1 = m1.nelem / dt1 cost2 = m2.nelem / dt2 print('nelem1: ', m1.nelem) print('nelem2: ', m2.nelem) print('speedup: {:.2f}'.format(cost2 / cost1))

nelem1: 273 nelem2: 480 speedup: 5.74

```
m1
```

VBox(children=(HTML(value='\n <p><i>This is only a preview. Use ParaView to actually investigate\n … <salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x10446f400>

```
m2
```

VBox(children=(HTML(value='\n <p><i>This is only a preview. Use ParaView to actually investigate\n … <salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x102ac7da0>

### Outlook: Time step optimization¶

The mesh can be further improved by smoothing with a criterion that directly optimizes the time step. Note that the method is work in progress and currently does not ensure that the resolution criterion remains fulfilled and hence should be used with caution.

from salvus_mesh.optimize_dt import optimize_dt # collect side sets where nodes should not be moved side_sets = [m1.get_side_set(ssn) for ssn in m1.side_sets.keys()] # actual optimization m3 = optimize_dt(m1, [side_sets] * 2, maxiter=50) # compute new dt and speedup compared to unsmoothed mesh dt3, dt3_elem = m3.compute_dt(vp, fast=False) cost3 = m3.nelem / dt3 print('dt3: {:10.4f}'.format(dt3)) print('speedup by smoothing: {:10.4f}'.format(cost1 / cost3)) print('total speedup: {:10.4f}'.format(cost2 / cost3)) # plot dt over the mesh m3.attach_field('dt', dt3_elem) m3

dt3: 0.0324 speedup by smoothing: 1.6269 total speedup: 9.3382 VBox(children=(HTML(value='\n <p><i>This is only a preview. Use ParaView to actually investigate\n … <salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x11742a710>