External 2D Meshes

This tutorial assumes that you’ve already read the introduction for this external mesh series.

Constructing the geometry

First, let’s start by creating a rectangle and a circle. In Coreform Cubit, this can be achieved using:

create surface rectangle width 2 height 1.5 zplane
create surface circle radius 0.25 zplane

Note that currently the circle and the rectangle are completely disjoint from each other. In order to remove this overlap between these two entities, the boolean operator can be used to subtract the circle from the rectangle:

subtract surface 2 from surface 1 keep_tool

While the rectangle now has a circular cutout in its center, the boundaries between the yellow circle and the green rectangle are not shared. The merge tool can be used to consolidate these curves into one:

merge curve 5 6

Generating the mesh

We will use the following material properties for this example:

As a general rule of thumb, using approximately 2.0 elements per wavelength tends to yield minimal numerical dispersion for most applications. For this example we’ll assume that we’re using a source with a maximum resolved frequency of 30 kHz, meaning that we can compute the approximate element size $\Delta x$ in each region as

and

In order to mitigate spatial dispersion errors, it is advantageous to begin meshing from the region with the finest discretization. Thus, first the circular region is meshed with

surface 2 scheme pave
surface 2 size 0.025
mesh surface 2

Next, the rectangular region is meshed using

surface 3 scheme pave
surface 3 size 0.05
surface 3 sizing function interval
mesh surface 3

Mesh smoothing can be used to improve the uniformity of the element sizes throughout the mesh using

surface all smooth scheme edge length
smooth surface all

In order to make it easier to select the individual discrete regions of the mesh, we’ll assign the circular inclusion as belonging to block 1, while the surrounding part of the domain will belong to block 2:

block 1 add surface 2
block 2 add surface 3

Finally, the mesh can be exported to Exodus such that it can be imported into Salvus. It is important to export the mesh with the following settings:

• Dimensions: 2D
• Format: Use NetCDF4/HDF5 Format

This can be achieved using the Journal commands

set exodus netcdf4 on
export mesh "/path/to/file/sample_2D.e" dimension 2

Complete journal file

Here’s a complete list of the journal commands necessary for creating the mesh:

# Create the background rectangle and circular inclusion
create surface rectangle width 2 height 1.5 zplane
create surface circle radius 0.25 zplane

# Subtract the circle from the rectangle
subtract surface 2 from surface 1 keep_tool

# Merge overlapping curves
merge curve 5 6

# Mesh the circular region
surface 2 scheme pave
surface 2 size 0.025
mesh surface 2

# Mesh the rectangular region
surface 3 scheme pave
surface 3 size 0.05
surface 3 sizing function interval
mesh surface 3

# Smooth the mesh
surface all smooth scheme edge length
smooth surface all

# Assign the individual surfaces to blocks 1 and 2
block 1 add surface 2
block 2 add surface 3

# Export to Exodus
set exodus netcdf4 on
export mesh "/path/to/file/sample_2D.e" dimension 2

Mesh import and material assignment

As usual, we’ll start off by importing a few useful modules.

Let’s begin by first importing the Exodus mesh. We can do this using the from_exodus constructor that is part of the UnstructuredMesh class. Here’s an example of how we can do that:

Currently the mesh only has the block indices attached; we haven’t specified any other material properties yet. The block indices let us flag parts of the mesh as belonging together and serve as a convenient way for distinguishing between different discrete regions of the domain.

We can check for potential issues with the mesh using the qc_test method:

Uh oh — since we don’t have any material properties attached yet, we wouldn’t be able to run a simulation just yet. Let’s attach the material properties that we defined earlier in the tutorial.

We’re going to want to assign all of the elements in each region of the mesh to have these material properties. To do this, we can use the following general recipe:

1. Initialize each elemental field
2. Create a mask that selects the elements for a particular block
3. Assign the material properties from the block_props dictionary
4. Depending on whether the material is acoustic or elastic (i.e., whether there is a VS present), set the fluid flag for each region

Let’s run the qc_test again to see if there are any remaining issues with the mesh.

It looks like we’re good to go!

Resolved frequencies and global time step

We checked before whether we’ll be able to execute the simulation with the mesh we’ve imported. While the qc_test tells us if there are issues with the mesh that will prevent us from running the simulation, it is still up to the user to make sure that the external mesh has been constructed in such a way that it will perform well for the simulation:

• Having elements that are too large will introduce numerical dispersion errors into the simulation
• Having elements that are too small will limit the global time step under the Courant (or CFL) criterion, which will make the simulation much more computationally expensive to run.

To make sure that our simulation will perform well, there are two things we want to check here:

1. Compute the resolved frequencies of each element in the mesh to make sure that they’re (mostly) around the desired maximum frequency of the source.
2. Compute the time step $\Delta t$ to see where the limiting elements in the mesh are.

In the above plot, we can see that the range of resolved frequencies seems to mostly be in the same range. Let’s visualize these mesh statistics with a histogram to see the frequencies that most of the elements are resolving:

We can see here that most of our elements are resolving frequencies that are around 30kHz, which is the approximate maximum frequencies that our source will have.

While there are some elements that are somewhat too large, they are sufficiently few in number that the additional dispersion introduced by them for the sake of this example will likely be minimal. However, depending on the application, one may wish to go back and remesh this domain to eliminate these large elements.

As noted before, the abnormally small elements will limit the global time step. We can nicely visualize this by creating a similar histogram of the time steps for each element:

The spread of time steps is relatively narrow, meaning there is no singular element that is significantly smaller than its surrounding elments. Thus, while it may be possible to further optimize the time step for this mesh, it is of satisfactory quality.

Finally, we can visualize the time steps of each element using the mesh visualization widget:

Setting up and running a simulation

To demonstrate that this mesh actually works in a forward simulation, let’s set up a very minimalistic project and output a wavefield for it.

Before running the simulation, let’s check that the setup looks OK:

Looks like we’re good to go!

Here we’ll output a few wavefield snapshots for the simulation. In order to easily compare the outputs for the acoustic and elastic parts, we’ll output the fields displacement and gradient-of-phi. We can convert the gradient-of-phi in the acoustic part to displacement using the conversion of

where $u (\mathbf{x}, t)$ is the displacement wavefield, $\rho (\mathbf{x})$ is the mass density (i.e., the RHO field in the acoustic part of the mesh), and $\nabla \varphi (\mathbf{x}, t)$ is the gradient-of-phi wavefield that we output in the acoustic part of the mesh.

To visualize this data, one can use an application such as ParaView. Included in this tutorial is a ParaView state file that can be used to load in the necessary data for this. Here’s a snapshot of what the wavefield looks like at the final time step: