Solution algorithm

The implicit discretization of the action balance equation (2.19) as described in Section 3.2 yields a system of linear equations that need to be solved. The corresponding matrix structure can take different forms, mainly depending on the propagation of wave energy in the geographic space. For instance, suppose that $c_x > 0$ and $c_y > 0$, everywhere. Then, the matrix structure has the following form:

  $$\left[ \begin{array}{l} \left( \begin{array}{l} \mathrm{.\, x... ... x\, .\, x\, x\, x\, .\, x\,} \end{array} \right)\\ \end{array} \right]\, .$$ (3.31)

One recognizes that the subblocks on the main diagonal express coupling among the unknowns in the $(\sigma,\theta)-$space for each geographic grid point, whereas the off-diagonal subblocks represent coupling across geographical grid points. The ordering of grid points is determined by the direction of wave propagation. This system can be solved with a Gauss-Seidel technique in one step, if the wave characteristic is a straight line and constant everywhere (Wesseling, 1992). In addition, this number is independent of grid size. Hence, the complexity of this algorithm is ${\cal O} (M)$ for a total of $M$ grid points.


After each propagation update at geographic grid point, an update in the spectral space must be made. We divide all wave directions into a number of groups according to their directions, and order the grid points accordingly for the update. This is called a sweep. Since every internal grid point has the same number of edges with fixed directions, same division can be made within the $(\sigma,\theta)$-space, resulting in four quadrants of each $90^o$, as illustrated in Figure 3.1. In this case we have four sweeps and the selected wave directions in each sweep form the domain of dependence appropriate for update. This immediately satisfies the CFL criterion. This criterion is related to the causality principle. In general, a numerical scheme must look for information by following characteristics in an upwind fashion. Causality can be preserved in the iteration process by means of ordering grid points according to the propagation direction, which guarantees convergence in a finite number of iterations. We will come back to that later in Section 3.4.

Figure 3.1: The solution procedure for wave energy propagation in geographical space with the appropriate directional quadrant (indicated by shaded area) for each of four sweeps.

The Gauss-Seidel iteration process is done as follows. For each iteration, sweeping through grid rows and columns in geographical domain are carried out, starting from each of the four corners of the computational grid. After four sweeps, wave energy has been propagated over the entire geographical domain. During each sweep, only a subset of the unknown values of $N$ are updated depending on the sign of $c_x$ and $c_y$. For instance, the first sweep starts at the lower left hand corner and all grid points with $c_x > 0$ and $c_y > 0$ are updated. Because of the causality principle these transport velocities must be positive in those ordered grid points along the wave ray, in order to make a stable iterative update. Moreover, adapting the ordering of updates of the unknowns $N$ in geographical space to the propagation direction improves the rate of convergence of the Gauss-Seidel iterative procedure (Wesseling, 1992). For an illustrative explanation of this technique, see Section 3.5.


Due to the implicit nature of the spectral propagation terms in Eq. (3.2), a system of equations must be formed (i.e. one of the main diagonal of the matrix Eq. (3.31)). Furthermore, due to the fact that the source term $S_{\rm tot}$ in Eq. (3.2) is nonlinear in $N$, linearization is required in order to find a solution. Generally, the term $S_{\rm tot}$ in each bin $(l,m)$ is treated by distinguishing between positive and negative contributions and arranging these in the linear form (Ferziger and Perić, 1999):

  $$S_{\rm tot} = S_{\rm tot}^p + S_{\rm tot}^n N \, ,$$ (3.32)

where $S_{\rm tot}^p$ consists of positive contributions and $S_{\rm tot}^n$ of negative ones. Both contributions are independent of the solution $N$ at the corresponding bin $(l,m)$. Any negative term that does not contain $N$ as a multiplier is first divided by $N$ obtained from the previous iteration level and then added to $S_{\rm tot}^n$. This stabilizes the iteration process. Details on the application of this principle to each source term in SWAN can be found in Booij et al. (1999).


The strongly nonlinear source term of depth-induced wave breaking is linearized by means of the Newton-Raphson iteration, as follows

  $$S^n \approx \phi^{n-1} E^{n} + \left( \frac{\partial S}{\partial E} \right)^{n-1} (E^n - E^{n-1})$$ (3.33)

Since this process of depth-induced wave breaking has been formulated such that $S = aS_{\rm tot}$ and $E = aE_{\rm tot}$, the derivative $\partial S/\partial E$ is analytically determined as $\partial S_{\rm tot}/\partial E_{\rm tot}$. Here, $a$ is identical in both expressions and the total energy $E_{\rm tot}$ and total source $S_{\rm tot}$ are the integrals over all frequencies and directions of $E(\sigma,\theta)$ and $S_{\rm {ds,br}}(\sigma,\theta)$, respectively.


As such, each difference equation (3.2) using expressions (3.20), (3.21) and (3.32) provides an algebraic relation between $N$ at the corresponding bin and its nearest neighbours:

  $$a_{\rm P} N_{\rm P} = a_{\rm L} N_{\rm L}+ a_{\rm R} N_{\rm R}+ a_{\rm B} N_{\rm B}+ a_{\rm T} N_{\rm T}+ b_{\rm P} \, ,$$ (3.34)

where P corresponds to central bin $(l,m)$ and L(eft), R(ight), B(ottom) and T(op) correspond to $(l-1,m)$, $(l+1,m)$, $(l,m-1)$ and $(l,m+1)$, respectively. Furthermore, the coefficients $a_k$, $k \in \{ \rm P,\rm L, \rm R, \rm B, \rm T \}$ arise from the discretizations of the fluxes $c_{\sigma} N$ and $c_{\theta} N$ and $b_{\mathrm{P}}$ contains the positive contributions of the source term $S_{\rm tot}^p$ in (3.32) and the updated fluxes $c_x N$ (3.3) and $c_y N$ (3.4). Note that coefficient $a_{\rm P}$ includes $-S_{\mathrm{tot}}^n$.


The linear system of equations (3.34) for all bins within a directional quadrant at a particular geographical point is denoted by

  $$A \, \vec{N} = \vec{b} \, ,$$ (3.35)

where $A \in {\mbox{$I\!\!R$}}^{K \times K}$ contains the coefficients $a_k$, $k \in \{ \mathrm{P}, \mathrm{L}, \mathrm{R}, \mathrm{B}, \mathrm{T} \}$ (and corresponds to a subblock on the main diagonal of (3.31)), $\vec{b} \in {\mbox{$I\!\!R$}}^{K}$ contains the coefficient $b_{\rm P}$ and boundary values and $\vec{N} \in {\mbox{$I\!\!R$}}^{K}$ denotes an algebraic vector containing the unknown action density values. Matrix $A$ is non-symmetric. The dimension $K$ of a directional quadrant equals $N_{\sigma} \times 1/4 N_{\theta}$. Note that linearization of the source term (3.32) enhances diagonal dominance of $A$, thereby improving numerical stability. Also note that neither $A$ nor $\vec{b}$ depends on the unknowns. Each row in the matrix $A$ corresponds to a bin $(l,m)$. The main diagonal contains the coefficients $a_{\rm P}$ and directly to the left and right are the coefficients $-a_{\rm B}$ and $-a_{\mathrm{T}}$, respectively. The coefficients $-a_{\rm L}$ and $-a_{\rm R}$ are on the diagonals that are $N_{\theta}$ positions to the left and right of the main diagonal, respectively.


The solution $\vec{N}$ is given by $A^{-1} \vec{b}$. Since, the only non-zero matrix elements are situated in five diagonals, iterative solution methods that utilize the sparsity of $A$ optimally are very attractive. In SWAN, the solution of Eq. (3.35) is found by means of an incomplete lower-upper decomposition method followed by an iteration process called the Strongly Implicit Procedure (SIP) (Ferziger and Perić, 1999). This procedure is specifically designed for (non-symmetric) penta-diagonal systems and is relatively fast. Note that in the absence of mean current there are no shifts in the frequency, and consequently the structure of $A$ reduces to a tri-diagonal one, i.e. $a_{\rm L} = a_{\rm R} = 0$, which can be inverted efficiently with the Thomas algorithm (Press et al., 1993; Ferziger and Perić, 1999).