A GPU Accelerated Method for 3-D Nonlinear Kelvin Ship Wave Patterns Simulation

1. Introduction

The focus of this research is on innovation of highly parallel algorithms to simulate the contours of a three-dimensional free surface that appears to be stationary at the stern of a moving vessel, known as “Kelvin ship waves” [1]. Research on the kelvin wave shape has been continuously put to practical use in hull design, ship detection and environmentally friendly shipping policies [2].

Froude [3], a famous naval architect, first comprehensively described the morphology and main characteristics of ship waves. Under the assumption of infinite water depth, Kelvin [4] replaced a moving ship with a pressure disturbance point moving in a constant velocity straight line on the water surface, proposed the famous Kelvin angle of $19.{47}^{\circ }$. In recent years, with the further study of ship wave characteristics, Rabaud [5] noted that the wake angle will be less than the well-known Kelvin angle if the vessel speed is sufficiently large. Subsequently, various effect factors for the Kelvin wake form were discussed in plenty of papers, e.g., Froude number [6], non-axisymmetric and interference effects, shear current,surface tension, the bottom topography, submergence depth, finite water depth and viscosity [7], etc. Accordingly, the research method of ship waves has gradually shifted from the previous analytical algorithms to numerical simulation.

The overwhelming majority of analytical algorithms of ship wave patterns concerns linear theories. Havelock [8] provided a linear solution for the problem of flow under a pressure distribution. Such ideal perturbations can also be replaced by a single submerged point source singularity [8] and submerged bodies [2]. Moreover, thin ship theory was alse used in the study of the ship wave pattern [9]. As the development of computer technology, numerical simulation methods are becoming increasing popular, the research focus shifted from linear problems to nonlinear problems. Nowadays, there are three numerical methods widely used to solve surface wave problems: boundary integral method, finite-difference method and finite-element method. In particular, Forbes [10] apply boundary integral method to build a series of integro-differential equation, the full nonlinear free surface flow problem was solved with moderate efficiencies. In more recent times, according to this method, many papers solve fully 3-D nonlinear ship waves with meshes between $60\times 20$ and $181\times 61$[11,12]. And Pethiyagoda [6] noted that the points used along the $x-$direction should be more than 100 to make a sufficient standard about grid-independence.

With increasing mesh size, however, the computation time increases exponentially using only Central Processing Unit (CPU) computation power. As the rapid improvement of the electronics industry, the Graphics Processing Unit (GPU) has become another method of acceleration for optimizing the execution of large numbers of threads. Currently, the powerful GPU parallel computing ability has been used to improve the studies on ocean engineering. Hori [13] simulated 2-D dam-break flow by developing a GPU-based MPS code and achieved 7 times speedup. As for 3-D nonlinear free surface problem, Pethiyagoda [6] combined the GPU acceleration technique with the boundary integral method and LU [14] developed a GPU-accelerated high-order spectral solver. Xie [15] developed the MPSGPU-SJTU solver with GPU acceleration technique for the liquid sloshing simulation.

This paper presents a parallel solution framework based on GPU for nonlinear ship wave problem, in which almost all operations are performed in GPU device. Since the nonlinear boundary integral equation on each node is independent of the synchronous equations on other nodes, plenty of threads on GPU can be used to complete the integration operation for each node simultaneously. In addition, the parallel computing method can be used for the calculation of the large-scale linear sparse system, the complex inversion process is quickly finished by using Compute Unified Architecture (CUDA) language. According to this framework, a highly-paralleled GPU solver is proposed to simulate 3-D nonlinear Kelvin ship waves. The computation speed for the 3-D nonlinear ship waves simulation can be significantly increased, it is convenient to study the larger scale problems. On the other hand, the size of Random-Access Memory limits grid growth, the application of the banded preconditioner method can greatly save running memory to break through this limitation. The banded preconditioner method helps to achieve the standard for the grid-independence.

The rest of the paper is as follows. A brief introduction of the problem formulation is given in Section 2. In Section 3, the banded preconditioner JFNK algorithm is described. In Section 4, the theory and implementation of the GPU acceleration technique are presented. The accuracy, efficiency and capability of the GPU solver are verified in Section 5, and a summary in Section 6 concludes the paper.

2. Numerical Model

This paper supposes that a flow is moving at a uniform speed U along the positive x-axis direction. Considering the inviscid incompressible fluid of infinite depth without rotational flow, ignoring the influence of surface tension, the potential flow theory is applied. Therefore, a source singularity of strength m is introduced at a distance L below the surface, as illustrated in Figure 1. The transient waves can be generated with the disturbance of source. Free surface wave height and flow field velocity potential can be expressed as $z=\zeta (x,y)$ and $\Phi (x,y,z)$.

Dimensionless analysis is performed with fluid velocity U and distance L. The velocity potential $\Phi (x,y,z)$ satisfies Laplace’s equation, the free surface kinematic and dynamic boundary condition, the radiation condition and the limiting behavior of source singularity. With $\varphi \left(x,y\right)=\Phi \left(x,y,\zeta \left(x,y\right)\right)$, the boundary integral equation is written:

$$\begin{array}{ccc}\hfill & 2\pi (\varphi \left(q\right)-x)={\int }_0^{\infty }{\int }_{-\infty }^{\infty }[\varphi \left(p\right)-\varphi \left(q\right)+x-\rho ]K_{1}d\sigma d\rho \hfill & \hfill \\ & +{\int }_0^{\infty }{\int }_{-\infty }^{\infty }{\zeta }_{\rho }\left(P\right)K_{2}d\sigma d\rho -\frac{ϵ}{{[y^2+x^2+{\left(\zeta \right(q)+1)}^2]}^{\frac{1}{2}}}\hfill \end{array}$$

(1)

where the $K_1$ and $K_2$ are kernel functions [12].

Figure 1. Flow field diagram.

Figure 1. Flow field diagram.

Moreover, the free surface conditions can be simplified by the symbol $\varphi \left(x,y\right)$. Then the kinematic and dynamic boundary conditions of the free surface are combined to be

$$\begin{array}{cc}\hfill & \frac{{(1+\zeta }_x^2){\varphi }_y^2+\left(1+{\zeta }_y^2\right){\varphi }_x^2-{2{\zeta }_x{\zeta }_y\varphi }_x{\varphi }_y}{{2(1+\zeta }_x^2+{\zeta }_y^2)}+\frac{\zeta }{F^2}=\frac{1}{2}\hfill \end{array}$$

(2)

To solve the above nonlinear problem numerically, the $N\times M$ mesh is established on the free surface (N and M represent the number of longitude and latitude lines of the mesh, respectively). The x-coordinates and y-coordinates of nodes are $x_1,x_2,\dots ,x_N$ and $y_1,y_2,\dots ,y_N$ with regular intervals in the coordinate system, thus the vector $\mathit{u}$ of $2\left(N+1\right)M$ unknowns is

$$\begin{array}{ccc}\hfill \mathit{u}=& [{\varphi }_{1,1},{\left({\varphi }_x\right)}_{1,1},\dots ,{\left({\varphi }_x\right)}_{N,1},\dots ,{\varphi }_{1,M},{\left({\varphi }_x\right)}_{1,M},\dots ,{\left({\varphi }_x\right)}_{N,M}\hfill & \hfill \\ & {\zeta }_{1,1},{\left({\zeta }_x\right)}_{1,1},\dots ,{\left({\zeta }_x\right)}_{N,1},\dots ,{\zeta }_{1,M},{\left({\zeta }_x\right)}_{1,M},\dots ,{\left({\zeta }_x\right)}_{N,M}{]}^T.\hfill \end{array}$$

(3)

With more $4M$ equations provided by applying the radiation condition as follows:

$$x_1({\left({\varphi }_x\right)}_{1,l}-1)+\gamma ({\varphi }_{1,l}-x_1)=0$$

(4)

$$x_1({\left({\varphi }_{xx}\right)}_{1,l}-1)+\gamma ({\left({\varphi }_x\right)}_{1,l}-1)=0$$

(5)

$$x_1{\left({\zeta }_x\right)}_{1,l}+\gamma {\zeta }_{1,l}=0$$

(6)

$$x_1{\left({\zeta }_{xx}\right)}_{1,l}+\gamma {\left({\zeta }_x\right)}_{1,l}=0$$

(7)

where $\gamma$ is the decay coefficient.

Furthermore, more details about the governing equations, the boundary integral method and numerical discretization are provided by Sun[12].

3. Banded Proconditioner JFNK Algorithm

3.1. Jacobian-free Newton-Krylov method

JFNK method combines inexact Newton iteration method with Krylov subspace method. Its core content is the Generalized Minimum Residual (GMRES) algorithm, according to the matrix free idea, uses the finite difference form to approximate the product of coefficient matrix and vector, avoiding the Jacobian matrix calculation and storage alone. JFNK method mainly has two processes, namely external and internal iterations. The external iteration is the inexact Newton iteration method, and the damping parameter ${\lambda }_k$ is used to ensure that the nonlinear residual decreases significantly in each iteration for $t=0,1,2,\dots$, as follows:

$$\begin{array}{cc}\hfill & {\mathit{u}}_{\mathit{t}+\mathbf{1}}={\lambda }_k\mathit{\delta }{\mathit{u}}_{\mathit{t}}+{\mathit{u}}_{\mathit{t}},{\lambda }_k\in (0,1]\hfill \end{array}$$

(8)

Its internal iteration is GMRES algorithm [16], which efficiently solves the correction in inexact Newton iteration, that is, computes large-scale linear equations as follows:

$$\begin{array}{cc}\hfill & \mathrm{J}\left({\mathit{u}}_{\mathit{t}}\right)\mathit{\delta }{\mathit{u}}_{\mathit{t}}=-\mathit{F}\left({\mathit{u}}_{\mathit{t}}\right)\hfill \end{array}$$

(9)

where $\mathrm{J}\left({\mathit{u}}_{\mathit{t}}\right)=\mathbf{\partial }\mathit{F}\left({\mathit{u}}_{\mathit{t}}\right)/\mathbf{\partial }{\mathit{u}}_{\mathit{t}}$ is the Jacobian matrix [17].

Firstly, the approximate solution of $\mathit{\delta }{u}_t$ is found by projecting obliquely onto the Krylov subspace

$$\begin{array}{cc}& K_m({\mathrm{J}}_{\mathit{t}}{\mathrm{P}}^{-\mathbf{1}},{\mathit{F}}_{\mathit{t}})=\mathrm{span}\{{\mathit{F}}_{\mathit{t}},{\mathrm{J}}_{\mathit{t}}{\mathrm{P}}^{-\mathbf{1}}{\mathit{F}}_{\mathit{t}},\dots ,{\left({\mathrm{J}}_{\mathit{t}}{\mathrm{P}}^{-\mathbf{1}}\right)}^{\mathbf{m}-\mathbf{1}}{\mathit{F}}_{\mathit{t}}\}\hfill \end{array}$$

(10)

where m is the value of the subspace dimension.

An initial linear residual ${\mathit{r}}_{\mathbf{0}}$ is defined, given an initial guess ${\mathit{u}}_0$, for the Newton correction,

$$\begin{array}{cc}& {\mathit{r}}_{\mathbf{0}}=-\mathit{F}\left({\mathit{u}}_{\mathbf{0}}\right)-{\mathrm{J}}_{\mathbf{0}}{\mathrm{P}}^{-\mathbf{1}}\mathit{\delta }{\mathit{u}}_{\mathbf{0}}\hfill \end{array}$$

(11)

Subsequently, $∥{\mathit{r}}_{\mathit{t}}∥$ is minimized to a suitable value by the GMRES iteration. Wherein Jacobian-vector products are approximated with finite difference:

$$\begin{array}{cc}\hfill & {\mathrm{J}}_{\mathit{t}}{\mathrm{P}}^{-\mathbf{1}}\mathit{v}\approx \frac{\mathit{F}({\mathit{u}}_{\mathit{t}}+\mathit{h}{\mathrm{P}}^{-1}\mathit{v})-\mathit{F}\left({\mathit{u}}_{\mathit{t}}\right)}{h}\hfill \end{array}$$

(12)

where v represents an arbitrary vector used in building the Krylov subspace [18], and the h is a small perturbation

$$\begin{array}{cc}& h=\frac{\sqrt{(1+∥{\mathit{u}}_{\mathit{t}}∥)h_{mach}}}{∥\mathit{v}∥}\hfill \end{array}$$

(13)

Finally, the initial guess ${\mathit{u}}_{\mathbf{0}}$ can be defined as below:

$$\begin{array}{c}\hfill {\zeta }_{1,l}=0,{\left({\zeta }_x\right)}_{k,l}=0,{\varphi }_{1,l}=x_0,{\left({\varphi }_x\right)}_{k,l}=1.\end{array}$$

(14)

The nonlinear equations are solved according to the calculation flow of baned preconditioner JFNK method, as shown in Figure 2. Note that $\mathit{v}$ in the figure is a unit orthogonal vector in the orthonormal basis of Krylov subspace.

Figure 2. Calculation flow chart of the banded preconditioner JFNK method.

Figure 2. Calculation flow chart of the banded preconditioner JFNK method.

3.2. Banded preconditioner method

Iterative methods, e.g. GMRES method etc., are currently most popular choices for solving large sparse linear systems of equations. However, this process of prcconditioning is essential to most successful application of iterative methods, since the convergence of a matrix iteration depends on the properties of the matrix, e.g. the eigenvalue, etc., [19]. Generally, the methods for choosing the appropriate preconditioner are different for the specific problems. In this section, a banded preconditioner method for solving the nonlinear ship wave problem is proposed.

3.2.1. Building preconditioner matrix

For a good preconditioner $\mathit{P}$, it should be cheap to form and to factorizeis. Meanwhile, the preconditioned Jacobian ${\mathrm{J}}_{\mathit{t}}{\mathit{P}}^{-\mathbf{1}}$ should be easier to solve, which means the eigenvalues are more concentrated. In general, it is feasible to consider a matrix constructed from the same problem under simplified physics. This paper applies the numerical scheme to the linearized governing equations which apply formally in the limit $\mathit{ϵ}\to \mathbf{0}.$

The equations of the linear free surface boundary condition are described [20]:

$$\begin{array}{ccccc}& {\mathit{\zeta }}_{\mathit{x}}={\mathit{\varphi }}_{\mathit{z}}\hfill & & \mathrm{on}\mathit{z}=\mathit{\zeta }(\mathit{x},\mathit{y})\hfill & \end{array}$$

(15)

$$\begin{array}{cccc}& {\mathit{\varphi }}_{\mathit{x}}-\mathbf{1}+\frac{\mathit{\zeta }}{{\mathit{F}}^{\mathbf{2}}}=\mathbf{0}\hfill & & \mathrm{on}\mathbf{z}=\mathit{\zeta }(\mathit{x},\mathit{y})\hfill \end{array}$$

(16)

According to the linear free surface boundary condition, the boundary integral equation is described:

$$\begin{array}{cc}& \mathbf{2}\mathit{\pi }(\mathit{\varphi }\left(\mathit{q}\right)-\mathit{x})=-\frac{\mathit{ϵ}}{{({\mathit{x}}^{\mathbf{2}}+{\mathit{y}}^{\mathbf{2}}+\mathbf{1})}^{\frac{\mathbf{1}}{\mathbf{2}}}}+{\int }_{\mathbf{0}}^{\infty }{\int }_{-\infty }^{\infty }{\mathit{\varphi }}_{\mathit{\rho }}\left(\mathit{p}\right){\mathit{K}}_{\mathbf{3}}(\mathit{\rho },\mathit{\sigma },\mathit{x},\mathit{y})\mathit{d}\mathit{\sigma }\mathit{d}\mathit{\rho }\hfill \end{array}$$

(17)

After numerical discretization, the linear system can be described, as follows:

$$\begin{array}{cc}& \mathit{F}{\mathbf{1}}_{\mathbf{k},\mathbf{l}}={\mathit{\varphi }}_{\mathit{k},\mathit{l}}\left(\mathit{q}\right)+\frac{{\mathit{\zeta }}_{\mathit{k},\mathit{l}}\left(\mathit{q}\right)}{{\mathit{F}}^{\mathbf{2}}}-\mathbf{1}\hfill \end{array}$$

(18)

$$\begin{array}{c}\hfill \mathit{F}{\mathbf{2}}_{\mathit{k},\mathit{l}}=\mathbf{2}\mathit{\pi }({\mathit{\varphi }}_{\mathit{k},\mathit{l}}\left(\mathit{q}\right)-{\mathit{x}}_{\mathit{k}})+\frac{\mathit{ϵ}}{{[{{\mathit{x}}^{\mathbf{2}}}_{\mathit{k}}\left(\mathit{q}\right)-{{\mathit{y}}^{\mathbf{2}}}_{\mathbf{l}}\left(\mathit{q}\right)+\mathbf{1}]}^{\frac{\mathbf{1}}{\mathbf{2}}}}-\sum _{\mathit{i}=\mathit{1}}^{\mathit{N}}\sum _{\mathit{j}=\mathit{1}}^{\mathit{M}}\mathit{w}(\mathit{i},\mathit{j})[{\left({\mathit{\zeta }}_{\mathit{\rho }}\right)}_{\mathit{i},\mathit{j}}-{\left({\mathit{\zeta }}_{\mathit{x}}\right)}_{\mathit{i},\mathit{j}}]{\mathit{K}}_{\mathit{i},\mathit{j},\mathit{k},\mathit{l}}^{\left(\mathbf{3}\right)}-{\left({\mathit{\zeta }}_{\mathit{x}}\right)}_{\mathit{i},\mathit{j}}\mathit{I}\end{array}$$

(19)

where $\mathit{w}(\mathit{i},\mathit{j})$ is the weighting function for numerical integration, for $\mathit{k}=\mathbf{1}\dots (\mathit{N}-\mathbf{1}),\mathit{l}=\mathbf{1}\dots \mathit{M}.$ Then the linear Jacobian can be calculated directly, by differentiating the linear system with respect to ${\mathit{\varphi }}_{\mathbf{1},\mathit{m}},{\left({\mathit{\varphi }}_{\mathit{x}}\right)}_{\mathit{n},\mathit{m}},{\mathit{\zeta }}_{\mathbf{1},\mathit{m}}\mathrm{and}{\left({\mathit{\zeta }}_{\mathit{x}}\right)}_{\mathit{n},\mathit{m}}$. Therefore, the preconditioner matrix $\mathit{P}$ can be formed cheaply, and the eigenvalues of ${\mathbf{J}}_{\mathit{t}}{\mathbf{P}}^{-\mathbf{1}}$ obviously cluster as shown in Figure 3.

Figure 3. The distribution of eigenvalues of ${\mathbf{J}}_{\mathit{t}}$ and ${\mathbf{J}}_{\mathit{t}}{\mathbf{P}}^{-\mathbf{1}}$ on a $\mathbf{31}\times \mathbf{11}$ mesh.

Figure 3. The distribution of eigenvalues of ${\mathbf{J}}_{\mathit{t}}$ and ${\mathbf{J}}_{\mathit{t}}{\mathbf{P}}^{-\mathbf{1}}$ on a $\mathbf{31}\times \mathbf{11}$ mesh.

3.2.2. Preconditioner factorisation and storage

The JFNK method requires the result of the product of the inverse preconditioner matrix and vector, ${\mathbf{P}}^{-\mathbf{1}}\mathit{v}$. In general, the operation of inverting a matrix should be converted to solving a system of linear equations, $\mathbf{P}\mathit{r}=\mathit{v}$. Find the solution $\mathit{r}$, the result of ${\mathbf{P}}^{-\mathbf{1}}\mathit{v}$ will be got. In order to calculate this linear system rapidly, the following block matrix method is used to process the preprocessing matrix, where $\mathit{I},\mathit{A},\mathit{B},\mathit{C}$ and $\mathit{D}$ are the unit matrix and the four submatrices.

Accordingly, the vector $\mathit{v}$ can be divided into upper and lower parts ${\left[{\mathit{v}}_{\mathbf{1}}{\mathit{v}}_{\mathbf{2}}\right]}^{\mathbf{T}}$, and then the solution $\mathit{r}$ can be got after three cheap steps, as follows:

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}{\mathit{o}}_{\mathbf{1}}\\ {\mathit{o}}_{\mathbf{2}}\end{array}\right]=\left[\begin{array}{c}{\mathit{v}}_{\mathbf{1}}\\ {\mathit{v}}_{\mathbf{2}}-\mathit{C}{\mathit{A}}^{-\mathbf{1}}{\mathit{v}}_{\mathbf{1}}\end{array}\right]\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}{\mathit{s}}_{\mathbf{1}}\\ {\mathit{s}}_{\mathbf{2}}\end{array}\right]=\left[\begin{array}{c}{\mathit{A}}^{-\mathbf{1}}{\mathit{o}}_{\mathbf{1}}\\ {(\mathit{D}-\mathit{C}{\mathit{A}}^{-\mathbf{1}}\mathit{B})}^{-\mathbf{1}}{\mathit{o}}_{\mathbf{2}}\end{array}\right]\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}{\mathit{r}}_{\mathbf{1}}\\ {\mathit{r}}_{\mathbf{2}}\end{array}\right]=\left[\begin{array}{c}{\mathit{s}}_{\mathbf{1}}-{\mathit{A}}^{-\mathbf{1}}\mathit{B}{\mathit{s}}_{\mathbf{2}}\\ {\mathit{s}}_{\mathbf{2}}\end{array}\right]\hfill \end{array}$$

(22)

The calculation of ${\mathbf{P}}^{-\mathbf{1}}\mathit{v}$ in the Eq. (12) can be facilitated according to the progressive order from Eqs. (20) - (22).

3.2.3. The banded preconditioner

After the factorisation operation, the calculation and storage of preconditioner matrix $\mathit{P}$ are optimized. However, the size of submatrix $\mathbf{D}$ is $(\mathit{N}+\mathbf{1})\mathit{M}\times (\mathit{N}+\mathbf{1})\mathit{M}$, it will increase dramatically as the size of mesh increases. Consequently, there will be two problems when the preconditioner matrix size is large. One is a memory problem, the running memory of this computer cannot accommodate this preconditioner matrix; the other is an efficiency problem, inverting the preconditioner matrix will take much time.

By observing the preconditioner matrix, it can be found that the values decay with distance from the main block diagonal. This observation suggests using a banded approximation to the matrix for our preconditioner, as shown in Figure 4. Moreover, batch construction avoids the problem of insufficient running memory due to the large size of the submatrix $\mathit{D}$. The compressed sparse row (CSR) data format is used to save this matrix. Hence, a lot of memory can be saved.

Figure 4. Construction of the banded preconditioner.

Figure 4. Construction of the banded preconditioner.

The feasibility of the banded preconditioner matrix method is verified, as shown in Figure 5. The tightness of clustering can be further improved by increasing the bandwidth. When the $\mathit{band}=\mathbf{21}$, the eigenvalues of ${\mathbf{J}}_{\mathit{t}}{\mathbf{P}}^{-\mathbf{1}}$ have been clustered, satisfying the requirement of the GMRES method.

For certain bandwidth values, the computing speed of GMRES will not be significantly improved by increasing the bandwidth further. However, the time required for inverse operation will increase in these cases as the banded preconditioner matrix size increases. The bandwidth regulates the runtime of inverting banded preconditioner matrix and the number of the inner iterations of the GMRES method. The runtime of inverting banded preconditioner matrix increases with the bandwidth, while the inner iterations decrease with the bandwidth. Therefore, the total runtime will decrease first and then increase with the bandwidth, as shown in Figure 6. The case is $\mathit{ϵ}=\mathbf{0}.\mathbf{4}$, $\mathbf{121}\times \mathbf{41}$ mesh and $\mathit{F}=\mathbf{0}.\mathbf{7}$, when ${\mathit{b}}^{\mathbf{\prime }}$ (bandwidth $\mathit{band}={\mathit{b}}^{\mathbf{\prime }}\times (\mathit{N}+\mathbf{1})$) is less than 14, an ill-conditioned coefficient matrix is formed, the accuracy of the solution is low. The runtime decreases with ${\mathit{b}}^{\mathbf{\prime }}$ ranges from 14 to 16, then the runtime increases monotonically with ${\mathit{b}}^{\prime }$ ranges from 16 to 20. For the case of $\mathbf{121}\times \mathbf{41}$ mesh, the shortest running time is 5.6s with the optimal bandwidth $\mathit{band}=\mathbf{16}\times (\mathit{N}+\mathbf{1})$. Therefore, provided that the appropriate bandwidth is selected, not only can save memory, but also can improve the computational efficiency.

Figure 5. The distribution of eigenvalues of ${\mathbf{J}}_{\mathit{t}}$ and ${\mathbf{J}}_{\mathit{t}}{\mathbf{P}}^{-\mathbf{1}}$ on a $\mathbf{31}\times \mathbf{11}$ mesh for: $\mathit{band}=\mathbf{1},\mathit{band}=\mathbf{11},\mathit{band}=\mathbf{21}$.

Figure 5. The distribution of eigenvalues of ${\mathbf{J}}_{\mathit{t}}$ and ${\mathbf{J}}_{\mathit{t}}{\mathbf{P}}^{-\mathbf{1}}$ on a $\mathbf{31}\times \mathbf{11}$ mesh for: $\mathit{band}=\mathbf{1},\mathit{band}=\mathbf{11},\mathit{band}=\mathbf{21}$.

Figure 6. The plot of runtime against the bandwidth, computed on a $\mathbf{121}\times \mathbf{41}$ mesh, $\mathit{band}={\mathit{b}}^{\prime }\times (\mathit{N}+\mathbf{1})$.

Figure 6. The plot of runtime against the bandwidth, computed on a $\mathbf{121}\times \mathbf{41}$ mesh, $\mathit{band}={\mathit{b}}^{\prime }\times (\mathit{N}+\mathbf{1})$.

4. GPU Parallel Computing Framework

Although the banded preconditioner JFNK algorithm can improve the computational efficiency of nonlinear ship wave problem, the running time of the program will increase significantly with the increase of the mesh size, which is very unfavorable to the further study of nonlinear ship wave. The reason is that the CPU is not good at handling such large-scale nonlinear equations. Compared with CPU, GPU possesses more arithmetic logic units in the same chip area [21]. The computational efficiency of nonlinear ship wave can be greatly improved by utilizing the GPU acceleration technique.

4.1. Parallel computing framework design

Compute Unified Device Architecture (CUDA) language is used to develop the numerical scheme for computing ship wave patterns. CUDA is a parallel computer platform and programming model developed by NVIDIA, powered by the GPU [21]. CUDA toolkit includes abundant GPU accelerated libraries, tools and runtime library, which can be compiled in C language, C++ language and Fortran language. In addition, CUDA source program can be executed on multiple GPUs. By applying the hybrid programming model, the parallel computing process consists of kernel function on device and serial code on host CPU. Figure 7 shows the CUDA execution mode and thread hierarchy.

Figure 7. Illustration of CUDA execution mode and thread organization hierarchy [21].

Figure 7. Illustration of CUDA execution mode and thread organization hierarchy [21].

As for the solver of ship wave pattern, as described above, there are four main parts: building preconditioner matrix, creating nonlinear system, inverting preconditioner matrix and solving linear equations by GMRES algorithm. Simulation results of the CPU solver proposed by Sun[12] show that the parts of creating nonlinear system and inverting preconditioner matrix take up most of the time, as shown in Figure 8. This figure shows the computation time distribution of the CPU solver on a $\mathbf{151}\times \mathbf{51}$ mesh case. The total runtime is 185.8 seconds, in which the runtime of inverting preconditioner matrix and creating nonlinear system is 95.8 and 80.4 seconds respectively. Each of them takes up nearly half the total runtime. Therefore, calculations on these two parts parallelly are vital for improving computational efficiency. And the part of building preconditioner matrix and solving linear equations will also be executed in GPU to further shorten the program running time.

Figure 8. The computation time distribution of ship wave solver. The alphabet I represents the part of inverting preconditioner matrix, the alphabet C represents the part of creating nonlinear system, the alphabet B represents the part of building preconditioner matrix, the alphabet S represents the part of solving linear equations by GMRES algorithm and the alphabet O represents the part of other code in the solver.

Figure 8. The computation time distribution of ship wave solver. The alphabet I represents the part of inverting preconditioner matrix, the alphabet C represents the part of creating nonlinear system, the alphabet B represents the part of building preconditioner matrix, the alphabet S represents the part of solving linear equations by GMRES algorithm and the alphabet O represents the part of other code in the solver.

Based on the above analysis, the GPU solver of Kelvin ship waves adopts a hybrid programming model. The entire parallel computing procedure is shown as follows :

Step 1: Input calculation parameters including the initial guess, the data is transferred from CPU to GPU;

Step 2: According to the calculation parameters, the nonlinear equations are created in the GPU device;

Step 3: The banded preconditioner method is applied to build the banded preconditioner matrix in GPU;

Step 4: QR decomposition algorithm is used to invert the preconditioner matrix, and saving the decomposition results outside the loop body to avoid repeated QR decomposition of preprocessing;

Step 5: The result of ${\mathbf{P}}^{-\mathbf{1}}\mathit{v}$ is calculated directly using the QR decomposition results, by combining the result of ${\mathbf{P}}^{-\mathbf{1}}\mathit{v}$ with the approximate solution $\mathit{u}$ of the nonlinear equations, the finite difference approximation is carried out to obtain the linear equations;

Step 6: The GMRES algorithm is used to calculate the linear equations, obtain the correction values and update the approximate solutions $\mathit{u}$;

Step 7: Check of the approximate solutions of the nonlinear equations: If the accuracy requirement is not met, back to Step 5; if the accuracy requirement is met, the result is transferred from the GPU to the CPU.

The corresponding calculation flow chart is shown in Figure 9 which shows the calculation procedure more clearly.

Figure 9. The computational flow chart of GPU implementation.

Figure 9. The computational flow chart of GPU implementation.

4.2. GPU solver implementation

4.2.1. Creating nonlinear system

The programming for creating nonlinear system on GPU by using CUDA language is briefly shown in Table 1. On the whole, the dimension of the GPU grid is the equivalent of the size of mesh, which means that one block can complete the relevant equations of one node in mesh. One block has 1024 threads, these plenty of threads can calculate Eq. (1), Eq. (2) and Eqs. (4)-(7) simultaneously.

Table 1. The programming on GPU.

Device Part
Device function 1
1	$\_\_$device$\_\_$double Integral(double a,double b,double c,double d,double e)
2	{
3	double val= b/sqrt(c)*log(2*c*a + d*b+ 2*sqrt(c*(c*a*a+ d*b*b +e*b*b)));
4	return val;
5	}
Device function 2
1	$\_\_$device$\_\_$double Integrall(double a,double b,double c,double d,double e)
2	{
3	double val = a/sqrt(e)*log(2*d*b + d*s+ 2*sqrt(e*(c*a*a+ d*a*b +e*b*b)));
4	return val;
5	}
Kernel function
1	$\_\_$global$\_\_$void nonlinear(const long int N, const long int M, … )
2	{
3	$\_\_$shared$\_\_$double A, B, C, … ;
4	$\_\_$shared$\_\_$long int K, l, blockpos, … ;
5	long int i, j, threadPos, … ;
6	threadPos = threadIdx.x; k = blockIdx.x; l = blockIdx.y ;
7	// Calculate necessary values
8	while(threadPos<(M*N)) { … }
9	// Calculation of the 16 parts to the closed integral
10	if(threadIdx.x==blockDim.x-1) { … }
11	…
12	if(threadIdx.x==blockDim.x-16) { … }
13	// Sum up all thread contributions
14	for(i=blocDim/2;i>0;i=i/2){ … }
15	//Split the free surface condition and radiation conditions between 5 blocks
16	if(k==0&&l==0) { … } … if(k==0&&l==5) { … }
17	}
Host Part
1	int main
2	{
3	// Initialize data on CPU
4	InitialData(double* cpuData, double* gpuData);
5	long int i, j, threadPos, … ;
6	// Copy data from CPU to GPU
7	cudaMemcpy(cpuData, gpuData, Datasize, cudaMemcpyHostToDevice);
8	// Set dimensions of the grid and thread
9	dim3 block(1024, 1); dim3 grid(M, N - 1);
10	// Start Kernel function
11	Nonlinear≪ grid, block≪(gpuData);
12	// Copy data from GPU to CPU
13	cudaMemcpy(gpuData, cpuData, Datasize, cudaMemcpyDeviceToHost);
14	return 0;
15	}

Table 1. The programming on GPU.

Device Part
Device function 1
1	$\_\_$device$\_\_$double Integral(double a,double b,double c,double d,double e)
2	{
3	double val= b/sqrt(c)*log(2*c*a + d*b+ 2*sqrt(c*(c*a*a+ d*b*b +e*b*b)));
4	return val;
5	}
Device function 2
1	$\_\_$device$\_\_$double Integrall(double a,double b,double c,double d,double e)
2	{
3	double val = a/sqrt(e)*log(2*d*b + d*s+ 2*sqrt(e*(c*a*a+ d*a*b +e*b*b)));
4	return val;
5	}
Kernel function
1	$\_\_$global$\_\_$void nonlinear(const long int N, const long int M, … )
2	{
3	$\_\_$shared$\_\_$double A, B, C, … ;
4	$\_\_$shared$\_\_$long int K, l, blockpos, … ;
5	long int i, j, threadPos, … ;
6	threadPos = threadIdx.x; k = blockIdx.x; l = blockIdx.y ;
7	// Calculate necessary values
8	while(threadPos<(M*N)) { … }
9	// Calculation of the 16 parts to the closed integral
10	if(threadIdx.x==blockDim.x-1) { … }
11	…
12	if(threadIdx.x==blockDim.x-16) { … }
13	// Sum up all thread contributions
14	for(i=blocDim/2;i>0;i=i/2){ … }
15	//Split the free surface condition and radiation conditions between 5 blocks
16	if(k==0&&l==0) { … } … if(k==0&&l==5) { … }
17	}
Host Part
1	int main
2	{
3	// Initialize data on CPU
4	InitialData(double* cpuData, double* gpuData);
5	long int i, j, threadPos, … ;
6	// Copy data from CPU to GPU
7	cudaMemcpy(cpuData, gpuData, Datasize, cudaMemcpyHostToDevice);
8	// Set dimensions of the grid and thread
9	dim3 block(1024, 1); dim3 grid(M, N - 1);
10	// Start Kernel function
11	Nonlinear≪ grid, block≪(gpuData);
12	// Copy data from GPU to CPU
13	cudaMemcpy(gpuData, cpuData, Datasize, cudaMemcpyDeviceToHost);
14	return 0;
15	}

In the device part, there are two device functions which are called by the kernel function multiple times. These two device functions are formed to solve the singularity in the second integral of the boundary integral equation. The 16 special threads are set separately for fast computation. Other threads with the same CUDA code are used to complete the calculation of the remaining parts of the boundary integral equation. After threads have finished computing, all thread contributions are summed up, and $(\mathit{N}-\mathbf{1})\mathit{>M}$ nonlinear equations are built. Then arbitrarily choose 5 blocks to calculate the free surface condition and the radiation condition, $(\mathit{>N}+\mathbf{3})\mathit{>M}$ nonlinear equations can be obtained. Therefore, these equations are formed by using the $\mathbf{1024}\times (\mathit{>N}-\mathbf{1})\mathit{>M}$ threads on GPU.

In the host part, the environment variables are configurated firstly. Then the data is transferred from CPU to GPU and the parallel instruction is sent to GPU. Finally, CPU gets the computation results from GPU.

4.2.2. Building preconditioner matrix

The building preconditioner matrix is decomposed to several tasks that can be operated in parallel by corresponding kernel functions in GPU blocks. The parallel idea and program structure are roughly similar to the part of creating nonlinear equations.

In the device part, in order to avoid data storage conflicts in GPU, three kernel functions are used to construct the preconditioner matrix in turn. As mentioned above, the preconditioner matrix size is $\mathbf{2}(\mathit{N}+\mathbf{1})\mathit{M}\times \mathbf{2}(\mathit{N}+\mathbf{1})\mathit{M}$, and four submatrices are formed by block decomposition method and banded preconditioner method: Submatrix $\mathit{A}$ is a tridiagonal matrix of size $\mathbf{3}\times (\mathit{N}+\mathbf{1})\mathit{M}$; the submatrix $\mathit{B}$ and $\mathit{C}$ only differ between coefficients, and the base matrix ${\mathbf{B}}_{\mathbf{0}}$ can be constructed to represent them respectively with a size of $(\mathit{N}+\mathbf{1})\times (\mathit{N}+\mathbf{1})$, $\mathit{B}=\mathbf{1}/{\mathit{F}}^{\mathbf{2}}·{\mathbf{B}}_{\mathbf{0}}$ and $\mathit{C}=\mathbf{2}\mathit{\pi }·{\mathbf{B}}_{\mathbf{0}}$; the submatrix $\mathit{D}$ is a sparse matrix with a size of $\mathit{band}\times (\mathit{N}+\mathbf{1})\mathit{M}/\mathbf{3}$. Firstly, $\mathit{M}\times (\mathit{N}+\mathbf{1})$ thread blocks are called in the GPU to fulfill the parallel construction of the four submatrices by the kernel function precondition(), and the two device functions mentioned above are also used to eliminate the singularity of the linear boundary integral equation. Then the kernel function matrix() is written to call $\mathit{N}+\mathbf{1}$ thread blocks for parallel computation of ${\mathit{CA}}^{-\mathbf{1}}\mathit{B}$, which involves solving multiple right-handed linear systems and matrix multiplication. Finally, the subtraction operation between matrices is completed by kernel function Schur(), and $\mathit{M}\times (\mathit{N}+\mathbf{1})$ thread blocks are called to perform parallel operation of $\mathit{D}-{\mathit{CA}}^{-\mathbf{1}}\mathit{B}$. The programming for building preconditioner matrix on GPU by using CUDA language is briefly shown in Table 2.

In the host part, the variables are first defined according to calculation parameters, then data is transferred from the CPU to the GPU, then the dimension of the thread blocks and thread grid is specified, and finally the kernel functions precondition(), matrix(), and Schur() are successively released. This part of the host side code is similar to the establishment of nonlinear equations and will not be repeated here.

Table 2. The programming on GPU.

Device Part
Kernel functions
1	$\_\_$global$\_\_$void precondition(const long int N, const long int M, … )
2	{
3	// Calculate submatrix $\mathit{A}$
4	while(threadPos<N+3) { … };
5	if(blockIdx.x==0) { … };
6	// Calculate basis matrix ${\mathrm{B}}_{\mathbf{0}}$
7	if(blockIdx.x==0) { while(threadPos<N+1) { … } };
8	// Calculate submatrix $\mathbf{D}$
9	while(threadPos<N*M) { … };
10	// Sum up all thread contributions
11	for(i=blocDim/2;i>0;i=i/2){ … }
12	}
13	$\_\_$global$\_\_$void matrix(const long int N, const long int M, … )
14	{
15	// Calculate ${\mathbf{CA}}^{-\mathbf{1}}\mathit{B}$
16	double Bcoef = 1/($\mathit{F}\times \mathit{F}$);
17	double Ccoef = 2*3.1415;
18	if(blockIdx.x==0) { while(threadPos<N+1) { … } };
19	}
20	$\_\_$global$\_\_$void Schur(const long int N, const long int M,double*_b,double*D )
21	{
22	// Calculate $\mathbf{D}-{\mathbf{CA}}^{-\mathbf{1}}\mathit{B}$
23	int l=blockIdx.x;
24	int k=blockIdx.y;
25	int pos=1*(N+1)+k;
26	while(threadPos<N+1)
27	{
27	D[(1*(1+N)+threadPos)*M*(1+N)+pos]=D[(1*(1+N)+
28	threadPos)*M*(1+N)+pos]-_b[threadPos*(1+N)+k];
29	threadPos+=blockDim.x;
30	}
31	}

Table 2. The programming on GPU.

Device Part
Kernel functions
1	$\_\_$global$\_\_$void precondition(const long int N, const long int M, … )
2	{
3	// Calculate submatrix $\mathit{A}$
4	while(threadPos<N+3) { … };
5	if(blockIdx.x==0) { … };
6	// Calculate basis matrix ${\mathrm{B}}_{\mathbf{0}}$
7	if(blockIdx.x==0) { while(threadPos<N+1) { … } };
8	// Calculate submatrix $\mathbf{D}$
9	while(threadPos<N*M) { … };
10	// Sum up all thread contributions
11	for(i=blocDim/2;i>0;i=i/2){ … }
12	}
13	$\_\_$global$\_\_$void matrix(const long int N, const long int M, … )
14	{
15	// Calculate ${\mathbf{CA}}^{-\mathbf{1}}\mathit{B}$
16	double Bcoef = 1/($\mathit{F}\times \mathit{F}$);
17	double Ccoef = 2*3.1415;
18	if(blockIdx.x==0) { while(threadPos<N+1) { … } };
19	}
20	$\_\_$global$\_\_$void Schur(const long int N, const long int M,double*_b,double*D )
21	{
22	// Calculate $\mathbf{D}-{\mathbf{CA}}^{-\mathbf{1}}\mathit{B}$
23	int l=blockIdx.x;
24	int k=blockIdx.y;
25	int pos=1*(N+1)+k;
26	while(threadPos<N+1)
27	{
27	D[(1*(1+N)+threadPos)*M*(1+N)+pos]=D[(1*(1+N)+
28	threadPos)*M*(1+N)+pos]-_b[threadPos*(1+N)+k];
29	threadPos+=blockDim.x;
30	}
31	}

4.2.3. Inverting preconditioner matrix

Comparing the Math Kernel Library which is famous for the computation of sparse linear algebra, the cuSolverSP library is generally faster for solving sparse linear systems [22]. In this paper, the cuSolverSP library is adopted to invert preconditioner matrix. The present sparse linear system is special, the right-hand term of the system $\mathit{v}$ changes continuously in the iteration whereas the left-hand term does not. The characteristic of the sparse linear system suggests using QR factorization to calculate ${\mathit{P}}^{-\mathbf{1}}\mathit{v}$[19]. By QR factorization, the sparse matrix is decomposed into an orthogonal matrix and an upper triangular matrix, which are saved in GPU memory and are directly used to solve linear equations in each iteration. Finally, the preconditioner-vector products ${\mathit{P}}^{-\mathbf{1}}\mathit{v}$ can be obtained.

Step1: Using CSR data format to save preconditioner matrix with an appropriate bandwidth;

Step2: In the analysis stage, cusolverSpXcsrqrAnalysis() function is used to analyze the sparsity of orthogonal matrix and upper triangular matrix in QR decomposition. This process may consume a large amount of memory. If the memory is insufficient to complete the analysis, the program will stop running and return the corresponding error message;

Step3: In the preparation stage, cusolverSpXcsrqrAnalysis() function is used to select the appropriate computing space to prepare for QR decomposition. Here, two memory blocks are prepared in the GPU, one to store the orthogonal matrix and the upper triangular matrix, and the other to perform QR decomposition;

Step4: The cusolverSpDcsrqrSetup() function is called to allocate storage space for the orthogonal and upper triangular matrices based on the results of the preparation stage. Then, cusolverSpDcsrqrFactor() function is used to complete the QR decomposition of coefficient matrix outside the cycle;

Step5: Using cusolverSpDcsrqrZeroPivot() function checks the singularity of the decomposition results, if the nearly singular the program terminates operation and error is given, return to step 1 to choose the bandwidth again;

Step6: In the loop body, the cusolverSpDcsrqrSolve() function is repeatedly called, and the solution of linear equations can be obtained directly by using the decomposition results stored in GPU;

The main CUDA functions are shown in Table 3.

Table 3. The list of CUDA functions for QR factorization.

No.	Function name	Goal
1	cusolverSpXcsrqrAnalysisHost();	Analyze structure
2	cusolverSpDcsrqrBufferInfoHost();	Set up workspace
3	cusolverSpDcsrqrSetupHost();	QR factorization
4	cusolverSpDcsrqrFactorHost();	QR factorization
5	cusolverSpDcsrqrZeroPivotHost();	Check singular
6	cusolverSpDcsrqrSolveHost();	Solve system

Table 3. The list of CUDA functions for QR factorization.

No.	Function name	Goal
1	cusolverSpXcsrqrAnalysisHost();	Analyze structure
2	cusolverSpDcsrqrBufferInfoHost();	Set up workspace
3	cusolverSpDcsrqrSetupHost();	QR factorization
4	cusolverSpDcsrqrFactorHost();	QR factorization
5	cusolverSpDcsrqrZeroPivotHost();	Check singular
6	cusolverSpDcsrqrSolveHost();	Solve system

4.2.4. Solving linear equations by GMRES algorithm

In the process of solving linear equations, because the matrix free idea is adopted to avoid the storage of coefficient matrix, there is no product operation of coefficient matrix and vector in GMRES algorithm, so the operations that can be parallel in this part are operations between vectors. Therefore, this paper mainly uses cuBLAS library to complete the CUDA programming of GMRES algorithm to solve linear equations.

The cuBlasDdot() function is used to realize the inner product of vectors in the GMRES algorithm; the vector subtraction is calculated using cublasDaxpy() function; cublasDnrm2() function is used to calculate the Euclidean norm of the vector; cublasDscal() function is used to divide vector and scalar. After obtaining the orthonormal basis of Krylov subspace and the upper Hessnberg matrix, cublasDrotg() function is used to perform Givens rotation transformation on the upper Hessnberg matrix in GPU device to obtain the upper triangular matrix. Then the solution of linear least squares problem in GMRES algorithm is obtained, and cublasDspmv() function is used to achieve orthonormal basis and vector multiplication to get the solution of linear equations.

5. Numerical Simulations and Discussion

In this section, numerical simualtion of ship wave in multiple cases are carried out using the CPU and GPU solvers, and the simulation results are discussed. The effectiveness of the developed banded preconditioner JFNK method is first verified. Then, comparisons between the proposed GPU solver and the CPU solver on accuracy and efficiency are performed. Finally, verifying the capability of GPU solver by comparing simulation results with real ship wakes. The parameters of the computing environment are listed in Table 4.

Table 4. The computing environment of high-performance computing cluster.

	CPU	GPU
Card	Intel Xeon Bronze 3204	NVIDIA Tesla A100
Memory	64 GB	40 GB
Max Cores	6 per node	6912
Programming language	C++	CUDA, C++

Table 4. The computing environment of high-performance computing cluster.

	CPU	GPU
Card	Intel Xeon Bronze 3204	NVIDIA Tesla A100
Memory	64 GB	40 GB
Max Cores	6 per node	6912
Programming language	C++	CUDA, C++

5.1. Verification of the banded preconditioner JFNK method

To reveal the effectiveness of the banded preconditioner method, numerical simulations on different mesh sizes,i.e.,$\mathbf{181}\times \mathbf{61}$, $\mathbf{241}\times \mathbf{81}$, $\mathbf{301}\times \mathbf{101}$ and $\mathbf{361}\times \mathbf{121}$ with $\mathbf{\Delta }\mathit{x}=\mathbf{0}.\mathbf{3},\mathbf{\Delta }\mathit{y}=\mathbf{0}.\mathbf{3}$ are carried out.

The overall runtimes against bandwidth on these four mesh sizes are illustrated in Figure 10. An optimal value of bandwidth ${\mathit{b}}^{\mathbf{\prime }}$ exists for a certain mesh size. Furthermore, the optimal value of ${\mathit{b}}^{\mathbf{\prime }}$ increases with the mesh size and approximately equals $\frac{\mathit{M}}{\mathbf{3}}$, in which $\mathit{M}$ means the number of latitude lines of the mesh. Therefore, the optimal bandwidth can be set to $\frac{\mathit{M}}{\mathbf{3}}\times (\mathit{N}+\mathbf{1})$ to get an optimal efficiency.

According to the optimal bandwidth selection rule, the running memories against the bandwidth are shown in Table 5. Correspondingly, the required running memory is drastically reduced by applying the banded preconditioner JFNK method. The mean reduction ratio is about 3.2, this means that the banded preconditioner JFNK method can save running memory by at least two-thirds.

Figure 10. Optimal values of bandwidth ${\mathit{b}}^{\prime }$ for different mesh sizes.

Figure 10. Optimal values of bandwidth ${\mathit{b}}^{\prime }$ for different mesh sizes.

Table 5. The running memory usage before and after applying the banded preconditioner method.

Mesh size	Before	${\mathit{b}}^{\prime }$	After	Reduction ratio
$\mathbf{181}\times \mathbf{61}$	0.91GB	19	0.28GB	3.2
$\mathbf{241}\times \mathbf{81}$	2.9GB	24	0.88GB	3.3
$\mathbf{301}\times \mathbf{101}$	6.9GB	33	2.3GB	3.0
$\mathbf{361}\times \mathbf{121}$	15GB	38	4.6GB	3.3

Table 5. The running memory usage before and after applying the banded preconditioner method.

Mesh size	Before	${\mathit{b}}^{\prime }$	After	Reduction ratio
$\mathbf{181}\times \mathbf{61}$	0.91GB	19	0.28GB	3.2
$\mathbf{241}\times \mathbf{81}$	2.9GB	24	0.88GB	3.3
$\mathbf{301}\times \mathbf{101}$	6.9GB	33	2.3GB	3.0
$\mathbf{361}\times \mathbf{121}$	15GB	38	4.6GB	3.3

5.2. Verification of the GPU solver

5.2.1. Accuracy

To verify the accuracy of the GPU solver, numerical simulations are conducted on $\mathit{F}=\mathbf{0}.\mathbf{7}$ and $\mathit{ϵ}=\mathbf{0}.\mathbf{4}$ with a $\mathbf{361}\times \mathbf{121}$ mesh and $\mathbf{\Delta }\mathit{x}=\mathbf{0}.\mathbf{3},\mathbf{\Delta }\mathit{y}=\mathbf{0}.\mathbf{3}$. The simulated wave heights on the centerline are compared with that of the CPU solver proposed by Sun [12], which is shown in Figure 11. The almost all points in the figure are traversed through the center by line, indicating that the calculation results of the GPU solver are very consistent with those of the CPU solver.

Furthermore, the MSE is used to further explain the error between them, as follows:

$$\begin{array}{cc}\hfill & \mathit{MSE}=\frac{\mathbf{1}}{\mathit{n}}\sum _{\mathbf{i}=\mathbf{1}}^{\mathit{n}}{({\mathit{Truch}}_{\mathbf{i}}-{\mathit{Value}}_{\mathit{i}})}^{\mathbf{2}}\hfill \end{array}$$

(23)

where $\mathit{n}$ is the amount of data, ${\mathit{Truch}}_{\mathit{i}}$ and ${\mathit{Value}}_{\mathit{i}}$ represent CPU results and GPU results respectively. According to Eqs.(23), the calculated MSE is 9.37E-8, indicating that the calculation error between the GPU and CPU solver is minimal. Since the CPU solver has been verified by Sun [12], the accuracy of the proposed GPU solver can also be acceptable.

Figure 11. A comparison of the centerline profiles for the simulation results of CPU solver and GPU solver, computed on a $\mathbf{361}\times \mathbf{121}$ mesh with $\mathbf{\Delta }\mathit{x}=\mathbf{0}.\mathbf{3},\mathbf{\Delta }\mathit{y}=\mathbf{0}.\mathbf{3},\mathit{F}=\mathbf{0}.\mathbf{7}$ and $\mathit{ϵ}=\mathbf{0}.\mathbf{4}$. The solid line represents the simulation result of the GPU solver, the solid circles represent the simulation result of the CPU solver.

Figure 11. A comparison of the centerline profiles for the simulation results of CPU solver and GPU solver, computed on a $\mathbf{361}\times \mathbf{121}$ mesh with $\mathbf{\Delta }\mathit{x}=\mathbf{0}.\mathbf{3},\mathbf{\Delta }\mathit{y}=\mathbf{0}.\mathbf{3},\mathit{F}=\mathbf{0}.\mathbf{7}$ and $\mathit{ϵ}=\mathbf{0}.\mathbf{4}$. The solid line represents the simulation result of the GPU solver, the solid circles represent the simulation result of the CPU solver.

5.2.2. Efficiency

To verify the efficiency of the GPU solver, numerical simulation are conducted on $\mathit{F}=\mathbf{0}.\mathbf{7}$ and $\mathit{ϵ}=\mathbf{0}.\mathbf{4}$ with five mesh sizes, namely $\mathbf{121}\times \mathbf{41}$, $\mathbf{181}\times \mathbf{61}$, $\mathbf{241}\times \mathbf{81}$, $\mathbf{301}\times \mathbf{101}$, $\mathbf{361}\times \mathbf{121}$ and $\mathbf{\Delta }\mathit{x}=\mathbf{0}.\mathbf{3},\mathbf{\Delta }\mathit{y}=\mathbf{0}.\mathbf{3}$. The overall runtimes of GPU solver are compared with that of the CPU solver proposed by Sun [12], as shown in Figure 12. It clearly shows that the overall runtimes of the GPU solver are much shorter than that of the CPU solver. The clear accelerated-up ratios between the GPU solver and the CPU solver are shown in Table 6. The accelerated-up ratio on all cases are around 20.0. Therefore, the computaion efficiency of the proposed GPU solver is much higher than the CPU solver.

Figure 12. The runtime of the GPU solver and CPU solver at different mesh sizes, red bars represent the GPU solver results and blue bars represent the CPU solver results.

Figure 12. The runtime of the GPU solver and CPU solver at different mesh sizes, red bars represent the GPU solver results and blue bars represent the CPU solver results.

The proposed GPU solver has also been compared with another GPU solver proposed by Pethiyagoda [6] on these cases. The comparison of computation time between them is shown in Table 6. Obviously, the efficiency of the GPU solver proposed in this paper is higher than the GPU solver proposed by Pethiyagod [6], and the advantage is more significant with the increase of the mesh size. The reason is that Pethiyagoda [6] only introduced the GPU acceleration technique in the boundary integral method not the whole processs, whereas this paper proposes a complete parallel computing framework including the parallel process for inverting preconditioner matrix. In the process of inverting, the larger mesh size the heavier computational load, hence the advantage of the GPU solver proposed in this paper can be more significant.

Table 6. The comparisons of runtime between CPU solver and GPU solvers on different mesh sizes, the CPU solver is proposed by Sun [12], one GPU solver is proposed in this paper, the other GPU solver is proposed by Pethiyagoda [6](results of Exp. ).

Mesh size	CPU solver	Exp.	GPU solver	Accelerated-up ratio
$\mathbf{121}\times \mathbf{41}$	8.96E+1 s	1.61E+1 s	5.60E+0 s	16.1
$\mathbf{181}\times \mathbf{61}$	2.70E+2 s	1.22E+2 s	1.13E+1 s	23.9
$\mathbf{241}\times \mathbf{81}$	8.70E+2 s	5.51E+2 s	5.37E+1 s	16.0
$\mathbf{301}\times \mathbf{101}$	2.56E+3 s	1.78E+3 s	1.32E+2 s	19.3
$\mathbf{361}\times \mathbf{121}$	6.08E+3 s	5.04E+3 s	2.95E+2 s	20.6

Table 6. The comparisons of runtime between CPU solver and GPU solvers on different mesh sizes, the CPU solver is proposed by Sun [12], one GPU solver is proposed in this paper, the other GPU solver is proposed by Pethiyagoda [6](results of Exp. ).

Mesh size	CPU solver	Exp.	GPU solver	Accelerated-up ratio
$\mathbf{121}\times \mathbf{41}$	8.96E+1 s	1.61E+1 s	5.60E+0 s	16.1
$\mathbf{181}\times \mathbf{61}$	2.70E+2 s	1.22E+2 s	1.13E+1 s	23.9
$\mathbf{241}\times \mathbf{81}$	8.70E+2 s	5.51E+2 s	5.37E+1 s	16.0
$\mathbf{301}\times \mathbf{101}$	2.56E+3 s	1.78E+3 s	1.32E+2 s	19.3
$\mathbf{361}\times \mathbf{121}$	6.08E+3 s	5.04E+3 s	2.95E+2 s	20.6

5.2.3. Capability

In the proposed GPU solver, like the CPU solver [12], three parameters can be used to regulate simulation results, namely source strength, source type and Froude number. The wake characteristics can be regulated by adjusting these parameters. For example, the wake waves of vessels of high speed and small overall length can be generated by using Rankine source, lower source strength and larger Froude number; oppositely, it is appropriate to choose small Froude number and a higher strength Kelvin source. As shown in Figure 13, the simulation patterns of GPU solver are compared with the real ship wave patterns. It is found that the simulation patterns are consistent with the real ship waves, the proposed GPU solver can also generate high quality simulation patterns for 3-D nonlinear ship wave.

Figure 13. Wake pattern of real ship waves and the GPU solver computational results. The picture of real speedboat wake pattern comes from internet https://www.quanjing.com, accessed on 1 september 2023; the picture of real fishing ship wake comes from https://www.shutterstock.com, accessed on 1 september 2023; the picture of real large vessel wake comes from internet https://blogs.worldbank.org, accessed on 6 september 2023.

Figure 13. Wake pattern of real ship waves and the GPU solver computational results. The picture of real speedboat wake pattern comes from internet https://www.quanjing.com, accessed on 1 september 2023; the picture of real fishing ship wake comes from https://www.shutterstock.com, accessed on 1 september 2023; the picture of real large vessel wake comes from internet https://blogs.worldbank.org, accessed on 6 september 2023.

6. Conclusions

The numerical simulation of ship wave is important for practical ocean engineering. This paper proposes a highly-paralleled numerical scheme for simulating three-dimensional (3-D) nonlinear Kelvin ship waves effectively, including a numerical model for nonlinear ship waves, a banded preconditioner JFNK method and a GPU based parallel computing framework. Numerical simulations show that the proposed GPU solver can save GPU memory and obtain high efficiency significantly. This highly-paralleled numerical scheme provides an opportunity for the further study of the nonlinear Kelvin ship waves on a large scale.

(1)

The bandwidth has an effect on the running memory and runtime of the GPU solver. Based on the mesh size, the value of the most appropriate bandwidth is around $\frac{\mathit{M}}{\mathbf{3}}\times (\mathit{N}+\mathbf{1})$, with more than 66% GPU memory can be saved.

(2)

The GPU solver can obtain an accurate numerical solution. The mean square error of GPU solver results and CPU solver results is $\mathit{MSE}$=9.37E-8, which is acceptable.

(3)

By designing the GPU parallel computing framework, the computation of ship wave simulation is accelerated up to 20 times.

Although an highly-paralleled numerical scheme for nonlinear ship wave is proposed in this paper, some assumptions are still made in the construction of the numerical model, such as infinite water depth and the steady motion of ship on calm water. It is of great significance to improve simulation results by further exploring the influence of finite water depth, tangential flow and unsteady ship motion on nonlinear ship waves.