South Atlantic Bight Coastal Mesh Generation Example



Introduction

Initial steps  
BatTri mesh generation

Introduction

The purpose of  this exercise is to resolve tidal and wind driven currents in the South Atlantic Bight (SAB) in a mesh that resolves the estuarine system sufficiently to allow correct tidal propagation. The size of this mesh (number of nodes, number of elements, number of boundary nodes) should be small enough to allow operational runs of the SABLAM operational system. This effort follows the work done by Keston Smith with BatTri on the generation of the "peach" series of meshes and previous meshes like "empanada" . BatTri, developed by Ata Bilgili and Keston Smith, is a graphical Matlab interface to the C language two-dimensional quality grid generator Triangle developed by Jonathan Richard Shewchuk. Triangle received the Wilkinson Prize for Numerical Software at ICIAM03.
 

Initial Steps

Matlab setup

Start Matlab and set the path:

>>addpath /usr/local/matlab/toolbox/BatTri
>>addpath /usr/local/matlab/toolbox/opnml_toolbox
>>start_opnml %or whatever startup script you have to add the full opnml path.

Bathymetry

The bottom topography used is the Coastal Relief Model bathymetry database (Volumes 2 and 3). The 3 seconds spacing data is subsample into 6 seconds and a basic matlab structure (x,y,z) is created using the function makesabbat.m . The resulting structure is ready to be used by BatTri, because it is in an adequate format and in cartesian coordinates.

 makesabbat
Top

BatTri Mesh Generation

We start running BatTri using the generate_mesh.m script.  For detailed explanations of every step, check the BatTri manual .

>>g=generate_mesh(bat);

Basic set up

We are asked for a series of inputs:  BatTri will display the defined Triangle path.  We should check that the triangle_path parameter is correct. If it is not, we have to change that manually inside generate_mesh.m . We first will be asked if we want to turn off the exact arithmetic for scattered bathymetry (This option is relevant for Intel architectures). We will then be prompted to choose between linear interpolation and some simple Objective Analysis and usually the interpolation is chosen for computational speed . After this we have to choose the filename of the final grid ('finalgrid') that will be generated in this BatTri session . We will then be prompted to enter the array of bathymetric contour values that will be drawn on the screen when the input .poly file is displayed and usually none is used because of  the computational cost of working with interfaces with complex graphics. BatTri will then ask us to enter the minimum depth for nodes in the grid. In the final mesh with interpolated depths, any depth larger than this value will be truncated to the value of the minimum depth. In this example we used 0 m as our minimum depth.

basic set up
Top

Boundary creation

The next step is to enter the name of the .poly file that you will edit or create a grid from. BatTri gives 3 options at this stage and in our example we are starting from scratch. The program will then ask us if we would like to plot bathymetric data point locations and for simplicity reasons we choose not to plot them.

We are presented with a 17 options menu that constitutes the first part of the mesh editing menu.
The first step is to choose the coastal resolution that we are going to use for the mesh. In our case we choose 2000 m resolution at the coast and 20,000 m at the open boundary at the 200 m isobath. We smooth the bathymetry/coastline in the same scale that we interpolate for boundary nodes (we use the same length scale for the filter that we use for the node spacing) so that the estuaries are as close to the real ones as possible. In other experiments a larger length scale for the filter is used to avoid sharp angles in the coastal region of the resulting mesh.

The boundary is constructed from the bathymetry segments provided by the matlab structure obtained from the  Coastal Relief Model bathymetry database using the function makesabbat.m . For each degree of latitude we have a segment  of coastline.

first segment
first segment plot
Top

After we load the two first segments (option 0, add contour line) we can start connecting them. The extremes of each segment are represented by green asterisks. We have to be careful that we make all the connections so that the boundary is continuous (option 4, add edge). We continue loading coastline segments and connecting them to the previous ones until we have the whole coastline.
Note: A useful option is "save current pslg" (option 13). It saves the current polygon as we are working on it, therefore allowing for posterior loading of the polygon if we want to make corrections after the fact.
 It is important to be certain that all the segments are connected. After we finish with the coastline, we use exactly the same procedure for the open boundary at the 200m isobath. Because of the size of the arrays we are dealing with, it is important to zoom in the regions we want to modify (option -1, zoom and pan), for instance when we are making boundary connections. These connections need to be as close to the real boundary as possible and for that sometimes we need to divide edges (option 6) or add vertices (option 1) and then make the connections adding edges (option 4).
Note: To make the characterization of island and other holes on the grid easier it is important that we use the option for including the holes directly when we are loading the boundary, so that BatTri determines whether it is a hole or not. In the regions where we made the boundary connection by hand we need to add the corresponding holes too (option 8).

connection 1
3a
3b
3c
Top

The next step is to connect the offshore boundary to the coastline. We make this connection adding an edge (option 4) from a vertex on the coastal boundary to a vertex on the offshore boundary. Afterwards we can remove portions of the coastline lying outside the mesh with options 2 and 12. After we have a close boundary we can start playing with it so that it includes the rest of necessary features of the grid, for instance zones and holes.

connect boundaries

We want to treat the resolution in the estuaries differently than the resolution on the shelf, because we need a grid in which the tides have the right characteristics on the whole mesh (or at least a decent approximation to them). If we use zones it is necessary that every part of the grid is in its own zone, so the first thing we do is to define a zone on the shelf (option 10). This is going to be our zone number one. After that we can start closing the edges that separate the different zones using option 4 for adding edges or option 7 for adding spline curves that we can then connect to the existing boundary. We need to make sure that the edges separate the zones completely. Then we use option 10 to add zone numbers to each zone we want to define. In our case we define a total of 9 zones that are mostly concentrated on the Ga/SC coast, that is the region where the estuarine system is more significant. Remember that because of the size of the arrays it is useful to zoom in the regions we want to modify (option -1, zoom and pan, option -2 to go back to normal view).

adding zones

Then, we designate which regions of the domain are outside of the finite element mesh (apart from the ones that we have designated while we were adding the boundary segments). We use option 8 to place a mark anywhere inside each island or hole we want excluded from the domain. It is convenient to also add a hole to the mainland, to guaranty that the mesh is closed.
 
with the zones
Top

Initial mesh creation

After we are done with the polygon specifications, we can create a first triangulation that can serve us to determine regions that need better resolution. First we need to impose a minimum angle and usually we choose 25 degrees. Then we choose a maximum number of nodes to add, that in our example is not a big limitation if we select the maximum area in the next step in a sensible way. Then we need to specify a maximum element area for each of the zones that we defined in the boundary and zones creation step. In our case the minimum resolutions are going to be related to the different kind of regions that we want to resolve. On the shelf the minimum resolution is going to be related to the fact that we have a 20,000 m offshore boundary resolution. So we can choose the elements maximum area to adjust to this resolution. We use the formula of the area of  an equilateral triangle of side length L;
area
The maximum area for the shelf zone (zone number 1) is sqrt(3)*(20000) 2/4.
For the estuarine zones (zone numbers 2-9), we choose a minimum resolution close to the resolution of the coastal boundary, so that these areas are better resolved. The maximum area chosen is sqrt(3)*(2500) 2 /4.
 
text initial mesh
Top

Mesh Refinement

At this point we should be able to see if we are going in the right direction or if we need to start over. The initial mesh should have reasonable characteristics and should be much more resolve on the estuaries than on the shelf (use plot option 6). We can see that the transition between the two regions is very abrupt. We could have avoided this problem if we would have created another zone running parallel to the coast and surrounding the previously defined zones (to try this suggestions you should go back to the boundary (polygon) creation menu).

initial mesh initial mesh zoom
Top

The first refinement we are going to use is related to wave propagation.  This constraint is h/area (refinement option 2) and it is equivalent to a wavelength to gridsize ratio. We choose a h/area >= alpha constraint related to the maximum area limitations that we imposed in the previous step. We apply this constraint only to elements which have at least one node deeper than 15 meters ([-200,-15]). The rest of the elements are refined in the next step.

refine1
1 refin text 1 ref plot 1 refin zoom
Top

For the second refinement we impose a constant area constraint for elements shallower than 15 meters (refinement option 6). To be consistent with the previous refinement  we match the resolution of our h/area constraint at the 15 meter isobath for the region between the 15 meter isobath and the shoreline  ([-15,0]), choosing:
 
refine5

2 ref text 2 refin plot 2 ref zoom
Top

 The final refinement we are going to use is to impose a  h/(grad(h)*area) constraint  (refinement option 1) shallower than 10 m ([-10,0]). This constrain helps resolve the steep topography of some regions of the mesh, like the mouth of some of the estuaries, so that we resolve the deep channels in the centers of the rivers. This constrain is somehow related to the delta(h)/h constraint.
 
refine0

3 ref text 3 refin plot 3 ref zoom
Top
Finally, we are going to use the algorithm to do the spring relaxation of interior node positions (option 0) to avoid large changes in angles between adjacent elements. These improves significantly the quality of the mesh and relaxes the transition between different refinement zones.
 
spring and final output
Top

 After we are finished with the refinements we take a last look at the mesh and we check that the different areas are appropriately resolved (plot options 3 and 6). This shows that the estuaries are resolved and that the steep bathymetry regions are resolved in more detail.
 
bathy
Top

Output mesh

If we are happy with the resulting mesh we can output the current mesh to NML standard mesh files. The resulting mesh has been bandwidth reduced.

So now you can start playing with your wonderful mesh!