GPU-acceleration of the Discontinuous Galerkin Shallow Water Equations Model (DG-SWEM) with OpenACC

TL;DR

This paper demonstrates how using OpenACC and Unified Memory simplifies porting the DG-SWEM storm surge model to GPUs, achieving significant performance improvements on NVIDIA hardware while maintaining a unified codebase.

Contribution

It introduces a GPU porting approach for DG-SWEM using OpenACC and Unified Memory, avoiding complex code conversions and enabling efficient storm surge simulations.

Findings

GPU implementation outperforms CPU version in storm surge scenarios

OpenACC simplifies GPU porting process for complex models

Unified Memory enables seamless CPU-GPU code compatibility

Abstract

This paper presents a porting of {DG-SWEM}, a first-order discontinuous Galerkin solver for storm surge based on the Advanced Circulation Model (ADCIRC), to NVIDIA GPUs. Time-explicit discontinuous Galerkin methods contain a large number of degrees of freedom but have been shown to exhibit a large amount of data parallelism due to the loose coupling between elements, and thus are naturally mapped to the GPU architecture. A previous framework in porting DG-SWEM to GPUs required converting subroutines from Fortran to C++ to be used with CUDA C++. By using OpenACC and Unified Memory, we simplify the porting process and maintain a single codebase for both CPU and GPU versions. We test the code using a large Hurricane Harvey scenario on NVIDIA's Grace Hopper chip, and compare the GPU code's performance on multiple H200 nodes to the CPU version on the same amount of Grace CPU nodes.

Tables

Table 1. Table 1 : A Grace Hopper (GH) GPU node

CPU	NVIDIA Grace, 72 cores @ 3.1 GHz
CPU memory	120 GB LPDDR5X
GPU	NVIDIA H100
No. of SMs	132
FP64	34 TFLOPs/s
FP32	67 TFLOPs/s
GPU memory	96 GB HBM3
GPU memory bandwidth	4 TB/s
NVLink-C2C CPU-GPU bandwith	900 GB/s

Table 2. Table 2 : A Grace Grace (GG) CPU node

CPU	NVIDIA Grace, 144 cores @ 3.4 GHz
FP64	7 TFLOPs/s
CPU memory	237 GB LPDDR5
CPU memory bandwidth	850 GiB/s

Table 3. Table 3 : Compiler flags used for the single-precision option.

Version	Flags
GPU	mpif90 -r4 -O3 -cuda -acc -Mextend -Mpreprocess -gpu=unified
CPU	mpif90 -r4 -O3 -Mextend -Mpreprocess

Table 4. Table 4 : Primary data arrays in the code.

Array	Meaning
ZE[k,e,irk]	$k^{th}$ mode of elevation solution at element $e$ on RK stage irk
QX[k,e,irk]	$k^{th}$ mode of x-velocity solution at element $e$ on RK stage irk
QY[k,e,irk]	$k^{th}$ mode of y-velocity solution at element $e$ on RK stage irk
RHS_ZE[k,e,irk]	Contribution to the $k^{th}$ mode of $d\zeta /dt$ at element $e$ on RK stage irk
RHS_QX[k,e,irk]	Contribution to the $k^{th}$ mode of $d{Q}_x/dt$ at element $e$ on RK stage irk
RHS_QY[k,e,irk]	Contribution to the $k^{th}$ mode of $d{Q}_y/dt$ at element $e$ on RK stage irk
PHI_AREA[k,i]	$k^{th}$ basis evaluated at area quadrature point i
PHI_EDGE[k,i]	$k^{th}$ basis evaluated at edge quadrature point i
XFAC[k,i,e]	Weight for the $k^{th}$ basis of element $e$ in x-direction at area quadrature point i
YFAC[k,i,e]	Weight for the $k^{th}$ basis of element $e$ in x-direction at area quadrature point i
SRFAC[k,i,e]	Weight for the $k^{th}$ basis of the source term on element $e$ at area quadrature point i
EDGEFAC[k,i,l]	Weight for the $k^{th}$ basis of the source term on element $e$ at area quadrature point i
EDGE_ze[k,l,e]	Numerical flux contribution to the $k^{th}$ basis of $\hat{\zeta }$ at local edge $l$ of element $e$
EDGE_qx[k,l,e]	Numerical flux contribution to the $k^{th}$ basis of $\widehat{Q_x}$ at local edge $l$ of element $e$
EDGE_qy[k,l,e]	Numerical flux contribution to the $k^{th}$ basis of $\widehat{Q_y}$ at local edge $l$ of element $e$

Table 5. Table 5 : Estimated and experimental DRAM arithmetic intensities α \alpha for several kernels. Values for GFLOPs and Gbytes represent the total count for all elements/nodes/edges for the Harvey test case.

		Estimated		Experimental
Kernel	GFLOPs	GBytes	$\alpha$	$\alpha$
EDGE	7.95	1.6	4.97	2.7
AREA	4.09	2.2	1.85	1.81
GATHER	0.24	1.92	0.125	0.09

Table 6. Table 6 : Occupancy and instruction breakdown of the integration kernels. Note that these are warp instructions , i.e. each one operates on 32 scalars.

Kernel	Registers per thread	Occupancy (%)	FP32 Inst.	INT Instr.
EDGE	249	9.4	197 $\times {10}^6$	$118\times {10}^6$
AREA	163	14.9	$110\times {10}^6$	$78\times {10}^6$
GATHER	72	40.6	$5.6\times {10}^6$	$110\times {10}^6$

Table 7. Table 7 : Line count of OpenACC directives.

Directive	Description	# of lines
!$acc parallel	Loop parallelism	29
!$acc routine	Convert of CPU to GPU routines	15
!$acc loop seq	Serialize loops	27
!$acc atomic	Atomic updates	15
!$acc wait	Barrier	8
Continuation lines		4
Total		98

Equations

$$\begin{array}[]{ll}\text{Find }(\zeta,\mathbf{u})\text{ in }\Omega\text{ such that:}\qquad\qquad\\ \\ \displaystyle\frac{\partial\zeta}{\partial t}+\nabla\cdot(H{\mathbf{u}})=R,\\ \\ \displaystyle\frac{\partial(Hu_{x})}{\partial t}+\nabla\cdot\left(Hu_{x}^{2}+\frac{g}{2}(H^{2}-h_{b}^{2}),Hu_{x}u_{y}\right)-g\zeta\frac{\partial h_{b}}{\partial x}+\tau_{b}Hu_{x}=F_{x},&\\ \\ \displaystyle\frac{\partial(Hu_{y})}{\partial t}+\nabla\cdot\left(Hu_{x}u_{y},Hu_{y}^{2}+\frac{g}{2}(H^{2}-h_{b}^{2})\right)-g\zeta\frac{\partial h_{b}}{\partial y}+\tau_{b}Hu_{y}=F_{y},\end{array}$$

$$\mathbf{w}={\left[\begin{array}{ccc}\zeta ,& uH,& vH\end{array}\right]}^T,$$

$$\mathbf{s}={\left[\begin{array}{ccc}R,& g\zeta \frac{\partial h_b}{\partial x}+F_x,& g\zeta \frac{\partial h_b}{\partial y}+F_y\end{array}\right]}^T,$$

$$\mathbf{F}(\mathbf{w})=\left[\begin{array}{cc}uH& vH\\ H{u}^2+g(H^2-h_b^2)& Huv\\ Huv& H{v}^2+g(H^2-h_b^2)\end{array}\right],$$

$$\frac{\partial \mathbf{w}}{\partial t}+\nabla \cdot \mathbf{F}(\mathbf{w})=\mathbf{s},i=1,2,3.$$

$$\Omega =\bigcup _e{\Omega }_e$$

$${\left(\frac{\partial \mathbf{w}}{\partial t},\mathbf{v}\right)}_{{\Omega }_e}=-{\left(\nabla \mathbf{v},\mathbf{F}\right)}_{{\Omega }_e}-⟨\hat{\mathbf{F}}\cdot \mathbf{n},\mathbf{v}{⟩}_{\partial {\Omega }_e}+{\left(\mathbf{s},\mathbf{v}\right)}_{{\Omega }_e}\text{for each element }e,$$

$$\tilde{\mathbf{w}}=\sum _i\sum _j{\tilde{\mathbf{w}}}_{ij}{\varphi }_{ij},$$

$$\frac{d}{dt}(\tilde{\mathbf{w}})=L_h(\tilde{\mathbf{w}})\text{on each element }e,$$

$${\left(\nabla v,F\right)}_{{\Omega }_e}+{\left(s,v\right)}_{{\Omega }_e}$$

$$F_x,F_y$$

$$G_x,G_y$$

$$H_x,H_y$$

$$⟨{\hat{\mathbf{F}}}_i\cdot \mathbf{n},\mathbf{v}{⟩}_{\partial {\Omega }_e},i=1,2,3,$$

$$\alpha =\frac{\text{Floating point operations}}{\text{Bytes read from DRAM }+\text{Bytes written to DRAM }}.$$

$$\text{Attainable performance (FLOPs/s)}=\min (\text{Peak FLOPs/s},\alpha \times \text{Peak DRAM bandwidth}),$$

Full text

1 Introduction

One of the primary causes of flooding on the Texas and US Gulf and Atlantic coasts is storm surge: the inundation of coastal regions by water in the ocean that is driven inland by hurricanes or tropical storms. This phenomenon is extremely devastating to life and property of those nearshore as well as coastal ecosystems. In mitigating damage and guiding relief and emergency efforts, computer simulation of storm surge is critical in both operational forecasting as well as in hindcasting past storms for validation, engineering planning, and design studies. Some well-established and commonly used models include, for storm surge, the ADvanced CIRCulation (ADCIRC) model [Pringle et al. (2020), Luettich et al. (1992)], the Sea, Lake, and Overland Surges from Hurricanes (SLOSH) [Jelesnianski (1992)], and Delft3D [Veeramony et al. (2017)]; they have all been extensively validated for past hurricane events.

In this paper, the focus is on a model closely related to ADCIRC, the Discontinuous Galerkin Shallow Water Equation Model (DG-SWEM) [Kubatko et al. (2006), Dawson et al. (2011), Wichitrnithed et al. (2024)] and aims to improve the main criticism of this model, i.e., the computational cost. The formulation of ADCIRC is based on the two-dimensional shallow water equations (SWE) which are obtained by depth-averaging the incompressible Navier-Stokes equations. The advantage of using finite elements is the ability to model complex coastal domains through unstructured meshes. At the same time, a straightforward usage of the continuous Bubnov-Galerkin (CG) method results in spurious oscillations for advection-dominated flows. ADCIRC handles this by using a reformulation of the continuity equation in the SWE called the Generalized Wave Continuity Equation (GWCE), but as a consequence of this and the CG methodology, local mass conservation is no longer guaranteed. Furthermore, in practice, it still necessary to turn the advective terms in the GWCE off for some cases to achieve numerical stability.

One critical type of flooding that presents difficulties for ADCIRC is compound flooding which arises when multiple water sources such as storm surge and rainfall with resulting river discharge occur simultaneously. In practice, ADCIRC has encountered instability when rainfall is added to the source term of the GWCE and its subsequent interaction with wetting and drying scheme. DG-SWEM has recently [Wichitrnithed et al. (2024)] been shown to be able to address many of these problems by using the discontinuous Galerkin (DG) finite element formulation on the primitive SWE [Dawson et al. (2011), Dawson and Proft (2002)]. The DG method can be viewed as a synthesis of the Galerkin method and the finite volume method, whereby the solution is approximated using polynomial basis functions local to each finite element (cell) which guarantees local conservation. Solutions need not have a single value at each node (vertex) and coupling between adjacent elements is done through a numerical flux. It has been widely applied to hyperbolic conservation laws in various domains such as electromagnetics, acoustics, plasma physics, and shallow water equations. We refer to [Hesthaven and Warburton (2008), Cockburn (1998)] for the history and details of the DG method.

A primary disadvantage of the DG method which has confined DG-SWEM to mainly research and not operational forecasting like ADCIRC is the increased computational cost associated with the computations required from the higher number of global degrees of freedom compared to the continuous formulation. At the same time, the localized nature of computing the solution for each element independently generates a high degree of data parallelism, and several GPU implementations of the DG method have been using different approaches. Early implementations using NVIDIA’s CUDA include [Klöckner et al. (2009), Siebenborn et al. (2012), Fuhry (2013), Brodtkorb et al. (2012)]. More recent implementations using different application programming interfaces (API) can be seen in [Abdi et al. (2019), Kirby and Mavriplis (2020)] which use OCCA, and [Xia et al. (2014), Dai et al. (2022)] which use OpenACC. Lastly, we mention the recent open-source FEniCSx SWE solverby Dawson et al. in [Dawson et al. (2024)] and GPU accelerated by Pachev in [Pachev (2024)], where implicit schemes are considered.

In Section 2, we introduce the governing model, i.e., the SWE and our DG discretization. Next, in Sections 3 and 4, we discuss the related work and the hardware setup we utilize. Section 5 details our implementation. In Section 6, we discuss multi-node performance results when running the code on a large Hurricane dataset. Finally, in Section  7 we draw conclusions and discuss potential future research directions.

2 Model

2.1 Shallow water equations

The governing model in DG-SWEM is the two-dimensional shallow water equations which consists of mass and momentum conservation in each direction. The problem can be visualized as in Figure 1 and formulated as follows [Tan (1992)]:

[TABLE]

where:

$\zeta$

is the free surface elevation (positive upwards from geoid),

$h_{b}$

the bathymetry (positive downwards from the geoid),

$H$

the total water column,

$\mathbf{u}=\{u_{x},u_{y}\}^{\text{T}}$

the depth-averaged velocity field,

$\tau_{b}$

the bottom friction factor,

$R$

the mass source term such as rain,

$F_{x},F_{y}$

the relevant sources which induce flow, including, e.g., Coriolis force, tidal potential forces, wind stresses, and wave radiation stresses, and

$\Omega$

the computational domain, e.g., the coastal ocean, and its boundary $\Gamma$ is identified by three distinctive sections $\Gamma=\Gamma_{ocean}\cup\Gamma_{land}\cup\Gamma_{river}$.

2.2 The Discontinuous Galerkin Formulation

The particular DG formulation and solution method used in this work is based on previous works [Kubatko et al. (2006), Dawson et al. (2011)]. Hence, we do not include a comprehensive introduction to the details of the formulation but rather a brief overview to make this presentation sufficiently self-contained. This description is an updated version from our recent work [Wichitrnithed et al. (2024)].

First, we define the vector of state variables to be solved, $\mathbf{w}=\{\zeta,Hu_{x},Hu_{y}\}^{\text{T}}$, the source vector $\mathbf{s}$ and the flux matrix $\mathbf{F}(\mathbf{w})$ as:

[TABLE]

[TABLE]

[TABLE]

we can write these three equations a conservative form

[TABLE]

As is standard in the discretization of (2) with a finite element formulation, we partition the domain into (triangular) non-overlapping elements,

[TABLE]

and the corresponding weak formulation is obtained by multiplying each equation by a test function $\mathbf{v}$ and integrating over each element and subsequently performing an integrating by parts to the term with the spatial derivative:

[TABLE]

where we use standard inner product notation for integration over elements, i.e., $\left(u,v\right)_{\Omega_{e}}=\int_{\Omega_{e}}u\,v$ and element edges $\langle u,v\rangle_{\partial\Omega_{e}}=\int_{\partial\Omega_{e}}u\,v$. $\hat{\mathbf{F}}$ represents the numerical flux, e.g., Local Lax-Friedrichs or Roe flux. The role of the numerical flux is to couple adjacent, discontinuous elements and correctly capture the flow characteristics across the edges; this is necessary for the solution to be stable [Hesthaven and Warburton (2008), Chapter 2]. Thus, aside from calculating $\hat{\mathbf{F}}$ (which depends on two neighboring elements), the computation can proceed on each element independently.

In accordance with the discontinuous Galerkin method, we next approximate the solution $\mathbf{w}$ in space as $\widetilde{\mathbf{w}}$ using a set of basis functions. We opt to use the orthogonal Dubiner basis [Dubiner (1991)] $\phi_{ij}$ as:

[TABLE]

where $\widetilde{\mathbf{w}}_{ij}$ are the modal degrees of freedom. The indices $i,j$ indicate the order of the polynomial: $\phi_{ij}$ is of order $i+j$. This basis is chosen as it produces a diagonal mass matrix, and higher-order approximations can be obtained by simply adding more terms. Note that while DG-SWEM supports higher-order polynomial orders, this work implements only the linear case, i.e. the first three basis functions are used. In the context of simulating storm surges and similar to ADCIRC, we found that using linear (or even constant) approximation is sufficient. Substituting (4) into the weak formulation (3) above, and requiring that it holds for all basis functions $\mathbf{v}\in\{\phi_{00},\phi_{01},\phi_{10}\}$ reduces the problem to a system of ODEs in the form:

[TABLE]

where the terms from (3) have been collected into the RHS $L_{h}$. Numerical integration of each term is performed using Gauss-Lobatto quadrature.

To integrate in time, we use an optimized version of the explicit Strong Stability Preserving Runge-Kutta scheme (SSPRK) based on [Shu and Osher (1988)] to discretize (5) in time and compute the difference approximation $\mathbf{w}_{h}^{n+1}$ of $\widetilde{\mathbf{w}}$, at the next timestep. This scheme has been specifically optimized with respect to DG spatial discretizations to enhance the CFL stability regions and therefore increase the possible time step size for a given order and stage [Kubatko et al. (2014)]. This is also combined with an optional slope limiter to suppress oscillations. Finally, to output the state variables, we convert the modal representation to a nodal representation by averaging the neighboring elements of each node.

In terms of execution, DG-SWEM supports parallel execution on CPUs through domain decomposition with METIS [Karypis and Kumar (1997)] using the Message-Passing Interface (MPI). As the scheme is explicit in both space and time, most computations are local to each element, except the calculation of the numerical flux $\hat{\textbf{F}}$ that is shared between two neighboring elements. We use the standard ghost-cell technique to transfer state variables of elements at the boundary of each domain. An evaluation of the parallel efficiency of the code is reported in [Kubatko et al. (2009)].

3 GPU Programming approach

As an alternative to purely distributed parallelism, graphics processing units (GPUs) can be viewed as a multithreaded vector processor [Hennessy and Patterson (2018)], where each unit executes the same instructions on multiple data (SIMD), but individual lanes can take different paths. To take advantage of a large number of such processors, the algorithm must exhibit a high degree of data parallelism, where multiple elements can be processed independently at the same time. This is particularly applicable to explicit DG methods due to the inherently local computation of the elements.

In the past, as GPUs have originally been designed for graphics programming, the available application programming interface (API) such as DirectX and OpenGL imposed limitations of using GPUs for general-purpose computing (GPGPU) [Kirk and Hwu (2010)]. In 2007, NVIDIA developed CUDA [NVIDIA (2024)] which provides a more generic API to GPU instructions. This significantly extended the capabilities of NVIDIA GPUs to scientific and engineering codes: the programmer gained control of individual GPU threads, synchronization, and shared memory. Using the NVIDIA HPC compiler stack (a successor of PGI compilers), CUDA is supported in C, C++, and Fortran [NVIDIA (2020)]; CUDA code can be written in the same source files and compiled with appropriate flags. While the low-level nature of CUDA allows more opportunities for control and optimization by the programmer, it typically results in having to maintain a separate code base from the CPU code. In the past, a framework has been created to generate GPU code from DG-SWEM through CUDA C++ [DuChene et al. (2011)]. The workflow consisted of scripts that transform Fortran subroutines to C++ and handle variable declaration and allocation. The work noted that other alternatives at the time, such as using the PGI accelerator or CUDA-Fortran, were lacking in features or accessibility.

A more modern approach to GPU programming is to use OpenACC [OpenACC-Standard.org (2021)] or OpenMP [Openmp.org (2021)] to write directives instead of manually managing individual GPU threads. This approach allows the programmer to insert high-level statements that generate GPU instructions and appear as comments in the source code (similar to OpenMP loop constructs, for instance). This can sometimes result in a single code base — only the compilation flags such as -acc need to be turned on to enable the GPU code generation. Of course, the programmer is still required to sufficiently expose the parallelism in the code, which often results in modifications of the CPU code as well. We choose to use OpenACC in this work as it currently has better support on the NVHPC compilers. In this approach, we aim to mimimize integration and maintenance efforts, as well as to take advantage of the new NVIDIA Grace-Hopper architecture [Singer (2024)] which provides a more efficient CPU-GPU data transfer.

4 Experimental setup

The primary hardware used in this work is the recent Vista Supercomputer [Singer (2024)] at the Texas Advanced Computing Center (TACC). Vista uses the ARM-based NVIDIA GH200 Grace Hopper architecture. Each GPU node consists of an NVIDIA H100 GPU connected to an NVIDIA Grace CPU directly through high-bandwidth NVLink-C2C interconnect. The Grace CPU contains 72 Arm Neoverse V2 cores and 120 GiB of LPDDR5X memory, while the Hopper GPU contains 96 GB of HBM3 memory. Each CPU node in the CPU-only subsystem contains 144 Grace cores on two sockets with a total memory of 237 GB LPDDR5. Different nodes are connected in a fat-tree topology, with GG nodes through an InfiniBand 200 GB/s and GH nodes an InfiniBand 400 GB/s. Specifications for these nodes are listed in Table 1 and 2.

The NVLink interconnect provides a faster data transfer between CPU and GPU. Previous generations of GPUs connect to CPUs through PCIe and trigger data migration on page faults. GH200 uses a hardware-based approach to achieve the same functionality by sharing the same page table and virtual address space which results in less page faults. In other words, CPU and GPU memories are now treated as NUMA nodes [Fusco et al. (2024)].

The original CPU code was compiled using the Intel oneAPI compiler stack. This work uses the NVIDIA SDK compiler stack 24.7 (nvfortran, nvc, nvc++) for both CPU and GPU versions. Compilation flags for the CPU and GPU codes are listed in Table 3 .

5 Implementation

5.1 Unified Memory

A significant part of porting code to GPUs is managing data movement. As a CPU (host) and a GPU (device) have separate physical memories, data must be transferred back and forth to perform the respective calculations. Traditionally, even if some arrays are only used in device calculations, they still need to be initialized or read from files through the host, and we would have to declare device copies of those arrays. For example, in CUDA Fortran, an existing array in the code

real, allocatable :: A(:)

will need an additional copy on the device

real, allocatable, device :: A_d(:)

and copied to and from the device as necessary:

A_d = A

! Use A_d on the GPU

A = A_d

To avoid having to do this for each array in the code, we exploit Unified Memory which was introduced in CUDA 6.0. This allows GPU kernels to access data on the host without needing to explicitly copy the data to device. Data migration is automatically performed once a page fault is triggered on the device. In the past, this requires a special allocation, e.g. cudaMallocManaged in CUDA C. For newer systems that support Heterogenous Memory Management (HMM), all device data declarations can be completely eliminated using the compilation flag -gpu=unified. This significantly simplifies programming effort and enhances maintainability. In both approaches, a synchronization statement like !$acc wait is needed for correctness before the CPU can access the computed data, since GPU kernels immediately return control to the host.

In Table 4 we summarize the important data arrays, some of which are referred to in the kernel descriptions in the next section.

5.2 GPU kernels

This section outlines the main subroutines and how their parallelization. For the integration kernels, we follow the one-thread-per-element and one-thread-per-edge approach as in [Fuhry (2013)].

5.2.1 Area integration

The area integration routine (Algorithm 1) is naturally data-parallel, as computations are done locally to each element $\Omega_{e}$. At a given Runge-Kutta stage irk, it computes the terms

[TABLE]

from (3) and accumulates them into the RHS vectors. Recall from Eq 2 that the flux terms (explicitly written in terms of the state variables) are

[TABLE]

and these are represented by the procedures Continuity_Flux, X_Flux, and Y_Flux in the code. The source terms $S_{R},S_{x},S_{y}$ are calculated in the same way but retrieve additional variables such as current rain intensity and wind speed.

There are technically three levels of loop parallelism: looping over the elements, looping over the quadrature points, and looping over the degrees of freedom. Here, the element-level loop is accelerated through the gang-vector parallelism over the number of elements $N_{e}$ (line 2), and the inner loops over quadrature points and degrees of freedom are enforced to be serial using !$acc seq. This needs to be explicitly inserted, as by default the inner loops will also be parallelized, and in this case the trip counts is too short to be profitable. We also add the private(...) clause to keep any local arrays private to each thread. Importantly, the async clause is added to eliminate kernel launch overhead, and the default(present) the overhead of checking if data is already on the device. This same structure is applied for the edge integration kernel. We have observed that using these clauses speeds up the program by almost 1.5$\times$.

5.2.2 Edge integration

Edge integration routines compute the line integrals

[TABLE]

and accumulate them to the RHS. Similar to the volume integration kernel, a natural parallelism is in the number of edges $N_{d}$. There are multiple ways to add contributions from inter-element fluxes. One is to let each edge compute the numerical flux between the two elements, and then atomically add back the contributions to each. Atomics are necessary since two edges can share the same element.

Another approach is to first store the flux contributions to each element separately for each of its edges, i.e., we create an additional array elem_edge(E, i) which stores the contribution at edges $i=1,2,3$ of element $E$. Subsequently, in each element we then sum the three contributions locally without needing atomics. We use this approach in our work as it leads to reproducible, deterministic results due to the fixed order of summation. The pseudocode for these two steps are shown in Algorithm 2 and Algorithm 3. Note that we have abstracted away the wetting/drying branching for brevity.

In Algorithm 2, each GPU thread is mapped to an interior edge, again through gang-vector parallelism (line 2). In line 4-9, the edge number is used to obtain the IDs of the two neighboring elements as well as their local edge numbers (1-3), with the interior side denoted with $+$ and exterior $-$. Line integration on the edge is done by looping through the edge quadrature points. First, at each point we compute the state variables (Line 14-23), use them to compute the numerical flux. Then, we add this flux contributions to both sides of the edge, indexed using the local edge number of each element (line 35-41). The subsequent elemental flux summation is shown in Algorithm 3, where each thread is mapped to an element and sums the contributions from its 3 edges. Note that for brevity, we have opted for the atomic approach for other less common boundary integration routines.

5.2.3 Other kernels

Other kernels like slope limiting (SLOPELIM) and wet/dry calculations (WETDRY) are similarly parallelized at the element-level loop, with small inner loops serialized. Wind forcing and pressure calculations (WIND) are done in a series of loops, each operating on the mesh nodes independently. We simply parallelize each loop with !$acc parallel loop.

5.2.4 Data transfers between host and device

At every fixed number of time steps specified by the user, the program writes the state variables to ASCII files. Thus data will have to be copied to the CPU. While it is possible to completely rely on HMM to initiate the transfers, they can incur large overhead as there are enough access counts to trigger page migration. To avoid this, we use cudaMemPrefetchAsync to transfer data at once, as shown in Listing 1.

5.2.5 Memory usage of automatic arrays

We sometimes observe runtime memory error when parallelizing loop with many automatic (private) arrays per thread. From our observation, this seems to be due to the allocation on these arrays on the GPU stack, a different behavior from using the CUDA API. We can avoid this by limiting the number of gangs to a fixed number, e.g.

!$acc parallel loop gang vector num_gangs(1024) private(<long list of arrays>).

as opposed to allowing an automatic scheduling of gangs.

5.3 Multiple GPUs

To exchange data among different domains when using multiple processors, DG-SWEM uses MPI persistent communication. At the start of the simulation, calls to MPI_SEND_INIT() and MPI_RECV_INIT() are made to specify the location of data to exchange to each neighbor. Then, as shown in Listing 2, whenever we need to exchange data, we call MPI_STARTALL() to initiate the communication in a non-blocking manner. At the end, we wait until all the messages have arrived. The packing of the elements to send is done through the subroutine packElems(), and the unpacking similarly through unpackElems().

The existing MPI routines need no further modification in the GPU code, as they accept unified memory addresses. However, we observe data being transferred to host and back to device during MPI communication. This appears to be a limitation of using unified memory with CUDA-aware MPI which prevents direct GPU-to-GPU data transfer [Chenbr (2025)]. The packing and unpacking have been parallelized to take advantage of the higher memory bandwidth of the GPU, as shown in Listing 3.

6 Results

6.1 Test Case

To measure performance improvements, we run the solver on a dataset for Hurricane Harvey (2017) which is representative of a typical simulation scenario. The mesh is shown in Figure 2 which contains 6,675,517 elements and 3,352,598 nodes with a maximum resolution of 20 meters. We impose a tidal (periodic) elevation boundary on the right side and use Harvey Best Track data from the National Hurricane Center HURDAT2 database. The timestep size is set to be $\Delta t=0.3$s, and computation is done with single-precision.

Figure 3 shows the performance comparison between the GPU implementation and the CPU implementation from 1 to 8 nodes, where we run the simulation for 28,800 time steps. As described in Section 4, each GPU node consists of a single NVIDIA H200, while each CPU node contains 144 NVIDIA Grace cores. Here, ”performance” is measured by the rate of work done, defined to be the number of degrees of freedom computed per second. In a node-by-node comparison, we observe a speedup of GPU code over CPU of around $4\times$ to $2.55\times$. Figure 4 further shows the breakdown of time spent in each kernel. In this case, we only consider the total time spent in the timestepping loop (i.e. exclude I/O or startup time) and run the simulation for 2,880 time steps. Here, we see that parallel communication time in the GPU case increases with more GPUs to around 6 s at 8 nodes, resulting in a scaling plateau. For the same amount of nodes in the CPU case case, the communication time appears to stabilize to around 18 s.

6.2 Performance Evaluation

To further evaluate performance of the time-consuming kernels, we first calculate the theoretical arithmetic intensity $\alpha$ at the DRAM level. This is the ratio of floating point operations (FLOPs) to total memory traffic in bytes from/to DRAM:

[TABLE]

We count each arithmetic operation like addition and division counts as one FLOP, including special operations like square root and power. In estimating $\alpha$ from the source code (by manually counting), we also assume that everything is loaded once from memory for simplicity, which is generally not valid for large enough datasets and when accesses are not fully coalesced. Thus we expect the estimated values to be higher than the true values. The estimated arithmetic intensities for the integration kernels are shown in Table 5, along with their experimental values obtained from the NVIDIA Nsight Compute profiler. We focus on these two kernels as they are representative of the DG method.

It can be seen that the arithmetic intensity of AREA is very close to the estimated value. This is expected, as most data are not shared between different elements, and are loaded in a (mostly) coalesced manner. A similar reasoning applies for GATHER. On the other hand, EDGE has a markedly lower intensity than estimated, a little over half. This can be attributed to the different memory access patterns. While area integration loops over elemental arrays sequentially, edge integration requires each edge to fetch data from neighboring elements which may not lie close together in memory. This poor spatial locality is a well-known issue in using unstructured meshes [Hadade et al. (2020)] and in our case results in uncoalesced accesses and therefore more data traffic.

We further evaluate the performance using the empirical Roofline which prescribes a performance bound through the relation:

[TABLE]

which makes some assumptions, e.g., latencies are ignored and computation perfectly overlaps with data transfer [Hager and Wellein (2010)]. In kernels with low $\alpha$ (like most in our code), performance will, at best, reach only a fraction of the peak FLOP/s of the GPU.

To plot the empirical roofline, we use the procedure described in [Yang (2020)]; this is shown in Figure 5. The AREA, EDGE, and GATHER kernels achieve compute throughputs of 1.9 TFLOPs/s, 1.3 TFLOPs/s, and 0.24 TFLOPs/s, respectively. Other relevant performance characteristics from Nsight Compute are shown in Table 6. A few observations can be made. First, the INT instruction count shows that a large portion of the computation are in fact array index calculations. While this is largely unavoidable, they do lower the measured FLOP rate, as we cannot assume they overlap with FP32 instructions. Second, the achieved occupancy for the first two kernels are quite low due to the high register usage, and this can also limit performance. While it is possible to tune register usage through the -maxregcount compiler flag, we found that the penalty from spilling largely offsets the speedup from higher occupancy. We also note that using linear ($p=1$) elements is typically insufficient to reach close to peak performance of the GPU [Fuhry (2013), Klöckner et al. (2009), Abdi et al. (2019)].

6.3 Lines of code changed

To get a rough estimate of the level of programming effort, we compare the total lines of code before and after incorporating OpenACC (we ignore blank lines). Lines of code increase by 147, from 48,674 to 48,821; this change includes the addition OpenACC directives and other modifications such as the flux-gathering routine and MPI packing/unpacking. Table 7 shows the breakdown of the specific OpenACC directives used.

7 Conclusion

We have ported DG-SWEM to NVIDIA GPUs using OpenACC Fortran in a basic manner by exploiting the data parallelism inherent in the explicit-time DG method through element and edge computations . We tested the code on multiple nodes on the NVIDIA Grace Hopper system using a large flooding test case and observed significant speedups over the original CPU implementation. When combined with the CUDA Unified Memory System, the porting requires very few code modifications and can be maintained in the same code base; this is a marked improvement from our previous porting framework [DuChene et al. (2011)]. Finally, we note that there are potential performance issues that can be further investigated on top of this basic parallelization. First, a more detailed benchmark on multi-node scaling can be done to observe the code’s relative performance, and whether further optimizations like overlapping computation and communication will be beneficial. Second, we can explore alternative ways to structure edge integration to exploit GPU shared memory and improve data locality among neighboring edges.

8 Data Availability

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

• •

Flood datasets used

• •

DG-SWEM code - https://github.com/foci/dgswem

9 Acknowledgements

This work has been supported by the United States National Science Foundation NSF PREEVENTS Track 2 Program, under NSF Grant Numbers 1855047. This material is also based on work supported by the US Department of Homeland Security under Grant No.2015-ST-061-ND0001-01. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the US Department of Homeland Security.

The authors also would like to gratefully acknowledge the use of the ”ADCIRC” allocation on the Vista supercomputer, and “ADCIRC”, “DMS23001”, and “DMS21031” allocations on the Frontera supercomputer at the Texas Advanced Computing Center at the University of Texas at Austin.