The Fast Discrete Interaction Approximation Concept

: Hasselmann and coauthors proposed the discrete interaction approximation (DIA) as the best tool replacing the nonlinear evolution term in a numerical wind–wave model. Much later, Polnikov and Farina radically improved the original DIA by means of location all the interacting four wave vectors, used in the DIA conﬁguration, exactly at the nodes of the numerical frequency–angular grid. This provides a nearly two-times enhancement of the speed of numerical calculation for the nonlinear evolution term in a wind–wave model. For this reason, the proposed version of the DIA was called as the fast DIA (FDIA). In this paper, we demonstrate all details of the FDIA concept for several frequency–angular numerical grids of high-resolution with the aim of active implementation of the FDIA in modern versions of world-wide used wind–wave models.


Introduction
Nonlinear interactions between waves play a very important role in description of wind-wave evolution governed by the equation [1] where D/Dt is the total derivative operator, N k ≡ N(k, x, t) is the wave-action spectrum in the wave vector k-space, at location x, and time t; IN, NL, DISS are the evolution terms responsible for the input, conservative nonlinear transfer among wave components, and dissipation of wave action, respectively. The nonlinear evolution term NL is described by the six-fold Hasselmann kinetic integral I NL with a very complicated integrand [2] ∂N(k 4 ) ∂t where M 2 1,2,3,4 ≡ M 2 (k 0 , k 1 , k 2 , k 3 ) is the second power of the matrix elements corresponding to the four-wave nonlinear interactions, δ 1+2−3−4 ≡δ(σ 1 + σ 2 − σ 3 − σ 4 )δ(k 1 + k 2 − k 3 − k 4 ) is the Dirac delta-function responsible for the resonant feature of the four-wave interactions, and σ i = σ(k i ) is the radian frequency of the wave component with wave vector k i . Due to this complicity, the calculation of integral I NL takes too much time; therefore, to be used in a practical numerical wind-wave model, the integral should be replaced by some theoretically justified approximation. The best approximation was proposed by Hasselmann et al. [3], based on replacing integral I NL by the only configuration of four interacting waves, located at a singular subsurface. This subsurface in the six-fold k-space is defined by the resonance conditions Fluids 2020, 5, 176 2 of 12 The wave vectors k i are usually represented in the frequency-angular space, (σ, θ), where the wave frequencies σ i are related to k i by the dispersion relation, in the deep-water case having the kind σ(k i ) = σ i = (gk i ) 1/2 . (3) The proposed approximation is named as the discrete interaction approximation (DIA). Example of the four vectors configuration (which is usually called as a quadruplet) used in the original DIA [3] is schematically shown in Figure 1. (3) The proposed approximation is named as the discrete interaction approximation (DIA). Example of the four vectors configuration (which is usually called as a quadruplet) used in the original DIA [3] is schematically shown in Figure 1. It is worthwhile to mention that there are a lot of papers devoted to optimization the performance of the NL-term calculations (e.g., [4][5][6][7][8][9][10][11], see also a list of special workshops on wind wave evolution [11]). Some certain approaches to simplify the NL-term calculations, different from the DIA (except the two-scale approximation (TSA) [5,6] derived later), are presented in our paper [12], where a radical preference of the DIA was shown. This is why we dwell on the features of the DIA.
Regarding the other methods, we would say the following. Tamura et al. [4] studied the impact of non-linear energy transfer on realistic wave fields of the Pacific Ocean using the Simplified Research Institute of Applied Mechanics (SRIAM) model, which was developed to accurately reproduce non-linear source terms with lower computational cost than more rigorous algorithms, and the widely used the DIA method. Comparison of the model with buoy observations revealed a negligible difference in significant wave heights but pronounced bias in peak frequency with the DIA. The analysis of spectral shape indicated that the SRIAM method can quantitatively capture the overshoot phenomena around the spectral peak during wave growth The TSA method [5] resides in shearing the whole wave spectrum under the kinetic integral in two parts: a high frequency and the energy containing ones. This allows excluding a reasonable value of the quadruplets from the NL-term calculations, saving accuracy. The TSA has recently been presented as a new method to estimate nonlinear transfer rates in wind waves, and has been tested for idealized spectral data, as well as for observed field measurements [6]. The TSA has been implemented in the wind-wave model WAVEWATCH III and shown to perform well for wave spectra from field measurements, even for cases with directional energy shearing, compared to the DIA. It is worthwhile to mention that there are a lot of papers devoted to optimization the performance of the NL-term calculations (e.g., [4][5][6][7][8][9][10][11], see also a list of special workshops on wind wave evolution [11]). Some certain approaches to simplify the NL-term calculations, different from the DIA (except the two-scale approximation (TSA) [5,6] derived later), are presented in our paper [12], where a radical preference of the DIA was shown. This is why we dwell on the features of the DIA.
Regarding the other methods, we would say the following. Tamura et al. [4] studied the impact of non-linear energy transfer on realistic wave fields of the Pacific Ocean using the Simplified Research Institute of Applied Mechanics (SRIAM) model, which was developed to accurately reproduce non-linear source terms with lower computational cost than more rigorous algorithms, and the widely used the DIA method. Comparison of the model with buoy observations revealed a negligible difference in significant wave heights but pronounced bias in peak frequency with the DIA. The analysis of spectral shape indicated that the SRIAM method can quantitatively capture the overshoot phenomena around the spectral peak during wave growth.
The TSA method [5] resides in shearing the whole wave spectrum under the kinetic integral in two parts: a high frequency and the energy containing ones. This allows excluding a reasonable value of the quadruplets from the NL-term calculations, saving accuracy. The TSA has recently been presented as a new method to estimate nonlinear transfer rates in wind waves, and has been tested for idealized spectral data, as well as for observed field measurements [6]. The TSA has been implemented in the wind-wave model WAVEWATCH III and shown to perform well for wave spectra from field measurements, even for cases with directional energy shearing, compared to the DIA.
Regarding the DIA, the papers by van Vledder [7] and Tolman [8,9], devoted to an extension of the DIA to multi-configuration versions, are the closest to our topic. Tolman [8] presented a generalized multiple DIA (GMD) which expands upon the DIA by (i) expanding the definition of the representative quadruplet, (ii) formulating the DIA for arbitrary water depths, (iii) providing complimentary deep and shallow water scaling terms, and (iv) allowing for multiple representative quadruplets. The GMD is rigorously derived to be an extension of the DIA, and is backward compatible with it. The free parameters of the GMD have been optimized holistically, by optimizing full model behavior in the WAVEWATCH III wave model [9]. The results showed that in deep water, GMD configurations can be found which remove most of the errors of the DIA. Most of these improvements were implemented in a new version (4.18) of the WAVEWATCH code.
In these papers [7][8][9] (and further elaborations of them), they did not avoid the main shortage of the DIA-the location of all the quadruplets exactly at a singular subsurface defined by the resonant conditions (2).
The shortages of the original DIA become clear if we consider details of the original DIA.

Details of the Discrete Interaction Approximation
The quadruplet used in the original DIA in the polar coordinates (σ, θ) has the following parameters (see Figure 1): In consistency with the resonant conditions (2), the parameters of the configuration are λ= 0.25, ∆θ + = 11.5 • , and ∆θ − = 33.6 • .
The nonlinear transfer at all the mentioned k-points takes the form [3] NL where In Equation (6), the fitting dimensional constant, C, proportional to the integration k-space, has value C = 3000 which is valid for the integration grid used by Hasselmann et al., [3]. (The value of C depends not only on the grid resolution parameters but on the kind of the source term of the model, as well). The net nonlinear transfer at any fixed (σ, θ)-point of the numerical grid is found by running external vector k in Equations (6) and (7) through all the points of the frequency-angle integration grid σ i , θ j .
Note that the set of Equation (6) provides the conservative feature for the DIA [3] regardless the location of the quadruplet with respect to the singular subsurface.
The main advantage of this approximation is its evident simplicity. As the analogue of the method of "the integration-in-mean", the DIA has certain accuracy for a certain initial spectrum [3]. The mean error of the DIA (relative to the exact calculations of I NL ) is about 60%, if estimated on the ensemble of different wave-spectra shapes [12]. This is the first shortage of the original DIA.
Nevertheless, due to the several orders increase of the speed of calculation for NL-term, the DIA is widely used in practical wind-wave modeling [1]. The third-generation wind-wave numerical models, WAM [13] and WAVEWATCH (WW) [14], are the examples of successful implementation of the DIA.
The other technical shortage of the DIA routine resides in a presence of intermediate and cumbersome interpolation procedures induced by the mismatch of the spectral grid nodes and vectors k + , k − (see Equation (4)). This leads to the time-consuming about 50% CPU for the nonlinear evolution-term calculation during the numerical simulation of wave-field evolution [12].
The radical improvement of the DIA was done by Polnikov and Farina [12], who located all the interacting wave vectors of the DIA configuration exactly at the nodes of the frequency-angular grid, σ i , θ j , used in both the kinetic integral and numerical model under application. This provides a nearly two-times enhancement of the speed of numerical calculation for the NL-term in a wind-wave model (with preservation of the accuracy). For this reason, the proposed version of the DIA was called as the fast DIA (FDIA).
Below, we demonstrate all details of the FDIA elements, based on several frequency-angular numerical grids of high-resolution. The aim of this demonstration was to stimulate an active implementation of the FDIA in modern wind-wave models.

The Concept of the FDIA
In the original version of DIA [3], two of four interacting vectors (i.e., k + , k − ) are not located at the nodes of integrating grid, what leads to the necessity of the spectrum interpolation. For this reason, the speed of numerical wave forecast calculations is remarkably reduced. The main idea of the FDIA is to use quadruplets which are adjusted to the integration grid for the kinetic integral. To specify this idea, first of all, one should introduce the principal parameters of the grid. Then, the features of configurations in FDIA could be described.

The Integration Grid Properties
Integration grid for kinetic integral will be considered in the polar co-ordinates where each of interacting wave vector k i is represented by the frequency-angular point (σ i , θ i ). Usually [1,3], the integration grid is given by the formulas Thus, parameters of the grid are as follows: • the lowest frequency, σ 0 ; • the frequency exponential increment, q; • the maximum number of frequencies, I; • the angle resolution in radians, ∆θ; • and the maximum number of angles, J.
To our aims, the principal parameters are q and ∆θ, as far as they define the resolution of the grid. The numbers I and J should be rather great (several tens), but for the concept under consideration their explicit values I and J are not principal. Note only that the FDIA concept is valid for the rather fine grid (to save an accuracy) when q ≤ 1.1 and ∆θ ≤ π/10.
Everywhere below, the restriction (8) is supposed to be met. Initially, the FDIA was proposed in [12] for the resolution parameters q = 1.05 and ∆θ= π/18, that is called as the "standard" integration grid which was used in [12] for the exact calculation of the kinetic integral I NL based on the author's method described in [15].

The Choice of Configuration
In the FDIA, the so called "basic configuration", that is the closest to the original DIA, is described by the following ratios. (Pay attention that in the FDIA, a choice of the independent vectors of a quadruplet is changed, accepting k 4 as the external variable). (1) First, we fix the external vector of integration which is located at a current grid node (σ 4 , θ 4 ) represented in polar co-ordinates. (2) Then, we fix vector which is also located at numerical grid node. Here, the new parameter of the configuration, ∆θ 34 = θ 3 − θ 4 , is the angle between vectors k 4 and k 3 . These two vectors make the summary vector k a = k 4 + k 3 as a benchmark for directions for the other two vectors, as far as all the vectors of a quadruple are to be allocated in the vicinity of the resonance "figure-of-eight" in the k-space (Figure 1) [3]. In the original DIA, vectors k 1 and k 2 are simply located on this vector k a = (k a , θ a ) making the basis for the DIA configuration. This is the principle difference between quadruplets in the DIA and the FDIA (3) Finally, we choose vectors k 1 and k 2 , which are also to be allocated at the nodes of the grid, to be directed closely to the direction of vector k a , what is defined by the ratios Vector k a plays the role of the reference direction along the angle θ a = θ 4 + ∆θ a4 , where parameters k a and σ a , in terms of the independent variables, σ 4 , σ 3 , and ∆θ 34 , have the kind and k a = σ 4 4 + σ 4 3 + 2σ 2 σ 2 3 cos(∆θ 34 ) Herewith, the difference between angles θ a and θ 4 is given by the ratio whilst the correspondence of the quadruplet location near the figure-of-eight is given by the ratio [16]: Thus, after fixing vectors k 4 and k 3 , and determining ∆θ 34 , Equations (11)- (14) determine the values k a and σ a for the given σ 4 and σ 3 . After that, the expression for θ a = θ 4 + ∆θ a4 finalize possibilities to choose vectors k 1 and k 2 . Varying independent parameters for k 4 , k 3 and ∆θ 34 (below they are called as "general"), one can vary the values for the dependent parameters θ a , k a , and σ a , determining possible positions for vectors k 1 and k 2 .
The main differences between the configurations used in the FDIA and the original DIA are as follows: (a) All wave vectors k 1 , k 2 , k 3 , and k 4 should be allocated at the nodes of the integration grid; (b) vectors k 1 and k 2 may be unequal, i.e., they may have some (but small) discrepancies in both values and directions ( Figure 2); (c) the resonance conditions (2) may be rather approximately met, and the quadruplet may be unclosed ( Figure 2). (a) All wave vectors 1 k , 2 k , 3 k , and 4 k should be allocated at the nodes of the integration grid; (b) vectors 1 k and 2 k may be unequal, i.e., they may have some (but small) discrepancies in both values and directions ( Figure 2); (c) the resonance conditions (2) may be rather approximately met, and the quadruplet may be unclosed ( Figure 2).

Specification of the Configuration Parameters
To specify the FDIA configuration, it needs to define several integer values corresponding to the abovementioned requirement (a) (allocation of the vectors on the grid). According to the grid (7a), this requirement can be expressed via a set of integer digits by the following equations: Here m1, m2, n3, and na are the integer values to be found for any given integer m3. The first two are found from requirement (10c), and the latter two do from formulas and the corresponding choice for the angle parameters of the vectors 1 k and 2 k : where n1, n2 are the angular parameters of the vectors 1 k and 2 k corresponding to Equations (11c) and (16b). Sign (±) means the permutation symmetry for vectors 1 k and 2 k .
The choice of (±1) means a possible inequality of vectors 1 k and 2 k . Equations (16b) and (18) mean that for a certain configuration, given by values m1, m2, m3, and n1, n2, n3, the angle parameters of interacting vectors have the values

Specification of the Configuration Parameters
To specify the FDIA configuration, it needs to define several integer values corresponding to the abovementioned requirement (a) (allocation of the vectors on the grid). According to the grid (7a), this requirement can be expressed via a set of integer digits by the following equations: and Here m1, m2, n3, and na are the integer values to be found for any given integer m3. The first two are found from requirement (10c), and the latter two do from formulas n3 = Int(∆θ 34 /∆θ), na = Int(∆θ a4 /∆θ) (16c) via the previously determined values for ∆θ 34 and ∆θ a4 (as it is described above). In (16c), the function Int( . . . ) means the integer number nearest to the value of the argument. Requirement (b) (inequality of vectors k 1 and k 2 ) means that one can use the following choice for modulus parameters of the vectors |k 1 | and |k 2 |, defined via σ 1 and σ 2 : and the corresponding choice for the angle parameters of the vectors k 1 and k 2 : where n1, n2 are the angular parameters of the vectors k 1 and k 2 corresponding to Equations (11c) and (16b). Sign (±) means the permutation symmetry for vectors k 1 and k 2 . The choice of (±1) means a possible inequality of vectors k 1 and k 2 . Equations (16b) and (18) mean that for a certain configuration, given by values m1, m2, m3, and n1, n2, n3, the angle parameters of interacting vectors have the values where sing (±) denotes a set of two mirror configurations due to the matrix M 1,2,3,4 symmetry (see [3]). Taking into account the change of the interacting wave vectors order, the net expression for NL-term in the FDIA (for the energy-spectrum representation: NL S (k) ≡ ∂S(k)/∂t) is given in the (σ, θ)-coordinates by the formulas NL S (σ 4 , θ 4 ) = I(σ 1 , θ 1 , σ 2 , θ 2 , σ 3 , θ 3 , σ 4 , θ 4 ), Fluids 2020, 5, 176 7 of 12 NL S (σ 1 , θ 1 ) = −I(σ 1 , θ 1 , σ 2 , θ 2 , σ 3 , θ 3 , σ 4 , θ 4 ), (20c) where [12], and the dimensional fitting constant C is depending on the grid parameters. The final 2D-function for NL-term is found by the running vector (σ 4 , θ 4 ) in Equations (20) and (21) through the whole integration grid σ i , θ j , similarly to the original DIA procedure.
After some numerical simulations for the grid (7a,b) with parameters q = 1.1, ∆θ = π/12 (typical for the WAM), the fitting constant C in (20) is found to be equal to 12,000 [17]. In our case, the change of C is related to the change of the quadruplet configuration (and due to other evolution terms in the wind-wave model used). (This fitting coefficient C is tuned to the total source function of the wind-wave model proposed in [17]).
Hereby, the algorithm of the FDIA configuration calculations is fully described. The certain set of configurations will be given in the next section. It needs only to add that effectiveness of the FDIA against DIA was numerously and successfully verified in comparison with the WAM [17][18][19][20] and the WW [21].

Parameters for Several Efficient FDIA Configurations
To compare the efficiency of any DIA configurations, a certain criterion was derived in [12], based on accuracy of an approximation (as a measure of discrepancy between approximated and exact estimations of I NL ) and the time of calculation for a tested DIA configuration. Owing to this criterion, we found several the most efficient FDIA configurations presented below. In the case of unequal k 1 and k 2 , the configuration has a symmetry with respect to permutation k 1 ↔ k 2 , that means the possible permutation (m2, n2) ↔ (m1, n1), making the new quadruplet.
Pay attention that from Table 1 some other configurations are seen which could be used for the DIA. As it was shown in [22], some of them are more efficient than that given by (23a,b) for the original DIA. The relative efficiency for these configurations was checked by means of the special method constructed to this task and presented in [12].
Case 4. For the grid parameters with more fine angular resolution as q = 1.1 and ∆θ= π/18 (or ∆θ = 10 • ), Table 1 has the same kind, whilst the proper configurations are as follows.

Single Configuration
For applications which can be applied in the future, the following very high-resolution grid is preferable q = 1.05 and ∆θ= π/18 (or ∆θ = 10 • ).
Principal and auxiliary parameters for this grid are presented in Table 2. Here within, in the case of the grid (31), the most efficient FDIA single configurations, which could be used in practice, are presented in Table 3 (for a proof of relative efficiency of these configurations among all other configurations, see [22]). Table 3. Principle parameters for the set of the most efficient single configurations for grid (31). 8  4  5  3  2  2  2  S2  8  4  5  3  2  3  2  S3  9  5  5  4  3  3  3  S4 *  9  4  5  4  3  3  2  S5  10  5  6  4  3  3  3  S6  10  6  6  4  3  3  3 Notes. 1. Index of configuration includes the symbol of the single configuration type, S, and the conventional number of configuration (for notations, see [14]). Supindex "*" means that the configuration has unequal n1 and n2; 2. Configuration S6 is marked as the closest one to the original DIA configuration; 3. Parameters m3, n3, na are marked as "general" as far as m3 is independent parameter, and n3, na are directly defined by formulas (15a,b) and constant for a given m3.

Multiple Constructions of Single Configurations
Finally, we add that some multiple configurations (i.e., constructions of several single configurations [22]) which are more efficient than the simple configurations mentioned in Table 3. These constructions are presented in Tables 4 and 5 given for auxiliary configurations.

Index of Construction Composition of Two Simple Configurations
M5 S1 + S8 M6 S1 + 0.7·S8 * M7 S1 + S10 M8 S1 + 0.7·S10 Note. The coefficient in front of configuration means the weight of a proper single configuration from Tables 3 and 5. Supindex "*" means that the configuration has unequal n1 and n2. Notes. For legend, see notes for Table 3. Super-index "*" means that the configuration has unequal n1 and n2.
These are the parameters of the auxiliary simple configurations used in Table 4.

Remarks on the Efficiency
For the completeness of the text, the digitized values of the conventional efficiency parameter, E f f, for the abovementioned single configurations and multiple constructions are presented in Table 6.

Discussion
The DIA was proposed in 1985 [3], and for a long time was unchanged for the reasons of complexity of the point. Some ideas of improving the DIA were declared by van Vledder in [7], but the radical step was made by Polnikov and Farina in [12]. This was possible due to owning the routine for the exact calculation of the kinetic integral [15], that allows formulating the criterion of comparing an efficiency of different versions for DIA and its modifications. Finally, the idea of locating the interacting wave vectors at the nodes of the numerical grid was proposed and realized in [12]. Despite of mismatch of the exact resonance conditions (2), the conservative feature of the NL-term is saved in the FDIA due to ratios (20) (analogous to ratios (6)).
It is found that this modification provides not only an enhancement of the speed of calculation of the NL-term but has better accuracy as well. The calculation speed is increased due to eliminating the interpolation procedures in the original DIA, whilst the better accuracy of FDIA is due to the better choice of the quadruplet configuration [12,22].
This double positive effect in calculation of the NL-term is due the fact of rather crude efficiency of the original DIA (the mean error is about 60% [12]), and better choice of the configuration (see Sections 4.1 and 4.2 above). For the NL-term, the FDIA provides the increase of accuracy in 20%, whilst the speed of calculation is enhanced nearly twice. The tables of comparison for the accuracy and time-consuming values of FDIA and DIA are not given here to save the room of this paper. They are presented in the numerous early papers, both for the net NL-term [12,22] and for the real wind-wave models WAM and WW [18][19][20][21].
Based on these results, the FDIA was implemented in the National Institute of Oceanography in India [20]. It is still left to spread this positive result to the modern versions of the world-wide used models: WAM and WW. The present paper aimed to prompt this implementation.

Conclusions
Details of the original discrete interaction approximation (DIA) are presented, and the concept of the fast DIA (FDIA) is comprehensively described.
Numerous versions of the FDIA configurations for different numerical grids are presented, including the single and multiple configurations in a high-resolution case.
The preference of the FDIA compared to the original DIA in accuracy and time-consuming are mentioned and explained. Some estimations of increased efficiency of the FDIA are shown.
Funding: This research received no external funding.