Mondaic
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.
# 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_blockSimulations 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 0x7fa66590e4d0>
PAGE CONTENTS