[book_gpu.git] / BookGPU / Chapters / chapter7 / ch7.tex

\chapterauthor{Allan P. Engsig-Karup}{Technical University of Denmark}
\chapterauthor{Stefan L. Glimberg}{Technical University of Denmark}
\chapterauthor{Allan S. Nielsen}{Technical University of Denmark}
\chapterauthor{Ole Lindberg}{Technical University of Denmark}

\chapter{Fast hydrodynamics on heterogenous many-core hardware}
\label{ch7}

\begin{figure}[!htb]
\centering
\includegraphics[width=0.95\textwidth]{Chapters/chapter7/figures/figSeries60CB06Type7StedaySnapshot-eps-converted-to.pdf}
%\caption{Snapshot of steady state wave field generated by a Series 60 ship hull.}
\end{figure}

\newpage
In this chapter, we present details of a heterogenous and massively parallel GPU library implementation in CUDA C/C++ of a  nonlinear free surface water wave model \cite{ch7:EngsigKarupEtAl2011}. We describe how flexible-order finite difference\index{finite difference} approximations to the partial differential equations of the model can be prototyped using library components provided in an in-house library. In this library hardware-specific implementation details are hidden via template-based components, as described in chapter \ref{ch5}. We provide details of the modelling basis and important unique numerical properties which has been made tuneable to create a powerful tool that can be tailored for specific purposes in engineering analysis. The model is based on unified potential flow theory, and can be applied in scientific applications related to maritime engineering. It can be applied for cost-efficient estimation of wave propagation and transformation of irregular multidirectional waves over variable depth, kinematics and structural wave loads over large areas and scales.

A main motivation of this work is to deliver exceptional performance to minimize calculation times, using modern parallel hardware technologies in combination with a proper choice of discretization methods and data-local algorithms. Data-local algorithms with optimal complexity, such that work and memory requirements grow (scale) linearly with problem size on any hardware system. For the wave model this is achieved by explicit time integration and iterative solution of a large non-symmetric and sparse linear $\sigma$-transformed Laplace problem. For the latter, we use an iterative Preconditioned Defect Correction (PDC) method, accelerated using a geometric multigrid preconditioning strategy. We use modern programming paradigms in the form of MPI and CUDA for development of a novel massively parallel wave modelling tool, executable on modern heterogenous many-core hardware.

One purpose of the developed numerical model is to perform hydrodynamic calculations in computationally intensive interactive real-time simulations. Realistic simulations are, with present technology and available computational resources, a tremendous challenge in this setting. Yet, our aim is to take a first step in this direction, and compute first-order accurate hydrodynamics for near-realistic simulations of unsteady ship-wave dynamics in a large ship simulator, used for training purposes in seakeeping operations. For this type of application, a mandatory ingredient for real-time and interactive simulation is a truly high-performance parallel implementation to ensure data processing in time for interactive visualization and responses. Details of the model properties, implementation, and promising novel combinations of techniques and algorithms for acceleration of performance are presented. Numerical experiments and benchmarks are provided to demonstrate the accuracy and efficiency of the model across recent generations of many-core CUDA-enabled hardware.

\section{On hardware trends and challenges in scientific applications}

During the last two decades we have seen how computer graphics hardware has been developed from fixed pipeline processors with no level of programmability, to flexible high-performance hardware platforms, suitable for general purpose scientific computations other than computer graphics. This trend has been disruptive high-performance computing on mass-produced commodity hardware and give new opportunities for computational science and engineering on work stations for a broad range of scientific applications. This emphasizes the increasingly important role of computers in simulation of real world dynamics \cite{ch7:Keyes201170}. In recent years, the Compute Unified Device Architecture (CUDA) programming model, based on the standard C/C++ programming language and introduced by Nvidia, has become popular as a proprietary and widely used standard in high performance communities. It is by design and supported functionality, easy and sufficient to be deployed for wide improvement of existing and new applications across science and engineering fields, that can benefit from the the use of heterogenous hardware.

We should be careful about speculating about the future and extrapolate from current trends. The TOP 500 list\footnote{\url{http://www.top500.org.}} of supercomputers shows that there are some general noticeable hardware trends and give indication of what to expect in the near future. First, since 2005 we have seen how power constraints and resulting heat dissipation problems forced chip producers to increase the number of cores rather than clock frequency.  Multi-core processors have become the new standard and many-core processors are becoming available as a standard in commodity hardware, from personal laptops to super-computing clusters.

This trend suggests, that there will be less fast low-latency memory available per core in the future, favoring data-locality. In addition, we have also seen how communication speed to computation speed ratio decreases, making it increasingly difficult to supply data to hungry floating point units. In addition, there will likely be increasing amounts of data to store as a result of increasing processing capabilities. The rapidly increasing floating point performance following Moore's law for transistor production has resulted in a significant memory gap which leaves most PDE-based scientific applications bandwidth bound rather than compute bound. This trend is driven by pure commercial needs and not the needs of high-performance computing. Roads to better performance includes standardization of software infrastructure, rethinking algorithms to better exploit memory hierarchies optimally to boost strong scaling properties, increase locality in algorithms and introduce as much concurrency and work as possible to both utilize and exploit the many cores. Also, software that can utilize many cores should be fault-tolerant to maximize time to solution for application users. We should also expect to see multiple layers of parallelism that will have to be exploited and possibly auto-tuned to optimally utilize hardware. This introduces new challenges in compilers, requires programming experts with hardware knowledge and introduces new trends in software developments to leverage productivity and utilize available computational resources in more optimal ways. We have observed a fundamental paradigm shift of underlying hardware design towards much more heterogeneity and parallelism.

A key problem is, that improvements in performance require porting legacy codes\footnote{In the worst case, a legacy code is an undocumented serial code developed by a developer no longer around and for a long time ago for execution on single core hardware.} to new hardware, and possibly changing algorithms which have been developed for the conventional single core processors decades ago. Without this, it may be impossible to utilize and scale algorithmic work optimally to achieve high performance on modern and emerging hardware. This problem is currently addressed with rapid progress by researchers and industry by development of new optimized libraries, that can utilize such new hardware at minimum effort. While we have seen significant improvements in such efforts, and today see much more rapid development of applications, there are still few applications running entirely on heterogenous hardware. However, increasing amounts of applications are utilizing accelerators to parts of their code to gain speedups albeit with less dramatic improvements of performance as one can potentially find by adapting most, if not all of the application to modern hardware.

In this work, we explore some of these trends by developing, by bottom-up-design, a water wave model which can be utilized in maritime engineering and with intended use on affordable office desktops as well as on more expensive modern compute clusters for analysis purposes.  
\section{On modeling paradigms for highly nonlinear and dispersive water waves}
\label{ch7:sec:modernwavemodellingparadigms}

We see development of new hardware technologies as a key driver for exploring new and revisiting existing approaches that can contribute to next-generation modelling techniques.

For instance, the dominant wave modelling paradigm today for numerical simulation in coastal engineering tools is the use of Boussinesq-type \index{Boussinesq models} formulations for approximate solution of unified potential flow\index{potential flow} equations over varying bathymetry \cite{ch7:MS98}. The use of Boussinesq-type models in engineering tools was pioneered by Abott et. al. (1978) \cite{ch7:AbottEtAl1978,ch7:AbottEtAl1984} based on the original Boussinesq equations due to Peregrine (1967) \cite{ch7:Peregrine1967} for calculations of waves in a harbour area. New formulations for highly nonlinear and dispersive water waves, useful for description of wave propagation in the important application range from deep to shallow areas, have been subject of intense research for more than two decades. Such higher order formulations can be derived by first introducing an infinite Mclaurin series solution to the Laplace equation as described in \cite{ch7:AMS99}. This technique was later generalized to arbitrary expansion levels \cite{ch7:MBL02}. By analytical truncation of such series solutions, a polynomial variation in the vertical is assumed, and provides the basis for efficient higher order Boussinesq-type formulations \cite{ch7:MBS03,ch7:Bingham2009467} for fully nonlinear and dispersive water waves. It is attractive, since it is then possible to eliminate the vertical coordinate in the analytical formulation of the Laplace problem. The resulting approximate model typically contains higher order derivatives that require treatment in numerical models. Thus, this truncation procedure inherently limits the practical application range, however, can be significantly improved by optimization via Pad\'e approximations together with introduction of free parameters for extending the application range via optimization of accuracy in dispersion, kinematics and shoaling characteristics.

Main challenges of Boussinesq-type models are accurate and large-scale simulation of waves propagating towards near-shore from deep to shallow waters through surf zones, while accounting for high-order dispersive and nonlinear effects \cite{ch7:Cavaleri2007603}. Within the last two decades, much research has focused on extending the application range through improved formulations in terms of both dispersion, shoaling, kinematic and nonlinear properties. The ultimate high-order Boussinesq-type model due to \cite{ch7:MBS03} was at the time considered a breakthrough in this direction, and since then new promising formulations have been proposed. For example, the methodology behind Boussinesq-type formulations can be extended via a multi-layer approach \cite{ch7:LynettEtAl2004a,ch7:LynettEtAl2004b,ch7:ChazelEtAl2010}, that makes it possible to achieve a similar range of application and levels of accuracy, but without higher derivatives in the formulation that can cause numerical difficulties.

Boussinesq-type formulations for free surface waves are conventionally evaluated against the unified potential flow theory to evaluate limits to application range and accuracy limits. The use of unified potential theory as a basis for numerical models has traditionally been perceived too expensive \cite{ch7:Lin2008} to solve in comparison with the typically more efficient Boussinesq-type models. This may be true in a strict comparison between the models, especially with respect to applications towards the more restricted shallow regions. However, this is despite that a numerical unified potential flow model can be used for a larger range of practical scientific applications. A unified potential flow model has at most second order derivatives in the formulation. In a numerical setting it has good opportunities for balancing accuracy and work effort by appropriate tuning of discrete parameters. This comes without a need for changing the underlying wave model to extend application range towards deep waters. Thus, the main problem related to the practical use of a unified model in maritime applications is arguably an issue of numerical efficiency.

To address this issue, we have recently proposed a proof-of-concept approach, that combines modern many-core hardware with appropriate numerical and parallel strategies to facilitate efficient, accurate and scalable solution of water wave problems~\cite{ch7:EngsigKarupEtAl2011}. The use of potential theory for unsteady water wave computations can be traced at least back to \cite{ch7:HausslingVanEseltine1975}, and the fully nonlinear potential equations have been solved using various numerical methods since then, e.g., see reviews \cite{ch7:Yeung1982,ch7:TsaiYue1996,ch7:DiasBridges2006,ch7:Lin2008}. In the context of the finite-difference method, an efficient and scalable second-order geometric multigrid approach was first proposed by Li \& Fleming (1997) \cite{ch7:LiFleming1997}. Since then, the numerical strategy has been significantly improved in several works \cite{ch7:BinghamZhang2007,ch7:EBL08} via more efficient discretization techniques, with the objective to develop an efficient general purpose strategy, that can be used for a broad range of practical maritime applications. Recently, a comparison with a High-Order Spectral (HOS) model \cite{ch7:DucrozetEtAl2011} was also reported to assess  accuracy and relative differences in efficiency on single-core hardware against a superior spectral modelling basis for a numerical wave tank setup in a structured domain with a flat sea bed.

\section{Governing equations}
\label{ch7:goveq}

We describe how, by physical principles via mathematical procedures and assumptions, it is possible to formulate a fully nonlinear and dispersive water wave model, describing incompressible, irrotational and inviscid fluid\index{fluid} flow above an uneven seabed.

Conservation of mass\index{mass conservation} for an infinitely small control volume can be stated as
\begin{align}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho {\bf u}) = 0,
\end{align}
where $\rho$ is fluid density, and ${\bf u}=(u,v,w)$ is the velocity field vector and $\nabla=(\partial/\partial x, \partial/\partial y, \partial/\partial z)$ a Cartesian gradient operator. If we assume that fluid density is constant, i.e. that the fluid is incompressible, the mass continuity equation simplifies to
\begin{align}
\label{ch7:conteq}
 \nabla \cdot {\bf u} = 0,
\end{align}
which shows that the divergence of the velocity field is zero everywhere.

Conservation of momentum\index{momentum conservation} for an infinitely small control volume is expressed by the famous Navier-Stokes equations
\begin{align}
\label{ch7:NSeq}
\rho \frac{D{\bf u}}{Dt} = -\nabla p + \mu \nabla^2{\bf u} + {\bf F},
\end{align}
where $p$ is pressure and ${\bf F}$ is the net force vector acting on the fluid volume, assumed to be of the form ${\bf F}=\rho{\bf g}$ with ${\bf g}=(0,0,-g_z)$ accounting for gravitational effects in the vertical direction. This implies that surface tension effects on the free surface is neglected. Exact solutions to the Navier-Stokes equations are in general difficult to obtain and this motivate the use of numerical methods for direct simulation of fluid flow and if necessary analytical simplifications to only account for physics of interest.

The material derivative for a co-moving coordinate system used in Lagrangian formulations
\begin{align}
\frac{D}{Dt}\equiv \frac{\partial }{\partial t} + ({\bf u}\cdot \nabla),
\end{align}
is defined as the sum of a time derivative and a convective term measured in a static (Eularian) coordinate system and accounts for the time rate of change following the motion. Thus, for the velocity vector ${\bf u}$, the total acceleration is defined as
\begin{align}
\label{ch7:momentum}
\frac{D{\bf u}}{Dt}\equiv \frac{\partial {\bf u}}{\partial t} + \frac{1}{2}\nabla {\bf u}^2 +{\boldsymbol \omega} \times {\bf u},
\end{align}
where the curl of the velocity field ${\boldsymbol \omega}\equiv\nabla\times {\bf u}$ is referred to as the vorticity vector field accounting for rotation of fluid particles. If we assume that the flow is irrotational
\begin{align}
\nabla\times {\bf u} = 0,
\end{align}
and make use of the following relationship known from vector calculus
\begin{align}
\frac{1}{2}\nabla({\bf u}\cdot {\bf u}) = ({\bf u}\cdot \nabla){\bf u} + {\bf u}\times (\nabla \times {\bf u}),
\end{align}
we can rewrite \eqref{ch7:momentum} into
\begin{align}
\label{ch7:momentum2}
\frac{D{\bf u}}{Dt} = \frac{\partial {\bf u}}{\partial t} +  \frac{1}{2}\nabla({\bf u}\cdot {\bf u}).
\end{align}
Since the Navier-Stokes equations \eqref{ch7:NSeq} can be interpreted as the application of Newton's second law to an infinitely small fluid volume, we can establish that the changes in momentum in an infinitely small control volume of a fluid is simply the sum of forces due to pressure gradients, dissipative forces, gravitational forces and possibly other forces acting inside the fluid volume.

To accurately simulate propagation of long gravity waves and high-reynolds number flows in the context of maritime applications it is for many applications acceptable to assume that viscous forces are small in comparison with inertial forces. In this case, it is reasonable to assume that the fluid is inviscid and we can neglect the viscous terms in \eqref{ch7:NSeq}. Then we obtain the set of Euler equations
\begin{align}
\label{ch7:momentum3}
\frac{D{\bf u}}{Dt} = - \frac{1}{\rho}\nabla p + \frac{1}{\rho} {\bf F},
\end{align}
which does not account for any losses in energy via dissipative physical mechanisms.

If we introduce a scalar velocity potential function $\phi$ that relates to the velocity components in the following way
\begin{align}
\label{ch7:vectorfield}
{\bf u}\equiv \nabla \phi = \left(\frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y}, \frac{\partial \phi}{\partial z} \right),
\end{align}
the number of unknowns can be lowered and we find that the scalar velocity potential function due to \eqref{ch7:conteq} must satisfy the Laplace equation
\begin{align}
\label{ch7:laplaceeq}
\nabla^2\phi=0,
\end{align}
Solutions to this equation are completely determined by the boundary conditions\index{boundary condition}. Thus, it is possible, given appropriate boundary conditions, to determine the scalar velocity potential $\phi$ in all of the domain by solving the resulting Laplace problem. When the scalar velocity potential function is known, detailed information of the kinematics can immediately be obtained.

With the definition of the vector field \eqref{ch7:vectorfield}, we can collect the momentum equations \eqref{ch7:momentum2} and \eqref{ch7:momentum3} and express the momentum equation as
\begin{align}
\rho\left( \frac{\partial \nabla \phi}{\partial t} + \nabla \left(\frac{1}{2} \nabla\phi \cdot \nabla\phi\right)  \right) = -\nabla p - \nabla(\rho g z),
\end{align}
having assumed that the net force can be decomposed into pressure and gravity forces only. This set of equations can be rewritten into
\begin{align}
\nabla \left[ \rho \frac{\partial \phi}{\partial t} + \rho \frac{1}{2} \nabla\phi \cdot \nabla\phi + p + \rho g z \right] = 0,
\end{align}
and by integration in space we arrive at the unsteady Bernoulli's Equation
\begin{align}
\rho \frac{\partial \phi}{\partial t} + \rho \frac{1}{2} \nabla\phi \cdot \nabla\phi + p + \rho g z = G(t),
\end{align}
where $G(t)$ is an arbitrary function of the integration that can be assumed to be zero as it only defines a reference value for the unphysical scalar velocity potential function. Bernoulli's equation is typically used as a dynamic condition for the free fluid surface, expressing that the fluid pressure at the free surface is equal to the pressure in the air above the free surface.

In the following, we assume that the displacement of the free surface $z=\eta({\bf x},t)$ is described in a Cartesian coordinate system with the $xy$-plane located at the still water level and the positive $z$-axis pointing upwards. It is typical to assume a constant atmospheric pressure at the free surface by defining $p=0$ as reference. This leaves us with a dynamic boundary condition for the free surface velocity potential function stated as
\begin{align}
\frac{\partial \phi}{\partial t} = -\frac{1}{2} \nabla\phi \cdot \nabla\phi - g z, \quad z=\eta.
\end{align}
At the free surface, we can determine a kinematic free surface condition by determining the rate of change of the streamline for the surface as
\begin{align}
\frac{\partial \eta}{\partial t} = -\frac{\partial \phi}{\partial x} \frac{\partial \eta}{\partial x} - \frac{\partial \phi}{\partial y} \frac{\partial \eta}{\partial y} +  \frac{\partial \phi}{\partial z}, \quad z=\eta.
\end{align}
Spatial and temporal differentiation of the free surface variables are related by the chain rule
%
\begin{align}
\boldsymbol{\nabla}\tilde{\phi} &= (\boldsymbol{\nabla}\phi)_{z=\eta} + \left(\frac{\partial \phi}{\partial z}\right)_{z=\eta}\boldsymbol{\nabla}\eta,  \\
\frac{\partial \tilde{\phi}}{\partial t}  &= \left( \frac{\partial \phi}{\partial t} \right)_{z=\eta} +\frac{\partial  \eta}{\partial t} \left(\frac{\partial \phi}{\partial z} \right)_{z=\eta},
\end{align}
%
and can be used to transform the free surface problem to variables defined solely at the free surface as
%
\begin{subequations}
\begin{align}
\frac{\partial}{\partial t} \eta &= -\boldsymbol{\nabla}\eta\cdot\boldsymbol{\nabla}\tilde{\phi}+\tilde{w}(1+\boldsymbol{\nabla}\eta\cdot\boldsymbol{\nabla}\eta), \label{ch7:FSeta} \\
\frac{\partial}{\partial t} \tilde{\phi} &= -g\eta - \frac{1}{2}\left(\boldsymbol{\nabla}\tilde{\phi}\cdot\boldsymbol{\nabla}\tilde{\phi}-\tilde{w}^2(1+\boldsymbol{\nabla}\eta\cdot\boldsymbol{\nabla}\eta)\right).
\label{ch7:FSphi}
\end{align}
\label{ch7:FSorigin}
\end{subequations}
%
with the $\boldsymbol{\nabla}$-operator from here and forward defined as a horizontal gradient operator $\boldsymbol{\nabla}=(\partial_x,\partial_y)$ and tilde's are used for free surface variables. To solve the set of unsteady free surface equations \eqref{ch7:FSorigin}, we need a closure between the horizontal and vertical free surface velocity variables. This can be established by solving the Laplace equation \eqref{ch7:laplaceeq} in the interior domain together with suitable boundary conditions.

A kinematic bottom condition can be derived by assuming that the fluid particles follow a streamline along the solid sea bed. Consider the rate of change of such a streamline at still-water depth $z=-h({\bf x},t)$ and we find
\begin{align}
\label{ch7:kinbot}
\frac{\partial z}{\partial t} = - \frac{\partial h}{\partial x}  \frac{\partial \phi}{\partial x} - \frac{\partial h}{\partial y}  \frac{\partial \phi}{\partial y} - \frac{\partial h}{\partial t}, \quad z=-h({\bf x},t).
\end{align}
We assume that the sea bed is static allowing us to neglect the last term. Thus, by specifying $\tilde\phi$ as a Dirichlet condition at the free surface together with a kinematic bottom boundary condition at the sea bed defines a Laplace problem
%
\begin{subequations}
\label{ch7:eq:laplaceproblem}
\begin{align}
\phi & =  \tilde{\phi}, \quad z = \eta({\bf x},t), \\
\boldsymbol{\nabla}^2\phi + \partial_{zz}\phi & = 0, \quad -h\leq z<\eta({\bf x},t), \label{ch7:Laplace} \\
\partial_z \phi + \boldsymbol{\nabla}h\cdot\boldsymbol{\nabla}\phi &= 0, \quad z=-h. \label{ch7:KB}
\end{align}
\end{subequations}
where we have used that $\partial_t z \equiv \partial_z \phi$ to rewrite the first term of the kinematic bottom condition \eqref{ch7:kinbot}.
%
The moving free surface makes the spatial fluid domain $\Omega$ vary in time. The main challenges in solving these equations numerically are to deal with the time-dependent fluid domain and nonlinearity of the equations.
%

\subsection{Boundary conditions}\index{boundary condition}

We consider three types of boundaries, namely, fully reflective boundaries, incident wave boundaries and absorbing boundaries. The fully reflective boundaries are handled through numerical approximations of the boundary conditions for solid walls and bottom surfaces stating that the velocity in the normal direction is zero
\begin{align}
{\bf n}\cdot \nabla \phi = 0, \quad {\bf x}\in\partial\Omega,
\end{align}
where ${\bf n}=(n_x,n_y)$ is a two-dimensional normal vector pointing outwards from the solid surface. A complementary condition for the free surface elevation variable is
\begin{align}
{\bf n} \cdot \nabla\eta = 0.
\end{align}

Incident wave and absorbing boundary conditions are imposed via an embedded penalty forcing technique as described in section \ref{ch7:wavegen}.

\section{The numerical model}
\label{ch7:sec:nummodel}

The unified potential flow model is attractive as a basis due to the underlying analytical properties.
From a numerical point of view, an efficient and scalable discretization strategy should be based on using a data-local method, e.g., a flexible-order finite difference method for discretely approximating the governing equations and imposing boundary conditions via fictitious ghost points techniques as described in \cite{ch7:BinghamZhang2007,ch7:EBL08}. Such an approach has several attractive features from a scientific computing perspective. For example, finite difference methods are among the simplest and most efficient methods due to the use of structured grids and data structures. This result in low implementation and computational complexity which maps efficiently to modern computer architectures. Formal accuracy and tuneable numerics are achieved by employing flexible-order finite difference\index{finite difference} (local stencil)\index{stencil} approximations. 

We present scalability and performance tests based on the same two test environments outlined in chapter \ref{ch5} section \ref{ch5:sec:testenvironments}, plus a a third test environment based on the most recent hardware generation:
\begin{description}
\item[Test environment 3.] Desktop with dual-socket Sandy Bridge Intel Xeon E5-2670 (2.60 GHz) processors, 64GB RAM, 2x Nvidia Tesla K20 GPUs.
\end{description}
Performance results can be used to predict actual runtimes as described in \cite{ch7:EngsigKarupEtAl2011}, e.g., for estimation of whether a real-time constraint for a given problem size can be met.

\subsection{Strategies for efficient solution of the Laplace problem}\index{Laplace problem}

As explained in section \ref{ch7:sec:modernwavemodellingparadigms}, for the formulation of potential flow problems there are two widely used paradigms for solving the Laplace problem efficiently. The most widely used approach is of Boussinesq-type where essentially the three-dimensional formulation is reduced to a two-dimensional formulation. The main argument for this type of model reduction procedure is the resulting efficiency in the numerical models. The price paid is typically high-order derivatives in the approximate formulation and is justified by efficient solution of an approximate Laplace problem.

A second approach, is to transform the equations at PDE level to provide a basis for efficient direct solution of the discrete Laplace problem for the entire volume. This strategy is based on a paradigm where approximations are done by discrete approximations rather than analytical manipulations of the equation at PDE level. This approach introduces at a first look more complexity in the formulation, e.g. by the introduction of mixed derivatives, however, essentially does not limit the application range beyond those resulting from numerical approximations and properties hereof. Using this second approach, it is standard to introduce a $\sigma$-transformation in the vertical coordinate
\begin{align}
\label{ch7:sigtrans}
\sigma \equiv \frac{z+h(\boldsymbol{x})}{d(\boldsymbol{x},t)}, \quad 0\leq \sigma \leq 1,
\end{align}
where $d(\boldsymbol{x},t)=\eta(\boldsymbol{x},t) + h(\boldsymbol{x})$ is the height of the water column above the bottom. This enables a transformation of the problem to a time-constant computational domain at the expense of time-varying coefficients.

The fluid domain is mapped to a time-invariant computational domain in which the Laplace problem \eqref{ch7:eq:laplaceproblem} is expressed as
%
\begin{subequations}
\label{ch7:TransformedLaplace}
\begin{align}
  \Phi & =  \tilde{\phi}, \quad \sigma = 1, \label{ch7:FSDirichlet} \\
  \boldsymbol{\nabla}^2\Phi+\boldsymbol{\nabla}^2\sigma(\partial_\sigma\Phi) +2\boldsymbol{\nabla}\sigma\cdot\boldsymbol{\nabla}(\partial_\sigma\Phi) + & \nonumber \\
  (\boldsymbol{\nabla}\sigma\cdot\boldsymbol{\nabla}\sigma+(\partial_z\sigma)^2)\partial_{\sigma\sigma}\Phi&=0, \quad 0\leq\sigma<1, \label{ch7:convlaplace}\\
{\bf n} \cdot ({\bf \nabla},\partial_z\sigma\partial_\sigma)\Phi & = 0, \quad ({\bf x},\sigma)\in\partial\Omega, \label{ch7:strucbcs}
%  (\partial_z\sigma+\boldsymbol{\nabla}h\cdot\boldsymbol{\nabla}\sigma)(\partial_\sigma\Phi) + \boldsymbol{\nabla}h\cdot\boldsymbol{\nabla}\Phi&=0, \quad \sigma=0, \label{eq:BBC}
\end{align}
\end{subequations}
%
where the scalar velocity function $\Phi(\boldsymbol{x},\sigma,t)=\phi(\boldsymbol{x},z,t)$ contains all information about the flow kinematics in the entire fluid volume. The spatial derivatives of the coordinate $\sigma$ appearing in \eqref{ch7:TransformedLaplace} are expressed as
%
\begin{subequations}
\label{ch7:timedepcoef}
\begin{align}
\boldsymbol{\nabla}\sigma & = \tfrac{1-\sigma}{d}\boldsymbol{\nabla}h
- \tfrac{\sigma}{d}\boldsymbol{\nabla}\eta, \\
\boldsymbol{\nabla}^2\sigma & = \tfrac{1-\sigma}{d}\left(\boldsymbol{\nabla}^2
h-\tfrac{\boldsymbol{\nabla} h\cdot\boldsymbol{\nabla} h}{d}\right)
-\tfrac{\sigma}{d}\left(\boldsymbol{\nabla}^2 \eta-\tfrac{\boldsymbol{\nabla}
\eta\cdot\boldsymbol{\nabla}\eta}{d}\right)  \nonumber \\
 &-\tfrac{1-2\sigma}{d^2}\nabla h\cdot\nabla\eta
-\tfrac{\boldsymbol{\nabla}\sigma}{d}\cdot\left(\boldsymbol{\nabla} h
+ \boldsymbol{\nabla}\eta\right), \\
\partial_z  \sigma & =
\tfrac{1}{d}.
\end{align}
\end{subequations}
%
All of these coefficients can be computed explicitly from the known two-dimensional free surface and bottom positions at given instants of time.

The velocity field can be determined from a known $\Phi$ using the relation
\begin{align}
({\bf u},w) = (\boldsymbol{\nabla}, \partial_z \sigma \partial_\sigma) \Phi.
\end{align}
The flow can be computed from the scalar velocity potential and used for estimating non-hydrostatic pressure and resulting wave loads. An exact expression for local pressure as a function of the vertical coordinate can be found by vertical integration of the vertical momentum equation to be of the form
\begin{align}
p(z) = \rho g (\eta -z ) + \int_{z}^\eta \partial_t w dz + \frac{1}{2}(\tilde{u}^2-u(z)^2 + \tilde{v}^2 - v(z)^2 + \tilde{w}^2 - w(z)^2).
\end{align}
The integral term can be estimated numerically using an accurate quadrature rule. Estimation of structural forces is determined by integration of pressure
\begin{align}
{\bf F} =  -\iint_S p{\bf n} dS,
\label{ch7:forcecalc}
\end{align}
where $S$ is a structural surface. 

\subsection{Finite difference approximations}\index{finite difference}

The numerical scheme is implemented as a flexible-order finite difference collocation scheme where all finite difference approximations of derivatives are constructed from one-dimensional approximations in a standard way each having the maximum possible accuracy. In explicit numerical schemes, finite difference approximations can be implemented using a matrix-free technique to exploit that only a few different stencils are in fact needed. This can significantly reduce memory requirements of the implemented model by exploiting that the same small set of stencils can be reused. See chapter \ref{ch5} for more details about matrix free stencil operations supported in our in-house library for heterogenous and massively parallel computing using GPUs.

%\newpage
\subsection{Time integration}\index{time integration}

For users of scientific applications robustness is of paramount importance for the solution of time-dependent PDEs. This makes stability considerations relevant in the context of both explicit and iterative numerical methods often considered most suitable for massively parallel applications. In the following, we address aspects of explicit time integration schemes which is associated with a stability\index{stability} requirement on time steps.

A Method of Lines\index{method of lines} approach is used for the discretization of the wave model. The spatial discretization yields a system of ordinary differential equations which can be expressed as a semi-discrete system.
%
We use the classical fourth order Explicit Runge-Kutta Method (ERK4). This algorithm is suitable for massive parallel computations via a data-parallel implementation where the spatial discretization terms are processed. As a means to introduce more concurrency into the time integration we consider for the first time the 'Parareal' algorithm as described in section \ref{ch7:parareal}.
For explicit time-integration schemes a Courant-Levy-Friedrichs (CFL) condition defines a necessary restriction for temporal stability of the form
\begin{align}
\Delta t\leq \frac{C}{\max_{n} |\lambda_n(\mathcal{J}_h)|},
\end{align}
with $C\in\mathbb{R}_+$ a CFL constant typically of size $\mathcal{O}(1)$ dependent on chosen scheme, and $\mathcal{J}_h\in\mathbb{R}^{2m\times2m}$, where $m=N_xN_y$, is a discrete Jacobian\index{Jacobian} matrix obtained by local linearization in time of \eqref{ch7:FSorigin}. For ERK4, $C=2\sqrt{2}$ if all eigenvalues are purely imaginary.

To gain insight into necessary conditions for stability, we employ a linear stability analysis based on the semi-discrete linear system
%
\begin{align}
\label{ch7:linstab}
\frac{d}{dt}\left[ \begin{array}{c} \eta \\ \tilde{\phi} \end{array} \right] = \mathcal{J} \left[ \begin{array}{c} \eta \\ \tilde{\phi} \end{array} \right], \quad \mathcal{J}=\left[ \begin{array}{cc} 0 & \mathcal{J}_{12} \\ -g & 0 \end{array} \right],
\end{align}
%
where $\mathcal{J}_{12} = k \tanh(kh)$ for the continuous problem which results in purely imaginary eigenvalues of the Jacobian $\lambda(\mathcal{J}) = \pm i \sqrt{gk\tanh(kh)}$ that grow unbounded in the limit $kh\to\infty$ corresponding to deep water relative to wave length. The waves are described in terms of wave number $k=2\pi/L$ where $L$ is the wave length.

For any numerical method, a numerical discretization of the same system of equations should mimic the continuous eigenspectrum in well-resolved parts of the spectrum. For a given discretization method, the stability of the discrete problem \eqref{ch7:linstab} is governed by the eigenvalues of the discrete Neumann-Dirichlet operator $\mathcal{J}_{12,h}\in\mathbb{R}^{m\times m}$, which connects the vertical velocity at the free surface with that of the velocity potential at the free surface. It is possible to partition the discrete $\sigma$-transformed Laplace problem and the differential in the $z$-direction in terms of equations for the interior (i) and those at the free surface boundary (b) such that
%
\begin{align}
\mathcal{A}\phi = \left[ \begin{array}{cc} \mathcal{A}_{bi} & \mathcal{A}_{bb} \\  \mathcal{A}_{ii} & \mathcal{A}_{ib} \end{array}  \right] \left[ \begin{array}{c}  \tilde{\phi} \\  \phi_{i} \end{array}  \right], \quad
\mathcal{D} \phi = \left[ \begin{array}{cc} \mathcal{D}_{bi} & \mathcal{D}_{bb} \\  \mathcal{D}_{ii} & \mathcal{D}_{ib} \end{array}  \right] \left[ \begin{array}{c}  \tilde{\phi} \\  \phi_{i} \end{array}  \right].
\end{align}
%
From these equations we find for the discrete block Jacobian operator
%
\begin{align}
\frac{\partial \phi}{\partial z}\Big|_{z=\eta} = \mathcal{J}_{12,h} \tilde{\phi}, \quad  \mathcal{J}_{12,h} = \mathcal{D}_{bb} - \mathcal{D}_{bi} \mathcal{A}_{ii}^{-1} \mathcal{A}_{ib}.
\end{align}
%
As shown in \cite{ch7:RobertsonSherwin1999} the eigenvalues of $\mathcal{J}_{12,h}$ is related to the eigenvalues of the discrete Jacobian $\mathcal{J}_h$ through the following relationship
%
\begin{align}
\lambda(\mathcal{J}_h) = \pm i \sqrt{ \lambda(\mathcal{J}_{12,h}) g }.
\end{align}
%
and is all imaginary confirming the hyperbolic (energy-conserving) nature of the potential flow formulation. Thus, for a given discretization of the linearized equations, it is possible to compute the eigenvalues of the discrete block operator to determine the eigenspectrum of the full operator. A discrete analysis of the eigenvalues is given in \cite[Section 4.1]{ch7:EBL08}, but it is not clearly pointed out that in fact the discrete eigenspectrum is compact (bounded) for a fixed polynomial order in the vertical, i.e. that for a constant depth $h$
\begin{align}
\max|\lambda(\mathcal{J}_h)| = \lim_{kh\to\infty}|\lambda(\mathcal{J}_h)|\leq C(N_z)\sqrt{\frac{g}{h}}.
\end{align}
Similar results were reported for the first time in the context of high-order Boussinesq-type equations in \cite{ch7:ENG06,ch7:EHBM06}. This is an important practical property of the discrete scheme as it is favourable to numerical stability. It implies that the linear model is not severely limited by the spatial resolution in the horizontal for a specific choice of the number  of collocation nodes ($N_z$) in the vertical. This suggests that the model is quite robust due to insensitivity in the choice of time step, with the implication that local grid adaptivity can be used for improving spatial accuracy. Interestingly, for the unified potential flow model we find that this also holds for nonlinear simulations. Large time steps can be chosen when using dense grids and high-order numerics without severely degrading overall numerical stability and efficiency. This is confirmed in numerical experiments and demonstrated in figure \ref{ch7:numexp}. However, for very steep nonlinear waves and very densely clustered non-uniform grids, stability is found to be compromised without filtering. A proper filtering strategy can be used to remedy stability problems without destroying accuracy. 
%
\begin{figure}[!htb]
\centering
\subfigure[Grid scaling, $x=(1-a)\xi^3+a\xi$.]{
% MainLaplace2D_ex03.m
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/scalingNx25-eps-converted-to.pdf}
}
\subfigure[High-order spatial discretisation and stable explicit time-stepping with large time steps for a nonlinear standing wave. Scaling based on $a=0$. ]{
% MainLaplace2D_ex03.m
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/standingwaveglozman-eps-converted-to.pdf}
}
\subfigure[Uniform grid ($a=1$).]{
% MainLaplace2D_ex035_nonlinearLaplace.m
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/SFwaves_snapshots_uniform-eps-converted-to.pdf}
}
\subfigure[Clustered grid ($a=0.05$).]{
% MainLaplace2D_ex035_nonlinearLaplace.m
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/SFwaves_snapshots_clustered-eps-converted-to.pdf}
}
\caption[Numerical experiments to assess stability properties of numerical wave model.]{Numerical experiments to assess stability properties of numerical wave model. In three cases, computed snapshots are taken of the wave elevation over one wave period of time. In a) the grid distribution of nodes in a one-parameter mapping for the grid is illustrated. Results from changes in wave elevation are illustrated for b) a mildly nonlinear standing wave on a highly clustered grid, c) regular stream function wave of medium steepness in shallow water $(kh,H/L)=(0.5,0.0292)$ on a uniform grid ($N_x=80$) and d) for a nonuniform grid with a minimal grid spacing 20 times smaller(!). In every case the step size remains fixed at $\Delta t = T/160$ s corresponding to a Courant number $C_r=c\tfrac{\Delta t}{\Delta x}=0.5$ for the uniform grid. A 6'$th$ order scheme and explicit EKR4 time-stepping is used in each test case.}
\label{ch7:numexp}
\end{figure}
%\newpage

\subsection{Wave generation and absorption}
\label{ch7:wavegen}

To simulate waves using a numerical model, a general purpose technique for both generating and absorbing waves inside the finite numerical domain is needed. It is preferable that the technique is suitable for easy integration in a software library component setup. One such technique is the line relaxation\index{relaxation} method attributed to \cite{ch7:LD83}. This is a simple ad hoc technique sufficiently accurate for engineering purposes. It modifies the computed solution every time step during simulation by a post-processing step where the relaxed solution becomes
%determined as
\begin{align}
\label{ch7:eq:discreteupdate}
g^*(x_i,t) = \Gamma(x_i)g(x_i,t) + (1-\Gamma(x_i))g_e(x_i,t), \quad x_i\in\Omega_\Gamma.
\end{align}
Here $g(x,t)$ is one of the free surface variables $\tilde{\phi},\eta$ at an instant in time and $0\leq \Gamma(x)\leq 1$ is a single-valued function within the relaxation region $x_i\in\Omega_\Gamma$. The first term acts as a sponge layer which is responsible for effectively dissipating energy inside a specified relaxation zone\index{relaxation!zone} $\Omega_\Gamma$. The terms containing $g_e(x,t)$, where $g_e$ is an analytical solution (e.g. such as stream function wave theory \cite{ch7:Dean1965}) act as source terms in the relaxation zone. This makes it possible to generate arbitrary waves accurately in the computational domain in accordance with an analytical representation of incident waves.

We can interpret \eqref{ch7:eq:discreteupdate} as a discrete update of the solution at an isolated spatial point inside a relaxation zone. We introduce the notation $g^n=g(x_i,t_n)$ and $g^{*,n+1}=g^*(x_i,t_{n+1})$ and assume that $t_{n+1}=t_n+\tau$ is an instant in time. Then we can rewrite \eqref{ch7:eq:discreteupdate} to motivate an analytical modification to time-dependent equations that can provide similar modification (forcing) in simulation.

Subtract $g^n$ in \eqref{ch7:eq:discreteupdate} and divide by a pseudo time step size $\tau$ to obtain the equivalent form

%%
\begin{align}
\frac{g^{*,n+1}-g^n}{\tau} =\frac{(1-\Gamma)}{\tau} (g_e^n-g^n).
\end{align}
%
The first term is similar to a first order accurate Forward Euler\index{forward Euler} approximation of a rate of change term. This motivates an {\em embedded penalty forcing technique} based on adding a correction term of the form
%
\begin{align}\label{ch7:eq:penalty}
\partial_t g = \mathcal{N}(g) + \frac{1-\Gamma(x)}{\tau} (g_e(t,x)-g(t,x)), \quad {\bf x}\in\Omega_\Gamma,
\end{align}
%
where $\mathcal{N}$ represents a general nonlinear operator for the right hand side. The immediate advantage is that a time stepping scheme can easily be interchanged in a model implementation. The added terms is a source term resulting in forcing inside relaxation zones when $g_e(t,x)\neq g(t,x)$ and $\Gamma(x)\neq1$ and otherwise has no effect. The strength of the forcing is influenced by the arbitrary parameter $\tau\in\mathbb{R}_+$ which can be defined to match the time scale of the dynamics. We have found that  $\tau\approx\Delta t$ works well, however, it is possible that a more optimal choice exist. Remark, a too small $\tau$ may degrade numerical stability of the model.

A simple validation of the zones is shown in figure \ref{ch7:figstandwave} where waves are generated at the left wall, and propagate to the right wall, where reflection occur leading formation of standing waves du to resulting interaction with the incident waves inside the numerical wave tank.

The following relaxation functions proposed in \cite{ch7:ENG06} guarantees continuity across interface of relaxation zone and computational domain and are used in simulations for, respectively, sponge layers and wave generation zones
%
\begin{align}
% \label{relaxfunc1}
\Gamma_{s}(x)  = 1-(1-x)^p, \quad \Gamma_{g}(x) = -2x^3+3x^2, \quad x\in[0,1].
\label{ch7:relaxfunc2}
\end{align}
%
The profiles can be reversed by a change of coordinate, i.e. $\Gamma(1-x)$, and scaled to interval sizes of interest. The first function satisfies the condition that any derivative at the left boundary vanishes at $x=1$. The first derivatives of the second function vanish at both ends. The relaxation zones are positioned appropriately where waves are to be both/either generated and/or absorbed. The rule of thumb is that a relaxation used for absorption has a spatial length of at least two wave lengths. For absorption zones, we find that this technique is more efficient in velocity formulation of the free surface equations often used in Boussinesq-type formulations, e.g., see \cite{ch7:ENG06,ch7:EHBM06,ch7:EHBW08} in comparison with scalar potential formulations \eqref{ch7:FSorigin}. However, similar performance can be achieved by merely increasing the length of relaxation zones in such regions. Demonstrations of the technique are seen in figure \ref{ch7:figstandwave} where vertical dashed lines indicates interfaces between relaxation zones and the computational region. Incident waves propagate from left to right in both examples.
%
\begin{figure}[!htb]
\centering
\subfigure[Wave generation, reflection and absorption of small-amplitude waves.]{
% Script : MainLaplace2D_ex03penalityLINEAR_REFLECTEDWAVES.m
\includegraphics[width=0.98\textwidth]{Chapters/chapter7/figures/standingwavespenalty-eps-converted-to.pdf}
% Nx = 480, 6th order, vertical clustering, Nz=6;
}
\subfigure[Wave generation and absorption of steep finite-amplitude waves.]{
% Script : MainLaplace2D_ex03penalityNONLINEAR_GENERATEWAVES.m
\includegraphics[width=0.98\textwidth]{Chapters/chapter7/figures/nonlinearwavespenalty-eps-converted-to.pdf}
% Nx = 540, 6th order, vertical clustering, Nz=6;
}
\caption[Snapshots at intervals $T/8$ over one wave period in time.]{Snapshots at intervals $T/8$ over one wave period in time of computed a) small-amplitude $(kh,kH)=(0.63,0.005)$ and b) finite-amplitude $(kh,kH)=(1,0.41)$ stream function waves elevations having reached a steady state after transient startup. Combined wave generation and absorption zones in the western relaxation zone of both a) and b). In b) an absorption zone is positioned next to the eastern boundary and causes minor visible reflections. }
\label{ch7:figstandwave}
\end{figure}

\subsection{Preconditioned Defect Correction (PDC) method}\index{preconditioning}\index{defect correction}
\label{ch7:PDCmethod}

For the solution of sparse linear systems
\begin{align}
\label{ch7:linsyscompact}
\mathcal{A}{\bf \Phi} = {\bf b}, \quad \mathcal{A}\in\mathbb{R}^{n\times n}, \quad {\bf b}\in\mathbb{R}^{n},
\end{align}
it is attractive to use iterative methods for large system sizes $n=N_xN_yN_z$ and for parallel implementations. Acceleration of suitable iterative methods can be done, e.g., by instead solving a left-preconditioned\index{preconditioning!left-preconditioning} system of the form %typical %A key challenge is preconditioning for acceleration of convergence of iterative methods.
\begin{align}
\label{ch7:eq:linsys}
\mathcal{M}^{-1} ( \mathcal{A} {\bf \Phi} = {\bf b}), \quad \mathcal{M}\in\mathbb{R}^{n \times n},
\end{align}
where $\mathcal{M}$ is a preconditioned with the property that $\mathcal{M}^{-1}\approx \mathcal{A}^{-1}$ can be computed at low cost.

The bottleneck problem in a unified potential flow model is the solution of a discrete $\sigma$-transformed Laplace problem stated in the compact forms \eqref{ch7:linsyscompact} or \eqref{ch7:eq:linsys}. It is attractive to find an efficient iterative strategy where convergence\index{convergence} is understood via a convergence theory, have modest storage requirements, have minimal global communication requirements (in the form of global inner products) and is suitable for flexible-order\index{flexible order} discretizations. The class of geometric Multigrid Methods\index{multigrid} fulfils these requirements and have shown to be among the most efficient iterative strategies for a wide class of problems \cite{ch7:Trottenberg01}. In particular, the time required to solve a system of linear equations to a given accuracy level can be made to scale proportional to the number of unknowns.

There are several known approaches to multigrid methods \cite{ch7:MR744926} for high-order discretizations. Among these, Defect Correction\index{defect correction} Methods (DCMs) \cite{ch7:Stetter1978,ch7:AuzingerStetter1982,ch7:Hackbusch1982,ch7:Auzinger1987,ch7:Trottenberg01} have been employed successfully, e.g., in Computational Fluid Dynamics \cite{ch7:LaytonEtAl2002}, in numerical simulations since the early 1970s.  
The fundamental idea of DCMs is to combine the good stability properties of low-order discretizations with high-order accuracy discretizations for explicit residual evaluations. These iterative methods impose low storage requirements and have scalable work effort under suitable choices of preconditioning strategies and may be accelerated using a multigrid method based on low-order discretizations, while still achieving high-order accuracy.
Furthermore, it has been shown, that the rate of convergence of DCM combined with standard multigrid methods can achieve rates of convergence corresponding to the most efficient multigrid methods~\cite{ch7:Auzinger1987}.

Therefore, for the efficient and scalable solution of the unified potential flow model, we have recently \cite{ch7:EngsigKarupEtAl2011} proposed a Preconditioned Defect Correction (PDC) Method for efficient iterative low-storage solution of high-order accurate discretization of the $\sigma$-transformed Laplace problem \eqref{ch7:TransformedLaplace}. The proposed strategy can be seen as a generalization of the multigrid strategy proposed by \cite{ch7:LiFleming1997}. The PDC method enables significant improvement of overall efficiency and accuracy with the preconditioning based on a second-order linearized version of the full coefficient matrix $\mathcal{A}$ as described in \cite{ch7:EBL08}.

 Starting from some initial guess ${\bf \Phi}^{[0]}\in\mathbb{R}^n$, the PDC method for solving \eqref{ch7:eq:linsys} can be stated compactly as a three-step recurrence for $k=1,2,...$
\begin{align}
\label{ch7:eq:defectcorrectionprocedure}
\Phi^{[k]} = \Phi^{[k-1]} + \delta^{[k-1]}, \;\; \delta^{[k-1]}=\mathcal{M}^{-1} {\bf r}^{[k-1]}, \;\; {\bf r}^{[k-1]}={\bf  b} - \mathcal{A}\Phi^{[k-1]},
\end{align}
where $\Phi^{[k]},{\boldsymbol \delta}^{[k]},{\bf r}^{[k]}\in\mathbb{R}^n$ are, respectively, the approximate solution, the defect (preconditioned residual) and the residual of \eqref{ch7:eq:linsys} at the $k$'th iteration. The algorithm can be speedup by using mixed-precision calculations on modern many-core hardware as demonstrated in \cite{ch7:Glimberg2011}.  

\subsection{Distributed computations via domain decomposition}\label{ch7:sec:dd}\index{domain decomposition}

Numerical modelling of large ocean areas to account for nonlinear wave-wave interactions and wave-structure interactions require large degrees of spatial resolution, significant computational resources and parallel computations to be practical. The recent generations of programmable GPUs are heavily optimized for on-chip bandwidth performance but not capacity. This implies that for the solution of large-scale PDE problems, distributed computing on multiple GPU devices is required due to limited capacity in the global memory space of current GPUs. Via a combination of MPI and CUDA we have recently demonstrated how both small and large systems can be solved efficiently by heterogenous computations using a data domain decomposition technique in parallel \cite{ch7:GlimbergEtAl2012}. The idea is to distribute the computational tasks to multiple GPUs, to enable reduced computational times and increased problem sizes. A homogenous partitioning of the data is used to ensure that the load balance across the GPUs is uniform. Data distribution and message passing introduce a data transfer bottleneck in the form of the PCIe link and network interconnection. Thus, if the computational intensity of the local problem is not large enough to enable sufficient latency hiding of this bottleneck, the whole application are likely to be (severely) limited by the PCIe link or network bandwidth performance rather than the high on-chip bandwidth of the individual GPUs.

The ratio between necessary data transfers and computational work for the proposed numerical model for free surface water waves is high enough to expect reasonable latency hiding. The data domain decomposition method consists of a logically structured division of the computational domain into multiple subdomains. Each of these subdomains are connected via fictitious ghost layers at the artificial boundaries of width corresponding to the half-width of the finite difference stencils employed. This results in a favourable volume-to-boundary ratio as the problem size increases, diminishing communication overhead for message passing. Information between subdomains are exchanged through ghost layers at every step of the iterative PDC method, in connection with the matrix-vector evaluation for the $\sigma$-transformed Laplace problem, and before relaxation steps in the multigrid method. A single global synchronization point occur at most once each iteration, if convergence is monitored, where a global reduction step (inner product) between all processor nodes takes place. The main advantage of this decomposition strategy is, that the decomposition into multiple subdomains is straightforward. However, it comes with the cost of extra data transfers to update the set of fictitious ghost layers.

\begin{figure}[!htb]
    \setlength\figureheight{0.30\textwidth}
    \setlength\figurewidth{0.33\textwidth}
    \begin{center}
    \subfigure[Test environment 3.]{
    {\scriptsize\input{Chapters/chapter7/figures/TeslaK20PerformanceScaling3D.tikz}}
    }
    \subfigure[Test environment 3.]{
    {\scriptsize\input{Chapters/chapter7/figures/TeslaK20SpeedupGPUvsCPU3D.tikz}}
    }
    \end{center}
    \caption[Performance timings per PDC iteration as a function of increasing problem size $N$, for single, mixed and double precision arithmetics.]{Performance timings per PDC iteration as a function of increasing problem size $N$, for single, mixed and double precision arithmetics. Three dimensional nonlinear waves, using $6^{th}$ order finite difference approximations, preconditioned with one multigrid V-cycle and one pre- and post- Red-black Gauss-Seidel smoothing. Speedup compared to fastest known serial implementation. Using Test environment 3. CPU timings represent starting point for our investigations and has been obtained using Fortran 90 code and is based on a single-core run on a Intel Core i7, 2.80GHz processor.}\label{ch7:fig:perftimings}
\end{figure}

The parallel domain decomposition solver has been validated against the sequential solvers with respect to algorithmic efficiency to establish that the code produce correct results. An analysis of the numerical efficiency have also been carried out on different GPU systems to identify comparative behaviors as both the problems sizes and number of compute nodes vary. For example, performance scalings on Test environment 1 and Test environment 2 are presented in figure \ref{ch7:fig:multigpuperformance}. The figure confirms that there is only a limited benefit from using multiple GPUs for small problem sizes, since the computational intensity is simply too low to efficiently hide the latency of message passing. A substantial speedup is achieved compared to the single GPU version, while being able to solve even larger systems.
With the linear scaling of memory requirements and improved computational speed, the methodology based on multiple GPUs makes it possible to simulate water waves in very large numerical wave tanks with improved performance.





\begin{figure}[!htb]
    \setlength\figureheight{0.4\textwidth}
    \setlength\figurewidth{0.68\textwidth}
    \begin{center}
    \subfigure[Test environment 1.]{
    {\scriptsize\input{Chapters/chapter7/figures/GTX590MultiGPUScaling3D.tikz}}
    }
    \subfigure[Test environment 2.]{
    {\scriptsize\input{Chapters/chapter7/figures/TeslaM2050MultiGPUScaling3D.tikz}}
    }
    \end{center}
    \caption[Domain decomposition performance on multi-GPU systems.]{Domain decomposition performance on multi-GPU systems. Performance timings per PDC iteration as a function of increasing problem sizes using single precision. Same setup as in figure \ref{ch7:fig:perftimings}.}
    \label{ch7:fig:multigpuperformance}
\end{figure}

Future work involves extending the domain decomposition method to include support for more general curvilinear boundary-fitted meshes for representing the free surface plane and include bottom-mounted structures which is typically encountered in near-costal areas. This will enable opportunities to resolve wave-structure interactions such as those encountered in large harbour regions and isolated human-made structures in offshore regions, e.g., legs of oil rigs or offshore windmill foundations.

\subsection{Assembling the wave model from library components}

It is described in chapter \ref{ch5} how we have developed a heterogenous library has for fast prototyping of PDE solvers, utilizing the massively parallel architecture of CUDA-enabled GPUs. The objective is to provide a set of generic components within a single framework, such that software developers can assemble application specific solvers efficiently at a high abstraction level, requiring a minimum of CUDA specific kernel implementations and parameter tuning.

The CUDA-based numerical wave model has been developed based on all the numerical techniques described in preceding sections. These techniques are a part of the library implemented as generic components, which makes them useful for the numerical solutions of partial differential equation (PDE) systems. The components of the numeral model as described in section \ref{ch7:sec:nummodel} is an ERK4 time integrator, flexible-order finite difference approximations on regular grids, and an iterative multigrid PDC solver for the $\sigma$-transformed Laplace equation \eqref{ch7:TransformedLaplace}. Application developers can either use these components directly from the library or they are free to combine their own implementations with library components, if they need alternative strategies that are not present in the library.

For the unified potential flow model the user will need to provide implementations of the following components; the right hand side operator for the semi-discrete free surface variables \eqref{ch7:FSorigin}, the matrix-vector operator for the discretized $\sigma$-transformed Laplace equation \eqref{ch7:TransformedLaplace}, a smoother for the multigrid relaxation step, and the potential flow solver itself, that reads initial data and advance the solution in time. In order to make the library as generic as possible, all components are template-based, which makes it possible to assemble the PDE solver by combining type definitions in the preamble of the application. An excerpt of the potential flow assembling is given in listing \ref{ch7:lst:solversetup}.

\lstset{label=ch7:lst:solversetup,caption={Generic assembling of the potential flow solver for fully nonlinear free surface water waves.}
%,basicstyle=\scriptsize
}
\begin{lstlisting}
// Basics
typedef double value_type;
typedef gpulab::grid<value_type> vector_type;
typedef free_surface::laplace_sigma_stencil_3d<vector_type> matrix_type;

// Multigrid setup
typedef gpulab::solvers::multigrid_types<
    vector_type
  , matrix_type
  , free_surface::gs_rb_low_order_3d
  , gpulab::solvers::grid_handler_3d_boundary> multigrid_types;
typedef gpulab::solvers::multigrid<multigrid_types> preconditioner_type;

// Defect Correction setup
typedef gpulab::solvers::defect_correction_types<
    vector_type
  , matrix_type
  , monitor_type
  , preconditioner_type> dc_types;
typedef gpulab::solvers::defect_correction<dc_types> laplace_solver_type;

// Potential flow solver setup
typedef free_surface::potential_flow_solver_types<
    vector_type
  , laplace_solver_type
  , gpulab::integration::ERK4> potential_flow_types;
typedef free_surface::potential_flow_solver_3d<potential_flow_types> potential_flow_solver_type;
\end{lstlisting}

Hereafter, the potential flow solver is aware of all component types that should be used to solve the entire PDE system, and it will be easy for developers to exchange parts at later times. The \texttt{laplace\_sigma\_stencil\_3d} class implements both the matrix-vector and right hand side operator. The flexible-order finite difference kernel for the matrix-free matrix-vector product for the two-dimensional Laplace problem is presented in listing \ref{ch7:lst:fd2d}. Library macros and reusable kernel routines are used throughout the implementations to enhance developer productivity and hide hardware specific details. This kernel can be used both for matrix-vector products for the original system and for the preconditioning.

\lstset{label=ch7:lst:fd2d,caption={CUDA kernel implementation for the two dimensional finite difference approximation to the transformed Laplace equation.}
%,basicstyle=\scriptsize\ttfamily
}
\begin{lstlisting}
template <typename value_type, typename size_type>
__global__ void laplace_sigma_transformed(
      value_type* out            , value_type const* p
    , value_type const* eta      , value_type const* etax
    , value_type const* etaxx    , value_type const* h
    , value_type const* hx       , value_type const* hxx
    , value_type dx              , value_type ds
    , size_type Nx               , size_type Ns
    , size_type xstart           , size_type alpha
    , value_type const* stencil1 , value_type const* stencil2)
{
        size_type j    = IDX1Dx;
        size_type i    = IDX1Dy;
        size_type rank = alpha*2+1; // Total stencil size

        // Get shared memory pointers
        value_type* ss1 = device::shared_memory<value_type>::get_pointer();
        value_type* ss2 = ss1 + rank*rank;
        
        // Load stencils into ss1 and ss2
        gpulab::device::load_shared_memory(ss1,stencil1,rank*rank);
        gpulab::device::load_shared_memory(ss2,stencil2,rank*rank);
        
        // Only internal points
        if(j>=xstart && j<Nx-xstart && i>0 && i<Ns-1)
        {                       
                size_type offset_i = i < alpha ? 2*alpha-i : i >= Ns-alpha ? Ns-1-i : alpha;
                size_type row_i    = offset_i*rank;
    // Always centered stencils in x-dir
                size_type offset_j = alpha;  
                size_type row_j    = alpha*rank;
                        
                value_type dhdx    = hx[j];
                value_type dhdxx   = hxx[j];
                value_type detadx  = etax[j];
                value_type detadxx = etaxx[j];

                // Calculate the sigma transformation variables
                value_type sig   = kernel::sigma(i,ds);
                value_type d     = h[j]+eta[j];
                value_type dsdz  = kernel::dsdz(d);
                value_type dsdx  = kernel::dsdx (d,sig,detadx,dhdx);
                value_type dsdxx = kernel::dsdxx(d,sig,detadx,detadxx,dhdx,dhdxx,dsdx);

                // dp/dxx
                value_type dpdxx = values<value_type>::zero();
                for(size_type s = 0; s<rank; ++s)
                {
                        dpdxx += ss2[row_j+s]*p[i*Nx+j-alpha+s];
                }
                dpdxx /= (dx*dx);

                // Calculate dp/dss and dp/ds
                value_type dpdss = values<value_type>::zero();
                value_type dpds  = values<value_type>::zero();
                value_type ps;
                for(size_type s = 0; s<rank; ++s)
                {
                        ps = p[(i+offset_i-s)*Nx+j];
                        dpds  += ss1[row_i+s]*ps;
                        dpdss += ss2[row_i+s]*ps;
                }
                dpds  /= ds;
                dpdss /= (ds*ds);

                // Calculate dpdxds
                value_type dpdxds = values<value_type>::zero();
                for(size_type ss = 0; ss<rank; ++ss)
                {
                        value_type dpdx   = values<value_type>::zero();
                        for(size_type sx = 0; sx<rank; ++sx)
                        {
                                dpdx += ss1[row_j+sx]*p[(i+offset_i-ss)*Nx+j-offset_j+sx];
                        }
                        dpdx /= dx;
                        dpdxds += ss1[row_i+ss]*dpdx;
                }
                dpdxds /= ds;

                // Calculate total
                out[i*Nx+j] = dpdxx + dsdxx*dpds + 2.0*dsdx*dpdxds + (dsdx*dsdx + dsdz*dsdz)*dpdss;
        }
}
\end{lstlisting}

In a similar template-based approach, the kernel for the right hand side operator of the two dimensional problem is implemented and listed in Listing \ref{ch7:lst:rhs2d}. The kernel computes the right hand side updates for both surface variables, $\eta$ and $\tilde{\phi}$, and applies an embedded penalty forcing \eqref{ch7:eq:penalty}, for all nodes within generation or absorption zones. The penalty forcing functions are computed based on linear or non-linear wave theory in a separate device function.

\lstset{label=ch7:lst:rhs2d,caption={CUDA kernel implementation for the 2D right hand side.}%,basicstyle=\scriptsize\ttfamily
}
\begin{lstlisting}
template <typename value_type, typename size_type>
__global__ void rhs(value_type const* p    , value_type const* p_surf
                  , value_type const* eta  , value_type* dp_surf_dt
                  , value_type* deta_dt    , value_type const* etax
                  , value_type const* h    , value_type dx
                  , value_type ds          , size_type Nx
                  , size_type xstart       , size_type xend
                  , unsigned int alpha     , value_type const* stencil
                  , value_type t           , value_type dt
                  , bool generate          , bool absorb)
{
        size_type j    = IDX1D;
        size_type rank = alpha*2+1; // Total stencil size

    // Load shared memory for stencil into ss
        value_type* ss = gpulab::device::shared_memory<value_type>::get_pointer();
    gpulab::device::load_shared_memory(ss,stencil,rank*rank);

    // Only interior points
        if(j>=xstart && j<Nx-xstart)
        {
                value_type g     = values<value_type>::gravity();
        value_type detax = etax[j];

                // Calculate the sigma transformation variables
                value_type d     = h[j]+eta[j];
                value_type dsdz  = kernel::dsdz(d);

                // Calculate dp/ds on surface
                value_type dpds  = values<value_type>::zero();
                for(size_type s=0; s<rank; ++s)
                {
                        dpds += p[j+s*Nx]*ss[rank*rank-1-s];
                }
                dpds /= ds;
                value_type w     = dsdz * dpds;

                // Calculate dp/dx on surface
                value_type dpdx  = values<value_type>::zero();
                for(size_type s=0; s<rank; ++s)
                {
                        dpdx += p_surf[j-alpha+s]*ss[rank*alpha+s];
                }
                dpdx /= dx;

        // Update surface variables rhs
                dp_surf_dt[j] = - g*eta[j] - 0.5*(dpdx*dpdx - w*w*(1.0 + (detax*detax)));
                deta_dt[j]    = - detax*dpdx + w*(1.0 + (detax*detax));

                if(generate || absorb)
                {
                        // Relaxation terms
                        value_type reta  = values<value_type>::zero();
                        value_type rp    = values<value_type>::zero();
            free_surface::kernel::rhs_relax(reta, rp, eta[j], p_surf[j], (((int)j)-(int)xstart)*dx, (xend-xstart-1)*dx, t+dt, dt, false, generate, absorb);
                        
                        deta_dt[j]    = deta_dt[j]    + reta;
                        dp_surf_dt[j] = dp_surf_dt[j] + rp;
                }
        }
}
\end{lstlisting}

\section{Properties of the numerical model}

We now consider different properties of the numerical model in order to shed light on unique features and limits of the model with respect to maritime engineering applications. The presented results extend and complements earlier studies \cite{ch7:BinghamZhang2007,ch7:EBL08} for the same model. In particular, we seek to highlight that the properties are tuneable to the practical applications of interest through proper choice of discretization parameters and therefore also provide details of numerical properties.

\subsection{Dispersion and kinematic properties}\index{dispersion}\index{kinematic}
\label{ch7:dispkin}

The dispersion and kinematic properties of the unified model is determined by the tuneable discretization parameters and should in general be chosen for specific problems. For assessment of errors, we introduce the metrics proposed in \cite{ch7:MBS03}
\begin{subequations}
\begin{align}
E_c(N_x,N_z,kh,h/L) & = \frac{c^2}{c_s^2},
\label{ch7:errdisp} \\ % \left( \frac{c-c_s}{c_s}  \right)^2, \\
\label{ch7:errkin}
E_m(N_x,N_z,kh,h/L) &= \frac{1}{h}\int_{-h}^\eta \left( \frac{\phi(z)-\phi_s(z)}{\phi_s(z)} \right)^2 dz,
\end{align}
\end{subequations}
where $m$ is one of the scalar functions $\phi,u,w$ describing kinematics, $c$ is the numerical phase celerity of regular waves and $c_s=\sqrt{g\tanh(kh)/k}$ is the exact phase celerity according to linear Stokes Theory \cite{ch7:SJ01}.  Measurements of the error are taken in the vertical below the crest of a wave which is well-resolved in the horizontal direction.

% kinematicAnalysis.m
\begin{figure}[!htb]
\begin{center}
\subfigure[Uniform vertical grid.]{
%\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/lineardispersion_Nx30-HL90-p6-vergrid0_Linear.eps}
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/lineardispersion_Nx30-HL90-p6-vergrid0_Linear-eps-converted-to.pdf}
}
\subfigure[Cosine-clustered vertical grid.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/lineardispersion_Nx30-HL90-p6_Linear-eps-converted-to.pdf}
}
\end{center}
\caption{The accuracy in phase celerity $c$ determined by \eqref{ch7:errdisp} for small-amplitude (linear) wave.
$N_z\in[6,12]$. Sixth order scheme.}
\label{ch7:figlinear}
\end{figure}

% kinematicAnalysis.m
\begin{figure}[!htb]
\begin{center}
\subfigure[Linear]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/kinematicsPHI_Nx30-HL90-p6_Linear-eps-converted-to.pdf}
}
\subfigure[Linear]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/kinematicsW_Nx30-HL90-p6_Linear-eps-converted-to.pdf}
}
\subfigure[Nonlinear]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/kinematicsPHI_Nx30-HL90-p6_Nonlinear-eps-converted-to.pdf}
}
\subfigure[Nonlinear]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/kinematicsW_Nx30-HL90-p6_Nonlinear-eps-converted-to.pdf}
}
\end{center}
\caption[Assessment of kinematic error is presented in terms of the depth-averaged error.]{Assessment of kinematic error is presented in terms of the depth-averaged error determined by \eqref{ch7:errkin} for a) scalar velocity potential and b) vertical velocity for a small-amplitude (linear) wave, and c) scalar velocity potential and d) vertical velocity for a finite-amplitude (nonlinear) wave with wave height $H/L=90\%(H/L)_\textrm{max}$.
$N_z\in[6,12]$. Sixth order scheme. Clustered vertical grid. }
\label{ch7:figlinear2}
\end{figure}

The accuracy of the dispersion and kinematic properties are found to be excellent with small differences between the accuracy of computed profiles for  linear and nonlinear waves. Figure \ref{ch7:figlinear} shows curves from a discrete analysis of how the accuracy can be improved by increasing the number of nodes in the vertical for a) uniform grids and b) cosine-clustered grids. Similarly, figure \ref{ch7:figlinear2} shows how the accuracy of the kinematic properties of the model can be controlled by choosing an appropriate number of co-sine clustered vertical collocation points.

%\section{Verification and Validation}
\section{Numerical experiments}

The numerical model detailed has been subject to careful verification and validation utilizing a range of standard benchmarks, cf. \cite{ch7:EBL08,ch7:EngsigKarupEtAl2011}. Here we exclusively focus on properties and performance of the numerical wave model. We provide several new results that highlights possibilities for acceleration of the wave model via simple and readily applicable techniques that works well on massively parallel hardware. Finally, we describe how we can extend the implementation of the wave model into a novel GPU implementation of a linear ship-wave model for fast hydrodynamics calculations.

\subsection{On precision requirements in engineering applications}\index{precision}

Practical engineering applications are widely used for analysis purposes and give support to decision making in engineering design. For engineering purposes the turn-around time for producing analysis results is of crucial importance as it affects cost-benefit of work efforts. The key interest is often just 'engineering accuracy' in computed end results which suggest that we can do with less precision in calculations. One may ask what is the precision requirements for engineering applications? In a recent study \cite{ch7:EngsigKarupEtAl2011,ch7:Glimberg2011}, it is shown that the PDC method when executed on GPUs can be utilized to efficiently solve water wave problems. This was done by trading accuracy for speed in parts of the PDC algorithm, e.g. by using either single or mixed-precision\index{precision!mixed} computations. Without preconditioning the PDC method reduces to a classical iterative refinement technique, which is known to be fault tolerant \cite{ch7:Higham:2002:ASN}.

Previously reported performance results for the wave model can be taken a step further. We seek to demonstrate how single precision computations can be used for engineering analysis without significantly affecting accuracy in final computational results. At the same time improvements in computational speed can be as much as a factor of two for large problems as a direct result of reduced data transfer, cf. figure \ref{ch7:fig:perftimings}. Therefore, in pursue of high performance, it is of interest to exploit the reduced data transfers associated with replacing double precision with single precision floating point calculations. In a well organized code this step can be taken with minimal programming effort.
%
\begin{figure}[!htb]
\begin{center}
% MainLaplace2D_ex025nonlinearLaplaceSINGLE.m
\subfigure[Single precision.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/PrecisionSINGLE-eps-converted-to.pdf}
}
\subfigure[Double precision.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/PrecisionDOUBLE-eps-converted-to.pdf}
}
\end{center}
\caption[Comparison between convergence histories for single and double precision computations using a PDC method for the solution of the transformed Laplace problem.]{Comparison between convergence histories for single and double precision computations using a PDC method for the solution of the transformed Laplace problem. Very steep nonlinear stream function wave in intermediate water $(kh,H/L)=(1,0.0903)$. Discretizaiton based on $(N_x,N_z)=(15,9)$ with 6'$th$ order stencils.}
\label{ch7:convhist}
\end{figure}

Most scientific applications \cite{ch7:lessismore} use double precision calculations to minimize accumulation of round-off errors, to employ higher precision for ill-conditioned problems, or to stabilize critical sections in the code that requires higher precision. The main reason is that roundoff-errors tend to accumulate slower when higher precision is used, thereby avoiding significant losses of accuracy due to round-off errors. Paradoxically, for many computational tasks, the need for high precision connected with the above-mentioned restrictions do not apply at all, or only apply to a portion of the tasks. In addition, on modern hardware there can be relative differences in peak floating point performance of up to about one order of magnitude in favor of single precision over double precision math calculations depending on choice of hardware architecture. However, for bandwidth bound applications, the key performance metric is not floating point performance, but rather bandwidth performance. In both cases, data transfer can be effectively halved by switching to single precision storage, in which case bandwidth performance increases and at the same time makes it possible to feed floating point units with data at effectively twice the speed. If maximizing performance is an ultimate goal, such considerations suggests that it can be possible to compute faster by using single over double precision arithmetics fi accuracy requirements can be fulfilled.
%

As a simple illustrative numerical experiment, we can consider the iterative solution of \eqref{ch7:eq:linsys} using the PDC method using, respectively, single and double precision math. As a simple test case, we consider the solution of periodic stream function waves in two spatial dimensions. Computed convergence histories are presented in figure \ref{ch7:convhist} where it is clear that the main difference is in the attainable accuracy level achievable before stagnation. In both cases, the attainable accuracy, defined in terms of the absolute error for the exact stream function solution to the governing equations, is associated with accuracy of approximately $10^{-4}$ (solid line). In single precision math, the algebraic and estimated errors measured in the 2-norm can reach the level of machine precision. These results suggest that single precision math is sufficient for calculating accurate solutions at the chosen spatial resolution. We find that the iterative solution of the $\sigma$-transformed Laplace problem by the PDC method does not immediately lead to significant accumulation of roundoff-errors and further investigations are warranted for unsteady computations.

Elaborating on this example, we examine the propagation of a regular stream function wave in time. We consider the errors in wave elevation as a function of time, respectively, with and without a filtering\index{filtering} strategy for single precision calculations in comparison with double precision calculations. With an objective to exert control on accumulation of round-off errors that appear as high-frequency noise, the idea is to employ an inexpensive stencil-based filtering strategy. For example, a central filter in one spatial dimension
\begin{align}
\label{ch7:filter}
\mathcal{F}u(x_i) = \sum_{n=-\alpha}^{\alpha} c_n u(x_{i+n}),
\end{align}
where $c_n\in\mathbb{R}$ are the stencil coefficients and $\alpha\in\mathbb{Z}_+$ is the stencil half-width . An active filter can, e.g., be based on employing a Savitzky-Golay smoothening filter \cite{ch7:PT90}, e.g. the mild 7-point SG(6,10) filter, and applying it after every 10th time step to each of the collocation nodes defining the free surface variables. The same procedure can be used for stabilization of nonlinear simulations to remove high-frequency 'saw-tooth' instabilities as shown in \cite{ch7:EBL08}. This filtering technique can also remove high-frequency noise resulting from roundoff-errors in computations that would otherwise potentially pollute the computational results and in the worst case leave them useless. The effect of this type of filtering on the numerical efficiency of the model is insignificant.

Results from numerical experiments are presented in figure \ref{ch7:filtering} and most of the errors can be attributed to phase errors resulting from difference in exact versus numerical phase speed. In numerical experiments, we find that while results computed in double precision are not significantly affected by accumulation of round-off errors, the single precision results are. In figures \ref{ch7:filtering} a) and b), a direct solver based on sparse Gaussian Elimination within Matlab\footnote{\url{http://www.mathworks.com}.} is used to solve the linear system every stage and a comparison is made between single and unfiltered double precision calculations. It is shown in figure \ref{ch7:filtering} a) that without a filter, the single precision calculations result in 'blow-up' after which the solver fails just before 50 wave periods of calculation time. However, in figure \ref{ch7:filtering} b) it is demonstrated that invoking a smoothening filter, cf. \eqref{ch7:filter}, stabilizes the accumulation of round-off errors and the calculations continue and achieves reduced accuracy compared to the computed double precision results. Thus, it is confirmed that such a filter can be used to control and suppress high-frequency oscillations that results from accumulation of round-off errors. In contrast, replacing the direct solver with an iterative PDC method using the GPU-accelerated wave model appears to be much more attractive upon inspection of figures \ref{ch7:filtering} c) and d). The single precision results are found to be stable with and {\em without} the filter-based strategy for this problem. The calculations shows that single precision math leads to slightly faster error accumulation for this choice of resolution, however, with only small differences in error level during long time integration. This highlights that fault-tolerance of the iterative PDC method which contributes to securing robustness of the calculations.

\begin{figure}[!htb]
\begin{center}
% DriverWavemodelDecomposition.m
\subfigure[Direct solve without filter.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/ComparisonLUNoFiltering-eps-converted-to.pdf}
}
\subfigure[Direct solve with filter.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/ComparisonLUWithFiltering-eps-converted-to.pdf}
} \\
\subfigure[Iterative PDC solve without filter.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/ComparisonDCNoFiltering-eps-converted-to.pdf}
}
\subfigure[Iterative PDC solve with filter.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/ComparisonDCWithFiltering-eps-converted-to.pdf}
}
\end{center}
\caption[Comparison between accuracy as a function of time for double precision calculations vs. single precision with and without filtering.]{Comparison between accuracy as a function of time for double precision calculations vs. single precision with and without filtering. The double precision result are unfiltered in each comparison and shows to be less sensitive to roundoff-errors. Medium steep nonlinear stream function wave in intermediate water $(kh,H/L)=(1,0.0502)$. Discretization is based on $(N_x,N_z)=(30,6)$, A courant number of $C_r=0.5$ and 6'$th$ order stencils.}
\label{ch7:filtering}
\end{figure}

Last, we demonstrate using a classical benchmark for propagation of nonlinear waves over a semi-circular shoal that single precision math is likely to be sufficient for achieving engineering accuracy. The benchmark is based on Whalin's experiment~\cite{ch7:Whalin1971} which is often used in validation of dispersive water wave models for coastal engineering applications, e.g., see provious work \cite{ch7:EBL08}. Experimental results exists for incident waves with wave periods $T=1,2,3\,s$ and wave heights $H=0.0390, 0.0150, 0.0136\,m$. All three test cases have been discretized with a computational grid of size ($257 \times 41 \times 7$) to resolve the physical dimensions of $L_x=35\,m$, $L_y=6.096\,m$. The still water depth decreases in the direction of the incident waves as a semi-circular shoal from $0.4572\,m$ to $0.1524\,m$ with an illustration of a snapshot of the free surface given in figure \ref{ch7:fig:whalinsetup}. The time step $\Delta t$ is computed based on a constant Courant number of $Cr=0.8$, where $c$ is the incident wave speed and $\Delta x$ is the grid spacing. Waves are generated in the generation zone $0 \leq x/L \leq 1.5$, where $L$ is the wave length of incident waves, and absorbed again in the zone $35 - 2L \leq x \leq 35\, m$.
%

A harmonic analysis of the wave spectrum at the shoal center line is computed and plotted in figure \ref{ch7:whalinresults} for comparison with the analogous results obtained from the experiments data. The three harmonic amplitudes are computed via a Fast Fourier Transform (FFT) method using the last three wave periods up till $t=50\, s$. There is a satisfactory agreement between the computed and experimental results and no noticeable loss in accuracy resulting from the use of single precision math.
%
\begin{figure}[!htb]
    \setlength\figureheight{0.3\textwidth}
    \setlength\figurewidth{0.32\textwidth}
%    \begin{center}
    \subfigure[Whalin test case at $t=50s$, wave period $T=2s$, and grid dimensions ($257 \times 41 \times 7$)]{\label{ch7:fig:whalinsetup}
    \includegraphics[width=0.5\textwidth]{Chapters/chapter7/figures/WhalinWaveSol_t50_T2_single.pdf}
    }
    \subfigure[$T=1s$]{
    {\scriptsize\input{Chapters/chapter7/figures/WhalinWaveHarmonics_T1_single.tikz}}
    }
    \subfigure[$T=2s$]{
    {\scriptsize\input{Chapters/chapter7/figures/WhalinWaveHarmonics_T2_single.tikz}}
    }
    \subfigure[$T=3s$]{
    {\scriptsize\input{Chapters/chapter7/figures/WhalinWaveHarmonics_T3_single.tikz}}
    }
%    \end{center}
    \caption[Harmonic analysis for the experiment of Whalin for $T=1,2,3\,s$ respectively.]{Harmonic analysis for the experiment of Whalin for $T=1,2,3\,s$ respectively. Measured experimental and computed results (single precision) are in good agreement. Test environment 1.}\label{ch7:whalinresults}
\end{figure}

\subsection{Acceleration via parallelism in time using 'Parareal'}\label{ch7:parareal}\index{parareal}

% Performance is no longer free  with new hardware
With modern many-core architectures, performance is no longer intrinsic and free on new generations of hardware. Added performance with new hardware now comes with the requirement of sufficient parallelism in the application to be accelerated. Methods, tricks and techniques for extracting parallelism in scientific applications are thus becoming increasingly relevant to enable added numerical accuracy as well as minimization of time to solution in the pursuit of faster and better analysis for engineering applications.

% Parareal as component
The Parareal algorithm has been introduced as a component in our in-house GPU library as described in section \ref{ch5:parareal}. The parareal library component makes it possible to easily investigate potential opportunities for further acceleration of the water wave model on a heterogeneous system and to assess practical feasibility of this algorithmic strategy for various wave types. We omit a detailed review of the Parareal algorithm and refer to details given in section \ref{ch5:sec:parareal} together with recent reviews \cite{ch7:LMT01,ch7:MS07,ch7:ASNP12}.

% Explain key parameters shortly
In section \ref{ch5:parareal} it is assumed that communication costs can be neglected and a simple model for the algorithmic work complexity is derived. It is found that there are four key discretization parameters for parareal that needs to be balanced appropriately in order to achieve high parallel efficiency. It is the number of coarse-grained time intervals $N$, the number of iterations $K$, the ratio between the computational cost of the coarse to the fine propagator $\mathcal{C}_\mathcal{G}/\mathcal{C}_\mathcal{F}$ and the ratio between fine and coarse time step sizes $\delta t/\delta T$.

% How to obtain speed-up
Ideally, the ratio $\mathcal{C}_\mathcal{G}/\mathcal{C}_\mathcal{F}$ is small and convergence happens in $k=1$ iteration. This is seldom the case though, as it requires the coarse propagator to achieve accuracy close to that of the fine propagator while at the same time being substantially cheaper computationally, these two objectives obviously being conflicting. Obtaining the highest possible speed-up is a matter of trade-off, typically, the more GPUs used, the faster the coarse propagator should be. The performance of parareal is problem and discretization dependent and as such one would suspect that different wave parameters influence the suitability of the method. This was investigated in \cite{ch7:ASNP12} and indeed the performance does change with wave parameters. Typically the method work better for deep water waves with low to medium wave amplitude.
%
\begin{figure}[!htb]
    \begin{center}
    \setlength\figureheight{0.35\textwidth}
    \setlength\figurewidth{0.37\textwidth}
    \subfigure[Performance scaling]{
%        {\small\input{Chapters/chapter7/figures/PararealScaletestGTX590.tikz}}
      \includegraphics[width=0.5\textwidth]{Chapters/chapter7/figures/PararealScaletestGTX590_conv.pdf}
    }
    \subfigure[Speedup]{
       % {\small\input{Chapters/chapter7/figures/PararealSpeedupGTX590.tikz}}
 \includegraphics[width=0.5\textwidth]{Chapters/chapter7/figures/PararealSpeedupGTX590_conv.pdf}
    }
    \end{center}
    \caption[Parareal absolute timings and parareal speedup.]{(a) Parareal absolute timings for an increasingly number of water waves traveling one wave length, each wave resolution is ($33\times 9$). (b) Parareal speedup for two to sixteen compute nodes compared to the purely sequential single GPU solver. Notice how insensitive the parareal scheme is to the size of the problem solved. Test environment 2.}\label{ch7:fig:DDPA_SPEEDUP}
\end{figure}
%

% What did we do and what are the results
We have performed a scalability study for parareal using 2D nonlinear stream function waves based on a discretization with $(N_x,N_z)=(33,9)$ collocation points, cf. figure \ref{ch7:fig:DDPA_SPEEDUP}. The study shows that moderate speedup are possible for this hyperbolic system. Using four GPU nodes a speedup of slightly more than two was achieved while using sixteen GPU nodes resulted in a speedup of slightly less than five. As noticed in figure \ref{ch7:fig:DDPA_SPEEDUP}, parallel efficiency decrease quite fast when using more GPUs. This limitation is due to the usage a fairly slow and accurate coarse propagator and linked to a known difficulty with parareal applied to hyperbolic systems. For hyperbolic systems, instabilities tend to arise when using a very inaccurate coarse propagator. This prevents using a large number of time sub-domains, as this by Amdahl's law also requires a very fast coarse propagator. The numbers are still impressive though, considering that it comes as additional speedup to an already efficient and fast code.
Performance results for the Whalin test case have also been reported in figure \ref{ch7:fig:whalinparareal}. There is a natural limitation to how much we can increase $R$ (the ratio between the complexity of the fine and coarse propagators), because of stability issues with the coarse propagator. In this test case we simulate from $t=[0,1]s$, using up to $32$ GPUs. For low $R$ and only two GPUs, there is no speedup gain, but for configuration with eight or more GPUs and $R\geq6$, we are able to get more than $2$ times speedup. Though these hyperbolic systems are not optimal for performance tuning using the parareal method, results still confirm that reasonable speedups are in fact possible on heterogenous systems.

\begin{figure}[!htb]
    \setlength\figureheight{0.3\textwidth}
    \setlength\figurewidth{0.32\textwidth}
%    \begin{center}
    \subfigure[Speedup]{
    {\small\input{Chapters/chapter7/figures/WhalinPararealSpeedup.tikz}}
    }
    \subfigure[Efficiency]{
    {\small\input{Chapters/chapter7/figures/WhalinPararealEfficiency.tikz}}
    }
%    \end{center}
    \caption[Parallel time integration using the parareal method.]{Parallel time integration using the parareal method. $R$ is the ratio between the complexity of the fine and coarse propagators. Test environment 2.}\label{ch7:fig:whalinparareal}
\end{figure}

% Comparison with DD
The parareal method is observed to be the most viable approach at speeding up small-scale problems due to the reduced communication and overhead involved. For sufficiently large problems, where sufficient work is available to hide the latency in data communication, we find the spatial domain decomposition method to be more favorable as it does not involve the addition of computational work and thereby allows for ideal speed-up, something usually out of reach for the parareal algorithm. An important thing to note here is that it is technically possible to extend the work and wrap the parareal method around the domain decomposition method, thereby obtaining a multiplication of the combined speed-up of both methods. This is of great interest in the sense that for any problem size, increasing the number of spatial sub-domains will eventually degrade speed-up due to the latency in communication of boundaries. By exploiting the latency robustness of parareal in conjunction with domain decomposition parallelism, it may be possible to go large scale for problems that would otherwise be to small to exploit a large number of GPUs. These investigations are subject to future work.

% Final remarks
Finally, we remark that the Parareal algorithm is also a fault tolerant algorithm. This property follows from the iterative nature of the algorithm and implies that a process can be lost during computations and regenerated without restarting the computations. This can be exploited to minimize total run time in case of such failures. 

%\newpage
\subsection{Towards real-time interactive ship simulation}\index{real-time simulation}

A fast GPU-accelerated ship hydrodynamic model is developed for real-time interactive ship simulation by modification of the unified potential flow model presented in section \ref{ch7:goveq}. The target scientific application is an interactive full mission marine simulator, where multiple ships controlled by naval officers can navigate in a near-realistic virtual marine environment. Full mission simulators are used for education and training of naval officers in critical manoeuvring operations and for evaluation of ship and marine infrastructure designs. To predict the motion of ships, a hydrodynamics model is required for prediction of forces by \eqref{ch7:forcecalc} which is affected by the kinematic properties of the model, cf. section \ref{ch7:dispkin}. The state-of-the-art for such a hydrodynamic model in todays real-time ship simulators is based on fast interpolation and proper scaling of experimental model data. The amount of experimental model data are limited with respect to hull forms and configurations, requiring the need for extrapolation that compromises the accuracy.

The objective of current and ongoing work is aimed at removing these limitations by replacing the existing hydrodynamic model and instead calculate at full-scale the flow field, wave field, ship-structure and ship-ship interaction forces in real-time using massive parallel computation technology. The potential flow (OceanWave3D) model presented in section \ref{ch7:goveq} is suitable as the modeling basis for this purpose since it is robust, accurate, efficient and scalable to arbitrary large domains. Furthermore, it can accurately account for dispersive waves in the range from shallow to deep waters in marine settings where the sea bed may be uneven.

The inclusion of ships in the wave model requires an approximate representation of such ships in the model. This ship approximations have to be chosen carefully with consideration to the computational performance of the numerical model to enable interactive real-time computing on today's modern hardware. For a first proof-of-concept linear wave and ship models is used as the model basis. This implied that wave heights and ship draft are linearized around the mean sea level $z=0$ m.

The physical domain for the computation is bounded from below by the seabed and from above by the free surface of the sea or the hull of the ship. If the ship is navigating in open water the ship physical spatial domain is unbounded in the horizontal direction and in confined waters it is bounded by harbour structures, etc. The representation of the physical domain surrounding the ship is done by finite truncation in the horizontal directions. The resulting time-varying finite physical domain $\Omega(t)$ is a box fixed to follow the ship motion, with the ship in the middle of the top face and with Cartesian coordinate axes aligned with the horizontal components of the forward and sideward ship directions and the upward is the opposite of the direction of gravitational acceleration.

A linearized model can be formulated in terms of kinematic and dynamic boundary conditions at the mean sea level \eqref{ch7:FSorigin} together with a Laplace problem subject to a variable depth kinematic boundary condition \eqref{ch7:eq:laplaceproblem}. The effects of ship hulls can be accounted for by splitting the scalar velocity potential function into steady $\phi_0$ and unsteady $\phi_1$ potentials such that $\phi = \phi_0 + \phi_1$ and with a quasi-static approximation of the pressure acting on the ship hull as suggested in \cite{ch7:LindbergEtAlIWWWFB2012}. This leads to a relatively simple ship model that enables a flexible and computationally efficient approximation of the ship geometry.
The steady potential $\phi_0$ is calculated using a double body approximation \cite{ch7:RavenJMST2010} of the ship and the unsteady potential $\phi_1$ is calculated by a linear free surface flow model. The resulting double-body flat-ship problem becomes
\begin{subequations}
\begin{align}
 \nabla^2 \phi_0 &= 0,  \quad -h \leq z \leq 0 \\
 \frac{\partial \phi_0}{\partial n} &= 0, \quad z = -h, \\
 \frac{\partial \phi_0}{\partial z} &= -U \frac{\partial \eta_0}{\partial x}, \quad z = 0,
 \end{align}
 \end{subequations}
where $U$ is the velocity of the ship. The unsteady linear water problem is used to calculate the unsteady wave evolution around and away from the ship
\begin{subequations}
\begin{align}
 \nabla^2 \phi_1 &= 0, \quad -h \leq z \leq 0, \\
 \frac{\partial \phi_1}{\partial n} &= 0, \quad z = -h, \\
\frac{\partial \phi_1}{\partial t}
+ (u_0 - U) \frac{\partial \phi_1}{\partial x}
+ v_0 \frac{\partial \phi_1}{\partial y}
 &= - \left(\frac{1}{2} \mathbf{u}_0 \cdot \mathbf{u}_0 + g \eta_0 + \frac{p_{ship}}{\rho} \right), \quad z=0, \label{ch7:quasistatic} \\
\frac{\partial \eta_1}{\partial t}
+ (u_0 - U) \frac{\partial \eta_1}{\partial x}
+ v_0 \frac{\partial \eta_1}{\partial y} &= \frac{\partial \phi_1}{\partial z}, \quad z=0,
 \end{align}
 \end{subequations}
where the pressure on the ship hull $p_{ship}$ is calculated explicitly based on a quasi-static approximation which is determined by assuming $\partial_t\phi_1\approx0$ and rewriting \eqref{ch7:quasistatic}. In general, a ship hull is a complex surface in three-dimensional space, but its draft can be approximated by a single valued function of the horizontal coordinates $\eta_0 = \eta_0(x,y)$ and the no-flux condition on the ship hull is approximated by a flat-ship approximation. Radiation boundary conditions are approximated by a Sommerfelt absorbing boundary condition \cite{ch7:DgayguiJolySJAM1994} on the vertical sides of the physical domain to let waves escape the domain.

The modified numerical model can still be based on flexible-order finite difference method as discussed in section \ref{ch7:sec:nummodel}. The computational bottleneck problem is the efficient solution of the Laplace problem twice which can be done efficiently by the GPU-accelerated iterative PDC method as explained in section \ref{ch7:PDCmethod}. A snapshot of the steady state wave field is provided the introduction to this chapter. Computed resistance curves for a Series 60 hull moving at forward speed corresponding to Froude number $F_n=0.316$ knots in calm water are compared to experimental data \cite{ch7:TodaEtAl1992} in figure \ref{ch7:fig:shiphydro} a). The computed Kelvin wave system is shown in figure \ref{ch7:fig:shiphydro} b). The computed results compare well with experiments at moderate ship Froudes numbers $F_n=U/\sqrt{gh}$ in the range 0.1-0.25 as expected for a linear model. The real-time constraint required to fulfil the interactive and visualization requirements can be met with the GPU-accelerated hydrodynamics model for problem sizes of approximately $10^6$ for ship Froudes numbers in the range 0.1-0.3. The modelling and real-time aspects will be addressed in more detail in ongoing work.

%
\begin{figure}[!htb]
\begin{center}
\subfigure[Hydrodynamic force calculations.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/figSeries60CB06Type7Resistance-eps-converted-to.pdf}
}
\subfigure[Kelvin pattern.]{
\includegraphics[width=0.45\textwidth]{Chapters/chapter7/figures/figSeries60CB06Type7kelvin-eps-converted-to.pdf}
}
\end{center}
\caption{Computed results. Comparison with experiments for hydrodynamics force calculations confirming engineering accuracy for low Froudes numbers.}
\label{ch7:fig:shiphydro}
\end{figure}
%

%\newpage
\section{Conclusion and future work}

We have presented implementation details together with several novel results on development of a new massively parallel and scalable tool for simulation of nonlinear free surface water waves on heterogenous hardware. The tool is based on the unified potential flow model referred to as OceanWave3D \cite{ch7:EBL08} which provide the basis for efficient and scalable simulation of water waves over uneven bottoms on arbitrary domain sizes. We have demonstrated in a few examples how we can accelerate performance by using single precision math without comprising accuracy. We have shown that performance can be accelerated by introducing concurrency in the time integration using the Parareal algorithm and for the first time in a heterogenous setup based on the use of multiple GPUs. Interestingly, we find that parallel computations using parareal may be more efficient than using conventional data-parallel distributed computations in a multi-GPU setup for moderate problem sizes. We have measured absolute performance and scalability on several of the most recent generations of Nvidia GPUs to detail the efficiency of the current code. This is useful to predict time to results as explained in \cite{ch7:EngsigKarupEtAl2011} and may be compared against other wave models.

Work in progress focus on extending the governing equations to account for lack of physics such as wave runup and wave breaking. Also, we plan to extend the domain decomposition method to unstructured grids of blocks that can be boundary-fitted to more general bottom-mounted structures to be able to address wave-structure problems, cf. \cite{ch7:EHBM06,ch7:EHBW08}. For example, this will provide the basis for simulations of wave transformations in large harbour areas or predict wave climates in near-coastal areas.

We anticipate that a tool based on the proposed parallel solution strategies will be useful for further advancement in fast and robust analysis techniques and large-scale simulation of free surface wave simulation (e.g. for use as an efficient far-field solver at large scales) and be a basis for next-generation wave models. We also expect that the tool can be useful for hybrid-solution strategies with local flow features possibly resolved by other models and for advancing state-of-the-art in fast physics-based wave-body simulations, e.g., ship-wave interactions in ship simulation where real-time constraints are imposed due to visualization. These subjects will be part of ongoing work addressing application aspects. 

\section{Acknowledgement}

This work was supported by grant no. 09-070032 from the Danish Research Council for Technology and Production Sciences. A special thank goes to Professor Jan S. Hesthaven for supporting parts of this work. Scalability and performance tests was done in GPUlab at DTU Informatics, Technical University of Denmark and using the GPU-cluster at Center for Computing and Visualization, Brown University, USA. Nvidia Corporation is acknowledged for generous hardware donations to facilities of GPUlab.

\putbib[Chapters/chapter7/biblio7]