1. Introduction

2. Numerical Model

Figure 1. Flow field diagram.

Figure 1. Flow field diagram.

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

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.

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.

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.

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

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

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.

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. ).

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

(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.

Abbreviations

References

1. Dias, F. Ship Waves and Kelvin. j.Fluid Mech. 2014, 746, 1–4. [Google Scholar] [CrossRef]
2. Tuck, E.; Scullen, D. A comparison of linear and nonlinear computations of waves made by slender submerged bodies. J Eng Math 2002, 42, 255–264. [Google Scholar] [CrossRef]
3. Froude, w. Experiments upon the effect produced on the wave-making resistance of ships by length of parallel middle body; Institution of Naval Architects, 1877.
4. Kelvin, L. On Ship Waves. Proc. Inst. Mech. Engrs. 1887, 38, 409–434. [Google Scholar] [CrossRef]
5. Rabaud, M.; Moisy, F. Ship Wakes: Kelvin or Mach Angle? Phys. Rev. Lett 2013, 110, 214503.1–214503.5. [Google Scholar] [CrossRef] [PubMed]
6. Pethiyagoda, R.; Moroney, T.; Lustri, C.; McCue, S. Kelvin Wake Pattern at Small Froude Numbers. J. Fluid Mech. 2021, 915, A126. [Google Scholar] [CrossRef]
7. Ma, C.; Zhu, Y.; Wu, H.; He, J.; Zhang, C.; Li, W.; Noblesse, F. Wavelengths of the Highest Waves Created by Fast Monohull Ships or Catamarans. Ocean Eng. 2016, 113, 208–214. [Google Scholar] [CrossRef]
8. Havelock, T. Wave resistance: Some cases of three-dimensional fluid motion. Proc. R. Soc. London. Ser. A, Contain. Pap. a Math. Phys. 1919, 95, 354–365. [Google Scholar] [CrossRef]
9. J.H., M. The wave resistance of a ship. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 1898, 45, 106–123. [Google Scholar] [CrossRef]
10. Forbes, L. An algorithm for 3-dimensional free-surface problems in hydrodynamics. J. Comput. Phys. 1989, 82, 330–347. [Google Scholar] [CrossRef]
11. Parau, E.; Vanden-Broeck, J. Three-dimensional waves beneath an ice sheet due to a steadily moving pressure. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 2011, 369, 2973–2988. [Google Scholar] [CrossRef]
12. Sun, X.; Cai, M.; Wang, J.; Liu, C. Numerical Simulation of the Kelvin Wake Patterns. Appl. Sci. 2022, 12, 6265. [Google Scholar] [CrossRef]
13. Hori, C.; Gotoh, H.; Ikari, H.; Khayyer, A. GPU-Acceleration for Moving Particle Semi-Implicit Method. Comput. Fluids. 2011, 51, 174–183. [Google Scholar] [CrossRef]
14. Lu, X.; Dao, M.H.; Le, Q.T. A GPU-accelerated domain decomposition method for numerical analysis of nonlinear waves-current-structure interactions. Ocean Eng. 2022, 259, 111901. [Google Scholar] [CrossRef]
15. Xie, F.; Zhao, W.; Wan, D. CFD Simulations of Three-Dimensional Violent Sloshing Flows in Tanks Based on MPS and GPU. J. Hydrodyn. 2020, 32, 672–683. [Google Scholar] [CrossRef]
16. Saad, Y.; Schultz, M. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM J. Sci. Stat. Comput. 1986, 7, 856–869. [Google Scholar] [CrossRef]
17. Brown, P.; Saad, Y. Hybrid Krylov Methods for Nonlinear Systems of Equations. SIAM J. Sci. Stat. Comput. 1990, 11, 450–481. [Google Scholar] [CrossRef]
18. Dembo, R.; Eisenstat, S.; Steihaug, T. Inexact Newton Methods. SIAM J. Numer. Anal. 1982, 19, 400–408. [Google Scholar] [CrossRef]
19. Trefethen, L.; Bau, D. Numerical linear algebra; Siam, 1997.
20. Lustri, C.J.; Chapman, S.J. Steady Gravity Waves Due to a Submerged Source. J. Fluid Mech. 2013, 732, 400–408. [Google Scholar] [CrossRef]
21. NVIDIA. CUDA Toolkit Documentation v11.7.1., 2022.
22. Grossman, M.; Mckercher, T. Professional CUDA C programming; China Machine Press: Beijing, China, 2017. [Google Scholar]