Get SWAN at Fast, secure and Free Open Source software downloads

Unstructured grids in SWAN - a tutorial

This tutorial gives you an opportunity to learn how to setup and to use unstructured meshes in SWAN. By following several steps in this site, you will learn to:
  • motivate the use of flexible meshes in SWAN
  • generate an unstructured grid with Triangle
  • specify grid layout in SWAN
  • impose boundary conditions on specific side(s) or segment(s) of computational grid
  • get acquaint with numerical integration and convergence check
  • visualize wave parameters as (colored) maps in Matlab
Unstructured meshes provide much better representation of complex boundaries such as coastlines and areas around islands than do conventional regular grids and also provide the opportunity to concentrate mesh resolution in areas of interest e.g., regions of strong bathymetry variations in estuaries and fjords, to a degree not possible on a curvilinear grid. Hence, there is no need for nesting. The use of unstructured grids facilitates to resolve the model area with a relative high accuracy but with much fewer grid points than with regular grids. Although, the CPU cost per iteration is relative higher than cases with structured grids (as is often the case), this effect is more than offset by the reduction in the number of grid points.

Grid generation
The first step you have to do is generating a grid with a variable resolution. There are, however, a huge number of mesh generators for grids containing triangles, quadrilaterals, or a combination of these types of grids. To give you an idea, take a look at a list of grid generation software. SWAN uses triangular cells only.

At the moment, SWAN supports three mesh generators:
SMS (Surface-water Modeling System) is a pre- and post-processor for surface water modeling developed by Environmental Modeling Research Laboratory. SMS supports models like ADCIRC (Advance Circulation Model originally developed by Rick Luettich and Joannes Westerink) that calculates shallow water flows and transport problems in two and three dimensions. It employs the finite element technique in geographical space allowing the use of unstructured grids. A dynamic coupling of ADCIRC with SWAN is available.

SMS provides data for ADCIRC by creating a triangular mesh and assigning boundary conditions at appropriate boundaries. This data is saved in a file called fort.14 as part of the input to ADCIRC. This file can be interpret by SWAN as well. Hence, SWAN can apply the same grids used by ADCIRC.

This well-known public-domain software was developed by Jonathan Shewchuk. Triangle is a C program for 2D mesh generation through construction of Delaunay triangulation. (The Delaunay criteria: any vertex must not be contained within the circumsphere of any triangle within the mesh. This is not a restriction to SWAN, i.e. SWAN may allow non-Delaunay triangulation as well.)

Just another 2D Delaunay-based mesh generator developed by Bojan Niceno (freeware).
As a starting point, a flexible mesh for SWAN can be prepared using Triangle. If you don't have a C compiler, you may download the executable that can be run on Windows. The user has the possibility to specify the user-specified constraints on angles and triangle areas and user-specified holes and concavities. Triangle's input is a planar straight line graph (PSLG) defined to be a collection of vertices and segments, where the endpoints of every segment are included in the list of vertices. An example:

# 4 vertices of a rectangle of 5m by 1m
4 2 0 0 # no attributes and markers here
1 0. 0.
2 5. 0.
3 5. 1.
4 0. 1.
# 4 segments of rectangle with boundary markers
4 1
1 1 2 1 # south
2 2 3 1 # east
3 3 4 1 # north
4 4 1 2 # west, another marker for this side
# no holes

Save this to a file called rectangle.poly (extension .poly is obliged). Next, type the following at the command prompt:

triangle -pq30a0.05I rectangle

where the following options have been used:

-p: reads a .poly file
-q30: imposes a minimum angle of 30 degrees
-a0.05: imposes a maximum triangle area of 0.05 m2
-I: supresses mesh iteration numbers

The execution of the triangle command will produces two more files: rectangle.node and rectangle.ele. The first file contains the (x,y)-coordinates of the vertices (or nodes) including boundary markers. The second file is a connectivity table for triangles (or elements). This table consists of three corner vertices around each triangle in counterclockwise order. These files can be read by SWAN if the following command lines are specified in the SWAN command file:

CGRID UNSTRUCTURED CIRCLE [mdc] [flow] [fhigh] [msc]

After setting up the vertices and the connectivity tables for cells and faces (automatically done in SWAN), SWAN will print some information concerning the used mesh in PRINT file, among others, number of vertices, cells and faces and minimum and maximum gridsizes. Furthermore, SWAN will check at two levels for a possible occurence of badly shaped triangles. Firstly, the number of triangles that meet at each vertex inside the mesh should not be smaller than 4 or larger than 10. Secondly, the angles inside each triangle should not be higher than 143 degrees. If, at least, one of these two situations occur, SWAN will print an error message.

You may visualize the grid by typing the following command:

showme rectangle

and you will see something like this
Note that showme can only be executed on a platform supporting X11 windows system (usually Unix-like ones). Alternatively, you may use the Matlab script plotgrid.m to plot the grid, as follows (to be executed in Matlab):


where the parameter 'rectangle' refers to the basename of the Triangle files.

It is also quite easy to put a number of holes in a computational domain as shown below. Usually, these holes represent islands or obstacles. For instructions how to setup holes in your domain, you may consult the online manual of Triangle.

Generally, the user would like to have an efficient grid in which areas where the bathymetry or evolution of the waves change rapidly require a higher resolution than areas where the physics or depth changes less. It seems not so easy to get such a mesh using Triangle only. Usually, the modeler need to have an indication how to determine the refinement based on bathymetry or geometric variations through preliminary evaluations. Such an iterative mesh editing procedure, including checking and improving grid quality, may be controlled by BatTri, a graphical Matlab public-domain interface to Triangle. (Note that this GUI is not supported anymore.) The nice thing of this interface is that boundary nodes, segments and holes can be created from only bathymetry data with the use of the editing options. This process may require manual deleting of unnecessary segments and nodes, closing of islands by segment adding, addition of an open ocean boundary segment, etc. The final information on nodes and segments is forced into the triangulation of the domain. Triangle is called within BatTri to generate the actual grid.

BatTri provides many pre-defined depth-dependent contraints for further mesh refinement. An example is the h refinement which relates the maximum triangle area A to depth h via a user-defined constant alpha, as follows:


resulting in a new (preliminary) mesh that is more finer in the relative shallow parts of the domain compared to the deep parts. This refinement process is then repeated iteratively until the user is happy with the final grid.

An example is the Haringvliet estuary as part of the Rhine river in the south-west of the Netherlands. Initially, a grid is generated with approximately 1700 triangular elements (maximum triangle area is 300,000 m2) as depicted in the next figure; the colored map shows the depth ranging from -25 meter to 5 meter (positively upwards from still water level).

To optimize the grid, the size of these elements is made proportional to the depth, i.e. the triangles are small in regions where the depth is shallow and large in deeper water. This is done in regions where the depth is at most 10 meter, i.e. in the range [-10 m, 5 m]. At the same time, 5000 grid points have been added as well. The next picture shows the final result.

The number of grid cells is 11,395. The minimum grid size is roughly 25 meters and the maximum one about 700 meters. Suppose, a regular mesh should be used and the same level of accuracy should be obtained without nesting, this would require a rectangular grid containing more than 500,000 grid cells covering the same part of the area.

Simulation with SWAN can be done in the usual way. However, two aspects need to be considered:
  • imposing spectra at boundaries
  • numerical integration e.g., sweeping and convergence check
There are two ways to define the part of the boundary at which the spectra are imposed. The first way is using the command SIDE and is easiest if the boundary is one full side of the computational grid. The other is realized with the command SEGMENT and can be used if the boundary segment goes around the corner of the grid, or if the segment is only part of one side of the grid. With respect to the latter approach, we may still specify the (x,y)-coordinates (indicate as problem coordinates) as a series of points forming a segment together. (You can not indicate the segment points by means of two-pairs grid indices!)

When applying the command SIDE, the so-called boundary markers must be employed. Boundary markers are tags to identify which vertices occur on a boundary of the mesh. With respect to the SWAN command file, they are associated with boundary conditions. In the above example with the rectangle, we would like to impose a Jonswap spectrum (Hs=1.2 m, Tp=7.3 s and cos2-directional spreading) at the west side. At the other sides no boundary conditions will be given. As indicated in the second part of the file rectangle.poly (see above, the last column), the west side (4th segment) is marked 2 while the other sides all have the same boundary marker 1. By indicating a boundary marker in the command SIDE, we can now specify on which side the boundary condition is applied, as follows:

BOU SIDE 2 CCW CON PAR 1.2 7.3 90 2

After finish editing the SWAN command file, you may run SWAN on your computer in the usual way.

Contrary to the most third-generation spectral wave models, the numerical propagation scheme in SWAN is unconditionally stable i.e., it is not subjected to a CFL stability criterion. It made use of the four-direction (symmetric) Gauss-Seidel iteration technique. This method appears to be very effective on regular meshes. Furthermore, it is also very robust in practical shallow coastal applications since, the permitted grid resolution and time steps are mutually independent, which is a very important advantage for flexible meshes. Because of these nice properties, this solution technique is tailored to unstructured grids. A vertex-based algorithm is employed in which the variables are stored at the vertices of the mesh and the wave balance equation is solved in each vertex assuming a constant spectral grid resolution in all vertices. Next, an ordering of vertices is established in such a way that sweeping through each vertex can be made by using the updated information from the surrounding vertices as soon as it is available. This ordering is based on spherical wavefronts around a reference point that is chosen at the boundary where boundary conditions are imposed. The two upwave faces connecting the vertex to be updated with its two updated neighbours encloses those wave energy propagation directions that can be processed in the spectral space without having stability problems. See figure below.

After updating all vertices, the process continues by sweeping again through the vertices. This is repeated until all vertices are updated in the sense that wave energy from all directions has been transmitted through geographical space. It is very important to realize that the number of sweeps per iteration is exactly 1 (not 4 as in regular grids)! In addition, the average number of iterations in a stationary calculation is more or less equal to regular grids.

Graphical presentation of wave parameters
To display the wave parameters e.g., significant wave height, peak period, directional spreading, etc. as (colored) maps (i.e. as function of geographical coordinates) they need to be stored in a Matlab binary format as specified in the SWAN command file:


Of course, you are free to choose whatever wave parameters you like. However, the two parameters XP and YP must be stored in the Matlab file.

With the Matlab script plotunswan.m, a plot of the significant wave height can be made as follows (to be executed in Matlab):


The first parameter of command plotunswan indicates the name of the Matlab file containing the wave parameter(s), the second parameter indicates the basename of the Triangle files and the last parameter indicates the wave parameter to be plotted. If you want to know which other wave parameters are stored in the Matlab file waveparam.mat, just type

who -file waveparam

Help about plotunswan can be accessed by running help plotunswan in Matlab.

The next figures serve as an example of plotting with plotunswan, and show detailed snapshots of the significant wave height (in meters) and the peak period (in seconds), respectively, of the Haringvliet case.

What's next?
This tutorial should be enough to perform some basic steps for a successful simulation of SWAN using unstructured grids. Also, you should be able to understand in that respect most of the User Manual. You may start with the Haringvliet case immediately to examine it further. Try things out. Next, we highly recommend you to explore more thoroughly the Matlab interface BatTri, or some other advanced GUI grid generator since, it is the most essential part of the whole work.

If you have any questions, ideas or suggestions, please feel free to contact us!

Good luck!