arXiv CSopen access

An efficient multi-GPU implementation for the Discontinuous Galerkin ocean model SLIM

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

An efficient multi-GPU implementation for the Discontinuous Galerkin ocean model SLIM M. De Le Court1 , V. Legat1 , A. P. Ishimwe3 , C. Scherpereel2 , E. Hanert1,2 , and J. Lambrechts1 1 Institute of Mechanics, Materials and Civil Engineering, UCLouvain, Louvain-la-Neuve, Belgium 2 Earth and Life Institute, UCLouvain, Louvain-la-Neuve, Belgium

arXiv:2605.16082v1 [cs.DC] 15 May 2026

3 Department of Ecoscience, Aarhus University, Roskilde, Denmark

Abstract Unstructured-mesh ocean models are increasingly used for coastal applications due to their ability to represent complex geometries and apply local grid refinement where needed. However, their broader use has been hindered by their high computational cost, particularly for models based on the Discontinuous Galerkin finite element (DG-FE) method, which involves significantly more degrees of freedom than traditional finite volume or continuous finite element approaches. The rapid emergence of GPU-based high-performance computing architectures now offers a pathway to address this limitation, as DG-FE formulations are inherently well suited to massively parallel, element-wise computations. Here, we present a full 3D DG-FE ocean model implementation optimized for both single- and multi-GPU systems, with support for both NVIDIA and AMD architectures. We detail the computational strategies employed to achieve high performance, including memory layout optimization, kernel-level parallelization, and matrix-free solvers for key vertical processes. Benchmark results demonstrate that a single HPC-grade GPU (e.g. NVIDIA A100) delivers performance equivalent to approximately 1500 CPU cores, while replacing a 128-core CPU node with a 4×A100 GPU node yields a speedup of around 50×. Weak-scaling efficiency is maintained up to 1024 GPUs. We further demonstrate the model’s capabilities on a real-world application in the Great Barrier Reef, achieving a spatial resolution five times finer than the most accurate existing model while maintaining a physical-to-numerical time ratio of 100. These results highlight how GPU-accelerated DG-FE methods can dramatically advance the capabilities of unstructured-mesh ocean modeling, enabling ultra-high-resolution coastal simulations that were previously infeasible.

Contents 1 Model 1.1 Primitive equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Split momentum equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 4 5

2 Efficient GPU implementation of the model 6 2.1 Memory layout and thread assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Cell structure for linear system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Going to and from the cell layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 3D hydrodynamical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Matrix-free solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Fully-assembled column solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Data management . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3 Multi-GPU 3.1 Structure of a compute kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 MPI communications for the 3D kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Challenges with the short 2D kernels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15 15 15 15

4 Performance and scaling of the model 18 4.1 Single-GPU performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.2 Multi-GPU scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5 Application to the Great Barrier Reef

6 Conclusion

Introduction Unstructured-mesh ocean models, once considered a niche approach, have progressively gained prominence and now represent a mainstream modelling approach, particularly for coastal ocean studies. Their growing adoption in coastal applications is largely due to their geometrical flexibility, allowing them to accurately represent intricate bathymetric and topographic features through local mesh refinement, as well as to resolve multiscale processes along the land-sea continuum (Deleersnijder et al., 2010; Lermusiaux et al., 2013). Unstructured-mesh ocean models can be based on different numerical methods, including the finite volume (FV) method (Chen et al., 2003; Danilov et al., 2017; Korn et al., 2022; Zhang et al., 2016), the continuous Galerkin finite element (CG-FE) method (Umgiesser et al., 2004; Wang et al., 2014; Westerink et al., 2008) and the Discontinuous Galerkin finite element (DG-FE) method (Kärnä et al., 2013, 2018). However, despite their advantages, unstructured-mesh models are often limited by a higher computational cost per degree of freedom than structured-grid models. As a result, they remain less commonly employed for global ocean simulations as compared to structured-grid models (Danilov et al., 2017). The higher computational demand of unstructured-mesh ocean models, and particularly DG-FE models, necessitates the implementation of cost-reduction strategies to help make them viable for largescale simulations. The DG-FE formulation combines inter-element fluxes (typical of FV methods) with intra-element calculations (typical of CG-FE methods), hence implying more calculations than these two approaches. Additionally, the data organization on an unstructured mesh with irregular connections between elements requires non-trivial indexing and neighbour-finding algorithms due to its non-systematic structure. Consequently, neighbour relationships lack implicit connectivity information and must be explicitly stored and managed. This explicit storage necessitates more sophisticated data handling and memory management methods, ultimately leading to higher computational costs. On the other hand, various strategies have been proposed to address these computational challenges, including local mesh refinement (Danilov & Wang, 2015), adaptive time stepping (Dawson et al., 2013; Seny et al., 2014), or mode-splitting procedures (Shchepetkin & McWilliams, 2005). The first two strategies complement each other by focusing computing resources where they are most needed and employing smaller time steps exclusively within regions of high mesh resolution. In contrast, the latter optimizes the use of computational resources by splitting the treatment of barotropic and baroclinic processes, using small explicit

timesteps for fast processes while employing a single, larger timestep for slower processes. However, these strategies alone are not sufficient to make the routine application of unstructured-mesh models straightforward, as their inherent computational burden remains significant, particularly for high-resolution or long-duration simulations. Over the last 10-15 years, significant improvements in model efficiency have also been facilitated by the evolution of high-performance computing (HPC) infrastructures, particularly the widespread transition from CPU-based to GPU-based architectures. Graphics Processing Units (GPUs), characterized by their substantial computational power per chip, have notably accelerated HPC capabilities, reaching exascale performances (Chang et al., 2023). This GPU revolution, supported by the development of programming systems such as CUDA, HIP, or OpenCL (Dally et al., 2021), is clearly illustrated by the fact that GPU-based supercomputers now dominate the list of the world’s fastest HPC systems. To date, nine of the top ten supercomputers in the TOP500 list (Strohmaier et al., 2025) are accelerated by either INTEL, AMD or NVIDIA GPUs. However, leveraging GPU performance necessitates a complete redesign of computational codes, representing a non-trivial challenge, especially for mature models. Achieving performance, portability and code maintainability have now become three major conflicting objectives preventing well-established modelling frameworks from evolving as fast as their computing clusters. Despite this rapid hardware evolution, many of the most widely used structured-grid ocean models remain predominantly CPU-based. Systems such as NEMO (Madec et al., 2022), HYCOM (Bleck, 2002), ROMS (Shchepetkin & McWilliams, 2005), and MITgcm (Marshall et al., 1997) are still implemented mainly in Fortran and parallelized with MPI or OpenMP. While partial GPU adaptations exist such as for ROMS (Panzer et al., 2013), these are generally limited to specific components and do not deliver full model acceleration. Recent progress in global circulation modelling, however, illustrates what is possible when GPU architectures are fully embraced. LICOM3-HIP (Wei et al., 2024) demonstrates efficient multi-GPU scaling for eddy-resolving global simulations, while Veros (Häfner et al., 2021) shows how a structured finite-difference model rewritten in Python/JAX can achieve strong performance gains through just-intime compilation. Oceananigans.jl (Wagner et al., 2025), developed natively in Julia using kernel-fusion techniques that make highly effective use of GPU memory bandwidth, represents a new generation of GPU-first models, reaching 10 SYPD at 8 km resolution on 64 GPUs (Silvestri et al., 2024). These examples underscore how models purpose-built for GPUs can outperform legacy codes even at global scales. Regional ocean modelling still has comparatively few GPU-ready options. In this paper, we present the latest iteration of the three-dimensional unstructured-mesh coastal ocean model SLIM, fully redesigned and optimized to efficiently utilize multiple GPUs. A light abstraction layer allows switching the target architecture between CPUs or NVIDIA or AMD GPUs using the same source code. The use of the Discontinuous Galerkin method in the SLIM model is especially interesting for the GPU implementation. Compared to other unstructured approaches, the DG method enjoys a high data locality and arithmetic intensity. This lets it use more of the GPUs’ computational power. Its block-based discretization is well suited to their high throughput and parallel capabilities. SLIM runs on triangular meshes extruded vertically into columns of prisms. This vertical structure is excellent for the computational efficiency of the model as the penalty of horizontally unstructured accesses can be significantly mitigated. The other major benefit of vertical extrusion is that the implicit treatment of vertical quantities only couples together nodes of the same column of prisms. Therefore, all columns can be processed in parallel, taking full advantage of the GPU. The first section describes the model’s DG-FE formulation. Section 2 details the GPU-focused implementation on a single GPU. Section 3 presents the multi-GPU strategy. Section 4 evaluates performance through comprehensive benchmarking. Finally, we demonstrate the model’s scientific capability using a realistic case study of the Great Barrier Reef.

Model

SLIM is a regional 2D and 3D ocean model. It uses the Discontinuous Galerkin (DG) method with linear Lagrange basis functions. 2D equations are discretized on an unstructured triangular mesh (Figure 1a), while 3D equations use a prismatic mesh made by extruding the 2D mesh vertically (Figure 1b). SLIM 3D solves the hydrostatic Boussinesq equations. This section introduces the numerical formulation of SLIM, which serves as the foundation for the optimizations described in the remainder of the paper.

(a) 2D mesh around the Whitsundays islands.

(b) Simplified 3D mesh.

Figure 1: Examples of meshes used in SLIM. The unstructured 2D triangular mesh is extruded vertically in columns of prisms.

1.1

Primitive equations

The primitive hydrostatic Boussinesq equations being solved are ∂z w = − ∇ h · u

(1)

∂t H = ∂t η = − ∇h · Q + s

(2)

∂t u + ∇h · (u ⊗ u) + ∂z (wu) =∇ · (ν∇u) − f ez × u −

1 ∇h p ρ0

(3)

with the following boundary conditions on the surface and bottom of the ocean: w − u · ∇h η = ∂t η − s w − u · ∇h b = 0

at the ocean surface,

(4)

on the ocean floor.

(5)

In the equations above, H denotes the height the water column, b the bathymetry and η the deviation from the reference level. u = [u, v] is the 3D horizontal velocity, and w the vertical velocity. Q is the horizontal transport defined as Z η

u dx .

The viscosity tensor is written as ν, f is the Coriolis parameter and ez is the vertical unit vector. The source term s represents rain (s > 0) or evaporation (s < 0). Additionally, temperature and salinity are modelled as tracers, following the advection-diffusion equation ∂t T + ∇h · (uT ) + ∂z (wT ) = ∇ · (κ∇T )

(6)

for a generic tracer T and associated diffusivity tensor κ. The parameterization of the eddy viscosity and diffusivity tensors is split between the horizontal and vertical directions. Vertically, the 2-equation GLS turbulence closure model (Umlauf & Burchard, 2003) discretized as in (Kärnä, 2020) is used for the eddy coefficients. Horizontally, the viscosity is parameterized by the Smagorinsky model (Smagorinsky, 1963), and the diffusivity uses the Okubo model (Okubo, 1971).

1.2

Split momentum equation

Under the hydrostatic assumption and accounting for the atmospheric pressure patm , the pressure p can be written in terms of the density ρ(x) = ρ0 + ρ′ (x) integrated over the water column. Z η Z η p(z) = patm + g ρ(z̃) dz̃ = patm + gρ0 (η − z) + g ρ′ (z̃) dz̃ z

Using this separation, we express the pressure gradient as 1 1 1 ∇h p = ∇h patm + g∇h η + r ρ0 ρ0 ρ0

(7)

with the horizontal pressure gradient r obtained by vertical integration  r = gρ′ (η)∇h η on Γs

(8)

∂ r = g∇ ρ′ z

from top to bottom. The momentum equation is thus split as S

}|2 { z r 1 ∂t u + ∇h · (u ⊗ u) + ∂z (wu) = ∇ · (ν∇u) − f ez × u − −g∇h η − ∇h patm ρ0 ρ0 Z η Z η H r ∂t Q = −gH∇h η − ∇h patm + ∇h · −u ⊗ u dz + ∇ · (ν∇u) − f ez × u − dz ρ0 ρ 0 b b {z } |

(9) (10)

where S2 are the 2D terms of the 3D momentum equation and S3 are the inegrated 3D terms in the 2D momentum equation. In combination with the 2D continuity equation (2), the term gH∇h η of the 2D momentum equation (10) creates the dynamics of fast-moving gravity waves. Like most ocean models (Kärnä et al., 2018; Madec et al., 2022; Shchepetkin & McWilliams, 2005), SLIM separates the dynamics into a fast barotropic (external) mode and a slower baroclinic (internal) mode. The barotropic mode represents depth-averaged gravity waves and is integrated on the two-dimensional mesh with a small time step. The baroclinic mode captures the fully three-dimensional internal motions and is advanced with a larger time step. A split-IMEX Runge-Kutta scheme (Ishimwe et al., 2023, 2025) couples the two: for each baroclinic step, several barotropic iterations are performed. This mode-splitting approach greatly reduces computational cost compared to applying the small barotropic time step to the full 3D system. The internal mode is further separated between its vertical and horizontal components. The Horizontal components of the advection, viscosity and diffusion for tracers are computed explicitly, but their vertical counterparts can be implicit. This mitigates the otherwise strict CFL condition arising from the vertical viscosity and diffusion. More details on the temporal and spatial discretizations can be found in the supporting information.

Efficient GPU implementation of the model

Having established the mathematical formulation and time-stepping structure, we now turn to the practical question of how these equations can be mapped efficiently to GPUs. We start by giving an overview of SLIM’s two-step split-IMEX scheme where the first step is vertically implicit, and the second is explicit. Owing to the mode-splitting approach, each of these steps of the internal mode is accompanied by several iterations of the external mode. The external mode use a three-step explicit Runge-Kutta method. Each step of the internal scheme consists of 5 main components, as illustrated in Figure 2a. Figure 2b presents a timeline of a full iteration, highlighting the 5 components within each step. Those will serve as the guiding structure for the core of this section where we discuss how each compnent is implemented.

2.1

Memory layout and thread assignment

In finite element computations, efficiently mapping elements to threads is a crucial consideration. In most kernels, we assign one thread per mesh element, rather than the other common approach of one thread per node. This choice increases the computational workload per thread, which is beneficial for low-order methods like ours, where the GPU often lacks sufficient computations to stay fully utilized, and the number of nodes per element is relatively small. Given this one-thread-per-element approach, we organize element data in a structure-of-arrays (SoA) format to optimize memory coalescence (Figure 3). This differs from the conventional array-of-structures format, where data for each element is stored contiguously. However, the structure-of-arrays approach comes with a drawback: reduced cache efficiency, particularly when dealing with an unstructured mesh. To address this, we reorder the 2D mesh following a Hilbert curve. The 3D mesh, in contrast, is less affected by this issue due to its vertical structure. Neighboring accesses either target the top or bottom elements (likely already in cache) or lateral elements, which are stored in contiguous memory. Consequently, as the number of vertical layers increases, the performance gap between structured and horizontally unstructured models narrows, a property we aim to leverage. 2.1.1

Cell structure for linear system

The primary drawback of the SoA data layout arises when solving the linear systems associated with each column of prisms. The sequential nature of solving a banded system would ideally suit a one-thread-percolumn approach rather than the one-thread-per-prism strategy. With the SoA layout, using one thread per column would result in large strides when accessing memory, leading to poor cache utilization and performance degradation. Distributing memory loads across multiple threads would add synchronization overheads and complexity. To address this, we maintain a single-thread-per-column approach but introduce what we call the cell layout (Figure 4), a data structure designed to optimize memory access patterns for linear solves while enabling efficient conversion to and from the SoA layout. The cell layout groups sets of prism columns (typically 128) and stores their data in a matrix (Figure 5). In the following, the term cell refers to such a group of 128 columns. Each column of the matrix correspond to one column of prisms. Each row corresponds to a data entry from the column of prisms. For example, a scalar field would be unrolled to 6Nlayers rows, 6 being the number of nodes per prism. When the number of layers varies, padding is added to match the deepest column of the cell. This ensures that a block of 128 threads solving the linear systems for 128 columns of prisms achieves perfectly coalesced memory access. The layout is akin to an array-of-structure-of-array (AoSoA), which is well suited for GPU computing. 2.1.2

Going to and from the cell layout

We assemble the system with one thread per element. A cell is divided into multiple GPU blocks of threads, each responsible for a subset of the cell’s columns. If a block of threads handles n columns, it writes to the cell in contiguous chunks of n values and reads from the SoA layout in chunks of ⌊128/n⌋ values. The choice of n is a trade-off. A larger n improves write efficiency to the cell and benefits kernels assembling larger matrices. Conversely, a smaller n enhances coalescence when reading data and is preferable when no linear system needs to be solved. In practice, n is determined per cell based on the maximum number of layers in the cell and an estimate of the average execution time across all kernels. In our case, prioritizing vertical size (⌊128/n⌋) over horizontal size (n) is generally more effective, as most kernels perform more reads than writes, and many do not require solving a system. With 128

1. First, a prediction for the 3D horizontal momentum fluxes is computed. This step is always explicit. The fluxes are then summed vertically and passed to the external 2D mode as a source term.

3D Horizontal

2. The 2D external mode is advanced in time with many RK iterations. It is responsible for the free surface movement and providing a correction to the horizontal transport.

3. The turbulent variables are evolved in time before the rest of the simulation. This enables the last two components to use the vertical viscosity and diffusivity coefficients at the end of the step.

Turbulence

1 4. The horizontal momentum fluxes are recomputed to include the 2D correction to the transport, as well as vertical momentum fluxes and the 2D contribution to update the velocity field. This step may be vertically implicit.

3D Momentum

5. The tracer fluxes (horizontal and vertical) are computed and the tracer field is advanced in time. This step may be also vertically implicit.

Tracers

(a) Schematic view of the five main components of a time step tal on

riz

m ntu

len rbu Tu

me Mo

l nta

m ntu me

o riz

r ace Tr

ace Tr

nce

ule

rb Tu

20 ms

40 ms

60 ms

80 ms

100 ms

120 ms

134 ms

(b) Timeline of a full iteration

Figure 2: (a) Schematic view of the five main components of a time step. The ordering shown corresponds to a vertically implicit step. For vertically explicit steps, the turbulence update is performed last, after the momentum and tracer updates. (b) Timeline of a full iteration of the 3D scheme. Black stripes indicate memory operations (e.g., vector initialization). The 3D momentum and tracer components are more expensive in the first step due to the vertically implicit solve. In contrast, the 2D component is more costly in the second step, as it integrates over a larger time interval (∆t vs. ∆t/2), leading to more 2D iterations. All timelines in this work were obtained on one or two A100 GPUs, using single precision, on a mesh with 210k triangles and 32 vertical layers.

18 1

27 10

20 3

29 12

23 6

32 15

25 8

34 17

21 4

30 13

Node 0

Node 1

Node 2

Node 3

Figure 3: Structure-of-Arrays (SoA) memory layout used in our implementation. Prisms within a column are ordered from top to bottom, with columns arranged sequentially. Nodal values are stored first by node order, followed by element order. For vector fields, each component field is stored contiguously rather than being interleaved. The complete hierarchy is thus field → node → column → layer.

4 8

6 10

5 9

7 11

28 32

30 34

29 33

31 35

12 16

14 18

13 17

15 19

21 {z Cell 0

}

{z Cell 1

}

Figure 4: Example of a cell layout for the same mesh as Figure 3. In this minimal example, a cell is made of two columns instead of the typical 128. This figure shows the ordering of the nodes in memory and corresponds to the memory layout of a scalar field. Vector fields are also possible with the complete memory hierarchy being cell → layer → node → field → column.

Column: 0 1 2 3 4 . . . . . .

. . . . . . 123 124 125 126 127

node 0 node 1 node 2

Layer 0

node 3 node 4

node 5 node 0 node 1 node 2

Layer 1

node 3 node 4

}

node 5

1 Figure 5: Example of a cell matrix with two layers and 128 columns for a scalar field (6 values per layer). The arrows indicate the order in which the data is stored in memory.

Column: 0 1 . . . . . . 6 7 8 . . . . . .

.. .

Layer 15 Layer 16

Second chunk of 16 layers

}

First chunk of 16 layers

Layer 0

. . . . . . 125 126 127

.. .

}

Layer 31 | {z } 8 columns handled by a given block of 128 threads

1 Figure 6: Subset of a cell processed by a block of 128 threads. This example shows a cell with 128 columns and 32 layers being processed by blocks of 128 threads. Here, the block reads from the SoA layout in chunks of 16 values and writes to the cell in chunks of 8 values. Since the cell has 32 layers, the block processes its assigned columns in two passes: first for layers 0 to 15, then for layers 16 to 31. threads per block, a typical choice is n = 8, meaning each block processes chunks of 8 columns × 16 layers at a time (Figure 6). If a column has more than 16 layers, the block iterates over the remaining layers in steps of 16, and so on. After assembling a matrix and its right-hand side (RHS) in the cell format, the system is solved while maintaining the same layout. Since all columns are independent, the solving step remains parallelizable over the individual columns. Finally, a transposition step is required to convert the results back to the SoA layout. This is performed by a kernel with the same structure as the assembly kernels. Thanks to the block-based nature of reads and writes between the cell and SoA layouts, this kernel nearly achieves peak memory bandwidth. Any imperfections in access sizes and alignment are fully masked by the cache hierarchy.

2.2

3D hydrodynamical scheme

The tracer and the hydrodynamics schemes are very similar in structure, thus we will only detail the hydrodynamics scheme here. From a high level, a step of the hydrodynamics scheme involves solving the following equation: M 1 u1 − M 0 u0 h v = F3D (u, q̄, r) + M1 (F2D /H1 ) + F3D (u, w̃, u1 ) ∆t

(11)

where the unknown is u1 , the velocity field at the end of the step, u0 is the velocity field at the start of the step, u is current the velocity field from the Runge-Kutta scheme and q̄ is the horizontal h transport used to advect u. F3D represents the discrete horizontal advection and viscosity fluxes, F2D denotes the horizontal momentum input from the external 2D mode, and the vertical fluxes are given v by F3D (u, w̃, u1 ). M0 and M1 are the mass matrices at the beginning and end of the step. The implicit v vertical fluxes have been linearized such that F3D (u, w̃, u1 ) = A(u, w̃) u1 . Using that property, equation (11) can be rewritten as h (M1 − ∆t A(u, w̃)) u1 = M0 u0 + ∆t F3D (u, q̄, r) + M1 (F2D /H1 ) (12)

where the linear system to be solved appears explicitly. We will now focus on each of these three h v , F2D and F3D (or A)), and the kernels involved in their computation. contributions (F3D 3D horizontal flux prediction

h As detailed in Section S3.2 of the supporting information, the horizontal fluxes F3D are computed twice. Once with a transport q as a prediction step, and later with a transport q̄ to advance the momentum in

h time. The prediction of F3D involves two main steps. First, the density is updated and used to compute h the hydrostatic horizontal pressure gradient r. Then, r is used to evaluate F3D (u, q, r). A detailed derivation of the discretization of these terms can be found in Sections S2.2 and S3.2 of the supporting information. The complete timeline for the assembly of the horizontal terms is illustrated in figure 7. y sit

ρ( 0

bl sem As

o e isc fre ev ix atr o SoA ute q mput m mp reco ll t lve Co So Ce P

dra m tto D o b 2 rf/ D→ Su 3

h Assemble F3D (u, q, r)

2 ms

4 ms

6 ms

8 ms

10 ms

12 ms

13.1 ms

Figure 7: Kernels used for the computation of the horizontal terms of the momentum equation. h The dominant kernels are as expected the assembly of the RHS for r and the computation of F3D . Computing the horizontal pressure gradient r requires the solution of a linear system. Fortunately, the matrix has a peculiar structure (See Section 2.3) that allows solving the system without ever assembling it, hence why the solving step is so short. Computing the horizontal transport q and the horizontal h viscosity fields are operations that could be integrated in the assembly of F3D , but it is both simpler and faster to precomputed them at all nodes beforehand. The two small kernels at the end are the surface h h and bottom drag, and vertically summing F3D to produce F3D→2D for the coupling with the external mode. The reason for both of these is explained Section S3.2 of the supporting information, but we won’t go over the details since they remain negligible in the timings.

Fluxes from the 2D external mode The external mode produces the fluxes F2D due to the gravity wave, and also accumulate the mean transport Q̄ during all the 2D iterations. The momentum change F2D is computed based on the difference of 2D transport before and after 2D iterations, as detailed in Section 1.2 of the supporting information. h F2D = (Q1 − (Q0 + ∆tF3D→2D ))/∆t

In order to get a conservative and consistent advection scheme, the vertical sum of the 3D transport q̄ must exactly match the averaged 2D transport Q̄. Therefore, at the end of each Runge-Kutta step of the 2D mode, a kernel adds the current solution to the build the time-averaged transport Q̄. This explains the 3+1 kernels for each iteration of the 2D mode illustrated in Figure 8. 2D

rs py ax

ite

rs ite py ax

ite

rs py ax

rs ite py ax

rs py ax

ite

rs py ax

ite

rs ite py ax

250 µs

500 µs

750 µs

1 ms

1.25 ms

1.5 ms

1.75 ms

2 ms

2.23 ms

Figure 8: Kernels of the 2D external mode during a full step of the scheme. In this example, the external mode performs 10 iterations of a 3-stage Runge-Kutta scheme.

3D momentum update with implicit time-stepping ee -fr w̃ atrix oA q̄ e e l S t b u m o mp sem olve ell t Co As S C

yst

s lve So

h Assemble M − ∆tA(u, w̃), F3D (u, q̄, r)

ll Ce

A So

5 ms

10 ms

15 ms

20 ms

25 ms

30 ms

35 ms 36.2 ms

Figure 9: Kernels used for the computation of the vertical terms of the momentum equation during an implicit timestep. With implicit time-stepping, the momentum update is divided into four steps. First, a corrected transport q̄ is computed from the vertically averaged 2D transport Q̄ and the velocity u. This transport 10

is then used to compute the vertical velocity w̃. Solving for w̃ requires the solution of a linear system with a structure similar to that of the horizontal pressure gradient r, which explains the short solution times observed in Figure 9. Next, a single kernel performs the assembly of both the matrix M − ∆t A(u, w̃) h and the corresponding right-hand side. The latter involves recomputing F3D (u, q̄, r) using the corrected transport q̄, and adding the contributions from M0 u0 and the 2D terms. Finally, the resulting linear system is solved. In contrast to the system for w̃, this system has no exploitable structure beyond its banded form, leading to the longer assembly and solution times shown in Figure 9. 3D momentum update with explicit time-stepping During fully explicit substeps, the temporal scheme remains almost identical. The only difference is that v F3D no longer depends on u1 , so equation (11) becomes M1 u1 − M0 u0 h v = F3D (u, r) + M1 (F2D /H1 ) + F3D (u, w̃, u) ∆t

(13)

or equivalently h v u1 = M1−1 M0 u0 + ∆t F3D (u, r) + M1 (F2D /H1 ) + F3D (u, w̃, u)

(14)

Since the mass matrix M1 is block-diagonal, solving a full linear system per column is no longer required, and the values of u1 can be computed independently for each element. Avoiding matrix assembly and u mp

q̄ te

ble sem As

w̃

e fre ixatr o SoA m t lve ell So C

v h (u, q̄, r) and solve for u1 Assemble F3D (u, w̃, u), F3D

2 ms

4 ms

6 ms

8 ms

10 ms

12 ms

14 ms

16 ms

17.3 ms

Figure 10: Kernels used for the computation of the vertical terms of the momentum equation during an explicit step. eliminating the need to write it to global memory frees up significant computational resources and memory bandwidth. This is why the explicit step is considerably faster than the implicit one.

2.3

Matrix-free solvers

The equations for the horizontal pressure gradient r and the vertical velocity w share a common structure: ∂(•) =f . ∂z In such cases, the discrete equations simplify significantly, allowing the resulting linear systems to be solved without explicitly assembling the matrix. Horizontal pressure gradient The discrete form of the horizontal pressure gradient (8) is detailed in Section S2.2 of the supporting information and essentially takes the form φr ext Jh

T̂top

+ φr int Jh

T̂bot

− ⟨r∂ẑ (φ) Jh ⟩P̂ = F .

(15)

The left-hand side of this equation couples the degrees of freedom within each vertical column, resulting in one linear system per column. Consequently, the computation of r is divided into two stages: the assembly of the right-hand side of equation (15), and solving the system. Due to the tensor product formulation of the mesh (i.e. the 3D mesh is made by extruding the 2D mesh in layers of prisms), all the terms that depend on the vertical coordinate cancel out and the system on each column can be written as Dvu r = F (16)

with Dvu a constant matrix per column that only depends on the 2D mass matrix Mh . For example, on a 3-layer column we get   −Mh −Mh  Mh −Mh      2 1 1   J 1 2Mh −Mh −Mh h  with Mh = 1 2 1 Dvu =   Mh −Mh 24 2  1 1 2  2Mh −Mh −Mh  Mh −Mh

and the pattern continues with more layers. Knowing the structure of the linear system in advance removes the need to assemble an expensive matrix. Solving the system is done via a single-pass uplooking solver that encodes the matrix structure in a recursion formula. This leads to leading to the matrix-free solver described in algorithm 1. Since the RHS of the system is stored and solved in the cell Algorithm 1 Matrix free solver for r on a single column of prisms r ← RHS M ← diag(Mh ) ∈ R6×6 s ← 0 ∈ R3×2 r ∈ R6×2 for l = 1, . . . , Nlayers do r ← LoadLayer(r, l) r ← M −1 r for i = 1, . . . , 3 do si ← si + ri + ri+3 ri ← −si + 2ri+3 ri+3 ← −si end for r ← WriteToLayer(r, r, l) end for

▷ Initialize r to solve in-place ▷ Bloc-diagonal mass matrix ▷ Initialize local accumulator ▷ Allocate the local r ▷ Loop over all layers ▷ Load layer in local memory ▷ Inverse the 2D mass matrix ▷ Loop over the 3 nodes of a face ▷ Update the accumulator ▷ Compute r on the top face ▷ Compute r on the bottom face ▷ Write the local result to global memory

format, as described in Section 2.1.1, an additional step is required after solving to convert the solution to the AoS layout. Vertical velocity Similarly, the equation for the vertical velocity (1), whose discretization is described in Section S2.1 of the supporting information, takes the form φwint n̂z Jh

T̂top

+ φwext n̂z Jh

T̂bot

− ⟨wJh ∂ẑ φ⟩P̂ = F .

As with r, the resulting linear system for each vertical column can be written as Dvd w = F . For a 3-layer column, for instance, the matrix Dvd takes the form:  Mh −Mh Mh Mh −2Mh  1 Mh −Mh Dvd =   M Mh −2Mh 2 h  Mh Mh

(17) 

      −Mh  Mh

using the same 2D mass matrix Mh as in the r system. Thanks to the predictable structure of Dvd , this system can also be solved using a matrix-free solver analogous to the one described in Algorithm 1.

2.4

Fully-assembled column solvers

Beyond the equations for r and w, where matrix-free solvers are applicable, other parts of the model require solving linear systems that couple all nodes within a column and cannot rely on a-priori know 12

matrices. This is the case in the turbulence closure model, and in the computation of vertical fluxes during the vertically implicit substeps for both tracers and hydrodynamics. The matrices involved in the momentum equation (12) and the tracer equation share a common sparsity structure. Each element has six nodes, with the top three nodes coupled to the six nodes of the element above, and the bottom three nodes coupled to the six nodes of the element below. For a 3-layer column of prisms, this results in the following sparsity pattern:               

      M − ∆tA(u, w) =       

Because memory access is expensive, only the nonzero structure of the matrix is stored during assembly. The linear systems are then solved in-place using Gaussian elimination, with one thread assigned per system. Leveraging the sparsity pattern, a local buffer of size 36 is sufficient to store the relevant portion of the matrix at any time. The right-hand side only requires 6 values (or 12 for a 2-component vector like u), all of which fit within GPU registers. This eliminates register spilling and ensures that global memory bandwidth is used efficiently for meaningful data transfers. The turbulence closure model is comparatively much simpler. It involves only a single degree of freedom per element, resulting in tridiagonal linear systems. These systems are stored using a diagonal matrix format and are solved via Gaussian elimination.

2.5

Data management

Data movement is a primary performance bottleneck on GPUs. This includes transfers within device memory and, more critically, transfers from host memory or from disk. As a regional model, SLIM relies on time-dependent external forcing (e.g., boundary currents, temperature, salinity, wind). Updating forcing fields on the CPU and transferring them to the GPU at every timestep is infeasible given the execution rate of GPU kernels, where multiple updates may be required within tens of microseconds. Consequently, forcing data must be updated directly on the device whenever possible. To this end, forcing fields are chosen to vary linearly in time between two precomputed input states, typically separated by one hour. Temporal interpolation is performed within compute kernels, thereby minimizing data transfers and avoiding additional kernel launches. When a kernel requires data outside the available temporal window, new data are loaded from disk, interpolated in space, and transferred asynchronously to the GPU. Because this loading is performed on demand, short wait times may occur; however, these delays are negligible in practice. Although not currently necessary, data prefetching may become beneficial as dataset sizes increase. Linear transport while advancing the external mode As noted in Section 2.2, tracer consistency requires the time-averaged solution of the external (2D) mode over a full step. This average is obtained by accumulating the 2D solution at each external iteration (the axpy kernels in Figure 8) and is subsequently used in the 3D tracer kernels. The same requirement applies to external forcing: the tracer must be driven by a velocity corresponding to the time average of the external transport seen by the 2D mode. The tracer requires a 3D boundary condition, whereas the external mode uses a 2D one. At the discrete level, for a full step in which the external mode performs m iterations, consistency requires m−1

1 X Qd (t0 + (i + 1/2)∆t/m) = m i=0

qd (t0 + ∆t/2) ,

Vertical DOFs

where Qd (t) and qd (t) denote the 2D and 3D external transports at time t, respectively. A consistent definition would be to set the 2D transport as the vertical sum of the 3D transport, X Qd (t) ··= qd (t), (18) Vertical DOFs

and to explicitly time-average qd so that the tracer uses m−1

qd,tracer (t0 + ∆t/2) ··=

1 X qd (t0 + (i + 1/2)∆t/m) . m i=0

However, this approach would require evaluating vertical sums of 3D data at every 2D iteration and storing the accumulated values, which is computationally expensive. Instead, we interpolate qd linearly in time over [t0 , t0 + ∆t]. Under this assumption, Qd (t0 ) and Qd (t0 + ∆t) are computed using (18), and Qd (t) is obtained by linear interpolation between these two states. Since ∆t is much smaller than the temporal resolution of the forcing data, this approximation introduces no significant loss of accuracy.

Multi-GPU

While single-GPU performance is critical, large-scale coastal simulations demand the use of many GPUs working in parallel. Multi-GPU parallelism within SLIM is handled through domain decomposition and distributed memory with MPI, with one GPU per MPI rank. The domain is split horizontally, such that each GPU owns a subset of the triangles of the 2D mesh and its associated columns. Each partition also stores a layer of ghosts triangles from the neighboring partitions to be used in the computations that require neighboring accesses. Such computations are then followed by a halo exchange, to update the ghost elements.

3.1

Structure of a compute kernel

Since GPUs significantly accelerate computations compared to CPUs, the relative cost and frequency of communications increase accordingly. Overlapping computation and communication therefore becomes essential. To achieve this, the computation within each partition is split into two parts: (i) boundary elements and elements adjacent to ghost cells, and (ii) interior elements. Boundary computations are launched first, and as soon as they complete, halo exchanges are initiated while interior computations continue. This decomposition also enables kernel specialization. Interior kernels can assume the presence of all neighbors and avoid conditional logic required at boundaries. While this distinction has limited impact on CPUs, it improves GPU performance by reducing control flow divergence and lowering resource usage (e.g., registers and shared memory). To enable overlap on the GPU, two streams are used. The compute stream handles the bulk of the computations, while the communication stream is reserved for boundary computations and halo exchange operations. The communication stream is assigned the highest priority. A full iteration proceeds as follows: 1. Boundary computations are launched on the communication stream. 2. Interior computations are launched asynchronously on the compute stream. 3. A packing kernel is launched on the communication stream to gather data for halo exchange, as the memory layout is not contiguous per element. 4. The CPU synchronizes with the communication stream to ensure packing completion, then initiates MPI communications. 5. Upon completion of MPI communications, an unpacking kernel is launched on the communication stream. GPU events are then used to ensure that the compute stream waits for communication to complete before proceeding to the next step. Since MPI predates widespread GPU usage, not all implementations are GPU-aware, i.e., capable of directly accessing device memory. In such cases, halo exchanges require additional copies between host and device memory.

3.2

MPI communications for the 3D kernel

The 3D kernels, as shown in Figures 7, 9 and 10, are quite computationally heavy, and are thus not bound by latency. As a result, the communications easily overlap with the rest of the computations and the GPU is kept busy, as seen in Figure 11. Compute kernels taking much longer than the synchronous MPI communications also means that the next kernels will already been enqueued by the time the computation is over.

3.3

Challenges with the short 2D kernels

Kernels used for the 2D mode execute in much shorter time, and thus the latency of the MPI communications starts to become a bottleneck. Figure 12 shows a typical timeline for a few iterations of the 2D mode where the white gaps between the colored blocks indicate an idling GPU. Since communications now take a significant amount of time, it is often no longer possible to launch the next kernels in advance,

ρ( 0

l mb

ity cos ree vis xf A e i t q r t u o ute mp ma o S lve ll t omp reco C So Ce P

dra m tto D o b 2 rf/ → Su 3D

h Assemble F3D (u, q, r)

Compute stream (0) .

Communications stream (1) 0s

1 ms

Pa MP Un ck I p

2 ms

3 ms

4 ms

Pa MP Un ck I p

ack

5 ms

6 ms

6.68 ms

ack

ee -fr q̄ le w̃ atrix oA S te b u m o mp ssem olve ell t Co A S C

lve So

h Assemble M − ∆tA(u, w̃), F3D (u, q̄, r)

tem SoA o ll t Ce

sys

Compute stream (0) .

Communications stream (1) 0s

2.5 ms

ute

q̄ As

5 ms

ble sem

7.5 ms

ree x-f tri SoA a m o lve ll t So Ce

10 ms

12.5 ms

15 ms

17.5 ms 18.3 ms

w̃

h v (u, w̃, u), F3D (u, q̄, r) and solve for u1 Assemble F3D

Compute stream (0) .

Communications stream (1) 0s

1 ms

2 ms

3 ms

4 ms

5 ms

6 ms

7 ms

8 ms

8.61 ms

Figure 11: Timeline of both the Compute stream and Communications stream for the two main phases of internal mode of the hydrodynamics computations. On top are the horizontal fluxes, in the middle the vertical fluxes during an implicit step, and on the bottom the vertical fluxes during an explicit step.

r ite

ite

er py D it 2

ite

e py D it 2

ite

r ite

er py D it 2

Compute stream (0) .

Communications stream (1) 0s

MP Un I pa

100 µs

200 µs

300 µs

400 µs

500 µs 518 µs

Figure 12: Timeline of both the Compute stream and Communications stream for 3 iterations of the 2D external mode. Contrary to the 3D kernels shown in Figure 11, the packing is done by the computation kernel and not separately.

and the kernel launch latency cannot be hidden. To partially mitigate this problem, the boundary computation kernel is also responsible for the packing in the MPI buffers, which saves a few microseconds of latency. This duality in MPI regimes, with heavy compute kernels for the 3D components, and much lighter 2D kernels where the communication cost is significant, will be the determinant factor in scaling the code to multiple GPUs. At the limit, when the strong scaling starts to weaken, the 2D mode takes an almost constant amount of time, due to the latency of everything, while the 3D continues to scale. This is also whey our timings are a very good fit for Amdahl’s law.

Performance and scaling of the model

On CPUs, explicit models typically exhibit predictable performance, with iteration time scaling approximately linearly with problem size. In contrast, GPU performance depends strongly on workload characteristics due to the massively parallel architecture. In particular, small problem sizes lead to underutilization of the device. This effect is further exacerbated in multi-GPU configurations with MPI, even before accounting for communication overheads. This sensitivity to workload size motivates a detailed examination of single- and multi-GPU performance.

4.1

Single-GPU performance

We first examine single-device performance across architectures, focusing on both scaling behavior and effective hardware utilization. A large GPU–CPU speedup alone does not necessarily indicate an efficient GPU implementation, as it may reflect limitations of the CPU baseline. In this work, the same codebase is used for both CPU and GPU targets, although the optimizations are primarily GPU-oriented. A comparison of single-device performance across several architectures is shown in Figure 13. The figure includes results for a 16-core AMD R9 5950X CPU (using 16 threads), consumer GPUs from NVIDIA (RTX 3060) and AMD (RX 6800XT), and data-center GPUs from NVIDIA (A100) and AMD (MI250X). The MI250X consists of two Graphics Compute Dies (GCDs), each exposed as a separate device at the programming level; results reported here correspond to a single GCD. Single precision

Time per iteration

Double precision

AMD R9 5950X CPU (18×)

AMD R9 5950X CPU (16×)

NVIDIA RTX3060 GPU (431×)

NVIDIA RTX3060 GPU (28×)

AMD 6800XT GPU (721×)

AMD 6800XT GPU (132×)

AMD MI250X GPU (510×)

AMD MI250X GPU (264×)

NVIDIA A100 GPU (1477×)

NVIDIA A100 GPU (515×)

100 ms

10 ms

1 ms 10 k

100 k

1M DG Nodes

10 M

100 M

10 k

100 k

1M DG Nodes

10 M

100 M

Figure 13: Performance of the 3D model on various hardware platforms with 32 layers and increasing horizontal resolution. The multiplier reported in the legend denotes the speedup relative to a hypothetical single-core CPU operating in double precision for a sufficiently large problem. In practice, the 16-thread execution on the AMD R9 5950X is used as the reference and assigned a multiplier of 16×. While double-precision performance is consistently lower than single precision, three distinct categories emerge. On CPUs, the performance degradation is modest, as expected for a code dominated by scalar arithmetic, for which the cost of single and double precision is similar. Data-center GPUs typically exhibit FP32/FP64 throughput ratios of 1/2 (e.g., A100) or even 1 (MI250X), resulting in double-precision performance between one-third and one-half of FP32. Memory bandwidth also plays a significant role, which explains why the MI250X is approximately twice as slow in double precision despite having identical peak throughput in FP32 and FP64. In contrast, consumer GPUs have FP32/FP64 ratios ranging from 1/16 to 1/64, resulting in a pronounced degradation in double-precision performance. Nevertheless, GPUs provide substantial overall performance gains: a laptop RTX 3060 can achieve performance comparable to approximately 430 CPU cores, an MI250X GCD exceeds 500 cores, and an A100 approaches 1500 cores. The resulting implementation makes effective use of hardware resources. Memory-bound kernels sustain up to 80% of peak bandwidth, while compute-bound kernels reach approximately 60% of peak floating-point throughput. Many kernels exhibit mixed behavior due to kernel fusion, combining memorybound and compute-intensive operations. Figure 14 shows the compute and bandwidth utilization during a typical time step. Sustained averages of approximately 30% of peak for both compute and bandwidth are notable for a low-order method,

m tal e ntu c on e z n i m r le Ho Mo rbu 3D 2D Tu 3D

r ace Tr

o riz Ho 3D 2D

l nta

nt me

ce rs bulen e c r a Tr Tu

FLOPs [% of peak]

80 60 40 20

Memory throughput [% of peak]

0 80 60 40 20 0 0s

20 ms

40 ms

60 ms

80 ms

100 ms

120 ms

134 ms

Figure 14: Global memory bandwidth and floating-point throughput as a percent of peak over a complete time-step. The sustained average is about 30% for both. Measurements performed on an A100 GPU using single precision. and reflect a high degree of efficiency despite the combination of memory-bound and compute-intensive kernels. Non-linear scaling with horizontal resolution While the code exhibits near-linear scaling at large mesh sizes, reducing the horizontal resolution does not necessarily lead to proportional performance gains. As shown in Figure 13, the iteration time on GPUs becomes nearly constant below a threshold between 105 and 106 nodes (corresponding to approximately 500 to 5000 triangles with 32 layers). A 3D step involves many 2D computations, which are significantly faster than the 3D kernels. Consequently, the 2D component is the first to underutilize the GPU as the resolution decreases. In this regime, the 2D execution time is dominated by kernel launch latency and therefore remains effectively constant. Since the 3D kernels are substantially more expensive, their latency overhead is negligible in comparison. Non-linear scaling with vertical resolution Due to the block-structured execution, the time per step does not scale linearly with the number of layers, as shown in Figure 15. As described in Section 2.1.2, each thread block processes all layers of its assigned columns. This can lead to partial thread utilization. For instance, with 20 layers, an efficient configuration assigns 6 columns per block, yielding 120 active threads out of 128. Moreover, memory writes occur in chunks equal to the number of columns, here 6, which is suboptimal for memory coalescence. For larger layer counts, columns are split vertically and processed sequentially by threads within the block from top to bottom. As one can expect, a number of layers that perfectly divides the block size will often result in better performance. Because each block processes all layers of its assigned columns, sharing 2D data among threads within a column is advantageous. This is implemented using shared memory. To reduce complexity and maximize performance, only statically allocated shared memory is used, with its size fixed at compile time. Consequently, an upper bound on the number of columns per block must be specified in advance, and shared memory for 2D data is allocated accordingly. In SLIM, this maximum is set to 32 columns 19

1×A100

1×6800XT

Normalized time per iteration

1.0

0.8

0.6

0.4

0.2

0.0 1 4

21 25 32 36 Number of layers

1 4

21 25 32 36 Number of layers

Figure 15: Normalized time per iteration of the 3D scheme as a function of the number of layers for an NVIDIA A100 GPU (left) and AMD 6800XT GPU (right) using single precision. The visible dips at values like 16, 32 and 64 are numbers that allow all threads from a block to be used, and also have a good coalescence in the memory accesses. per block, with a block size of 128 threads. As a result, when fewer than 4 layers are present, the total number of elements per block is below 128, leaving some threads idle. This accounts for the nearly constant runtime observed for 1 to 4 layers in Figure 15.

4.2

Multi-GPU scaling

SLIM is designed for distributed-memory execution using MPI, with one GPU per rank. The scaling behavior closely mirrors that observed on a single GPU, with latency-dominated components limiting performance at small workloads. On a single GPU, this latency is primarily due to kernel launches; in the multi-GPU setting, it additionally includes MPI communication overheads and the cost of packing and unpacking data. This overhead is partially mitigated through the overlap strategy described in Section 3. However, due to the relatively low computational cost of the 2D component, its associated latency cannot be fully hidden and remains the dominant limitation for strong scaling. As a result, the overall scaling of the 3D model is well described by Amdahl’s law. The latency associated with the 2D component effectively acts as the sequential fraction of the computation, while the 3D component exhibits near-ideal strong scaling. Benchmarks were conducted on two EuroHPC systems. The first is MeluXina, equipped with 4 NVIDIA A100 (40 GB) GPUs per node, using a non GPU-aware OpenMPI implementation. Within each node, GPUs are interconnected via NVLink, while communication with the host and network relies on PCIe. The second system is LUMI (GPU partition LUMI-G), equipped with 4 AMD MI250X GPUs per node, each providing 128 GB of memory. Each MI250X consists of two Graphics Compute Dies (GCDs), exposed as separate devices at the programming level; results are therefore reported in terms of the number of GCDs. Within a node, MI250X modules are connected through high-speed GPU–GPU links, and each GPU is directly connected to the Slingshot-11 interconnect, providing 200 Gb/s per network endpoint. GPU-aware Cray MPICH is used on this system. Figure 16 illustrates the scaling on the MeluXina cluster with A100 GPUs. Two distinct latency regimes can be observed: on a single GPU, the minimum time per iteration is approximately 3 ms, while it increases to about 6 ms as soon as multiple GPUs are used. Importantly, this latency remains essentially constant as the number of GPUs increases, indicating that network contention is not a limiting factor.

Figure 16: Scaling of the 3D model with 32 layers on the MeluXina cluster with A100 GPUs. This configuration uses a time step ratio of 20 between the internal and external modes, i.e., one internal iteration corresponds to 20 external iterations, which is representative of typical applications. Each internal iteration involves approximately 100 halo exchanges, about 90% of which originate from the 2D mode. This corresponds to roughly 100 stream synchronizations, 100 MPI send/receive operations, and 200 additional kernel launches. Assuming comparable costs for these operations, the additional 3 ms

per iteration corresponds to an average overhead of approximately 7.5 µs per synchronization, communication, or kernel launch, consistent with the expected order of magnitude. Figure 17 shows the efficiency as a function of the number of elements per GPU. For two or more GPUs, the efficiency depends primarily on this quantity, indicating near-ideal weak scaling. The observed behavior is also well described by Amdahl’s law, as expected. A100-MLUX w/32 layers : Efficiency 1.0

Efficiency

0.8

0.6

A100-MLUX 2×A100-MLUX 4×A100-MLUX 8×A100-MLUX 16×A100-MLUX 32×A100-MLUX 64×A100-MLUX 128×A100-MLUX Ahmdal’s law

0.4

0.2

0.0

102

103

104 Number of triangles per GPU

105

106

Figure 17: Efficiency of the 3D model with 32 layers on the MeluXina cluster with A100 GPUs. Figure 18 illustrates the scaling on the LUMI cluster with MI250X GPUs. Owing to the larger scale of the system, experiments were conducted with significantly larger configurations, reaching up to 1024 Graphics Compute Dies (GCDs), corresponding to 512 GPUs across 128 nodes. Although efficiency decreases at large scale, the model continues to scale provided the problem size per GPU remains sufficiently large. As in many parallel applications, SLIM involves a trade-off between time to solution and resource efficiency. An efficiency of approximately 80% is often considered a practical target. For SLIM, this corresponds to about 4 × 104 triangles per GPU, or approximately 5 × 106 nodes globally. Given that the intended use case involves simulations with at most a few million triangles distributed over 32 to 64 GPUs, the current implementation achieves satisfactory performance.

Figure 18: Scaling of the 3D model with 32 layers on the LUMI cluster equipped with MI250X GPUs.

Application to the Great Barrier Reef

To demonstrate that the model is applicable beyond idealized benchmarks and can handle realistic coastal domains, we consider a high-resolution simulation of the Great Barrier Reef (GBR). The circulation in the GBR is driven by a combination of large-scale currents, winds, tides, and complex bathymetry. In particular, the South Equatorial Current (SEC) bifurcates near 14–18◦ S into the northward Coral Sea Coastal Current and the southward East Australian Current (Lambrechts et al., 2008). These currents, modulated by local forcing, control reef flushing, larval connectivity, and the transport of nutrients and pollutants. In addition, small-scale features such as tidal jets and eddies (from 100 m to a few kilometers) play a key role in reef-scale processes. Modeling these interacting processes across a wide range of spatial scales remains challenging. Uniformresolution models in the GBR, typically ranging from 1.5 km to 4 km, are unable to resolve these fine-scale dynamics. The present approach addresses this limitation using a highly non-uniform mesh, refined near reefs and progressively coarsened offshore, enabling the representation of both basin-scale circulation and reef-scale processes. Figure 19 shows the computational domain and horizontal mesh. The mesh consists of 3.3 million triangles, with resolution ranging from 200 m near reefs to 10 km offshore. The 2D mesh is extruded vertically into a prismatic grid with 10 layers in shallow regions and up to 29 layers in deeper areas, resulting in approximately 34 million elements. The resulting resolution enables the model to capture a wide range of flow features, as illustrated below. Figure 20 shows the surface vertical vorticity at increasing levels of zoom, highlighting flow structures down to scales of a few hundred meters. Mesh resolution has a strong impact on coral connectivity, as demonstrated by (Saint-Amand et al., 2023). Figure 21 illustrates the fine-scale variability in the surface temperature field, capturing features associated with complex reef topography. In particular, resolving tidal jets and their influence on local flow patterns requires sufficiently fine spatial resolution. The simulation was performed on 32 AMD MI250X GPUs (64 GCDs) in double precision, achieving a throughput of approximately 100 simulated days per day of wall-clock time. This demonstrates that high-resolution coastal simulations of this scale are computationally feasible on modern GPU systems. While the present configuration is not yet calibrated for scientific analysis, it provides a representative test of model performance and scalability in a realistic setting. The model setup relies on multiple data sources. Bathymetry is derived from 30 m datasets for the Great Barrier Reef and Torres Strait (Beaman, 2020a, 2023), complemented by 100 m datasets for the Coral Sea (Beaman, 2020b) and the Gulf of Papua (Daniell, 2020). Coastlines are obtained from OpenStreetMap (OpenStreetMap contributors, 2015). Open boundary conditions (temperature, salinity, and currents) are derived from the BRAN2023 reanalysis dataset (Chamberlain & Commonwealth Scientific and Industrial Research Organisation, 2024). Tidal forcing is provided by TPXO10v2 (Egbert & Erofeeva, 2002), while atmospheric forcing (wind, precipitation, and heat fluxes) is obtained from the BARRA2-C2 regional atmospheric reanalysis dataset (Bureau of Meteorology, 2023), described by (Su et al., 2022, 2024). Reef extent is based on (“CoralMapping/AllenCoralAtlas”, 2020; UNEP-WCMC et al., 2021). The density ρ is computed from the equation of state ρ(S, T, p) following (Jackett et al., 2006).

Figure 19: Computational mesh for the Great Barrier Reef configuration. The horizontal resolution varies from 200 m in reef regions to 10 km in the open ocean, enabling the resolution of fine-scale coastal dynamics while maintaining tractable computational cost offshore. Insets highlight the progressive refinement near complex topographical structures.

Figure 20: Surface vertical vorticity at increasing levels of zoom in the Great Barrier Reef domain. The solution reveals fine-scale flow structures down to scales of a few hundred meters, illustrating the ability of the model and mesh to capture small-scale dynamics induced by the topography.

Figure 21: Modelled sea surface temperature on October 31, 2024 at 10:00, after two months of simulation. The solution resolves fine-scale spatial variability associated with reef topography, including signatures of tidal jets and localized mixing processes.

Conclusion

We demonstrated that unstructured-mesh ocean models based on the Discontinuous Galerkin finite element method can effectively leverage modern GPU architectures. The DG-FE formulation is inherently well suited to GPU execution, enabling hardware utilization levels that can rival structured-grid models. The resulting multi-GPU implementation scales efficiently from consumer laptops to large HPC clusters with hundreds of devices. This performance gain makes it feasible to perform three-dimensional coastal simulations at resolutions and scales that were previously difficult to attain. The DG-FE formulation offers properties that align naturally with GPU architectures: high data locality, element-wise independence in explicit computations, and a high arithmetic intensity relative to other unstructured-mesh approaches. These properties allowed us to sustain up to 80% of peak memory bandwidth for memory-bound kernels and approximately 60% of peak floating-point throughput for compute-bound kernels, with a sustained average of about 30% of peak compute and memory throughput over a complete time step. These figures are notable for a low-order method on an unstructured mesh, and are comparable to utilization rates reported for structured-grid GPU models. For example, Veros achieves strong performance gains through JAX just-in-time compilation (Häfner et al., 2021) but on a regular grid, while Oceananigans.jl reaches high efficiency through kernel-fusion on a structured finite-difference grid (Silvestri et al., 2025). Our results demonstrate that comparable efficiency can be achieved on unstructured meshes when the DG formulation is exploited through appropriate data layouts (structureof-arrays with Hilbert curve reordering), matrix-free solvers for vertically structured operators, and a dedicated cell layout for column-wise implicit solves. Our implementation achieved efficient scaling from a single consumer GPU to 1024 HPC-grade devices across two state-of-the-art EuroHPC systems (MeluXina and LUMI). A single NVIDIA A100 delivers performance equivalent to approximately 1500 CPU cores, and replacing a 128-core CPU node with a 4×A100 GPU node yields a speedup of approximately 50. These gains are preserved as the number of devices increases: weak-scaling efficiency remains high, and strong scaling follows Amdahl’s law closely, with the latency-dominated 2D external mode acting as the effectively sequential fraction. This performance is obtained from a single codebase targeting CPU, CUDA, and HIP backends through a lightweight abstraction layer, ensuring portability across architectures with minimal overhead. This scaling behavior reflects the dual nature of the mode-splitting approach common to most ocean models (Kärnä et al., 2018; Madec et al., 2022; Shchepetkin & McWilliams, 2005). The 3D baroclinic component, which dominates the computational cost, exhibits near-ideal scaling owing to its high arithmetic intensity and effective computation-communication overlap. In contrast, the 2D barotropic mode, while computationally cheap, involves many short kernels and frequent halo exchanges, making it sensitive to MPI latency and kernel launch overhead. This dichotomy has also been observed in GPU implementations of structured-grid models (Xu et al., 2014). Similar challenges are reported with the barotropic mode in LICOM3-HIP (Wei et al., 2024). The overlap strategy employed here, where boundary elements are processed first to initiate asynchronous halo exchanges while interior computations proceed, partially mitigates this issue but cannot fully hide the latency at extreme strong-scaling limits. GPU acceleration makes high-resolution, three-dimensional coastal simulations over large domains computationally feasible. The Great Barrier Reef application, with mesh resolution down to 200 m near the reefs and 34 million prisms, ran on 32 MI250X GPUs (64 GCDs) while maintaining a physical-tonumerical time ratio of approximately 100. This represents the first application of a 3D model to the entire GBR with a sub-reef scale resolution. Such a model can resolve fine-scale features such as tidal jets and reef-scale eddies that are critical for larval connectivity, nutrient transport, and sediment dynamics. This capability addresses a longstanding bottleneck in coastal ocean modeling. While unstructuredmesh models have long been recognized for their geometrical flexibility, their higher computational cost per degree of freedom compared to structured-grid models has historically limited their practical application (Danilov et al., 2017). The GPU acceleration presented here effectively removes this constraint for regional and coastal domains. The present implementation brings the cost of high-resolution unstructured simulations in line with what is routinely achieved on structured grids at coarser resolution, potentially opening the door to operational applications in environmental management. Some limitations should nonetheless be acknowledged. The performance scaling with mesh resolution is not strictly linear. GPU utilization saturates below approximately 106 DG nodes, corresponding to roughly 5000 triangles with 32 layers, below which the device is underutilized and iteration time becomes dominated by kernel launch latency. Similarly, performance exhibits non-monotonic behavior with the number of vertical layers due to thread-block occupancy constraints: layer counts that do not evenly divide the block size result in idle threads and suboptimal memory coalescence. While these effects are

well understood and primarily affect small simulations, they imply that the model is most efficient for medium-to-large problems with layer counts aligned to powers of two. Furthermore, the Great Barrier Reef simulation, while demonstrating the model’s scalability and resolution capabilities, has not been calibrated or validated against observational data. The results should therefore be interpreted as a computational feasibility demonstration rather than a scientifically validated study of GBR dynamics. More broadly, this work shows that the DG-FE approach, long considered too expensive for routine use in ocean modeling, can become computationally competitive with structured-grid methods when GPU architectures are effectively exploited. As exascale computing systems become more widely available and GPU architectures continue to evolve, the inherent advantages of unstructured meshes, geometrical flexibility, local refinement, and multi-scale resolution, may increasingly be leveraged without the traditional computational overhead. In this context, GPU-accelerated DG models emerge as a promising tool for the next generation of coastal and regional ocean simulations, where high resolution, complex topography, and multi-physics coupling are simultaneously required.

Open Research The SLIM model source code used in this study is openly developed at https://git.immc.ucl.ac. be/slim/slim4 and is distributed under the GNU General Public License v3.0 or later (GPLv3+). Documentation, installation instructions, and examples are available at https://slim.git-page.immc. ucl.ac.be/slim4/v-0.9/index.html, while an overview of the model and publications using SLIM are provided at https://www.slim-ocean.be/. A permanent, citable archive of the version of SLIM used for this manuscript, together with the scripts and input files needed to reproduce the numerical experiments and figures, will be deposited in a public repository such as Zenodo and assigned a DOI before publication. The synthetic benchmark configurations used to evaluate the multi-GPU implementation are available in the SLIM benchmark example at https://slim.git-page.immc.ucl.ac.be/slim4/v-0.9/ examples/Benchmark.html. These benchmarks are designed to assess computational performance and do not rely on external geophysical data. The Great Barrier Reef (GBR) model configuration used for the realistic performance test is available in the SLIM example repository at https://slim.git-page.immc.ucl.ac.be/slim4/v-0.9/examples/ Great_Barrier_Reef_testcase.html. This configuration is provided to reproduce the performanceoriented test case reported in this manuscript. It was not calibrated or validated for scientific analysis of GBR hydrodynamics, and should therefore not be interpreted as a research-grade GBR circulation simulation. The public GBR example uses bathymetry, coastline, atmospheric forcing, ocean boundary conditions, and coral reef distribution data obtained from public data providers. The bathymetric data are from the Geoscience Australia Great Barrier Reef and Coral Sea depth products (Beaman, 2020a, 2020b, 2023) and the Gulf of Papua bathymetry data set (Daniell, 2020). Coral reef distribution data are from the Allen Coral Atlas (“CoralMapping/AllenCoralAtlas”, 2020) and the global coral reef distribution data set compiled by UNEP-WCMC, WorldFish Centre, WRI, and TNC (UNEP-WCMC et al., 2021). Coastline data are from OpenStreetMap (OpenStreetMap contributors, 2015). Atmospheric forcing is from BARRA2, the Australian Regional Atmospheric Reanalysis (Bureau of Meteorology, 2023), and ocean boundary and initial conditions are from BRAN2023, the Bluelink Ocean Reanalysis (Chamberlain & Commonwealth Scientific and Industrial Research Organisation, 2024). Access conditions and licenses are those specified by the respective data providers. The original TPXO tidal forcing data used in the internal GBR performance tests cannot be redistributed, because access to TPXO products must be requested directly from the data provider. The archived and public GBR example therefore uses simplified tidal forcing with realistic magnitudes. This substitution does not affect the performance conclusions of the manuscript, because the GBR case is used only to benchmark the implementation under realistic mesh, bathymetry, and forcing complexity, not to draw scientific conclusions about tides or circulation in the Great Barrier Reef.

Acknowledgments This work was supported by the Special Research Fund (FSR) - UCLouvain. Computational resources have been provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI) funded by the Fond de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under convention 2.5020.11 and by the Walloon Region. The present research also benefited from computational resources made available on Lucia, the Tier1 supercomputer of the Walloon Region, infrastructure funded by the Walloon Region under the grant agreement n°1910247. The simulations on multiple NVIDIA GPUS were performed on the Luxembourg national supercomputer MeluXina. The authors gratefully acknowledge the LuxProvide teams for their expert support. The authors acknowledge LUMI-BE for awarding this project access to the LUMI supercomputer, owned by the EuroHPC Joint Undertaking, hosted by CSC (Finland) and the LUMI consortium through a LUMI-BE Regular Access call. LUMI-BE is joint effort from BELSPO (federal), SPW Économie, Emploi, Recherche (Wallonia), Department of Economy, Science & Innovation (Flanders) and Innoviris (Brussels). The authors declare there are no conflicts of interest for this manuscript.

References Allen coral atlas [dataset]. (2020). https://doi.org/10.5281/zenodo.3833242 Beaman, R. J. (2020a). Ausbathytopo (great barrier reef) 30m 2017 – a regional-scale depth model (20170025c) [dataset]. https://doi.org/10.4225/25/5a207b36022d2 Beaman, R. J. (2020b). High resolution depth model for the great barrier reef and coral sea 100 m [dataset]. https://doi.org/10.26186/5E2F8BB629D07 Beaman, R. J. (2023). Torres strait bathymetry 30m 2020 – a high-resolution depth model (20200021c) [dataset]. https://doi.org/10.26186/144348 Bleck, R. (2002). An oceanic general circulation model framed in hybrid isopycnic-cartesian coordinates. Ocean Modelling, 4 (1), 55–88. https://doi.org/10.1016/S1463-5003(01)00012-9 Bureau of Meteorology. (2023). Bureau of meteorology atmospheric high-resolution regional reanalysis for australia – version 2 (barra2) [dataset]. https://doi.org/10.25914/1X6G-2V48 Chamberlain, M., & Commonwealth Scientific and Industrial Research Organisation. (2024). Bluelink Ocean Reanalysis - BRAN2023 [dataset]. https://doi.org/10.25914/2WXJ-VT48 Chang, C., Deringer, V. L., Katti, K. S., Van Speybroeck, V., & Wolverton, C. M. (2023). Simulations in the era of exascale computing. Nature Reviews Materials, 8 (5), 309–313. https://doi.org/ 10.1038/s41578-023-00540-6 Chen, C., Liu, H., & Beardsley, R. C. (2003). An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries. Journal of Atmospheric and Oceanic Technology, 20 (1), 159–186. https : / / doi . org / 10 . 1175 / 1520 0426(2003)020<0159:AUGFVT>2.0.CO;2 Dally, W. J., Keckler, S. W., & Kirk, D. B. (2021). Evolution of the graphics processing unit (gpu). IEEE Micro, 41 (6), 42–51. https://doi.org/10.1109/MM.2021.3113475 Daniell, J. (2020, March). Gulf of papua bathymetry raster dataset [dataset]. https://doi.org/10. 6084/m9.figshare.11986797.v1 Danilov, S., & Wang, Q. (2015). Resolving eddies by local mesh refinement. Ocean Modelling, 93, 75–83. https://doi.org/10.1016/j.ocemod.2015.07.006 Danilov, S., Sidorenko, D., Wang, Q., & Jung, T. (2017). The finite-volume sea ice-ocean model (fesom2). Geoscientific Model Development, 10 (2), 765–789. https://doi.org/10.5194/gmd-10-7652017 Dawson, C., Trahan, C. J., Kubatko, E. J., & Westerink, J. J. (2013). A parallel local timestepping runge– kutta discontinuous galerkin method with applications to coastal ocean modeling. Computer Methods in Applied Mechanics and Engineering, 259, 154–165. https://doi.org/10.1016/j. cma.2013.03.015 Deleersnijder, E., Legat, V., & Lermusiaux, P. F. J. (2010). Multi-scale modelling of coastal, shelf, and global ocean dynamics. Ocean Dynamics, 274. https://doi.org/10.1007/s10236-010-0363-6 Egbert, G. D., & Erofeeva, S. Y. (2002). Efficient inverse modeling of barotropic ocean tides. Journal of Atmospheric and Oceanic Technology, 19 (2), 183–204. https://doi.org/10.1175/15200426(2002)019<0183:EIMOBO>2.0.CO;2 Häfner, D., Nuterman, R., & Jochum, M. (2021). Fast, cheap, and turbulent—global ocean modeling with gpu acceleration in python. Journal of Advances in Modeling Earth Systems, 13 (12), e2021MS002717. https://doi.org/10.1029/2021MS002717 Ishimwe, A. P., Deleersnijder, E., Legat, V., & Lambrechts, J. (2023). A split-explicit second-order runge–kutta method for solving 3d hydrodynamic equations. Ocean Modelling, 186, 102273. https://doi.org/10.1016/j.ocemod.2023.102273 Ishimwe, A. P., Deleersnijder, E., Legat, V., & Lambrechts, J. (2025). A multi-scale imex second-order runge–kutta method for 3d hydrodynamic ocean models. Journal of Computational Physics, 520, 113482. https://doi.org/10.1016/j.jcp.2024.113482 Jackett, D. R., McDougall, T. J., Feistel, R., Wright, D. G., & Griffies, S. M. (2006). Algorithms for density, potential temperature, conservative temperature, and the freezing temperature of seawater. Journal of Atmospheric and Oceanic Technology, 23 (12), 1709–1728. https://doi.org/ 10.1175/JTECH1946.1 Kärnä, T. (2020). Discontinuous galerkin discretization for two-equation turbulence closure models. Ocean Modelling, 150, 101619. https://doi.org/10.1016/j.ocemod.2020.101619 Kärnä, T., Kramer, S., Mitchell, L., Ham, D., Piggott, M., & Baptista, A. (2018). Thetis coastal ocean model: Discontinuous galerkin discretization for the three-dimensional hydrostatic equations.

Geoscientific Model Development, 11 (11), 4359–4382. https://doi.org/10.5194/gmd- 114359-2018 Kärnä, T., Legat, V., & Deleersnijder, E. (2013). A baroclinic discontinuous galerkin finite element model for coastal flows. Ocean Modelling, 61, 1–20. https://doi.org/10.1016/j.ocemod.2012.09. 009 Korn, P., Brüggemann, N., Jungclaus, J. H., Lorenz, S. J., Gutjahr, O., Haak, H., Linardakis, L., Mehlmann, C., Mikolajewicz, U., Notz, D., Putrasahan, D. A., Singh, V., von Storch, J.-S., Zhu, X., & Marotzke, J. (2022). Icon-o: The ocean component of the icon earth system model—global simulation characteristics and local telescoping capability. Journal of Advances in Modeling Earth Systems, 14 (10), e2021MS002952. https://doi.org/10.1029/2021MS002952 Lambrechts, J., Hanert, E., Deleersnijder, E., Bernard, P.-E., Legat, V., Remacle, J.-F., & Wolanski, E. (2008). A multi-scale model of the hydrodynamics of the whole great barrier reef. Estuarine, Coastal and Shelf Science, 79 (1), 143–151. https://doi.org/10.1016/j.ecss.2008.03.016 Lermusiaux, P. F. J., Schröter, J., Danilov, S., Iskandarani, M., Pinardi, N., & Westerink, J. J. (2013). Multiscale modeling of coastal, shelf, and global ocean dynamics. Ocean Dynamics, 63 (11), 1341–1344. https://doi.org/10.1007/s10236-013-0655-8 Madec, G., Bourdallé-Badie, R., Chanut, J., Clementi, E., Coward, A., Ethé, C., Iovino, D., Lea, D., Lévy, C., Lovato, T., Martin, N., Masson, S., Mocavero, S., Rousset, C., Storkey, D., Müeller, S., Nurser, G., Bell, M., Samson, G., . . . Moulin, A. (2022, March). Nemo ocean engine [software]. https://doi.org/10.5281/zenodo.6334656 Marshall, J., Adcroft, A., Hill, C., Perelman, L., & Heisey, C. (1997). A finite-volume, incompressible navier-stokes model for studies of the ocean on parallel computers. Journal of Geophysical Research: Oceans, 102 (C3), 5753–5766. https://doi.org/10.1029/96JC02775 Okubo, A. (1971). Oceanic diffusion diagrams. Deep Sea Research and Oceanographic Abstracts, 18 (8), 789–802. https://doi.org/10.1016/0011-7471(71)90046-5 OpenStreetMap contributors. (2015). Planet dump retrieved from https://planet.osm.org [dataset]. Panzer, I., Lines, S., Mak, J., Choboter, P., & Lupo, C. (2013). High performance regional ocean modeling with gpu acceleration. OCEANS 2013 - San Diego, 1–4. https://doi.org/10.23919/OCEANS. 2013.6741366 Saint-Amand, A., Lambrechts, J., & Hanert, E. (2023). Biophysical model resolution affects coral connectivity estimates. Scientific Reports, 13 (1), 9414. https://doi.org/10.1038/s41598-02336158-5 Seny, B., Lambrechts, J., Toulorge, T., Legat, V., & Remacle, J.-F. (2014). An efficient parallel implementation of explicit multirate runge–kutta schemes for discontinuous galerkin computations. Journal of Computational Physics, 256, 135–160. https://doi.org/10.1016/j.jcp.2013.07.041 Shchepetkin, A. F., & McWilliams, J. C. (2005). The regional oceanic modeling system (roms): A splitexplicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling, 9 (4), 347–404. https://doi.org/10.1016/j.ocemod.2004.08.002 Silvestri, S., Wagner, G. L., Constantinou, N. C., Hill, C. N., Campin, J.-M., Souza, A. N., Bishnu, S., Churavy, V., Marshall, J., & Ferrari, R. (2025). A gpu-based ocean dynamical core for routine mesoscale-resolving climate simulations. Journal of Advances in Modeling Earth Systems, 17 (4), e2024MS004465. https://doi.org/10.1029/2024MS004465 Silvestri, S., Wagner, G. L., Hill, C., Ardakani, M. R., Blaschke, J., Campin, J.-M., Churavy, V., Constantinou, N. C., Edelman, A., Marshall, J., Ramadhan, A., Souza, A., & Ferrari, R. (2024, October). Oceananigans.jl: A julia library that achieves breakthrough resolution, memory, and energy efficiency in global ocean simulations [software]. https://doi.org/10.48550/arXiv. 2309.06662 Smagorinsky, J. (1963). General circulation experiments with the primitive equations: I. the basic experiment. Monthly Weather Review, 91 (3), 99–164. https://doi.org/10.1175/1520-0493(1963) 091<0099:GCEWTP>2.3.CO;2 Strohmaier, E., Meuer, H. W., Dongarra, J., & Simon, H. D. (2025). The top500 list and progress in high-performance computing. Proceedings of the ACM/IEEE Supercomputing Conference (SC). Su, C.-H., Rennie, S., Dharssi, I., Torrance, J., Smith, A., Le, T., Steinle, P., Stassen, C., Warren, R. A., Wang, C., & Le Marshall, J. (2022). Barra2: Development of the next-generation australian regional atmospheric reanalysis (Bureau Research Report No. 067). Bureau of Meteorology. http://www.bom.gov.au/research/publications/researchreports/BRR-067.pdf Su, C.-H., Rennie, S., Torrance, J., Howard, E., Stassen, C., Lipson, M., Warren, R., Pepler, A., Dharssi, I., & Franklin, C. (2024). Barra-c2: Development of the kilometre-scale downscaled atmospheric

reanalysis over australia (Bureau Research Report No. 097). Bureau of Meteorology. http : //www.bom.gov.au/research/publications/researchreports/BRR-097.pdf Umgiesser, G., Canu, D. M., Cucco, A., & Solidoro, C. (2004). A finite element model for the venice lagoon: Development, setup, calibration, and validation. Journal of Marine Systems, 51 (1), 123– 145. https://doi.org/10.1016/j.jmarsys.2004.05.009 Umlauf, L., & Burchard, H. (2003). A generic length-scale equation for geophysical turbulence models. Journal of Marine Research, 61 (2). https://doi.org/10.1357/002224003322005087 UNEP-WCMC, WorldFish Centre, WRI, & TNC. (2021). Global distribution of coral reefs, compiled from multiple sources including the millennium coral reef mapping project [dataset]. https : //doi.org/10.34892/t2wk-5t34 Includes contributions from IMaRS-USF and IRD (2005); IMaRS-USF (2005); Spalding et al. (2001). Wagner, G. L., Silvestri, S., Constantinou, N. C., Ramadhan, A., Campin, J.-M., Hill, C., Chor, T., Strong-Wright, J., Lee, X. K., Poulin, F., Souza, A., Burns, K. J., Marshall, J., & Ferrari, R. (2025, February). High-level, high-resolution ocean modeling at all scales with oceananigans. https://doi.org/10.48550/arXiv.2502.14148 Wang, Q., Danilov, S., Sidorenko, D., Timmermann, R., Wekerle, C., Wang, X., Jung, T., & Schröter, J. (2014). The finite element sea ice-ocean model (fesom) v1.4: Formulation of an ocean general circulation model. Geoscientific Model Development, 7 (2), 663–693. https : / / doi . org / 10 . 5194/gmd-7-663-2014 Wei, J., Lin, P., Jiang, J., Liu, H., Zhao, L., Zhang, Y., Han, X., Zhang, F., Huang, J., Wang, Y., Li, Y., Yu, Y., & Chi, X. (2024). Accelerating lasg/iap climate system ocean model version 3 for performance portability using kokkos. Future Generation Computer Systems, 160, 901–917. https://doi.org/10.1016/j.future.2024.06.029 Westerink, J. J., Luettich, R. A., Feyen, J. C., Atkinson, J. H., Dawson, C., Roberts, H. J., Powell, M. D., Dunion, J. P., Kubatko, E. J., & Pourtaheri, H. (2008). A basin- to channel-scale unstructured grid hurricane storm surge model applied to southern louisiana. Monthly Weather Review, 136 (3), 833–864. https://doi.org/10.1175/2007MWR1946.1 Xu, S., Huang, X., Zhang, Y., Fu, H., Oey, L.-Y., Xu, F., & Yang, G. (2014). Gpupom: A gpu-based princeton ocean model. Geoscientific Model Development Discussions, 7, 7651–7691. https : //doi.org/10.5194/gmdd-7-7651-2014 Zhang, Y. J., Ye, F., Stanev, E. V., & Grashorn, S. (2016). Seamless cross-scale modeling with schism. Ocean Modelling, 102, 64–81. https://doi.org/10.1016/j.ocemod.2016.05.002

Supporting Information for “An efficient multi-GPU implementation for the Discontinuous Galerkin ocean model SLIM” M. De Le Court1 , V. Legat1 , A. P. Ishimwe3 , C. Scherpereel2 , E. Hanert1,2 , and J. Lambrechts1 1 Institute of Mechanics, Materials and Civil Engineering, UCLouvain, Louvain-la-Neuve, Belgium 2 Earth and Life Institute, UCLouvain, Louvain-la-Neuve, Belgium 3 Department of Ecoscience, Aarhus University, Roskilde, Denmark

Spatial Discretization

is slightly non-conformal along lateral faces. To avoid ambiguity in the integration domain, all integrals are expressed in the parent element. This change of variables introduces a jacobian, denoted as J when integrating inside an element, or J∂ when integrating on its boundary. Depending on the context, J and J∂ are sometimes expressed in a different form. Due to the tensor product formulation, the volume Jacobian J is the product of the 2D Jacobian (Jh ) and the vertical 1D Jacobian (Jz ). The same is true for the lateral jacobian where JQ = Jz Jl , with Jl the jacobian of the edge of the triangle. On the top and bottom faces of the prisms, the jacobian can also be split as JT = Jh /|nz |. Finally, to simplify notation, we denote volume integrals over a domain Ω with single angled brackets and boundary integrals over ∂Ω with double brackets: Z Z f dΩ = ⟨f ⟩ , f dσ(Ω) = ⟨⟨f ⟩⟩ .

The spatial discretization used in this article is based on the Discontinuous Galerkin (DG) method with linear nodal basis functions. Since the 3D mesh is constructed by vertically extruding the 2D mesh, the corresponding basis functions φ are defined as the tensor product of the 2D horizontal basis functions φh and a vertical 1D basis function φz : φ(ξ, η, ζ) = φh (ξ, η) · φz (ζ) . When using nodal basis functions, we sometimes refer to a single 3D index i, or to a pair of indices (ih , iz ) corresponding to the horizontal and vertical indices, respectively: φi ··= φih ,iz (ξ, η, ζ) ··= φihh (ξ, η) φizz (ζ) . Since the top and bottom boundaries of the 3D elements may not be horizontal, the vertical coordinate ζ is a function of the physical coordinates (x, y, z), while ξ and η do not depend on z. Using the DG method introduces discontinuities (jumps) in discretized fields at element interfaces. We therefore define the following operators for the mean value, mean difference, maximum, and upwind value on an interface: {α} =

αint + αext , 2

JαK =

⌈ α⌉⌉ = max(αint , αext ),

αup = int

2D equations

All discrete equations are derived from the weak form of their corresponding continuous equations, using the same set of basis functions φ for both the test and trial spaces. During the derivation, interface terms that have not yet been explicitly defined are marked with an asterisk (∗ ).

αint − αext , 2 ( αint , if n · u∗ ≥ 0, ext

∂Ω

Ω

S1.1

∗

if n · u < 0,

Free Surface

We start with the free-surface equation (1), which forms the first component of the external mode:

ext

for any field α, where α and α denote the interior and exterior values, respectively. The vector u∗ is the discrete advection velocity at the interface, and n is the unit outward normal vector. When a quantity is taken in its upwind form, it is always multiplied by the velocity field from which it originates. Since SLIM uses triangular meshes in 2D, the 3D elements are prisms, denoted by P, with the corresponding triangle denoted T2D . The top and bottom faces of a prism are referred to as Ttop and Tbot , collectively called the horizontal boundaries T = Ttop ∪ Tbot . The lateral boundaries are denoted Q, and the full prism boundary is ∂P. In 2D, the edges of the triangle T2D are denoted L. When integrating over the parent element, a “hat” is used to distinguish it: P̂ denotes the parent prism, while T̂ and Q̂ refer to its horizontal and lateral faces, respectively. Because the free-surface elevation is discretized discontinuously and defines the shape of the 3D mesh, the mesh

∂t η = −∇h · Q + s

(1)

Its weak form, after integrating by parts, is ⟨φ∂t η⟩ = ⟨∇h φ · Q⟩ − ⟨⟨φn · Q∗ ⟩⟩ + ⟨φs⟩ , where Q∗ denotes the Lax–Friedrichs numerical flux across the interface, defined as Q∗ = {Q}+n ⌈ c⌉⌉ JηK. Expressed in the parent element, this yields the discrete form of the free-surface equation: ⟨φJh ∂t η⟩T̂2D = ⟨Jh ∇h φ · Q⟩T̂2D − ⟨⟨φ(n · {Q} + ⌈ c⌉⌉ JηK)Jl ⟩⟩ L̂ + ⟨φsJh ⟩T̂2D

(2)

S1.2

Depth-averaged momentum

over the m iterations of the external mode performed during the current step of the internal mode. While F2D can build similarly, it is more efficient to derive it from the 2D momentum before and after the 2D iterations (Q0 and Q1 ).

The depth-averaged momentum equation (3) constitutes the second component of the external mode: H ∇h patm + S3 ρ0 H2 H = −g∇h − gH∇h b − ∇h patm + S3 2 ρ0

∂t Q = −gH∇h η −

(3)

F2D =

1 m

h ⟨φ∂t QJh ⟩T̂2D − F3D→2D

2D iterations

(6)

h = (Q1 − (Q0 + ∆tF3D→2D ))/∆t

Neglecting the source terms from the 3D mode and the atmospheric pressure for now, the weak form becomes

⟨φ∂t Q⟩ = g(H 2 /2)∇h φ − nφg(H 2 /2)∗ − ⟨φgH∇h b⟩

3D Diagnostic Equations

This section focuses on the two main diagnostic equations of the 3D model: the continuity equation (7) and the equation for the horizontal pressure gradient (10).

where (H 2 /2)∗ is also evaluated using a Lax–Friedrichs flux: (H 2 /2)∗ = H 2 /2 + n · ⌈ c⌉⌉ JQK . This leads to the discrete equation:

S2.1

Continuity Equation

⟨φ∂t Q⟩ = g(H 2 /2)∇h φ − nφg H 2 /2

The continuity equation (7) is solved on the 3D mesh, which is slightly non-conformal. In general, a lateral interface between two prisms does not coincide exactly on However, this formulation involves subtracting two both sides, creating an ambiguity in the integration doO(H 2 ) terms, introducing numerical errors of order main for the lateral terms. To avoid this, all integrals are O(H 2 εmachine ). This is excessive since the total flux is expressed in the parent element, and variations in the inonly O(Hη), and H is typically several orders of magnitegration domain are represented through a discontinuous tude η. qTo mitigate this, we rewrite H 2 as 2 larger than y Jacobian across the lateral interfaces. H = (H 2 )int − H 2 and reverse the integration by parts, yielding ∂z w = −∇h · u (7) q 2 y ⟨φ∂t Q⟩ = − ⟨gH∇h η⟩ + nφg H /2 − ⟨⟨⌈⌈c⌉⌉ JQK⟩⟩ . The discrete form of equation (7) is based on the folq 2y lowing weak formulation: Since H represents a difference of two squares, we can − ⟨φgH∇h b⟩ − ⟨⟨⌈⌈c⌉⌉ JQK⟩⟩ .

express it as

⟨⟨φnz w∗ ⟩⟩ T − ⟨w∂z φ⟩P = ⟨u · ∇h φ⟩P − ⟨⟨φnh · u∗ ⟩⟩ T

q 2 y (H int + H ext )(η int − η ext ) = 2 {H} JηK H = 2

− ⟨⟨φnh · u∗ ⟩⟩ Q ,

which, when expressed in the parent element, becomes which significantly reduces numerical noise, yielding an error on the order of O(Hηεmachine ). This “reverse integra⟨⟨φnz w∗ JT ⟩⟩ T̂ − ⟨w∂z φJz Jh ⟩P̂ tion by parts” procedure will be reused for other equations, = ⟨Jz u · ∇h φJh ⟩P̂ − ⟨⟨JT φnh · u∗ ⟩⟩ T̂ though it is most impactful here. − ⟨⟨φnh · u∗ Jz Jl ⟩⟩ Q̂ . Reintroducing the atmospheric pressure and the 3Dh mode contribution F3D→2D (the discrete equivalent of S3 , Here, we use the tensor-product structure of the mesh, discussed later), we obtain where the volume Jacobian satisfies J = Jz Jh and the ⟨φ∂t Q⟩ = − ⟨gφH∇h η⟩ + ⟨⟨nφg {H} JηK⟩⟩ lateral Jacobian JQ = Jz Jl . A key requirement of the continuity equation is its conh − ⟨⟨⌈⌈c⌉⌉ JQK⟩⟩ − ⟨φH/ρ0 ∇h patm ⟩ + F3D→2D . sistency with the free-surface equation (1) in some parts Finally, expressing the integrals in the parent triangle of the temporal scheme. In practice, this means that the yields the final discrete 2D momentum equation: vertical integral of (7), combined with the impermeability condition at the ocean floor, must reproduce the free⟨φ∂t QJh ⟩T̂2D = − ⟨gφH∇h ηJh ⟩T̂2D surface equation (1). This property must also hold in the + ⟨⟨nφg {H} JηK Jl ⟩⟩ L̂ discrete formulation. To ensure this, the discrete continuity equation employs a linearized transport q, consis− ⟨⟨⌈⌈c⌉⌉ JQK Jl ⟩⟩ L̂ tent with the averaged 2D linear transport Q̄ from equah − ⟨φH/ρ0 ∇h patm Jh ⟩T̂2D + F3D→2D . tion (5). (4) In the continuous limit, we have q = Jz u. In practice, however, q is always projected onto the linear basis When coupling the internal and external modes in the functions from Jz u. When consistency is required, the 3D temporal scheme, the average of the 2D fluxes (F2D ) in (4) transport is build to match the mean 2D transport Q̄ and is required, together with and their corresponding transdenoted as q̄. port Q. Both are used an input for the vertical part of X q̄ = Q̄ the 3D time step. Thus, during the iterations of the 2D Vertical DOFs mode, the average transport Q̄ is accumulated as X 1 Q̄ = Q (5) m 2D iterations

— In the term ⟨⟨φnz r ∗ JT ⟩⟩ T̂ , r ∗ refers to the value of r above the interface. Thus, at the top of an element we use r from the exterior, while at the bottom we use r from the interior. As in equation (9), this choice follows from the vertical nature of the equation and the boundary condition being applied at the ocean surface.

Otherwise, q is simply the projection of Jz u to the linear basis functions. To express the continuity equation in terms of q, we replace Jz u by q wherever possible, and use q/Jz in place of u otherwise: ⟨⟨φnz w∗ JT ⟩⟩ T̂ − ⟨w∂z φJz Jh ⟩P̂ = ⟨q · ∇h φJh ⟩P̂ − ⟨⟨φnh · (q/Jz )∗ JT ⟩⟩ T̂ − ⟨⟨φnh · q ∗ Jl ⟩⟩ Q̂ .

— In the term ⟨⟨φnh ρ′∗ JT ⟩⟩ T̂ , ρ′∗ also refers to the value above the interface, for the same reason.

This equation can be simplified by noting that JT = |Jh /nz |

∂ẑ =

and

∂z ∂z = Jz ∂z , ∂ ẑ

— In ⟨⟨φnρ′∗ Jz∗ Jl ⟩⟩ Q̂ , we discretize ρ′∗ Jz∗ as {ρ′ } {Jz }. We do not use {ρ′ Jz }, since for constant density fields with slightly discontinuous interfaces, this alternative would spuriously modify the flow through artificial density gradients.

(8)

which gives ⟨⟨φn̂z w∗ Jh ⟩⟩ T̂ − ⟨w∂ẑ φJh ⟩P̂ = ⟨q · ∇h φJh ⟩P̂ − ⟨⟨φnh · (q/Jz )∗ |Jh /nz |⟩⟩ T̂

With these definitions and applying the simplifications from equation (8), we rewrite the equation as

− ⟨⟨φnh · q ∗ Jl ⟩⟩ Q̂ . We now define the interface variables explicitly: ∗

φr ext n̂z Jh

∗

T̂top ′

+ φr int n̂z Jh

− ⟨r∂ẑ (φ) Jh ⟩P̂

T̂bot ′ext

— In the term ⟨⟨φn̂z w Jh ⟩⟩ T̂ , w corresponds to the = −g ⟨∇h φρ Jh Jz ⟩P̂ + g φnh ρ |Jh /nz | T̂top value of w below the interface. Thus, at the top of + g φnh ρ′int |Jh /nz | T̂ + g ⟨⟨φn {ρ′ } {Jz } Jl ⟩⟩ Q̂ . an element, w is taken from the interior, while at the bot bottom it is taken from the exterior. This choice is consistent with the vertical nature of the equation and As in Section S1.2, we integrate the right-hand side backwards, leading to the impermeability condition at the ocean floor. + φr int n̂z Jh T̂ − ⟨r∂ẑ (φ) Jh ⟩P̂ T̂top bot ′ ′ext = g ⟨φ∇h ρ Jh Jz ⟩P̂ + g φnh (ρ − ρ′int )|Jh /nz | T̂

φr ext n̂z Jh

— In the term ⟨⟨φnh · (q/Jz )∗ |Jh /nz |⟩⟩ T̂ , the interface value is computed as the simple average (q/Jz )∗ = {q/Jz }.

top

+ g φn({ρ′ } − ρ′int ) {Jz } Jl

— The term ⟨⟨φnh · q ∗ Jl ⟩⟩ Q̂ requires particular care to ensure consistency with the boundary terms of equation (4). To achieve this, we define the lateral flux as q ∗ = {q} + nh {Jz /H} ⌈ c⌉⌉ JηK .

The integral over the bottom of each element cancels out during this backward integration. Since n̂z appears only on the top face, it simplifies to n̂z = 1. Finally, using (ρ′ext − ρ′int ) = −2 Jρ′ K and ({ρ′ } − ρ′int ) = − Jρ′ K, we obtain the final discrete form of equation (10):

This yields the final discrete continuity equation: φwint n̂z Jh

T̂top

+ φwext n̂z Jh

T̂bot

+ φr int Jh T̂ − ⟨r∂ẑ (φ) Jh ⟩P̂ T̂top bot = g ⟨φ∇h ρ′ Jh Jz ⟩P̂ − g ⟨⟨2φnh Jρ′ K |Jh /nz |⟩⟩ T̂top

φr ext Jh

− ⟨wJh ∂ẑ φ⟩P̂

= ⟨q · ∇h φJh ⟩P̂ − ⟨⟨nh · {q/Jz } φ|Jh /nz |⟩⟩ T̂

. Q̂

(9)

(11)

′

− g ⟨⟨φnh Jρ K {Jz } Jl ⟩⟩ Q̂ .

− ⟨⟨φ (nh · {q} + {Jz /H} ⌈ c⌉⌉ JηK) Jl ⟩⟩ Q̂ .

As mentioned previously, q should be understood as the S3 Prognostic Equations linearized 3D horizontal transport Jz u, i.e., a projection of the exact (second-degree) transport field onto the linear In the 2D equations, only the spatial components are basis functions. discretized, resulting in a system of ordinary differential equations (ODEs) in time. This system can then be integrated using any of a wide range of implicit or explicit S2.2 Horizontal Pressure Gradient time-stepping methods. In contrast, the discretization of The continuous equation for the horizontal pressure gra- the 3D equations follows the temporal scheme introduced in (Ishimwe et al., 2025), which couples the spatial and dient r reads ′ ∂z r = g∇h ρ , (10) temporal discretizations while treating the vertical terms implicitly in time. with its corresponding weak form Compared to the original formulation in (Ishimwe et al., 2025), this work adopts a slightly modified im∗ ′ ′∗ ⟨⟨φnz r ⟩⟩ − ⟨r∂z φ⟩ = −g ⟨∇h φρ ⟩ + g ⟨⟨φnh ρ ⟩⟩ , plicit/explicit splitting strategy. Rather than decomposing the 3D velocity field v into horizontal (u) and vertical which, when expressed in the parent element, becomes (w) components, we split it into a mesh-aligned compo⟨⟨φnz r ∗ JT ⟩⟩ T̂ − ⟨r∂z (φ) Jh Jz ⟩P̂ nent ũ and a vertical component w̃, as illustrated in Figure S1. The vertical component w̃ is treated implicitly, = − g ⟨∇h φρ′ Jh Jz ⟩P̂ + g ⟨⟨φnh ρ′∗ JT ⟩⟩ T̂ while the mesh-aligned component ũ is handled explicitly. + g ⟨⟨φnρ′∗ J ∗ Jl ⟩⟩ . z

Q̂

On the interfaces, the discontinuous variables are defined as follows: 3

Because ũ is aligned with the mesh, it is orthogonal to the normal vector on the top and bottom faces. Consequently, the integrals over T̂ vanish. On the lateral faces, the normal vector lies in the horizontal plane, so n · ũ = nh · u, resulting in the same term as in the standard continuity u equation. ũ To simplify the volume term, we decompose the gradient 1 operator as ˜h + m ∂ , ∇=∇ Figure S1: Left: original splitting, where the horizontal ∂ζ velocity u is treated explicitly and the vertical velocity w implicitly. Right: modified splitting used in this work, where where the mesh-aligned velocity component ũ is treated ∂ξ ∂ ∂η ∂ ∂ξ ∂ ∂η ∂ ˜ · + , + ,0 . ∇h ·= explicitly and the vertical component w̃ implicitly. ∂x ∂ξ ∂x ∂η ∂y ∂ξ ∂y ∂η w

w̃

˜ h represents the horizontal gradient along The operator ∇ iso-ζ slices of a prism. For 2D fields, or for fields that ˜ h reduces to ∇h . Applying this are vertically constant, ∇ decomposition to ∇φ yields

Note that this does not introduce any additional spatial error. To explicitly define ũ in terms of u, let us introduce the vector m as ∂ζ ∂ζ ∂ζ ∂ζ = , , , m= ∂x ∂x ∂y ∂z

˜ h φ + m ∂φ = φz ∇h φh + m ∂φ . ∇φ = ∇ ∂ζ ∂ζ

where ζ denotes the vertical parametric coordinate inside Since ũ is perpendicular to m, the second term from ũ·∇φ ˜ hφ = u · ∇ ˜ h φ. Substituting this the prism, ranging from ζ = −1 (bottom) to ζ = 1 (top). vanishes, leaving ũ · ∇ By construction, m is aligned with the normal vector on result, the integrals simplify to the top and bottom faces of each element, and orthogonal ⟨⟨φnz w̃∗ JT ⟩⟩ T̂ − ⟨w̃∂z φJz Jh ⟩P̂ to ũ. We can then express = ⟨Jh Jz u · φz ∇h φh ⟩P̂ − ⟨⟨φnh · u∗ Jz Jl ⟩⟩ Q̂ . ũ = (u, −mh · u/mz ) = (u, v, −mh · u/mz ), From this point, the derivation proceeds identically to that leading to equation (9), resulting in the following discrete form for w̃:

while the corresponding vertical component is w̃ = w + mh · u/mz .

φw̃int n̂z Jh

The purpose of this decomposition is to minimize the explicit fluxes across the top and bottom faces. By construction, ũ is tangent to these faces, thereby eliminating explicit advective fluxes. To further reduce diffusive and viscous contributions, a similar splitting is applied to the diffusivity and viscosity tensors. Let D denote either of these tensors. We decompose it as D = Di + De , where Di and De represent the implicit and explicit parts, respectively. The implicit component is defined as   0 m·D·m  . 0 Di = (12) ez ⊗ ez =  m2z Di

T̂top

+ φw̃ext n̂z Jh

T̂bot

− ⟨w̃Jh ∂ẑ φ⟩P̂

= ⟨Jh q · φz ∇h φh ⟩P̂ − ⟨⟨φ (nh · {q} + {Jz /H} ⌈ c⌉⌉ JηK) Jl ⟩⟩ Q̂ .

(13)

The modified continuity equation thus provides a consistent definition of w̃ that is fully compatible with the implicit/explicit velocity decomposition introduced earlier. This expression of w̃ is then used in the implicit treatment of the vertical fluxes within the 3D momentum and tracer equations, discussed below.

Since De = D − Di , this ensures that m · De · m = 0, S3.2 3D Momentum and thus n · De · n = 0 on the top and bottom faces. The discretization of equation (14) introduces additional complexity due to its two-way coupling with the external mode: S3.1 Modified Continuity Equation ∂t u + ∇h · (u ⊗ u) + ∂z (wu)

Although w̃ can, in principle, be computed from w, u, and m, it is more efficient and numerically stable to solve for w̃ directly. The governing equation for w̃ has the same form as equation (7):

= ∇h · (κh ∇h u) + ∂z (κv ∂z u) 1 − f ez × u − r + S2 ρ0

∂z w̃ = −∇ · ũ .

(14)

To handle this coupled system, SLIM reformulates the discrete version of (14) as

Following the same procedure as before, we express its weak form in the parent domain as

M1 u 1 − M 0 u 0 M1 F2D h v = F3D (u, q̄, r) + + F3D (u, w̃, u1 ) ∆t H1 (15)

⟨⟨φnz w̃∗ JT ⟩⟩ T̂ − ⟨w̃∂z φJz Jh ⟩P̂ = ⟨Jz ũ · ∇φJh ⟩P̂ − ⟨⟨JT φn · ũ∗ ⟩⟩ T̂ − ⟨⟨φn · ũ∗ Jz Jl ⟩⟩ Q̂ . 4

h The discrete form of F3D is

when vertical quantities are treated implicitly, or M1 u1 − M0 u0 M1 F2D v h = F3D (u, q̄, r) + + F3D (u, w̃, u) ∆t H1 (16) when everything is explicit. Here, M0 and M1 are the (time-varying) mass matrices at the beginning and end of the time step, respectively. u0 is the velocity field at the beginning of the step, and u1 is the unknown at its end. u represents the intermediate solution from the Runge–Kutta scheme. The transport q̄ is derived from u h and the mean transport Q̄. The term F3D includes the horizontal advection and viscosity, the Coriolis force and v the horizontal pressure gradient. F3D accounts for the vertical advection and viscosity, and wind stresses and bottom friction as viscous boundary conditions. Finally, F2D represents the horizontal momentum change from the external mode, i.e. the change induced by S2 and computed via (6). Since F2D is vertically integrated, a division by H is required to obtain a velocity increment. In practice, H is taken at the end of the step (H1 ). The multiplication by the mass matrix ensures that the velocity increment is uniform vertically. With this splitting, the interdependence between (14) and (3) becomes manageable. The 3D contribution to the 2D mode (S3 ) corresponds to the vertical integral of (14) without S2 , which discretely translates to the sum over all v h : and F3D vertical degrees of freedom of F3D X h h v F3D→2D = (F3D + F3D ).

h F3D (u, q, r) = ⟨Jh u(q · φz ∇h φh )⟩P̂

− ⟨⟨φuup (nh · {q} + {Jz /H} ⌈ c⌉⌉ JηK) Jl ⟩⟩ Q̂ − ⟨J(∇φ · κe · ∇)u⟩P̂

+ ⟨⟨φ {J∂ (n · κe · ∇)u}⟩⟩ ∂ P̂ − ⟨⟨σ3 ⌈ n · κe · n⌉⌉ {JQ } JuK⟩⟩ Q̂

− ⟨Jφf ez × u⟩P̂ − ⟨Jφr/ρ0 ⟩P̂ . (17) v The corresponding expression for F3D is v F3D (u, w̃, u1 ) = ⟨J∂z (φ) (w̃ − wm )u1 ⟩P̂ bot − Jh n̂z φuup − wm ) 1 (w̃

T̂

− ⟨J∂z (φ) κi ∂z (u1 )⟩P̂

(18)

+ ⟨⟨φ {JT nz κi ∂z u1 }⟩⟩ T̂ − σ3 κi n2z JT Ju1 K T̂ . In both equations, the viscosity tensor κ is parameterized using a horizontal viscosity κh based on Smagorinsky’s formulation, and a vertical viscosity κv provided by the turbulence closure model. The complete tensor κ is split into explicit and implicit parts following (12): In the interface integrals, the upwind direction is determined from the advection velocity immediately adjacent to uup . The mesh used in (17) and (18) always corresponds to that of u, meaning u and the mesh come from the same Runge–Kutta stage. The interior penalty parameter (Rivière, 2008), (Riviere & Sardar, 2014) σ3 is defined as N0 (o + 1)(o + d) (19) σd = 2d · min(Lint , Lext )

Vertical DOFs

However, neglecting wind drag and bottom friction for v integrates to zero over the vertical and the moment, F3D h is therefore not required for computing F3D→2D . This v allows equation (4) to be solved before F3D is known. Wind drag and bottom friction must still be included in h . During explicit steps, their contribution can be F3D→2D added directly. For implicit steps, however, the bottom drag depends on the unknown velocity u1 and is therefore approximated using a prediction, rather than computed exactly. Wind drag, on the other hand, is treated explich h is itly. Since consistency is not relevant for F3D→2D , F3D evaluated with q rather than q̄. With this strategy, the discretization of (14) proceeds as follows:

where d = 3 is the problem dimension, o = 1 is the polynomial degree of the basis functions, and N0 = 5 is the average number of neighboring elements. L approximates the effective length scale of an element. On horizontal interfaces, L is taken as the average prism height, and on lateral faces, L = A/l with A the triangle area and l the edge length. Although originally derived for simplices, this definition performs well in practice for prisms.

S3.3

Tracers

1. Solve for r using (11).

The discrete tracer equation is formulated on the same 2. Compute q by projecting Jz u to the linear basis func- moving mesh as the momentum equation. Its discretization of advection and diffusion terms follows an almost tions. identical structure, differing only by the absence of source h 3. Compute F3D (u, q, r). terms such as the external mode, the wind, or the bottom 4. Compute the wind drag τ and compute or predict friction. The resulting discrete form reads: s

the bottom friction τb∗ . P h h 5. Compute F3D→2D = Vertical DOFs F3D + τs + τb∗ .

⟨φT1 ⟩P1 − ⟨φT0 ⟩P0 = ⟨Jh T (q̄ · φz ∇h φh )⟩P̂ ∆t up − ⟨⟨φT (nh · {q̄} + {Jz /H} ⌈ c⌉⌉ JηK) Jl ⟩⟩ Q̂

6. Advance (4), accumulating F2D and Q̄ following (6) and (5).

− ⟨J(∇φ · νe · ∇)T ⟩P̂ + ⟨⟨φ {J∂ (n · νe · ∇)T }⟩⟩ ∂ P̂ − ⟨⟨σ3 ⌈ n · νe · n⌉⌉ {JQ } JT K⟩⟩ Q̂

7. Update the vertical mesh coordinate (z) and velocity (wm ).

+ ⟨J∂z (φ) (w̃ − wm )T1 ⟩P̂ − Jh n̂z φT1up (w̃bot − wm ) − ⟨J∂z (φ) νi ∂z (T1 )⟩P̂ + ⟨⟨φ {JT nz νi ∂z T1 }⟩⟩ T̂ − σ3 νi n2z JT JT1 K T̂ .

8. Compute q̄ from Jz u and Q̄. 9. Compute the modified vertical velocity w̃ from Q̄. 10. Solve equation (15) or (16) for u1 . 5

(20)

T̂

Although the corrected transport q̄ is strictly required ocean models. Journal of Computational Physics, only for tracer advection, it is also applied to the advec520, 113482. https://doi.org/10.1016/j.jcp. tion of momentum. Since q̄ must be computed for tracers 2024.113482 in any case, using it consistently for velocity ensures a uni- Rivière, B. (2008, January). Discontinuous galerkin methfied and coherent treatment of advection with acceptable ods for solving elliptic and parabolic equations. Soadditional costs. ciety for Industrial; Applied Mathematics. https: //doi.org/10.1137/1.9780898717440 Riviere, B., & Sardar, S. (2014). Penalty-free discontinReferences uous galerkin methods for incompressible navier– stokes equations. Mathematical Models and MethIshimwe, A. P., Deleersnijder, E., Legat, V., & Lamods in Applied Sciences, 24 (6), 1217–1236. https: brechts, J. (2025). A multi-scale imex second//doi.org/10.1142/S0218202513500826 order runge–kutta method for 3d hydrodynamic

Record · ID 194284 · SHA-256 b92de89bea606ae8
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.