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