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.
In this notebook you build a 2D mesh with alternating solid and fluid blocks of elements.
Copy
# initialize notebook
%matplotlib inline
%config Completer.use_jedi = False
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["figure.figsize"] = (10, 10)

from salvus.mesh import mesh_block
Simulations in coupled solid-fluid media require two different kinds of wave equations, the acoustic wave equation in the fluid part of the domain and the elastic wave equation in the solid part, as well as some interface conditions to ensure the continuity of traction.
Imposing the interface conditions on unstructured meshes can be challenging, because it requires some book-keeping which elements are facing each other along the interface.
Fortunately, this is automatically done within Salvus_Compute_ and we only need to flag elements as being either solid and fluid. In this example, we create a simple quadratic domain that consists of small patches of acoustic and elastic material. Of course, this is not a realistic physical domain, but it shows how easy it is to create meshes for coupled domains with Salvus.
# input in SI units
solid_vp = 1500.0  # m/s
solid_vs = 1000.0  # m/s
solid_rho = 1000.0  # kg / m^3
fluid_vp = 1000.0  # m/s
fluid_rho = 1200.0  # kg / m^3
lx = 10000  # domain size in x direction in m
ly = 10000  # domain size in y direction in m
fmax = 1.0  # maximum frequency in Hz
block_size = 3  # size of the solid and fluid blocks in
# number of elements

courant_number = 0.5
elements_per_wavelength = 2
# maximum element size
hmax = min(solid_vs, fluid_vp) / fmax / elements_per_wavelength

# compute number of elements
nelemx = int(lx / hmax) + 1
nelemy = int(ly / hmax) + 1

# generate rectilinear grid
sg = mesh_block.generators.cartesian.rectangle_2d(
nelemx, nelemy, max_x=lx, max_y=ly
)

# create unstructured mesh
m = sg.get_unstructured_mesh()

# make fluid blocks
fluid = np.zeros((nelemx, nelemy), dtype="bool")
for i in range(block_size):
for j in range(block_size):
fluid[i :: block_size * 2, j :: block_size * 2] = True
fluid[
i + block_size :: block_size * 2, j + block_size :: block_size * 2
] = True

fluid = fluid.ravel()

# a mask for the solid elements
nfluid = np.logical_not(fluid)

# attach material properties
vpa = nfluid * solid_vp + fluid * fluid_vp
vsa = nfluid * solid_vs
rhoa = nfluid * solid_rho + fluid * fluid_rho

m.attach_field("fluid", fluid.astype("int"))
m.attach_field("VP", np.outer(vpa, np.ones(m.nodes_per_element)))
m.attach_field("VS", np.outer(vsa, np.ones(m.nodes_per_element)))
m.attach_field("RHO", np.outer(rhoa, np.ones(m.nodes_per_element)))

# define outer boundaries
m.find_side_sets(mode="cartesian")

# compute time step
dt, dt_elem = m.compute_dt(vpa, courant_number, fast=False)
m.attach_global_variable("dt", dt)
m.attach_field("dt", dt_elem)  # this is mainly for visualization

print("================")
print("nelem = %d" % m.nelem)
print("dt    = %e" % dt)
print("================")
================
nelem = 441
dt    = 1.122392e-01
================

We have created a simple mesh with patches of acoustic and elastic elements. Now, let's verify that we indeed attached the coorect physics by visualizing the mesh and checking the fluid flag. You also want to ensure that the shear wave velocity VS is zero in the fluid domain.
m
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fd39411d490>
PAGE CONTENTS