# Real-World Oceanographic Simulations on the GPU
using a Two-Dimensional Finite-Volume Scheme

André R. Brodtkorb111Corresponding author<EMAIL_ADDRESS>Norwegian
Meteorological Institute, P.O. Box 43 Blindern, NO-0313 Oslo, Norway.  Oslo
Metropolitan University, Department of Computer Science, P.O. Box 4 St. Olavs
plass, NO-0130 Oslo, Norway.  Håvard Heitlo Holm  SINTEF Digital, Mathematics
and Cybernetics, P.O. Box 124 Blindern, NO-0314 Oslo, Norway.  Norwegian
University of Science and Technology, Department of Mathematical Sciences,
NO-7491 Trondheim, Norway.

###### Abstract

In this work, we take a modern high-resolution finite-volume scheme for
solving the rotational shallow-water equations and extend it with features
required to run real-world ocean simulations. Our contributions include a
spatially varying north vector and Coriolis term required for large scale
domains, moving wet-dry fronts, a static land mask, bottom shear stress, wind
forcing, boundary conditions for nesting in a global model, and an efficient
model reformulation that makes it well-suited for massively parallel
implementations. Our model order is verified using a grid convergence test,
and we show numerical experiments using three different sections along the
coast of Norway based on data originating from operational forecasts run at
the Norwegian Meteorological Institute. Our simulation framework shows perfect
weak scaling on a modern P100 GPU, and is capable of providing tidal wave
forecasts that are very close to the operational model at a fraction of the
cost. All source code and data used in this work are publicly available under
open licenses.

##  1 Introduction

The aim of this paper is to simulate realistic oceanographic scenarios using a
modern finite-volume scheme on GPUs. Modern operational ocean models, such as
the Regional Ocean Modeling System (ROMS) [1] and the Nucleus for European
Modelling of the Ocean (NEMO) [2], are based on solving the primitive
equations that describe conservation of mass, momentum, temperature, and
salinity. Even though these ocean models are heavily optimized, their complex
physical and mathematical properties make them very computationally demanding.
The NorKyst-800 model system [3], for example, is based on the ocean model
ROMS and covers the Norwegian coast with 800 m horizontal resolution. It is
run operationally by the Norwegian Meteorological Institute on a daily basis.
The complete model consists of 42 vertical layers with $2600\times 900$ grid
cells each, and a 24 hour forecast takes 45 minutes when using 256 CPUs [4].
It is unfeasible to use such models for flexible on-demand simulations during
search-and-rescue operations, oil spills, heavy lifting at sea, etc. In these
cases, fast on-demand simulations based on the shallow-water equations can
have a significant value for decision makers.

We simulate the shallow-water equations in a rotating frame of reference,
targeting efficient simulation of ocean currents based on the operational
NorKyst-800 data. Our model equations are essentially two-dimensional
simplifications of the primitive equations and assume constant density and a
vertically integrated momentum (see more detailed derivations in, e.g., Røed
[5]). Our work starts with the recent high-resolution finite-volume scheme by
Chertock et al. [6], and extends this scheme to support bottom shear stress,
wind forcing, projection-free spatial variations in the north vector,
spatially varying Coriolis parameter, land mask, and a moving wet-dry
boundary. We present how to reconstruct a piecewise bilinear bathymetry from
the piecewise constant bathymetry in the NorKyst-800 data, nesting of the
model into the NorKyst-800 data using linearly interpolated boundary
conditions, and generation of both higher and lower spatial resolution of the
initial conditions.

We simulate increasingly challenging scenarios along the Norwegian coast, and
compare with reference results from the operational model. As the shallow-
water equations are particularly well suited for tidal and storm-surge
predictions, we also compare our tidal-wave predictions at several locations
with the reference dataset and observations. All datasets and source code used
in this work are publicly available under open licenses.222The source code
used to create plots and results in this article is available in revision
7281301a3a286351b6e857add12c2fbd0165b229 of the gpu-ocean repository,
https://github.com/metno/gpu-ocean, released under GPL version 3. The plots
are created using Jupyter Notebooks [7] in the folder
demos/realisticSimulations, and the simulator code itself is under the
directory SWESimulators as a Python module. We have used the NorKyst-800
forecast for 2019-07-16 from the Norwegian Meteorological Institute in this
work, available on
https://thredds.met.no/thredds/dodsC/fou-
hi/norkyst800m-1h/NorKyst-800m_ZDEPTHS_his.an.2019071600.nc. The observed sea-
surface levels are obtained from Se Havnivå,
https://www.kartverket.no/sehavniva/, provided by the Norwegian Mapping
Authority, Hydrographic Service. The data sets are available under the
Creative Commons 4.0 BY International license as described on
https://www.met.no/en/free-meteorological-data/Licensing-and-crediting and
https://www.kartverket.no/en/data/Terms-of-use/. color=yellow!40]Håvard: Cite
the future zenodo release of the code

The rest of the paper is organized as follows: Section 2 gives a high-level
description of the shallow-water equations and how they can be solved
numerically. In Section 3, we present all the extensions that are required to
tackle relevant oceanographic problems. The performance and numerical accuracy
of our scheme are assessed in Section 4, where we show perfect weak scaling of
the GPU implementation and run a grid convergence test with the new terms
introduced in the scheme. We present operational-grade simulations in Section
5, and finally summarize the work in Section 6.

##  2 Mathematical formulation

We use the shallow-water equations with source terms to simulate ocean
dynamics. Since the equations are purely barotropic, they are not expected to
generate or preserve features such as mesoscale eddies in the long term. These
features arise from baroclinic instabilities, typically caused by variations
in temperature and salinity that are not included in the mathematical model,
but the equations are nevertheless capable of capturing the most important
physics required to simulate short-term ocean dynamics. We can write the
equations with source terms as

$\begin{split}\left[\begin{array}[]{c}h\\\ hu\\\
hv\end{array}\right]_{t}+\left[\begin{array}[]{c}hu\\\
hu^{2}+\tfrac{1}{2}gh^{2}\\\
huv\end{array}\right]_{x}&+\left[\begin{array}[]{c}hv\\\ huv\\\
hv^{2}+\tfrac{1}{2}gh^{2}\end{array}\right]_{y}\\\
\,=&\left[\begin{array}[]{c}0\\\ ghH_{x}\\\
ghH_{y}\end{array}\right]+\left[\begin{array}[]{c}0\\\
-ru\sqrt{u^{2}+v^{2}}/h\\\
-rv\sqrt{u^{2}+v^{2}}/h\end{array}\right]+\left[\begin{array}[]{c}0\\\ fhv\\\
-fhu\end{array}\right]+\left[\begin{array}[]{c}0\\\ \tau_{x}\\\
\tau_{y}\end{array}\right].\end{split}$ (1)

Here, $h$ is water depth, $hu$ momentum along the $x$-axis, and $hv$ momentum
along the $y$-axis (see also Figure 1). Furthermore, $g$ is the gravitational
acceleration, $H$ is the equilibrium ocean depth measured from a reference
sea-surface level, $r$ is the friction coefficient, $f$ is the Coriolis
parameter, and $\tau_{x}$ and $\tau_{y}$ are the wind stress along the $x$\-
and $y$-axis, respectively. In vector form, we can write the equations as

$Q_{t}+F(Q)_{x}+G(Q)_{y}=B(Q)+S(Q)+C(Q)+W(Q).$ (2)

Here, $Q=[h,hu,hv]^{T}$ is the vector of conserved variables and $B$, $S$,
$C$, and $W$ are source terms representing varying bathymetry, bed shear
stress, Coriolis force, and wind drag, respectively.

By using a finite-volume discretization on a Cartesian grid to solve (2), we
get the following semi-discrete formulation

$\displaystyle\frac{\partial Q}{\partial t}$ $\displaystyle=-\frac{1}{\Delta
x}\bigl{[}F_{i+1/2,j}-F_{i-1/2,j}\bigr{]}-\frac{1}{\Delta
y}\bigl{[}G_{i,j+1/2}-G_{i,j-1/2}\bigr{]}+B_{i,j}+S_{i,j}+C_{i,j}+W_{i,j}.$
(3)

To evolve the solution in time, we use a second-order strong stability-
preserving Runge-Kutta scheme [8],

$\begin{split}Q^{*}_{i,j}&=Q^{n}_{i,j}+\Delta tM(Q^{n})_{i,j}\\\
Q^{n+1}_{i,j}&=\tfrac{1}{2}\Bigl{(}Q^{n}_{i,j}+\left(Q^{*}_{i,j}+\Delta
tM(Q^{*})_{i,j}\right)\Bigr{)},\end{split}$ (4)

in which $M$ represents the right hand side of (3). The time step is
restricted by the CFL condition given by

$\Delta t\leq\frac{1}{4}\text{min}\left\\{\frac{\Delta
x}{\text{max}\left|u\pm\sqrt{gh}\right|},\frac{\Delta
y}{\text{max}\left|v\pm\sqrt{gh}\right|}\right\\},$ (5)

if we assume no Coriolis and wind drag. The dominant term for oceanographic
applications is typically the speed of gravity waves $\sqrt{gh}$, and we use a
conservative Courant number less than one in our simulations to account for
Coriolis and wind drag.

###  2.1 Numerical scheme

The starting point for our simulator is the scheme proposed by Chertock et al.
[6], from now on referred to as CDKLM. It is a modern central-upwind scheme
that uses a reconstruction of the physical variables to be able to capture
steady-state solutions important for physical oceanography, and includes
source terms for varying bathymetry and Coriolis. The original paper includes
several idealized test cases and examples demonstrating that the steady states
in geostrophic balance are well captured. However, there are several
shortcomings that need to be addressed before the scheme can be used for real-
world simulations.

The CDKLM scheme is based on the same principles as the MUSCL scheme [9] to
obtain a high-order, total variation diminishing discretization. For each cell
in the domain, we first reconstruct a piecewise (linear) function for our
physical variables (see Figure 1). At the interface between two cells, we
evaluate the slope limited function from the right and left cell, and use
these two values to compute the flux across the interface. The scheme uses the
Riemann-solver-free central-upwind flux function,

$\displaystyle
F_{i+1/2}=\frac{a^{+}F(Q_{i+1/2}^{l})-a^{-}F(Q_{i+1/2}^{r})}{a^{+}-a^{-}}+\frac{a^{+}a^{-}}{a^{+}-a^{-}}\left[Q_{i+1/2}^{r}-Q_{i+1/2}^{l}\right],$
(6)

to evaluate the numerical flux [10]. Here, $a^{+}$ and $a^{-}$ are the largest
positive and negative wave speeds at the interface, respectively.

When adding source terms, it becomes challenging to preserve steady states
such as “lake at rest”, i.e.,

$\frac{\partial Q}{\partial
t}=0\qquad\mathrm{if}\qquad\frac{\partial\eta}{\partial
x}=\frac{\partial\eta}{\partial y}=hu=hv=0,$ (7)

in which $\eta$ is the sea-surface deviation from a mean equilibrium level,
given by

$\eta=h-H.$ (8)

Kurganov and Levy [11] presented a discretization of the bottom slope source
term $B(Q)$ and a matching reconstruction based on water elevation to capture
such steady states. For oceanographic simulation, an important steady state is
the so-called geostrophic balance, in which the angular momentum caused by the
rotational reference frame is balanced by the pressure gradient:

$\begin{split}fhv-gh\frac{\partial\eta}{\partial x}=0,\qquad-fhu-
gh\frac{\partial\eta}{\partial y}=0.\end{split}$ (9)

The novel idea in the CDKLM scheme is to reconstruct the face values
$Q^{l}_{i+1/2}$ and $Q^{r}_{i+1/2}$ based on this geostrophic balance so that
the resulting flux $F_{i+1/2}$ also balances the Coriolis force for steady-
state solutions. This makes the scheme suitable for simulating ocean currents
in a rotating frame of reference.

(a)

(b)

(c)

(d)

Figure 1: a) The relationship between the physical variables $h$, $H$, $\eta$,
and $hu$ used in the shallow-water equations, here shown in one dimension. b)
The ocean state is discretized in terms of cell average values, whereas the
bathymetry is defined at the cell intersections. c) We find the slopes of the
ocean state within each cell. d) The flux between the cells is found from two
one-sided point values at each interface, as seen from the two adjacent cells.

##  3 Efficient simulation of real-world ocean currents

Our simulator is implemented on the GPU for efficiency. Explicit finite-volume
schemes normally have an embarrassingly parallel nature, and it is therefore
natural to implement them on such massively parallel accelerator hardware (see
e.g., [12, 13, 14, 15]). We use CUDA [16] coupled with PyCUDA [17] to access
the GPU from Python in this work, and we have previously shown that using
Python instead of C++ has a negligible performance impact when used like this
[18]. Using the Python ecosystem for tasks such as pre-processing input data
and post-processing results has significantly increased our productivity,
while still maintaining high performance.

A simulation starts by reading initial conditions, forcing terms, and other
grid data from NetCDF files hosted by _thredds_ servers at the Norwegian
Meteorological Institute. We continue by initializing the simulator, which
internally allocates memory and uploads initial conditions on the GPU.
Finally, the main loop steps the simulator forward in time, and occasionally
(e.g., every simulated hour) downloads and writes results to NetCDF files on
disk.

The rest of this section covers additions and improvements to the CDKLM scheme
required for simulating real-world scenarios.

###  3.1 Efficient parallel formulation of CDKLM

The original CDKLM scheme uses a recursive formulation of the potential
energies used for the well-balanced flux reconstruction. This turns the
reconstruction procedure into a global operation with a data-dependency
spanning the whole domain, making the scheme prohibitively expensive on
parallel architectures. The recursive expression can be calculated somewhat
efficiently using a prefix sum, as suggested by the authors themselves, but by
carefully reformulating the recursive terms it can also be rephrased as a
local operation.

The scheme reconstructs $h$, $hu$, and $hv$ on each side of a face to
calculate the fluxes in (6). However, to capture the rotating steady-state
solutions in geostrophic balance, the reconstruction is based on the potential
energies,

$K=g(\eta-V),\qquad L=g(\eta+U),$ (10)

instead of physical variables. Here, $U$ and $V$ are the primitives of the
Coriolis force, given by

$V_{x}=\frac{f}{g}v,\qquad U_{y}=\frac{f}{g}u.$ (11)

Note that setting $K_{x}=0$ and $L_{y}=0$ is equivalent to the geostrophic
balance from (9). In the following, we detail the reconstruction along the
abscissa, but the same derivations apply in the ordinate. The recursive terms
appear from the discretization of $U$ and $V$, which are defined on the cell
faces as

$V_{i+1/2,j}=V_{i-1/2,j}+\frac{\Delta x}{g}f_{i,j}v_{i,j},\qquad
V_{-1/2,j}=0.$ (12)

The Coriolis parameter $f_{k,j}$ is allowed to change both along the $x$ and
$y$ axis, even though $f$ only varies with latitude. This is because we will
use a spatially varying latitude, as described in more detail in Section 3.7.
Values of $V$ in the cell centers are obtained by taking the mean of the face
values, so that the cell-average value of $K$ from (10) becomes

$K_{i,j}=g\left(\eta_{i,j}-\frac{1}{2}\left(V_{i+1/2,j}+V_{i-1/2,j}\right)\right).$
(13)

By reconstructing the (limited) slopes of $K$ within each cell, the face
values of $h$ are found by combining (13) and (8), so that

$h_{i,j}^{E}=\frac{1}{g}\left(K_{i,j}+\frac{\Delta
x}{2}(K_{x})_{i,j}\right)+V_{i+1/2,j}+H_{i+1/2,j},$ (14)

and

$h_{i,j}^{W}=\frac{1}{g}\left(K_{i,j}-\frac{\Delta
x}{2}(K_{x})_{i,j}\right)+V_{i-1/2,j}+H_{i-1/2,j}.$ (15)

The reconstructed values (14) and (15) of $h$ at the cell faces quickly become
a performance bottleneck as long as they depend on calculating the recursive
relation in (12). As it turns out, however, we never need to explicitly
compute the recursive terms. First of all, notice that (14) and (15) only
depend on derivatives of $K$ and not on $L$. Similarly, the reconstructions of
$h_{i,j}^{N}$ and $h_{i,j}^{S}$ depend on $L$ and $L_{y}$. This means that we
are only interested in derivatives of $K$ (and $L$) in the same direction as
the recursive terms for $V$ (and $U$). The derivatives are limited using the
generalized minmod function using the backward, central, and forward
differences. The backward difference is given by

$\frac{K_{i,j}-K_{i-1,j}}{\Delta x}=\frac{g}{\Delta
x}\left(\eta_{i,j}-\frac{1}{2}\left(V_{i+1/2,j}+V_{i-1/2,j}\right)-\eta_{i-1,j}+\frac{1}{2}\left(V_{i-1/2,j}+V_{i-3/2,j}\right)\right).$
(16)

We see that $V_{i-1/2,j}$ cancels, and focus on the remaining $V$ values. By
applying the recursive expression in (12) twice on $V_{i+1/2,j}$ we get

$\begin{split}\frac{1}{2}\left(-V_{i+1/2,j}+V_{i-3/2,j}\right)=&\frac{1}{2}\left(-V_{i-1/2,j}-\frac{\Delta
x}{g}f_{i,j}v_{i,j}+V_{i-3/2,j}\right)\\\
=&\frac{1}{2}\left(-V_{i-3/2,j}-\frac{\Delta
x}{g}f_{i-1,j}v_{i-1,j}-\frac{\Delta x}{g}f_{i,j}v_{i,j}+V_{i-3/2,j}\right)\\\
=\,&-\frac{\Delta
x}{2g}\Big{(}f_{i-1,j}v_{i-1,j}+f_{i,j}v_{i,j}\Big{)}.\end{split}$ (17)

Inserted back into (16), the backward difference needed to evaluate
$(K_{i,j})_{x}$ can be written as

$\frac{K_{i,j}-K_{i-1,j}}{\Delta x}=\frac{g}{\Delta
x}\left(\eta_{i,j}-\eta_{i-1,j}-\frac{\Delta
x}{2g}\Big{(}f_{i-1,j}v_{i-1,j}+f_{i,j}v_{i,j}\Big{)}\right),$ (18)

which no longer contains any recursive terms. Similar derivations can be used
for the forward and central differences, and equivalently to obtain
$(L_{i,j})_{y}$.

After obtaining $(K_{x})_{i,j}$ we look at the reconstruction of $h_{i,j}^{E}$
from (14). We start by inserting the expression for $K$ from (13) into (14)
and gather all $V$-terms,

$\begin{split}h_{i,j}^{E}&=\eta_{i,j}-\frac{1}{2}\left(V_{i+1/2,j}+V_{i-1/2,j}\right)+\frac{\Delta
x}{2g}(K_{x})_{i,j}+V_{i+1/2,j}+H_{i+1/2,j}\\\
&=\eta_{i,j}+H_{i+1/2,j}+\frac{\Delta
x}{2g}(K_{x})_{i,j}+\frac{1}{2}\left(V_{i+1/2,j}-V_{i-1/2,j}\right).\end{split}$
(19)

Here, the values for $\eta_{i,j}$, $H_{i,j}$ and $(K_{x})_{i,j}$ are known,
and by using the recursive expression in (12) we get that

$\begin{split}\frac{1}{2}\left(V_{i+1/2,j}-V_{i-1/2,j}\right)&=\frac{1}{2}\left(V_{i-1/2,j}+\frac{\Delta
x}{g}v_{i,j}f_{i,j}-V_{i-1/2,j}\right)\\\ &=\frac{\Delta
x}{2g}f_{i,j}v_{i,j}.\end{split}$ (20)

Inserting (20) back into (19) we see that we get

$h_{i,j}^{E}=\eta_{i,j}+H_{i+1/2,j}+\frac{\Delta
x}{2g}\Big{[}(K_{x})_{i,j}+f_{i,j}v_{i,j}\Big{]},$ (21)

Again, the same derivation can be applied on the other three face values as
well. This enables us to express the reconstruction step in terms of only
local variables, which unlocks an embarrassingly parallel numerical algorithm.

###  3.2 Accuracy, precision and GPUs

Today’s GPUs can be categorized in two major classes: GPUs meant for the
consumer market, and GPUs meant for supercomputers and other professional
users. For the most part, these GPUs are close to identical, but there are
some key differences. One is price: professional GPUs can cost over ten times
more than equivalent consumer market GPUs, making the latter tractable from a
financial stand point. Another difference is in the feature set. Professional
GPUs offer more memory, error-correcting code (ECC) memory, and significantly
higher performance for features such as double precision.

Double precision is not an important feature on consumer market GPUs, and
typically offer just 3% of the performance of single-precision arithmetics.
The NVIDIA TITAN RTX, for example, has a performance of 16.3 teraFLOPS in
single precision, and 0.51 teraFLOPS in double precision. On professional
GPUs, on the other hand, the performance is 50% of single precision, which is
the same ratio as on CPUs. If we are able to utilize single-precision
calculations, we should therefore be able to get the highest performance at
minimum monetary cost. We have previously shown that using single precision is
sufficiently accurate for most simulation scenarios for the shallow-water
equations [19], and recent studies have shown a 40% increase in performance
for complex forecasting models when using single precision [20].

It is well known that floating-point numbers have a round-off error that
increases for larger numbers. A single-precision floating-point number is
represented as

$\displaystyle(-1)^{s}\cdot 2^{e-127}\cdot(1.0+f)$ (22)

in which $s$ is the sign bit, $e$ is represented using 8 bits, and $f$ is the
fractional part, represented by 23 bits. This formulation means that the
distance between two floating-point numbers is smallest close to zero, and
increases as $e$ increases. This is important to note in oceanographic
simulations, as the water depth, $h$, can often be thousands of meters,
whereas the relevant tidal wave height is perhaps 0.5 meters. The loss in
precision due to the water depth makes floating point prohibitively inaccurate
if we implement the shallow-water equations as formulated in (1). To counter
this loss in precision, we can possibly use double precision (at twice the
computational and memory cost), but a better alternative is to reformulate the
problem to represent the relevant physical quantities. This can be done by
basing our simulation on $\eta$ from (8) instead of $h$ and use
$[\eta,hu,hv]^{T}$ as our vector of conserved variables. By combining this
approach with single-precision floating-point numbers, we get more than
sufficient accuracy for the relevant oceanographic simulations.

When the water depth becomes close to zero, the calculation of the particle
velocity, $u=hu/h$ can become severely ill conditioned as the expression is
prone to large round-off errors. To counter this, Kurganov [21] and Kurganov
and Petrova [22] have proposed to _desingularize_ the term using expressions
such as

$\displaystyle u^{*}_{1}$
$\displaystyle=\frac{\sqrt{2}h(hu)}{\sqrt{h^{4}+\max(h^{4},\kappa^{4})}}$ (23)
$\displaystyle u^{*}_{2}$ $\displaystyle=\frac{h\cdot
hu}{h^{2}+\kappa^{2}},\qquad\text{or}$ (24) $\displaystyle u^{*}_{3}$
$\displaystyle=\frac{2h\cdot hu}{h^{2}+\max(h^{2},\kappa^{2})}.$ (25)

The suggestion of using $\kappa=\max(\Delta x^{4},\Delta y^{4})$ has some
obvious problems for real-world simulations with $\Delta x>1$, as pointed out
in [13]. Furthermore, the desingularization of $u$ here is not very well-
behaved. Instead, we propose to use the following formulation

$\begin{split}h^{*}&=\text{sign}(h)\max\Bigl{(}|h|,\min(h^{2}/(2\kappa)+\kappa/2,\kappa)\Bigr{)}\\\
u^{*}&=hu/h^{*}.\end{split}$ (26)

The advantages of this formulation is that (i) $\kappa$ now controls directly
what magnitudes of $h$ to desingularize, (ii) it improves the numerical
stability as it avoids the square root and the problematic $h^{4}$ term, (iii)
it yields a smooth transition between normal and desingularized values, and
(iv) the effective magnitude of $h^{*}$ is controlled. Figure 2 shows the
three desingularization strategies and clearly indicates that the
desingularized quantity $h$ is better behaved with our proposed approach.

(a)  (b)

Figure 2: Desingularization of the quantity $u=h/hu$ using existing and our
proposed method with $\kappa=1.0\cdot 10^{-11}$ for all cases. We compute the
value of $hu/h$ using the different approaches for $u=1$ and different values
of $h$. Notice how the single precision floating-point errors are clearly
visible for $u_{1}^{*}$, that $u_{2}^{*}$ has a distinct kink, and that
$u_{3}^{*}$ significantly underestimates the true magnitude of $u$.

###  3.3 Global wall boundary conditions

We have implemented wall and periodic boundary conditions to be able to run
test cases to assess the numerical correctness of our simulator. For
reflective von Neumann-type wall boundary conditions, we use so-called ghost
cells that mirror the ocean state across the boundary to enforce a zero flux
out of the domain. In the absence of Coriolis forces, wall boundary conditions
at $x=0$ are constructed by setting the ghost cell values according to

$\eta(-x,y)=\eta(x,y),\qquad hu(-x,y)=-hu(x,y),\qquad hv(-x,y)=hv(x,y).$ (27)

Naïve implementation of these conditions with the CDKLM scheme leads to
gravitational waves along the boundary. Special care is therefore required to
balance the Coriolis potential energies in the reconstruction of $h$, so that

$\begin{split}0&=h^{W}_{0,j}-h^{E}_{-1,j}\\\
&=\left(\eta_{0,j}+H_{-1/2,j}-\frac{\Delta
x}{2g}\Big{[}(K_{x})_{0,j}+f_{0,j}v_{0,j}\Big{]}\right)-\left(\eta_{-1,j}+H_{-1/2,j}+\frac{\Delta
x}{2g}\Big{[}(K_{x})_{-1,j}+f_{-1,j}v_{-1,j}\Big{]}\right),\end{split}$ (28)

using (21). Here, we see that $H$ cancels immediately, whereas $\eta$ cancels
by the use of (27). To make the last terms cancel, we need to set
$f_{-1,j}v_{-1,j}=-f_{0,j}v_{0,j}$, which means that in addition to (27), we
need to enforce

$f(-x,y)=-f(x,y).$ (29)

Note also that this results in $(K_{x})=0$ across the boundary from (18).

###  3.4 Moving wet-dry boundary and land mask

The original CDKLM scheme breaks down on land as the gravitational wave speed
$\sqrt{gh}$ becomes imaginary with negative depths. To support a moving wet-
dry boundary, we use a similar approach as Kurganov and Petrova [22] by
adjusting the slope of $h$ so that all the reconstructed point values become
non-negative. For example, if we get $h_{i,j}^{E}<0$ from (21), we adjust the
slope to become

$(K_{x})_{i,j}=-f_{i,j}v_{i,j}-\frac{2g}{\Delta
x}\Big{[}\eta_{i,j}+H_{i+1/2,j}\Big{]},$ (30)

and thus forcing $h_{i,j}^{E}=0$ instead.

The moving wet-dry boundary is required for simulating, e.g., tsunamis and
other phenomena in which the run-up is important. The approach nonetheless
violates the steady-state reconstruction, causing non-physical waves and often
also small time step sizes. For static wet-dry boundaries, we have therefore
implemented wall boundaries following a land mask, i.e., ensuring a zero flux
for the shore-line. In our approach, we explicitly set the flux between two
neighboring cells to zero if any one of them is masked. The mask is
represented by a unique no-data value in the bathymetry, thereby having no
increase in the memory footprint. It should be noted that both of these
approaches reduce the reconstruction to first order along the wet-dry
boundary.

###  3.5 Wind forcing and bed friction source terms

The original formulation of CDKLM includes the bed slope source term, $B(Q)$,
and Coriolis force, $C(Q)$, and we extend the scheme to also include the bed
shear stress, $S(Q)$, and wind forcing, $W(Q)$. Bed shear stress is
discretized using a semi-implicit formulation in the strong stability
preserving Runge-Kutta scheme,

$\begin{split}Q_{i,j}^{*}&=\left(Q^{n}_{i,j}+\Delta
tM(Q^{n})_{i,j}\right)/\left(1+\Delta tS(Q_{i,j}^{n})\right)\\\
Q^{n+1}_{i,j}&=\tfrac{1}{2}\Bigl{(}Q^{n}_{i,j}+\left(Q^{*}_{i,j}+\Delta
tM(Q^{*})_{i,j}\right)\Bigr{)}/\Bigl{(}1+\tfrac{1}{2}\Delta
tS(Q_{i,j}^{*})\Bigr{)},\end{split}$ (31)

instead of (4). This essentially means that we apply half of the bed friction
using $Q^{n}$, and half of the friction using the predictive step $Q^{*}$.

The wind stress comes from the atmospheric wind transferring momentum to the
water column. We take the approach of Large and Pond [23] and use

$\displaystyle\vec{\tau}=\frac{\rho_{a}}{\rho_{w}}C_{D}\,|W_{10}|\,W_{10},$
(32)

in which $\rho_{a}=1.225$ is the specific density of air, $\rho_{w}=1025$ is
the specific density of sea water, $W_{10}$ is the wind speed at 10 meters,
and $C_{D}$ is the drag coefficient computed as

$C_{D}=10^{-3}\begin{cases}1.2&\text{if}\quad|W_{10}|<11~{}\mathrm{m/s},\\\
0.49+0.065\,|W_{10}|&\text{if}\quad|W_{10}|\geq 11~{}\mathrm{m/s}.\end{cases}$
(33)

The NorKyst-800 model stores its output, including the wind forcing, every
hour. Because the time step of our simulator is significantly smaller than one
hour, we use linear interpolation in time to get the most accurate wind stress
source term for each time step. Given the wind stresses $W_{10}^{a}$ at time
$t_{a}$ and $W_{10}^{b}$ at time $t_{b}$, we compute the wind stress at time
$t\in[t_{a},t_{b}]$ as a convex combination of the two

$W_{10}^{n}=(1-s)W_{10}^{a}+sW_{10}^{b},\qquad s=(t-t_{a})/(t_{b}-t_{a}).$
(34)

In space, we use specialized GPU texture hardware for bilinear interpolation
(see e.g., [24]), resulting in a trilinearly interpolated wind stress. This is
highly efficient, as the GPU has dedicated hardware for bilinear spatial
interpolation that also includes caching. An added benefit is that the grid
size of the wind stress data does not need to match up with the underlying
simulation grid.

###  3.6 Nesting the model into NorKyst-800

A common technique in operational oceanography is to nest high-resolution
local models within lower-resolution models that cover a larger domain. For
instance, the NorKyst-800 model is nested within TOPAZ, which is a 12-16 km
resolution HYCOM-based coupled ocean and sea ice model with data assimilation
[25]. In this work, we nest our simulations into the NorKyst-800 model by
setting initial ocean state and boundary conditions appropriately. We use the
so-called flow relaxation scheme [26], in which we have a relaxation zone
between the internal computational domain and the external boundary
conditions. In this region, we solve the shallow-water equations as normal to
obtain an internal solution $Q^{\mathit{int}}_{i,j}$, and then gradually nudge
it towards an external solution $Q^{\mathit{ext}}$ using

$\begin{split}Q_{i,j}&=(1-\alpha)Q^{\mathit{int}}_{i,j}+\alpha
Q^{\mathit{ext}},\qquad\alpha=1-\tanh(d_{i,j}/d_{0}).\end{split}$ (35)

Here, $d_{i,j}$ is distance in number of cells to the external boundary, and
the parameter $d_{0}$ is typically set to 2 or 3 to control the fall-off of
the $\tanh$ relaxation function.

The implementation of the boundary conditions follows the same ideas as for
the wind stress, using textures for linear interpolation in time and space. We
use float4333Textures on GPUs are originally designed to hold the color of one
pixel as the color channels red, green, blue, and alpha, hence float4.
textures to hold our physical variables, $\eta,hu,hv$, because float3 is not
supported as a texture format. We furthermore only use two textures, as we
pack the north and south boundaries into one texture, and similarly for the
east and west.

###  3.7 Coriolis force

The original CDKLM scheme requires that the grid is aligned with the cardinal
directions, which significantly restricts the use of the scheme. For higher
latitudes this becomes especially pronounced, as both the latitude and
direction towards north will vary across the domain. We have therefore
extended the scheme to support both a spatially varying north-vector and a
spatially varying latitude, which is required to be able to run realistic
simulations for the areas covered by the NorKyst-800. To account for the
varying north vector, we project the momentum onto the local north and east
vectors,

$\displaystyle\begin{split}hu^{e}&=\vec{e}\cdot[hu,hv]^{T},\qquad\vec{e}=[\cos(\theta),-\sin(\theta)]\\\
hv^{n}&=\vec{n}\cdot[hu,hv]^{T},\qquad\vec{n}=[\sin(\theta),\cos(\theta)],\end{split}$
(36)

and use these vectors when computing the source term $C(Q)$ so that

$\displaystyle C(Q)=f\cdot\left[\begin{array}[]{c}0\\\
\vec{x}\cdot[hv^{n},-hu^{e}]^{T}\\\
\vec{y}\cdot[hv^{n},-hu^{e}]^{T}\end{array}\right].$ (40)

in which the vectors $\vec{x}=[\cos(\theta),\sin(\theta)]$ and
$\vec{y}=[-\sin(\theta),\cos(\theta)]$ project the result back into the
$(x,y)$-coordinate system. Note that we also need to use the same approach in
the reconstruction of $h$ as discussed in Section 3.1 to maintain the
rotational steady states.

Our spatially varying Coriolis parameter is computed from the latitude,
$l_{i,j}$, of a cell

$\displaystyle f_{i,j}=2\omega\sin(l_{i,j}),\qquad\omega=7.2921\cdot
10^{-5}\text{rad/s},$ (41)

or using a beta-plane model,

$\displaystyle
f_{i,j}=f_{\mathit{ref}}+\beta\cdot\Big{(}\vec{n}\cdot\big{[}(x_{i}-x_{\mathit{ref}}),(y_{j}-y_{\mathit{ref}})\big{]}^{T}\Big{)}.$
(42)

Here, $f$ is linearized around a reference point,
$(x_{\mathit{ref}},y_{\mathit{ref}})$, with a slope $\beta$ in the direction
of $\vec{n}$. Setting $\beta$ to zero yields a constant Coriolis parameter
throughout the domain.

Figure 3 demonstrates how the angle to north plays an important role in the
generation of planetary Rossby waves, caused by variations in the Coriolis
force using (42). The figure shows the resulting $\eta$ after long simulations
initialized with a rotating Gaussian bump in geostrophic balance in the middle
of the domain, similar to the case presented in Holm et al. [27]. Here, we
have used a domain consisting of $350\times 350$ cells with $\Delta x=\Delta
y=20\;\mathrm{km}$, a depth of $H=50\;\mathrm{m}$, and the center of the
domain corresponding to approximately $33^{\circ}$ north, with $f=8\cdot
10^{-4}$ and $\beta=2\cdot 10^{-11}$. The simulation ends after approximately
35 days. The figure clearly shows the effect of varying the north vector and
hence the Coriolis source term.

Figure 3: Planetary Rossby waves generated by a beta plane model for the
Coriolis force and with different directions for north. The simulations are
initialized by a rotating bump in the middle of the domain, which develops
into Rossby waves that propagates westwards with wave energy slowly
propagating eastwards. The different figures show that we get the same Rossby
waves for arbitrary orientation of our domain. Axes are given in 1000 km, and
the colormap shows $\eta$ in m.

###  3.8 Bathymetry

The bathymetry in the NorKyst-800 model is given as cell midpoint values, and
can not be used directly in our framework as CDKLM requires the bathymetry to
be defined as a piecewise bilinear surface specified by point values at the
cell intersections. It is still desirable that the averages of the resulting
four intersection values end up to be equal to the NorKyst-800 cell average
values, as large deviations can disturb the balanced relations in the provided
ocean state. Constructing such a piecewise bilinear surface from the cell
averages is an ill-posed problem, and a simple averaging of cell values
excessively smears important features such as deep fjords and straits.

We have therefore devised a pragmatic iterative algorithm, which is based on
reconstruction of slopes for each cell, similar to the reconstruction
procedure in our numerical shallow-water scheme. Our initial guess is computed
as follows:

  1. 1.

Compute the local gradient in each cell so that

$\displaystyle[H_{i,j}^{x},H_{i,j}^{y}]=\nabla H_{i,j}.$ (43)

  2. 2.

Evaluate the piecewise planar function at the four corners of each cell,

$\displaystyle\begin{split}H_{i-1/2,j-1/2}^{i,j}&=H_{i,j}-\tfrac{1}{2}H_{i,j}^{x}-\tfrac{1}{2}H_{i,j}^{y},\\\
H_{i+1/2,j-1/2}^{i,j}&=H_{i,j}+\tfrac{1}{2}H_{i,j}^{x}-\tfrac{1}{2}H_{i,j}^{y},\\\
H_{i-1/2,j+1/2}^{i,j}&=H_{i,j}-\tfrac{1}{2}H_{i,j}^{x}+\tfrac{1}{2}H_{i,j}^{y},\\\
H_{i+1/2,j+1/2}^{i,j}&=H_{i,j}+\tfrac{1}{2}H_{i,j}^{x}+\tfrac{1}{2}H_{i,j}^{y}.\end{split}$
(44)

  3. 3.

Average the estimate from the four cells meeting at one intersection point so
that

$\displaystyle
H^{*}_{i+1/2,j+1/2}=\tfrac{1}{4}\left(H_{i+1/2,j+1/2}^{i,j}+H_{i-1/2,j+1/2}^{i+1,j}+H_{i-1/2,j-1/2}^{i+1,j+1}+H_{i+1/2,j-1/2}^{i,j+1}\right)$
(45)

After obtaining the initial guess, we start an iterative procedure consisting
of two stages. The first stage minimizes the difference between computed and
true midpoint values, and the second dampens oscillations created by the first
stage. We start by computing the difference between our estimate and the
target bathymetry as

$\displaystyle\Delta H_{i,j}=H_{i,j}-H^{*}(i,j),$ (46)

in which $H^{*}(i,j)$ means that we evaluate the piecewise bilinear estimate
at midpoints. We then apply (43)–(45) to compute an estimate of the difference
at intersections, $\Delta H_{i+1/2,j+1/2}^{*}$, and subtract this difference
from our initial guess,

$\displaystyle H^{*}_{i+1/2,j+1/2}:=H^{*}_{i+1/2,j+1/2}+\Delta
H_{i+1/2,j+1/2}^{*}.$ (47)

This brings our updated midpoints of $H^{*}(i,j)$ closer to the target
$H_{i,j}$, but it also results in unwanted oscillations, as there is nothing
that limits the slopes of this piecewise bilinear surface. We therefore smooth
$H^{*}$, and repeat the procedure until the update between two subsequent
iterations is sufficiently small. Around ten iterations is typically
sufficient to create a piecewise bilinear surface that lies close to the
target. Figure 4 shows the result of our approach applied to the bathymetry
from NorKyst-800 in the innermost parts of Vestfjorden (see also Section 5.2).

Figure 4: Reconstruction of the bathymetry at intersection points. Starting
with the initial bathymetry $H_{m}$ defined at cell midpoints (upper left), we
compute the $H_{i}$ values at intersections (upper right). The lower-left
figure shows the convergence of our reconstruction algorithm, and the
reconstructed bathymetry is plotted together with the original in the lower-
right figure.

Figure 5: Widening of narrow straits, in which the land mask is eroded using
binary erosion. The left figure shows particle velocities in the Tjeldsund
strait in Lofoten/Vesterålen between the main land and Hinnøya (the bathymetry
is shown in Figure 4). Left shows the NorKyst-800 reference, center shows
original (reconstructed) bathymetry, and the right figure shows the result
with one grid cell land erosion. The strait width and geometry limits the
amount of water passing through it in the original bathymetry, but by eroding
one cell of the land mask the amount of water is much closer to the reference.

CDKLM is a full two-dimensional scheme, yet it exhibits a grid orientation
bias close to the land mask. For example, we have experienced that the scheme
has difficulties in propagating the momentum through narrow straits and canals
that are only one or two cells wide. To ameliorate this, we can use binary
erosion on the land mask to erode one row of cells at a time as shown in
Figure 5. We extrapolate $\eta$ into these areas using gray dilation (maximum
of neighbors), and set the water depth equal the land value, which in the
NorKyst-800 model is 5 meters. The process can be repeated to erode the land
mask even more, and is required for tidal waves to properly propagate in
features such as narrow fjords and straits. We have found that using single
iteration is sufficient to propagate tides in and out of the fjords in the
NorKyst-800 model.

###  3.9 Refinement and coarsening of grid

The possibility to run simulations at different resolutions is attractive to
quickly obtain low-resolution preliminary results, and also to resolve more of
the dynamics for more detailed and accurate results with higher resolution. We
have included functionality in our simulation framework to increase and
decrease the grid resolution with factors of two to support both of these
scenarios. We use a static grid refinement/coarsening performed on the initial
conditions themselves, but the same approach can be used for efficient
adaptive mesh refinement following the same approach as Sætra, Brodtkorb and
Lie [28].

The physical variables are represented by cell average values, and refinement
is based on a slope-limited reconstruction within each cell, which we evaluate
at the refined grid cell centers. As the bathymetry is specified by a
piecewise bilinear surface, the high-resolution representation of the
bathymetry evaluates the surface on the refined cell intersections. Low-
resolution representation of the physical variables on a coarsened grid is
simply the average of the original cell values, and intersection values for
the bathymetry are sampled from the corresponding original intersection
values. With a factor two coarsening, both these operations are well-defined.

Note that the boundary conditions and wind forcing terms are independent of
the mesh resolution due to our use of GPU textures.

##  4 Performance and accuracy of GPU implementation

We have previously evaluated the standard CDKLM scheme on several test cases
which target oceanographic dynamics relevant for the shallow-water equations
when implemented on a Cartesian grid [27]. In this work, however, we have
reformulated the scheme in terms of $\eta$ instead of $h$, added several new
source terms, and enabled varying latitude and north vector. It is therefore
important to revisit the grid convergence test. We also measure computational
performance of our simulation framework through assessing its weak scaling.

###  4.1 Convergence of the numerical scheme

To evaluate numerical convergence, we define a benchmark case covering
$500~{}\mathrm{km}\times 500~{}\mathrm{km}$. We model the bathymetry after the
Matlab peaks function,

$\displaystyle\text{peaks}(s,t)=3(1-s)^{2}e^{-s^{2}-(t+1)^{2}}-10(s/5-s^{3}-t^{5})e^{-s^{2}-t^{2}}-1/3e^{-(s+1)^{2}-t^{2}},\qquad
s,t\in[-3,3]$ (48)

which yields a non-symmetric smooth surface. We then define the equilibrium
ocean depth as

$H(x,y)=100+10\cdot\mathrm{peaks}\left(\frac{6x}{500}-3,\frac{6y}{500}-3\right),$
(49)

with $x$ and $y$ given in km, so that the depth varies between 34 and 181
meters (see Figure 6). We initialize the sea-surface level with a bump
centered in the domain,

$\displaystyle\eta(r)=\begin{cases}\tfrac{1}{2}\left(1+\cos\left(\tfrac{r}{c}\pi\right)\right),\quad&\mathrm{if}\;r\leq
c,\\\ 0,&\mathrm{if}\;r>c,\end{cases}$ (50)

in which $r=\sqrt{(x-250)^{2}+(y-250)^{2}}$ is the distance from the center of
the domain, and $c$ controls the size of the bump, set to 300 km in our case.
We assume that the North pole is located at
$(3000~{}\mathrm{km},3000~{}\mathrm{km})$, which means that the angle between
the $y$-axis and the north vector varies from $33.7^{\circ}$ to $56.3^{\circ}$
across the domain and that the latitude varies from $56.0^{\circ}$ to
$67.4^{\circ}$ north. These parameters roughly correspond to the North Sea.
Including land inevitably yield only first-order accurate fluxes, and we have
thus deliberately not included such areas in our grid convergence test.

The simulation is run for 800 seconds, after which we compute the difference
for the whole domain between the solution and a reference computed on a
$1024\times 1024$ grid. Figure 6 and Table 1 show that we achieve second-order
accuracy in both the $L_{1}$ and the $L_{2}$ norm, whilst the convergence in
$L_{\mathrm{inf}}$ is somewhat slower.

Figure 6: Convergence plot of our numerical scheme. The numerical case has a water disturbance consisting of a cosine bump at the center of the domain, and the bathymetry consists of the smooth Matlab peaks function. The Coriolis force varies spatially across the domain according to variations in latitude and direction towards north. The left figure shows the self convergence of the numerical error compared to the reference solution with $1024^{2}$ cells. To the right we show the bathymetry with the sea-surface level initially (top) and at the end of the simulation (bottom). Table 1: Error of computed solutions at different resolutions and numerical reduction order for each refinement. The time step is fixed throughout these simulations.  Resolution | $\Delta t$ | $L_{1}$ | Order | $L_{2}$ | Order | $L_{\mathrm{inf}}$ | Order
---|---|---|---|---|---|---|---
$16^{2}$ | $1600/16$ | 0.002898 | - | 0.006072 | - | 0.037367 | -
$32^{2}$ | $1600/32$ | 0.000923 | 1.65 | 0.001877 | 1.69 | 0.012651 | 1.56
$64^{2}$ | $1600/64$ | 0.000256 | 1.85 | 0.000502 | 1.90 | 0.002621 | 2.27
$128^{2}$ | $1600/128$ | 0.000058 | 2.14 | 0.000122 | 2.04 | 0.000741 | 1.82
$256^{2}$ | $1600/256$ | 0.000015 | 1.96 | 0.000033 | 1.87 | 0.000322 | 1.20
$512^{2}$ | $1600/512$ | 0.000005 | 1.58 | 0.000009 | 1.89 | 0.000136 | 1.25

###  4.2 Computational performance

To assess the computational performance, we continue to use the same
benchmark, but run the simulation twice as long. We also measure the cost of
downloading the results from the GPU to the CPU. Figure 7 shows the results
for domains with $16^{2}$ to $16~{}384^{2}$ cells, running on a Tesla P100
GPU. We see that for the larger domain sizes, we obtain perfect weak scaling,
as a factor two grid refinement leads to a factor eight increase in
computational complexity (the timestep size is cut in half due to the CFL
condition, and we have four times the number of cells). Cases with less than
$512^{2}$ cells are dominated by overheads and hence, the cost of a single
timestep is close to constant. Time for downloading data scales linearly with
the problem size, as is expected. Again, we see the same behavior for small
domains as for the computational performance.

Figure 7: Weak scaling of our simulator. For small domain sizes, the runtime is dominated by overheads, but for domain sizes larger than approximately $512\times 512$, we see that the computational time scales with the number of cells. For each doubling of our domain, we perform approximately eight times as many operations (four times as many cells, and twice the number of time steps). Table 2: Performance for our simulator. The simulation time is approximately eight times larger when we double the resolution. This is because we have four times as many cells and twice as many timesteps to reach the same end time. The download increases by a factor four, which is to be expected as we transfer four times as many data elements. Resolution | Simulation time | Factor | Download time | Factor
---|---|---|---|---
$512^{2}$ | 0.169710 | - | 0.000929 | -
$1024^{2}$ | 1.170814 | 6.9 | 0.002773 | 3.0
$2048^{2}$ | 9.133564 | 7.8 | 0.019642 | 7.1
$4096^{2}$ | 72.832028 | 8.0 | 0.077541 | 3.9
$8192^{2}$ | 582.599290 | 8.0 | 0.305405 | 3.9
$16384^{2}$ | 4695.904546 | 8.0 | 1.218880 | 4.0

##  5 Real-world simulations

Using the NorKyst-800 model system, the Norwegian Meteorological Institute
issues daily ocean forecasts for the complete Norwegian coast and makes them
publicly available through thredds servers. The results contain full 3D values
as well as depth-averaged values, wind forcing, bathymetry, longitudes,
latitudes, and angle towards north for each cell. This means that it is simply
a matter of accessing the file and cutting out the correct sections to be able
to initialize our simulator. To start a simulation within our framework, we
simply need the full url to the specific forecast we want to nest our
simulation within, and specify a desirable sub-domain.

We run our simulations on the three different sub-domains shown in Figure 8.
Our first simulation is in the Norwegian Sea with no land values within the
domain, and is therefore a fairly simple test for using real initial and
boundary conditions. The second case is centered on the Lofoten archipelago in
Northern Norway, which is dominated by a complex coastline consisting of
narrow fjords and a large number of islands, and serves as a challenging test
for our reconstruction of the bathymetry and land mask. Finally, we run almost
the complete domain, leaving only 25 cells in each direction to serve as
boundary conditions. With the entire Norwegian coastline, it becomes important
to account for the difference in the north vector throughout the domain (see
the longitudinal lines and north arrows in Figure 8), and we should see the
tides following the coast as Kelvin waves. The last case is also used to make
sea-level predictions at five different locations, indicated as stars in
Figure 8.

All three cases are run with three different resolutions; a high resolution
with a grid refined to 400 m cells, the original 800 m resolution of
NorKyst-800, and a low resolution grid consisting of 1600 m cells. We use an
adaptive time step, and update the time step every 20 minutes according to the
CFL condition in (5) with a Courant number 0.8. Each dataset provided from the
NorKyst-800 forecasts covers 23 hours, which means that we only have boundary
conditions to run simulations for the same time range.

Figure 8:  The domain covered by the NorKyst-800 model showing the equilibrium
depth off shore. The location of Case 1 in the Norwegian Sea (yellow
rectangle), Case 2 around the Lofoten archipelago (red rectangle), and Case 3
(beige rectangle) within the NorKyst-800 model domain, with red stars marking
the locations for sea-level predictions. The values on the $x$\- and $y$-axes
are given in km, and the latitude-longitude grid shows how the direction
towards north changes throughout the domain.

###  5.1 Case 1: Norwegian Sea

The first case runs an open water domain in the Norwegian Sea and tests our
ability to use the NorKyst data as initial and boundary conditions. The sea-
surface level changes significantly during the simulation time range as the
tidal waves enter and exit the domain, and this dynamic is largely dependent
on the boundary conditions. The velocity field on the other hand, is dominated
by already present features of slowly rotating dynamics, and will not change
much during the simulation time range. This enables us to test both that the
initial conditions are correctly captured, and how well our simulator
maintains complex rotating structures over time.

Figure 9 shows simulation results of the Norwegian Sea for all three different
grid resolutions after 23 hours, compared to the reference solution from
NorKyst-800. First, we see that the contour lines for the sea-surface level in
the left column show the same values at similar places for all simulations.
This shows that the boundary conditions allow tidal waves to propagate
correctly in and out of the domain as intended. The right column shows the
corresponding particle velocities. Note that the initial state for both the
high-resolution and original runs are visually inseparable from the
NorKyst-800 model state, whereas the low-resolution initial conditions
slightly smear out the sharper features for the velocity. At 23 hours, we see
that the high-resolution model to a large degree maintains the sharp features
from NorKyst-800. With the original resolution, however, the particle
velocities resemble a slightly blurred image of the reference solution, and in
the low-resolution simulation, the features are blurred even more. Also, we
are starting to see clear signs of grid effects both internally and from the
boundary.

(a)

(b)

(c)

(d)

Figure 9: Case 1 of the Norwegian Sea after 23 hours, comparing the reference
solution from NorKyst-800 with our simulations using three different grid
resolutions. All simulations obtain sea-surface level (left) similar to the
reference solution, meaning that the boundary conditions allow the tides to
correctly enter and exit our domain. From the particle velocities (right), we
see that only the high-resolution model is able to maintain sharp features,
whereas the low-resolution simulation has strong signs of grid effects. The
values on the $x$\- and $y$-axes are in km relative to the location in the
complete domain.

###  5.2 Case 2: Lofoten

(a)

(b)

(c)

(d)

Figure 10: Case 2 around the Lofoten archipelago after 23 hours, comparing the
reference solution from NorKyst-800 with our simulations using three different
grid resolutions. The values on the $x$\- and $y$-axes are in km relative to
the location in the complete domain. Our sea-surface levels (left column) are
similar but slightly higher within Vestfjorden compared to the reference
solution. From the particle velocities (right), we see that we get stronger
currents than the reference. The level of details in our simulations are
stronger with higher grid resolution, whereas the low-resolution simulation
blurs the structures and are influenced by grid effects.

The second case is centered around the Lofoten archipelago in Northern Norway.
This region has a complex coast line consisting of a large number of fjords,
straits and islands. In addition to the aspects tested in the first case, this
case tests the reconstruction of the bathymetry at intersections and the
handling of dry zones through the complex land mask. We expect to see a Kelvin
wave following the coast, but because the Lofoten archipelago protrudes into
the Atlantic, some of the wave will be partially trapped and amplified in
Vestfjorden between Lofoten and the main land. Furthermore, as can be seen
from the depth map in Figure 8, the edge of the continental shelf crosses
through the upper right corner of the sub-domain, and we can here see a line
of eddies driven by baroclinic instabilities.

We again run simulations with the three different grid resolutions, and Figure
10 shows the final simulation state after 23 hours, with a high tide in
Vestfjorden. All three simulations show that the contour lines for the sea-
surface level are slightly off compared to the reference solution, with
marginally more water into the fjord. Similarly as in Case 1, the low-
resolution simulation is unable to preserve the fine structures in the
currents and shows grid effects. The original and high grid resolutions are
also unable to preserve the exact eddy structure as found in the reference
solution, but we can instead more clearly see a consistent coastal current
with natural irregularities. This is consistent with what we can expect using
a barotropic model. These results also show a larger area with high particle
velocity around the smallest islands in our simulations, possibly due to the
reconstructed bathymetry, which widens some straits by the landmask erosion.

###  5.3 Case 3: Norway

Our third case consists of almost the complete domain covered by the
NorKyst-800 model. We leave out the 25 outermost cells and use these as
boundary conditions. With the complete domain, we cover a large enough area to
see the Kelvin waves forming in the middle of the domain and travelling north.

Figure 11 shows the simulation results after 23 hours in the same manner as
the earlier figures. We now see from the sea-surface level (left column) that
there is a low tide in the middle of the domain, whereas there are high tides
both in the north at the Russian border and in the Oslo fjord in the lower
left corner. The contour lines are similar for all three grid-resolutions, and
the particle velocities (right column) are in general consistent with the
reference solution from NorKyst-800. Due to the scale of the domain, it is
harder to compare the fine details in Figure 11, but the large scale features
from NorKyst-800 are found in both the high and original grid resolutions. As
we have seen in the other cases, the low-resolution simulation maintains less
details. It should be noted that by zooming into the same areas simulated in
Case 1 and 2, we have seen that the results from running the complete domain
are similar within each resolution as to the results shown in Figures 9 and
10.

(a)

(b)

(c)

(d)

Figure 11: Simulation result for Case 3 containing almost the complete domain used by NorKyst-800 after 23 hours. The top row shows the reference solution from NorKyst-800, followed by our simulation results using three different grid resolutions below. The values on the $x$\- and $y$-axes are in km relative to the location in the complete domain. The sea-surface levels (left) show how the tide varies along the coast, with low tide in the Oslofjord (lower left corner), high tide in the middle of the domain, and low tide again towards the border with Russia (lower right corner). Our results are consistent with the reference solution for all grid resolutions. The particle velocities (right column) are qualitatively similar for the reference solution (NorKyst-800), and our simulations with high and original resolution. Even at this scale, it is apparent that the low resolution grid smears the features in the ocean currents.  Table 3: Timing results different resolutions of running the complete Norwegian coast (Case 3) with the original, high and low resolution on different computer types. Resolution | Num cells | Average $\Delta t$ |  Laptop GeForce 840M Maxwell (2014) |  Desktop GeForce GTX780 Kepler (2013) |  Server Tesla P100 Pascal (2016)
---|---|---|---|---|---
Low | 541 875 | 1.85 s |  9m 22s |  1m 29s |  47s
Original | 2 167 500 | 0.92 s |  53m 49s |  10m 42s |  3m 13s
High | 8 670 000 | 0.46 s |  6h 10m 09s |  1h 23m 22s |  23m 21s

Table 3 reports the run time for simulating this case on a laptop, a desktop,
and a server GPU. The laptop contains a dedicated GeForce 840M GPU, which is
powerful for a laptop but on the low end of the GPU performance spectrum. It
completes the simulation in slightly less than 54 minutes for the original
grid resolution of NorKyst-800, which is comparable to the time required for
256 CPUs to run the NorKyst-800 simulation. On a desktop with a slightly old
GeForce 780GTX gaming GPU, the same simulation is carried out in approximately
10 minutes. Tesla P100 is a GPU often found in supercomputers, thus belonging
to the very high-end of the performance (and price) spectrum, and runs the
23-hour forecast in approximately 3 minutes. Note that for the high-resolution
scheme, the time step is on average just 0.46 seconds, meaning that the
results shown in Figures 9–11 are obtained after 180 000 time steps, using
only single-precision floating-point operations.

###  5.4 Tidal forecast verification and validation

The oceanographic phenomena best captured and maintained by the shallow-water
equations are arguably tides and storm surges. We therefore use Case 3 to
generate tidal forecasts for five locations along the Norwegian coast (see
Figure 8) and compare these with the hourly values from the NorKyst-800
dataset. The locations have been chosen based on the availability of actual
sea-level observation stations, and we use these to show the realism in the
forecasts generated by both NorKyst-800 and the current model. It should be
noted that the sea level variation depends on more physical processes than
those modeled by NorKyst-800 and our code, and we therefore expect a
discrepancy between the forecasts and the observations. For us, it is
therefore most relevant to use the NorKyst-800 result as a reference solution.
The five locations (see Figure 8) are Bergen and Honningsvåg close to the open
sea; Oslo and Narvik located in the innermost parts of long fjords; and Bodø,
which is affected by how the Lofoten archipelago catches the tidal waves.

(a)

(b)

(c)

(d)

(e)

Figure 12: Tidal predictions, reference NorKyst-800 results, and official
gauge measurements at selected locations generated by the simulations
discussed in Case 3. All grid resolutions capture the tide as predicted by
NorKyst-800 with only minor discrepancies at coastal locations, such as
Bergen, Bodø and Honningsvåg. The grid resolution comes much more into play
deep inside the fjords, such as at Oslo and Narvik. Here, the low-resolutions
grid no longer manages to represent the bathymetry well enough for the tide to
enter the fjord correctly (especially in the very narrow Oslo fjord), but we
still get fair results with the high-resolution grid. The weakly dotted line
is actual observed sea-surface level at the respective locations in the time
range of the forecast. It should be noted that official tidal forecasts have
additional physical terms, which is the main reason why the NorKyst-800 is
off.

Figure 12 shows the tidal forecasts generated from Case 3 for each of the five
locations, sorted from south to north. First of all, we see that our forecasts
at Honningsvåg are very well in accordance with NorKyst-800. The same is the
case in Bergen, even though our simulations here consistently give slightly
too low values for $\eta$. At Bodø, we see that NorKyst-800 gives a delay in
the transition from high to low tide between 14 and 16 hours, and we can see
weak indications of this in our high-resolution simulation as well. All these
locations are at the coast, which makes the tides here relatively easy to
capture, but it should still be noted that they are all protected by islands.
The tides in Oslo and Narvik are seen to be harder to capture, as these two
cities are located at the innermost parts of two long fjords, and we get
larger differences between the different grid resolutions. In Narvik, we get a
large improvement by using the high-resolution grid over the original and low
resolutions, but overall, our results are not too far from the reference
solution. Oslo is a particularly hard case, as it is positioned within a very
narrow and shallow fjord. This leads to larger relative discrepancies here,
even though the absolute difference is roughly the same as in Narvik. We see
that the shape of our high-resolution result is a bit different than for the
reference solution. The tidal signal predicted by the low-resolution grid in
Oslo is very weak, as the fjord almost becomes closed off at this resolution.
Finally, note that our solutions are well behaved already from the start,
meaning that we manage to initialize our simulations in a balanced state.

##  6 Summary

We have presented a GPU simulation framework suitable for real-world
oceanographic applications. The framework is based on a modern high-resolution
finite-volume method for the shallow-water equations [6] with new adaptations
and extensions to handle moving wet-dry fronts, land mask, bottom friction,
wind forcing, varying north vector, varying latitude, and external boundary
conditions. The numerical algorithm is also improved to facilitate for
efficient implementation on GPUs using single-precision floating-point
operations and includes an efficient reformulation of a problematic recursive
term. The framework is designed to be initialized from operational ocean
forecasts issued by the Norwegian Meteorological Institute from a three-
dimensional ocean model.

We have presented second-order grid convergence for a synthetic but
challenging benchmark containing complex topography and varying north vector.
We have also shown that the computational performance follows the expected
weak scaling.

Finally, we have validated the simulation framework through three different
real-world cases along the Norwegian coastline, initialized from ocean
forecasts produced by the NorKyst-800 operational model with an
$800~{}\mathrm{m}\times 800~{}\mathrm{m}$ horizontal resolution. Our results
demonstrate that the framework manages to maintain fine rotational structures,
especially when increasing the grid resolution to $400~{}\mathrm{m}\times
400~{}\mathrm{m}$. On a coarse grid (1600 m grid cells), however, the features
are lost and the simulations get dominated by grid effects. The tidal
forecasts show that we are able to predict the observed tides relatively well
at coastal locations on any of the three grid resolutions used, when compared
to the NorKyst-800 reference data. We see, however, that within long fjords,
our results are improved when increasing the grid resolution, compared to
using the original horizontal resolution of the NorKyst-800 model.

### Acknowledgements

This research has been funded by the Research Council of Norway under grant
number 250935 (GPU Ocean). The GPU Ocean core project team consists of Göran
Broström, Kai Christensen, Knut-Andreas Lie, and Martin Lilleeng Sætra, and
the authors are sincerely grateful for their collaboration and discussions,
which have influenced and enabled this work. The authors also thank Øyvind
Sætra and Jon Albretsen for feedback on drafts of the manuscript. The GPU
Ocean project has received support in form of compute time on UNINETT Sigma2 -
the National Infrastructure for High Performance Computing and Data Storage in
Norway under project number nn9550k. The authors declare that we have no
competing interests related to this work.

### CRedIT author statement

A. R. Brodtkorb: Conceptualization, Methodology, Software, Investigation,
Writing - Original Draft, Writing - Review & Editing, Visualization, Project
administration, Supervision, H. H. Holm: Conceptualization, Methodology,
Software, Investigation, Writing - Original Draft, Writing - Review & Editing,
Visualization.

## References

  * [1] A. Shchepetkin and J. McWilliams.  The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model.  Ocean Modelling, 9(4):347–404, 2005.
  * [2] G. Madec and the NEMO team.  NEMO ocean engine.  Note du Pôle de modélisation, Institut Pierre-Simon Laplace (IPSL), France, No 27, ISSN No 1288-1619, 2008.
  * [3] J. Albretsen, A. Sperrevik, A. Staalstrøm, A. Sandvik, and F. Vikebø.  NorKyst-800 report no. 1: User manual and technical descriptions.  Technical Report 2, Fisken og Havet, Institute of Marine Research, 2011.
  * [4] K. Christensen.  Head of Division for Ocean and Ice, Norwegian Meteorological Institute.  [personal email communication].
  * [5] L. Røed.  Atmospheres and Oceans on Computers.  Springer International Publishing, 2019.
  * [6] A. Chertock, M. Dudzinski, A. Kurganov, and M. Lukácová-Medvidová.  Well-balanced schemes for the shallow water equations with Coriolis forces.  Numerische Mathematik, Dec 2017.
  * [7] T. Kluyver, B. Ragan-Kelley, F. Pérez, B. Granger, M. Bussonnier, J. Frederic, K. Kelley, J. Hamrick, J. Grout, S. Corlay, P. Ivanov, D. Avila, S. Abdalla, and C. Willing.  Jupyter notebooks - a publishing format for reproducible computational workflows.  In F. Loizides and B. Schmidt, editors, Positioning and Power in Academic Publishing: Players, Agents and Agendas, pages 87–90. IOS Press, 2016\.
  * [8] S. Gottlieb, C.-W. Shu, and E. Tadmor.  Strong stability-preserving high-order time discretization methods.  SIAM Review, 43(1):89–112, 2001.
  * [9] B. van Leer.  Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method.  Journal of Computational Physics, 32(1):101 – 136, 1979.
  * [10] A. Kurganov, S. Noelle, and G. Petrova.  Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton–Jacobi equations.  SIAM Journal on Scientific Computing, pages 707–740, 2001.
  * [11] A. Kurganov and D. Levy.  Central-upwind schemes for the Saint-Venant system.  Mathematical Modelling and Numerical Analysis, 36:397–425, 2002.
  * [12] M. de la Asunción, J. Mantas, and M. Castro.  Simulation of one-layer shallow water systems on multicore and CUDA architectures.  The Journal of Supercomputing, 58(2):206–214, Nov 2011.
  * [13] A. Brodtkorb, M. Sætra, and M. Altinakar.  Efficient shallow water simulations on GPUs: Implementation, visualization, verification, and validation.  Computers & Fluids, 55(0):1–12, 2012.
  * [14] P. Parna, K. Meyer, and R. Falconer.  GPU driven finite difference WENO scheme for real time solution of the shallow water equations.  Computers & Fluids, 161:107 – 120, 2018.
  * [15] X. Qin, R. LeVeque, and M. Motley.  Accelerating an adaptive mesh refinement code for depth-averaged flows using GPUs.  Journal of Advances in Modeling Earth Systems, 11(8):2606–2628, 2019.
  * [16] NVIDIA.  NVIDIA CUDA C programming guide version 10.1, 2019.
  * [17] A. Klöckner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov, and A. Fasih.  PyCUDA and PyOpenCL: A scripting-based approach to GPU run-time code generation.  Parallel Computing, 38(3):157–174, 2012.
  * [18] H. Holm, A. Brodtkorb, and M. Sætra.  GPU computing with Python: Performance, energy efficiency and usability.  2019\.  [in review].
  * [19] A. Brodtkorb, T. Hagen, K.-A. Lie, and J. Natvig.  Simulation and visualization of the Saint-Venant system using GPUs.  Computing and Visualization in Science, 13(7):341–353, 2010.
  * [20] F. Váňa, P. Düben, S. Lang, T. Palmer, M. Leutbecher, D. Salmond, and G. Carver.  Single precision in weather forecasting models: An evaluation with the ifs.  Monthly Weather Review, 145(2):495–502, 2017.
  * [21] A. Kurganov.  Finite-volume schemes for shallow-water equations.  Acta Numerica, 27:289–351, 2018.
  * [22] A. Kurganov and G. Petrova.  A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system.  Communications in Mathematical Sciences, 5(1):133–160, 03 2007\.
  * [23] W. Large and S. Pond.  Open ocean momentum flux measurements in moderate to strong winds.  Journal of Physical Oceanography, 11(3):324–336, 1981.
  * [24] A. Brodtkorb, C. Dyken, T. Hagen, J. Hjelmervik, and O. Storaasli.  State-of-the-art in heterogeneous computing.  Journal of Scientific Programming, 18(1):1–33, May 2010.
  * [25] J. Xie, L. Bertino, F. Counillon, K. Lisæter, and P. Sakov.  Quality assessment of the TOPAZ4 reanalysis in the Arctic over the period 1991–2013.  Ocean Science, 13(1):123–144, 2017.
  * [26] H. Davies.  A lateral boundary formulation for multi-level prediction models.  Quarterly Journal of the Royal Meteorological Society, 102(432):405–418, 1976.
  * [27] H. Holm, A. Brodtkorb, K. Christensen, G. Broström, and M. Sætra.  Evaluation of selected finite-difference and finite-volume approaches to rotational shallow-water flow.  Communications in Computational Physics, 2019.  [to appear].
  * [28] M. Sætra, A. Brodtkorb, and K.-A. Lie.  Efficient GPU-implementation of adaptive mesh refinement for the shallow-water equations.  Journal of Scientific Computing, 63(1):23–48, April 2015.