Prévision Numérique du Temps

Prévision Numérique du Temps

Soit N un opérateur non-linéaire régissant le comportement du système atmosphérique, s’appliquant à un vecteur d’état Ψ. On appelle Prévision Numérique du Temps (PNT) la discipline qui intègre temporellement de t = 0 à t = tf le système dynamique

∂Ψ/∂t = N (Ψ), (1.1)

par des moyens numériques et qui trouve le meilleur compromis parmi des composantes clés identifiées, pour répondre à un objectif donné. Ces composantes (C1)–(C5) ont été listées par Mengaldo et al. (2018) :
— le temps à calculer la solution, (C1)
— l’énergie nécessaire à calculer la solution, (C2)
— la qualité de la solution, (C3)
— la robustesse du modèle, (C4)
— la quantification des incertitudes. (C5)

La connaissance de l’état initial Ψ(t = 0,~x) est nécessaire et estimée par une discipline nommée assimilation de données (Carrassi et al., 2018). La discipline prévision d’ensemble valorise les incertitudes portant sur l’état initial et l’opérateur N , et traite notamment la composante (C5). Lorsque le domaine géographique est la planète entière, le modèle est dit global. A contrario, lorsqu’il est restreint à une partie seulement, le modèle est dit à aire limitée (LAM) et nécessite l’imposition de conditions aux limites latérales fournies par un modèle global. Un modèle de surface impose la condition à la limite inférieure.

Moyens technologiques

Les moyens technologiques utilisés pour résoudre le problème de PNT ont fortement progressé depuis la première expérience concluante menée par Charney et al. (1950). Depuis l’apparition des premiers ordinateurs et pendant les plusieurs décennies qui s’ensuivirent, la puissance de calcul augmentait rapidement suivant la loi de Moore. Depuis quelques années néanmoins, la loi de Moore est de moins en moins respectée (Asanović et al., 2006) et l’augmentation de la puissance de calcul passe d’une part par une augmentation du nombre de cœurs par processeur (CPU), et d’autre part par une augmentation du nombre de CPU travaillant en parallèles, ou par l’utilisation de nouvelles technologies pour le calcul scientifique, comme les cartes graphiques (GPU) a priori plus économes en énergie (C2) (Michalakes et Vachharajani, 2008, Müller et al., 2013).

Dans cette étude, on modélise un supercalculateur comme un ensemble de nœuds de calculs reliés en réseau. On appelle coût en communication le coût d’un échange de données entre les différents nœuds de calcul. On considère que ce coût est régi : par la bande passante, c’est-àdire le volume de communication qui peut transiter via le réseau de communication et par le temps de latence, c’est-à-dire le temps que prend le nœud pour recevoir une information depuis un autre nœud. Chaque nœud de calcul contient un certain nombre de processeurs reliés par un réseau interne et pour lequel le coût en communication est cette fois-ci considéré comme négligeable. Ces notions sont abordées plus en détail par Zheng et Marguinaud (2018). Une tâche est assignée à chaque nœud de calcul avec le protocole de parallélisation MPI (Barros et al., 1995 donnent des détails de l’implémentation de ce protocole dans le modèle IFS). Le paradigme de parallélisation en point de grille actuel, découpe le domaine géographique de sorte que chaque noeud de calcul traite une partie du domaine géographique horizontal et les colonnes verticales correspondantes (Figure 1.1a). Sur les machines actuelles, la performance totale d’un modèle dépend particulièrement de sa capacité à tourner dans un environnement très parallèle. On distingue deux types de scalabilité : la faible et la forte. La scalabilité faible consiste à mesurer le temps d’exécution d’un modèle lorsqu’on le fait tourner avec une configuration plus avancée (une résolution plus fine par exemple), et ce avec le même nombre de nœuds de calcul. La scalabilité forte consiste à mesurer le temps d’exécution du modèle pour une version donnée, en augmentant le nombre de nœuds de calcul. Des expériences de scalabilité sont régulièrement menées par de nombreux centres météorologiques en utilisant des supercalculateurs différents, comme par exemple par Müller et Scheichl (2014) et Müller et al. (2018).

Les modèles de PNT se complexifient régulièrement au gré des avancées de la recherche. Les plus fortes demandes en puissance de calcul sont généralement dues à un accroissement de la résolution, doublant environ tous les 5 ans (Wedi, 2014) et correspondant aussi à la période caractéristique du renouvellement d’un supercalculateur (Figure 1.2). Cela pose alternativement des problèmes de scalabilité faible (lorsque l’on double la résolution du modèle sur une machine donnée) et forte (lorsque l’on change de supercalculateur pour une version donnée du modèle). Par conséquent, nous visons à améliorer ces deux types de scalabilité et par la suite, on utilisera le terme scalabilité sans préciser s’il s’agit de la faible ou de la forte. La plupart des centres météorologiques actuels disposent de puissants supercalculateurs. Par exemple, en 2020, Météo-France se dote de deux nouveaux supercalculateurs nommés Belenos et Taranis, disposant chacun d’un peu moins de 300 000 cœurs de calcul et développant une puissance maximale de plus de 10 PétaFlops. Cela représente un gain d’un facteur 5 en puissance de calcul par rapport aux précédents supercalculateurs (Figure 1.2). En juin 2020, Belenos occupe la 29ème place dans le Top500 qui répertorie les 500 supercalculateurs les plus puissants au monde (https://top500.org/).

Problèmes scientifiques 

Problème multi-échelle, multi-physique

L’étude du système atmosphérique fait intervenir de nombreuses disciplines de la physique à de multiples échelles : mécanique des fluides, thermodynamique, micro-physique, rayonnement, etc. Bauer et al. (2015) montrent que les modèles de PNT se sont améliorés par l’étude de chacune des composantes atmosphériques dont nous donnons ici un bref aperçu.

Comme nous le verrons ensuite, les équations sont résolues à une certaine résolution spatiale. Cette résolution sépare la partie des phénomène résolus par le modèle, qui est appelée « dynamique », de la partie des phénomènes sous-maille qui sont paramétrisés (c’est-à-dire dont seul leur effet moyen dans la maille est considéré) via les « paramétrisations physiques ». Le raccord entre ces deux parties est fait pour garantir la « cascade d’énergie » discutée par exemple par Skamarock (2004). L’opérateur non-linéaire N du système général (1.1) est alors décomposé en la somme d’une partie dynamique D et d’une partie de paramétrisation physique P :

∂Ψ/∂t = D(Ψ) + P(Ψ). (1.2)

L’amélioration des descriptions de D et de P ou de leur couplage permet d’améliorer la qualité (C3). Or, certains phénomènes, comme la convection profonde, sont difficiles à paramétriser. Une alternative est de les traiter dans la partie dynamique en affinant la résolution horizontale. Ainsi, pour éviter le délicat exercice de paramétrisation de la convection profonde, la plupart des centres météorologiques affinent fortement la résolution pour que ces phénomènes soient traités dans la partie dynamique. Cette raison a été l’une des motivations à l’émergence du projet AROME (Seity et al., 2011) dont la résolution de 2,5 km, lors de son lancement en 2008, a permis de se passer de cette paramétrisation. Récemment, pour les mêmes raisons, des expériences à des résolutions comparables ont été menées sur certains modèles globaux (Stevens et al., 2019), même si à l’heure actuelle, compte tenu des ressources de calcul disponibles, aucun modèle global en fonctionnement opérationnel n’atteint des résolutions inférieures à 5 km. Les ressources en calcul limitées actuelles (et très probablement pour de nombreuses années encore), ne nous permettent pas de résoudre les équations de Navier-Stokes dans la partie dynamique D et les petits tourbillons qui assurent la cascade d’énergie sont paramétrisés. L’enjeu est donc de choisir un système d’équations en phase avec la technologie de son époque. La première prévision numérique réussie menée par Charney et al. (1950) utilisait le système de conservation du tourbillon barotrope qui filtrait notamment les ondes de gravité et ne modélisait que les ondes de Rossby, en adéquation avec la technologie de l’époque. Depuis, pour éviter de filtrer des ondes d’intérêt météorologique, les équations d’Euler sont désormais résolues. Néanmoins, elles autorisent la propagation des ondes acoustiques, peu utiles en météorologie et de plus, fortement pénalisantes numériquement. Ainsi, des systèmes d’équations intermédiaires ont émergé, pour filtrer ces ondes, comme le système anélastique (une formulation est par exemple donnée par Lipps et Hemler (1982)). Cette approximation ralentit les ondes de Rossby (Davies et al., 2003) et dégrade la qualité de la simulation à l’échelle synoptique mais elle est cependant parfois encore utilisée dans les modèles de mésoéchelle (Lac et al., 2018). Depuis, des systèmes d’équations plus sophistiqués, filtrant seulement les ondes acoustiques sans subir les effets indésirables de l’approximation anélastique, sont étudiés (voir par exemple Voitus et al. (2019)). Que le système soit filtré ou non, il est possible d’effectuer l’approximation hydrostatique (H) pour les phénomènes dont l’échelle caractéristique horizontale est supérieure à la hauteur caractéristique de l’atmosphère (soit environ 10 km). Lorsque cette approximation est effectuée, la vitesse verticale est déduite des autres variables du modèle. La plupart des modèles globaux actuels dont la résolution horizontale est de cet ordre de grandeur, résolvent les équations d’Euler en faisant cette approximation. L’approximation hydrostatique filtre notamment les ondes acoustiques (à l’exception de l’onde de Lamb qui se propage uniquement suivant l’horizontale) et simplifie considérablement la résolution numérique de la partie dynamique. Pour une résolution horizontale plus fine que 10 km, quelques différences apparaissent entre un modèle faisant l’approximation hydrostatique et un modèle ne la faisant pas. Par exemple, à une résolution de 8 km, quelques différences sur les champs de pluie sont perceptibles suite aux expériences menées par Janjic et al. (2001). Pour des résolutions encore plus fines, de l’ordre du kilomètre, les vitesses verticales sont sur-estimées dans les zones de convection profonde lorsque cette approximation est effectuée en l’absence de paramétrisation (Wedi et Malardel, 2010). D’autres expériences sur des cas idéalisés ont montré que cette approximation filtrait une partie des ondes orographiques dans les atmosphères stratifiées (Wedi et Smolarkiewicz, 2009). Toutefois, comme nous le montrerons dans le paragraphe 1.4.11, la qualité de cette représentation dépend également de la résolution temporelle : un pas de temps trop grand conduit à ne représenter que partiellement la contribution linéaire Non-Hydrostatique (NH).

La mauvaise représentation des ondes de gravité d’une part, et des vitesses verticales d’autre part, dégradent la qualité des prévisions des phénomènes à enjeux comme les tempêtes de pentes aval ou la convection profonde, pour lesquels la contribution non-hydrostatique est majeure. Ainsi, contrairement aux modèles globaux, les modèles à aire limitée (LAM) sont majoritairement non-hydrostatiques. Le passage des équations, faisant l’approximation hydrostatique, appelées aussi équations primitives (HPE en anglais), aux équations NH a probablement été l’un des enjeux principal de l’amélioration de la partie dynamique ces dernières années pour les LAM et le sera aussi probablement pour les modèles globaux, pour lesquels les premières expériences en NH ont été menées par exemple lors du projet DYAMOND (Stevens et al., 2019).

Table des matières

Introduction
1 État de l’art
1.1 Prévision Numérique du Temps
1.2 Moyens technologiques
1.3 Problèmes scientifiques
1.3.1 Problème multi-échelle, multi-physique
1.3.2 Stabilité et instabilité physique
1.4 Numérisation
1.4.1 Convergence numérique
1.4.2 Maillage et coordonnée
1.4.3 Grille et discrétisation
1.4.4 Filtrage spatial et diffusion numérique
1.4.5 Discrétisation explicite ou implicite
1.4.6 Modèle eulérien ou semi-lagrangien
1.4.7 Discrétisation temporelle des modèles Eulériens
1.4.8 Discrétisation temporelle des modèles semi-lagrangiens
1.4.9 Opérateur linéaire
1.4.10 Discrétisation spectrale
1.4.11 Pas de temps
1.4.12 Exemple de modèles utilisés
1.4.13 Enjeux des prochaines années
1.5 Conclusion
2 Noyau dynamique du modèle AROME
2.1 Partie numérique du modèle AROME
2.1.1 Présentation générale
2.1.2 Coordonnées et grilles
2.1.3 Équations du modèle
2.1.4 Tangent linéaire
2.1.5 Opérateur linéaire considéré
2.1.6 Variables pronostiques
2.1.7 Discrétisation spatiale
2.1.8 Discrétisation temporelle ICI
2.1.9 Inversion du problème implicite
2.1.10 Schéma semi-Lagrangien
2.1.11 Fonctionnement général du code actuel
2.1.12 Pas de temps
2.1.13 Diffusion et filtrage
2.2 Limites de la version actuelle
2.2.1 Scalabilité
2.2.2 Fortes pentes
2.2.3 Origine des problèmes
2.3 Solutions étudiées dans la thèse
2.3.1 Scalabilité
2.3.2 Fortes pentes
3 Scalabilité d’AROME points de grille
3.1 Version points de grille
3.1.1 Problème implicite
3.1.2 Dérivées horizontales
3.1.3 Choix du solveur
3.1.4 Solveur de Krylov
3.1.5 Fonctionnement
3.2 Méthodologie d’estimation de la scalabilité
3.2.1 Mesure directe de la scalabilité
3.2.2 Mesure indirecte de la scalabilité
3.3 AROME 2D points de grille
3.4 Article : évaluation de la scalabilité du modèle AROME 2D PG
3.5 AROME 3D partiellement points de grille
3.5.1 Présentation générale
3.5.2 Critère d’arrêt
3.5.3 Estimation de la qualité
3.6 Évaluation de la scalabilité du modèle AROME 3D PG
3.6.1 Paramétrage du solveur de Krylov
3.6.2 Validation sur plusieurs champs et niveaux verticaux
3.6.3 Comparaison sur cas réel fortement précipitant
3.6.4 Scores sur un mois
3.7 Conclusion
Conclusion

Lire le rapport complet