 Home Engineering  # Numerical efficiency

An important question is the computational efficiency of the model. In the implicit time marching scheme implemented in the Lagrangian solver, the most time consuming element is the inversion of a Jacobi matrix used by the Newton iterations for solving nonlinear grid equations. The required calculation time grows fast with increased matrix dimension, which equals to the number of mesh equat ions and is proportional to the product of the dimensions of a numerical mesh. For a matrix inversion algorithm used in this work, the inversion time is approximately proportional to the square of the dimension of the matrix. For a 2D problem with a constant value of Nz much smaller than Nx, the matrix dimension is proportional to Nx. To provide a uniform discretisation error in space and time, the value of time step should be proportional to 8x which implies that the number of steps is proportional to Nx. Thus, the overall calculation time grows with a rate proportional to A'3, as illustrated in Figure 1 for the cases presented in Table 1.

High demand for computational resources for large scale problems is a well recognised disadvantage of implicit schemes, which often outweighs their advantages in numerical stability. The radical method of increasing computational efficiency is using parallel computing. Implicit solvers have a single standard time-consuming operation and are, therefore, suitable for efficient parallelisation. A parallel version of an implicit solver can be created with minimal changes to the original code by replacing a matrix inversion subroutine with a parallel analogue. This feature is particularly useful in light of recent advances in GPU-based matrix inversion algorithms, which are much faster than conventional parallelisation using multiple CPUs (e.g., Sharma et al., 2013). Figure 1: Computation time for modelling of 20 sec of wave propagation by the Lagrangian solver on a single 2.4 GHz CPU for different mesh sizes: Nz = 11 (solid); Nz = 16 (dashed); Nz =21 (dash-dotted).

Table 1: Numerical cases for convergence and validation tests of the Lagrangian model.

 Nx Nz St, sec Nx Nz St, sec 101 11 16 21 0.010 251 11 16 21 0.004 126 11 16 21 0.008 401 11 16 0.0025 201 11 16 21 0.005 501 11 0.002

# Model validation

We validate numerical results against a set of experimental data on propagation of focussed wave groups obtained in a wave flume of the Civil Engineering department at UCL. The flume has a width of 45 cm and the length of the working section between two piston wavemakers is 12.5m (see Figure 2). A paddle on the right end of the flume is used as a wave generator and the opposite paddle as an absorber. Water depth over the horizontal bed of the flume was set to h = 40 cm. We use the centre of the flume as the origin of the coordinate system with the ж-axis directed towards the wave generator positioned at x = 6.25 nr. The vertical г-axis with the origin at the mean water surface is directed upwards. The wavemaker uses a control system with force feedback operating in a frequency domain, which allows precise control and partial absorption of incident waves to reduce reflections. The input of the control system is the desired linear amplitude spectrum of the wave at the centre of the flume. The control system uses discrete spectrum and generates periodic paddle motions. For our experiments we used an overall return period of 128 sec, which is the time between repeating identical events produced by the paddle. This means that the wavemaker generates the discrete spectrum with frequencies n/128Hz, where n is an integer. We use n = 1... 512 to generate 512 frequency components in the range from 1/128 to 4 Hz. The waves generated were monitored by a series of resistance wave probes measuring surface elevation, and an ultrasonic sensor was used to record paddle motion.

We apply an iterative procedure (Buldakov et ah, 2017) to generate a Gaussian wave group with peak frequency of fp = 1 Hz focussed at the centre of the flume with the linear Figure 2: Wave flume layout and positions of wave probes. Figure 3: Linearised spectra of experimental wave groups at the focus point (x = 0). Amplitude (left) and phase (right). Solid- A = 2.5cm; dashed- A = 5cm; dash-dotted A = 7.5 cm.

focus amplitude of 2.5 cm. Then, we use the resulting input spectrum to generate higher amplitude waves by multiplying the input by factors 2 and 3 without further focussing or spectrum corrections. We, therefore, obtain additional waves with linear focus amplitudes A = 5 cm and A = 7.5 cm. The linearised spectra of generated waves at x = 0 are shown in Figure 3. As can be expected, non-linear defocussing and transformation of the spectrum can be observed for higher amplitude waves. The linearisation is done using the spectral decomposition method described in Buldakov et al. (2017), which requires generation of waves with 4 constant phase shifts Дф = 0, jt/2, jt, 3tt/2. In this section we are using results with Дф = 0, which corresponds to a peak-focussed linear wave. The resulting waves are of three distinct qualitative types. The small amplitude wave (A = 2.5 cm) has weakly non-linear features. The wave with A = 5 cm can be described as a strongly nonlinear nonbreaking wave. And the high amplitude wave (A = 7.5 cm) exhibits two intensive breaking events, while it propagates along the flume.

We reconstruct the experimental setup using a numerical wave tank based on the previously described Lagrangian numerical model. The dimensions of the NWT are the same as the dimensions of the experimental tank. The wave is generated by implementing the boundary condition (6) on the right boundary of the computational domain with Xn{z. t) specified by the paddle displacement recorded in the experiment. To account for gaps between the paddle and walls and bottom of the experimental wave tank and for frictional losses at the wavemaker the amplitude of the numerical wave generator is reduced by 18.5%. The surface displacement damping term in (9) is activated near the tank wall opposite the wavemaker to absorb the reflections. Calculations are performed for all experimental cases, including 3 amplitude values and 4 phase shifts, and repeated with different time steps and horizontal and vertical numbers of mesh points. The time step has been reduced with the increasing number of horizontal mesh points to maintain the dispersion error levels due to temporal and spatial discretisation as given by (15), approximately equal to each other. This provides uniform convergence by both parameters. The summary of parameters for the numerical cases is given in Table 1.

The convergence is tested using an L-2 norm of the difference between the experimental and calculated spectral components of the surface elevation at the linear focus point x = 0. The norms for spectral amplitudes and phases are calculated as where the sum is taken over discrete spectral components in the range 0.5 Hz < / < 1.5 Hz from the set generated by the wavemaker (n = 64... 192), where the amplitude components are large enough not to be affected by experimental errors. It should be noted that full convergence of the calculated results to the experimental measurements can not be expected. The numerical model is based on a set of assumptions that are satisfied with limited precision. In addition, the measurements include some experimental errors. Therefore, the difference between the experimental and numerical results converges to a certain small value and does not change with an additional increase in the resolution of the numerical model.

The selected results of the convergence tests are presented in Figures 4-6. As can be seen in the top row of Figure 4, the results with and without dispersion correction converge to the Figure 4: Z/2-norm of difference between amplitude (left) and phase (right) spectra of a wave group at x = 0 measured in experiment and calculated by the Lagrangian solver. Top - A = 5 cm; without (solid) and with (dashed) dispersion correction. Bottom - A = 7.5 cm; breaking control without (solid) and with (dashed) crest correction. Nz = 11. Figure 5: Convergence of time history of surface elevation for a wave group at ж = 0. Experiment (thick dashed) and Lagrangian solver: Nx = 101 (dashed); Nx = 201 (dash- dotted); Nx = 401 (solid); Nz = 11. Top - A = 5cm without (left) and with (right) dispersion correction. Bottom A = 7.5cm, breaking control without (left) and with (right) crest correction. Figure 6: Convergence of the normalised total energy for the Lagrangian solver with dispersion correction. A = 5 cm. Left - increasing horizontal mesh resolution: Nx = 101 (dashed); Nx =201 (dash-dotted); Nx =401 (solid); Nz = 11. Right - increasing vertical mesh resolution: Nz = 11 (dashed); Nz = 16 (dash-dotted); Nz = 21 (solid); Nx = 251.

same solution. However, the introduction of a dispersion correction considerably increases the speed of convergence for both amplitudes and phases. The experimental results can be reproduced with sufficient accuracy for a relatively small number of horizontal mesh points and a relatively large time step. For a breaking wave (Figure 4, bottom row), the shape correction of the crest introduces a new physical process. For this reason, the numerical results with and without crest correction converge to different solutions, and the numerical result with the correction shows a much better comparison with the experiment. A general impression of convergence and accuracy of the different versions of the numerical model can be obtained from the graphs of the time history of surface elevation presented in Figure 5.

Figure 6 shows the behaviour of the total wave energy in the wave tank calculated for different temporal and horizontal resolutions and different resolutions of the vertical mesh. Wave absorption is disabled for energy tests. The energy is normalised by the energy of one wave length Л of a linear regular wave with the frequency and the amplitude equal to the peak spectral frequency the linear focus amplitude of the wave group. The kinetic energy is calculated by numerical integration over the entire Lagrangian fluid domain using a bilinear interpolation of the velocities of the fluid particles inside mesh cells. This provides the second order approximation of the integral with respect to the mesh resolution. According to (13), the error of the velocity profiles within the fluid domain, and thus the kinetic energy, is determined by the horizontal and vertical resolution of the mesh. The potential energy is calculated as an integral of the potential energy density at the free surface and is not directly affected by the discretisation of the vertical mesh. As can be seen in Figure 6, the total energy in the tank increases up to t яз 4 sec due to the energy generated by the wave paddle. For low horizontal mesh resolution, the dissipative term in the numerical dispersion relation leads to the reduction of the energy of the propagating wave. For higher resolutions, energy conservation is satisfied with high accuracy. The right graph of Figure 6 shows the rapid convergence of energy with increased vertical resolution. This includes both the convergence of the numerical integral used to calculate kinetic energy and the convergence of the numerical solution for wave kinematics, as indicated by expansion (13). Overall, the convergence tests show that the numerical results converge towards a solution, which approximates the experiment with a good accuracy. The implementation of the dispersion correction term greatly increases the convergence rate and provides accurate results with a smaller number of mesh points and a larger time step.

 Related topics