User's Manual

  1. Definition of terms
  2. Theory
  3. Input and output files
  4. Link to ADCIRC-2DDI

1. Definition of terms

Node: A vertex in a grid.
Element: a collection of 3 nodes.
Edge: a side of a triangular element.
Foreground grid: also referred to as edit grid or simply grid, it is an unstructured 2D grid that is to be optimized. Initially it is stored in fort.14.
Background grid: sometimes also referred to as 1st background grid, it contains the errors at each node in it (fort.13). It is also an unstructured grid, but unlike the foreground grid, it remains the same throughout the optimization.
Search grid: a structured grid used to speed up the interpolation, which is fequently necessary during the optimization process. See the illustration in Theory on how it is used. It also remains unchanged throughout the optimization.
Internal edge: any edge that is shared by 2 elements.
Boundary edge (or interval): an edge that can only be found in one element.
Negative element area: element area < 0. But in practice, a small positive number is used in lieu of zero to prevent underflow problem, i.e., if an element's area falls below this positive number, the operation in the 5 optimization modules is then aborted.
Skewness of an element: defined as the ratio between the longest edge of the triangular element and its equivalent radius (=sqrt(area/pi)).

2. Theory

AGrid is a 2D grid optimization code that seeks to equi-distribute edge errors in the foreground grid, given an error field associated with every node in the background grid. The edge error is defined as the integration of nodal errors along an edge. It utilizes five core modules to achieve the error equi-distribution: move node, refine, coarsen and swap edge, and snap node to an edge. These 5 modules are iterated in the way illustrated below:

* *

The clean-up loop consists of move-node, swap-edge, coarsen-edge and snap-node modules written in the Euclidean space to eliminate skew elements (in Euclidean space), which are detrimental to a number of circulation packages.

AGrid is also iterated with the solver to produce a final optimized grid in the following fashion:

*

2.1 Five optimization modules
1) Move-node module:
In this module, a node (i_0) is moved in order to better equidistribute edge errors for the edges around it.
*

2) Swap-edge module:
* In this module, an internal edge, shared by 2 elements, is swapped to the other diagonal of the quadrilateral formed by the 2 elements, in order to improve aspect ratios (or skewness) of the 2 elements in the Riemann space.

3) Refine-edge module:
* In this module, an edge is refined whenever its error exceeds a pre-given threshold lambda_1.

4) Coarsen-edge module:
* In this module, an edge is removed from the grid if its error is smaller than a pre-given threshold lambda_2.

5) Snap-node module:
* In this module, a node is snapped to its opposing edge if the "snapping distance" in Riemann space falls below a certain threshold lambda_3. There are two distinctive cases in this module: snapping to a boundary edge or an internal edge, which are shown respectively above and below this sentence. *

2.2 Second background grid ("search grid")

The code has been made very efficient to handle very large grids, but it was found that the interpolation of errors on the background grid took considerable amount of time (more than 90%). Therefore to maximize efficiency, it is highly recommended to create a search grid, which is a structured grid.

The way the search grid is constructed and utilized is as follows. The entire domain is divided into small rectangles (which may not be equal in size), and each rectangle intersects some elements in the 1st background grid (see the diagram below). Note that this information can be computed and stored once and for all as the two background grids do not change during the optimization.

*

During the process of optimization, it is frequently necessary to find error for a given point in space via interpolation from its parent triangle in the 1st background grid. Instead of searching through all elements in the 1st background grid to find the parent triagnle, we first find its parent rectangle in the search grid (which is trivial), and then search through all the intersecting triangles of this rectangle.

AGrid supports three options for the construction of search grid, but only two are most useful. The easiest way is to evenly divide the domain into nxx times nyy small rectangles (nxx and nyy are specified in agrid.h), and this corresponds to a flag 1 for iflag2bg on the 6th line of fort.15. On the other hand, a flag of 2 entails a user-defined search grid and the x,y values used for the division in x,y directions are specified in fort.18. This option is useful for highly uneven 1st background grids. In any case, there is a cap on the maximum number of intersecting elements allowed for a rectangle (mie in agrid.h), and the smaller this parameter is, the faster the code will run, but possibly at the expense of larger memory required. A good choice of this parameter is around 70. Should the code fail to construct the search grid due to too many intersecting elements, the error message contains the required mie.

2.3 The issue of reserving boundary nodes

It's essential to reserve some of the boundary nodes to preserve its integrity. AGrid automatically reserves those boundary nodes near which the curvature is high. Cubic spline is used to compute the curvature. Then starting from the maximum curvature point on a boundary, the code will go through the boundary and select reserved nodes in such a way that the change of slope angle/tangent (which is the integration of curvature along the boundary) is always below a pre-given constant (degrm on the 8th line of fort.15). Reserved nodes are never moved, deleted or snapped in the optimization process. However, in order to reduce the number of skew elements left after the geometric run, the reserved nodes are allowed to be deleted in this run. This may cause the boundary to be self-looped. Therefore a check is made at the end of this run to indicate problematic boundary nodes and the user can easily fix it with Gredit.

2.4 Automatic search for thresholds lambda_1-3

A bi-section method is used to automaticlly search for proper thresholds lambda_1-3 in order that the number of nodes after the first refinement falls into the prescribed interval [ndmin,ndmax] (the 11th line of fort.15). Currently, the upper bound used in the bi-section method is the larger of 5*A.E. and A.E.+2*S.D. (A.E. is the initial averaged error and S.D. is the standard deviation), and the lower bound is either A.E.-S.D. (if it's positive) or 0.2*A.E. Therefore this option may not be robust and the bi-section method may fail. Also the final number of nodes is less than ndmax, but may be less than ndmin due to later iterations. Obviously, the larger the interval [ndmin,ndmax] is, the faster the search becomes.

Automatic search option should not be used from the 2nd solver-adaptor iteration onwards; instead, set lambda_1-3 to be the values found in the first iteration, in order to guarentee the exactly the same set of parameters is used for every iteration for convergence check purpose.

2.5 Hot start

The user can hot start AGrid anywhere except within the geometric run (the clean-up loop). Also the hot start cannot be made when the following flags are set as:

  1. icstflag=0;
  2. ibifur=1;
  3. icsflag=2.

For hot start, first move save.dat to save0.dat, turn on the hot start option in fort.15, and then place statement 999 to where you want to restart from.

3. Input and output files

3.1 Input files:
Grid file (fort.14), background grid file (fort.13), boundary nodes file (fort.16), main parameter input file (fort.15), search grid file (fort.18).

Description of input files

fort.13 and fort.14:
standard grid file in Gredit format, with the depth being the error at the node in fort.13 (for fort.14, the depth is not used).
fort.16:
in Gredit format (edit boundary nodes).
fort.18 (enabled when iflag2bg>=2 in fort.15):
ndivx: no. of divisions in x;
xmin: minimum x (leave as it is);
xmid(i), i=1,ndivx: x values defining each division;
xmax: maximum x (leave as it is);
mx(i), i=1,ndivx: no. of subdivisions in each interval. The sum of mx(i) must be equal to nxx-2!
ndivy: no. of divisions in y;
ymin: minimum y (leave as it is);
ymid(i), i=1,ndivy: y values defining each division;
ymax: maximum y (leave as it is);
my(i), i=1,ndivy: no. of subdivisions in each interval. The sum of my(i) must be equal to nyy-2!
fort.15 (parameter file):
  1. warning message;
  2. ics, slam0,sfea0: coorsinates system. ics=1: cartesian; =2: lat/long. slam0 and sfea0 are center of projection.
  3. ihot, nstep: hot start option; =0: cold start; =1: hot start and nstep is the new starting step. The code cannot be restarted within the geometric run;
  4. icoflag, emin, emax: flag; =0: error cutoff disabled; =1: enabled, and emin and emax are lower and upper error cutoff values (ignored when icoflag=0);
  5. ismflag, it_sm: flag; =0: no error field smoothing; =1: error smoothing using least square method, and and it_sm is the no. of iterations in smoothing (ignored when ismflag=0);
  6. iflag2bg: flag; =0: no search grid; =1: uniform structured grid as search grid; =2: user-defined search grid (fort.18); =3: use strips in x and y directions defined in fort.18 (user needs to comment out two lines and uncomment two lines in agrid.h). The search grid determines efficiency.
  7. icstflag, cst_neg: flag and a constant to determine negative area during the optimization process. If icstflag=0, cst_neg is a ratio to the smallest element area in the inital grid; if icstflag=1, use absolute value for cst_neg.
  8. degrm: maximum change of slope angle (in degrees) within a boundary interval. Used to reserve boundary nodes.
  9. ivsflag: flag; =0: virus scan of the grid disabled; =1: enabled.
  10. ibifur: flag to use automatic search for thresholds rlam1-3. =0: disabled; =1: enabled.
    if ibifur=0, next line will be:
  11. icsflag, rlam1, rlam2, rlam3: flag; =0: relative thresholds rlam1-3 for refine, coarsen and snap modules; =1: absolute thresholds; =2: absolute thresholds as ratios to the averaged error for the initial grid (the user is not allowed to hot start with flag 2);
    if ibifur=1, next line will be:
  12. ndmin, ndmax: target interval for # of nodes.
  13. iconflag: flag for convergence criterion used in optimization; =0: absolute convergence; =1: relative convergence.
  14. dn_r, an_r, se_r: used when iconflag=1. The three tolerances give the maximum ratios of deleted nodes, added nodes and swapped edges allowed for convergence.
  15. iflaggeo: flag for geometric optimization to get rid of very skew elements; =0: disabled (then the remaining parameters are all ignored); =1: enabled.
  16. askew: acceptable skewness.
  17. it_geo: no. of iterations for geometric run.

3.2 Output files:
screen output, main output file (fort.10), edge error distribution file (err_dis.dat), reserved nodes (rsv.dat), intermediate grid and output files (debug.dat and save.dat), optimized grid (op.gr and op.ll), optimized grid after geometric correction (converged.gr and converged.ll), and error message output files (fort.11).

Appendix: Compressed sample input and output files

4. Trouble shooting

It is quite common to encounter fatal errors during the run, but most of those errors are caused by improper dimensioning in the header file (agrid.h) or improper threshold used in the refinement (rlam1), and can be easily rectified. Below is a list of commonly encountered error messages and the remedy.
  1. Errors in generating search grid
    Solution: if iflag2bg=1, a uniform grid is generated, and you may need to increase nxx and nyy (no. of divisions in x and y) or increase mie (maximum no. of intersecting triangles for any rectangle). Try to keep mie as small as possible to maximize the speed. For highly irregular domains, uniform grid may not be sufficient, and you may consider using a non-uniform grid by changing the flag to 2 and prepare a fort.18. Put more divisions around "crowed" areas. The error message helps you estimate mie required.

  2. Too many nodes/elements during the refinement
    Solution: The 3 thresholds are related by rlam1=4*rlam2, and rlam3=0.8*rlam2. Since most added/deleted nodes are done during the very first iteration, the number of nodes in the final grid is smaller than that after the first iteration. If the refinement cannot go through due to too many added nodes/elements, raise rlam1 (and adjust rlam2-3 accordingly) and retry. It is a good idea to start with rlam1= 2 times the initial average edge error (there is an option to specify rlam1-3 as ratios).
    An easier way is to use the automatic search utility to automatically adjust rlam1-3.

  3. How to adjust emin, emax and cst_neg to obtain a desired grid?
    Solution: Set cst_neg to be of the same order as the smallest area in the initial foreground grid so that most elements will have a chance to be altered by the optimizer. For example, you may set the flag=0, cst_neg=2.0 etc. Remember the smallest element size in the final grid is also related to the time step admissible to a solver.

    The 2 cut-off values emin and emax serves to change the smoothness or contrast of the optimized grid. Often it's necessary to impose these 2 cut-offs for the initial error field because otherwise overflow or over-coarsening would occur due to the "noise" in the field. If the low-error region is over-coarsened after the first iteration, raise emin; if the high-error region is over-refined, lower emax.

  4. Flow solver blows up even if a smaller time step is used.
    Solution: For Adcirc, it should be able to run if the maximum skewness in the grid is less than 6. However, the smallest element size also affects the run. The user can try to raise the minimum element size in the final grid by increasing cst_neg.

  5. Too many neighbors during refinement when using auto search for rlam1-3.
    Solution: Automatic search is not robust, and it can fail due to the following reasons: # of iterations in bi-section method exceeding 20; inadequate maximum and minimum rlam1 used in bi-section method, which are currently being set at 5 and 0.2 times the initial averaged error; # of neighbors exceeding maxnei during the refinement. The last problem is often an indication that the initial error field is too uneven and thus closing up the gap between emin and emax would help. But the user may also try to manually set rlam1-3 to avoid the problem.

5. Link to ADCIRC-2DDI

The Columbia River estuary is used here as an example to illustrate the solver-optimizer loop.

Starting with an initial grid, ADCIRC-2DDI (2D depth integrated model) is used to calculate an error norm (e.g., mass conservation errors). After the ADCIRC run, fort.13 and fort.14 are therefore eventually generated (and they are actually identical in this case). Fort.16 can be easily generated by using Gredit.

After setting all parameters in fort.15, agrid.f can be compiled and run. A new optimized grid is generated at the end of the geometric run (which is needed in this case as ADCIRC is sensitive to skew elements). Most likely there are still a few skew elements in the grid, and Gredit can be used to fix them directly. The resultant grid is then used as new grid in ADCIRC.

To generate fort.15 in ADCIRC run, the user needs to re-compute the open boundary condition of tidal forcing. A series of Perl and FORTRAN codes need to be run to achieve this (spcs.pl, ecp.f, genbcs.f, intel_deg.f). The grid file in ADCIRC, fort.14 is used as input to these programs. The bottom friction file (fort.21) can be generated by a simple FORTRAN code friction.f, the wind file (fort.22) can be generated by data obtained from the web, and the river open-boundary forcing file (fort.17) can be generated from a database (LOADMAX). Then ADCIRC is run and the process repeated until converged.