next up previous index
Next: Output Up: Model description Previous: Physics   Index

Numerics


\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

NONHYDrostatic  /    STAndard \  [theta]                                   &
                \ -> BOX      /

                SUBGrid [pmax]  REDuced [qlay]                             &

                SOLVer [rhsaccur] [initaccur] [maxiter] [relax]            &

                PREConditioner  ILUD|ILU

\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

With this optional command the user can include the non-hydrostatic pressure in the shallow water equations. If this command is not used, SWASH will not account for non-hydrostatic pressure, i.e. pressure is assumed to be hydrostatic.


Hydrostatic pressure assumption can be made in case of propagation of long waves, such as large-scale ocean circulations, tides and storm surges. This assumption does not hold in case of propagation of short waves, flows over a steep bottom, unstable stratified flows, and other small-scale applications where vertical acceleration is dominant.


In SWASH two different schemes for the vertical pressure gradient are available, i.e., the classical central differencing (option STANDARD) and the Keller-box scheme (option BOX). The former approximation is particularly meant for applications where vertical structures are important, e.g. stratified flows with density currents, undertow and flows over steep and rapidly varying bottoms, while the latter will be mainly used for accurate short wave propagation with two or three layers employed. In addition, when applying the Keller-box scheme, it is advised to choose an equidistant layer distribution (variable thicknesses only) in order to get the optimal linear frequency dispersion accuracy (i.e. the relative error of the phase speed predicted by the model compared to that from the linear dispersion theory is minimal).


The time integration of the vertical pressure gradient is the so-called $ \theta$ -scheme (a mix of explicit and implicit Euler schemes). With [theta] = 0.5 we have the well-known second order accurate Crank-Nicolson scheme with the smallest truncation error, while [theta] = 1 indicates the first order implicit Euler scheme. Note that only values of [theta] in the range [0.5,1] are allowed for stability reasons.


The number of vertical layers is chosen sufficiently large to resolve the vertical structure of the flow field (see command VERTICAL). However, this grid resolution in the vertical is also employed in the discretization of the non-hydrostatic pressure field, whereas a few number of layers usually suffices to compute wave dispersion without deteriorating accuracy. With a subgrid approach the user has an option to resolve the flow and pressure field with a different resolution in the vertical. The vertical is equidistantly divided into a few number of layers for the non-hydrostatic pressure and a relative large number of (non-equidistant) layers for the horizontal velocities and turbulent stresses. In this way, the solution for the non-hydrostatic pressure and vertical acceleration can be obtained on a relatively coarse vertical grid, while on a subgrid with a high vertical resolution the vertical structure of turbulent flow is resolved.


Since most of the computational effort is devoted to inverting the Poisson pressure matrix, an effective way to minimize this effort is by reducing the rank of the Poisson matrix. This leads to less computational cost and memory. With the keyword REDUCED the dimension of the Poisson equation, i.e. the number of pressure unknowns per water column, can be diminished. In spite of this reduction it can still provide an accurate description of dispersive waves. For instance, a model with two layers but one reduced pressure layer, the so-called reduced two-layer model, has more or less the same linear dispersion accuracy as the full two-layer model but saves about 30% CPU time. The stencil of the pressure equation can be reduced by assuming a constant pressure in the vertical near the bottom. This is determined by the parameter [qlay] which indicates the number of layers, started from the bottom, where the pressure is constant. In this way, it is possible to eliminate the pressure in these layers from the set of equations. As an example, we have a model with five pressure layers ([pmax] = 5) and we assume that the pressure is constant in the two lowest layers ([qlay] = 2). So there are effectively three (upper) layers in which the pressure need to be computed. An important remark needs to be made. Although this reduced pressure equation method is more efficient in terms of CPU time, it is less so concerning the dispersion accuracy. To optimize this accuracy often a different than equidistant layer distribution needs to be sought. For instance, the optimal linear dispersion relation of the reduced two-layer model ([kmax] = [pmax] = 2, [qlay] = 1) is obtained by setting the following layer distribution: the top layer has a relative thickness of 84% whereas the thickness of the bottom layer is 16% of the total depth.


By inclusion of the non-hydrostatic pressure, a solution of the Poisson pressure equation is required. In SWASH this equation is solved by an iterative solution method and the user may controls this by means of the keyword SOLVER. Two linear solvers are adopted: Strongly Implicit Procedure (SIP) and BiCGSTAB preconditioned with an incomplete LU factorization. The former one is particularly meant for the depth-averaged case, while the latter one can only be applied for the multi-layered case. The incomplete LU factorization is either ILU or ILUD. The latter is restricted to the main diagonal of the matrix, i.e. the off-diagonals remain unchanged. The choice of the preconditioner is indicated with the keyword PRECONDITIONER. For parallel computing, the ILUD preconditioner is a good choice. However, the ILU preconditioner is more robust. For instance, when high waves or very short waves are involved, or when the bottom topography exhibits steep slopes, or when a considerable number of layers (> 20) is involved, it may be wise to choose the ILU preconditioner. The weighting parameter $ \alpha$, as indicated by [relax], may improve the rate of convergence. With this parameter a combination of the classical ILU and the modified ILU (MILU) can be given. This combination is given by (1-$ \alpha$)ILU + $ \alpha$MILU, and is also hold for ILUD and its modified variant (MILUD). Based on several numerical experiments, an optimum in the convergence rate is found by taking 55% of MILU and 45% of ILU in case of the Keller-box scheme, and 99% of MILUD and 1% ILUD for the same scheme, while for the standard discretization of the vertical pressure gradient, a combination of 90% of MILU(D) and 10% of ILU(D) is chosen.


It is common to use the reduction of the residual as a stopping criterion, because the BiCGSTAB method requires calculation of the residual. When solving the system Ax = b, after m iterations we have an approximate solution xm and the residual rm = b - Axm is related to the convergence error em = x - xm by Aem = rm, so the reduction of the residual results in the reduction of the convergence error. This does not necessarily mean that the relative error of the solver is identical to the decrease of the residual. The iteration process stops at each time step if the ratio of the 2-norm of the residual and of the right-hand side or initial residual is less than a given accuracy: | rm|2/| b|2 < $ \epsilon$ and | rm|2/| r0|2 < $ \epsilon$, respectively. They are indicated by the parameters [rhsaccur] and [initaccur], respectively. If both these accuracies are given, the sum of the two is used as termination criterion. Often, the stopping criterion for the iterative methods is basically a compromise between efficiency and accuracy. Decreasing the required accuracy can save a considerable amount of CPU time. Numerical experiments showed $ \epsilon$ = 0.01 gives the optimum.


The default option is the Keller-box scheme with [theta] = 1.0, while the Poisson pressure equation is solved with either SIP with [rhsaccur]=0.01, in the case of depth-averaged mode, or ILUD-BiCGSTAB with [rhsaccur]=0.01 and [initaccur]=0, in the case of multi-layered mode.

STANDARD this option indicates that the classical central differencing for the  
  vertical pressure gradient will be employed.  
BOX this option indicates that the Keller-box scheme will be adopted.  
  This is the default.  
[theta] parameter for the $ \theta$ -method. [theta] = 0.5 corresponds to the  
  Crank-Nicolson scheme and [theta] = 1 to the implicit Euler  
  scheme.  
  Default: [theta] = 1.0  
SUBGRID this option indicates that the subgrid approach will be adopted.  
[pmax] number of vertical pressure layers and 1 $ \leq$ [pmax] $ \leq$ kmax.  
  Default: [pmax] = kmax  
REDUCED this option indicates that the reduced pressure equation method will  
  be adopted.  
[qlay] number of layers, started from the bottom layer, for which pressure is  
  constant in the vertical to reduce the pressure equation.  
  Default: [qlay] = 1  
SOLVER With this option the user can influence some of the numerical parameters  
  of iterative solvers.  
[rhsaccur] relative accuracy with respect to the right-hand side.  
  Default: [rhsaccur] = 0.01  
[initaccur] Relative accuracy with respect to the initial residual.  
  Default: [initaccur] = 0.0  
[maxiter] maximum number of iterations.  
  Default: [maxiter] = 500  
[relax] parameter for a linear solver:  
  $ \bullet$ relaxation parameter for the ILU(D) preconditioning. [relax] = 0  
  corresponds to ILU(D) and [relax] = 1 corresponds to MILU(D).  
  Default: model dependent (see above).  
  $ \bullet$ under relaxation parameter for the SIP solver.  
  Default: [relax] = 0.91  
PRECONDI With this option the user can choose a preconditioner.  
ILUD indicates the ILUD preconditioner.  
  This is the default.  
ILU indicates the ILU preconditioner.  


\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

                |         | UMOM  MOMentum|HEAD  / -> Horizontal   |
                |         |                      \    Vertical     |
                | UPWind <                                         |
                |         | WMOM  / -> Horizontal                  |
DISCRETization <          |       \    Vertical                     >      &
                |                                                  |
                | CORRdep                                          |
                |                                                  |
                | TRANSPort  / -> Horizontal                       |
                |            \    Vertical                         |


                   | NONe                                 |
                   |                                      |
                   | FIRstorder                           |
                   |                                      |
                   | HIGherorder [kappa]                  |
                   |                                      |
                   |          | -> SWEBy [phi]            |
                   |          |                           |
                   | LIMiter <  RKAPpa [kappa]            |
                   |          |                           |
                   |          | PLKAPpa [kappa] [mbound]  |
                   | FROmm                                |
                   |                                      |
                   | -> BDF|LUDs                          |
                   |                                      |
                  <  QUIck                                 >
                   |                                      |
                   | CUI                                  |
                   |                                      |
                   | MINMod                               |
                   |                                      |
                   | SUPerbee                             |
                   |                                      |
                   | VANLeer                              |
                   |                                      |
                   | MUScl                                |
                   |                                      |
                   | KORen                                |
                   |                                      |
                   | SMArt                                |

\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

With this optional command the user can influence the space discretization.


For convection-dominated flows it is possible that wiggles in the solution arise. In that case upwind discretization may be necessary. At this moment, three types of upwind schemes are implemented:

Schemes of up to third order accuracy can be constructed by piecewise polynomial interpolation, the so-called $ \kappa$ -formulation. For all values of $ \kappa$ $ \in$ [- 1, 1], a blended form arises between second order backward difference scheme (BDF) and second order central differencing. The schemes BDF, QUICK and CUI are obtained by setting $ \kappa$ = -1, 1/2 and 1/3, respectively. The value $ \kappa$ = 0 gives the Fromm's scheme, while $ \kappa$ = 1 corresponds to central differencing. For $ \kappa$ $ \neq$1/3 the local truncation error is of second order; for $ \kappa$ =1/3 it is of third order. All upwind schemes are applied in each computational direction. Hence no streamline upwinding is used.


In case of transport equations for salinity, temperature, or suspended sediment (keyword TRANSPORT) a TVD scheme must be applied to prevent wiggles in the solution or to avoid negative concentration values.


The water depth in velocity points is not uniquely defined. An appropriate approximation is based on first order upwinding instead of the usual interpolation. To achieve second order accuracy in space, we add a higher order interpolation augmented with a flux limiter. See keyword CORRDEP.


Conservation properties become crucial for rapidly varied flows. These properties are often sufficient to get solutions that are acceptable in terms of local energy losses, location of incipient wave breaking, propagation speed of a bore, etc. In flow expansions, the horizontal advective terms in u - and v -momentum equations are approximated such that they are consistent with momentum conservation; see option MOMENTUM. In flow contractions, the approximation is such that constant energy head is preserved along a streamline, i.e. the Bernoulli equation; see option HEAD.


The default option is the second order backward difference (BDF) scheme ($ \kappa$ = - 1) for all horizontal advective terms in both u/v -momentum and w -momentum equations, and first order upwinding for the vertical term in the u/v -momentum equation and in the w -momentum equation. The exception is when the BREAK command is employed, which in that case the central differences ($ \kappa$ = 1) are applied to all horizontal advective terms in the u/v -momentum equations. In addition, the water depth in velocity points is approximated with the MUSCL limiter. All the advection terms in any transport equation are approximated with the second order Van Leer limiter. By default, SWASH decides whether energy head or momentum conservation is to be applied. In principle, energy head conservation will be applied only in strong flow contractions, while elsewhere the momentum conservation is applied.

UPWIND indicates the type of discretization for momentum equations.  
UMOM indicates the discretization employed for u/v -momentum equation.  
MOMENTUM indicates that momentum must be conserved everywhere.  
HEAD indicates that energy head must be conserved everywhere.  
HORIZONTAL indicates the discretization of the horizontal advective terms.  
VERTICAL indicates the discretization of the vertical advective terms.  
WMOM indicates the discretization employed for w -momentum equation.  
CORRDEP indicates the type of discretization for water depth in velocity points.  
TRANSPORT indicates the discretization employed for transport equation(s).  
NONE indicates that no upwinding is applied, hence the standard  
  central difference scheme is used.  
  This is the default in case of vertical advective terms.  
FIRSTORDER indicates that the standard first order upwind scheme is used.  
HIGHERORDE a higher order upwind $ \kappa$ -scheme is activated.  
[kappa] parameter defining the type of $ \kappa$ -scheme used.  
  For instance, [kappa] = 0.5 gives the QUICK scheme  
  and [kappa] = 0 the Fromm's scheme.  
  Note: -1 $ \leq$ [kappa] $ \leq$ 1.  
LIMITER indicates that a TVD scheme with flux limiter must be applied.  
SWEBY a Sweby limiter is activated.  
[phi] parameter defining the type of Sweby limiter used. For instance,  
  [phi] = 1 defines the minmod limiter.  
  Note: 1 $ \leq$ [phi] $ \leq$ 2.  
RKAPPA an R- $ \kappa$ limiter is activated.  
PLKAPPA a PL- $ \kappa$ limiter is activated.  
[mbound] parameter indicating the upper bound of the PL- $ \kappa$ limiter.  
  Note: [mbound] > 0.  
FROMM indicates the Fromm's scheme.  
BDF indicates the second order backward upwind scheme.  
  This is the default in case of horizontal advective terms.  
QUICK indicates the QUICK scheme.  
CUI indicates the CUI scheme.  
MINMOD indicates the minmod limiter.  
SUPERBEE indicates the superbee limiter.  
VANLEER indicates the Van Leer limiter.  
MUSCL indicates the MUSCL limiter.  
  This is the default in case of CORRDEP.  
KOREN indicates the Koren limiter.  
SMART indicates the SMART limiter.  


\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

           |    MIN
           |
           | -> MEAN
BOTCel    <
           |    MAX
           |
           |    SHIFt

\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

With this optional command the user can determine how the bottom levels need to be computed in cell centers.


The bottom levels are usually defined in the corners of computational cells. These bottom levels have been extracted from the input grid (see command INPGRID BOTTOM). Nevertheless, SWASH uses a staggered grid and to determine the total water depth at water level points, a bottom level in the cell center is then required. As a rule, we take the average bottom level from the surrounding bottom corners to determine the total depth at cell center. This is fine for many cases. However, in the vicinity of steep bottom slopes, use of the average bottom level may lead to an inaccurate flooding and drying process (e.g. too much or too less volume of water left on tidal flats or flood plains) or inaccurate tide propagation. Another example is the inaccurate computation of pressure in front of the vertical wall which may negatively affect wave overtopping. For those situations it is more appropriate to take the bottom level at cell center being equal to the minimum of the surrounding bottom levels at the corner points. This so-called tiled bottom approach is also suitable to represent a small channel of one grid cell width.


Alternatively, it is possible to specify the bottom level at cell center as the maximum of the surrounding bottom levels or being shifted from the upper-right corner of the computational cell. This latter implies that, if the bottom input grid is identical to the computational grid in terms of resolution and orientation, the user-defined bottom values are specified on input at cell centers (instead of upper-right corners).


Note that for determining the bottom level in cell centers the positive downwards orientation is considered. So the minimum operation results into the shallowest bottom level. As a consequence, when the bottom level is above the reference level, i.e. a negative value, the maximum operation should be chosen, instead of the minimum one, in the context of the tiled bottom approach.


The default option is MEAN.

MIN the bottom level in cell center is the minimum of the four surrounding  
  bottom corner points.  
MEAN the bottom level in cell center is the average of the four surrounding  
  bottom corner points.  
MAX the bottom level in cell center is the maximum of the four surrounding  
  bottom corner points.  
SHIFT the bottom level in upper-right corner is shifted to the cell center.  


\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

            | -> EXPL [cfllow] [cflhig]                                 |
TIMEI METH <                                                             > &
            |    IMPL [thetac] [thetas] SOLVer [tol] [maxiter] [weight] |

      VERTical [thetau] [thetaw] [thetat]

\begin{picture}(18,0.12)
\thicklines
\put(0,0){\line(1,0){17}}
\put(0,0.1){\line(1,0){17}}
\end{picture}

With this optional command the user can influence the time integration.


If time integration is explicit, a time step restriction must be applied based on a Courant number associated with the long wave speed. For definition, see Appendix A. Note that a maximum Courant number of 0.5 is advised in case of high waves, nonlinearities (e.g. wave breaking, wave-wave interactions), and wave interaction with structures with steep slopes (e.g. quays, piers).


An automatic time step control is implemented as follows. The actual maximum of the Courant number over all wet grid points is determined. The time step is halved when this number becomes larger than a preset constant [cflhig] < 1, and the time step is doubled when this number is smaller than another constant [cfllow], which is small enough to be sure the time step can be doubled. Information on the actual time step changes is provided in the PRINT file.


If time integration is semi-implicit, then the gradient of the water level in the momentum equations and the velocity divergence in the (global) continuity equation are discretized implicitly by means of the $ \theta$ -method, while both the horizontal advective and viscosity terms are discretized explicitly. As a consequence, the stability of the method will not depend upon the long wave speed. However, the time step will be restricted owing to the explicit treatment of the horizontal advective terms, although this restriction is mild. This method is particular useful and efficient for the simulation of three-dimensional circulation driven by buoyancy, tides and winds, as the Courant number can easily be larger than 1 while still providing sufficiently accurate solution.


This semi-implicit time stepping requires the solution of a system of linear equations to obtain the water levels. This system is symmetric and positive definite and can be solved efficiently by using a preconditioned conjugate gradient method. The keyword SOLVER controls the use of this iterative method. The iteration process stops if the norm of the residual falls below a small fraction of the initial residual; this small fraction is the user prescribed error tolerance [tol]. The preconditioner is a weighted combination of ILUD and its modified variant MILUD. For details on the weighting parameter $ \alpha$, which may improve the convergence, see command NONHYDROSTATIC. This parameter is indicated by [weight].


In the multi-layered mode, the vertical advective and viscosity terms in the momentum and transport equations are discretized implicitly by means of the $ \theta$ -method.

METHOD with this option the type of time integration is indicated.  
EXPLICIT indicates that the time integration is explicit.  
  This is the default.  
[cfllow] minimum Courant number to be used for automatic time step control.  
  Default: [cfllow] = 0.4  
[cflhig] maximum Courant number to be used for automatic time step control.  
  Default: [cflhig] = 0.8  
IMPLICIT indicates that the time integration is semi-implicit.  
[thetac] parameter for the $ \theta$ -method applied to the continuity equation.  
  [thetac] = 0.5 corresponds to the Crank-Nicolson scheme and  
  [thetac] = 1 to the implicit Euler scheme.  
  Default: [thetac] = 0.5  
[thetas] parameter for the $ \theta$ -method applied to the water level gradient.  
  [thetas] = 0.5 corresponds to the Crank-Nicolson scheme and  
  [thetas] = 1 to the implicit Euler scheme.  
  Default: [thetas] = 0.5  
SOLVER With this option the user can influence some of the parameters of  
  the preconditioned conjugate gradient method.  
[tol] error tolerance to terminate the iteration process.  
  Default: [tol] = 0.01  
[maxiter] maximum number of iterations.  
  Default: [maxiter] = 100  
[weight] weighting parameter for the (M)ILUD preconditioning. [weight] = 0  
  corresponds to ILUD and [weight] = 1 corresponds to MILUD.  
  Default: [weight] = 0.9  
VERTICAL with this option the time integration of vertical terms is indicated.  
[thetau] parameter for the $ \theta$ -method in case of the u/v -momentum equation.  
  [thetau] = 0.5 corresponds to the Crank-Nicolson scheme and  
  [thetau] = 1 to the implicit Euler scheme.  
  Default: [thetau] = 0.5  
[thetaw] parameter for the $ \theta$ -method in case of the w -momentum equation.  
  [thetaw] = 0.5 corresponds to the Crank-Nicolson scheme and  
  [thetaw] = 1 to the implicit Euler scheme.  
  Default: [thetaw] = 0.5  
[thetat] parameter for the $ \theta$ -method in case of the transport equation.  
  [thetat] = 0.5 corresponds to the Crank-Nicolson scheme and  
  [thetat] = 1 to the implicit Euler scheme.  
  Default: [thetat] = 0.5  


next up previous index
Next: Output Up: Model description Previous: Physics   Index
The SWASH team 2017-04-06