# Portage GPU d'un solveur éléments finis discontinus hybridisé pour les problèmes d'ondes en fréquence

Ahmed Chabib, Roland Greffe, Christophe Geuzaine, Axel Modave

► To cite this version:

Ahmed Chabib, Roland Greffe, Christophe Geuzaine, Axel Modave. Portage GPU d'un solveur éléments finis discontinus hybridisé pour les problèmes d'ondes en fréquence. CFM 2025 - 26e Congrès Français de Mécanique, Aug 2025, Metz, France. hal-05235205

HAL Id: hal-05235205

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

Submitted on 2 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# Portage GPU d'un solveur éléments finis discontinus hybridisé pour les problèmes d'ondes en fréquence

A. Chabibª, R. Greffeᵇ, C. Geuzaineᵇ et A. Modaveª

a. POEMS, CNRS, Inria, ENSTA, Institut Polytechnique de Paris, 91120 Palaiseau, France

b. Institut Montefiore, Université de Liège, Belgique

## Résumé :

Dans ce travail, nous nous intéressons à la résolution par éléments finis de problèmes de propagation d'ondes en régime harmonique de très grande taille. L'utilisation de cartes graphiques (GPU) permet d'accélérer les calculs, mais il est difficile d'en exploiter pleinement la puissance. Nous considérons une méthode d'éléments finis discontinus de type Galerkin (DG) avec des flux amonts, hybridisée utilisant des variables de transmission définies aux faces des éléments. L'élimination des variables physiques conduit à un système linéaire adapté pour une résolution itérative et une implémentation parallèle efficace sur GPU. Après une description de la méthode, appelée CHDG, nous présentons quelques stratégies de mise en œuvre sur GPU et nous comparons et discutons leurs performances.

## Abstract :

In this work, we deal with the finite element solution of large-scale time-harmonic wave propagation problems. Using graphics processing units (GPUs) can speed up the computations, but it is difficult to take full advantage of their computing power. We consider a discontinuous Galerkin (DG) finite element method with upwind fluxes hybridized with transmission variables defined at the element faces. Eliminating the physical variables yields a linear system whose structure lends itself to efficient parallel iterative solution on GPUs. After describing the method, called CHDG, we present some implementation strategies on GPUs, and we compare and discuss their performance.

## Mots clefs : Éléments finis, HDG, Ondes, Helmholtz, GPU, HPC

## 1 Introduction

La simulation numérique de phénomènes de propagation des ondes est largement utilisée dans l'industrie et le monde académique pour répondre à d'importants défis, comme la réduction du bruit des avions ou l'évaluation des risques sismiques. Lorsqu'un nombre limité de fréquences doit être étudié, des approches basées sur la résolution d'équations en régime harmonique, dont l'équation de Helmholtz est l'exemple typique, sont généralement envisagées.

Les approches basées sur des méthodes d'éléments finis ont démontré leur capacité à résoudre de tels problèmes, notamment grâce à leur robustesse et leur flexibilité pour traiter des configurations géométriques et physiques complexes. Dans ce travail, nous considérons des méthodes d'éléments finis
discontinus de type Galerkin (DG), qui permettent de représenter des solutions discontinues et offre plus de flexibilité que les méthodes continues standards (utilisation de maillage non-conformes, fonctions de base non-standards, etc). Dans le cadre des problèmes en régime instationnaire, ces méthodes se prettent bien mieux au calcul sur GPU, que les méthodes continues. Par exemple, dans [5, 7], une approche basée sur des fonctions bases nodales de Lagrange facilite l'implémentation sur GPU, et permet d'exploiter ces architecture massivement parallèle. Des stratégies sont utilisées pour paralléliser les calculs correspondants à différents éléments, et pour gérer les échanges de données entre les éléments voisins.

L'hybridisation consiste à introduire des variables additionnelles aux interfaces entre les éléments, et à reformuler le problème global en termes de ces variables. Le sytème obtenu après élimination des variables physiques possède alors moins de degrés de liberté. Avec les méthodes hybridisées classiques, les variables de transmission correspondent à des traces numériques, voir [2]. La méthode CHDG [8] est une alternative pour les problèmes de propagation d'ondes en régime harmonique, qui repose sur l'utilisation de variables de transmission aux interfaces. Par rapport aux approches classiques, le nombre d'inconnues est doublé, mais la résolution itérative avec des solveurs classiques est plus rapide en termes de nombre d'itérations. C'est dans ce contexte que s'inscrit ce travail, qui vise à implémenter cette méthode sur GPU, en tirant parti de la nature intrinsèquement parallèle des méthodes DG hybridisées.

Après une présentation du modèle et des méthodes, nous détaillerons différentes stratégies d'implémentation sur GPU, et nous les étudions sur un cas de référence 3D.

## 2 Modèle mathématique, méthode numérique et schéma discret

Pour présenter notre approche, nous considérons le problème acoustique en régime harmonique suivant :

$$ egin{cases} -\iota \kappa p + \nabla \cdot \mathbf{u} = 0, & \text{dans } \Omega, \\ -\iota \kappa \mathbf{u} + \nabla p = 0, & \text{dans } \Omega, \\ p - \mathbf{n} \cdot \mathbf{u} = s_R, & \text{sur } \partial \Omega, \end{cases} \quad (1) $$

où $p$ and $\mathbf{u}$ sont les champs de pression et de vitesse respectivement, $\kappa$ est le nombre d'onde, $n$ est le vecteur unitaire normal à la surface, dirigé vers l'extérieur, et $s_R$ représente une source surfacique imposée sur $\partial\Omega$, bord du domaine $\Omega$.

### 2.1 Méthode des éléments finis de Galerkin discontinu

Nous considérons un maillage du domaine $\Omega$, composé de tétraèdres, dont les éléments et les faces sont notés $K$ et $F$, respectivement. Sur chaque élément, les champs physiques sont représentés par des polynômes de degré $n$. Pour un élément $K$, l'espace d'interpolation est noté $\mathbb{P}_n(K)$, et les champs approchés sont notés $p_K$ et $u_K$. Nous introduisons les produits scalaires

$$ (x, y)_K = \int_K x \bar{y} \, d\mathbf{x} \quad \text{et} \quad \langle x, y \rangle_F = \int_F x_{|F} \bar{y}_{|F} \, d\mathbf{x}. $$

Les mêmes notations sont utilisées pour les produits scalaires de champs vectoriels.

Avec la méthode DG (Discontinuous Galerkin) considérée, la solution numérique vérifie le problème local suivant sur chaque élément $K$ :
**Problème 1.** Trouver $p_K \in \mathbb{P}_n(K)$ et $\mathbf{u}_K \in [\mathbb{P}_n(K)]^3$, tel que :

$$ egin{cases} -\iota\kappa(p_K, q_K)_K - (\mathbf{u}_K, \nabla q_K)_K + \sum_F \langle \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F, q_K \rangle_F = 0, \\ -\iota\kappa(\mathbf{u}_K, \mathbf{v}_K)_K - (p_K, \nabla \cdot \mathbf{v}_K)_K + \sum_F \langle \hat{p}_F, \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F = 0, \end{cases} \quad (2) $$

$\forall (q_K, \mathbf{v}_K) \in \mathbb{P}_n(K) \times [\mathbb{P}_n(K)]^3$, où $\mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F$ et $\hat{p}_F$ représentent les flux numériques et la somme $\sum_F$ est appliquée sur toutes les faces $F$ de l'élément $K$.

Les flux numériques imposent la continuité de la solution aux interfaces entre les éléments, ainsi que les conditions aux limites au bord du domaine. Dans le cadre de ce travail, nous utilisons le flux amont (ou “upwind flux”) défini sur chaque face $F$ par

$$ \left\{ \begin{array}{ll} \hat{p}_F = (p_K + p_{K'})/2 + \mathbf{n}_{K,F} \cdot (\mathbf{u}_K - \mathbf{u}_{K'})/2 & \text{si } F \not\subset \partial\Omega, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F = \mathbf{n}_{K,F} \cdot (\mathbf{u}_K + \mathbf{u}_{K'})/2 + (p_K - p_{K'})/2 & \end{array} \right. \quad (3) $$

$$ \left\{ \begin{array}{ll} \hat{p}_F = (p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K + s_R)/2 & \text{si } F \subset \partial\Omega, \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F = (p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K - s_R)/2 & \end{array} \right. \quad (4) $$

où $\mathbf{n}_{K,F}$ le vecteur unitaire normal extérieur à $K$ sur la face $F$, et $K'$ désigne un élément voisin de $K$ si $F \not\subset \partial\Omega$; voir [4] pour plus de détails.

## 2.2 Hybridation avec le variable de transmission entrante

L'hybridation consiste à introduire une variable supplémentaire définie sur les faces des éléments, puis à éliminer les variables physiques afin d'obtenir un système défini uniquement sur les faces. Ce processus modifie les propriétés du système sans altérer la solution physique récupérable par post-traitement. Avec la méthode CHDG, la variable supplémentaire correspond à une variable de transmission, voir [8].

Sur chaque face $F$ d'un élément $K$, nous introduisons la variable de transmission sortante $g_{K,F}^{\oplus}$ et la variable de transmission entrante $g_{K,F}^{\ominus}$ définies par

$$ g_{K,F}^{\oplus} := p_K + \mathbf{n}_K \cdot \mathbf{u}_K \quad \text{et} \quad g_{K,F}^{\ominus} := \begin{cases} g_{K',F}^{\oplus} & \text{si } F \not\subset \partial\Omega, \\ s_R & \text{si } F \subset \partial\Omega. \end{cases} $$

Ces variables sont interprétées comme des données échangées à travers les interfaces entre les éléments, ainsi qu'avec les frontières du domaine. Pour une face donnée, la variable sortante d'une maille correspond à la variable entrante de la maille voisine.

Le problème local peut être réécrit en utilisant la variable de transmission entrante comme donnée :

**Problème 2.** Trouver $p_K \in \mathbb{P}_n(K)$ et $\mathbf{u}_K \in [\mathbb{P}_n(K)]^3$, tel que :

$$ \left\{ \begin{array}{ll} -\iota\kappa(p_K, q_K)_K - (\mathbf{u}_K, \nabla q_K)_K + \sum_F \frac{1}{2} \langle p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K, q_K \rangle_F = \sum_F \frac{1}{2} \langle g_{K,F}^{\ominus}, q_K \rangle_F, & (5) \\ -\iota\kappa(\mathbf{u}_K, \mathbf{v}_K)_K - (p_K, \nabla \cdot \mathbf{v}_K)_K + \sum_F \frac{1}{2} \langle p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K, \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F = -\sum_F \frac{1}{2} \langle g_{K,F}^{\ominus}, \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F, & \end{array} \right. $$

$\forall (q_K, \mathbf{v}_K) \in \mathbb{P}_n(K) \times [\mathbb{P}_n(K)]^3$.
En introduisant la variable sortante comme une variable inconnue sur toutes les faces, nous pouvons
découpler les problèmes physiques. L’élimination des variables physiques, et l’utilisation des relations
entre les variables sortantes, conduit au problème global qui s’écrit de façon abstraite comme suit :

**Problème 3.** Trouver $g_h^\ominus \in \bigcup_K \bigcup_F \mathbb{P}_n(F)$, tel que

$$
\sum_K \sum_F \langle (I - \Pi S) g_h^\ominus, \xi_h \rangle_F = \sum_K \sum_F \langle b, \xi_h \rangle_F, \qquad (6)
$$

$$
\forall \xi_h \in \bigcup_K \bigcup_F \mathbb{P}_n(F).
$$

Dans ce problème, nous avons introduit :

* l'opérateur d'échange $\Pi : \bigcup_K \bigcup_F \mathbb{P}_n(F) \to \bigcup_K \bigcup_F \mathbb{P}_n(F)$, qui permet de communiquer les données entre les éléments (en exploitant le fait que la variable sortante d'une maille correspond à la variable entrante de la maille voisine) et d'imposer les conditions aux limites :

$$
\Pi(g^{\oplus})|_{K,F} = \begin{cases} g_{K',F}^{\oplus} & \text{si } F \not\subset \partial\Omega \text{ et partagée entre } K \text{ et } K', \\ 0 & \text{si } F \subset \partial\Omega; \end{cases} \tag{7}
$$

* l'opérateur de diffraction $S : \bigcup_K \bigcup_F \mathbb{P}_n(F) \to \bigcup_K \bigcup_F \mathbb{P}_n(F)$, permet le passage d'une variable entrante vers une variable sortante :

$$
S(g_h^{\ominus})|_{K,F} := p_K(g_h^{\ominus}) + n_{K,F} \cdot u_K(g_h^{\ominus}), \quad (8)
$$

où les variables $(p_K, u_K)$ sont obtenues en résolvant le problème local 2 sur tous les éléments,
en utilisant la variable de transmission entrante $g_h^\ominus$ comme donnée;

* le membre de droite global $b$, qui dépend de la source $s_R$ du problème :

$$
b|_{K,F} = \begin{cases} 0 & \text{si } F \not\subset \partial\Omega, \\ s_R & \text{si } F \subset \partial\Omega; \end{cases} \tag{9}
$$

* le vecteur inconnu $g_h^\ominus$, avec $g_h^\ominus|_{K,F} := g_{K,F}^\ominus$, qui regroupe toutes les variables sortantes définies sur le squelette du maillage.

2.3 Schéma discret

Dans notre implémentation, les champs numériques sont représentés avec de fonctions de base nodales
de Lagrange, en suivant ce qui est utilisé pour des implémentations DG en temps [4, 5]. Le système
linéaire obtenu prend alors la forme suivante :

$$
M(I - \Pi S)g = Mb \tag{10}
$$

où *M* est la matrice de masse, les matrices **Π** et **S** représentent les versions discrètes des opérateurs **Π**
et **S**, et le vecteur **b** correspond au second membre **b**. Le vecteur **g** contient les valeurs nodales de $g_h^\ominus$
de toutes les faces du maillage. Pour plus de détails sur la structure du système linéaire, voir [6]. Vu les
propriétés des opérateurs, ce système peut être résolu avec l'algorithme du point fixe :

$$g^{\ell+1} = \Pi_S g^\ell + b, \quad \ell = 0, 1, \dots \qquad (11)$$

où $g^0$ est une solution de départ et $\ell$ est le pas d'itération. En pratique, d'autres méthodes itératives plus performantes devraient être utilisées (méthodes de Krylov, préconditionnement, etc). Dans ce travail, focalisé sur le portage sur GPU de la méthode CHDG, nous nous limiterons au point fixe.

## 3 Stratégies d'implémentation sur GPU

L'application de la matrice $\Pi_S$ constitue l'étape crutiale de la procédure de résolution. L'application de la matrice d'échange $\Pi$ ne demande que des transferts de données. En revanche, appliquer $S$ demande la résolution des systèmes locaux sur l'ensemble des éléments du maillage, ce qui est l'étape la plus coûteuse, et doit faire l'objet d'efforts d'optimisation.

Deux stratégies se posent :

1. La première consiste à construire et résoudre les systèmes locaux à la volée "On the Fly" (OTF) à chaque application de $S$. La résolution sera faite à l'aide de solveurs directs de la librairie MAGMA [1].

2. La deuxième consiste à inverser les matrices des systèmes locaux de tous les éléments, et à construire explicitement les blocs diagonaux de la matrice $S$ dans une étape d'initialisation. Cette matrice sera conservée dès sa création jusqu'à la fin de l'exécution. Son application se simplifie alors à un produit matrice-vecteur. Nous appelons cette stratégie "Input To Output" (I2O).

Dans le modèle de programmation CUDA, les calculs sont exécutés par un grand nombre d'unités de calcul qui s'exécutent simultanément, appelées threads. Ces dernières sont regroupées en blocs (thread blocks), unités d'exécution qui partagent une mémoire rapide et peuvent se synchroniser entre elles. Plussieurs blocs de threads sont ensuite organisés en une grille (grid). Chaque thread possède un identifiant unique qui lui permet de gérer ses propres données, ce qui permet d'exploiter massivement la puissance des GPU, en maximisant le parallélisme, en minimisant les accès mémoires coûteux et en réduisant également les copies entre le CPU et le GPU.

Nous avons adopté une stratégie qui consiste à attribuer un nœud du maillage à chaque "thread" de la grille et un élément à chaque "block". Toutes les données associées au maillage et nécessaires pour calculer les matrices éléments finis sont copiées sur le GPU avant le début des calculs. Cela permet d'exécuter la quasi-totalité des opérations directement sur le GPU, tout en minimisant les transferts de données vers le CPU.

<table><caption>**Algorithm 1: Implémentation OTF**</caption><tr><td>1</td><td>Application de la matrice S</td><td></td></tr><tr><td>2</td><td>Input : Variable de transmission sortante g<sup>&ominus;</sup></td><td></td></tr><tr><td>3</td><td>Output : Variable sortante g<sup>&oplus;</sup> = Sg<sup>&ominus;</sup></td><td></td></tr><tr><td>4</td><td>Assemblage des systèmes locaux</td><td>(Noyaux Volume et Surface)</td></tr><tr><td>5</td><td>Résolution des systèmes locaux</td><td>(Batched MAGMA)</td></tr><tr><td>6</td><td>Calcul de g<sup>&oplus;</sup><sub>K,F</sub> &larr; p<sub>K</sub> + n<sub>K,F</sub> &sdot; u<sub>K</sub> pour chaque F et K</td><td>(Noyau Transmission)</td></tr></table>
<table><caption>**Algorithm 2: Implémentation I2O**</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Précalcul</td><td></td></tr><tr><td>2</td><td>Assemblage des systèmes locaux</td><td>(Noyaux <i>Volume</i> et <i>Surface</i>)</td></tr><tr><td>3</td><td>Inversion des matrices locales</td><td>(Batched MAGMA)</td></tr><tr><td>4</td><td>Création de la matrice <i>S</i></td><td>(Deux noyaux)</td></tr><tr><td>5</td><td>Application de la matrice <i>S</i></td><td></td></tr><tr><td>6</td><td><i>Input</i> : <i>Variable de transmission sortante</i> <i>g</i><sup>⊕</sup></td><td></td></tr><tr><td>7</td><td><i>Output</i> : <i>Variable sortante</i> <i>g</i><sup>⊕</sup> = <i>Sg</i><sup>⊕</sup></td><td></td></tr><tr><td>8</td><td><i>g</i><sup>⊕</sup> ← <i>Sg</i><sup>⊕</sup></td><td>(Noyau <i>MatVec</i>)</td></tr></tbody></table>

La construction des systèmes locaux est réalisée par deux noyaux CUDA spécialisés. Un noyau "Volume" s'occupe des termes volumiques du problème 2, tandis qu'un noyau de "Surface" est dédié aux termes surfaciques de la matrice et du second membre. Le passage des variables physiques vers la variable sortante $g^{\oplus}$ est implémenté à l'aide d'un noyau *Transmission* peu coûteux. Ces noyaux sont appelés au cours des itérations dans le cadre de l'approche OTF, tandis qu'ils ne sont exécutés qu'une seule fois lors de l'initialisation dans le cas de la stratégie I2O. Pour cette dernière, les lignes associées au second membre sont également éliminées et l'application de la matrice de diffraction nécessite ensuite un appel à un noyau *MatVec* chargé d'exécuter le produit matrice-vecteur. Les pseudocodes 1 et 2 permettent d'avoir une idée de la manière dont ces noyaux sont utilisés dans les deux approches.

Pour la résolution des problèmes locaux, généralement de petite taille et en grand nombre, nous avons adopté une approche basée sur l'utilisation des opérations par lot "batch operation" préimplémentées dans MAGMA. Cette stratégie permet de traiter simultanément plusieurs petits calculs indépendants en regroupant leurs exécutions. En exploitant les capacités de parallélisme massif des GPU, cette méthode optimise l'utilisation des ressources matérielles tout en réduisant les surcoûts. Nous avons utilisé dans notre code C++ les fonctions "LU batched" de cette bibliothèque pour le calcul de l'inverse des matrices dans le cas I2O, et pour la résolution des systèmes dans le cas d'OTF.

# 4 Étude de performance

## 4.1 Configuration et problème de référence

Les expérimentations ont été réalisées sur un noeud de la machine LUCIA¹, disposant d'un GPU NVIDIA A100 doté de 40 Go de mémoire et d'un processeur AMD EPYC 7513 (32 Cores). Le GPU utilisé possède une performance maximale de 9.7 TFLOP/s (opérations en virgule flottante par seconde) et une bande passante maximale de 1.5 TB/s (octets par seconde).

Pour les études, nous considérons un cas de référence défini sur le cube $\Omega = ]0, 1[ \times ]0, 1[ \times ]0, 1[$ avec une condition limite de Robin sur la frontière $\partial\Omega$. Nous utilisons la solution $p(x) = e^{\iota\kappa d \cdot x}$, avec $d = \sqrt{1/3} [1, 1, 1]^T$, qui représente la propagation d'une onde plane dans le cube unité. Cette solution est utilisée dans le second membre de la condition aux limites. La partie réelle de le solution est illustrée dans la figure 1. Le logiciel Gmsh [3] est utilisé pour la génération de maillages et la visualisation.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Partie réelle de l’onde plane $κ = 2π$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – (a) Le débit arithmétique en fonction du degré d’interpolation polynomial pour les différents noyaux, (b) Analyse Roofline des noyaux CUDA. Les chiffres sur la figure indiquent le degré polynomial $n$ à chaque nœud. Les valeurs intermédiaires ne sont pas présentées : seuls $n = 1$ et $n = 7$ sont indiqués.</figcaption></figure>

## 4.2 Évaluation de la performance des noyaux

D'une manière générale, la performance des noyaux est limitée par le nombre d'opérations à virgule flottantes (FLOP) effectuées (*limitation due à la capacité de calcul du GPU*) et/ou par le volume de données à transférer aux unités de calcul (*limitation due à la bande passante maximale*). Deux grandeurs sont donc particulièrement importantes : le **débit arithmétique** (c'est-à-dire le nombre d'opérations *par seconde*) et l'utilisation de la **bande passante** (c'est-à-dire la quantité de données transférées en *octet par seconde*). Pour identifier la cause principale de limitation de performance, nous analysons également l'**intensité opérationnelle**, définie comme le rapport entre le nombre d'opérations flottantes et la quantité de données transférées. L'étude de performance de type Roofline permet de mettre en relation ces quantiées, voir [11].

Nous étudions la peformance des noyaux appelés au cours des itérations pour les deux stratégies d'implémentation : les noyaux *Volume*, *Surface* et de la résolution *LU* pour OTF, et le noyau de produit *MatVec*

<footnote n="1">1. https://www.oci-hpc.be/clusters.html</footnote>
pour I2O. Sur la figure 2a, le débit arithmétique est présenté en fonction du degré polynomial (de $n = 1$ à $n = 7$) pour un nombre d'éléments fixe (8106). Sur la figure 2b, nous présentons l'analyse Roofline : le débit arithmétique des différent noyaux est représenté en fonction de l'intensité opérationnelle. Les noyaux à gauche de la ligne verticale pointillée sont limités par les transferts mémoire (*memory-bound*), alors qu'ils sont limités par le calcul (*compute-bound*) à droite de cette ligne. Pour chaque zone, la ligne noire épaisse correspond au plafond de performance théorique du matériel. Cette démarche permet d'obtenir une compréhension plus fine des facteurs limitant la performance des noyaux.

Pour la stratégie I2O, nous constatons que le noyau principal *MatVec* opère dans une zone *memory-bound*. Il atteint environ 57 – 80% de la performance maximale. Le débit arithmétique s'améliore quand nous augmentons le degré polynomial, ce qui correspond à augmenter la taille des matrices. Toutefois, cet approche, qui demande de stocker les blocs diagonaux de la matrice *S*, et limité par l'espace mémoire, ce qui empêche d'augmenter *n* au-delà de 5 pour le maillage considéré.

Pour la stratégie OTF, les noyaux *Volume* et *Surface*, qui sont également *memory-bound*, sont relativement éloignés du plafond, n'exploitant respectivement qu'environ 5 – 30% et 30 – 40% de la puissance maximale du GPU. Les performances du noyau *Surface* s'améliorent légèrement avec l'augmentation du degré polynomial *n*, alors que celles du *Volume* se dégradent. Le noyau *LU* est *memory-bound* pour $n \le 3$, et ensuite *compute-bound* lorsque nous augmentons *n*. Ses performances atteignent 10 – 30%, et ensuite 45%, des capacités maximales. Le changement du comportement dans la région *compute-bound* s'explique alors par le coût de calcul qui devient suffisamment gros pour compenser les limitations dues aux transferts mémoire.

## 4.3 Comparaison des stratégies d'implémentation

Pour comparer les stratégies d'implémentations, nous mesurons le temps de calcul par itération de l'al goritme du point fixe et par élément. Ces temps sont présentés en fonction du nombre d'éléments sur la figure 3 pour différents degrés polynômiaux (de $n = 1$ à $n = 7$).

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Temps moyen par itération du point fixe et par élément pour les deux implémentations (OTF et I2O). L'échelle de l'axe des ordonnées diffère d'une figure à l'autre.</figcaption></figure>
Nous pouvons faire plusieurs observations :

* L'implémentation I2O est toujours plus rapide que OTF en terme de temps de calcul par itération. Pour un petit nombre d'éléments, le rapport entre des temps des deux implémentations est d'environ 30 pour des degrés faibles, et peut atteindre environ 100 pour $n = 7$. Ce rapport devient encore plus important lorsque nous augmentons le nombre d'éléments.

* Pour les deux implémentations et $n = 1$, le temps par élément décroit linéairement lorsque nous augmentons le nombre d'éléments, et ce temps diminue plus lentement lorsque nous augmentons $n$. Pour OTF, le coût ne diminue plus et atteint un plateau pour un nombre d'éléments suffisamment grand, contrairement à l'approche I2O dont les performances continuent à s'améliorer.

Les résultats présentés portent uniquement sur le coût de calcul des noyaux principaux utilisés au cours des itérations. Le coût du calcul de l'étape d'initialisation et l'empreinte mémoire, qui sont des aspects importants à prendre en compte, surtout pour l'approche I2O, feront l'objet d'une étude ultérieure.

# 5 Conclusion

Dans ce travail, nous nous intéressons au portage sur GPU d'un solveur éléments finis pour la résolution de problèmes de propagations d'ondes en régime harmonique. Le solveur repose sur une méthode d'éléments finis discontinus hybride qui introduit de variables de transmission aux faces des éléments, et qui élimine les varibles physiques. Le système linéaire résultant possède de meilleures propriétés pour des résolutions que les méthodes standards, mais le portage sur GPU reste une étape délicate.

Nous avons présenté et étudié deux stratégies d'implémentation pour du calcul sur un seul GPU : l'approche OTF (On-The-Fly) résout le système à chaque appel, et l'approche I2O (Input To Output) consiste à créer la matrice du système hybridisé dans une étape de calcul préliminaire, puis à effectuer des produits matrice-vecteur dans le solveur du système réduit.

En utilisant un problème de référence, nous avons constaté que l'approche I2O était en moyenne 30 fois plus rapide que l'approche OTF pour les degrés polynomiaux faibles, et qu'il pouvait avoir un facteur 100 pour les degrés élevés. Cette différence s'explique par la complexité différente des deux approches ainsi que par les performances des noyaux appelés au cours des itérations pour les deux méthodes. Bien que la stratégie I2O soit beaucoup plus rapide qu'OTF, elle nécessite une étape de pré-calcul et davantage de mémoire pour le stockage de la matrice de diffraction, ce qui peut devenir un facteur limitant.

La méthode CHDG est étendue à plusieurs physiques, telles que l'électromagnétisme [10] et les milieux hétérogènes [9], ainsi qu'à l'utilisation de différents solveurs. Le portage sur GPU est alors envisagé pour ces types de problèmes, de même que l'exploitation de configurations multi-GPU pour traiter des cas complexes de très grande taille.

# Remerciements

Ce travail a été soutenu par le projet ANR JCJC WavesDG (ANR-21-CE46-0010). Dans le cadre de France 2030, ce travail a bénéficié d'une aide de l'Etat gérée par l'Agence Nationale de la Recherche et attribuée au projet Exa-MA (qui fait partie du programme prioritaire de recherche NumPEx) sous la référence ANR-22-EXNU-0002. Les moyens de calcul ont été fournis par le Consortium des Équipements de Calcul Intensif (CÉCI), financé par le Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) sous la convention n° 2.5020.11 et par la Région Wallone. Les présents travaux ont bénéficié de moyens
de calcul mis à disposition sur Lucia, le supercalculateur Tier-1 de la Région Wallonne, infrastructure financée par la Région wallonne sous la convention n°1910247.

Références

[1] A. Abdelfattah, H. Azzam, T. Stanimire & J. Dongarra, Batched one-sided factorizations of tiny matrices using GPUs : Challenges and countermeasures, *J. Comput. Sci.* 26, 226-236 (2018).

[2] B. Cockburn, J. Gopalakrishnan & R. Lazarov, Unified Hybridization of Discontinuous Galerkin, Mixed, and Continuous Galerkin Methods for Second Order Elliptic Problems, *SIAM J. Numer. Anal.* 47, 1319-1365 (2009)

[3] C. Geuzaine & J.-F. Remacle, Gmsh : A 3-D finite element mesh generator with built-in pre- and post-processing facilities, *Int. J. Numer. Methods Eng.* 79, 1309-1331 (2009).

[4] J. S. Hesthaven & T. Warburton, *Nodal discontinuous Galerkin methods : algorithms, analysis, and applications*, Springer Science & Business Media (2007).

[5] A. Klöckner, T. Warburton, J. Bridge & J. S. Hesthaven, Nodal discontinuous Galerkin methods on graphics processors, *J. Comput. Phys.* 228, 7863–7882 (2009)

[6] A. Modave, Accelerated iterative DG finite element solvers for large-scale time-harmonic acoustic problems. *INTER-NOISE and NOISE-CON Cong. and Conf. Proc.*, Nantes, France, 1301-1307

[7] A. Modave, A. St-Cyr & T. Warburton, GPU performance analysis of a nodal discontinuous Galerkin method for acoustic and elastic models, *Computers & Geosciences* 91, 64-76 (2016).

[8] A. Modave & T. Chaumont-Frelet, A hybridizable discontinuous Galerkin method with characteristic variables for Helmholtz problems, *J. Comput. Phys.*, 493 (2023) 112459.

[9] S. Pescuma, G. Gabard, T. Chaumont-Frelet & A. Modave, A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic acoustic problems in heterogeneous media, *J. Comput. Phys.* 534 (2025)

[10] A. E. Rappaport, T. Chaumont-Frelet & A. Modave, A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic electromagnetic problems - Preprint(2025)

[11] S. Williams, A. Waterman & D. Patterson, Roofline : an insightful visual performance model for multicore architectures, *Communications of the ACM* 52, 65–76(2009).