Convergence-enhancing measures

As explained in Section 3.7.1, many time scales are involved in the evolution of wind waves. The high-frequency waves have much shorter time scales than the low-frequency waves, rendering the system of equations (3.35) stiff. If no special measures are taken, the need to resolve high-frequency waves at very short time scales would result in extreme computational time. For economy, it is desirable that a numerical technique can be used with a large, fixed time step. Moreover, we are mainly interested in the evolution of slowly changing low-frequency waves. For stationary problems, we are interested in obtaining the steady-state solution. Unfortunately, the convergence to the steady state is dominated by the smallest time scale and, in the absence of remedial measures, destabilizing over- and undershoots will prevent solution from converging monotonically during the iteration process. These oscillations arise because of the off-diagonal terms in matrix $A$, which can be dominant over the main diagonal, particularly when the ratio $\sigma_{\rm max}/\sigma_{\rm min}$ is substantially larger than one. As a consequence, convergence is slowed down and divergence often occurs. To accelerate the iteration process without generating instabilities, appropriately small updates must be made to the level of action density.


With the development of the WAM model, a so-called action density limiter was introduced as a remedy to the abovementioned problem. This action limiter restricts the net growth or decay of action density to a maximum change at each geographic grid point and spectral bin per time step. This maximum change corresponds to a fraction of the omni-directional Phillips equilibrium level (Hersbach and Janssen, 1999). In the context of SWAN (Booij et al., 1999), this is

  $$\Delta N \equiv \gamma \, \frac{\alpha_{\rm PM}}{2 \sigma k^3 c_g} \, ,$$ (3.42)

where $\gamma \geq 0$ denotes the limitation factor, $k$ is the wave number and $\alpha_{\rm PM} = 8.1 \times 10^{-3}$ is the Phillips constant for a Pierson-Moskowitz spectrum (Komen et al., 1994). Usually, $\gamma = 0.1$ (Tolman, 1992)^3.5. Note that when the physical wind formulation of Janssen (1989,1991a) is applied in SWAN, the original limiter of Hersbach and Janssen (1999) is employed. Denoting the total change in $N_{i,j,l,m}$ from one iteration to the next after Eq. (3.2) by $\Delta N_{i,j,l,m}$, the action density at the new iteration level is given by

  $$N_{i,j,l,m}^s = N_{i,j,l,m}^{s-1} + \frac{\Delta N_{i,j,l,m}}{\ve... ...N_{i,j,l,m}\vert} \, \min \{ \vert \Delta N_{i,j,l,m} \vert, \Delta N \} \, .$$ (3.43)

For wave components at relatively low frequencies, Eq. (3.43) yields the pre-limitation outcome of Eq. (3.2), because, for these components, the pseudo time step matches the time scale of their evolution. For high-frequency waves, however, Eq. (3.43) gives the upper limit for the spectrum to change per iteration due to the limiter, Eq. (3.42). For typical coastal engineering applications, it is sufficient to compute the energy-containing part of the wave spectrum accurately. In other words, action densities near and below the spectral peak should not be imposed by the limiter (3.42). However, our experiences with SWAN have shown that the limiter is active even close to the peak. Furthermore, during the entire iteration process, the limiter is typically active at almost every geographic grid point.


The alternative measure to enhance the convergence of the stable iteration process considered here is so-called false time stepping (Ferziger and Perić, 1999). Under-relaxation terms representing the rate of change are introduced to enhance the main diagonal of $A$ and thus stabilize the iteration process. The system of equations (3.35) is replaced by the following, iteration-dependent system

  $$\frac{{\vec{N}}^s - {\vec{N}}^{s-1}}{\tau} + A\,{\vec{N}}^s = \vec{b}$$ (3.44)

with $\tau$ a pseudo time step. The first term of Eq. (3.44) controls the rate of convergence of the iteration process in the sense that smaller updates are made due to decreasing $\tau$, usually at the cost of increased computational time. To deal with decreasing time scales at increasing wave frequency, the amount of under-relaxation is enlarged in proportion to frequency. This allows a decrease in the computational cost of under-relaxation, because at lower frequencies larger updates are made. This frequency-dependent under-relaxation can be achieved by setting ${\tau}^{-1} = \alpha \sigma$, where $\alpha$ is a dimensionless parameter. The parameter $\alpha$ will play an important role in determining the convergence rate and stability of the iteration process. Substitution in Eq. (3.44) gives

  $$(A + \alpha \sigma I)\, {\vec{N}}^s = \vec{b} + \alpha \sigma {\vec{N}}^{s-1} \, .$$ (3.45)

When the steady state is reached (i.e. $s \rightarrow \infty$), system (3.45) solves $A\,{\vec{N}}^{\infty} = \vec{b}$ since, ${\vec{N}}^{\infty}$ is a fixed point of (3.45).


Suitable values for $\alpha$ must be determined empirically and thus robustness is impaired. For increasing values of $\alpha$, the change in action density per iteration will decrease in the whole spectrum. The consequence of this is twofold. Firstly, it allows a much broader frequency range in which the action balance equation (3.2) is actually solved without distorting convergence properties. Secondly, the use of the limiter will be reduced because more density changes will not exceed the maximum change due to Eq. (3.42). Clearly, this effect may be augmented by increasing the value of $\gamma$ in Eq. (3.42).


To allow proper calculation of the second-generation first guess of the wave field (see Section 3.3), under-relaxation is temporarily disabled ($\alpha = 0$) during the first iteration. Whereas this measure is important in achieving fast convergence, it does not affect stability, since the second-generation formulations do not require stabilization.