Skip to content

Coupled medium

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.

Mesh generation for coupled solid-fluid domain

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

The Checker and the Board:

# initialize notebook
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
plt.style.use('ggplot')
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 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.

Inputs

# input in SI units
solid_vp = 1500.                # m/s
solid_vs = 1000.                # m/s
solid_rho = 1000.               # kg / m^3
fluid_vp = 1000.                # m/s
fluid_rho = 1200.               # kg / m^3
lx = 10000                      # domain size in x direction in m
ly = 10000                      # domain size in y direction in m
fmax = 1.                       # 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.

m
VBox(children=(HTML(value='\n            <p><i>This is only a preview. Use ParaView to actually investigate\n …





<salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x111bf3860>