2 · Calcul scientifique¶

Le premier domaine à avoir adopté le GPU hors du graphique, et celui où les contraintes de précision sont les plus sérieuses.

2.1 Les trois structures de problème¶

La quasi-totalité de la simulation numérique se ramène à trois structures, dont les caractéristiques GPU diffèrent radicalement.

Structure	Exemple	Régularité	Accélération GPU
Stencil	équation de la chaleur, ondes, fluides explicites	très régulière	excellente
Particules	dynamique moléculaire, N-corps, astrophysique	moyenne	très bonne
Systèmes linéaires	éléments finis, fluides implicites	dépend du solveur	variable

2.2 Les stencils¶

Un stencil met à jour chaque point d'une grille en fonction de ses voisins :

\[ u^{n+1}_{i,j} = u^n_{i,j} + \alpha \left( u^n_{i+1,j} + u^n_{i-1,j} + u^n_{i,j+1} + u^n_{i,j-1} - 4u^n_{i,j} \right) \]

C'est le cas GPU idéal : parfaitement parallèle, accès réguliers, arithmétique simple.

__global__ void laplacien(const float* u, float* u_new,
                          int nx, int ny, float alpha) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    if (i < 1 || i >= nx-1 || j < 1 || j >= ny-1) return;

    int idx = j * nx + i;
    u_new[idx] = u[idx] + alpha * (u[idx+1] + u[idx-1]
                                 + u[idx+nx] + u[idx-nx] - 4.0f*u[idx]);
}

L'intensité arithmétique d'un stencil à 5 points en float :

opérations : ~6 par point ;
octets : 4 en lecture (chaque point est lu 5 fois, mais avec une bonne réutilisation en cache ou en mémoire partagée, on approche 1 lecture par point) + 4 en écriture ;
$I \approx 6/8 = 0{,}75$.

Contre $I_{\text{crit}} = 20$ sur H100 en FP32 : massivement limité par la mémoire, d'un facteur 27.

Ce que cela implique

Optimiser un stencil, c'est exclusivement optimiser les accès mémoire.

Les techniques rentables, par ordre :

le pavage temporel (temporal blocking) : effectuer plusieurs pas de temps sur une tuile gardée en mémoire partagée, avant de la réécrire. Cela multiplie $I$ par le nombre de pas fusionnés ;
le pavage spatial en mémoire partagée avec halo ;
le stockage en 2,5D (garder trois plans en registres pour un stencil 3D) ;
la précision réduite quand la physique le permet.

Le pavage temporel est le levier le plus puissant : il transforme un problème limité par la mémoire en problème limité par le calcul. Il complique le traitement des bords et la parallélisation en domaine, ce qui explique qu'il ne soit pas systématique.

2.3 Les particules et la dynamique moléculaire¶

Le problème à $N$ corps : chaque particule interagit avec les autres.

Interactions à courte portée (van der Waals, liaisons) : on utilise une liste de voisins ou une grille de cellules, et le coût devient $O(N)$. C'est très bien adapté au GPU.

Interactions à longue portée (électrostatique) : le coût naïf est $O(N^2)$. On utilise :

PME (Particle Mesh Ewald) : décompose en une partie courte portée dans l'espace réel et une partie longue portée dans l'espace de Fourier, résolue par FFT. Coût $O(N \log N)$ ;
les méthodes multipolaires rapides (FMM) : $O(N)$, mais avec une constante élevée et une structure d'arbre irrégulière.

Le noyau de PME contient donc des FFT (cuFFT), ce qui en fait un cas intéressant : la partie longue portée est limitée par la bande passante et la communication, la partie courte portée par le calcul.

Les codes de référence — GROMACS, AMBER, NAMD, LAMMPS, OpenMM — sont tous accélérés GPU depuis une décennie. Les accélérations rapportées vont de 10× à 100× contre un nœud CPU, la variabilité venant surtout de la qualité de la référence CPU.

Le point critique : la précision

Une trajectoire de dynamique moléculaire intègre des équations différentielles sur des millions de pas. Les erreurs s'accumulent.

La pratique courante est la précision mixte : forces en simple précision, accumulation et intégration en double. GROMACS et AMBER le font depuis longtemps, avec validation croisée contre des références en double précision.

Ce n'est pas une optimisation qu'on improvise : il faut valider sur des grandeurs conservées (énergie totale, quantité de mouvement) sur des trajectoires longues.

2.4 La mécanique des fluides¶

Deux familles, avec des profils GPU opposés.

Solveurs explicites¶

Le pas de temps est calculé directement depuis l'état courant. Structure de stencil, très bien adapté au GPU.

Contrainte : la condition CFL impose un pas de temps petit, donc beaucoup d'itérations. Mais chaque itération est rapide.

Codes : MFC/MFC 5.0, AMR-Wind, et de nombreux solveurs de recherche. Le solveur multi-physique MFC 5.0 est explicitement conçu pour l'exascale sur GPU.

Solveurs implicites¶

Chaque pas de temps demande de résoudre un grand système linéaire creux.

C'est là que ça se complique :

Composant	Adaptation GPU
Produit matrice-vecteur creux (SpMV)	correcte, limitée par la mémoire et les accès irréguliers
Gradient conjugué, GMRES	correcte (composé de SpMV et de réductions)
Préconditionneur	c'est le problème
Factorisation LU/ILU incomplète	intrinsèquement séquentielle
Multigrille algébrique	partiellement parallélisable, phase de construction difficile

Le préconditionneur est le nœud. Un ILU est séquentiel par nature (dépendances en cascade) ; les variantes parallèles (Jacobi par blocs, polynomial, Chebyshev) convergent moins bien, ce qui annule une partie du gain.

Résultat pratique : les accélérations sur solveurs implicites sont de 2× à 10×, contre 20× à 50× sur les explicites. C'est le meilleur exemple du principe selon lequel l'algorithme optimal sur CPU n'est pas l'algorithme optimal sur GPU : sur GPU, un solveur qui converge moins vite mais se parallélise mieux peut gagner.

OpenFOAM avance vers l'exascale par accélération GPU des solveurs linéaires creux itératifs et vectorisation des solveurs d'écoulement réactif. Ansys Fluent a développé des architectures de solveur natives GPU précisément pour contourner les limites imposées par la loi d'Amdahl sur les portages partiels.

2.5 Le climat et la météorologie¶

Un domaine où le GPU a mis longtemps à s'imposer, pour deux raisons :

les codes sont en Fortran, écrits sur des décennies ;
ils sont dominés par des paramétrisations physiques (nuages, rayonnement, turbulence) pleines de branchements.

Les voies de portage :

Approche	Effort	Performance
OpenACC	directives, faible	60-80 % du natif
OpenMP target	directives, faible	similaire
DSL dédiés (GridTools, PSyclone)	moyen	bonne
Réécriture (CUDA/HIP)	élevé	maximale

OpenACC a été conçu pour ce cas d'usage : ajouter des directives à du Fortran existant.

!$acc parallel loop collapse(2)
do j = 2, ny-1
  do i = 2, nx-1
    u_new(i,j) = u(i,j) + alpha * (u(i+1,j) + u(i-1,j) &
                                 + u(i,j+1) + u(i,j-1) - 4*u(i,j))
  end do
end do
!$acc end parallel loop

Une tendance forte de 2024-2026 : le remplacement de paramétrisations physiques par des réseaux de neurones, entraînés sur des simulations à haute résolution. Cela déplace le problème vers un régime où le GPU excelle.

2.6 Le problème du FP64¶

C'est la contrainte la plus concrète du domaine.

Carte	FP64 (TFLOPS)	Rapport FP64/FP32
MI300X	81,7	1/1
B200	~40	~1/2
H100	34	1/2
A100	9,7	1/2
RTX 5090	~1,6	1/64
RTX 4090	1,3	1/64

Le piège de la carte de station de travail

Une RTX 5090 est plus rapide qu'un A100 en BF16 et 25 fois plus lente qu'un H100 en FP64. Le bridage est délibéré et commercial.

Conséquence pratique : vérifiez le débit FP64 de votre cible avant de dimensionner un projet scientifique. Un code de simulation en double précision développé sur une RTX puis déployé sur H100 verra un gain énorme ; l'inverse est une catastrophe.

Les stratégies de contournement :

Précision mixte : calculer en simple, accumuler en double, corriger itérativement (iterative refinement). Peut donner la précision du FP64 à un coût proche du FP32.
Compensation de Kahan : réduire l'erreur d'accumulation sans changer de type.
Analyse d'erreur : déterminer si le FP64 est réellement nécessaire. Il l'est souvent moins qu'on ne le croit, mais il faut le démontrer, pas le supposer.

Et le rappel de la partie 1 : désactivez TF32 si vous partez de FP32 et que la précision compte.

torch.backends.cuda.matmul.allow_tf32 = False

Résumé du chapitre¶

À retenir

Trois structures : stencils (idéal), particules (très bon), systèmes linéaires (variable).
Un stencil a $I \approx 0{,}75$ : il est limité par la mémoire d'un facteur ~27 sur H100. Le seul levier qui change le régime est le pavage temporel.
En dynamique moléculaire, les interactions à longue portée passent par PME (donc par des FFT) ou par les méthodes multipolaires.
En fluides, les solveurs explicites accélèrent 20-50×, les implicites 2-10× — le préconditionneur est le goulot.
Sur GPU, un algorithme qui converge moins bien mais se parallélise mieux peut gagner. L'optimum CPU n'est pas l'optimum GPU.
Vérifiez le débit FP64 de votre cible : facteur 25 entre carte grand public et carte de centre de données.

Vérifiez que vous avez compris¶

Pourquoi le pavage temporel change-t-il le régime d'un stencil ?

Parce qu'il fait $k$ pas de temps sur une tuile chargée une seule fois.

Sans pavage temporel : $I \approx 0{,}75$ (6 opérations pour 8 octets).

Avec $k$ pas fusionnés : $6k$ opérations pour les mêmes 8 octets (plus le coût des halos), donc $I \approx 0{,}75k$. Avec $k = 30$, on atteint $I \approx 22$ et on franchit le seuil FP32 de H100.

Le prix : il faut charger un halo de $k$ points de large autour de la tuile (les valeurs dont dépend le résultat après $k$ pas), ce qui ajoute du calcul redondant. L'optimum est un compromis, généralement $k$ entre 2 et 8.

Pourquoi un préconditionneur ILU est-il mauvais sur GPU ?

Parce que sa résolution triangulaire est intrinsèquement séquentielle : calculer $x_i$ exige de connaître $x_1, \dots, x_{i-1}$. Le parallélisme disponible est limité par la structure de dépendance de la matrice.

Les techniques d'ordonnancement par niveaux (level scheduling) extraient du parallélisme en identifiant les inconnues indépendantes, mais le nombre de niveaux est souvent grand et le parallélisme par niveau faible.

Les alternatives GPU : Jacobi par blocs, préconditionneurs polynomiaux, Chebyshev, multigrille algébrique. Toutes convergent moins bien qu'un ILU, et il faut mesurer le compromis « itérations supplémentaires » contre « itérations plus rapides ».

Un chercheur veut porter un code Fortran de 200 000 lignes sur GPU. Que lui conseiller ?

OpenACC ou OpenMP target, sans hésiter, pour trois raisons :

l'effort est en directives, pas en réécriture ;
le code reste compilable et exécutable sur CPU (validation croisée possible à tout moment) ;
on atteint typiquement 60 à 80 % de la performance d'un portage natif.

La méthode : profiler d'abord pour identifier les 5 à 10 boucles qui dominent, les annoter, puis traiter les mouvements de données — c'est là que se joue le résultat. Un portage OpenACC naïf transfère les tableaux à chaque boucle et est plus lent que le CPU. Les directives !$acc data et !$acc enter data qui gardent les données sur le GPU font toute la différence.

Chapitre suivant : 3 · Données et bases de données