CMSFlow Numerical Methods:General Transport Equation
3 Numerical Methods
Overview
CMSFlow has both implicit and explicit solution schemes. The explicit solver is designed for dynamic problems with extensive wetting and drying that require small computational time steps, while the implicit solver is intended for simulating tidal and waveinduced circulation at tidal inlets, navigation channels, and adjacent beaches. A detailed description of the numerical formulation of the explicit solver of CMSFlow can be found in Buttolph et al. (2006) and is not repeated here. The sections below specifically refer to the implicit solver of CMSFlow.
The implicit solver uses the SIMPLEC algorithm (Van Doormal and Raithby 1984) on a nonstaggered grid to handle the coupling of water level and velocity. Primary variables u, vvelocity, and water level are stored on the same set of grid points, and fluxes at cell faces are determined using a Rhie and Chow (1983) type momentum interpolation method (Wu et al. 2011). The explicit solver uses a staggered grid with velocities at the cell faces and the water levels and water depths at the cell centers (Buttolph et al. 2006a). CMSFlow also calculates salinity, sediment transport, and morphology change.
CMSFlow Computational Grid
The implicit version of CMSFlow uses a generic Cartesian grid which can be regular, nonuniform, or locally refined by splitting a cell into four subcells. Three requirements are necessary for the input grid:
 1. Cells must have a rectangular shape. Irregularly shaped cells are not allowed.
 2. Cells may have a total of four to six neighboring cells (faces).
 3. Only two neighboring cells are allowed in the same direction (i.e. North, South, East, West).
Mesh refinement can be achieved either by locally decreasing the grid spacing (nonuniform Cartesian grid) as shown in the left panel of Figure 3.1 or by subdividing or splitting a cell into multiple cells as shown in the right panel of Figure 3.1. The refined mesh can be further refined and split in multiple levels, as needed. This mesh is referred to here as a telescoping mesh.
Figure 3.1. Examples of a Cartesian grids allowed in CMS: Stretched Cartesian (left) and telescoping grid (right)
One important aspect of incompressible flow models is the location of primary variables: velocity and pressure (water level). On a staggered grid, the pressure (water level) is located at the center of cells and the u and vvelocities are located along the faces of cells (Harlow and Welsh 1965; Patankar 1980). On a nonstaggered grid, all of the primary variables are located at the center of cells. A staggered grid can more easily eliminate the checkerboard oscillations when compared to a nonstaggered grid; however, a nonstaggered grid involves a simpler source code and can minimize the number of coefficients that must be computed and stored during a simulation because many of the terms in the equations are approximately equal. In particular, a staggered grid is more complicated when handling the interface between coarse and fine cells where five or sixface control volumes are used. Therefore, a nonstaggered (collocated) grid approach is adopted for CMSFlow, with a Rhie and Chow (1983) momentum interpolation technique used to eliminate the checkerboard oscillations. Figure 3.1 shows the location of primary variables and the 5 and 7point stencils (computational molecule) used in the calculations for CMSFlow.
The data structure for a grid can be approached in three ways: 1) blockstructured, 2) hierarchical tree, and 3) unstructured. The blockstructured approach divides the domain into multiple blocks and each block is treated as structured. A special treatment is applied between blocks to ensure mass and momentum balance using this approach. The hierarchical tree approach is memory intensive and requires parentchild relationships and a tree traverse to determine mesh connectivity. For the unstructured approach, all cells are numbered in a onedimensional sequence and tables are used to determine the connectivity of neighboring cells. Among these three approaches, the unstructured approach is the most simple and is therefore applied in CMSFlow.
Computational cells are numbered in an unstructured manner via a 1D index array. Inactive cells (permanently dry) are not included in the 1D index array to save memory and computational time. All active computational nodes are numbered sequentially. For convenience with handling boundary conditions, each boundary cell has a neighboring ghost cell outside of the computational domain. Each ghost cell corresponds to a boundary face of the boundary cell. The ghost cells are stored at the end of the 1D index array.
General Transport Equation
In order to avoid redundant derivations of discretized equations, discretization of a general transport equation is outlined and described below. All of the governing equations are some form of this transport equation; therefore, the same discretization may be applied to all of the governing equations. The general transport equation is given by
 (31)
where is a general scalar, t is time, h is the total water depth, is the transport velocity, is the diffusion coefficient for , and includes all remaining terms. Note that in the case of the continuity and momentum equations, is equal to 1 and , respectively.
Spatial Discretization
A controlvolume technique is used in which the governing equations are integrated over a control volume to obtain an algebraic equation that can be solved numerically. Integration of Equation (31) over a control volume yields:
 (32)
 (33)
 (34)
where
 = outward unit vector normal to cell face f
 = outward cell face velocity
 h_f = linearly interpolated total water depth at the cell face f
 is the outward normal gradient of at cell face f
The above equations, the GreenGauss Theorem has been used to convert the area integral to a boundary integral. The symbol () indicates the cell face linear interpolation operator. The cell face velocity, , is calculated using a momentum interpolation method similar to that of Rhie and Chow (1983) and is described in a subsequent section.
Temporal Discretization
The general transport equation is rewritten as
 (35)
where includes all the remaining terms, such as …. For stability and efficiency, a fully implicit timestepping scheme is used in the form
 (36)
where is a weighting factor between 0 and 1. For , the scheme reduces to the firstorder backward Euler scheme, and with , the scheme reduces to the secondorder backward scheme (Ferziger and Peric 1997). The superscripts indicate the time step levels, with n+1 being the current timestep and n being the previous timestep.
Cellface interpolation operator
The general formula for estimating the cellface interpolation operator (?) value of is given by
 (37)
where is a linear interpolation factor, is the gradient operator in the direction parallel to face f, and is the distance from the cell center to the ghost point O parallel to the cell face f (see Figure 3.2).
Figure 3.2. Schematic showing two types of refined cells.
By definition, and . Note that for neighboring cells without any refinement, is equal to zero and the above equation is consistent with nonrefined cell faces. The linear interpolation factor is defined as
 (38)
where is the coordinate of f perpendicular to the face and is the cell dimension perpendicular to the face f.
Advection Schemes
Hybrid Scheme
The hybrid scheme is a composed of a firstorder upwind scheme and a secondorder central difference scheme. When the Peclet number, , is larger than 2, the firstorder upwind scheme is used; otherwise, the secondorder central difference scheme is used:
 (39)
where the subscripts D and C indicate the first downstream and first upstream nodes and is the Peclet number at the cell face in which .
Exponential Scheme
The exponential scheme interpolates the face value using an exact solution to the 1D steady advectiondiffusion equation:
 (310)
where . The exponential scheme has automatic upwinding and is stable, but is usually less than second order.
Hybrid Linear/Parabolic Scheme
The Hybrid Linear/Parabolic Approximation scheme of Zhu (1991) may be written as
 (311)
where the subscripts D, C and indicate the first downstream and first and second upstream cells, respectively. The normalized variable, , is determined based on the formulation of Jasak et al. (1999)
 (312)
where . The Hybrid Linear/Parabolic Approximation scheme is second order.
Cellface gradient operator
A linearly exact secondorder approximation for the normal gradient at cell face f is calculated using the auxiliary node concept of Ferziger and Peric (1997)
 (313)
where the subscripts P and N refer to two neighboring cells, is the distance between cells P and N, normal to the cell face (see Figure 3.2), and is the gradient operator in the direction parallel to face f. Ham et al. (2002) compared the auxiliary node formulation to the fullyunstructured discretization proposed by Zwart et al. (1998) for the viscous terms and found that the auxiliary node formulation is significantly more stable.
Cellcentered gradient operator
The cellcentered gradient operator is calculated using the GreenGauss Theorem as
 (314)
Equation 314 above is second order and conservative for regular and nonuniform grids.
Reconstruction, Monotonicity, and Slope Limiters
Given a conservative average value within a cell, a linear reconstruction of the variable is used within that cell. This can be expressed as
 (315)
The reconstruction is conservative in the sense that . If the reconstruction satisfies the local maximum principle
 (316)
then no new extrema are created within the cell and the solution is monotonic. Figure x shows two examples of linear reconstruction with and without slope limiters to ensure monotonicity.
then no new extrema are created within the cell and the solution is monotonic. Figure x shows two examples of linear reconstruction with and without slope limiters to ensure monotonicity.
Figure 3.3. Schematics showing examples of nonlimited (a) and limited (b) linear reconstructions.
For nontelescoping grids, the slope limiter is calculated as
where r is the ratio between two consecutive slopes
 (318)
Here the secondorder van Leer (1979) limiter is used because of its smoothness. Figure X shows a comparison of three different common limiters. The slope limiter is applied in each direction separately.
Figure 3.4. Comparison of three different slope limiters.
For joint cells, the standard slope limiters described above are difficult to implement because of the difficulty in defining forward and backward differences. Therefore for joint cells, a variation of the Limited Central Difference (LCD) slope limiting procedure of Hubbard (1999)
where . In the procedure outlined by Hubbard (1999), a scalar limiter is then calculated as . For telescoping grids a directional limiter can be calculated as , which is less dissipative. Look at this again
Source/sink term
The source/sink term is linearized as (Patankar 1980)
 (320)
where is the cell area, and is approximated as the cellaverage source/sink term. The coefficient is required to be nonpositive for stability.
Assembly of Algebraic Equations
Assembly refers to the process of combining terms to create a linear system of algebraic equations. The algebraic equation for each cell is obtained by first combining or assembling all of the terms. Then, the continuity equation is multiplied by and is subtracted from the transport equation. The resulting discretized equation for cell P is
 (321)
where the subscript N refers to the neighboring cell sharing cell face, and are linear coefficients for and . The last term, , contains all the remaining terms (such as…name the specific terms here). Applying a similar equation for all of the internal cells of a grid results in a system of algebraic equations. This set of equations are referred to as the discretized governing equations.
Implicit Relaxation
Underrelaxation stabilizes the convergence of the outer nonlinear iteration loop by introducing a relaxation parameter in the discretized equations (Patankar 1980) as
 (322)
where is an underrelaxation parameter and is the value of from the previous iteration. An effect of underrelaxation is to make the coefficient matrix more diagonally dominant.
Iterative Solvers
The selection of an iterative solver is one key issue concerning the overall performance of the model. CMS has three iteration solvers available and they are described in more detail below: 1) GMRES variation, 2) BiCGStab, and 3) GaussSeidel. The default iterative solver for CMS is a variation of the GMRES (\underline{G}eneralized \underline{M}inimum \underline{RES}idual) method (Saad 1993) and is used to solve the algebraic equations. The original GMRES method (Saad and Schultz 1986) utilizes the Arnoldi process to reduce the coefficient matrix to the Hessenburg form and minimizes the norm of the residual vector over a Krylov subspace at each iterative step. The variation of the GMRES method recommended by Saad (1993) allows changes in preconditioning at every iteration step. An Incomplete LU Factorization (ILUT; Saad, 1994) is used as the preconditioner to speedup convergence. The GMRES solver is applicable to symmetric and nonsymmetric matrices and leads to the smallest residual for a fixed number of iterations. However, the memory requirements and computational costs become increasingly expensive for larger systems.
The BiCGStab (BiConjugate Gradient Stabilized) iterative solver is also a Krylov subspace solver and is applicable to symmetric and nonsymmetrix matrices (Saad 1996). BiCGStab also uses ILUT as a preconditioner (Saad 1994). The BiCGStab method can be viewed as a combination of the standard Biconjugate Gradient solver where each iterative step is followed by a restarted GMRES iterative step. One advantage of the BiCGStab iterative solver is that the memory requirements are constant for each iteration and there are less computational costs when compared to the GMRES method (for GMRES restart numbers larger than 4).
The simplest iterative solver implemented in CMS is the pointimplict GaussSeidel solver. This method may be applied in CMS with or without SuccesiveOverRelaxation to speedup convergence (Patankar 1980). Even though the GaussSeidel method requires more iterations for convergence, the overall efficiency may be higher than the GMRES and BiCGStab because each iteration is computationally inexpensive and the code is parallelized. However, the GMRES and BiCGStab methods are more robust and perform better for large time steps.
Convergence and TimeStepping
During the iterative solution process, error is calculated and used to determine if a solution has converged, diverged, or stalled at an error below a set tolerance threshold. Anestimate of the error in solving the general algebraic equation is given by
 (323)
Statistics (such as…) can be defined based on normalized errors. For example, the norm is given by
 (324)
Since this value depends on the total number of cells, the final statistic (referred to as the residual) that is used for estimating the model convergence is obtained by dividing the norm by the number of cells
 (325)
is referred to as the “normalized residual error” and the superscript refers to the iteration number. is calculated for each variable that is solved at each iteration step of the solution process. Each equation has default maximum tolerances for determining if the solution has converged, diverged, or stalled. The maximum number of iterations that is imposed is set equal to M. A minimum of 5 iterations are required for the hydrodynamic equations, and a minimum of M/2 iterations are required for the sediment transport equations. Table 31 lists the default criteria to determine whether the iterative solution procedure has converged, diverged, or requires a reduced time step.
Table 31. Default criteria to determine whether the iterative solution procedure has converged, diverged, or requires a reduced time step.
Current velocity, m/s 
If R^{m} < 1x10^{7} or  R^{m}R^{m2} < 1x10^{7} 
If R^{M} > 1.0x10^{2} or U_{i} > 10 
If R^{M} > 1.0x10^{3} 
Pressure correction, m^{2}/s^{2} 
If R^{M} < 1x10^{8} or R^{M}  R^{M2} < 1x10^{8} 
If R^{M} > 1.0x10^{3} or p >50 
If R^{M}>1.0x10^{4} 
Totalload concentration, kg/m^{} 
If R^{M}<1x10^{8} or R^{m}R^{m2}<1x10^{8} 
If R^{m}>1.0x10^{3} or C_{tk}<0 
None 
Salinity, ppt  If R^{m}<1x10^{6}  If S < 0  None 
For the implicit model, the time steps for the hydrodynamics, sediment and salinity transport are the same in order to avoid mass conservation problems and for simplicity. If any of the time step reduction criteria are met, then the time step is reduced by half and a minimum number of 3 time steps are calculated at the newly reduced time step. If the last time step converged properly, then the time step is increased. The maximum time step allowed is equal to the userspecified initial time step.
Ramp Period
For most coastal applications, the model is initialized from a “cold start”, which means that the water level and current velocities are initially set to zero. The ramp period allows the model to slowly transition from the initial condition without “shocking” the system. In CMS, the ramp function is defined as
 (326)
where t is the simulation time and is the ramp duration. The ramp function provides a smooth function for transitioning from the initial condition and is plotted in Figure 3.5.
Figure 3.5. Ramp function used in CMS
The ramp function is applied to the model forcing conditions, including the wave forcing , surface wind , sediment concentration capacity , and significant wave height , by direct multiplication of these parameters by the ramp function at each time step during the ramp period. Boundary conditions in CMS are specified without consideration of this ramp period; therefore, the boundary conditions are also slowly transitioned from the initial condition by direct multiplication of the boundary conditions by the ramp function at each time step during the ramp period.
Spatial Discretization
Temporal Discretization
The temporal term of the momentum equations is discretized using a first order implicit Euler scheme

(9) 
where is the cell area, and is the hydrodynamic time step.
Cellface interpolation operator
The general formula for estimating the cellface value of is given by

(12) 
where is a linear interpolation factor given by and is the gradient operator in the direction parallel to face f. By definition . Note that for neighboring cells without any refinement and are zero and thus the above equation is consistent with nonrefined cell faces.
Advection Schemes
The advection scheme obtained using the divergence theorem as where is the outward unit normal on cell face f, is the cell face length and is the total water depth linearly interpolated to the cell face. Here the overbar indicates a cell face interpolation operator described in the following section. For Cartesian grids the cell face unit vector is always aligned with one of the Cartesian coordinates which simplifies the calculation. Defining the cell face normal velocity as the above equation simplifies to

(10) 
where is the outward unit normal on cell face f, is the cell face length and is the total water depth linearly interpolated to the cell face. Here the overbar indicates a cell face interpolation operator described in the following section. For Cartesian grids the cell face unit vector is always aligned with one of the Cartesian coordinates which simplifies the calculation. Defining the cell face normal velocity as the above equation simplifies to

(11) 
where , , with being the basis vector. is equal to 1 for West and South faces and equal to 1 for North and East cell faces. Lastly, is the advective value of on cell face f, and is calculated using either the Hybrid, Exponential, HLPA (Zhu 1991) schemes. The cell face velocities are calculated using the momentum interpolation method of Rhie and Chow (1983) described in the subsequent section. The advection value is calculated as , where the superscripts and indicate low and high order approximations and the superscripts and indicate either explicit and implicit treatment. The explicit term is solved directly while the implicit term is implemented through a deferred correction in which the terms are approximated using the values from the previous iteration step.
Diffusion term
The diffusion term is discretized in general form using the divergence theorem

(13) 
The discritization of the cellface gradient is described in the next section. On a Cartesian grid the above expression may be further simplified as

(14) 
where is gradient in the direction perpendicular to the cell face and .