%  LaTeX support: latex@mdpi.com 
%  For support, please attach all files needed for compiling as well as the log file, and specify your operating system, LaTeX version, and LaTeX editor.

%=================================================================
\documentclass[journal,article,submit,pdftex,moreauthors]{Definitions/mdpi} 

 \makeatletter
 \let\c@lofdepth\relax
 \let\c@lotdepth\relax
 \makeatother
\usepackage{subfigure}

\usepackage{float}
\usepackage{makecell}
\usepackage{caption}
\usepackage{tabularx}

%\usepackage[authoryear]{natbib}
%\usepackage[round, sort]{natbib}
%--------------------
% Class Options:
%--------------------
%----------
% journal
%----------
% Choose between the following MDPI journals:
% acoustics, actuators, addictions, admsci, adolescents, aerobiology, aerospace, agriculture, agriengineering, agrochemicals, agronomy, ai, air, algorithms, allergies, alloys, analytica, analytics, anatomia, animals, antibiotics, antibodies, antioxidants, applbiosci, appliedchem, appliedmath, applmech, applmicrobiol, applnano, applsci, aquacj, architecture, arm, arthropoda, arts, asc, asi, astronomy, atmosphere, atoms, audiolres, automation, axioms, bacteria, batteries, bdcc, behavsci, beverages, biochem, bioengineering, biologics, biology, biomass, biomechanics, biomed, biomedicines, biomedinformatics, biomimetics, biomolecules, biophysica, biosensors, biotech, birds, bloods, blsf, brainsci, breath, buildings, businesses, cancers, carbon, cardiogenetics, catalysts, cells, ceramics, challenges, chemengineering, chemistry, chemosensors, chemproc, children, chips, cimb, civileng, cleantechnol, climate, clinpract, clockssleep, cmd, coasts, coatings, colloids, colorants, commodities, compounds, computation, computers, condensedmatter, conservation, constrmater, cosmetics, covid, crops, cryptography, crystals, csmf, ctn, curroncol, cyber, dairy, data, ddc, dentistry, dermato, dermatopathology, designs, devices, diabetology, diagnostics, dietetics, digital, disabilities, diseases, diversity, dna, drones, dynamics, earth, ebj, ecologies, econometrics, economies, education, ejihpe, electricity, electrochem, electronicmat, electronics, encyclopedia, endocrines, energies, eng, engproc, entomology, entropy, environments, environsciproc, epidemiologia, epigenomes, est, fermentation, fibers, fintech, fire, fishes, fluids, foods, forecasting, forensicsci, forests, foundations, fractalfract, fuels, future, futureinternet, futurepharmacol, futurephys, futuretransp, galaxies, games, gases, gastroent, gastrointestdisord, gels, genealogy, genes, geographies, geohazards, geomatics, geosciences, geotechnics, geriatrics, grasses, gucdd, hazardousmatters, healthcare, hearts, hemato, hematolrep, heritage, higheredu, highthroughput, histories, horticulturae, hospitals, humanities, humans, hydrobiology, hydrogen, hydrology, hygiene, idr, ijerph, ijfs, ijgi, ijms, ijns, ijpb, ijtm, ijtpp, ime, immuno, informatics, information, infrastructures, inorganics, insects, instruments, inventions, iot, j, jal, jcdd, jcm, jcp, jcs, jcto, jdb, jeta, jfb, jfmk, jimaging, jintelligence, jlpea, jmmp, jmp, jmse, jne, jnt, jof, joitmc, jor, journalmedia, jox, jpm, jrfm, jsan, jtaer, jvd, jzbg, kidneydial, kinasesphosphatases, knowledge, land, languages, laws, life, liquids, literature, livers, logics, logistics, lubricants, lymphatics, machines, macromol, magnetism, magnetochemistry, make, marinedrugs, materials, materproc, mathematics, mca, measurements, medicina, medicines, medsci, membranes, merits, metabolites, metals, meteorology, methane, metrology, micro, microarrays, microbiolres, micromachines, microorganisms, microplastics, minerals, mining, modelling, molbank, molecules, mps, msf, mti, muscles, nanoenergyadv, nanomanufacturing,\gdef\@continuouspages{yes}} nanomaterials, ncrna, ndt, network, neuroglia, neurolint, neurosci, nitrogen, notspecified, %%nri, nursrep, nutraceuticals, nutrients, obesities, oceans, ohbm, onco, %oncopathology, optics, oral, organics, organoids, osteology, oxygen, parasites, parasitologia, particles, pathogens, pathophysiology, pediatrrep, pharmaceuticals, pharmaceutics, pharmacoepidemiology,\gdef\@ISSN{2813-0618}\gdef\@continuous pharmacy, philosophies, photochem, photonics, phycology, physchem, physics, physiologia, plants, plasma, platforms, pollutants, polymers, polysaccharides, poultry, powders, preprints, proceedings, processes, prosthesis, proteomes, psf, psych, psychiatryint, psychoactives, publications, quantumrep, quaternary, qubs, radiation, reactions, receptors, recycling, regeneration, religions, remotesensing, reports, reprodmed, resources, rheumato, risks, robotics, ruminants, safety, sci, scipharm, sclerosis, seeds, sensors, separations, sexes, signals, sinusitis, skins, smartcities, sna, societies, socsci, software, soilsystems, solar, solids, spectroscj, sports, standards, stats, std, stresses, surfaces, surgeries, suschem, sustainability, symmetry, synbio, systems, targets, taxonomy, technologies, telecom, test, textiles, thalassrep, thermo, tomography, tourismhosp, toxics, toxins, transplantology, transportation, traumacare, traumas, tropicalmed, universe, urbansci, uro, vaccines, vehicles, venereology, vetsci, vibration, virtualworlds, viruses, vision, waste, water, wem, wevj, wind, women, world, youth, zoonoticdis 
% For posting an early version of this manuscript as a preprint, you may use "preprints" as the journal. Changing "submit" to "accept" before posting will remove line numbers.

%---------
% article
%---------
% The default type of manuscript is "article", but can be replaced by: 
% abstract, addendum, article, book, bookreview, briefreport, casereport, comment, commentary, communication, conferenceproceedings, correction, conferencereport, entry, expressionofconcern, extendedabstract, datadescriptor, editorial, essay, erratum, hypothesis, interestingimage, obituary, opinion, projectreport, reply, retraction, review, perspective, protocol, shortnote, studyprotocol, systematicreview, supfile, technicalnote, viewpoint, guidelines, registeredreport, tutorial
% supfile = supplementary materials

%----------
% submit
%----------
% The class option "submit" will be changed to "accept" by the Editorial Office when the paper is accepted. This will only make changes to the frontpage (e.g., the logo of the journal will get visible), the headings, and the copyright information. Also, line numbering will be removed. Journal info and pagination for accepted papers will also be assigned by the Editorial Office.

%------------------
% moreauthors
%------------------
% If there is only one author the class option oneauthor should be used. Otherwise use the class option moreauthors.

%---------
% pdftex
%---------
% The option pdftex is for use with pdfLaTeX. Remove "pdftex" for (1) compiling with LaTeX & dvi2pdf (if eps figures are used) or for (2) compiling with XeLaTeX.

%=================================================================
% MDPI internal commands - do not modify
\firstpage{1} 
\makeatletter 
\setcounter{page}{\@firstpage} 
\makeatother
\pubvolume{1}
\issuenum{1}
\articlenumber{0}
\pubyear{2023}
\copyrightyear{2023}
%\externaleditor{Academic Editor: Firstname Lastname}
\datereceived{ } 
\daterevised{ } % Comment out if no revised date
\dateaccepted{ } 
\datepublished{ } 
%\datecorrected{} % For corrected papers: "Corrected: XXX" date in the original paper.
%\dateretracted{} % For corrected papers: "Retracted: XXX" date in the original paper.
\hreflink{https://doi.org/} % If needed use \linebreak
%\doinum{}
%\pdfoutput=1 % Uncommented for upload to arXiv.org

%=================================================================
% Add packages and commands here. The following packages are loaded in our class file: fontenc, inputenc, calc, indentfirst, fancyhdr, graphicx, epstopdf, lastpage, ifthen, float, amsmath, amssymb, lineno, setspace, enumitem, mathpazo, booktabs, titlesec, etoolbox, tabto, xcolor, colortbl, soul, multirow, microtype, tikz, totcount, changepage, attrib, upgreek, array, tabularx, pbox, ragged2e, tocloft, marginnote, marginfix, enotez, amsthm, natbib, hyperref, cleveref, scrextend, url, geometry, newfloat, caption, draftwatermark, seqsplit
% cleveref: load \crefname definitions after \begin{document}

%=================================================================
% Please use the following mathematics environments: Theorem, Lemma, Corollary, Proposition, Characterization, Property, Problem, Example, ExamplesandDefinitions, Hypothesis, Remark, Definition, Notation, Assumption
%% For proofs, please use the proof environment (the amsthm package is loaded by the MDPI class).

%=================================================================
% Full title of the paper (Capitalized)
\Title{A GPU accelerated method for 3-D nonlinear Kelvin ship wave patterns simulation}

% MDPI internal command: Title for citation in the left column
\TitleCitation{Title}

% Author Orchid ID: enter ID or remove command
\newcommand{\orcidauthorA}{0000-0001-5926-0382} % Add \orcidA{} behind the author's name
%\newcommand{\orcidauthorB}{0000-0000-0000-000X} % Add \orcidB{} behind the author's name

% Authors, for the paper (add full first names)
\Author{Xiaofeng Sun* $^{1}$\orcidA{}, Miaoyu Cai $^{2}$ and Junchen Ding $^{1}$ }

%\longauthorlist{yes}

% MDPI internal command: Authors, for metadata in PDF
\AuthorNames{Xiaofeng Sun, Miaoyu Cai and Junchen Ding}

% MDPI internal command: Authors, for citation in the left column
\AuthorCitation{Sun, X.; Cai, M.; Ding, J.}
% If this is a Chicago style journal: Lastname, Firstname, Firstname Lastname, and Firstname Lastname.

% Affiliations / Addresses (Add [1] after \address if there is only one affiliation.)
\address{%
$^{1}$ \quad Navigation College, Dalian Maritime University, No.1 Linghai Road, Dalian 116026, China;xfsun\_dlmu@163.com; 2398719038@qq.com\\
$^{2}$ \quad School of Naval Architecture, Ocean \& Civil Engineering, Shanghai Jiao Tong University, No.800 Dongchuan Road, Shanghai 200000, China; caimiaoyu@outlook.com}


% Contact information of the corresponding author
\corres{Correspondence: xfsun\_dlmu@163.com}

% Current address and/or shared authorship
%\firstnote{Current address: Affiliation 3.} 
%\secondnote{These authors contributed equally to this work.}
% The commands \thirdnote{} till \eighthnote{} are available for further notes

%\simplesumm{} % Simple summary

%\conference{} % An extended version of a conference paper

% Abstract (Do not insert blank lines, i.e. \\) 
\abstract{Ship wave is of great interest for wave drag, coastal erosion and ship detection.
In this paper, a highly-paralleled numerical scheme is proposed for simulating three-dimensional (3-D) nonlinear Kelvin ship waves effectively. 
First, a numerical model for nonlinear ship waves is established based on potential flow theory, boundary integral method and Jacobian-free Newton-Krylov (JFNK) method. 
In order to improve computational efficiency and reduce data storage of JFNK method, a banded preconditioner method is then developed by formulating the optimal bandwidth selection rule. 
After that, a Graphics Process Unit (GPU) based parallel computing framework is designed,  and a GPU solver is developed by using Compute Unified Device Architecture (CUDA) language.
Finally, numerical simulations of 3-D nonlinear ship waves under multiple scales are performed by using the GPU and CPU solvers.
Simulation results show that the proposed GPU solver is more efficient than the CPU solver with the same accuracy. More than 66\% GPU memory can be saved and the computational speed can be accelerated up to 20 times. 
Hence, the computation time for Kelvin ship waves simulation can be significantly reduced by applying the GPU parallel numerical scheme, which lays a solid foundation for practical ocean engineering.}

% Keywords
\keyword{Kelvin wake pattern;GPU acceleration;Boundary integral method;JFNK method;Banded preconditioner method}  

% The fields PACS, MSC, and JEL may be left empty or commented out if not applicable
%\PACS{J0101}
%\MSC{}
%\JEL{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Diversity
%\LSID{\url{http://}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Applied Sciences
%\featuredapplication{Authors are encouraged to provide a concise description of the specific application or a potential application of the work. This section is not mandatory.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Data
%\dataset{DOI number or link to the deposited data set if the data set is published separately. If the data set shall be published as a supplement to this paper, this field will be filled by the journal editors. In this case, please submit the data set as a supplement.}
%\datasetlicense{License under which the data set is made available (CC0, CC-BY, CC-BY-SA, CC-BY-NC, etc.)}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Toxins
%\keycontribution{The breakthroughs or highlights of the manuscript. Authors can write one or two sentences to describe the most important part of the paper.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Encyclopedia
%\encyclopediadef{For entry manuscripts only: please provide a brief overview of the entry title instead of an abstract.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Advances in Respiratory Medicine
%\addhighlights{yes}
%\renewcommand{\addhighlights}{%

%\noindent This is an obligatory section in “Advances in Respiratory Medicine”, whose goal is to increase the discoverability and readability of the article via search engines and other scholars. Highlights should not be a copy of the abstract, but a simple text allowing the reader to quickly and simplified find out what the article is about and what can be cited from it. Each of these parts should be devoted up to 2~bullet points.\vspace{3pt}\\
%\textbf{What are the main findings?}
% \begin{itemize}[labelsep=2.5mm,topsep=-3pt]
% \item First bullet.
% \item Second bullet.
% \end{itemize}\vspace{3pt}
%\textbf{What is the implication of the main finding?}
% \begin{itemize}[labelsep=2.5mm,topsep=-3pt]
% \item First bullet.
% \item Second bullet.
% \end{itemize}
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \setcounter{section}{-1} %% Remove this when starting to work on the template.
% \section{How to Use this Template}

% The template details the sections that can be used in a manuscript. Note that the order and names of article sections may differ from the requirements of the journal (e.g., the positioning of the Materials and Methods section). Please check the instructions on the authors' page of the journal to verify the correct order and names. For any questions, please contact the editorial office of the journal or support@mdpi.com. For LaTeX-related questions please contact latex@mdpi.com.%\endnote{This is an endnote.} % To use endnotes, please un-comment \printendnotes below (before References). Only journal Laws uses \footnote.

% The order of the section titles is different for some journals. Please refer to the "Instructions for Authors” on the journal homepage.

\section{Introduction}
This study is dedicated to studying highly-parallel algorithms for steady three-dimensional free surface profiles that are caused by a disturbance to a free stream. These profiles will appear stationary in the reference frame of the moving ship, referred to as “Kelvin ship waves” \citep{Dias2014}. Researches on Kelvin ship wave patterns have ongoing practical applications to ship hull design, ship detection and environmentally friendly shipping policies \citep{Michell1989,Sheremet2013,Wang2021,  Pethiyagoda2021, Luo2022}. 

Froude \cite{Froude1877}, a famous naval architect, first comprehensively described the morphology and main characteristics of ship waves. Under the assumption of infinite water depth, Kelvin \cite{Kelvin1887} 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 \cite{Rabaud2013} noted that for sufficiently fast-moving ships, contrary to commonly held views, the wake angle that is observed behind a steadily moving ship is less than the well-known Kelvin angle. This finding aroused the interest of many academics. Subsequently, various effect factors for the Kelvin wake form were discussed in plenty of papers, e.g., Froude number \citep{Pethiyagoda2021b, Verberck2013}, non-axisymmetric simplified ship models and interference effects \citep{Benzaquen2014, Miao2015, Ma2016, Zhu2017, Wu2019a}, shear current, submergence depth and finite depth \citep{Ellingsen2014, Li2016, Li2018, Wu2019b}, surface tension and the bottom topography \citep{Liang2018, Grue2017}, and viscosity \citep{Liang2019}, etc. Accordingly, the research method of ship waves has gradually shifted from the previous analytical algorithms to numerical simulation.

The mathematical analysis of ship wave patterns has a very long history, the overwhelming majority of which concerns linear theories. Havelock \cite{Havelock1919}  firstly presented the linear solution for the classic problem of flow past a pressure distribution. Such ideal perturbations can also be replaced by a single submerged point source singularity \cite{Havelock1932} and submerged bodies \cite{Tuck2002,Reed2002}. Another approach is to consider the ship wave pattern due to a thin ship \cite{Michell1898,Nobless2009}. With the development of numerical methods and computer technology, numerical simulation methods become more and more popular, and the study of ship waves had been extended from two-dimensional (2-D) linear problems to three-dimensional (3-D) nonlinear problems. Nowadays, there are three numerical methods widely used to solve surface wave problems: finite-difference method, finite-element method and boundary integral method \cite{Barbosa2005,Kan2022,Bettess1983,Ma2022,Bystricky2021,Wang2022,Gu2022}. In particular, Forbes \cite{Forbes1985, FORBES1989} apply boundary integral method to build a series of integro-differential equation by placing a mesh of N$\times$M grid points over the truncated plane, the full nonlinear free surface flow problem was solved with moderate efficiencies. In more recent times, many papers have applied Forbes’s formulation to solve fully 3-D nonlinear ship waves with meshes between $60\times20$ and $181\times61$ \cite{Parau2011,Sun2022,Cai2022}. And Pethiyagoda \cite{Pethiyagoda2016} noted that with less than 100 points used along the $x-$direction, the resolution over each wavelength is simply not of a sufficient standard for any claims about grid-independence to be made. 

With increasing mesh size, however, the computation time increases exponentially using only Central Processing Unit (CPU) computation power. Therefore, people are seeking a viable, cheap and portable alternative approach.  As the electric industry develops, the Graphics Processing Unit (GPU), a multi-processor designed to optimize for the execution of a massive number of threads, has become an alternative acceleration approach. Currently, the powerful GPU parallel computing ability has been used to improve the studies on ocean engineering. Hérault \cite{Herault2010} introduced the GPU acceleration technique in Smoothed Particle Hydrodynamics (SPH) to simulate the dam-break flow. Hori \cite{Hori2011} developed a GPU-based MPS code and achieved 7 times speedup in simulating 2-D dam-break flow with 110 thousand particles. Pethiyagoda \cite{Pethiyagoda2016}  introduced the GPU acceleration technique in  the boundary integral method to improve the computational efficiency of 3-D nonlinear ship wave problem. Xie \cite{Xie2020} employed an in-house solver MPSGPU-SJTU coupled with GPU techniques for the liquid sloshing simulation, to study the factors leading to the 3-D effect. LU \cite{LU2022} developed a GPU-accelerated high-order spectral solver to simulate the wave propagation and interaction with current in far field.

In this paper, a parallel solution framework based on GPU is designed 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 remainder of this paper is organized as follows. A brief introduction of the problem formulation is given in Section \ref{Section 2}. In Section \ref{Section 3}, the banded preconditioner JFNK algorithm is described. In Section \ref{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 \ref{Section 5}, and a summary in Section \ref{Section 6} concludes the paper.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Model}\label{Section 2}

This paper supposes that a flow is directed along the positive x-axis with uniform speed $U$. Considering the irrotational flow of an inviscid, incompressible fluid of infinite depth with ignoring the effects 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  \autoref{Fig: 1}. The transient waves will be generated on the free surface, on account of the disturbance caused by a source. Free surface wave height and flow field velocity potential can be expressed as $z=\zeta(x,y)$ and $\Phi(x,y,z)$. This paper is devoted to the steady-state problem that arises in the long time limit of this flow.
% Figure
\begin{figure}[htp]
	\centering
		\includegraphics[width=9 cm]{1.png}
	  \caption{Illustration of the disposition of the fluid.}\label{Fig: 1}
\end{figure}

Dimensionless analysis is performed with fluid velocity $U$ and distance $L$. The dimensionless 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. Then the problem is solved by the boundary integral method, the velocity potential can be regarded as a function simply of the two independent variables $x^\prime$ and $y^\prime$, and will be denoted by the symbol $\phi^\prime\left(x^\prime,y^\prime\right)=\Phi^\prime\left(x^\prime,y^\prime,\zeta^\prime\left(x^\prime,y^\prime\right)\right)$. The final boundary integral equation is described, as follows:

%\begin{adjustwidth}{-\extralength}{0cm}
\begin{flalign}
\label{Eq 8}
&2\pi(\phi '(Q)-x ')=\int_0^{\infty}\int_{-\infty}^{\infty}[\phi '(P)-\phi '(Q) -\rho '+x ']K^{(1)}(\rho ',\sigma ',x ',y ')d\sigma 'd\rho '&\notag\\
&+\int_0^{\infty}\int_{-\infty}^{\infty}\zeta '_{\rho}(P)K^{(2)}(\rho ',\sigma ',x ',y ')d\sigma 'd\rho '-\frac{\epsilon^\prime}{[{x^\prime}^2 +{x^\prime}^2 +{(\zeta '(Q)+1)}^2]^{\frac{1}{2}}}&
\end{flalign}
%\end{adjustwidth}

where the kernel functions $K^{\left(1\right)}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)$ and $K^{\left(2\right)}\left(\rho^\prime,\sigma^\prime, x^\prime,y^\prime\right)$ are described, as follows: 
%\begin{adjustwidth}{-\extralength}{0cm}
\begin{flalign*}
K^{\left(1\right)}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)=\frac{\zeta^\prime\left(P\right)-\zeta^\prime\left(Q\right)-\left(\rho^\prime-x^\prime\right){\zeta^\prime}_\rho\left(P\right)-\left(\sigma^\prime+y^\prime\right){\zeta^\prime}_\sigma\left(P\right)}{\left[\left(\rho^\prime-x^\prime\right)^2+\left(\sigma^\prime+y^\prime\right)^2+\left(\zeta^\prime\left(P\right)-\zeta^\prime\left(Q\right)\right)^2\right]^\frac{3}{2}}\notag\\
+\frac{\zeta^\prime\left(P\right)-\zeta^\prime\left(Q\right)-\left(\rho^\prime-x^\prime\right){\zeta^\prime}_\rho\left(P\right)-\left(\sigma^\prime-y^\prime\right){\zeta^\prime}_\sigma\left(P\right)}{\left[\left(\rho^\prime-x^\prime\right)^2+\left(\sigma^\prime-y^\prime\right)^2+\left(\zeta^\prime\left(P\right)-\zeta^\prime\left(Q\right)\right)^2\right]^\frac{3}{2}}
\end{flalign*}
%\end{adjustwidth}

%\begin{adjustwidth}{-\extralength}{0cm}
\begin{flalign*}
K^{\left(2\right)}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)=\frac{1}{\left[\left(\rho^\prime-x^\prime\right)^2+\left(\sigma^\prime+y^\prime\right)^2+\left(\zeta^\prime\left(P\right)-\zeta^\prime\left(Q\right)\right)^2\right]^\frac{3}{2}}\notag\\
+\frac{1}{\left[\left(\rho^\prime-x^\prime\right)^2+\left(\sigma^\prime-y^\prime\right)^2+\left(\zeta^\prime\left(P\right)-\zeta^\prime\left(Q\right)\right)^2\right]^\frac{3}{2}}
\end{flalign*}
%\end{adjustwidth}

There is a singularity in the second integral of boundary integral equation Eq. \eqref{Eq 8}, it can be solved as follows:
\begin{flalign}
\int_{y '_1}^{y '_M}\int_{x '_1}^{x '_N}{{[\zeta^\prime}_\rho\left(P\right)K^{\left(2\right)}}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)-{{\zeta^\prime}_x\left(Q\right)\ S}^{\left(2\right)}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)]d\sigma^\prime d\rho^\prime
+{\zeta^\prime}_x\left(Q\right)I^\prime\left(Q\right)
\end{flalign}
and
\begin{flalign}
&I^\prime\left(Q\right)=\int_{y^\prime_1}^{y^\prime_M}{\int_{{x^\prime}_1}^{x^\prime_ N}}{{S}^{\left(2\right)}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)d\sigma^\prime d\rho^\prime}&
\end{flalign}
where the kernel function $S^{\left(2\right)}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)$ can be described, as follows: 
\begin{flalign}
S^{\left(2\right)}\left(\rho^\prime,\sigma^\prime,x^\prime,y^\prime\right)&=\frac{1}{\sqrt{A^\prime{\left(\rho^\prime-x^\prime\right)}^2+B^\prime\left(\rho^\prime-x^\prime\right)\left(\sigma^\prime-y^\prime\right)+C^\prime{\left(\sigma^\prime-y^\prime\right)}^2}}&\notag\\
&+\frac{1}{\sqrt{A^\prime\left(\rho^\prime-x^\prime\right)^2-B^\prime\left(\rho^\prime-x^\prime\right)\left(\sigma^\prime+y^\prime\right)+C^\prime\left(\sigma^\prime+y^\prime\right)^2}}&\notag\\
A^\prime&=1+{\zeta^\prime_x}^2\left(Q\right)&\notag\\
B^\prime&=2{\zeta^\prime}_x^2\left(Q\right){\zeta^\prime_y}^2\left(Q\right)&\notag\\
C^\prime&=1+{\zeta^\prime_y}^2\left(Q\right)&
\end{flalign}

Now, the integral $I^\prime\left(Q\right)$ contains the singularity and is computed exactly in terms of logarithms:
\begin{flalign}
\label{Eq 16}
&\iint\frac{dsdt}{[As^2+Bst+Ct^2]}=\frac{t}{A^{\frac{1}{2}}}ln(2As+Bt+2[A(As^2+Bst+Ct^2)]^{\frac{1}{2}}&\notag\\
&+\frac{s}{C^{\frac{1}{2}}}ln(2Ct+Bs+2[C(As^2+Bst+Ct^2)]^{\frac{1}{2}}&\notag\\
&+h_1(s)+h_1(t)&
\end{flalign}

Moreover, the free surface conditions can be simplified by the symbol $\phi^\prime\left(x^\prime,y^\prime\right)$. Then the kinematic and dynamic boundary conditions of the free surface are combined to be
\begin{flalign}
\label{Eq 9}
&\frac{{(1+\zeta^\prime}_x^2){\phi^\prime}_y^2+\left(1+{\zeta^\prime}_y^2\right){\phi^\prime}_x^2-{2{\zeta^\prime}_x{\zeta^\prime}_y\phi^\prime}_x{\phi^\prime}_y}{{2(1+\zeta^\prime}_x^2+{\zeta^\prime}_y^2)}+\frac{\zeta^\prime}{{F^\prime}^2}=\frac{1}{2}&
\end{flalign}


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, \ldots, x'_N$ and  $y'_1,y'_2, \ldots, y'_N$ with intervals $\Delta x '$ and $\Delta y '$ in the $x$- and $y$- directions. The free surface elevation $\zeta^\prime\left(x^\prime,y^\prime\right)$ and the velocity potential $\phi^\prime\left(x^\prime,y^\prime\right)$ are represented by discrete point values $\zeta'_{k,l}$ and $\phi'_{k,l}$ at the nodes $(x^\prime_k, y^\prime_l), k=1,\ldots,N,l=1,\ldots,M$. And this paper chooses the $x$- derivatives of the functions $\phi^\prime$ and $\zeta^\prime$ as the basis for the solutions, together with the values of $\phi^\prime$ and $\zeta^\prime$ at the upstream boundary of the truncated domain, resulting in the vector of $2\left(N+1\right)M$ unknowns $\boldsymbol{u}$ to be

\begin{flalign}
\boldsymbol{u}=&[\phi '_{1,1},{(\phi '_x)}_{1,1}, \ldots ,\left(\phi^\prime_x\right)_{N,1}, \ldots ,
\phi^\prime_{1,M},\left(\phi^\prime_x\right)_{1,M},\ldots,\left(\phi^\prime_x\right)_{N,M}&\notag\\
&\zeta^\prime_{1,1},\left(\zeta^\prime_x\right)_{1,1},\ldots,\left(\zeta^\prime_x\right)_{N,1},\ldots,
\zeta^\prime_{1,M},\left(\zeta^\prime_x\right)_{1,M},\ldots,\left(\zeta^\prime_x\right)_{N,M}]^T.&
\end{flalign}

Therefore, The boundary integral equation can be discretized to the $(N-1)M$ nonlinear equations, and an additional $(N-1)M$ nonlinear equations are given by evaluating the free surface condition at the half mesh points. Moreover, $4M$ equations are provided by applying the radiation condition as follows. There are total $2(N+1)M$ equations to be solved for the $2(N+1)M$ unknowns in the nonlinear ship wave problem.
\begin{flalign}
\label{Eq 17}
&x_1^\prime((\phi_x^\prime)_{1,l}-1)+\gamma(\phi_{1,l}^\prime - x_1^\prime)=0&\\
\label{Eq 18}
&x_1^\prime((\phi_{xx}^\prime)_{1,l}-1)+\gamma((\phi_x^\prime)_{1,l}^\prime -1)=0&\\
\label{Eq 19}
&x_1^\prime(\zeta_x^\prime)_{1,l}+\gamma\zeta_{1,l}^\prime =0&\\
\label{Eq 20}
&x_1^\prime(\zeta_{xx}^\prime)_{1,l}+\gamma(\zeta_x^\prime)_{1,l}  =0&
\end{flalign}
where $\gamma$ is the decay coefficient.

Furthermore, more details about the governing equations, the boundary integral method and numerical discretization are provided by Sun\cite{Sun2022}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Banded Proconditioner JFNK Algorithm}\label{Section 3}
\subsection{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.   

After numerical discretization, a nonlinear system of equations could be obtained, as follows:
\begin{flalign}
&\boldsymbol{ F(u)}=0&
\end{flalign}
where $\boldsymbol u$ is the vector of unknowns of the length $2(N+1)M$. 

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{flalign}
\label{Eq 221}
&\boldsymbol{u_{t+1}} =\lambda_k \boldsymbol{\delta u_t+u_t }, \lambda_k\in(0,1]&
\end{flalign}

Its internal iteration is GMRES algorithm \citep{Saad1986}, which efficiently solves the correction in inexact Newton iteration, that is, computes large-scale linear equations as follows:
\begin{flalign}
\label{Eq 222}
&\boldsymbol{{\rm J}(u_t )\delta u_t=- F(u_t )} &
\end{flalign}
where $\boldsymbol{{\rm J}(u_t )= \partial F(u_t )/\partial u_t}$ is the Jacobian matrix \citep{Brown1990}.

%However, Newton’s method cannot solve large nonlinear systems efficiently. Therefore, based on Newton’s method, this paper adopts the Jacobian-free Newton-Krylov method with a preconditioner which is carefully constructed with entries taken from the Jacobian of the linearized problem. Eq. \eqref{Eq 22} is solved by using the iterative Generalized Minimum Residual (GMRES) algorithm \citep{Saad1986}, called the internal iteration. 
The GMRES method is one of the Krylov subspace methods which are attractive as linear solvers in the context of nonlinear Newton iteration. Since each iteration does not require the exact value of $\boldsymbol \delta u_t$ in solving nonlinear equations, the GMRES algorithm fits this need well and thus can increase total computation speed. 

Firstly, the approximate solution of $\boldsymbol \delta u_t$ is found by projecting obliquely onto the Krylov subspace
\begin{flalign}
& K_m \boldsymbol{({\rm J}_t {\rm P}^{-1}, F_t )}={\rm span}\boldsymbol{\{ F_t,{\rm J}_t {\rm P}^{-1}  F_t}, \ldots,\boldsymbol{({\rm J}_t {\rm P}^{-1} )^{m-1}  F_t\}}&
\end{flalign}
where $m$ is the value of the subspace dimension and the accuracy of the solution increases with the subspace dimension, $\boldsymbol{{\rm J}_t={\rm J}(u_t )}, \boldsymbol{ F_t= F(u_t )}$. The matrix $\boldsymbol{{\rm P}\approx{\rm J}_t}$ is the preconditioner matrix, whose purpose is to construct an approximation to the Jacobian $\boldsymbol{{\rm J}_t}$ which is cheap to form and to factorize. The calculation speed of the GMRES method can be significantly improved with a preconditioner matrix, because the spectrum of the preconditioned Jacobian $\boldsymbol{{\rm J}_t{\rm P}}^{-1}$ exhibits a clustering of eigenvalues \citep{Knoll2004}. 

An initial linear residual $\boldsymbol{{r}_0}$ is defined, given an initial guess $\boldsymbol u_0$, for the Newton correction,
\begin{flalign}
&  \boldsymbol{{ r}_0=-F(u_0)-{\rm J}_0{\rm P}^{-1}\delta u_0}&
\end{flalign}

Subsequently, the GMRES iteration minimizes $\left \| \boldsymbol{{ r}_t} \right \|$ to a suitable value. These Jacobian-vector products can be approximated by applying first-order difference quotients:
\begin{flalign}
\label{Eq 26}
&  \boldsymbol{{\rm J}_t{\rm P}^{-1}v }\approx \frac{\boldsymbol{F(u_t+h {\rm P^{-1}}v)-F(u_t)}}{h}&
\end{flalign}
where $v$ represents an arbitrary vector used in building the Krylov subspace \citep{Dembo1982}, and the $h$ is a small perturbation
\begin{flalign}
&  h=\frac{\sqrt{(1+\left \| \boldsymbol{{u}_t} \right \|)h_{mach}}}{\left \| \boldsymbol{v} \right \|}&
\end{flalign}

 Finally, as for the nonlinear ship wave problem, the initial guess $\boldsymbol{u_0}$ can be defined as below. The nonlinear equations are solved according to the calculation flow of baned preconditioner JFNK method, as shown in \autoref{Fig: JFNK}. Note that $\boldsymbol{v}$ in the figure is a unit orthogonal vector in the orthonormal basis of Krylov subspace.
\begin{flalign}
&\phi'_{1,l}=x'_0, 
(\phi'_{x})_{k,l}=1,
\zeta'_{1,l}=0 ,
(\zeta'_x)_{k,l}=0,&
&k=1,\dots,N  \; \text{and} \; l=1,\dots,M.&
\end{flalign}
%for $k=1,\dots,N$ and $l=1,\dots,M$.

\begin{figure}[H]
	\centering
		\includegraphics[width=4.3 cm]{JFNK.png}
	  \caption{Calculation flow chart of the banded preconditioner JFNK method}
	  \label{Fig: JFNK}
\end{figure}

\subsection{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., \citep{trefethen2022}.  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. 

\subsubsection{Building preconditioner matrix}
For a good preconditioner $\bold P$, it should be cheap to form and to factorizeis. Meanwhile, the preconditioned Jacobian $\boldsymbol{{\rm J}_t}\bold P^{-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 \citep{Knoll2004}. This paper applies the numerical scheme to the linearized governing equations which apply formally in the limit $\epsilon'\rightarrow0.$

The equations of the linear problem of computing the Havelock potential for flow past a submerged point source can be described, as follows \citep{Havelock1932, Lustri2013, Noblesse1978, Noblesse1981}:
\begin{flalign}
&{\zeta^\prime}_{x}={\phi^\prime}_{z}& &\text{on } z'={\zeta^\prime}(x',y')&\\
&{\phi^\prime}_{x}-1+\frac{\zeta^\prime}{{F'}^2}=0& &\text{on } z'={\zeta^\prime}(x',y')&
\end{flalign}

According to the linear free surface boundary condition, the boundary integral equation can be described, as follows: 
\begin{flalign}
&2{\pi}({\phi^\prime}(Q)-x')=-\frac{\epsilon^\prime}{(x'^2+y'^2+1)^\frac{1}{2}}
+\int_{0}^{\infty}\int_{-\infty}^{\infty}{{\phi^\prime}_{\rho}(P)K^{(3)}({\rho^\prime},{\sigma^\prime},x',y')}d\sigma^\prime d\rho^\prime&
\end{flalign}
where 
\begin{flalign*}
&K^{(3)}({\rho^\prime},{\sigma^\prime},x',y')=
\frac{1}{[({\rho^\prime}-x'^2)^2+({\sigma^\prime}+y')^2]^\frac{1}{2}}+\frac{1}{[({\rho^\prime}-x'^2)^2+({\sigma^\prime}-y')^2]^\frac{1}{2}}&
\end{flalign*}

After numerical discretization, the linear system can be described, as follows:
\begin{flalign}
&\boldsymbol{F1}_{k,l}={\phi}_{k,l}(Q)+\frac{{\zeta^\prime}_{k,l}(Q)}{{F'}^2}-1&
\end{flalign}

\begin{adjustwidth}{-\extralength}{0cm} 
\begin{flalign}
\boldsymbol{F2}_{k,l}=2{\pi}({\phi'}_{k,l}(Q)-x'_{k})+\frac{\epsilon^\prime}{[{x'^2}_{k}(Q)-{y'^2}_{l}(Q)+1]^\frac{1}{2}}-\sum_{i=1}^{N}\sum_{j=1}^{M}{w(i,j)[(\zeta^\prime_\rho)_{i,j}-{(\zeta^\prime_x)}_{i,j}]K^{(3)}_{i,j,k,l}-{(\zeta^\prime_x)}_{i,j}I}
\end{flalign}
\end{adjustwidth}
where ${w(i,j)}$ is the weighting function for numerical integration, for $k=1\dots(N-1), l=1\dots M.$ Then the linear Jacobian can be calculated directly, by differentiating the linear system with respect to ${\phi^\prime}_{1,m},{(\phi^\prime_x)}_{n,m},{\zeta^\prime}_{1,m}\text{ and }{(\zeta^\prime_x)}_{n,m}$. Therefore, the preconditioner matrix $\bold P$ can be formed cheaply, and the eigenvalues of $\boldsymbol{{\rm J}_t{\rm P}^{-1}}$ obviously cluster as shown in \autoref{Fig: 3}.

\begin{figure}[h] 
 \centering 
 \subfigure[$\boldsymbol{{\rm J}_t}$ ]{  \includegraphics[width=6 cm]{3.png}  } 
 \subfigure[$\boldsymbol{{\rm J}_t{\rm P}^{-1}}$ ]{  \includegraphics[width=6.7 cm]{4.png} } 
 \caption{A plot of the spectrum for $\boldsymbol{{\rm J}_t}$ and $\boldsymbol{{\rm J}_t{\rm P}^{-1}}$ on a $31\times11$ mesh.} 
\label{Fig: 3}
\end{figure}




\subsubsection{Preconditioner factorisation and storage}
The JFNK method requires the result of the product of the inverse preconditioner matrix and vector, $\boldsymbol{{\rm P}^{-1}v}$. In general, the operation of inverting a matrix should be converted to solving a system of linear equations, $\boldsymbol{{\rm P}r=v}$. Find the solution $\boldsymbol{r}$, the result of  $\boldsymbol{{\rm P}^{-1}v}$ will be got. In order to calculate this linear system conveniently, the preconditioner can be divided up into four submatrices and factorized using the block decomposition,
\begin{flalign}
&  \bold P=\begin{bmatrix}\bold A & \bold B \\ \bold C & \bold D \end{bmatrix}=\begin{bmatrix}\bold I & \bold 0 \\ \bold C \bold A^{-1} & \bold I \end{bmatrix}\begin{bmatrix}\bold A & \bold 0 \\ \bold 0 & \bold D-\bold C \bold A^{-1} \bold B \end{bmatrix}\begin{bmatrix}\bold I & \bold A^{-1} \bold B \\ \bold 0 & \bold I \end{bmatrix}&
\end{flalign}
where $\bold I$ is the unit matrix, $\bold A, \bold B, \bold C$ and $\bold D$ are the four submatrices which are constructed on the base of Jacobian $\boldsymbol{{\rm J}_t}$. 

Accordingly, the vector $\boldsymbol{v}$ can be divided into upper and lower parts  $[\boldsymbol{ v_1}\quad\boldsymbol{ v_2}]^{\rm T}$, and then the  solution $\boldsymbol{r}$ can be got after three cheap steps, as follows:
%It is possible to solve the system $\bold P \boldsymbol{ r}=\boldsymbol{v}$ by performing the following operations,
\begin{flalign}
\label{Eq 29}
&  \begin{bmatrix}\boldsymbol{o_1}  \\ \boldsymbol{o_2} \end{bmatrix}=\begin{bmatrix}\boldsymbol{v}_1  \\ \boldsymbol{v}_2-\bold C \bold A^{-1}\boldsymbol{v}_1 \end{bmatrix}&
\end{flalign}

\begin{flalign}
\label{Eq 30}
&  \begin{bmatrix}\boldsymbol{ s_1}  \\ \boldsymbol{ s_2} \end{bmatrix}=\begin{bmatrix}\bold A^{-1}\boldsymbol{ o_1}  \\ (\bold D-\bold C \bold A^{-1} \bold B)^{-1}\boldsymbol{ o_2} \end{bmatrix}&
\end{flalign}

\begin{flalign}
\label{Eq 31}
&  \begin{bmatrix}\boldsymbol{ r_1}  \\ \boldsymbol{ r_2} \end{bmatrix}=\begin{bmatrix}\boldsymbol{s_1} -\bold A^{-1} \bold B\boldsymbol{ s_2 } \\ \boldsymbol{s_2}\end{bmatrix}&
\end{flalign}

The calculation of $\boldsymbol{{\rm P}^{-1}v}$ in the Eq. \eqref{Eq 26} can be facilitated according to the progressive order from Eqs. \eqref{Eq 29} - \eqref{Eq 31}. The reasons are as follows: the submatrix $\bold A$ is tridiagonal, allowing for easy storage and fast factorization; the submatrices $\bold B$ and $\bold C$ are only used in matrix vector multiplication operations and thus can be implemented as functions that perform these operations rather than stored as matrices.






\subsubsection{The banded preconditioner}
After the factorisation operation, the calculation and storage of preconditioner matrix $\bold P$ are optimized. However, the size of submatrix $\bold D$ is $(N+1)M\times(N+1)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. 
%This paper takes the banded preconditioner matrix method and GPU acceleration technique to overcome these problems.

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 \autoref{Fig: bandedpreconditioner}. Moreover, batch construction avoids the problem of insufficient running memory due to the large size of the submatrix $\bold D$. The compressed sparse row (CSR) data format is used to save this matrix. Hence, a lot of memory can be saved.



 \begin{figure}[H]
	\centering
		\includegraphics[width=12cm]{bandedpre.png}
	  \caption{Construction of the banded preconditioner}
	  \label{Fig: bandedpreconditioner}
\end{figure}

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

\begin{figure}[h] 
\begin{adjustwidth}{-\extralength}{0cm}
\centering 
\subfigure[$band=1$ ]{  \includegraphics[width=5.4 cm]{5.png}  } 
\subfigure[$band=11$ ]{  \includegraphics[width=5.1cm]{6.png} } 
\subfigure[$band=21$ ]{  \includegraphics[width=5.1cm]{7.png} } 
\caption{A plot of the spectrum for $\boldsymbol{{\rm J}_t}$ and $\boldsymbol{{\rm J}_t{\rm P}^{-1}}$ on a $31\times11$ mesh for: $band=1,band=11,band=21$.}  
\label{verified}
\end{adjustwidth}
\end{figure}
 
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 \autoref{Fig: band}. The case is $F'=0.7$ and  $\epsilon'=0.4$, computed on a $121\times41$ mesh, when $b'$ (bandwidth $band=b'\times(N+1)$) is less than 14, an ill-conditioned coefficient matrix is formed, the accuracy of the solution is low. The runtime decreases with $b'$ ranges from 14 to 16, then the runtime increases monotonically with $b'$ ranges from 16 to 20. For the case of $121\times41$ mesh, the shortest running time is 5.6s with the optimal bandwidth $band=16\times(N+1)$. Therefore, provided that the appropriate bandwidth is selected, not only can save memory, but also can improve the computational efficiency.  
\begin{figure}[h]
	\centering
		\includegraphics[width=8cm]{band.png}
	  \caption{The plot of runtime against the bandwidth, computed on a $121\times41$ mesh with ${\Delta}x'=0.3,{\Delta}y'=0.3,band=b'\times(N+1)$.}
\label{Fig: band}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{GPU Parallel Computing Framework}\label{Section 4}


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 \citep{Chen2019}. This hardware framework makes GPU own plenty of threads naturally to process large amounts of data simultaneously. The computational efficiency of nonlinear ship wave can be greatly improved by  utilizing the GPU acceleration technique . 
%Therefore, this paper adopts the GPU acceleration technique to improve the computational efficiency of nonlinear ship waves simulation.
\subsection{Parallel computing framework design}
In this paper, Compute Unified Device Architecture (CUDA) language is used to develop the numerical scheme for computing ship wave patterns. CUDA is a parallel computing platform and programming model created by NVIDIA and implemented by GPU \citep{NVIDIA2022}. CUDA toolkit includes abundant GPU accelerated libraries, optimization tools and a runtime library, which can be compiled in C language, $\textup{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. Therefore, a typical CUDA program is consisted of two parts: a host part runs on CPU and a device part runs on GPU. The host code is responsible for the logic work, environment configuration, instructions to set up parallelism, and data communication between the host and the device. The main responsibility of device code is to process tasks in parallel, and the dimensions of the grid and thread are specified by configuration parameters in the kernel function. \autoref{Fig: 5} shows the CUDA execution mode and thread hierarchy.

\begin{figure}[h]
	\centering
		\includegraphics[width=8 cm]{8.png}
	  \caption{ Illustration of CUDA execution mode and thread organization hierarchy \citep{NVIDIA2022}}
\label{Fig: 5}
\end{figure}

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\cite{Sun2022} show that the parts of creating nonlinear system and inverting preconditioner matrix take up most of the time, as shown in \autoref{Fig: 6}. This figure shows the computation time distribution of the CPU solver on a $151\times51$ 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.
%In contrast, the parts of building preconditioner matrix and solving linear equations by Givens rotation take very little time relatively, this paper does not focus on parallel computing in these two parts. 

\begin{figure}[h]
\centering
\includegraphics[width=9 cm]{9.png}
  \caption{ The computation time distribution of ship wave solver. The alphabet \textup{I} represents the part of inverting preconditioner matrix, the alphabet \textup{C} represents the part of creating nonlinear system, the alphabet \textup{B} represents the part of building preconditioner matrix, the alphabet \textup{S} represents the part of solving linear equations by GMRES algorithm and the alphabet \textup{O} represents the part of other code in the solver.}
\label{Fig: 6}
\end{figure}


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 $\boldsymbol{{\rm P}^{-1}v}$ is calculated directly using the QR decomposition results, by combining the result of $\boldsymbol{{\rm P}^{-1}v}$ with the approximate solution $\boldsymbol{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 $\boldsymbol{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 \autoref{Fig: 7}, which shows the calculation procedure more clearly. 
%In the whole calculation, the exchange of data between CPU and GPU only represents a low percentage of the total runtime.

\begin{figure}[H]
	\centering
		\includegraphics[width=10 cm]{GPU.png}
	  \caption{The computational flow chart of GPU implementation.}
\label{Fig: 7}
\end{figure}



\subsection{GPU solver implementation}
\subsubsection{Creating nonlinear system }
The programming for creating nonlinear system on GPU by using CUDA language is briefly shown in \autoref{table2}. 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. \eqref{Eq 8}, Eq. \eqref{Eq 9} and Eqs. \eqref{Eq 17}-\eqref{Eq 20} simultaneously.

In the device part, there are two device functions which are called by the kernel function multiple times. According to Eq. \eqref{Eq 16}, these two device functions are formed to solve the singularity in the second integral of the boundary integral equation. Eq. \eqref{Eq 16} is distributed into 16 special threads 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 $(N-1)M$ nonlinear equations are built. Then arbitrarily choose 5 blocks to calculate the free surface condition and the radiation condition,  $(N+3)M$ nonlinear equations can be obtained. Therefore, $2(N+1)M$ nonlinear equations are formed by using the $1024\times(N-1)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.



\begin{table}[H]
  \centering
   \caption{The programming on GPU. }
     \label{table2}
   \begin{tabular}{l|l}
   \Xhline{1pt}
     \multicolumn{2}{l}{\textbf{Device Part}}  \\
     \Xhline{0.4pt}
    \multicolumn{2}{l}{\textit{Device function 1}} \\
     \Xhline{0.4pt}
     1   &  $\_\_$device$\_\_$double Integral(double s,double t,doubleA,double B,double C)\\
     2   &\{ \\
     3   &\quad double val=  t/sqrt(A)*log(2*A*s + B*t+ 2*sqrt(A*(A*s*s+ B*s*t +C*t*t))); \\
     4   &\quad return val; \\ 
     5   & \} \\
  \Xhline{0.4pt}
  \multicolumn{2}{l}{\textit{Device function 2}} \\
   \Xhline{0.4pt}
     1   & $\_\_$device$\_\_$double Integrall(double s,double t,doubleA,double B,double C)\\
     2   &\{ \\
     3   &\quad double val = s/sqrt(C)*log(2*C*t + B*s+ 2*sqrt(C*(A*s*s+ B*s*t +C*t*t))); \\
     4   &\quad return val; \\ 
     5   & \} \\
     \Xhline{0.4pt}
    \multicolumn{2}{l}{\textit{ Kernel function}}\\
     \Xhline{0.4pt}
     1   & $\_\_$global$\_\_$void nonlinear(const long int N, const long int M, $\dots$ )\\
     2   &\{ \\
     3   &\quad $\_\_$shared$\_\_$double A, B, C, $\dots$ ;\\
     4   &\quad $\_\_$shared$\_\_$long int K, l, blockpos, $\dots$ ;\\
     5   &\quad long int i, j, threadPos, $\dots$ ;\\
     6   &\quad threadPos = threadIdx.x;  k = blockIdx.x;  l = blockIdx.y ;\\
     7   &\quad // Calculate necessary values\\
     8   &\quad while(threadPos<(M*N)) \{ $\dots$ \}\\
     9   &\quad // Calculation of the 16 parts to the closed integral \\
    10  &\quad  if(threadIdx.x==blockDim.x-1) \{ $\dots$ \}\\
    11  &\quad  $\dots$ \\
    12  &\quad  if(threadIdx.x==blockDim.x-16) \{ $\dots$ \}\\
    13  &\quad // Sum up all thread contributions\\
    14  &\quad for(i=blocDim/2;i>0;i=i/2)\{ $\dots$ \}\\
    15  &\quad //Split the free surface condition and radiation conditions between 5 blocks \\
    16  &\quad  if(k==0\&\&l==0) \{ $\dots$ \} $\dots$ if(k==0\&\&l==5) \{ $\dots$ \}\\
    17  & \} \\
    \Xhline{1pt}
    \multicolumn{2}{l}{ \textbf{Host Part}}  \\
     \Xhline{0.4pt}
     1   & int main \\
     2   &\{ \\
     3   &\quad  // Initialize data on CPU \\
     4   &\quad InitialData(double* cpuData, double* gpuData);\\
     5   &\quad long int i, j, threadPos, $\dots$ ;\\
     6   &\quad // Copy data from CPU to GPU\\
     7   &\quad cudaMemcpy(cpuData, gpuData, Datasize, cudaMemcpyHostToDevice);\\
     8   &\quad // Set dimensions of the grid and thread \\
     9   &\quad dim3 block(1024, 1); dim3 grid(M, N - 1); \\
    10  &\quad // Start Kernel function\\
    11  &\quad Nonlinear$\ll$ grid, block$\ll$(gpuData); \\
    12  &\quad // Copy data from GPU to CPU\\
    13  &\quad cudaMemcpy(gpuData, cpuData, Datasize, cudaMemcpyDeviceToHost);\\
    14  &\quad return 0;\\
    15  & \} \\
     \Xhline{1pt}
   \end{tabular}
\end{table}

\subsubsection{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 $2(N+1)M\times2(N+1)M$, and four submatrices are formed by block decomposition method and banded preconditioner method: Submatrix $\bold A$ is a tridiagonal matrix of size $3\times(N+1)M$; the submatrix $\bold B$ and $\bold C$ only differ between coefficients, and the base matrix $\boldsymbol{{\rm B}_0}$ can be constructed to represent them respectively with a size of $(N+1)\times(N+1)$, $\bold B=1/F'^2 \cdot\boldsymbol{{\rm B}_0}$ and $\bold C=2\pi \cdot \boldsymbol{{\rm B}_0}$; the submatrix ${\bold D}$ is a sparse matrix with a size of $band\times(N+1)M/3$. Firstly, $M\times(N+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 $N+1$ thread blocks for parallel computation of $\bold{ CA^{-1}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 $M\times(N+1)$ thread blocks are called to perform parallel operation of $\bold{D-CA^{-1}B}$. The programming for building preconditioner matrix on GPU by using CUDA language is briefly shown in \autoref{PreconditionerMatrix}.

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.

\begin{table}[H]
  \centering
   \caption{The programming on GPU. }
     \label{PreconditionerMatrix}
   \begin{tabular}{l|l}
   \Xhline{1pt}
     \multicolumn{2}{l}{\textbf{Device Part}}  \\
    \Xhline{0.4pt}
    \multicolumn{2}{l}{\textit{ Kernel functions}}\\
     \Xhline{0.4pt}
     1   & $\_\_$global$\_\_$void precondition(const long int N, const long int M, $\dots$ )\\
     2   &\{ \\
     3   &\quad // Calculate submatrix $\bold{A}$ \\
     4   &\quad while(threadPos<N+3) \{ $\dots$ \};\\
     5   &\quad if(blockIdx.x==0) \{ $\dots$ \};\\
     6   &\quad // Calculate basis matrix $\boldsymbol{{\rm B}_0}$\\
     7   &\quad if(blockIdx.x==0) \{ while(threadPos<N+1) \{ $\dots$ \} \};\\
     8   &\quad // Calculate submatrix $\bold{D}$ \\
     9   &\quad while(threadPos<N*M) \{ $\dots$ \}; \\
    10  &\quad  // Sum up all thread contributions\\
    11  &\quad  for(i=blocDim/2;i>0;i=i/2)\{ $\dots$ \} \\
    12  &\} \\
    13  & $\_\_$global$\_\_$void matrix(const long int N, const long int M, $\dots$ )\\
    14  &\{ \\
    15  & \quad // Calculate $\bold{CA^{-1}B}$ \\
    16  & \quad double Bcoef = 1/($F'\times F'$); \\
    17  & \quad double Ccoef = 2*3.1415;\\
    18  & \quad if(blockIdx.x==0) \{ while(threadPos<N+1) \{ $\dots$ \} \};\\
    19  & \} \\
    20  & \ $\_\_$global$\_\_$void Schur(const long int N, const long int M,double*\_b,double*D )\\
    21  &\{ \\
    22  & \quad // Calculate $\bold{D-CA^{-1}B}$ \\
    23  & \quad int l=blockIdx.x; \\
    24  & \quad int k=blockIdx.y; \\
    25  & \quad int pos=1*(N+1)+k; \\
    26  & \quad while(threadPos<N+1)\\
    27  & \quad \{ \\
    27  & \quad \quad D[(1*(N+1)+threadPos)*M*(N+1)+pos]=D[(1*(N+1)+ \\
    28  & \quad \quad threadPos)*M*(N+1)+pos]-\_b[threadPos*(N+1)+k]; \\
    29  & \quad \quad threadPos+=blockDim.x; \\
    30  & \quad  \} \\
    31  & \} \\
     \Xhline{1pt}
   \end{tabular}
\end{table}



\subsubsection{Inverting preconditioner matrix}
%The process of inverting preconditioner matrix is shown in Section \ref{Section 3}, which is transformed to solve a system of sparse linear equations. In the CUDA library, the cuSolverSP library was mainly designed for sparse linear systems. 
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 \citep{Grossman2017}. 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 $\boldsymbol{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   $\bold P^{-1}\boldsymbol{v}$ \citep{trefethen2022}. 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 $\bold P^{-1}\boldsymbol{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 \autoref{table1}.
\begin{table*}
\caption{The list of CUDA functions for QR factorization}
\label{table1}
\begin{tabular}{lll}
 \Xhline{1pt}
     \textbf{No.} & \textbf{Function name} & \textbf{Goal}\\
   \Xhline{0.4pt}
     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\\
    \Xhline{1pt}
\end{tabular}
\end{table*}


\subsubsection{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.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Simulations and Discussion}\label{Section 5}
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, the capability of GPU solver is verified by comparing  simulation results with real ship wakes. 

Both CPU and GPU solvers are executed on a high-performance computing cluster, the CPU Clock Speed is 1.90 GHz, the GPU card is NVIDIA Tesla A100, which includes 6912 CUDA cores and 40 GB graphics memory. The parameters of the computing environment are listed in \autoref{table3}. 

\begin{table*}
   \caption{The computing environment of high-performance computing cluster}
   \label{table3}
   \begin{tabular}{lll}
   \Xhline{1pt}
                        & \textbf{CPU}                 & \textbf{GPU}\\
   \Xhline{0.4pt}
     Card           &Intel Xeon Bronze 3204 &NVIDIA Tesla A100\\
     Memory     & 64 GB                             & 40 GB\\
     Max Cores & 6 per node                      &6912\\
     Programming language& C++             & CUDA, C++\\
    \Xhline{1pt}
   \end{tabular}
\end{table*}



\subsection{Verification of the banded preconditioner JFNK method}
To reveal the effectiveness of the banded preconditioner method, numerical simulations on different mesh sizes,i.e.,$181\times61$, $241\times81$, $301\times101$ and $361\times121$ with  $\Delta x'=0.3,\Delta y'=0.3$ are carried out.

The overall runtimes against bandwidth on these four mesh sizes  are illustrated in  \autoref{Fig: bs}. From this figure, it can be  seen that   a  optimal value of bandwidth $b'$ exists for a certain mesh size. Furthermore, the optimal value of $b'$ increases with the mesh size and approximately equals $\frac{M}{3}$, in which $M$ means the number of latitude lines of the mesh. Therefore, the optimal bandwidth can be set to $\frac{M}{3}\times(N+1)$ to get an optimal efficiency, which is called the optimal bandwidth selection rule in this paper.

According to the optimal bandwidth selection rule, the running memorys against the bandwidth are shown in \autoref{table5}. 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.

 
\begin{figure}[h] 
 \centering 
   \includegraphics[width=8.5cm]{bands.png}
 \caption{Optimal values of bandwidth $b'$ for different mesh sizes.} 
\label{Fig: bs}
\end{figure}


\begin{table*}
   \caption{The running memory usage before and after applying the banded preconditioner method}
   \label{table5}
   \begin{tabular}{ccccc}
   \Xhline{1pt}
      \textbf{Mesh size} &\textbf{Before} &\textbf{$b'$} & \textbf{After} & \textbf{Reduction ratio}\\
   \Xhline{0.4pt}
     $181\times61$    & 0.91GB  & 19  & 0.28GB &3.2 \\
     $241\times81$    & 2.9GB    & 24  & 0.88GB &3.3 \\
     $301\times101$  & 6.9GB    & 33  & 2.3GB   &3.0 \\
     $361\times121$  & 15GB     & 38  & 4.6GB   &3.3 \\
    \Xhline{1pt}
   \end{tabular}
\end{table*}



\subsection{Verification of the GPU solver}
\subsubsection{Accuracy}
To verify the accuracy of the GPU solver, numerical simulations are conducted on $F'=0.7$ and  $\epsilon'=0.4$ with a $361\times121$ mesh and $\Delta x'=0.3,\Delta y'=0.3$. The simulated wave heights on the centerline are compared with that of the  CPU solver proposed by Sun \cite{Sun2022}, which is shown in \autoref{Fig 8}.  From this figure, it can be seen that 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 Mean Square Error (MSE) is used to further explain the error between them, as follows:
\begin{flalign}
\label{MSE} 
&MSE=\frac{1}{n}\sum_{i=1}^{n}(Truch_i-Value_i)^2&
\end{flalign}
where $n$ is the amount of data, $Truch_i$ and $Value_i$ represent CPU results and GPU results respectively.
According to  Eqs.\eqref{MSE}, 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 \cite{Sun2022}, the accuracy of the proposed GPU solver can also be acceptable.


\begin{figure}[h]
\begin{adjustwidth}{-\extralength}{0cm}
	\centering
		\includegraphics[width=16cm]{accuracy.png}
	  \caption{A comparison of the centerline profiles for the simulation results of CPU solver and GPU solver, computed on a $361\times121$ mesh with $\Delta x'=0.3,\Delta y'=0.3, F'=0.7$ and  $\epsilon'=0.4$. The solid line represents the simulation result of the GPU solver, the solid circles represent the simulation result of the CPU solver.}
   \label{Fig 8}
\end{adjustwidth}
\end{figure}  






\subsubsection{Efficiency}
To verify the efficiency of the GPU solver, numerical simulation are conducted on $F'=0.7$ and  $\epsilon'=0.4$ with five mesh sizes, namely  $121\times41$, $181\times61$, $241\times81$, $301\times101$,  $361\times121$ and $\Delta x'=0.3,\Delta y'=0.3$. The overall runtimes of GPU solver are compared with that of the CPU solver proposed by Sun \cite{Sun2022}, as shown in  \autoref{Fig: 12}. 
From this figure, it can be seen clearly 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  \autoref{table6}. From this table,  it can be seen  that 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.
 

\begin{figure*}[h]
\begin{adjustwidth}{-\extralength}{0cm}
   	\centering
		\includegraphics[width=18cm]{runtime.png}
	  \caption{The comparison of runtime between 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.}    
   \label{Fig: 12}
\end{adjustwidth}
\end{figure*}


\begin{table*}
   \caption{The comparisons of runtime between CPU solver and GPU solvers on different mesh sizes, the CPU solver is proposed by Sun \cite{Sun2022}, one GPU solver is proposed in this paper, the other GPU solver is proposed by Pethiyagoda \cite{Pethiyagoda2016}(results of Exp. ). }
   \label{table6}
    \begin{tabular}{ccccc}
   \Xhline{1pt}
      \textbf{Mesh size} &\textbf{CPU solver} &\textbf{Exp.}&\textbf{GPU solver} & \textbf{Accelerated-up ratio}\\
   \Xhline{0.4pt}
     $121\times41$    & 8.96E+1 s  &1.61E+1 s    & 5.60E+0 s   &16.1\\
     $181\times61$    & 2.70E+2 s  & 1.22E+2 s   & 1.13E+1 s   &23.9\\
     $241\times81$    & 8.70E+2 s  & 5.51E+2 s   & 5.37E+1 s   &16.0\\
     $301\times101$  & 2.56E+3 s  & 1.78E+3 s   & 1.32E+2 s   &19.3\\
     $361\times121$  & 6.08E+3 s  & 5.04E+3 s   & 2.95E+2 s   &20.6\\
    \Xhline{1pt}
   \end{tabular}
\end{table*}


The proposed GPU solver has also been compared with another GPU solver proposed by Pethiyagoda \cite{Pethiyagoda2016} on these cases. The comparison of computation time between them is shown in \autoref{table6}. Obviously, the efficiency of the GPU solver proposed in this paper is higher than the GPU solver proposed by Pethiyagod \cite{Pethiyagoda2016}, and the advantage is more significant with the increase of the mesh size. The reason is that Pethiyagoda \cite{Pethiyagoda2016} 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.


\subsubsection{Capability}
In the proposed GPU solver, like the CPU solver \cite{Sun2022,Cai2022}, three parameters can be used to regulate simulation results, namely Froude number, source strength and source type. Numerical simulations show that the wake angle is tending to decrease with Froude number but increase with the source strength, and that the divergent waves will be prominent with the sufficiently large Froude number or the Rankine source, and that the wavelength dramatically increases with Froude number.  By adjusting these parameters, the wake characteristics including the wake angle, wave amplitude, prominent waves and wavelength can be changed, resulting in the polymorphic ship wave patterns. For example, the wake waves of high-speed boats can be generated by using large Froude number, low source strength and Rankine source. As the vessel speed decreases and overall length increases, it is appropriate to choose small Froude number and a higher strength Kelvin source.  Comparisons between the three types of real ship waves and the corresponding GPU solver simulation results are shown in \autoref{Fig: 10}.  It is found that the simulation results are in good agreement with the real ship waves, the proposed GPU solver can also generate high quality simulation patterns for 3-D nonlinear ship wave.

\begin{figure}[H] 
\begin{adjustwidth}{-\extralength}{0cm}
  \centering 
 \subfigure[Real speedboat wake]{  \includegraphics[scale=0.6]{boat1.png}  } 
 \subfigure[Simulation result of  speedboat wake]{  \includegraphics[scale=0.6]{boatresult1.png}  } 
  \subfigure[Real fishing ship wake]{  \includegraphics[scale=0.6]{fish.png}  } 
 \subfigure[Simulation result of  fishing ship wake]{  \includegraphics[scale=0.6]{fishresult.png}  } 
 \subfigure[Real large vessel wake]{  \includegraphics[scale=0.6]{vessel.png} } 
 \subfigure[Simulation result of  large vessel wake]{  \includegraphics[scale=0.6]{vessresult.png} }   
  \caption{Comparisons between the three types of real ship waves and the corresponding the GPU solver simulation 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.} 
  \label{Fig: 10}
\end{adjustwidth}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}\label{Section 6}


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,
the following conclusions can be drawn.
\begin{enumerate}
\item[(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{M}{3}\times(N+1)$, with more than 66\% GPU memory can be saved.
\item[(2)] The GPU solver can obtain an accurate numerical solution. The mean square error of GPU solver results and CPU solver results is $MSE$=9.37E-8, which is acceptable. 
\item[(3)] By designing the GPU parallel computing framework, the computation of ship wave simulation is accelerated up to 20 times.
\end{enumerate}

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. 
%Moreover, for the optimal bandwidth selection rule, it only draws conclusion from plenty of cases without mathematical inference, and there are deviations in few cases. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \section{Patents}

% This section is not mandatory, but may be added if there are patents resulting from the work reported in this manuscript.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{6pt} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% optional
%\supplementary{The following supporting information can be downloaded at:  \linksupplementary{s1}, Figure S1: title; Table S1: title; Video S1: title.}

% Only for journal Methods and Protocols:
% If you wish to submit a video article, please do so with any other supplementary material.
% \supplementary{The following supporting information can be downloaded at: \linksupplementary{s1}, Figure S1: title; Table S1: title; Video S1: title. A supporting video article is available at doi: link.}

% Only for journal Hardware:
% If you wish to submit a video article, please do so with any other supplementary material.
% \supplementary{The following supporting information can be downloaded at: \linksupplementary{s1}, Figure S1: title; Table S1: title; Video S1: title.\vspace{6pt}\\
%\begin{tabularx}{\textwidth}{lll}
%\toprule
%\textbf{Name} & \textbf{Type} & \textbf{Description} \\
%\midrule
%S1 & Python script (.py) & Script of python source code used in XX \\
%S2 & Text (.txt) & Script of modelling code used to make Figure X \\
%S3 & Text (.txt) & Raw data from experiment X \\
%S4 & Video (.mp4) & Video demonstrating the hardware in use \\
%... & ... & ... \\
%\bottomrule
%\end{tabularx}
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\authorcontributions{Conceptualization,  X.S. and M.C.; methodology, X.S.and M.C.; software, M.C.; validation, X.S., M.C. and J.D.; investigation, X.S.; resources, X.S.; data curation,J.D.; writing---original draft preparation, X.S. and M.C.; writing---review and editing,X.S. and M.C.; visualization, M.C; supervision, X.S. All authors have read and agreed to the published version of the manuscript.}

\funding{The work was supported by the National Key R\&D Program of China (No.2022YFB4300803, 2022YFB4301402), the Ministry of Industry and Information Technology of the People’s Republic of China(No. CBG3N21-3-3), and the National Science Foundation of Liaoning Province, China(No.2022-MS-159).The authors would like to express sincere thanks for their support.}

%\institutionalreview{In this section, you should add the Institutional Review Board Statement and approval number, if relevant to your study. You might choose to exclude this statement if the study did not require ethical approval. Please note that the Editorial Office might ask you for further information. Please add “The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board (or Ethics Committee) of NAME OF INSTITUTE (protocol code XXX and date of approval).” for studies involving humans. OR “The animal study protocol was approved by the Institutional Review Board (or Ethics Committee) of NAME OF INSTITUTE (protocol code XXX and date of approval).” for studies involving animals. OR “Ethical review and approval were waived for this study due to REASON (please provide a detailed justification).” OR “Not applicable” for studies not involving humans or animals.}

%\informedconsent{Any research article describing a study involving humans should contain this statement. Please add ``Informed consent was obtained from all subjects involved in the study.'' OR ``Patient consent was waived due to REASON (please provide a detailed justification).'' OR ``Not applicable'' for studies not involving humans. You might also choose to exclude this statement if the study did not involve humans.

%Written informed consent for publication must be obtained from participating patients who can be identified (including by the patients themselves). Please state ``Written informed consent has been obtained from the patient(s) to publish this paper'' if applicable.}

%\dataavailability{We encourage all authors of articles published in MDPI journals to share their research data. In this section, please provide details regarding where data supporting reported results can be found, including links to publicly archived datasets analyzed or generated during the study. Where no new data were created, or where data is unavailable due to privacy or ethical restrictions, a statement is still required. Suggested Data Availability Statements are available in section ``MDPI Research Data Policies'' at \url{https://www.mdpi.com/ethics}.} 

% Only for journal Nursing Reports
%\publicinvolvement{Please describe how the public (patients, consumers, carers) were involved in the research. Consider reporting against the GRIPP2 (Guidance for Reporting Involvement of Patients and the Public) checklist. If the public were not involved in any aspect of the research add: ``No public involvement in any aspect of this research''.}

% Only for journal Nursing Reports
%\guidelinesstandards{Please add a statement indicating which reporting guideline was used when drafting the report. For example, ``This manuscript was drafted against the XXX (the full name of reporting guidelines and citation) for XXX (type of research) research''. A complete list of reporting guidelines can be accessed via the equator network: \url{https://www.equator-network.org/}.}

%\acknowledgments{In this section you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).}

\conflictsofinterest{The authors declare no conflict of interest.} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Optional
%\sampleavailability{Samples of the compounds ... are available from the authors.}

%% Only for journal Encyclopedia
%\entrylink{The Link to this entry published on the encyclopedia platform.}

\abbreviations{Abbreviations}{
The following abbreviations are used in this manuscript:\\

\noindent 
\begin{tabular}{@{}ll}
JFNK & Jacobian-free Newton-Krylov \\
GMRES &  Generalized Minimum Residual \\
CUDA & Compute Unified Device Architecture\\
GPU & Graphics Process Unit\\
CPU & Central Processing Unit\\
MSE & Mean Square Error
\end{tabular}
}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{adjustwidth}{-\extralength}{0cm}
%\printendnotes[custom] % Un-comment to print a list of endnotes

\reftitle{References}

%\bibliographystyle{model5-names}
%\bibliographystyle{plainnat}

\bibliography{refs}

% For the MDPI journals use author-date citation, please follow the formatting guidelines on http://www.mdpi.com/authors/references
% To cite two works by the same author: \citeauthor{ref-journal-1a} (\citeyear{ref-journal-1a}, \citeyear{ref-journal-1b}). This produces: Whittaker (1967, 1975)
% To cite two works by the same author with specific pages: \citeauthor{ref-journal-3a} (\citeyear{ref-journal-3a}, p. 328; \citeyear{ref-journal-3b}, p.475). This produces: Wong (1999, p. 328; 2000, p. 475)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% for journal Sci
%\reviewreports{\\
%Reviewer 1 comments and authors’ response\\
%Reviewer 2 comments and authors’ response\\
%Reviewer 3 comments and authors’ response
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\PublishersNote{}
\end{adjustwidth}
\end{document}

