MondaicAbout UsProductsServicesDocs
Version:

This documentation is not for the latest stable Salvus version.

This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

Model Interpolation

In this notebook we'll explore the concept of "model interpolation order". The spectral-element method works by approximating the dynamic field variables as Lagrange polynomials with time-varying coefficients. These polynomials are defined using the Gauss-Lobatto-Legendre (GLL) collocation points defined within each element. For a polynomial interpolation order of degree 4 (the default in Salvus), there will be 5 GLL points along each edge of a given element, which results in 25 and 125 points each for 2-D and 3-D elements, respectively. In general, a dd-dimensional spectral element equipped with order nn interpolating polynomials will have (n+1)d(n+1)^d GLL points in its interior.

While the GLL basis is the canonical method for representing the wavefield, we can also use it to represent the material parameters and spatial geometry. This tutorial shows how one can use different order (1, 2, and 4) interpolating polynomials to represent the model in Salvus. Specifically, we'll:

  • Discretize a linearly-varying 2-D model onto 1st and 4th order polynomial bases
  • Compare the wavefield solutions generated in each of these cases.
  • Do the same in 3-D.
Copy
%matplotlib inline
%config Completer.use_jedi = False
import os
from functools import partial
from pathlib import Path

import h5py
import matplotlib.pyplot as plt
import numpy as np

import salvus.mesh.structured_grid_2D as sg2d
import salvus.mesh.structured_grid_3D as sg3d
import salvus.mesh.unstructured_mesh as um
import salvus.toolbox.toolbox as st
from salvus.flow import api
from salvus.flow import simple_config as sc

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
ranks = 2

Setting up the models

First we set up a simple function to return basic meshes for testing in both 2 and 3 dimensions. The size of the mesh is small as these examples will only run for a trivial duration, and are just meant to be used to explore the possible parameterizations. We'll also set up a constant density throughout the domain.

n = 25
x = 1.0
rho = 1000.0
mesh = sg2d.StructuredGrid2D.rectangle(
    nelem_x=n, nelem_y=n, max_x=x, max_y=x
).get_unstructured_mesh()
mesh.find_side_sets()

Now we copy the mesh skeleton we created above into two separate meshes, each of which we'll set to a different model interpolation order. For brevity we'll use here orders 1 and 4 (order 2 is also supported and the workflow is identical).

mesh_order_1 = mesh.copy()
mesh_order_1.change_tensor_order(1)

mesh_order_4 = mesh.copy()
mesh_order_4.change_tensor_order(4)

Now we'll copy the coordinate arrays for use with the interpolation function we will write below. Note that we randomly choose which axis to interpolate over.

axis = np.random.choice([0, 1])
coords_order_1 = mesh_order_1.get_element_nodes()[:, :, axis].flatten()
coords_order_4 = mesh_order_4.get_element_nodes()[:, :, axis].flatten()

The function below will linearly interpolate a range of model parameters over the coordinate range of a mesh.

def li(x: np.ndarray, min_val, max_val) -> np.ndarray:
    """Linearly interpolate a parameter over a coordinate array.

    Parameters
    ----------
    x : np.ndarray
        Flattened coordinate array (either x or y).
    min_val : float
        Value of the parameter at the minimum coordinate.
    max_val : float
        Value of the parameter at the maximum coordinate.

    Returns
    -------
    np.ndarray
        An array of values linearly interpolated over the coordinate range.

    """
    m = (max_val - min_val) / (np.max(x) - np.min(x))
    return m * x + min_val

We'll use the linear M0M0 and M1M1 to parameterize our acoustic domain. For more information on different supported parameterizations in Salvus, please refer to the relevant tutorial.

vp0, vp1 = 1000, 2000
m10 = m11 = 1 / rho
m00 = m10 / vp0 ** 2
m01 = m11 / vp1 ** 2

Now we interpolate and attach the desired parameter to our 1st order mesh...

m0 = (np.apply_along_axis(li, 0, coords_order_1, m00, m01)).reshape(
    mesh_order_1.nelem, mesh_order_1.nodes_per_element
) * 0.5

m1 = (np.apply_along_axis(li, 0, coords_order_1, m10, m11)).reshape(
    mesh_order_1.nelem, mesh_order_1.nodes_per_element
) * 0.5

mesh_order_1.attach_field("M0", m0)
mesh_order_1.attach_field("M1", m1)
mesh_order_1.attach_field("fluid", np.ones(mesh.nelem))
mesh_order_1
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f91faa2cf10>

...and our 4th order mesh.

m0 = (np.apply_along_axis(li, 0, coords_order_4, m00, m01)).reshape(
    mesh_order_4.nelem, mesh_order_4.nodes_per_element
) * 0.5

m1 = (np.apply_along_axis(li, 0, coords_order_4, m10, m11)).reshape(
    mesh_order_4.nelem, mesh_order_4.nodes_per_element
) * 0.5

mesh_order_4.attach_field("M0", m0)
mesh_order_4.attach_field("M1", m1)
mesh_order_4.attach_field("fluid", np.ones(mesh.nelem))

mesh_order_4