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

This paper presents a highly efficient, multi-GPU-accelerated implementation of the Discontinuous Galerkin ocean model SLIM that achieves massive speedups over CPU-based systems and enables ultra-high-resolution coastal simulations, such as a five-fold resolution improvement for the Great Barrier Reef.

The Gist

The Big Picture: Making Ocean Models "Super-Fast"

Imagine trying to simulate the ocean. For a long time, scientists used a "grid" like a chessboard to map the water. But the ocean isn't a chessboard; it has jagged coastlines, deep trenches, and shallow reefs. To make the chessboard fit, you either have to make the squares tiny everywhere (which takes forever to calculate) or accept that the edges look blocky and wrong.

The SLIM model described in this paper uses a different approach: an unstructured mesh. Think of this like a mosaic made of irregularly shaped tiles. You can use tiny, intricate tiles right next to a rocky reef and huge, simple tiles in the deep, open ocean. This is perfect for coastal areas, but it's computationally expensive. It's like trying to paint a masterpiece with a tiny brush; it takes a lot of time and effort.

The authors of this paper asked: "How can we make this detailed, mosaic-style ocean model run fast enough to be useful?" Their answer was to build a version specifically designed for GPUs (the powerful graphics chips found in gaming computers and supercomputers).

The Core Innovation: The "GPU-Ready" Ocean

The paper focuses on a specific mathematical method called Discontinuous Galerkin (DG).

• The Analogy: Imagine a classroom.
  • Old methods (Continuous): The students are holding hands in a giant circle. If one student moves, they have to tell everyone else in the circle. It's connected, but slow to coordinate.
  • DG Method: Each student sits at their own desk. They work independently on their own math problems. They only talk to their immediate neighbors when they need to pass a note.
• Why this helps: Because the students (data points) work independently, you can hire 1,000 teachers (GPU cores) to help them all at the same time without them getting in each other's way. This is exactly what GPUs love to do: massive parallel work.

How They Made It Fast (The "Secret Sauce")

The authors didn't just put the code on a GPU; they completely redesigned how the data is stored and moved, using three main tricks:

1. The "Library" Organization (Memory Layout)
GPUs are like super-fast librarians. If books are scattered randomly, the librarian wastes time running around. If they are organized perfectly, they can grab them instantly.

• The team reorganized the data so that related information sits right next to each other in memory. They even used a "Hilbert curve" (a specific winding path) to arrange the irregular tiles so that neighbors are physically close in the computer's memory. This keeps the GPU's "librarian" running at top speed.

2. The "Cell" Assembly Line
The ocean model is 3D, made of vertical columns of water. Some calculations need to solve a puzzle for the whole column at once.

• The Problem: Usually, solving these puzzles one by one is slow.
• The Fix: They created a special "Cell" layout. Imagine a factory assembly line where 128 workers (threads) are assigned to 128 columns. Instead of passing parts back and forth, they organize the parts into a neat grid (a matrix) so all 128 workers can grab what they need simultaneously. This turns a slow, sequential process into a fast, parallel one.

3. The "No-Blueprint" Solver (Matrix-Free)
In many math problems, you have to build a giant blueprint (a matrix) before you can solve the problem. Building the blueprint takes time.

• The Trick: For certain parts of the ocean model (like pressure and vertical movement), the authors realized the blueprint always follows a predictable pattern. Instead of building the blueprint, they wrote a recipe that calculates the answer directly on the fly. It's like knowing the answer to a math problem without needing to write out the long division steps.

The Results: A Speed Revolution

The paper presents benchmark results that show just how effective this is:

• One GPU vs. A Room of Computers: A single high-end GPU (like an NVIDIA A100) can do the work of about 1,500 standard computer processors.
• The "50x" Leap: If you replace a massive server with 128 CPU cores with a single server containing just 4 of these GPUs, the simulation runs 50 times faster.
• Scaling Up: They tested this on supercomputers with up to 1,024 GPUs. The system scaled beautifully, meaning adding more GPUs kept the simulation running efficiently, provided the ocean area being simulated was large enough to keep all those GPUs busy.

The Real-World Test: The Great Barrier Reef

To prove this wasn't just a theoretical speed test, they ran a simulation of the Great Barrier Reef.

• The Challenge: The reef has incredibly complex shapes. Previous models had to use a "blurry" resolution (about 1.5 km to 4 km per tile) to run in a reasonable time.
• The New Result: Using their new GPU-accelerated model, they simulated the entire reef with a resolution five times finer (down to 200 meters).
• The Outcome: They could see tiny details like "tidal jets" (fast streams of water) and small eddies that were previously invisible. They achieved a speed where the computer simulated 100 days of ocean time for every 1 day of real time.

Summary

This paper demonstrates that by rethinking how data is organized and leveraging the unique power of modern graphics chips, scientists can finally run highly detailed, 3D ocean models of complex coastlines. They turned a process that used to be too slow and expensive into a fast, efficient tool, opening the door to ultra-high-resolution simulations of places like the Great Barrier Reef.

Technical Summary

Technical Summary: An Efficient Multi-GPU Implementation for the Discontinuous Galerkin Ocean Model SLIM

Problem Statement
Unstructured-mesh ocean models are increasingly vital for coastal applications due to their ability to represent complex geometries and apply local grid refinement. However, their broader adoption, particularly for models based on the Discontinuous Galerkin finite element (DG-FE) method, has been hindered by high computational costs. The DG-FE approach inherently involves more degrees of freedom and calculations per element than traditional finite volume or continuous finite element methods. Furthermore, the irregular connectivity of unstructured meshes necessitates explicit storage and management of neighbor relationships, leading to sophisticated data handling requirements and increased computational overhead. While strategies like local mesh refinement and mode-splitting have been employed, they are often insufficient for high-resolution or long-duration simulations. The rapid emergence of GPU-based high-performance computing (HPC) architectures offers a potential solution, as DG-FE formulations are well-suited for massively parallel, element-wise computations, yet fully leveraging these architectures requires a complete redesign of computational codes.

Methodology
The authors present a full 3D DG-FE ocean model implementation of SLIM, optimized for both single- and multi-GPU systems with support for NVIDIA and AMD architectures. The implementation relies on a lightweight abstraction layer to switch between CPU, CUDA, and HIP backends.

• Memory Layout and Thread Assignment: The model employs a "one-thread-per-element" strategy rather than per-node to increase computational workload per thread. To optimize memory coalescence, the 2D mesh is reordered using a Hilbert curve. A hybrid data layout is introduced: a Structure-of-Arrays (SoA) format for general computations and a specialized "cell layout" for solving linear systems. The cell layout groups sets of prism columns (typically 128) into a matrix structure (Array-of-Structure-of-Arrays), enabling perfectly coalesced memory access during the solution of banded linear systems associated with vertical columns.
• Numerical Schemes: The model utilizes a split-IMEX Runge-Kutta scheme. The external (barotropic) mode uses a three-step explicit Runge-Kutta method on the 2D mesh, while the internal (baroclinic) mode advances the 3D equations. The internal mode separates horizontal and vertical processes; horizontal advection and viscosity are explicit, while vertical processes can be implicit to mitigate strict CFL conditions.
• Solvers:
  • Matrix-Free Solvers: For the horizontal pressure gradient and vertical velocity, the linear systems possess a predictable structure that allows for solution without explicit matrix assembly. A single-pass up-looking solver encodes the matrix structure in a recursion formula.
  • Fully-Assembled Column Solvers: For turbulence closure and implicit vertical fluxes, the model assembles sparse matrices coupling nodes within a column. These are solved in-place using Gaussian elimination with one thread per system, utilizing local buffers to avoid register spilling.
• Multi-GPU Strategy: The domain is decomposed horizontally using MPI, with one GPU per rank. To overlap computation and communication, the code splits each partition into boundary/ghost elements and interior elements. Boundary computations and halo exchanges are launched on a high-priority communication stream, while interior computations proceed on a compute stream. This overlap is critical, particularly for the 2D external mode where kernel launch latency and communication costs can become bottlenecks.

Key Results

• Single-GPU Performance: Benchmarks on an NVIDIA A100 GPU demonstrate performance equivalent to approximately 1500 CPU cores. Replacing a 128-core CPU node with a 4×A100 GPU node yields a speedup of roughly 50×. The implementation sustains up to 80% of peak memory bandwidth for memory-bound kernels and 60% of peak floating-point throughput for compute-bound kernels, with a sustained average of ~30% for both over a complete time step.
• Scaling: Weak-scaling efficiency is maintained up to 1024 GPUs (512 AMD MI250X GCDs) on the LUMI supercomputer. Strong scaling follows Amdahl's law, where the 3D baroclinic component scales nearly ideally, but the 2D barotropic mode acts as a sequential fraction due to its high frequency of short kernels and halo exchanges.
• Application to the Great Barrier Reef (GBR): The model was applied to a realistic GBR simulation with a mesh of 3.3 million triangles (34 million 3D elements), featuring resolutions from 200 m near reefs to 10 km offshore. Running on 32 AMD MI250X GPUs (64 GCDs) in double precision, the simulation achieved a throughput of approximately 100 simulated days per wall-clock day. This configuration resolved flow structures down to a few hundred meters, capturing tidal jets and eddies previously unresolvable in uniform-resolution models.

Significance and Claims
The paper claims that GPU-accelerated DG-FE methods can dramatically advance the capabilities of unstructured-mesh ocean modeling. By effectively exploiting GPU architectures through specific data layouts, matrix-free solvers, and communication overlap strategies, the authors demonstrate that the computational cost of high-resolution unstructured simulations can be brought in line with structured-grid models at coarser resolutions. This removes a historical constraint on the practical application of unstructured-mesh models for regional and coastal domains.

The authors emphasize that this work enables ultra-high-resolution coastal simulations that were previously infeasible, allowing for the resolution of multi-scale processes along the land-sea continuum. While the GBR application serves as a demonstration of scalability and resolution capabilities rather than a fully validated scientific study, it highlights the potential for operational applications in environmental management. The paper concludes that the inherent advantages of unstructured meshes—geometrical flexibility, local refinement, and multi-scale resolution—can now be leveraged without the traditional computational overhead, positioning GPU-accelerated DG models as a promising tool for the next generation of coastal ocean simulations.