11.7.5.3.1 : Un calcul pour l'astrophysique nucléaire théorique
La présentation classique du phénomène de supernova (de type II) repose sur "Quand la pression gravitationnelle au cœur de l'étoile dépasse la pression de dégénérescence des électrons, ceux-ci sont capturés par les protons" qui se transforment alors en neutrons en éjectant un neutrino... Si l'image est efficace, elle mérite d'être améliorée. En effet, il reste assez peu de protons libres au cœur d'une telle étoile. Ainsi des estimations précises requièrent de prendre en compte la structure nucléaire qui contient ces protons (typiquement un atome de fer). Le groupe d'astrophysique nucléaire de l'IPNO a donc couplé ses programmes de calcul de structure nucléaire par "Random phase approximation" à une intégration de lasection efficace $e^- + p \to n +\bar\nu_e$ ~[226]Supernovae theory: Study of electro-weak processes during gravitational collapse of massive stars, 2010, Fantina, Anthea Francesca[227]Stellar electron-capture rates calculated with the finite-temperature relativistic random-phase approximation, 2011, Y.F. Niu and N.Paar and D. Vretenar and J.Meng[228]Stellar electron-capture rates on nuclei based on a microscopic Skyrme functional, 2012, Fantina, A. F. and Khan, E. and Colo, G. and Paar, N. and Vretenar, D.[229]Stellar electron-capture rates on nuclei based on Skyrme functionals, 2014, Fantina, A. F. and Khan, E. and Colo, G. and Paar, N. and Vretenar, D..La méthode d'intégration est directe~: par parcours uniforme le long de chacune des quatre variables d'intégration / sommation. Le cœur du programme est donc un nid de six boucles dont les deux externes servent à la tabulation de lasection efficace selon l'énergie de l'électron, la parité et le numéro d'onde partielle et les quatre internes sont des intégrations sur l'angle de la réaction et la position radiale ainsi que des sommations sur les états d'énergie entrant et sortant du nucléon dans le noyau. Cette structure de boucles indépendantes se prête très bien au parallélisme, et le nombre total d'itérations des boucles internes étant entre cent millions et dix milliards (suivant la température de l'étoile), il y a assez de matière à paralléliser, par thread ou par GPU.
Après avoir profilé le programme et identifié que $92\,\%$ du temps était passé dans l'appel de fonctions de Bessel sphériques $j_\ell$ , recalculées de $150$ à $1\,500$ fois pour le même argument, le stockage des résultats intermédiaire a permis de gagner un speedup notefacteur d'accélération de $12$ . Un banc d'essai des différentes implémentations Fortran et C de $j_\ell$ a donné un autre speedup de $60\,\%$ ainsi qu'un gain de précision correspondant au passage de la simple précision à la double notecalculs initialement réalisés en double, mais avec des constantes en simple précision, entre autres.. Un travail de déclaration systématique des variables, de centralisation des common, d'uniformisation et de nommage des constantes physico-mathématiques et de documentation, a amélioré la qualité du code qui, depuis la version livrée par les physiciens, a été mis sous gestion de version distribuée mercurial. Le système de build a également été rationalisé par un Makefile paramétrable. Le choix de fonctions issues de la bibliothèque standard au lieu de fonctions bricolées par les développeurs initiaux et la vectorisation de certaines boucles ont également fourni un speedup. À l'issue de la première phase le speedup par rapport à la version initiale était de $56$ .
Le programme était alors mûr pour la parallélisation. Notre choix s'est porté sur OpenMP pour distribuer la principale boucle de tabulation sur l'énergie de l'électron. Après un travail d'identification des données communes à tous les threads, le programme compilait, mais échouait à l'exécution avec une erreur de segmentation à l'entrée de la zone parallélisée. Pour répondre à la demande des physiciens nous avions minimisé les modifications apportées au code, alors composé de $8\,000$ lignes de Fortran 77 et de $4\,000$ lignes de C. L’enquête sur ce problème passait naturellement par valgrind, mais ce dernier renonçait devant la taille du segment .bss notepartie du segment de données contenant les variables statiques représentées initialement (c'est-à-dire, quand l’exécution commence) uniquement par des bits à zéro. En Fortran, les variables de blocs COMMON sont affectées à ce segment excédant $2\,$ Go en raison de tableaux alloués statiquement avec des dimensions surévaluées codées en dur. Le débuggueur gdb ne nous en apprenait pas plus. Nos collègues de l'IDRIS nous ont recommandé de tout porter dans un unique langage. Nous avons choisi Fortran 90/95/03 pour capitaliser sur sa capacité à exprimer les expressions vectorielles. Là où la version originale ne compilait qu'avec le couple ifort/icc (Intel), ce portage fût également l'occasion de compiler le code avec trois compilateurs distincts~: gfortran, ifort et pgfortran (PGI), ce dernier en vue d'exploiter ses fonctionnalités GPU. Les avertissements des trois compilateurs nous ont donné une vision croisée des tournures à améliorer dans le code. Depuis le code a également été compilé, de manière épisodique, sous nagfor (NAG), pour chercher des avertissements instructifs, mais aussi avec le compilateur récent flang (qui repose sur llvm) notehttps://github.com/flang-compiler/flang\qquadhttps://github.com/flang-compiler/f18. Cependant l'erreur de segmentation OpenMP subsistait, et ce avec les trois compilateurs. Les trois comptes-rendus d'erreurs étaient différents, mais apparaissaient tous à l'entrée de la zone parallèle.
Nous avons alors décidé de capitaliser sur l'aspect massivement parallèle de ce code et de le vectoriser pour GPU, le principe étant de distribuer le nid de boucles sur chaque étape élémentaire.
En s'appuyant sur le portage en Fortran moderne, nous avons exploité les capacités fonctionnelles méconnues de Fortran 95 en vue de sa vectorisation~: les fonctions ont été purifiées autant que faire se peut, et le cas échéant ont été promues elemental (c'est-à-dire automappante sur les tableaux). Naturellement, les boucles ont été transformées autant que faire se peut en expressions-tableaux, allégeant considérablement le code.
La version séquentielle de ce code réécrit de manière vectorisée nous a encore offert un speedup de $2$ . Le profilage a radicalement changé par rapport à la version initiale, en révélant la dominance du produit de matrice qui était auparavant disséminé dans l'unique nid de boucle. Le recours (par le compilateur) à la fonction GEMM de BLAS pour ce produit a également fourni du speedup.
La possibilité d'écrire du code générique nous a également permis de paramétrer la précision des calculs et de déborder du cadre initial, purement en double précision. Outre l'extension à la simple précision, nous disposons également d'une version en précision étendue et en quadruple précision~: les quatre résultats sont numériquement compatibles, ce qui nous permet d'envisager le speedup bonus de la simple précision sur GPU. Cette branche du développement requiert encore un peu de mise en forme pour être ramenée dans le tronc.
Le passage sur GPU nous apparaissait (et s'est avéré) prometteur dans la mesure où il n'y a pratiquement pas d'échanges de données entre la GPU et la CPU durant le calcul, et que l'essentiel du transfert des résultats du calcul porte sur quelques kilo-octets de données~: le calcul déploie des tableaux de dimension et de taille croissantes au sein de la mémoire GPU, avant de les agréger vers des tableaux de dimension et de taille décroissantes. La partie calcul de structure nucléaire est restée intouchée (et tourne sur CPU faute d'avoir été reprise sous forme vectorielle) et fournit deux matrices complexes de quelques mégaoctets préalables au nid de boucles. C'est l'essentiel des transferts, qui n'a lieu qu'une fois. Le facteur limitant est donc cette empreinte mémoire sur la carte graphique. Des difficultés techniques pour effacer, en cours de calcul, les tableaux temporaires qui ne sont plus utiles, ajoutent malheureusement à cette empreinte. C'est une des nombreuses pistes d'optimisation.
La technologie retenue est OpenACC, qui a permis la transition du code séquentiel par addition de directives. Le compilateur Fortran de PGI, lié à NVidia, est le seul praticable pour l'exécution du code. Le compilateur gfortran de GNU s'est toutefois révélé utile à la compilation par son adhésion stricte à la norme.
Les premières exécutions sur la plate-forme ipngrid01 notePlate-forme de développement de l'IPNO ont montré une totale compatibilité des résultats, mais pas d'accélération par rapport à la CPU. À quoi il faut immédiatement ajouter que ce matériel Tesla M2090 acquis en 2013 datait de 2011. De plus, la tabulation d'une seule valeur de l'énergie de l'électron occupait déjà toute la mémoire disponible.
Nous avons bénéficié de l'appui de nos collègues de la collaboration ACP pour utiliser le serveur llracp01 dans l'environnement de développement PGI dans un container. Ainsi, non contents de bénéficier d'une technologie plus récente (architecture Volta plutôt que Fermi), nous avons disposé de trois fois plus de mémoire, qui nous a permis de tabuler simultanément cinq fois plus de valeurs, et d'enregistrer un speedup relatif de cinq sur la partie GPUifiée. Nous avons au passage fait face à une difficulté encore inexpliquée~: le programme se compile seulement dans l'environnement SLC6, et s'exécute seulement dans l'environnement CC7. Ces aspects nous confortent dans l'idée que le développement HPC a besoin d'équipes mettant en commun les compétences des administrateurs système et réseau et celles des développeurs.
Nous avons également construit un banc de test reposant sur une base de données PostgreSQL et des scripts de parsing des fichiers de résultat afin de mettre en perspective différentes versions du code, les différents compilateurs ainsi que leur version et les options de compilation.