iccsa-20-waves

main.tex (27724B)

---

 \documentclass[runningheads]{llncs}
 
 \usepackage{amsmath}
 \usepackage{amssymb}
 \usepackage{booktabs}
 \usepackage{cite}
 \usepackage{graphicx}
 \usepackage{url}
 
 \newcommand{\Jacobian}{\mathbb{J}}
 \newcommand{\Real}[1]{\operatorname{Re}#1}
 \newcommand{\Imag}[1]{\operatorname{Im}#1}
 \newcommand{\VectorL}[1]{\left[\begin{array}{l}#1\end{array}\right]}
 \newcommand{\VectorR}[1]{\left[\begin{array}{r}#1\end{array}\right]}
 \newcommand{\DerivativeT}[1]{\frac{\partial{#1}}{\partial{}t}}
 
 \begin{document}
 
 \title{Virtual testbed: Simulation of ocean wave reflection from the ship
     hull\thanks{Supported by Saint Petersburg State University (grants
     no.~51129371 and~51129725) and Council for grants of the President of the Russian
     Federation (grant no.~\mbox{MK-383.2020.9}).}}
 \author{%
     Ivan Petriakov\orcidID{0000-0001-5835-9313}\and\\
     Alexander Degtyarev\orcidID{0000-0003-0967-2949} \and\\
     Denis Egorov\orcidID{0000-0001-5478-2798}\and\\
     Ivan Gankevich\textsuperscript{*}\orcidID{0000-0001-7067-6928} \and\\
     Anton Gavrikov\orcidID{0000-0003-2128-8368} \and\\
     Artemii Grigorev\orcidID{0000-0003-0501-7656} \and\\
     Vasily Khramushin\orcidID{0000-0002-3357-169X}
 }
 
 \titlerunning{Simulation of ocean wave reflection from the ship hull}
 \authorrunning{I.\,Petriakov et al.}
 
 \institute{Saint Petersburg State University\\
     7-9 Universitetskaya Emb., St Petersburg 199034, Russia\\
     \email{st049350@student.spbu.ru},\\
     \email{a.degtyarev@spbu.ru},\\
     \email{st047824@student.spbu.ru},\\
     \email{i.gankevich@spbu.ru},\\
     \email{st047437@student.spbu.ru},\\
     \email{st016177@student.spbu.ru},\\
     \email{v.khramushin@spbu.ru}\\
     \url{https://spbu.ru/}}
 
 \maketitle
 
 \begin{abstract}
 
     Diffraction and radiation forces result from the interaction between the
     ship hull and the moving fluid. These forces are typically simulated using
     added masses, a method that uses mass to compensate for not computing these
     forces directly. In this paper we propose simple mathematical model to
     compute diffraction force. The model is based on Lagrangian description of
     the flow and uses law of reflection to include diffraction term in the
     solution. The solution satisfies continuity equation and equation of
     motion, but is restricted to the boundary of the ship hull. The solution
     was implemented in velocity potential solver in Virtual testbed~--- a
     programme for workstations that simulates ship motions in extreme
     conditions. Performance benchmarks of the solver showed that it is
     particularly efficient on graphical accelerators.
 
 \keywords{%
 ocean wave diffraction
 \and ocean wave radiation
 \and fluid velocity field
 \and law of reflection
 \and OpenCL
 \and OpenMP
 \and GPGPU.
 }
 \end{abstract}
 
 \section{Introduction}
 
 There are two mathematical models that describe rigid body motion and fluid
 particle motion: equations of translational and angular motion of the rigid
 body (Newton's Second Law) and Gerstner equations for ocean waves (which are
 solutions to linearised equations of motions for fluid particles). Usually, we
 use these models independently to generate incident ocean waves and then
 compute body motions caused by these waves. To measure the effect of still
 fluid on an oscillating rigid body (radiation forces) and the effect of fluid
 particles hitting the body (diffraction forces), we use added masses and
 damping coefficients~--- simplified formulae derived for small-amplitude
 oscillatory motion.
 
 But, what if we want to simulate large-amplitude rigid body motion with greater
 accuracy? There are two possible ways. First, we may use numerical methods such
 as Reynolds-Averaged Navier Stokes (RANS) method~\cite{wackers2011rans}.  This
 method is accurate, can be used for viscous fluid, but not the most
 computationally efficient. Second, we may solve Gerstner equations with
 appropriate boundary condition and use the solution to compute both rigid body
 and fluid particle motion around it. This paper explores this second option.
 Similar approach was followed by Fenton
 in~\cite{fenton1978vertical,fenton1993shoaling,fenton1994spectral}, but the
 distinctive feature of our approach is the use of law of reflection to derive
 analytic expressions for reflected waves and fluid particles.
 
 %and compares the results to the results of the method of added masses and numerical solver.
 
 
 
 \section{Methods}
 
 \subsection{Equations of motions with a moving surface boundary}
 
 An oscillating rigid body that floats in the water and experiences incident
 waves both reflects existing waves and generates new waves:
 \begin{itemize}
 \item fluid particles hit the body, causing it to move, and then reflect from it;
 \item moving body hits fluid particles and makes them move.
 \end{itemize}
 Both wave reflection and generation have the same nature~--- they are
 caused by the collision of the particles and the body~--- hence we
 describe them by the same set of formulae. Hereinafter we borrow the
 mathematical notation for Lagrangian description of the flow
 from~\cite{nouguier2015}.
 
 In Lagrangian description of flow instantaneous particle coordinates
 \(\vec{R}=(x,y,z)\) depend on particle positions at rest (independent initial
 coordinates) \(\vec{\zeta}=(\alpha,\beta,\delta)\) and time \(t\),
 i.e.~\(\vec{R}=\vec{R}(\alpha,\beta,\delta,t)\).  Using this notation
 equation of motion (conservation of momentum) is written as
 \begin{equation}
 \vec{R}_{tt} + g \hat{z} + \frac{1}{\rho} \nabla_{\vec{R}} p = 0,
 \label{eq-momentum}
 \end{equation}
 where \(\vec{R}\)~--- particle coordinates, \(\hat{z}\)~---
 unit vector in the direction of positive \(z\), \(p\)~--- pressure,
 \(\rho\)~--- fluid density, \(g\)~--- gravitational acceleration.
 Continuity equation (conservation of mass) is written as
 \begin{equation*}
 \left|\Jacobian\right| = 1, \qquad
 \frac{\partial}{\partial t}\left|\Jacobian\right| = 0, \qquad
 \Jacobian = \left[
 \begin{array}{lll}
 x_\alpha & y_\alpha & z_\alpha \\
 x_\beta & y_\beta & z_\beta \\
 x_\delta & y_\delta & z_\delta \\
 \end{array}
 \right].
 \end{equation*}
 Multiplying both sides of \eqref{eq-momentum} by \(\Jacobian\) and noting that
 \(\nabla_{\vec{\zeta}}p =\Jacobian\nabla_{\vec{R}}p\) gives
 \begin{equation*}
 \Jacobian\vec{R}_{tt} +
 g\vec\nabla\left( \vec{R} \cdot \hat{z} \right) +
 \frac{1}{\rho}\vec\nabla p
 \end{equation*}
 Following \cite{nouguier2015} we seek solution to this equation in the form
 of a simultaneous perturbation expansion for position, pressure, and the
 vorticity function:
 \begin{equation*}
 \begin{aligned}
 \vec{R} &= \vec{R}_0 + \vec{R}_1 + \vec{R}_2 + ... \\
 p &= p_a - \rho g \delta + p_1 + p_2 + ...
 \end{aligned}
 \end{equation*}
 Zeroth order terms are related to particles positions at rest:
 \begin{equation*}
 \begin{aligned}
 \vec{R}_0 &= \vec{\zeta} \\
 p_0 &= p_a - \rho g \delta
 \end{aligned}
 \end{equation*}
 First-order terms are solutions to linearised equation of motion and equation of
 continuity:
 \begin{equation*}
 \begin{aligned}
 & \vec{R}_{1tt} + g\vec\nabla\left( \vec{R}_1 \cdot \hat{z} \right) +
 \frac{1}{\rho}\vec\nabla p_1 = 0 \\
 & \vec\nabla \cdot \vec{R}_1 = 0
 \end{aligned}
 \end{equation*}
 We seek solutions of the form \(\vec{R}_1=\nabla w\) to make the flow
 irrotational. Plugging this form into the equations gives
 \begin{equation}
 \label{eq-mass-momentum}
 \begin{aligned}
 & \vec\nabla\left( w_{tt} + g w_\delta + p_1/\rho \right) = 0 \\
 & \Delta w = 0
 \end{aligned}
 \end{equation}
 The first equation denotes conservation of momentum (Newton's second law,
 equation of motion) and the second equation denotes conservation of mass
 (equation of continuity).
 
 When we have no boundary condition we seek solutions of the form
 \begin{equation}
 \label{eq-inverse-fourier}
 w\left(\alpha,\beta,\delta\right) =
 \Real{
 f(u,v)\exp\left(iu\alpha+iv\beta+k\delta-i\omega{}t\right)
 },
 \end{equation}
 where \(u\) and \(v\) are wave numbers.  We plug \eqref{eq-inverse-fourier} into
 continuity equation \eqref{eq-mass-momentum} where \(p_1\) is constant and get
 \(k=\sqrt{u^2+v^2}\).  That means that expression~\eqref{eq-inverse-fourier}
 is the solution to this equation when \(k\) is wave vector magnitude,
 i.e.~\(w\) decays exponentially with increasing water depth multiplied by
 wave vector magnitude.
 
 Then we plug \eqref{eq-inverse-fourier} into equation of motion
 \eqref{eq-mass-momentum} and get \(\omega^2=gk\), which is dispersion relation
 from classic linear wave theory.  That means that the expression is the solution
 to this equation when angular frequency depends on the wave number,
 i.e.~waves of different lengths have different phase velocities.
 
 Before solving this system of equations for an arbitrary moving surface
 boundary, we consider particular cases to substantiate the choice of the form of
 the solution.
 
 We use \(\Real{\exp\left(iu\alpha+iv\beta+k\delta-i\omega{}t\right)}\) to
 describe fluid particle potential. Here \(u\) and \(v\) are wave numbers,
 \(\omega\) is angular frequency, and \(k\) is wave vector.  This notation makes
 formulae short and is equivalent to the description that uses traditional
 harmonic functions.  This notation allows for easy transition to irregular
 waves via Fourier transforms which are essential for fast computations.  Such
 solutions will be studied in future work.
 
 \subsection{Stationary surface boundary}
 
 In this section we explore solutions stationary surface boundary in a form of
 infinite plane surface.  On such a boundary the projection of particle velocity
 to the surface normal is nought.  We write boundary conditions and
 corresponding solutions for different orientations of this boundary and then
 generalise these solutions to a parametric surface.
 
 \subsubsection{Infinite wall}
 
 On a vertical surface the boundary condition is written as
 \begin{equation*}
 \frac{d}{d t} \vec\nabla w \cdot \vec{n} =
 \frac{d}{d t} \frac{\partial}{\partial\alpha} w = 0;
 \qquad
 \alpha = \alpha_0;
 \qquad
 \vec{n} = \VectorR{1\\0\\0}.
 \end{equation*}
 Here we consider only \(\alpha\) coordinate, the derivations for \(\beta\) are
 similar. The potential of incident fluid particle has the form
 \begin{equation*}
 w\left(\alpha,\beta,\delta,t\right) =
 \exp\left(k\delta - i\omega t\right)
 \exp\left(iu\alpha + iv\beta\right).
 \end{equation*}
 Velocity vector of this particle is
 \begin{equation*}
 \begin{aligned}
 &
 \frac{d}{dt} \vec\nabla w =
 i\omega \left( \vec{d}_k+i\vec{d}_i \right)
 \exp\left(k\delta - i\omega t\right)
 \exp\left(iu\alpha + iv\beta\right);
 \\
 &
 \vec{d}_k = \VectorR{0\\0\\k};
 \qquad
 \vec{d}_i = \VectorR{u\\v\\0}.
 \end{aligned}
 \end{equation*}
 where \(\vec{d}_i\) is horizontal incident wave direction and \(\vec{d}_k\) is
 a vector that contains amplitude damping coefficient. (We use the vector
 instead of the scalar to shorten mathematical notation, otherwise we would have write a
 separate formula for vertical coordinate.)
 %In addition to this, it makes the solution for finite depth look the same
 %way as all other solutions.
 The law of reflection
 states that the angle of incidence equals the angle of reflection. Then
 the direction of reflected wave is%
 \footnote{Initially, we included the third
 component of incident wave direction making the
 vector complex-valued, however, the solution blew up as a result of mixing real
 and imaginary parts in dot products involving complex-valued vectors. The
 problem was solved by reflecting in two dimensions which is intuitive for
 ocean waves, but not for particles.}
 \begin{equation*}
 \vec{d}_r =
 \vec{d}_i-\vec{d}_s =
 \vec{d}_i-2\vec{n}\left(\vec{d}_i\cdot\vec{n}\right) =
 \VectorR{-u\\v\\0}.
 \end{equation*}
 We seek solutions of the form
 \begin{equation}
 \label{eq-w-alpha-0}
 w\left(\alpha,\beta,\delta,t\right) =
 \left[ C_1 \exp\left(iu\alpha\right) + C_2 \exp\left(-iu\alpha\right) \right]
 \exp\left(k\delta - i\omega t\right) \exp\left(iv\beta\right).
 \end{equation}
 We plug this expression into the boundary condition and get
 \begin{equation*}
 C_1\exp\left(i u\alpha_0 \right) -
 C_2\exp\left(-i u\alpha_0 \right) = 0,
 \end{equation*}
 hence \(C_1=C\exp(-iu\alpha_0)\) and \(C_2=-C\exp(iu\alpha_0)\).
 Constant \(C\) may take arbitrary values, here we set it to 1.
 Plugging \(C_1\) and \(C_2\) into~\eqref{eq-w-alpha-0} gives the final solution
 \begin{equation*}
 w\left(\alpha,\beta,\delta,t\right) =
 \cosh\left(iu\left(\alpha_0-\alpha\right)\right)
 \exp\left(k\delta - i\omega t\right)
 \exp\left(iv\beta\right).
 \end{equation*}
 There are two exponents in this solution with the opposite signs before
 horizontal coordinate \(\alpha\). These exponents denote
 incident and reflected wave respectively.
 The amplitude of the reflected wave does not decay as we go farther
 from the boundary, but decay only when we go deeper in the ocean. This behaviour
 corresponds to the real-world ocean waves.
 
 \subsubsection{Infinite plate}
 
 On a horizontal surface the boundary condition is written as
 \begin{equation*}
 \frac{d}{d t} \vec\nabla w \cdot \vec{n} =
 \frac{d}{d t} \frac{\partial}{\partial\delta} w = 0;
 \qquad
 \delta = \delta_0;
 \qquad
 \vec{n} = \VectorR{0\\0\\1}.
 \end{equation*}
 Analogously to wave direction we write vector form of the incident particle
 trajectory radius damping coefficient as \((0,0,k)\), hence vector form of the
 reflected coefficient is \((0,0,-k)\). We seek solutions of the
 form
 \begin{equation}
 \label{eq-w-delta-0}
 w\left(\alpha,\beta,\delta,t\right) =
 \left[ C_1\exp\left(k\delta\right) + C_2\exp\left(-k\delta\right) \right]
 \exp\left(- i\omega t\right) \exp\left(iu\alpha + iv\beta\right).
 \end{equation}
 We plug this expression into the boundary condition and get
 \begin{equation*}
 C_1\exp\left(k\delta_0\right) - C_2\exp\left(-k\delta_0\right) = 0.
 \end{equation*}
 Hence \(C_1=C\exp\left(-k\delta_0\right)\) and
 \(C_2=C\exp\left(k\delta_0\right)\).  Constant \(C\) may take arbitrary values,
 here we set it to \(1/2\). Plugging \(C_1\) and \(C_2\)
 into~\eqref{eq-w-delta-0} gives the final solution
 \begin{equation*}
 w\left(\alpha,\beta,\delta,t\right) =
 \cosh\left(k\left(\delta-\delta_0\right)\right)
 \exp\left(-i\omega t\right)
 \exp\left(iu\alpha + iv\beta\right).
 \end{equation*}
 There are two exponents in this solution with opposite signs before
 vertical coordinate \(\delta\). These exponents make the radius of
 the particle trajectories decay exponentially while approaching the boundary 
 \(\delta_0\) (i.e.~with increasing water depth). This is known
 solution from linear wave theory.
 
 \subsubsection{Infinite panel}
 
 On an arbitrary aligned infinite surface the boundary condition is written as
 \begin{equation*}
 \frac{d}{d t} \vec\nabla w \cdot \vec{n} = 0;
 \qquad
 \vec{n} \cdot \left(\vec\zeta - \vec\zeta_0\right) = 0,
 \end{equation*}
 where \(\vec{\zeta_0}\) is the point on the boundary plane and the third
 component of the normal vector is nought: \(\vec{n}=(n_1,n_2,0)\),
 \(\left|\vec{n}\right|=1\).
 The direction of incident wave is \(\vec{d}_i=(u,v,0)\) and the
 direction of reflected wave is \(\vec{d}_r\).  We seek solutions of the form
 \begin{equation}
 \label{eq-solution-diffraction}
 \begin{aligned}
 w\left(\alpha,\beta,\delta,t\right) = &\:
 C_1\exp\left( \left(i\vec{d}_i+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega_1
 t\right) \\ + &\:
 C_2\exp\left( \left(i\vec{d}_r+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega_2 t\right).
 \end{aligned}
 \end{equation}
 We plug this expression into the boundary condition and get
 \begin{equation*}
 \left( i\vec{d_i}\cdot\vec{n} \right) C_1
 +
 \left( i\vec{d_r}\cdot\vec{n} \right) C_2
 \exp\left( -i\vec{d}_s\cdot\vec\zeta_0 \right) = 0.
 \end{equation*}
 Here we substitute \(\vec{d}_r\cdot\vec{n}\) with \(-\vec{d}_i\cdot\vec{n}\)
 which is derived from the formula for \(\vec{d}_r\).
 Hence, the boundary condition reduces to
 \begin{equation*}
 C_1 - C_2 \exp\left( -i\vec{d}_s\cdot\vec\zeta_0 \right) = 0.
 \end{equation*}
 Hence
 \(C_1=\frac{1}{2}\exp\left(-\frac{1}{2}i\vec{d}_s\cdot\vec\zeta_0\right)\)
 and \(C_2=\frac{1}{2}\exp\left(\frac{1}{2}i\vec{d}_s\cdot\vec\zeta_0\right)\).
 This solution reduces to the solution for the wall when \(\vec{n}=(0,0,1)\).
 
 In a computer programme it is more practical to set \(C_1=1\) and
 \(C_2=\exp\left(i\vec{d}_s\cdot\zeta_0\right)\): that way you have to integrate
 only the second term in the solution over all ship hull panels
 (see~sec.~\ref{sec-results-diffraction}).
 
 %\input{stationary-surface.tex}
 %\input{progressively-moving-surface.tex}
 
 %\subsection{Non-stationary surface boundary}
 
 \subsection{OpenCL implementation}
  
 Solution for fluid velocity field was implemented in velocity potential solver
 in the framework of Virtual testbed.  Virtual testbed is a programme for
 workstations that simulates ship motions in extreme conditions and physical
 phenomena that causes them: ocean waves, wind, compartment flooding etc.  The
 main feature of this programme is to perform all calculations nearly in real
 time, paying attention to the high accuracy of calculations, which is partially
 achieved using graphical accelerators.
 
 Virtual testbed uses several solvers to simulate ship motions.
 The algorithm for velocity potential solver is the following.
 \begin{itemize}
     \item First of all, we generate wavy surface, according to our solution and using
         wetted ship panels from the previous time step (if any).
     \item Second, we compute wetted panels for the current time step, which are
         located under the surface calculated on the previous step.
     \item Finally, we calculate Froude---Krylov forces, acting on a ship hull.
 \end{itemize}
 These steps are repeated in infinite loop. Consequently, wavy surface is 
 always one time step behind the wetted panels. This inconsistency is a result
 of the decision not to solve ship motions and fluid motions in one system of
 equations, which would be too difficult to do.
 
 Let us consider process of computing wavy surface in more detail.  Since wavy
 surface grid is irregular (i.e.~we store a matrix of fluid particle positions
 that describe the surface), we compute the same formula for each point of the
 surface.  It is easy to do with C++ for CPU computation, but it takes some
 effort to efficiently run this algorithm with GPU acceleration.  Our first
 naive implementation was inefficient, but the second implementation that used
 local memory to optimise memory loads and stores proved to be much more
 performant.
 
 First, we optimised storage order of points making it fully sequential.
 Sequential storage order leads to sequential loads and stores from the global
 memory and greatly improves performance of the graphical accelerator.  Second,
 we use as many built-in vector functions as we can in our computations, since
 they are much more efficient than manually written ones and compiler knows how
 to optimise them. This also decreases code size and prevents possible mistakes
 in the manual implementation. Finally, we optimised how ship hull panels are
 read from the global memory. One way to think about panels is that they are
 coefficients in our model, as array of coefficients is typically read-only and
 constant. This type of array is best placed in the constant memory of the
 graphical accelerator that provides L2 cache for faster loads by parallel
 threads. However, our panel array is too large to fit in constant memory, so we
 simulated constant memory using local memory: we copied a small block of the
 array into local memory of the multiprocessor, computed sum using this block
 and then proceeded to the next block. This approach allowed to achieve almost
 200-fold speedup over CPU version of the solver.
 
 A distinctive feature of the local memory is that it has the smallest latency,
 at the same time sharing its contents between all computing units of the
 multiprocessor.  Using local memory we reduce the number of load/store
 operations to global memory, which has larger latency.  As far as global memory
 bandwidth remains a bottleneck, this kind of optimisation would improve
 performance.  To summarise, our approach to write code for graphical
 accelerators is the following:
 \begin{itemize}
     \item make storage order linear,
     \item use as many built-in vector operations as is possible,
     \item use local memory of the multiprocessor to optimise global memory
         load and stores.
 \end{itemize}
 Following these simple rules, we can easily implement efficient algorithms.
 
 \section{Results}
 
 \subsection{Diffraction}
 \label{sec-results-diffraction}
 
 In the first experiment we use solution for infinite panel to simulate wave
 diffraction around Aurora's ship hull (see~fig.~\ref{fig-ships}). In order to
 apply~\eqref{eq-solution-diffraction} to this problem we use smoothing kernel
 that accumulates influence of every panel on a particular point of ocean
 surface:
 \begin{equation*}
     w = \sum\limits_{j} K_j w_j;
     \qquad
     K_j = \frac{1}{1 + \left|\zeta-\zeta_0\right|^2}
 \end{equation*}
 Here \(w_j\) is solution~\eqref{eq-solution-diffraction} written for panel
 \(j\), \(K_j\) is smoothing kernel, \(\zeta_0\) is the centre of the panel.
 In the centre of the panel \(K_j=0\) and far from the ship hull
 \(K_j\rightarrow{}0\).
 
 \begin{figure}
 	\centering
 	\includegraphics{build/ships/diogen.eps}\hfill
 	\includegraphics{build/ships/aurora.eps}\hfill
 	\includegraphics{build/ships/micw.eps}
 	\caption{Diogen, Aurora and MICW three-dimensional ship hull
 	models.\label{fig-ships}}
 \end{figure}
 
 In the experiment waves with amplitude 1 approach the ship from the aft.  The
 results of the experiment are shown in fig.~\ref{fig-diffraction}.  Near the
 aft waves change their direction to be tangent to the waterline curve, follow
 the curve to the bow, and then restore their original direction.  The amplitude
 of waves near the hull is also increased.
 
 \begin{figure}
     \centering
     \includegraphics{build/gnuplot/surface.eps}
     \caption{Wave diffraction around Aurora's hull (the hull is not
     shown). Top view.\label{fig-diffraction}}
 \end{figure}
 
 
 \subsection{Performance benchmarks}
 
 We implemented velocity potential solver using OpenMP for parallel computations
 on a processor and OpenCL for graphical accelerator.  The solver uses single
 precision floating point numbers. Benchmark results are presented in
 tab.~\ref{tab-benchmark}.
 
 We performed benchmarks for three ships: Diogen, Aurora and MICW.  Diogen is a
 small-size fishing vessel, Aurora is mid-size cruiser and MICW is a large-size
 ship with small moment of inertia for the current waterline
 (fig.~\ref{fig-ships}). The main difference between the ships that affects
 benchmarks is the number of panels into which the hull is decomposed. These
 numbers are shown in tab.~\ref{tab-ships}.
 
 \begin{table}
 	\centering
 	\caption{Parameters of ship hulls that were used in the
 	benchmarks.\label{tab-ships}}
 	\begin{tabular}{lrrr}
 		\toprule
 					  & Diogen & Aurora & MICW  \\
 		\midrule               
 		Length, m     &   60   & 126.5  &   260 \\
 		Beam, m       &   15   &  16.8  &    32 \\
 		Depth, m      &   15   &  14.5  &    31 \\
 		No. of panels & 4346   &  6335  &  9252 \\
 		\bottomrule
 	\end{tabular}
 \end{table}
 
 Benchmarks were performed using three workstations: DarkwingDuck, GPUlab,
 Capybara. DarkwingDuck is a laptop, GPUlab is a desktop workstation, and
 Capybara is a desktop with professional graphical accelerator server-grade
 processor (tab.~\ref{tab-config}). 
 
 \begin{table}
     \centering
     \caption{Performance benchmarks results. Numbers represent average time in milliseconds
     that is needed to generate waves with reflection.\label{tab-benchmark}}
     \begin{tabular}{l@{\hskip 3mm}r@{\hskip 2mm}r@{\hskip 3mm}r@{\hskip 2mm}r@{\hskip 3mm}r@{\hskip 2mm}r}
 		\toprule
 				 & \multicolumn{2}{l}{Diogen}
 				 & \multicolumn{2}{l}{Aurora}
 				 & \multicolumn{2}{l}{MICW} \\
 		\cmidrule(r){2-3}
 		\cmidrule(r){4-5}
 		\cmidrule(r){6-7}
 		Node     & MP & CL 
                  & MP & CL 
                  & MP & CL  \\
 		\midrule
 		DarkwingDuck  
                  &    5462   & 48
                  &    7716   & 41
                  &    7725   & 11   
                  \\
 		GPUlab  
 		         &    5529   & 11                    
                  &    8222   & 10   
 		         &    6481   & 3
                  \\
 		Capybara &    2908   & 16                    
 		         &    2091   & 8     
 		         &    2786   & 4     
                  \\
 		\bottomrule
 	\end{tabular}
 \end{table}
 
 \begin{table}
 	\centering
 	\caption{Hardware configurations for benchmarks. For all benchmarks we
 	used GCC version 9.1.0 compiler and optimisation flags \texttt{-O3
 	-march=native}.\label{tab-config}}
 	\begin{tabular}{lllrr}
 		\toprule
 		     &     &     & \multicolumn{2}{c}{GPU GFLOPS} \\
 		\cmidrule(r){4-5}
 		Node           & CPU              & GPU            & Single & Double \\
 		\midrule
         DarkwingDuck   & Intel i7-3630QM  & NVIDIA GT740M  &  622   &     \\
 		GPUlab         & AMD FX-8370      & NVIDIA GTX1060 & 4375   & 137 \\
 		Capybara       & Intel E5-2630 v4 & NVIDIA P5000   & 8873   & 277 \\
 		\bottomrule
 	\end{tabular}
 \end{table}
 
 
 \section{Discussion}
 
 All the solutions obtained for various boundaries in this papre satisfy
 continuity equation and equation of motion, but they are all written for
 \emph{plane} surface boundaries with different orientations. Typical ship hull
 three-dimensional model is represented by triangulated surface, and in the
 centre of each triangular panel fluid particle velocity vector does not depend
 on the surface normal of the other panels. So, the solution for plane surface
 boundary is enough to compute fluid velocity field directly \emph{on} the
 surface boundary.
 
 In order to generalise the solution for fluid velocity field \emph{near} the
 surface boundary, we need to calculate weighted average of reflection terms of
 each underwater panel of the surface. Using inverse squared distance as the
 weight gives acceptable results in our experiments, but may not be appropriate
 in others.
 
 It is not clear how much we have to increase wave amplitude near the boundary.
 One way to control the amptlitude increase is to introduce the coefficient
 \(C\) into the reflection formula \(\vec{d}_r=\vec{d}_i-C\vec{d}_s\) and before
 reflection term in~\eqref{eq-solution-diffraction}.  When \(C=1\) the wave is
 fully reflected from the boundary and the amplitude is doubled, when \(C=0\) no
 reflection occurres and the amplitude does not change.
 
 Performance benchmarks showed that graphical accelerator greatly improves
 performance of velocity potential solver. Linear memory access patterns and
 large amount of floating point operations make this solver an ideal candidate
 for running on a graphical accelerator, and these features are the result of
 deriving explicit solution for fluid motions near the ship hull boundary.
 
 \section{Conclusion}
 
 This paper proposes a new model for ocean wave diffraction near ship hull.
 This model uses law of reflection to simulate incident and reflected waves and
 fluid particle motion. Although, the solutions are written for infinite plane
 surfaces, they can be used to compute fluid velocity at the centre of each
 triangle of the triangulated ship hull surface and can be generalised to
 compute fluid velocity near the ship hull using weighted sum over all panels.
 
 The model was implemented in Virtual testbed velocity potential solver and was
 found to be highly efficient on a graphical accelerator.  Future work is to
 incorporate radiation into the model and compare the solution to existing
 empirical approaches.
 
 \subsubsection*{Acknowledgements.}
 Research work is supported by Saint Petersburg State University (grants
 no.~51129371 and~51129725) and Council for grants of the President of the
 Russian Federation (grant no.~MK-383.2020.9).
 
 \bibliographystyle{splncs04}
 \bibliography{references}
 
 \end{document}