Skip to content

Basic api

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.

The most basic mesh to run with Salvus

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

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

In the mesher, we build unstructured meshes by combining multiple structured grids. For the simplest case here, we build a single structured grid and plot it.

from salvus_mesh.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)

Then, we 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 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')

Open a preview inline, showing all parameters stored on the mesh.

m

3D

The interfaces are 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 salvus_mesh.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

Meshes composed from multiple structured grids

The general strategy in the mesher is to compose unstructured meshes from multiple structured grids. For example, a mesh with a doubling can be thought of as two structured rectilinear grids connected by a doubling grid. The doubling is also considered as a structured grid, with some elements masked out and nodes moved. Nodes that appear multiple times are reduced to single node by lexicographic sorting when converting to an unstructured mesh.

from salvus_mesh.skeleton import Skeleton

2D

Start by creating three structured grids

max_x = 1000.
y = np.array([0., 500., 750., 1000.])
nelem_x = 5
nelem_y = np.array([2, 2])

sg1 = StructuredGrid2D.rectangle(nelem_x, nelem_y[0], min_x=0, max_x=max_x,
                                 min_y=y[0], max_y=y[1])

sg2 = StructuredGrid2D.cartesian_doubling_layer(
    nelem_x, min_x=0, max_x=max_x, min_y=y[1], max_y=y[2])

sg3 = StructuredGrid2D.rectangle(nelem_x * 2, nelem_y[1], min_x=0, max_x=max_x,
                                 min_y=y[2], max_y=y[3])


# plot using different colors
sg1.plot(edgecolor='r')
sg2.plot(edgecolor='b')
sg3.plot(edgecolor='g')

plt.show()

Collect the structured grids in a Skeleton and convert to an unstructured mesh.

sk = Skeleton([sg1, sg2, sg3])
m = sk.get_unstructured_mesh()
m

This is how the doubling layer is built from a rectilinear grid:

sg2mn = StructuredGrid2D.cartesian_doubling_layer(
    nelem_x, min_x=0, max_x=max_x, min_y=y[1], max_y=y[2], move_nodes=False,
    apply_mask=False)
sg2mn.plot(edgecolor='b', scatter=True)

plt.show()

The stuctured grid in 2D is defined by 2 2D arrays for the node locations. The 'elements' are only implicitly defined by the structure of these arrays.

print(sg2mn.x)
print(sg2mn.y)

Elements can be masked out:

sg2mn = StructuredGrid2D.cartesian_doubling_layer(
    nelem_x, min_x=0, max_x=max_x, min_y=y[1], max_y=y[2], move_nodes=False)
sg2mn.plot(edgecolor='b', scatter=True)

plt.show()

And nodes can be moved, making non-rectilinear elements:

sg2mn = StructuredGrid2D.cartesian_doubling_layer(
    nelem_x, min_x=0, max_x=max_x, min_y=y[1], max_y=y[2])
sg2mn.plot(edgecolor='b', scatter=True)

plt.show()

When converting to an unstructured mesh, nodes are made unique by lexicographic sorting. A structured grid is a collection of points and elements defined by connecting these points.

m = sg2mn.get_unstructured_mesh()
print(m.points)
print(m.connectivity)
m.plot(show=False, linewidths=2., scatter=True)
plt.show()

3D

Very similar to 2D for a single layer doubling:

from salvus_mesh.structured_grid_3D import StructuredGrid3D

sg = StructuredGrid3D.cartesian_doubling_layer_single(
    nelem_x=1, nelem_y=1, apply_mask=False, move_nodes=False)
sg.plot(show=False)
plt.axis('off')
plt.show()
sg = StructuredGrid3D.cartesian_doubling_layer_single(
    nelem_x=1, nelem_y=1, move_nodes=False)
sg.plot(show=False)
plt.axis('off')
plt.show()
sg = StructuredGrid3D.cartesian_doubling_layer_single(
    nelem_x=1, nelem_y=1, apply_mask=False)
sg.plot(show=False)
plt.axis('off')
plt.show()

1D Model File Format

The 1D model file format is adopted from AxiSEM and also used as a standard within the NASA Insight mission. It is based on KEY VALUE pairs as follows:

KEY VALUE definitions

The key value format can accomodate 3 different types of values: single values, list of values and tables, all separated by white space. Keys are case sensitive and must not have leading white space. Comments are indicated byt the # character. Continuation lines of tables are indented by at least one blank.

# single value
KEY VALUE

# list of values
KEY VALUE1 VALUE2 VALUE3

# table
KEY HEADER1 HEADER2
  VALUE11 VALUE21
  VALUE12 VALUE22

Available Keys for Salvus

For Salvus, the following keys can be used:

NAME Name of the model. type: string, default: filename without ending.

ANISOTROPIC Whether the model is anisotropic. type: boolean, default: false.

ANELASTIC Whether the model includes an attenuation model. type: boolean, default: false.

REFERENCE_FREQUENCY Frequency at which the seismic velocities are defined in Hz. type: float, default: 1.0.

UNITS Units used for density, seismic velocities and depth/radius, either m for SI units (m, m/s, kg/m 3 ) or km for (km, km/s, g/cm 3 ). type: string, default: m.

COLUMNS Table containing the depth dependent model parameters. type: float.

Depth dependent parameters

Discontinuities in the model (both first and second order) are detected based on repeated radius/depth values. In between these discontinuities, paramteres are assumed to be smooth and interpolated using splines (cubic by default). The table columns in the table COLUMNS have no particular order, the rows can be sorted either from top to bottom or vice versa. The table needs to contain at least these colums:

radius or depth

rho

vp and vs if ANISOTROPIC is false.

vpv, vph, vsv, vsh and eta if ANISOTROPIC is true.

QMU and QKAPPA if ANELASTIC is true

Sample File 1: a local isotropic elastic model with 3 layers

A three layer 1D model where the seismic velocities have gradients in the upper 2 layers and constant below.

NAME         true_model
UNITS        m
COLUMNS      depth rho vp vs
    0.0     2384.4 3500.0 2020.0
    2000.0  2441.9 3850.0 2223.0
    2000.0  2570.1 4725.0 2728.0
    7000.0  2780.8 6475.0 3738.0
    7000.0  2835.5 7000.0 4041.0
    10000.0 2835.5 7000.0 4041.0

Sample File 2: a global PREM model

some lines removed as indicated by [...]:

# Input file for Salvus
NAME         prem_ani
ANELASTIC       T
ANISOTROPIC     T
UNITS        m
COLUMNS       radius      rho      vpv      vsv      qka      qmu      vph      vsh      eta
            6371000.  2600.00  5800.00  3200.00    57827.0      600.0  5800.00  3200.00  1.00000
            6356000.  2600.00  5800.00  3200.00    57827.0      600.0  5800.00  3200.00  1.00000
#          Discontinuity   1, depth:      15.00 km < this just a comment and ignored by the software
            6356000.  2900.00  6800.00  3900.00    57827.0      600.0  6800.00  3900.00  1.00000
            6346600.  2900.00  6800.00  3900.00    57827.0      600.0  6800.00  3900.00  1.00000
#          Discontinuity   2, depth:      24.40 km
            6346600.  3380.75  8190.32  4396.02    57827.0      600.0  8190.32  4611.80  0.90039
            6335480.  3379.54  8182.26  4398.58    57827.0      600.0  8182.26  4601.82  0.90471
            [...]
            5771000.  3975.82 10157.83  5515.93    57827.0      143.0 10157.83  5515.93  1.00000
#          Discontinuity   6, depth:     600.00 km > second order discontinuity
            5771000.  3975.82 10157.76  5516.02    57827.0      143.0 10157.76  5516.02  1.00000
            [...
            3480000.  5566.46 13716.62  7264.65    57827.0      312.0 13716.62  7264.65  1.00000
#          Discontinuity  10, depth:    2891.00 km > fluid by vs=0
            3480000.  9903.44  8064.79     0.00    57827.0        0.0  8064.79     0.00  1.00000
            [...]
            1221500. 12166.33 10355.72     0.00    57827.0        0.0 10355.72     0.00  1.00000
#          Discontinuity  11, depth:    5149.50 km
            1221500. 12763.61 11028.26  3504.31     1327.7       84.6 11028.26  3504.31  1.00000
            [...]
                  0. 13088.50 11262.20  3667.80     1327.7       84.6 11262.20  3667.80  1.00000