# Portage sur GPU du calcul des forces internes d'un code basé sur la méthode des éléments finis spectraux

Tom Budon, Florent de Martin, Sylvain Jubertie, Sébastien Limet, Emmanuel Melin

► To cite this version:

Tom Budon, Florent de Martin, Sylvain Jubertie, Sébastien Limet, Emmanuel Melin. Portage sur GPU du calcul des forces internes d'un code basé sur la méthode des éléments finis spectraux. Conférence francophone d'informatique en Parallélisme, Architecture et Système (COMPAS 2025), Jun 2025, Bordeaux, France. hal-05127166

HAL Id: hal-05127166

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

Submitted on 24 Jun 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 sur GPU du calcul des forces internes d'un code basé sur la méthode des éléments finis spectraux

Tom Budon (1), Florent De Martin (2), Sylvain Jubertie (3), Sébastien Limet (1), Emmanuel Melin (1)

Université d'Orléans, INSA CVL, LIFO, UR 4022, Orléans, France. (1)
BRGM (Service Géologique National), 45060 Orléans, France. (2)
Intel Corporation, Paris, France. (3)
tom.budon@univ-orleans.fr

## Résumé

La méthode des éléments finis spectraux est utilisée dans de nombreux codes de simulations. Ces codes sont très coûteux en terme de calculs et nécessitent donc d'être parallélisés. Dans le cadre de la simulation de la propagation d'ondes sismiques, nous proposons de faciliter la parallélisation sur GPU du cœur de cette méthode en s'appuyant sur une approche existante vectorisée sur CPU. L'attention est portée sur l'utilisation de la mémoire partagée et les éventuels conflits d'accès. La méthode proposée pour le portage sur GPU permet d'obtenir un speedup de 7.7 par rapport à la version vectorisée sur CPU utilisant un cœur.

**Mots-clés** : Parallélisme, GPU, Éléments finis spectraux

## 1. Introduction

La simulation numérique est importante dans de nombreux domaines scientifiques pour la compréhension et la prédiction de phénomènes. En géosciences l'étude de phénomènes complexes comme les tremblements de terre nécessitent des simulations coûteuses en terme de calculs et de mémoire. Par exemple le code *EFISPEC3D*<fnref n="1" /> [1], utilisé par le BRGM (Bureau de Recherches Géologiques et Minières), est conçu pour calculer des sismogrammes dans des environnements géologiques 3D complexes afin de prédire l'aléa sismique local. Afin de supporter cette charge de calcul, le code s'exécute sur des supercalculateurs composés de nombreux noeuds. Toutefois ce type de configuration est coûteux, à la fois financièrement et en ressources énergétiques. Le code *EFISPEC3D* est un code développé pour exploiter des CPU et un portage vers des GPU (Graphics Processing Units) permettrait de réduire ces coûts, d'améliorer les performances et d'accroître la taille ou la résolution des problèmes traités. La figure 1 présente une visualisation à un instant donné des résultats d'une simulation d'*EFISPEC3D*.

<footnote n="1">1. https://gitlab.com/brgm/seismic-modelling/efispec3d</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Arrêt sur image de la prédiction des mouvements du sol après la rupture d'une hypothétique faille.</figcaption></figure>

## 2. Travaux antérieurs

En analyse numérique, la méthode des éléments finis (MEF), ou Finite Element Method (FEM), est employée pour résoudre numériquement des équations aux dérivées partielles (EDP). Ces équations peuvent servir à décrire le comportement de nombreux systèmes physiques. La méthode des éléments finis vise la résolution approchée d'un problème. Une solution approchée de l'équation du mouvement est trouvée grâce à une discrétisation [6] du domaine qui réduit l'espace des solutions.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – (Gauche) Représentation du maillage d'un milieu (Droite) Représentation des points GLL d'un élément</figcaption></figure>

La méthode des éléments finis spectraux est une formulation de la méthode des éléments finis utilisant des polynômes à degrés élevés. Cette méthode a l'avantage d'obtenir des résultats plus précis sans nécessiter un maillage plus fin. En plaçant les points d'intégrations aux points Gauss-Lobatto-Legendre (GLL), la méthode gagne en efficacité au travers d'une matrice de masse diagonale par construction. Augmenter l'ordre de la solution recherchée augmente également le nombre de points GLL, ce qui permet d'améliorer la qualité du résultat mais nécessite plus de temps de calcul. Pour cette étude et sans perte de généralité, nous nous concentrons sur des bases polynomiales d'ordre 4 utilisées dans *EFISPEC3D*, permettant un compromis entre la charge de calcul et la précision du résultat. Pour cette configuration, il y a 5 points GLL par dimension, et par conséquent 125 points GLL par élément en 3D (Figure 2 schéma de droite).
Il a été montré que les compileurs ne sont pas capables de vectoriser efficacement le noyau d'**EFISPEC3D** [4]. Ce travail de vectorisation d'**EFISPEC3D** a été réalisé manuellement en utilisant les *intrinsics* SIMD (Single Instruction Multiple Data). L'approche consiste à traiter plusieurs éléments simultanément suivant la largeur des registres SIMD (de 4 à 16). Cette approche pour la vectorisation a été choisie car elle permet de conserver la structure originelle du code. Les versions vectorisées du code **EFISPEC3D** utilisant AVX2 et AVX-512 permettent d'obtenir un speedup de 2 par rapport à une version vectorisée par le compilateur [10].

En ce qui concerne le GPU, plusieurs publications traitent de portage de code utilisant la méthode des éléments finis. Ces publications portent sur la modification du stockage des données en mémoire pour favoriser les accès coalescents [9]. Certaines proposent des algorithmes de coloriage et de partitionnement permettant d'éviter des synchronisations ([2], [3], [7], [8], [9]). Différents formats de stockage des matrices de rigidité sont aussi proposés [9]. Certaines approches évitent le stockage des matrices en les recalculant à chaque itération [5]. Toutefois, il est peu fait mention de la disposition des données dans la mémoire partagée GPU ainsi que des potentiels conflits liés à cette mémoire, et la majorité des publications portent sur la méthode des éléments finis et non la méthode des éléments finis spectraux.

# 3. Portage GPU du calcul des forces internes

## 3.1. Répartition du calcul

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – À gauche, organisation des *threads* en *blocks*. À droite, organisation des cœurs dans un SM</figcaption></figure>

Les GPU sont conçus pour exécuter des milliers de *threads* en parallèle. Pour tirer partie de la puissance de calcul des GPU, nous utilisons CUDA (language GPU de NVIDIA). Comme montré sur la figure 3 (schéma de gauche), en CUDA le parallélisme est divisé en trois parties, les *threads* qui sont l'unité de base d'exécution, les *blocks* qui sont des groupes de *threads* qui peuvent partager de la mémoire et la *grid* qui est un ensemble de *blocks*. L'intérêt de ce découpage est qu'il est possible d'établir une correspondance entre les *blocks* et les unités de calcul du GPU, appelées SM (Streaming Multi-processor), décrits à droite sur la figure 3. Pour exécuter un code GPU, on appelle une fonction codée en CUDA appelée *kernel*. Ce *kernel* est lancé avec un nombre défini de *threads* et de *blocks*. Le modèle d'exécution de CUDA repose sur la création de nombreux *threads* puis leurs exécution par décomposition en *warps*. Les *warps* sont des groupes de 32 *threads* synchrones qui sont exécutés sur les unités de calcul SIMD (les SM). Le modèle ne garantit pas la synchronisation entre les *warps*.

Les GPU disposent de différents types de mémoire en leur sein, chacune a des caractéristiques
spécifiques, la mémoire globale est accessible par tous les *threads* mais a une latence plus élevée, la mémoire partagée (*shared memory*) est partagée entre les *threads* d'un même *block*, elle a une latence plus faible que la mémoire globale.

## 3.2. Organisation mémoire

Chaque groupe de données, appartenant à un seul élément est stocké en *shared memory*. Il est uniquement utilisé par le *thread* qui travaille sur cet élément. En conséquence un *block* contenant un seul *warp* (groupe de 32 *threads*) contient les données relatives à 32 éléments. Si le nombre de *threads* par *block* augmente, alors les données stockées en *shared memory* augmentent. Ce paramètre est à prendre en compte dans le choix du stockage des données en *shared memory* à cause des limites de stockage de celle-ci.

Le tableau contenant les déplacements pour chaque point GLL utile au calcul des forces internes, qui doit être stocké en *shared memory*, contient les trois valeurs de déplacements des trois dimensions pour chaque point GLL. Un élément possédant 125 points GLL, 375 valeurs sont stockées par élément. Comme montré sur la figure 4, les valeurs sont rangées successivement par groupe de 375 dans la *shared memory*. Le *thread* s'occupant du premier élément accède aux 375 premières valeurs, le second *thread* accède à partir de la 376ème valeur et ainsi de suite. La *shared memory* est organisée en bancs mémoire, ces bancs sont tous de la même taille et il y en a 32 par SM, soit le même nombre que de *threads* dans un *warp*. Les données stockées dans la *shared memory* sont rangées successivement dans les bancs par groupes de 4 octets successifs, c'est-à-dire des flottants simple précision dans le cas d'**EFISPEC3D**. La figure 4 représente les bancs mémoire de la *shared memory* : le schéma de gauche montre le stockage par défaut des données d'un élément au travers de plusieurs bancs, et le schéma de droite montre le stockage réorganisé des données d'un élément par banc mémoire.

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Représentation des 2 stratégies de rangement sur les bancs de la *shared memory*</figcaption></figure>

Avec l'agencement mémoire par défaut de la *shared memory*, des problèmes de conflits d'accès aux bancs mémoire peuvent survenir. Si plusieurs *threads* accédaient à différentes adresses situées dans le même banc mémoire, les accès seraient séquentialisés, cela entraînerait une baisse de la bande passante de la *shared memory*, ce qui pourrait avoir des conséquences sur le temps d'exécution. A contrario, la bande passante est optimale si tous les *threads* accèdent à des bancs différents ou si tous les *threads* accèdent à la même adresse sur le même banc. Dans le cadre de la méthode des éléments finis spectraux à l'ordre 4 en 3D, chaque *thread* accède à un banc mémoire différent. La raison en est ici que les valeurs 375 et 32 (nombre de bancs mémoire) sont premières entre elles. Cette manière de stocker le tableau permet ici des accès optimaux, ce qui ne serait pas forcement le cas à d'autres ordres.
## 4. Expérimentation

Les tests de performance présentés dans cette section ont été réalisés sur AMD Ryzen Threadripper PRO 5945WX 12-Cores pour le CPU et sur NVIDIA RTX A4000 pour le GPU. Les compilateurs utilisés sont GFORTRAN 9.4.0 pour le code Fortran, GCC 14.1.0 pour le code C/C++ et NVCC 10.1 pour le code CUDA.

Pour la première stratégie (placement par défaut), l'efficacité de notre placement mémoire en *shared memory* est confirmée lors des mesures réalisées avec Nsight-compute, l'outil de profilage CUDA de Nvidia. La métrique *smsp_perNJwavefront*<fnref n="2" /> permet de mettre en évidence les problèmes de conflits d'accès aux bancs mémoire de la *shared memory*. Cette métrique contient 3 valeurs, la première valeur définit le nombre d'octets maximum pouvant être transférés par la *shared memory*. La valeur obtenue est de 128 octets, en effet si chacun des 32 *threads* accède à une valeur de 4 octets, on obtient 128 octets transférés depuis la *shared memory* (4 octets × 32 bancs = 128 octets). La deuxième valeur représente le pourcentage des 128 octets réellement envoyé en moyenne par transfert. Et la troisième valeur donne le nombre moyen d'octets envoyés par transfert. Dans le cas de cette version l'utilisation maximale des bancs est atteinte, 100% des 128 octets sont utilisés à chaque transfert mémoire. Si nous sortons du cadre de l'ordre 4 en 3D, la deuxième stratégie consiste à ranger les valeurs par banc mémoire (à droite sur la figure 4). Cette version plus générale permet d'obtenir les mêmes résultats que ceux de la première version au prix d'une réorganisation mémoire. L'ordre 4 étant figé pour *EFISPEC3D*, dans la suite nous utilisons la première version non réorganisée.

### 4.1. Extraction du kernel de calcul

Lors de précédents travaux, le noyau principal de calcul d'*EFISPEC3D* a été isolé au sein d'une Miniapp, notamment pour simplifier son étude et étudier sa vectorisation. Pour la suite, nous reprenons cette Miniapp et en proposons une version GPU.

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Temps d'exécution médians des versions du *kernel*</figcaption></figure>

La figure 5 montre les temps d'exécution médians, sur 33 exécutions, de différentes implémen-
tations du *kernel* de calcul des forces internes pour une seule itération. La version f90 est la
version originale du code, écrite en Fortran, et utilisant un unique cœur CPU, celle qui sert de
référence, la version AVX est la version vectorisée manuellement en C/C++ utilisant également
un cœur CPU et son unité AVX, et les deux dernières implémentations sont celles en CUDA,
l'une sans la *shared memory* et l'autre avec. L'implémentation en CUDA utilisant la *shared me-
mory* obtient un speedup de 16 par rapport à la version originale en Fortran, un speedup de 7.7

<footnote n="2">2. métrique complète : *smsp__sass_average_data__bytes_perNJwavefront_mem_shared*</footnote>
par rapport à la version vectorisée CPU et un speedup de 1.3 par rapport à la version CUDA sans *shared memory*.

Dans ces versions CUDA, seul le *kernel* du calcul des forces internes a été porté sur GPU, le reste de code s'exécute toujours sur CPU, ce qui implique des copies de données entre CPU et GPU à chaque itération. Ces copies ont été prises en compte dans le calcul du speedup. Ces résultats représentent ainsi le cas le plus défavorable, celui où tous les éléments sont transférés à chaque itération. Une perspective sera de limiter ces transferts entre CPU et GPU et de les recouvrir par des calculs.

## 4.2. Validité des résultats

Indépendamment de la performance, une attention particulière a été apportée à la validation sémantique du code GPU, en effet le changement d'architecture entre CPU et GPU aboutit à des résultats qui ne sont pas totalement identiques. Ceci est du aux instructions différentes générées par les compilateurs, notamment les FMA (Fused Multiply Add), à leur organisation, ainsi qu'à l'ordre de parcours des éléments. Du fait de l'arithmétique des flottants ces modifications dans l'ordre des instructions créent des différences dans les résultats ce qui complique la validation des versions GPU. En forçant la suppression des FMA dans les deux versions et en imposant un parcours équivalent des éléments, nous obtenons le même résultat, au bit près, sur GPU et CPU ce qui valide le code parallèle. De plus la vérification des sismogrammes obtenus après une simulation complète permet de valider, auprès des géologues, les résultats obtenus avec les versions CUDA.

## 5. Conclusion

La parallélisation sur GPU du *kernel* d'**EFISPEC3D** permet d'obtenir un speedup de 7.7 par rapport au *kernel* CPU vectorisé. Alors même qu'il reste des opportunités d'optimisations, cela montre l'intérêt d'un passage sur GPU pour un code vectorisé tel que celui d'**EFISPEC3D**.
Notre implémentation permet de tirer profit de la *shared memory* au maximum mais l'occupation des unités de calcul du GPU n'est pas optimale. Les SM (Streaming Multi-processor) peuvent exécuter jusqu'à 64 *warps* simultanément, mais ce nombre peut être inférieur à cause de contraintes sur la taille des *blocks* CUDA, de la quantité de *shared memory* utilisée et du nombre de registres utilisés par *thread*. Ici l'occupation des SM n'est pas optimale à cause du nombre de registres utilisés, cette implémentation nécessite trop de mémoire par *block* CUDA.
Nous avons d'autres approches à explorer pour exploiter au mieux les ressources du GPU. En effet, un parcours réalisé par point GLL ou par degré de liberté ([5]) pourrait permettre de lancer plus de *threads* et ainsi mieux répartir la mémoire entre les *threads*. Cela permettrait d'améliorer l'occupation des SM (Streaming Multi-processor) grâce à des *warps* CUDA nécessitant moins de mémoire.

Un autre point d'amélioration concerne les accès mémoires effectués depuis le *kernel* dans la mémoire globale du GPU. Ces accès n'utilisent pas les lignes de cache (*cache lines*) à leurs maximum. Cela est dû au fait que le traitement d'un point GLL nécessite l'accès à des données correspondant à d'autres points GLL qui sont potentiellement partagés entre plusieurs éléments, lesquels ne sont pas forcément contigus en mémoire. Pour une meilleure utilisation des *cache lines*, il est possible d'envisager un changement du stockage des données et d'organisation des éléments ([11], [12]) pour rapprocher au maximum les différentes données qui seront accédées par des points GLL adjacents.
## 6. Remerciements

Ce travail a bénéficié d'une aide de l'État gérée par l'Agence Nationale de la Recherche au titre du projet intitulé « Amélioration des performances des simulations de géosciences basées sur la méthode des éléments finis spectraux par la combinaison d'optimisations de bas niveau à l'aide d'abstractions de haut niveau » (projet INSEPTION - aide ANR-23-CE46-0003).

À des fins de libre accès, une licence CC-BY 4.0³ a été appliquée par les auteurs au présent document et le sera à toute version ultérieure jusqu'au manuscrit auteur accepté pour publication résultant de cette soumission.

## Bibliographie

1. De Martin (F.). – Verification of a Spectral-Element Method Code for the Southern California Earthquake Center LOH.3 Viscoelastic Case. Bulletin of the Seismological Society of America, vol. 101, n6, 12 2011, pp. 2855-2865.

2. Fu (Z.), James Lewis (T.), Kirby (R. M.) et Whitaker (R. T.). – Architecting the finite element method pipeline for the GPU. Journal of Computational and Applied Mathematics, vol. 257, 2014, pp. 195-211.

3. Huthwaite (P.). – Accelerated finite element elastodynamic simulations using the GPU. Journal of Computational Physics, vol. 257, 2014, pp. 687-707.

4. Jubertie (S.), Dupros (F.) et De Martin (F.). – Vectorization of a spectral finite-element numerical kernel. – In Proceedings of the 2018 4th Workshop on Programming Models for SIMD/Vector Processing, WPMVP'18, WPMVP'18, New York, NY, USA, 2018. Association for Computing Machinery.

5. Kiran (U.), Gautam (S. S.) et Sharma (D.). – GPU-based matrix-free finite element solver exploiting symmetry of elemental matrices. Computing, vol. 102, n9, sep 2020, p. 1941-1965.

6. Komatitsch (D.). – Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène. – Theses, Institut de physique du globe de paris - IPGP, mai 1997.

7. Komatitsch (D.), Erlebacher (G.), Göddeke (D.) et Michéa (D.). – High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster. Journal of Computational Physics, vol. 229, n20, 2010, pp. 7692-7714.

8. Martínez-Frutos (J.), Martínez-Castejón (P. J.) et Herrero-Pérez (D.). – Fine-grained GPU implementation of assembly-free iterative solver for finite element problems. Computers Structures, vol. 157, 2015, pp. 9-18.

9. Reguly (I. Z.) et Giles (M. B.). – Finite element algorithms and data structures on graphical processing units. Int. J. Parallel Program., vol. 43, n2, avril 2015, p. 203-239.

10. Sornet (G.), Jubertie (S.), Dupros (F.), De Martin (F.) et Limet (S.). – Performance analysis of SIMD vectorization of high-order finite-element kernels. – In 2018 International Conference on High Performance Computing Simulation (HPCS), pp. 423-430, 2018.

11. Sornet (G.), Jubertie (S.), Dupros (F.), De Martin (F.), Thierry (P.) et Limet (S.). – Data-layout reorganization for an efficient intra-node assembly of a spectral finite-element method. – In 2018 26th Euromicro International Conference on Parallel, Distributed and Network-based Processing (PDP), pp. 255-258, 2018.

12. Zegard (T.) et Paulino (G. H.). – Toward GPU accelerated topology optimization on unstructured meshes. Struct. Multidiscip. Optim., vol. 48, n3, septembre 2013, p. 473-485.

<footnote n="3">3. https://creativecommons.org/licenses/by/4.0/</footnote>