# Numerics-driven uplifting, Automatic Parallelization, and Performance Optimizations with Deep Kernel Fusion for Ocean Models on Heterogeneous Architectures

Julien Remy, Hugo Brunie, Martin Schreiber, Sergi Siso, Andrew Porter, Rupert Ford, Laurent Debreu, Florian Lemarié

► To cite this version:

Julien Remy, Hugo Brunie, Martin Schreiber, Sergi Siso, Andrew Porter, et al. Numerics-driven uplifting, Automatic Parallelization, and Performance Optimizations with Deep Kernel Fusion for Ocean Models on Heterogeneous Architectures. 2025. hal-05061001

HAL Id: hal-05061001

https://hal.science/hal-05061001v1

Preprint submitted on 8 May 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerics-driven uplifting, Automatic Parallelization, and Performance Optimizations with Deep Kernel Fusion for Ocean Models on Heterogeneous Architectures

Julien Remy

Hugo Brunie

Martin Schreiber

julien.remy@univ-grenoble-alpes.fr

hugo.brunie@univ-grenoble-alpes.fr

martin.schreiber@univ-grenoble-alpes.fr

Univ. Grenoble Alpes, Inria, CNRS, LJK, Grenoble INP, Grenoble, France

Sergi Siso

Andrew Porter

Rupert Ford†

sergi.siso@stfc.ac.uk

andrew.porter@stfc.ac.uk

† Post-mortem authorship

STFC - Hartree Centre, Daresbury Laboratory, United Kingdom

Laurent Debreu

Florian Lemarié

Arthur Vidard

Laurent.Debreu@inria.fr

Florian.Lemarie@inria.fr

Arthur.Vidard@inria.fr

Univ. Grenoble Alpes / Inria / CNRS / Lab. Jean Kuntzmann, Grenoble, France

## Abstract

Achieving peak performance for ocean models on modern High-Performance Computing (HPC) architectures requires extensive, costly code rewrites that are not only time-consuming and erroneous but also highly architecture-specific and require numerics experts to be proficient in parallel programming models or Domain-Specific Languages (DSLs).

In this paper we introduce Poseidon, a source-to-source code optimization tool that employs numerics & HPC co-design. We developed an uplifting approach to a hypergraph for the data flow and to our own Poseidon Intermediate Representation (PosIR) for the computations, which recovers for the first time high-level information and semantics about the computations and memory management that are typically lost during the conversion of numerical algorithms to source code. This representation is then employed for model-driven optimizations. In the backend, we inject source code back into the original code supporting different target HPC architectures, reusing the existing model as the runtime.

Our evaluation investigates OpenMP and OpenACC parallel programming models on CPUs and GPUs. We demonstrate Poseidon's capability of automatic parallelization and optimization of existing Fortran code, with comparisons to state-of-the-art parallelized code with performance improvements of up to ×2.60 and automatically reducing the required memory footprint by up to ×2.29.

## Keywords

Language uplifting, source-to-source compilation, deep kernel fusion, loop fusion, model-driven performance optimization

## ACM Reference Format:

Julien Remy, Hugo Brunie, Martin Schreiber, Sergi Siso, Andrew Porter, Rupert Ford†, Laurent Debreu, Florian Lemarié, and Arthur Vidard. 2025. Numerics-driven uplifting, Automatic Parallelization, and Performance Optimizations with Deep Kernel Fusion for Ocean Models on Heterogeneous Architectures. In *Proceedings of *-. ACM, New York, NY, USA, 13 pages. https://doi.org/XXXXXXX.XXXXXXX

## 1 Introduction

Ocean simulations rank among the most computationally demanding codes for extreme-scale computers. Achieving high efficiency in these simulations is crucial for reducing the carbon footprint associated with their execution, lowering procurement costs, and enabling higher-resolution simulations. Higher resolution, in turn, leads to more accurate weather forecasts [17] and climate reports [11].

These numerical models have been developed over decades and are predominantly developed in (modern) Fortran. The source code represents a vast repository of accumulated knowledge, which is the true asset of these developments. Porting such extensive and intricate codebases to other DSLs would require substantial financial investments, necessitate convincing a large community (often exceeding 1000 users) to adopt a different language, and carry the significant risk of introducing errors due to a lack of understanding of the existing code. These challenges make such a transition highly complex and potentially disadvantageous.

A compromise is often taken using code annotations with OpenMP or OpenACC. However, such code annotations often disturb unfamiliar developers, bugs can be introduced by model developers due to code that is not valid under existing annotations, and no standard best-performing parallel programming model exists for all architectures. Instead, source-to-source compilers such as Rose[32] and PSyclone [2, 39] could be used, but they have similar limitations to existing compilers.

One of the biggest challenges for automatic tools is that transferring the discrete mathematical equations into source code leads to a *loss of knowledge on the meaning of variables, loops, computations, etc.*, hindering the compiler from performing optimizations that would be possible once this deeper understanding could be

<footnote>Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org.</footnote>

<footnote>-, -</footnote>
regained. A brief discussion on this for ocean simulations on struc-
tured grids – the main focus of the present work – will be provided
in Sec. 2. This is followed by introducing for the first time Poseidon,
enabling us to finally reduce the gap between the high-level math-
ematical equations and the lower-level source code in Sec. 3 - 7
with results given in Sec. 8. Although our radical source-to-source
approach with deeper optimizations is novel, some concepts are
partly similar to related developments discussed in Sec. 9. Finally,
our work should not be seen as some free-lunch source-to-source
performance optimizations, but as a novel methodology to approach
portability of legacy code of which we will provide an overview of
pros and cons in Sec. 10 followed by a summary in Sec. 11.

## 2 Ocean models

This section provides the necessary background on ocean models on structured grids for an improved understanding of the presented approach. Ocean models are composed of multiple modules where the dynamical core is responsible for solving the governing equations of fluid dynamics, and modern dynamical cores use splitting methods into two main solver components: The baroclinic parts of three-dimensional fluid dynamics equations and the barotropic solver which is purely horizontal. These barotropic solvers pose particular challenges due to the numerical stiffness, requiring time steps 50 times smaller than the baroclinic one, and will be the focus of this work.

## 2.1 Barotropic solver in ocean models

This section briefly describes a barotropic solver, used, e.g., in NEMO [26] (for global ocean circulation and climate) and Croco [6] (regional ocean model for simulating coastal processes), to discuss the challenges it poses to compilers. In its easiest form, it is given by the following Partial Differential Equations (PDEs):

$$
\begin{align*}
\frac{\partial}{\partial t} D &= -\nabla \cdot (D\bar{u}) \\
\frac{\partial}{\partial t} (D\bar{u}) &= -\frac{1}{2} g h \nabla D - D f k \times \bar{u} - \nabla \cdot (D\bar{u} \bar{u}^T)
\end{align*}
$$

with $D = \zeta + h$ the water column height, $\bar{u} = (\bar{u}, \bar{v})^T$ the vertically integrated velocity, $\zeta$ the free surface elevation relative to average height, $h$ the distance from average sea level to the ocean's floor, $f$ the Coriolis effect accounting for the rotation of the earth and $g$ the gravitational acceleration. To solve this efficiently with finite differences in NEMO/Croco, on structured grids, the numerical operations are similar to stencil-like access: gathering data from one or more input arrays, followed by (non-)linear combinations of the input to compute one scalar output in the output array. For numerical reasons, Arakawa-C grid spaces [5] (aka. staggered grid) result in array iteration ranges determined by numerical properties.

To give a better impression, we
provide a simplified 1D example of
the placement of degrees of free-
dom providing the grid iteration
space with StrD/EndD and StrU/EndU referring respectively to
the start and end indices of D and u related variables where the
start/end indices might even change depending on the choice of
boundary conditions, see image. Supplemental arrays are then re-
quired to store not only the state required for the next time step (D

and $\bar{u}$ above), but also to buffer values obtained when computing
intermediate results of the equations above.

For both NEMO and Croco models, the barotropic solver is based on an explicit multi-stepping scheme [35] with further mathematical details skipped here for the sake of brevity. Its implementation is realized by extending an array with an additional dimension used in a round-robin fashion (modulo) to store the previously computed states of earlier time steps, and by overwriting results that are not further required.

## 2.2 Compiler restrictions

The motivation of our work is given by the significant loss of numerical knowledge during the discretization process of the governing equations, hindering compilers from performing optimizations.

**Observation #1 “Limited loop fusion”**: Discretization with Arakawa C-grids leads to different grid iteration ranges expressed with variables as iteration range bounds, making loop fusion more challenging for compilers (e.g., requiring multiple variants). The *numerical perspective* provides information that loops perform updates of grid data. All consecutive blocks of loops with the same number of nested loops have close iteration ranges extended/reduced by 1. *Consequently*, loop fusion will lead to only a few or no obsolete computations since there is sufficient overlap in iteration ranges.

**Observation #2 “Limited kernel fusion”**: We define a kernel as the computational body of perfectly nested loops, representing the actual operations performed within these loops. Such kernels are also called *stencil-kernels* in the literature, see, e.g., [14, 30]. Understanding data dependencies, which data is computed, and the origins of the data are highly challenging to compilers. For *explicit time integration schemes*, each output scalar represents the numerical evaluation of PDE operators (evaluation of finite difference, grid interpolation, etc.). Most importantly to us, there are no loop-carried dependencies. *Consequently*, loops can be interpreted as the application of *kernels* that can be applied in parallel without any dependencies. The outputs of PDE operators can serve as inputs to subsequent PDE operators, enabling us to leverage state-of-the-art Data Flow Graph analysis for kernel-oriented fusion of multiple consecutively executed kernels.

**Observation #3 “Identifying boundary conditions”**: Solver code on structured grids is composed of nested loops. *Numerically*, *D*-dimensional solvers represent iteration spaces over the *D* dimensional grids, and (*D*−1) nested loops represent boundary conditions being applied to different domain boundaries (left, top, bottom, or right boundary). *Consequently*, we can identify boundary conditions and even compute them in arbitrary order.

**Observation #4 “Relevance-lifetime of data”**: The relevance-lifetime of data in variables, in particular array data itself, is not always clear from the code. For instance, the data buffered to arrays might only be relevant to the current subroutine’s lifetime or time step. Compilers are required to update this data even if they are not further relevant for computations. The numerical perspective provides information about whether the data is only temporarily relevant or must persist after the current time step. Also, ocean models have various options to enable additional features in solvers, making it hard to determine the optimal array reuse for all possible

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Overview of Poseidon's stages. We make use of PSyclone in the uplifting and backend, which are part of the connector (box with dashed lines) specific to each PDE model. Internally (within the purple box), Poseidon works with a representation *independent of particular programming models or PDE solvers*.</figcaption></figure>

configurations in a programmatical way. Consequently, we can avoid
updating data in such an array and, in addition, avoid allocating
memory for it, or reuse arrays for other temporary results, reducing
the overall memory footprint.

**Observation #5 “Overall data flow”**: A lack of understanding of the overall data flow hinders the application of various compiler optimizations. A holistic compiler perspective on all numerically discretized equations is required for this. Consequently, we designed our approach as a holistic view to optimize the entire code beyond local routines.

**Observation #6 “Existing optimizations”**: Existing ocean models in Fortran are often optimized for one specific architecture. However, these existing optimizations might hinder further optimization steps. For example, vertical/horizontal loop fusion is trivial to apply manually and is often done in ocean models, but it can prevent further compiler passes. Consequently, our Poseidon development can first undo optimizations to obtain the most simplified and primitive form of the code, which is more suitable for further optimizations and automatic parallelization.

Without a doubt, some of the aforementioned problems could
be solved partly with changes in the code, but this would still not
provide a solution to all observed compiler limitations, and further
require hand-written parallelization and offloading to GPUs. Next,
we will discuss Poseidon as a potential way to overcome all of these
challenges.

3 Poseidon

Poseidon is a source-to-source compiler designed to optimize sim-
ulation models on structured grids written in Fortran. It assumes
array updates to be programmed with multiple nested loops iter-
ating over the array elements and updating them with stencil-like
access patterns from other arrays [13, 14, 20, 24, 25, 30, 33, 40]. For
ocean simulations used in production by hundreds of users, non-
standard large software packages (e.g., ROSE[32], LLVM extensions,
etc.) or software without long-life support is avoided. Poseidon is
solely based on Python with a single dependency to the open-source

Python package PSyclone, which has long-life support and is open
source. Our development can be used without requiring extensive
installation and configuration of complex software stacks, thereby
addressing the main concerns of the user community. Next, we will
discuss the different stages of Poseidon with an overview in Fig. 1.

4 Uplifting to control flow

We start by discussing the uplifting to the control flow. The Poseidon
uplifter parses and extracts the variable information from a PDE
model and characterizes the code for followup optimization steps.
We introduce here the notion of a *connector*, responsible for parsing
the relevant code, creating an intermediate representation suitable
for the next step. A connector also writes back generated code and
can modify existing code as discussed later as part of the backend.
Since each PDE model is programmed differently, each one has its
own connector *tailored* to how the PDE model is written. We will
focus on the main uplifting aspects shared among all connectors.

## 4.1 Variable uplifting & management

Poseidon uses a VVariable Management (VAM) where all used vari-
ables must be registered. During the uplifting, the connector regis-
ters all used scalars and arrays for variable management, including
their rank, dimension, and type. Also, the VAM retains (partly
implicitly) information on where they are defined, allocated, and
initialized, which becomes crucial in the backend, e.g., to remove
or allocate additional arrays in memory. The VAM distinguishes
between the following variable types where N refers to the number
of dimensions the solver was written for:

* **Scalar**: Integer or floating point scalar. E.g., for time step size, multistepping index, temporary variable within or outside a loop.

* **ArrayND**: Basic array type if the upliftler cannot characterize the array into one of the following types.

* **GridDataND(ArraND)**: Arrays of dimension *N* storing coefficients related to a particular *N*-dimensional grid space, e.g., the velocity components, momentum, or some constant metric terms.

* **MultiStepArray**: Array of dimension *N* + 1, hence, one dimension higher than GridDataND. Those arrays are required for multi-stepping schemes where one slice is given by a Multi-StepArrayElement, see following type.

* **MultiStepArrayElement GridDataND**: This variable represents a view of a slice of *MultiStepArray* indexed by an integer scalar variable. It inherits all properties from *GridDataND* and can be treated in the same way in Poseidon as a *GridDataND*, which again simplifies the followup steps.

Two crucial operations are supported for each variable in the
VAM for further processing: **Removing**: If a variable is not required
anymore, it can be flagged for *removal*. It is important to note that
the variable itself is not removed from the VAM, but the backend
must remove it from the source code. **Cloning**: In case we need
additional data storage, it can be derived from an existing variable
type of the same kind by *cloning*, potentially multiple times. Then,
a new variable of the same type will be registered to the VAM and
returned, keeping track of its origin variable. The information about
the origin variable is required in the backend to decide where to
define, allocate, and setup the *cloned* variables.
Note that these operations do not yet create any code or a code representation but are solely done in the VAM representation itself, leaving it up to the code generation in the backend to consider this, leading to a clear separation of concerns. Relevant for gaining back lost numerical information, the grid-related arrays are enriched with information about which grid they correspond to. Also, information on the variable lifetime is added, e.g., whether its data is still relevant after the time step is finished, allowing for further optimization of the memory footprint (see Obs. #4).

After the uplifting of variables, Poseidon uses these variables solely in its internal representation as a matter of separation of concerns, hiding where they are actually used in the code itself.

## 4.2 Control flow

After the registration of variables, the computations are extracted to a control flow with a macroscopic level (see, e.g., task-based control flows [38] for similar concepts) where a single node of the control flow can represent multiple operations.

### 4.2.1 Segmentation and Characterization.
For this, we first process the code using PSyclone to make it suitable for further processing, which – depending on the model – can include, e.g., rewriting array slicing syntax to loops, inlining functions to simplify interprocedural dependency analysis, eliminating dead if-branches and further features provided by PSyclone. The next task is then to further process the PSyclone IR (PSyIR) with Poseidon performing a segmentation into a coarse-grained control flow where we characterize each code block into four node types in the following order or matching:

* **“KernelLoop” nodes** are related to nested loops performing computations on grid arrays. For an N-D solver, these are composed of N perfectly nested loops (see Obs. #1). Loop bodies may only consist of scalar and array assignments, intrinsic functions, and if-else statements. KernelLoop nodes are also tagged if they are parallelizable using, e.g., PSyclone (see Obs. #2). Array iteration ranges are uplifted to an expression in PosIR discussed further below.

* **“BoundaryConditionLoop” kernel loop nodes** are related to (nested) loops applying boundary conditions on grid arrays. For an N-D model, boundary conditions are identified by (N - 1)-D nested loops. The same restrictions on loop bodies apply to kernel loop nodes (see Obs. #3).

* **“Computations” nodes** are related to computations performed outside kernel loops, frequently in conditional branches. These are, e.g., the multistep indices (see 2.1) computed relatively to the current time step, the friction coefficients of the model, some scalar variables precomputed before a kernel that uses them, etc. Consecutively occurring computations of this type are gathered in one computation node.

* **“Blackbox” nodes** are anything not characterized so far, such as I/O, calls to subroutines, and print statements. Code in these nodes is not further modified and treated as it is.

To segment and characterize a PSyIR source code representation given in Annotated Syntax Tree (AST), we use a context-free grammar given in Backus-Naur form in Fig. 2 with the starting symbol $\alpha$, the empty word $\epsilon$, and canonically given terminal symbols. The

derivation of the terminal blackbox node $B^*$ shall be defined as everything not matched by the grammar symbols $K_L$, $K_B$, and $C$.

**Starting symbol:** $\alpha$ $::=$ $K_L \alpha$ KernelLoop node
| $K_B \alpha$ BoundaryCondition node
| $C \alpha$ Computations node
| $B^* \alpha$ BlackBox node
| $\epsilon$ Empty word

**KernelLoop node:** $K_L$ $::=$ if (e) $L_N$ else
| $L_N$

**BoundaryCondition node:** $K_B$ $::=$ if (e) $L_{N-1}$ else
| $L_{N-1}$

**Perfectly nested loops:** $L_{i+1}$ $::=$ for ... do $L_i$ endfor
| $L_0$ $::=$ $C$

**ComputationBlock:** $C$ $::=$ if (e) $C$ else $C$ else
| if (e) $C$ else
| $(A \parallel C)$ +

**Assignment:** $A$ $::=$ $v = e$

**Expression:** $e$ $::=$ $v$
| $c$
| $e$ binop $e$
| unop($e$)
| intrinsic_call($e$)
| scalar variable
| scalar constant
| binary op (+, -, *, ...)
| unary op (-)
| min, abs, ...

**Scalar variable:** $v$ $::=$ $x$
| array(...) Value of scalar variable
| Value of array access

<caption>**Figure 2:** Grammar in Backus-Naur form for uplifting the source code to the control flow with starting symbol $\alpha$. See text for more details.</caption>

### 4.2.2 Variable closure.
For each node, including blackbox ones, we keep track of the variable closure, namely lists of which variables are read, written, and used within the nodes, see “input/output data” boxes in Fig. 3. All of these variables have been previously registered to the VAM (see Sec. 4.1).

Failing to determine the correct variable closure will render Poseidon either unable to perform optimizations (due to too many variables in the closure, hence obsolete dependencies) or lead to incorrect code generation (due to missing variables in the closure, hence missing dependencies).

### 4.2.3 Poseidon IR.
Computation blocks of the form $C$ by themselves or within kernel loops are further uplifted to Poseidon IR (PosIR). This representation is a list of assignments with a tree-like structure on the right-hand side, whose representation is bijective to the PSyIR concerning the subset of Fortran language supported by PosIR. We refer to compiler literature for such representations [4] and point out two main features of PosIR:

1. **Access** in the form of reading (right-hand side of an assignment) and writing (left-hand side of an assignment) from scalar and array variables is represented as leaf nodes of the same type “Access”, hiding the details of the used variable.

2. **Conditionals** can be specified for each assignment, obtained from converting (nested) if-else branches.

Both features avoid special handling of arrays and scalars and if-branches in the kernels, simplifying further analysis and rewriting of such code discussed later. After these steps, the macro-control flow of one solver step is entirely given by a list of nodes with examples in Fig. 3.
<figure><img src="image_4.png" /><figcaption>Figure 3: Two examples of the constructed control flow nodes. (top) Computation node with a source code segment uplifted to the PosIR. (bottom) Kernel loop node representing a part of the multi-stepping time integration.</figcaption></figure>

5 Data Flow Graph

With Poseidon’s separation of concerns, the control flow is now
independent of the PDE model code itself. Given the list of nodes
in the control flow, each node contains the variable closure on
the read and written variables it accesses. Based on this variable
closure, a Data Flow Graph (DFG) can be constructed (see, e.g.,
[38]) where Poseidon uses a hypergraph with an example given
in top image in Fig. 4 (including further transformations explained
later). In this graph, each node corresponds to one node of the
control flow and each hyperedge to a data flow dependency of
a variable. This includes also false dependencies on the data flow
level with respective false dependency hyperedges to account for
write-after-read and write-after-write dependencies (not displayed
in the Figure). Rather than using a regular graph, such a hypergraph
representation strongly simplifies the following DFG operations
and the realization of different backends for task-based scheduling
which is part of the following sections.

5.1 DFG operations

Next, we briefly describe the fundamental DFG operations sup-
ported by Poseidon. Rather than giving a detailed overview of all
supported operations (including validation of the DFG, debug fea-
tures, etc.), we focus on those relevant to rewriting the DFG during
the optimization process discussed in the following section.

5.1.1 *Normalization.* Before performing any HPC optimizations, we first undo existing ones and rewrite our code into what we call a *normalized form*, which forms the basis of further kernel-oriented optimizations detailed later on.

<figure><img src="image_5.png" /><figcaption>Figure 4: Exemplary DFG of the barotropic solver produced by Poseidon with data flow from top (source node) to bottom (sink node). *Ellipse-shaped nodes*: hyperedges accounting for data flow via arrays - flows of scalar variables are not rendered. *Box-shaped nodes*: data flow nodes of different types - kernel loops (yellow), computations (blue), boundary cond. (red), blackboxes (black), communication (pink - not part of this work).</figcaption></figure>

Due to the PosIR representation, with if-branches uplifted to
conditional expressions, the body of kernel loops and boundary
conditions consists only of a list of assignments. Based on this, we
further define a kernel/boundary loop with PosIR to be in **normal-**
**ized form** if and only if the following properties are fulfilled:

(a) Exactly one output variable array of dimension *N* is written to in the final assignment.

(b) No relative indices, e.g., ubar[i+1, j], are allowed in the written array in the final assignment.

(c) For all *N*-D kernel loops, the variable array written in the final assignment must not be read from in this kernel.

(d) (*N* − 1)D kernel loops can read & write variable arrays.

Property (a) relates the kernel to the evaluation of one partic-
ular finite difference/volume operator (*N*-D) or to one particular
boundary condition (($N$−1)-D), Property (c) avoids any false depen-
dencies leading again to further advantages used, e.g., for parallel
scheduling. All properties (a) - (d) result in a clearer algorithmic
implementation of DFG operations in our development.

To obtain a code in normalized form, Poseidon provides the
following DFG operations undoing common optimizations done by
hand:
<table><caption>**Table 1:** Overview of pseudo code for different fusion approaches. (Left) Vertical and horizontal fusion, often performed by hand-tuning code. (Right) Kernel-oriented fusion, requiring either synchronization (top) or duplication of operations (bottom).</caption><thead><tr><th colspan="2">Vertical Fusion* (VF)</th><th>Complex Kernel Fusion (CKF)</th></tr></thead><tbody><tr><td>for i=1...M:<br/>c[i] = a[i] * b[i+1]<br/>for i=1...M:<br/>d[i] = a[i] * c[i]<br/>⇔ result-equivalent to<br/>for i=1...M:<br/>c[i] = a[i] * b[i+1]<br/>d[i] = a[i] * c[i]</td><td>for i=1...M+1:<br/>c[i] = a[i] * b[i+1]<br/>for i=1...M:<br/>d[i] = c[i] + 2 * c[i+1]<br/>⇔ numerically-equivalent to<br/>for i=1...M+1:<br/>scr[i] = a[i] * b[i+1]<br/>--- GPU SYNC BARRIER ---<br/>if 1 &lt;= i &lt;= M:<br/>d[i] = scr[i] + 2*scr[i+1]</td></tr><tr><td>Horizontal Fusion* (HF)</td><td>Flopping Kernel Fusion* (FKF)</td></tr><tr><td>for i=1...M:<br/>c[i] = a[i] * b[i+1]<br/>for i=1...M:<br/>d[i] = a[i] * e[i]<br/>⇔ result-equivalent to<br/>for i=1...M:<br/>c[i] = a[i] * b[i+1]<br/>d[i] = a[i] * e[i]</td><td>for i=-1...M+1:<br/>c[i] = a[i] * b[i+1]<br/>for i=1...M:<br/>d[i] = c[i] + 2 * c[i+1]<br/>⇔ numerically-equivalent to<br/>for i=1...M:<br/>d[i] = (a[i] * b[i+1]) + 2 * (a[i+1] * b[i+2])</td></tr></tbody></table>

Fusion marked with * is used in the present work.

* **Single Static Assignment (SSA):** The DFG is converted to a form where each variable is assigned exactly once. Such an SSA form eliminates all false dependencies, resulting in property (c). While this leads to an increase in memory footprint, in some cases, it also allows for increased parallelism due to fewer dependencies, as discussed later.

* **Loop splitting:** Multiple array assignments can occur within a single kernel loop related to hand-tuned codes. This operation splits the kernel loop into multiple kernel loops, each with only one assignment to obtain property (a).

* **Unshifting of array outputs:** With hand-tuned loop fusion, working with loop iteration variables not intended for the 2nd merged loop requires working with relative output indices. After splitting the loop, we unshift them again to obtain property (b). If a kernel loop cannot be written in this form, Poseidon can still be used, but DFG operations operating solely on kernels in normalized form do not apply to this kernel loop. Despite using this normalized form often throughout the following optimization process, the kernel loop code that is finally produced is not necessarily in normalized form. What looks first as a de-optimization lifts us from being potentially stuck in a local extrema of an optimization process.

### 5.1.2 Fusion operations.
Over the last decades, various fusion methods have been developed [8, 15, 25, 28, 29, 31, 37, 40]. An overview of a selection of state-of-the-art fusion-based optimizations is presented in Table 1. For the sake of better understandability, we provide only simplified examples: They are based on single-dimension loops, whereas loop fusion for ocean models is applied consecutively on all perfectly nested loops. Also, scalar variables found inside each loop must be dealt with and have been skipped for the sake of brevity. Furthermore, loop iteration ranges can also mismatch for ocean models due to grid staggering, see 2.1.

We first discuss the two examples of vertical and horizontal fusion; see [25]: Vertical Fusion (VF) has an output-to-input dependency between kernel loops. Such fusion avoids reading back data from memory written before to an array. Horizontal Fusion (HF) fuses loops that share the same input dependencies to avoid reading data multiple times. Both loop fusion methods can be found in many hand-tuned codes (e.g., NEMO & Croco) and are often also supported by the state-of-the-art compilers under conditions discussed in the following section.

Complex Kernel Fusion (CKF) is a method used, e.g., in [40, 41] for GPUs, which can also be found in the Croco model with OpenMP parallelization. Applying a vertical loop fusion would lead to loop-carried dependencies, which is taken into account in multiple steps: By 1. working with subblocks, 2. extending the iteration range depending on the stencil-shapes and iteration range used in the 2nd loop, 3. writing intermediate results to GPU shared memory (scratch array local to thread group), 4. performing a barrier synchronization effectively resolving the dependency, 5. executing kernels of the 2nd loop. Such an approach requires synchronization and has idling threads, which becomes more extensive for a deeper fusion of multiple loops. Another issue arises with the requirements of direct access to GPU shared memory, specification of block size, and barrier synchronization in threads: such features are not available in all parallel programming models, e.g., OpenACC and OpenMP.

Incentivized by HPC trends of increasing gap between computational and data access performance[16, 27], we make use of a kernel fusion technique producing significantly more flops, and we will refer to it as Flopping Kernel Fusion (FKF), see right-lower code in Table 1. As a result, it leads to an increase of computational intensity and decrease of the number of required memory access aligning even better with future trends in HPC. We would like to remind of the normalized form introduced in Sec. 5.1.1, converting the kernel loops to the most suitable form for reducing the complexity of developing the FKF. We will apply this FKF independently to the iteration range of the 1st loop following its mathematical understanding (see Obs. #2). We emphasize the following points: 1. We generate code that is not necessarily valid from a compiler perspective due to the lack of updating array data. 2. Multiple vertical FKF can lead to exponential growth of code. Consequently, it is rarely done by hand since it rapidly leads to hard-to-maintain code. Additionally, if done by hand, it can introduce subtle bugs that are difficult to trace. Compared to CKF, we have advantages: (a) avoiding the synchronization barrier with multiple duplications of computations and (b) being less restrictive on resources required for the scratch array, which is not required anymore.

Concerning Poseidon, support for HF, VF, and FKF operations on the hypergraph is provided. Also, not further discussed here, all operations support fusion with assignments to temporary scalars within kernel loops.
5.1.3 DFG operation validation.

We strictly avoid trial-and-error ap-
proaches to detect invalid DFG operations
in Poseidon. Instead, all fusion operations
are validated before being applied to the
DFG with algorithms based on hyperedges
analysis (based on standard graph-analysis
algorithms).

An example of an invalid kernel fusion
DFG transformation is displayed on the
right. In the left DFG, “A” reads variable “Z” and “B” updates “Z”.
This situation is caused by reusing variables multiple times (e.g., to
reduce the memory footprint by hand).

Here, FKF of “A” and “C” would lead to an invalid DFG: After fu-
sion (depicted on the right DFG), the variable “Z” is still consumed
by the merged node “A+C”, but before this node is executed, “Z” is
updated/overwritten by node “B”. Poseidon checks for such invalid
transformations by first determining the set of all nodes on all pos-
sible paths between the two nodes to be fused. The transformation
is not allowed if one of these nodes on the path updates the variable
used as an input of the node “A” to be fused.

Since the SSA form avoids such cases by, e.g., using a different
variable for B's output, we emphasize that such situations do not
occur after normalization (see Sec. 5.1.1). Nevertheless, validation
checks are still required for such dependencies that potentially
occur due to blackbox nodes or code without normalized kernel
loops.

6 Optimization of data flow graph

Finding the optimal sequence of operations to be applied on a DFG
can be a complex process [25, 43] that depends not only on the
performance characteristics of the graph nodes but also on the
target architecture. For this work, we use a trivial, purely model-
driven optimization strategy, which, as shown later, is sufficient to
closely approach the theoretical maximum speedup.

All kernel loops have a low computational density and are con-
sequently memory bound [42], making reducing memory access a
top priority. For each kernel loop, the number of memory accesses
is given by the number of hyperedges related to array variables,
and from hereon, we will refer to them as *array hyperedges*.

Incentivised by the high computational performance of HPC
architectures, we will start by applying solely FKF. To apply FKF to
kernel loops, we make the following observations: for kernel loops
with multiple incoming array hyperedges and the output hyperedge
connected to multiple children, fusing a kernel loop to only one
other kernel loop can even lead to an actual increase in the number
of array hyperedges. Therefore, if we fuse this kernel loop, we fuse
it *to all or no children kernel loops*, consequently eliminating the
fused kernel loop. Furthermore, such a strategy also reduces the
complexity of the optimization space.

In each optimization iteration, one DFG operation can involve
multiple DFG nodes. For FKF, the optimizer chooses the kernel loop
with the maximum number of hyperedge arrays, which would be
removed by fusing it with all its children. This information can be
inferred efficiently from the information already encoded with the
hypergraph edges of all involved kernel loops. Using this optimiza-
tion criterion can lead to multiple candidates, and we use a 2nd

criterion by choosing the kernel loop with the smallest number of
arithmetic instructions. Although this is a relatively trivial opti-
mization criterion, we will observe it is sufficient to achieve results
close to the maximum achievable speedup.

7 Backend

The backend of Poseidon generates code that can be compiled within the existing PDE model. Similarly to the uplifting process, the backend is part of the *connector* and requires to be tailored to each PDE model. The following algorithms are strongly aligned with task-based execution/programming models (cf. [1, 7, 9, 18, 38]), and we will use when adequate the terminology “task” for “node”.

To determine the **order of the tasks** in the generated code, we use topological sorting [22] and refer to this as the source order. We generate code in different (parallel) variants from the task Directed Acyclic Graph (DAG) obtained from the DFG:

1. Purely sequential over the task: each task will be enqueued and executed only after the previous one in source order has finished.

2. Asynchronous execution of task (only for GPUs): each task will be directly enqueued but must wait for the previous one in source order to be finished.

3. Parallel DAG scheduling: tasks are still enqueued in source order but wait for task dependencies to be fulfilled. Boundary condition loops shown sequentially in the DFG (see red boxes in Fig. 4) are specially treated and executed in parallel alongside each other (See Obs. #3).

We support two different **parallel programming models**: OpenMP and OpenACC. For OpenMP with CPUs, dependencies of DAGs can be directly mapped to *OpenMP in/out/inout task dependencies*. In practice, we do not specify task dependencies using variables associated to DFG hyperedges, but express dependencies directly between DAG tasks with the use of dummy integer variables since OpenMP dependencies require a memory reference. Using OpenACC with Nvidia GPUs, a DAG can be expressed by submitting each task to its own queue using *async* and expressing input dependencies with *wait* on the respective queue. The scheduling of parallel DAG tasks is left to the respective parallel runtime.

The **code generation itself** is split into two parts: The first part is the generation of the Fortran code from the DFG, and the second part is related to the variable management. Code is generated for the obtained source-order and chosen parallel programming model. Special handling is required for blackbox nodes, which use an internal representation of code not fully comprehended by Poseidon, see Sec. 4. Here, we use the PSyIR to regenerate the original code as it was in the uplifting. So far, this part of the code generation is independent of the respective PDE model.

Concerning the variables, some variables might now be refer-
enced, that do not yet exist due to cloning (e.g., to convert the
DFG into the SSA, see Sec. 4.1) or variables might not be required
anymore (e.g., because they are not referenced anymore due to a
kernel fusion). In the case of cloning, this requires the declaration
of variables and also their initialization. In the case of deleted vari-
ables, their declaration, their initialization, and other references
to them must be removed. How these variables are handled can
strongly differ between PDE models, which is why the connector
must be tailored to it, and PSyclone is used for interacting with
Fortran code.

<figure><img src="image_6.png" /></figure>
# 8 Results

Benchmarks were conducted using Grid'5000: We used the “Intel Xeon Gold 5318Y” for the CPU experiments on the “spirou” partition in Louvain. GPU experiments are based on “Nvidia H100 NVL (94 GiB)” from the “musa” partition in Sophia. We also conducted benchmarks on the NVidia A100 GPU, which all looked similar to some scaling factors of the H100 results. Therefore, we only show plots for the H100, but discuss differences to the A100 if relevant.

Since the PDE solvers we are dealing with are memory-bound, we chose compilers providing the best performance on the target hardware based on a Fortran array stream benchmark, leading to the Intel LLVM-based compiler (2025.0.1.47, options -O3 -xHost) for CPUs and nvfortran (25.1-0, options -O3 -fast) with CUDA (12.2, driver 535.183.06) for GPUs, which provided close-to-peak performance on each vendor's respective architecture.

To show the performance of Poseidon we use the barotropic solver of Croco from the “soliton” test case in a research code, allowing us to focus on conducting detailed performance experiments with the same challenges during uplifting, but without the need to deal with the complexity of the full Croco development.

In the following plots, we will refer to **Hand-normalized** as one version of the code with all kernel loops given in the normalized form (see Sec. 5.1.1), hence without any hand-tuned optimization. This represents a form that can be 1:1 generated by transcribing the PDE operators into kernel loops and is, hence, also desirable by numerical experts. We will refer to **Hand-tuned** as the code version hand-tuned by HPC engineers (used in Croco), including existing fusion and OpenACC parallel loop sync annotations. For benchmarks with **Poseidon**, we use this already optimized code as a starting point for the optimization process. As a first step, we *undo all hand-made optimizations* with the DFG **Normalization** operation, see Sec. 5.1.1. Each optimization step then applies one DFG operation, which can involve multiple DFG nodes. First, we use **FKF optimizations** on kernel loops as described in Sec. 5.1.2 where each DFG operation eliminates one kernel loop. We also included results that applied a single **HF** after several FKF DFG operations - required if only horizontally fusible nodes are left.

We will refer to different problem sizes as *small* (128 × 128), *medium* (1024 × 1024), and *large* (8192 × 8192). All benchmarks run 100 time steps where the first 10 time steps have been excluded from the measurements to avoid the influence of JIT compilation (on GPUs) and other hardware-related warmup.

## 8.1 CPU benchmarks

Results for CPUs are presented in Fig. 5. We can observe Poseidon's ability to robustly parallelize and optimize the PDE model with FKF. There are certain overheads if DAG-based tasking is used over the DFG nodes. Applying one horizontal fusion after a number of FKF leads to significant performance drops. This is due to various if-branches required to cope with different iteration spaces, leaving the compiler unable to vectorize. Using task-based DAG scheduled execution, we can observe that the performance drops for the last fusion operations. This suggests this is due to concurrent scheduling of task loops within DAG tasks, leading to the competition for resources causing such overheads. For *small* problem sizes (not shown here), the task-based execution suffers from severe tasking overheads, leading to a performance drop of over ×10. An overview of speedups after applying all DFG operations is given in here:

<figure><img src="image_7.png" /><figcaption>**Figure 5:** Benchmarks conducted on CPU with medium problem size with min/max plotted as transparent filled area. We compare the *hand-normalized*, *hand-optimized*, and *Poseidon* with FKF, HF optimizations, and different parallel DAG tasking strategies.</figcaption></figure>

<table><thead><tr><td>Problem size</td><td>Hand-normalized</td><td>Hand-optimized</td></tr></thead><tbody><tr><td>Large</td><td>×1.87</td><td>×1.78</td></tr><tr><td>Medium</td><td>×2.04</td><td>×1.72</td></tr><tr><td>Small</td><td>×1.65</td><td>×1.25</td></tr></tbody></table>

## 8.2 GPU benchmarks

We continue with in-depth performance experiments on the GPU as the most relevant target platform for NEMO & Croco models. For this, we investigate three different offloading strategies of kernels to GPUs; all realized in Poseidon's backend (see Sec. 7): (a) ACC parallel loop only, (b) ACC parallel loop with *async* submission to a single queue and (c) ACC parallel loop with *async dispatching of kernel loops to different queues* for DAG tasking. We continue to point out the automatic parallel scheduling of boundary conditions to different queues to execute them in parallel.

We first investigate the results of a large problem size depicted in Fig. 6 with each array exceeding the L2 cache size. This provides the optimal case for fusion strategies - avoiding memory access - for which we optimize. We can observe the normalized version of the code to be significantly slower than the hand-optimized one due to the additional memory access being introduced due to the normalization. Starting from this normalized form, we apply FKF DFG operations by choosing kernel loops to be fused next as described in Sec. 6. We observe a robust performance optimization for each of the first 16 fusion operations with some plateauing effects afterward. Nevertheless, by applying FKF, we continue to obtain more performance improvements. Compared to the *hand-optimized* version, we can achieve significant speedups of about 2 (see table below for details). Note that there are four kernel loops left that cannot be fused by FKF any further.

Next, we investigate the performance of the FKF optimization on the *medium* and *small* problem sizes with results shown in Fig. 7. For the *medium* sized problem, we observe that parallel graph-based execution of the DFG nodes provides the best performance for many DFG nodes. The more FKF we apply, the more its runtime overheads, leading to a performance drop of over ×10. An overview of speedups after applying all DFG operations is given in here:
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Benchmarks conducted on GPU with *large* problem size and with min/max plotted as transparent filled area. For better comparison, baselines for hand-normalized / -tuned timings are drawn as a horizontal line, although no DFG operations are applied. (Left image) Benchmark comparisons of *hand-normalized*, *hand-tuned*, and *Poseidon* with DFG operations are given in parentheses. If specified, *HF* is applied only once after a certain number of FKF, counting as one more DFG operation. The fully FKF-optimized version, followed by one more HF, provides the best results. (Right image) Bandwidth information and array utilization for the *large* problem size. The modeled bandwidth is based on the number of memory accesses and the size of the arrays.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 7: Benchmarks conducted on GPU with *medium* (left) and *small* (right) problem size with min/max plotted as transparent filled area. For better comparison, baselines for hand-normalized / -tuned timings are drawn as a horizontal line, although no DFG operations are applied. For the *medium* problem size, we observe that DAG async execution provides the best performance for a low number of applied DFG operations. This is different for the *small* problem size, where we observe overheads of parallel execution (and synchronization). In both cases, the fully FKF-optimized version followed by one more HF provides the best results.</figcaption></figure>

approaches the async dispatching of nodes to a single queue since fewer and fewer nodes are available for parallel execution.

For the *small* problem size, we observe that the parallel dispatching no longer provides the best performance if no FKF is applied, but the async dispatching to a single queue performs best. We account for this relative slowdown by parallel dispatching and synchronization overheads.

Although this paper is more focused on the methodological side of source-to-source code transformation, we want to point out that the FKF followed by a single HF DFG optimization is not only a theoretical concept but also provides significant performance improvements in practice with speedups given here:

<table><thead><tr><td>Problem size</td><td>Hand-normalized</td><td>Hand-optimized</td></tr></thead><tbody><tr><td>Large</td><td>×2.10</td><td>×1.80</td></tr><tr><td>Medium</td><td>×2.19</td><td>×1.76</td></tr><tr><td>Small</td><td>×2.60</td><td>×2.16</td></tr></tbody></table>

## 8.3 Bandwidth, memory footprint, theoretical speedup and impact on math

Since the fusion optimizations are motivated mainly by the bandwidth limitation, we also investigate how much bandwidth each time step consumes for the *large* problem size. We model the utilized bandwidth as follows: For each kernel loop, the number of accessed arrays is given by the number of respective DFG hyperedges. We multiply this with the size of each array (inferred from the resolution) and the number of time steps, providing the total number of bytes accessed in each kernel loop. The total bandwidth is then given by summing the bytes read/written over all kernel loops and dividing by the time for one timestep. We plotted this upper bound in the right top image of Fig. 6 (top black line). We observe that, despite introducing more computations due to FKF, we stay close to the optimal bandwidth of the GPU for the first 16 FKF operations. Afterward, the memory bandwidth drops, which we account for by
the increase in register usage, resulting in a loss of occupancy and
therefore parallel requests to the memory. Nevertheless, despite
this drop from the maximum bandwidth, FKF still leads to further
speedups.

With GPUs becoming increasingly available in regular workstations and laptops, another desire is to run simulations of higher resolutions directly on a single GPU. Due to the VAM, Poseidon can free variables (scalsars and arrays) when they are no longer needed due to the FKF operation. We can observe in Fig. 6 a significant reduction of the required temporary arrays: Each FKF operation reduces the number of temporary arrays by one, and applying FKF to all nodes avoids allocating 22 temporary arrays. With $4 \times 4 = 16$ arrays for the multi-stepping arrays, and one temporary array left (output of upper kernel loop in Fig. 4), Poseidon can reduce the memory footprint by $(16+1+22)/(16+1) = 2.29$. Hence, Poseidon can also be used to optimize for reducing the memory footprint.

Assuming the bandwidth is the only limiting factor for an arbitrary number of FKF, we can also compute the theoretical maximum speedup that can be obtained after FKF: The DFG with normalized nodes used here as an example has 102 input/output arrays in total. After all FKFs are applied, this is reduced to 52 arrays, leading to a potential speedup of 102/52 ≈ 1.96. Comparing this to the FKF benchmark plots, we conclude to be *close to the optimal reachable speedup* compared to the Hand-normalized baseline.

After obtaining this deeper understanding of the bandwidth
limitation, we can also conclude that the multi-stepping method
itself can be identified as the main performance limitation after
FKF, since it requires reading in 3 times more arrays than a single-
stepping method would require, making single-stepping methods a
prime candidate for future bandwidth-limited HPC architectures.

9 Related Work

To our knowledge, it is the first time that a numerics-HPC co-
design has been applied to automatically transform existing non-
parallelized Fortran code into highly performing variants for GPUs
and CPUs and reusing the same Fortran code for the backend. In
addition, we address one of the operational requirements, which
is to depend only on a small software stack that can be installed
on any computer. Next, we will show relations of some subparts of
Poseidon to other developments where we will skip well-known
annotation-based approaches such as OpenMP and OpenACC for
brevity.

Optimization for stencil computations has been a longer-standing
research topic, and many approaches (see, e.g., [14, 20, 21, 30, 31,
33, 36, 40, 41, 43]) have been proposed to improve the performance
of stencil codes on modern architectures. Our work exploits the
increasing compute performance vs. bandwidth ratio of modern
computer architectures with FKF for the first time for solvers of such
a complexity, with results close to the max. available bandwidth.

The standard way to obtain performance portability is the utilization of **Domain Specific Languages** (e.g., ExaStencils [23], OPS (Oxford Parallel library for Structured mesh solvers) [34], xDSL [10], or GridTools [3], Stella [19]). However, these approaches require a complete rewrite of the code, which is not always feasible for existing PDE models. Instead, our approach allows developers to keep working within their programming language and existing

code but still requires using a specific coding style co-developed
with the Poseidon connector.

There have also been alternatives to DSL, e.g., based on the
ROSE source-to-source compiler [32], but it requires a larger soft-
ware stack of particular versions to be compiled and installed. ROSE
was used to perform a CUDA-to-CUDA compilation [41] with CKF,
which requires code already written in CUDA and using CUDA
as the backend. ClawDSL [12] allows annotating generic code and
transforming it to a parallel version but without further optimiza-
tions. PSyclone is both: a source-to-source and DSL compiler for
codes/DSL written in Fortran. On the one hand, it is used in the
LFRic [2] weather and climate model, using a DSL approach for
kernel invocation. While it can also be applied directly to source
codes such as the NEMO and Croco ocean models, it only provides
straightforward annotation-based optimizations that could be ap-
plied by hand. As presented in this work, our approach builds on
PSyclone, providing complementary features to allow for a holistic
optimization approach and performance modeling based on a DFG.

Other recent work [25] for GPUs required a substantial search
time of at least half an hour, whereas Poseidon achieves significant
performance improvements close to the maximum possible speedup
on GPUs within less than 8 seconds for the best-performing variant
(full FKF + 1 HF). Nevertheless, such an hybrid model-driven &
autotuning approach will become necessary once more DFG nodes
are dealt with and, even if only minor further speedups could be
gained, using in addition their approach would be highly relevant
for operational executions.

10 Discussion on Poseidon’s limitations

We demonstrated a new approach to cope with automatic parallel-
elization and optimization of existing Fortran code with a highly
efficient model-driven optimizer that applies deep kernel fusion
within a few seconds.

However, we do not want to raise the expectation that this will
be applicable out-of-the-box to all Fortran code. Since each PDE
model is written differently, Poseidon requires the adaptation of the
connector for each PDE model. Further, Poseidon can only optimize
what can be expressed in the VAM and the DFG. Also, optimizations
can only be performed if kernel loops are transformed e.g. into a
normalized form.

Some ocean models (e.g., NEMO) pose some particular challenges: E.g., modules using private array variables pose some limitations in the uplifting and backend parts in the connector to make them accessible with a rewritten DFG in the backend. We do not see this as a limitation of the core of Poseidon itself, but in the uplifting and backend that requires more extensive work in the future. An easy remedy would be to adapt existing PDE models to better fit the connectors requirements as part of a software co-engineering process.

11 Summary

This paper introduces Poseidon, a novel framework designed to
transform legacy Fortran code into a high-level representation
based on DFGs. By leveraging Poseidon IR, the framework supports
advanced DFG operations, including kernel fusion and performance-
oriented transformations, to optimize computations for modern
HPC architectures. The backend retains Fortran compatibility, en-
suring support for vendor-specific compilers and maximizing hard-
ware performance. Experimental results demonstrate significant
performance gains, with deep kernel fusion and loop fusion achiev-
ing up to 2.16× speedup compared to hand-tuned code and a compile-
time of less than 8 seconds. Future work aims to extend Poseidon
with features like automatic communication injection and efficient
automatic differentiation.

Acknowledgements

This work is based on funding from the IRGA project, from ADUM
and the Inria associate team "Crocodiles". We thank Peter Messmer
and Lukas Mosimann for their exceptional support with the H100.
We used Grammarly and Co-Pilot for improvements in language
and style. Experiments presented in this paper were carried out
using the Grid'5000 testbed, supported by a scientific interest group
hosted by Inria and including CNRS, RENATER and several Univ. as
well as other organizations (see https://www.grid5000.fr).
References

[1] 2018. OpenMP Application Programming Interface 5.0. November (2018). https://www.openmp.org/wp-content/uploads/OpenMP-API-Specification-5.0.pdf

[2] S.V. Adams, R.W. Ford, M. Hambley, J.M. Hobson, I. Kavcic, C.M. Maynard, T. Melvin, E.H. Mueller, S. Mullerworth, A.R. Porter, M. Rezny, B.J. Shipway, and R. Wong. 2019. LFRic: Meeting the challenges of scalability and performance portability in Weather and Climate models. *J. of Parallel and Distributed Computing* 132 (2019), 383–396.

[3] Anton Afanasyev, Mauro Bianco, Lukas Mosimann, Carlos Osuna, Florian Thaler, Hannes Vogt, Oliver Fuhrer, Joost VandeVondele, and Thomas C. Schulthess. 2021. GridTools: A framework for portable weather and climate applications. *Software X* 15 (2021), 100707. doi:10.1016/j.softx.2021.100707 Publisher: Elsevier B.V.

[4] V Aho Alfred, S Lam Monica, and D Ullman Jeffrey. 2007. *Compilers principles, techniques & tools*. pearson Education.

[5] Akio Arakawa and Vivian R. Lamb. 1981. A Potential Enstrophy and Energy Conserving Scheme for the Shallow Water Equations. doi:10.1175/1520-0493(1981)109<0018:APEAEC>2.0.CO;2 ISBN: 0027-0644 ISSN: 0027-0644 Issue: 1 Pages: 18–36 Publication Title: Monthly Weather Review Volume: 109.

[6] Francis Auclair, Rachid Benshila, Lucie Bordois, Martial Boutet, Maurice Brémond, Matthieu Caillaud, Gildas Cambon, Xavier Capet, Laurent Debreu, Nicolas Ducoussou, François Dufois, Franck Dumas, Christian Éthé, Jonathan Gula, Christophe Hourdin, Serena Illig, Swen Jullien, Matthieu Le Corre, Solène Le Gac, Sylvie Le Gentil, Florian Lemarié, Patrick Marchesiello, Camille Mazoyer, Guillaume Morvan, Cyril Nguyen, Pierrick Penven, Renaud Person, Joris Pianezze, Stéphane Pous, Lionel Renault, Laurent Roblon, Andres Sepulveda, and Sebastien Theetlen. 2024. *Coastal and Regional Ocean Community model*. doi:10.5281/zenodo.11036115

[7] Cédric Augonnet, Samuel Thibault, Raymond Namyst, and Pierre-André Wacrenier. 2009. StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures. (2009).

[8] David F. Bacon, Susan L. Graham, and Oliver J. Sharp. 1994. Compiler transformations for high-performance computing. *Comput. Surveys* 26, 4 (Dec. 1994), 345–420. doi:10.1145/197405.197406

[9] Michael Bauer, Sean Treichler, Elliott Slaughter, and Alex Aiken. 2012. Legion: Expressing locality and independence with logical regions. In *SC '12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis*. 1–11. doi:10.1109/SC.2012.71 ISSN: 2167-4337.

[10] Nick Brown, Maurice Jamieson, Anton Lydike, Julien Bauer, and Tobias Grosser. 2023. Fortran performance optimisation and auto-parallelisation by leveraging MLIR-based domain specific abstractions in Flang. In *Proceedings of the SC '23 Workshops of The International Conference on High Performance Computing, Network, Storage, and Analysis (SC-W '23)*. Association for Computing Machinery, New York, NY, USA, 904–913. doi:10.1145/3624062.3624167

[11] Katherine Calvin, Dipak Dasgupta, Gerhard Krinner, Aditi Mukherji, Peter W. Thorne, Christopher Trisos, José Romero, Paulina Aldunce, Ko Barrett, Gabriel Blanco, William W.L. Cheung, Sarah Connors, Fatima Denton, Aida Diongue-Niang, David Dodman, Matthias Garschagen, Oliver Geden, Bronwyn Hayward, Christopher Jones, Frank Jotzo, Thelma Krug, Rodel Lasco,UNE-Yi Lee, Valérie Masson-Delmotte, Malte Meinshausen, Katja Mintenbeck, Abdalah Mokssit, Friederike E.L. Otto, Minal Pathak, Anna Pirani, Elvira Polozcanska, Hans-Otto Pörtner, Aromar Revi, Debra C. Roberts, Joyashree Roy, Alex C. Ruane, Jim Skea, Priyadarshi R. Shukla, Raphael Slade, Aimée Slangen, Youba Sokona, Anna A. Sörensson, Melinda Tignor, Detlef Van Vuuren, Yi-Ming Wei, Harald Winkler, Pan-mao Zhai, Zinta Zommers, Jean-Charles Hourcade, Francis X. Johnson, Shonali Pachauri, Nicholas P. Simpson, Chandni Singh, Adelle Thomas, Edmond Totin, Paola Arias, Mercedes Bustamante, Ismail Elgizouli, Gregory Flato, Mark Howden, Carlos Méndez-Vallejo, Joy Jacqueline Pereira, Ramón Pichs-Madruga, Steven K. Rose, Yamina Saheb, Roberto Sánchez Rodríguez, Diana Urge Vorsatz, Cunde Xiao, Noureddine Yassaa, Andrés Alegria, Kyle Armour, Birgit Bednar-Friedl, Konrélis Blok, Guéldadio Cissé, Frank Dentener, Siri Eriksen, Erich Fischer, Gregory Garner, Céline Guivarch, Marjolijn Haasnoot, Gerrit Hansen, Mathias Hauser, Ed Hawkins, Tim Hermans, Robert Kopp, Noëmie Leprince-Ringuet, Jared Lewis, Debora Ley, Chloé Ludden, Leila Niamir, Zebedee Nicholls, Shreya Some, Sophie Szopa, Blair Trewin, Kaj-Ivar Van Der Wijst, Gundula Winter, Maximilian Witting, Arlene Birt, Meeyoung Ha, José Romero, Jinmi Kim, Erik F. Haites, Yonghun Jung, Robert Stavins, Arlene Birt, Meeyoung Ha, Dan Jezreel A. Orendain, Lance Ignon, Semin Park, Youngin Park, Andy Reisinger, Diego Cammaramo, Andreas Fischlin, Jan S. Fuglestvedt, Gerrit Hansen, Chloé Ludden, Valérie Masson-Delmotte, J.B. Robin Matthews, Katja Mintenbeck, Anna Pirani, Elvira Polozcanska, Noëmie Leprince-Ringuet, and Clotilde Péan. 2023. *IPCC, 2023: Climate Change 2023: Synthesis Report. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change*. IPCC, Geneva, Switzerland. Technical Report. Intergovernmental Panel on Climate Change (IPCC). doi:10.59327/IPCC/AR6-9789291691647 Edition: First.

[12] Valentin Clement, Sylvaine Ferrachat, Oliver Fuhrer, Xavier Papillonne, Carlos E. Osuna, Robert Pincus, Jon Rood, and William Sawyer. 2018. The CLAW DSL: Abstractions for performance portable weather and climate models. *Proceedings of the Platform for Advanced Scientific Computing Conference, PASC 2018* (2018). doi:10.1145/3218176.3218226 ISBN: 9781450358910.

[13] Kaushik Datta, Shoait Kamill, Samuel Williams, Leonid Oliker, John Shalf, and Katherine Yelick. 2009. Optimization and performance modeling of stencil computations on modern microprocessors. *SIAM Rev.* 51, 1 (2009), 129–159. doi:10.1137/070693199

[14] K. Datta, M. Murphy, V. Volkov, S. Williams, J. Carter, L. Oliker, D. Patterson, J. Shalf, and K. Yelick. 2008. Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures. In *2008 SC - International Conference for High Performance Computing, Networking, Storage and Analysis*. IEEE, Austin, TX, 1–12. doi:10.1109/SC.2008.5222004

[15] Chen Ding and Ken Kennedy. 2000. The Memory Bandwidth Bottleneck and its Amelioration by a Compiler. (2000).

[16] Jack Dongarra, John Gunnels, Harun Bayraktar, Azzam Haidar, and Dan Ernst. 2024. Hardware Trends Impacting Floating-Point Computations In Scientific Applications. doi:10.48550/arXiv.2411.12090 arXiv:2411.12090 [math].

[17] ECMWF. 2021. ECMWF Strategy 2021-2030. (2021). doi:10.21957/S21EC694KD Publisher: ECMWF.

[18] Alejandro Fernández, Vicenç Beltran, Xavier Martorell, Rosa M. Badia, Eduard Œguadé, and Jesús Labarta. 2014. Task-Based Programming with OmpSs and Its Application. In *Euro-Par 2014: Parallel Processing Workshops*, Luis Lopes, Julius Žilinskas, Alexandru Costan, Roberto G. Cascella, Gabor Kecskemeti, Emmanuel Jeanlot, Mario Cannataro, Laura Ricci, Siegfried Benkner, Salvador Petit, Vittorio Scarano, José Gracia, Sascha Hunold, Stephen L. Scott, Stefan Lankes, Christian Lengauer, Jesús Carretero, Jens BREITBART, and Michael Alexander (Eds.). Springer International Publishing, Cham, 601–612. doi:10.1007/978-3-319-14313-2 51

[19] Tobias Gysi, Christoph Müller, Oleksandr Zinenko, Stephan Herhut, Eddie Davis, Tobias Wicky, Oliver Fuhrer, Torsten Hoefler, and Tobias Grosser. 2021. Domain-Specific Multi-Level IR Rewriting for GPU: The Open Earth Compiler for GPU-accelerated Climate Simulation. *ACM Trans. Archit. Code Optim.* 18, 4 (Sept. 2021), 51:1–51:23. doi:10.1145/3469030

[20] Bastian Hagedorn, Larisa Stoltzfus, Michel Steuwer, Sergei Gorlatch, and Christophe Dubach. 2018. High performance stencil code generation with Lift. In *Proceedings of the 2018 International Symposium on Code Generation and Optimization*. ACM, Vienna Austria, 100–112. doi:10.1145/3168824

[21] Georg Hager and Gerhard Wellein. 2011. *Introduction to High Performance Computing for Scientists and Engineers*.

[22] A B Kahn. 1962. Topological sorting of large networks. (1962).

[23] Christian Lengauer, Sebastian Kuckuk, Ulrich Rüde, Armin Größlinger, Matthias Bolten, Harald Köstler, Hannah Rittich, Jonas Schmitt, Christian Schmitt, Alexander Grebhahn, Sven Apel, Lisa Claus, Stefan Kronawitter, Jürgen Teich, Frank Hannig, Shigeru Chiba, and Stefan Groth. 2020. *ExaStencils: Advanced Multigrid Solver Generation*. Technical Report FZJ-2021-00052. Springer. https://juser.fz-juelich.de/record/889127 ISBN: 9783030479558.

[24] Randall J. LeVeque. 2007. *Finite Difference Methods for Ordinary and Partial Differential Equations*. doi:10.1137/1.9780898718739 ISSN: 00465070_eprint: arXiv:1011.1669v3.

[25] Xiaoyan Liu, Xinyu Yang, Kejie Ma, Shanghao Liu, Kaige Zhang, Hailong Yang, Yi Liu, Zhongzhi Luan, and Depei Qian. 2024. Moirae: Generating High-Performance Composite Stencil Programs with Global Optimizations. In *Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis (SC '24)*. IEEE Press, Atlanta, GA, USA, 1–15. doi:10.1109/SC41406. 2024.00026

[26] Gurvan Madec, Mike Bell, Rachid Benshila, Adam Blaker, Romain Boudrallé-Badie, Clément Bricaud, Diego Bruciaferri, Davi Carneiro, Miguel Castrillo, Daley Calvert, Jérôme me Chanut, Emanuela Clementi, Andrew Coward, Casimir de Lavergne, Srdan Dobricic, Italo Epicoco, Christian Éthé, Emma Fiedler, David Ford, Rachel Furner, Jonas Ganderton, Tim Graham, James Harle, Katherine Hutchinson, Doroteaciro Iovino, Robert King, Dan Lea, Claire Levy, Tomas Lovato, Eric Maisonnave, Julian Mak, Juan Manuel Castillo Sanchez, Matt Martin, Nicolas Martin, Diana Martins, Sébastien Masson, Pierre Mathiot, Francesca Mele, Silvia Mocavero, Aimie Moulin, Simon Müller, George Nurser, Paolo Oddo, Stella Paronuzzi, Julien Paul, Mathieu Peltier, Renaud Person, Clement Rousset, Stefanie Rynders, Guillaume Samson, David Schroeder, Dave Storkey, Andrea Storto, Sibylle Téchené, Martin Vancoppenolle, and Chris Wilson. 2024. NEMO Ocean Engine Reference Manual. doi:10.5281/zenodo.14515373 Version Number: v5.0.

[27] John David McCalpin. 2016. Memory bandwidth and system balance in hpc systems. (2016).

[28] K.S. McKinley. 1998. A compiler optimization algorithm for shared-memory multiprocessors. *IEEE Transactions on Parallel and Distributed Systems* 9, 8 (Aug. 1998), 769–787. doi:10.1109/71.706049

[29] Kathryn S. McKinley, Steve Carr, and Chau-Wen Tseng. 1996. Improving data locality with loop transformations. *ACM Transactions on Programming Languages and Systems* 18, 4 (July 1996), 424–453. doi:10.1145/233561.233564
[30] Anthony Nguyen, Nadathur Satish, Jatin Chhugani, Changkyu Kim, and Pradeep Dubey. 2010. 3.5-D Blocking Optimization for Stencil Computations on Modern CPUs and GPUs. In *2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis*. IEEE, New Orleans, LA, USA, 1-13. doi:10.1109/SC.2010.2

[31] Bo Qiao, Oliver Reiche, Frank Hannig, and Jirgen Teich. 2019. From Loop Fusion to Kernel Fusion: A Domain-Specific Approach to Locality Optimization. In *2019 IEEE/ACM International Symposium on Code Generation and Optimization (CGO)*. IEEE, Washington, DC, USA, 242-253. doi:10.1109/CGO.2019.8661176

[32] Dan Quinlan and Chunhua Liao. 2011. The ROSE Source-to-Source Compiler Infrastructure. (2011).

[33] Prashant Singh Rawat, Miheer Vaidya, Aravind Sukumaran-Rajam, Atanas Rountev, Louis-Noël Pouchet, and P Sadayappan. 2019. On Optimizing Complex Stencils on GPUs. (2019).

[34] István Z. Reguly, Gihan R. Mudalige, and Michael B. Giles. 2018. Loop Tiling in Large-Scale Stencil Codes at Run-Time with OPS. *IEEE Transactions on Parallel and Distributed Systems* 29, 4 (April 2018), 873-886. doi:10.1109/TPDS.2017.2778161 Conference Name: IEEE Transactions on Parallel and Distributed Systems.

[35] Alexander F. Shchepetkin and James C. McWilliams. 2009. Computational Kernel Algorithms for Fine-Scale, Multiprocess, Longtime Oceanic Simulations. In *Handbook of Numerical Analysis*. Vol. 14. Elsevier, 121-183. doi:10.1016/S1570-8659(08)01202-0

[36] Prashant Singh, Aravind Sukumaran-Rajam, Atanas Rountev, Fabrice Rastello, Louis-Noël Pouchet, and Ponnuswamy Sadayappan. 2018. Register Optimizations for Stencils on GPUs. (2018).

[37] Sharad Singhai and Kathryn McKinley. 1996. Loop Fusion for Data Locality and Parallelism. (1996).

[38] Oliver Sinnen. 2007. *Task Scheduling for Parallel Systems* (1 ed.). Wiley. doi:10.1002/0470121173

[39] Sergi Siso, Andrew R. Porter, and Rupert W. Ford. 2023. Transforming Fortran weather and climate applications to OpenCL using PSyclone. In *Proceedings of the 2023 International Workshop on OpenCL (IWOCL '23)*. Association for Computing Machinery, New York, NY, USA, 1-8. doi:10.1145/3585341.3585360

[40] Mohamed Wahib and Naoya Maruyama. 2014. Scalable Kernel Fusion for Memory-Bound GPU Applications. In *SC14: International Conference for High Performance Computing, Networking, Storage and Analysis*. IEEE, New Orleans, LA, USA, 191-202. doi:10.1109/SC.2014.21

[41] Mohamed Wahib and Naoya Maruyama. 2015. Automated GPU Kernel Transformations in Large-Scale Production Stencil Applications. In *Proceedings of the 24th International Symposium on High-Performance Parallel and Distributed Computing*. ACM, Portland Oregon USA, 259-270. doi:10.1145/2749246.2749255

[42] Samuel Williams, Andrew Waterman, David Patterson, Soda Hall, and U C Berkeley. 2009. Roofline: An Insightful Visual Performance Model for Floating-Point Programs and Multicore Architectures. (2009).

[43] Yongpeng Zhang and Frank Mueller. 2013. Autogeneration and Autotuning of 3D Stencil Codes on Homogeneous and Heterogeneous GPU Clusters. *IEEE Transactions on Parallel and Distributed Systems* 24, 3 (March 2013), 417-427. doi:10.1109/TPDS.2012.160