Defining Solid and Fluid Domains

Mesh generation for coupled solid-fluid domain

In this notebook you build a 2D mesh with alternating solid and fluid blocks of elements.

Copy
# initialize notebook
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

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

from salvus_mesh.structured_grid_2D import StructuredGrid2D
from salvus_mesh import Skeleton

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 SalvusCompute 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.

Inputs

# 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 = StructuredGrid2D.rectangle(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", vpa)
m.attach_field("VS", vsa)
m.attach_field("RHO", rhoa)

# 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.

When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x7f8c739cce80>
PAGE CONTENTS