Skip to content

# 1 Basic API

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.

## The most basic mesh to run with Salvus¶

This notebooks shows the most basic mesh with all necessary commands to put material properties and boundaries ready to run with Salvus.

# 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'] = (10, 8)


### 2D¶

Set Basic input parameters for a rectangular domain with homogeneous material properties.

vp = 1700     # p-wave velocity in m/s
vs = 1000     # s-wave velocity in m/s
rho = 2700    # density in kg / m^3

max_x = 10000 # domain size in x direction
max_y = 5000  # domain size in y direction

fmax = 1.     # maximum frequency in Hz
elements_per_wavelength = 2.  # resolution criterion


Compute element size and number of elements needed according to the resolution criterion.

hmax = vs / fmax / elements_per_wavelength
nelem_x = int(max_x / hmax)
nelem_y = int(max_y / hmax)


Build a structured Grid and plot it.

from pymesher.structured_grid_2D import StructuredGrid2D

sg = StructuredGrid2D.rectangle(nelem_x=nelem_x, nelem_y=nelem_y, max_x=max_x, max_y=max_y)
sg.plot(show=True)


Convert the structured grid to an unstructured mesh (trivial in this special case, but uses the same general routines).

m = sg.get_unstructured_mesh()


Attach the material properties, constant here, but using the general array storage on element nodes.

m.attach_field('VP', vp * np.ones(m.npoint))
m.attach_field('VS', vs * np.ones(m.npoint))
m.attach_field('RHO', rho * np.ones(m.npoint))
# this is currently necessary to tell salvus that this is elastic material everywhere. Set to 1 to make it fluid.
m.attach_field('fluid', np.zeros(m.nelem))


Find outer Boundaries of the domain assuming the domain is rectangular.

m.find_side_sets(mode='cartesian')


Write to an exodus file that can be read by paraview and salvus.

m.write_exodus('basic2D.e')


Open with paraview. You propably need to enable the advanced mode clicking on the gear symbol to see the boundaries.

!~/ParaView/bin/paraview basic2D.e


### 3D¶

Very much the same as in 2D, so less comments.

vp = 1700     # p-wave velocity in m/s
vs = 1000     # s-wave velocity in m/s
rho = 2700    # density in kg / m^3

max_x = 10000 # domain size in x direction
max_y = 5000  # domain size in y direction
max_z = 3000  # domain size in y direction

fmax = 1.     # maximum frequency in Hz
elements_per_wavelength = 2.  # resolution criterion

hmax = vs / fmax / elements_per_wavelength
nelem_x = int(max_x / hmax)
nelem_y = int(max_y / hmax)
nelem_z = int(max_z / hmax)


Build a structured Grid. Don't plot it, because in 3D matplotlib is to slow.

from pymesher.structured_grid_3D import StructuredGrid3D

sg = StructuredGrid3D.cube(nelem_x=nelem_x, nelem_y=nelem_y, nelem_z=nelem_z,
max_x=max_x, max_y=max_y, max_z=max_z)

m = sg.get_unstructured_mesh()

m.attach_field('VP', vp * np.ones(m.npoint))
m.attach_field('VS', vs * np.ones(m.npoint))
m.attach_field('RHO', rho * np.ones(m.npoint))

m.find_side_sets(mode='cartesian')

m.write_exodus('basic3D.e')

!~/ParaView/bin/paraview basic3D.e