Compte rendu support avancé MAR

Le code MAR (Modèle Atmosphérique Régional) est un code de calcul résolvant les équations des fluides parfaits par une méthode de différences finies.

Le but de ce projet de support avancé était d'améliorer les performances du code MAR après portage sur la machine Ada. La version de travail est une version sans toute la physique de sous-maille.

Ce projet de support avancé s'est déroulé de février à septembre 2014.

Avant d'optimiser MAR, nous avons voulu vérifier avec quelques outils qu'aucun problème ou bug n'apparaissait.

Une exécution avec Valgrind n'a pas fait apparaître d'erreurs dans l'utilisation de la mémoire (pas de débordements de tableau ou autre erreur).

Le positionnement de toutes les options de débogage du compilateur Intel ne nous a pas non plus permis d'identifier de problèmes (en dehors de quelques fausses alertes liées à des underflows).

Lors des phases d'optimisation, nous nous sommes rendu compte que pour une précision spatiale d'ordre 4, la solution n'était pas indépendante du nombre de processus. Il s'agit probablement d'une erreur dans le code qui devra être corrigée.

Dès les premières exécutions, nous nous sommes rendu compte que MAR avait besoin de beaucoup de mémoire. Par exemple, pour une exécution avec 32 processus MPI et un domaine de 578x578x30 points et pour le processus le plus gourmand, la mémoire virtuelle maximale est d'environ 5 Gio (voir le tableau TAB.1). Pour rappel, la mémoire virtuelle correspond à la mémoire allouée (réservée) par le processus mais pas nécessairement utilisée. C'est celle-ci qui est prise en compte pour les limites mémoire sur la machine Ada. La mémoire RSS (Resident Set Size) est la mémoire RAM réellement utilisée par le processus. La mémoire stack (pile en français) est la mémoire utilisée pour toutes les variables locales créées lors de l'entrée dans un sous-programme.

Ces besoins élevés en mémoire (particulièrement pour la mémoire virtuelle) étaient très contraignants sur Ada et obligeaient à s'exécuter dans les classes grosses mémoires (> 3,5 Gio) limitées à un nombre de noeuds de calcul assez restreint (environ 10 % de la machine).

Une partie non négligeable de la mémoire virtuelle provient de la bibliothèque MPI. Nous l'estimons à environ 1250 Mio. Il s'agit d'ailleurs d'un problème identifié à l'IDRIS et qui est en voie d'amélioration par IBM.

L'écart important entre la mémoire virtuelle maximale et la mémoire RSS maximale montre qu'une grande partie de la mémoire allouée n'est pas réellement utilisée par MAR. Potentiellement, une grande fraction de cet écart doit pouvoir être récupérée (en dehors de la partie liée à la bibliothèque MPI).

La diminution de l'empreinte mémoire s'est faite en suivant plusieurs pistes :

• Suppression de tous les tableaux de taille globale (en mxtt/mytt ). Chaque processus n'alloue plus que des tableaux de la taille de son sous-domaine avec ses mailles fantômes. Cette étape est essentielle car la présence de tableaux globaux est extrêmement pénalisante pour de gros maillages sur les architectures modernes ayant peu de mémoire disponible par coeur.
• Suppression de variables tableaux non-utilisées. Leur identification s'est en partie faite en utilisant l'outil DHAT de Valgrind.
• Allocation dynamique de certains tableaux seulement lorsqu'ils sont nécessaires.
• Ajustement au plus près de la taille de certains tableaux.
• Réutilisation de certains tableaux.
• Suppression de certains tableaux intermédiaires (avec comme autre avantage de ne plus devoir faire de recopies).

La consommation obtenue après optimisation pour 32 processus et un domaine 578x578x30 est présentée dans le tableau I.

Tab. 1 : Consommation mémoire en Mio de MAR pour un domaine de 578x578x30 points sur 32 processus (processus le plus consommateur, mesures obtenues en utilisant l'outil IdrisMemMPI).

L'évolution des besoins en mémoire en fonction du nombre de processus pour la même taille globale de domaine est visible sur la figure FIG.1. Pour chaque nombre de processus, l'occupation mémoire est donnée respectivement pour la version originale et la version optimisée. On vérifie bien qu'une grande partie de la mémoire qui n'était pas utilisée a été récupérée (écart entre mémoire RSS et virtuelle), la majorité de la différence restante étant due à la bibliothèque MPI. Les mémoires RSS et stack ont également été fortement réduites. On voit également que MAR peut maintenant exécuter ce cas test sans devoir passer dans les classes grosse mémoire d'Ada (> 3,5 Gio) sauf pour moins de 4 processus. La consommation mémoire de la version originale n'a d'ailleurs pas pu être mesurée au-delà de 64 processus car les classes grosse mémoire sont limitées à cette valeur.

Des réductions supplémentaires sont probablement encore obtenables mais les gains potentiels sont plus réduits. Par exemple, une fusion de la variable var_NF avec ps_dyn, ua_dyn, va_dyn… semble possible mais nécessiterait de nombreuses modifications dans les sources.

FIG.1 : Consommation mémoire en Mio de MAR pour un domaine de 578x578x30 points (processus le plus consommateur, mesures obtenues en utilisant l'outil IdrisMemeMPI).

L'optimisation séquentielle consiste à réduire le temps d'exécution du code en s'intéressant aux performances de chaque processus sans tenir compte des phases de communications entre processus.

Avant de modifier l'application proprement dite, le choix des options de compilation est important et peut déjà permettre d'obtenir des gains significatifs. Sur Ada avec le compilateur Intel (version 13.0.1), les options qui nous ont donné les meilleurs résultats sont -03 -xAVX(-O3 : optimisations aggressives, -xAVX : pour utiliser les instructions vectorielles AVX disponibles sur les processeurs Sandy Bridge d'Ada).

Pour décider dans quels sous-programmes il est le plus judicieux d'essayer d'améliorer les performances, il est nécessaire de réaliser un profilage qui permet de classer les sous-programmes en fonction du temps qu'ils consomment. Cela se fait sur un cas test le plus réaliste possible et sur une durée suffisante pour que les phases d'initialisation et de fin ne soient pas surestimées. Voici un extrait du profilage fait sur le cas test 578x578x30 exécuté sur 64 processus (résultats ici pour le processus de rang 10) :

  %   cumulative   self       self     total  
   time   seconds   seconds    calls   s/call   s/call  name 
 26.33    124.03   124.03     1200     0.10     0.12  mod_mardyn__fil_2d_mp_mardyn__fil_2d_ 
 11.28    177.17    53.14     1202     0.04     0.04  mod_marout__cdf_mp_marout__cdf_  
 11.11    229.49    52.32   288000     0.00     0.00  mod_mardyn__hyd_advhb_a_mp_mardyn__hyd_advhb_a_x_ 
 10.38    278.38    48.89   288000     0.00     0.00  mod_mardyn__hyd_advhb_a_mp_mardyn__hyd_advhb_a_y_  
 10.00    325.47    47.09    13200     0.00     0.02  mod_mardyn__hyd_ff_mp_mardyn__hyd_ff_  
  7.99    363.10    37.63    13200     0.00     0.00  mod_mardyn__hyd_advv_ec_mp_mardyn__hyd_advv_ec_ 
  4.50    384.29    21.19   288000     0.00     0.00  mod_mardyn__hyd_advhb_a_mp_mardyn__hyd_advhb_a_ 
  3.65    401.46    17.17   177600     0.00     0.00  mod_mpi_mp_solve_tri_  
  2.70    414.19    12.73   396000     0.00     0.00  mod_mardyn__hyd_ff_rk2_mp_mardyn__hyd_ff_rk2_  
  1.93    423.29     9.10     1200     0.01     0.02  mod_mardyn__hyd_ps_mp_mardyn__hyd_ps_  
  1.61    430.89     7.60     3600     0.00     0.00  mod_mardyn__dif_hv_mp_mardyn__dif_hv_  
  1.48    437.87     6.98     1200     0.01     0.01  mod_mardyn__dif_h_mp_mardyn__dif_h_  
  1.39    444.44     6.57    36000     0.00     0.00  mod_mardyn__hyd_advhb2a_mp_mardyn__hyd_advhb2a_x_ 
  1.35    450.78     6.34    36000     0.00     0.00  mod_mardyn__hyd_advhb2a_mp_mardyn__hyd_advhb2a_y_  
  1.04    455.68     4.90     1200     0.00     0.00  mod_mardyn__hyd_dgzdxy_mp_mardyn__hyd_dgzdxy_
  ...  

Au vu de ces résultats, les efforts principaux d'optimisation se sont portés sur les sous-programmes mardyn__fil_2d, mardyn__hyd_advhb_a, mardyn__hyd_ff, solve_tri et hyd__advv_ec. mardyn__hyd_advhl_a a également été optimisé car il intervenait beaucoup dans les phases de communications (voir section suivante). D'autres interventions ponctuelles ont été réalisées et ont permis de gagner quelques pourcents supplémentaires.

Le sous-programme marout__cdf était appelé à chaque itération. Son utilité est de fournir à la bibliothèque XIOS les données qu'elle pourrait avoir à écrire et il effectue des opérations de préparation des données. Or, actuellement, il ne sert qu'à écrire les champs à certains moments (pas de moyennes sur la durée). Il n'est donc pas nécessaire de l'appeller à chaque pas de temps, mais seulement lorsqu'une sortie est programmée.

Les sous-programmes mardyn__hyd_advhb_a et mardyn__hyd_advhb2a étaient très proches. La transformation consistant à travailler sur tous les plans à la fois les a rendu encore plus similaires et a ainsi rendu possible la fusion des 2.

Pour réduire le temps consommé par les différentes unités, les stratégies d'optimisation suivantes ont été appliquées (liste non exhaustive) :

• Réduction du nombre d'opérations.
• Précalcul de certaines valeurs réutilisées plusieurs fois (d'une itération à l'autre par exemple).
• Fusion de boucles pour améliorer la réutilisation de valeurs présentes dans les caches.
• Suppression de variables intermédiaires (moins de recopies).
• Sortie des invariants de boucles (même si le compilateur fait déjà une partie de ce travail).
• Suppression de certains tableaux dont les valeurs n'étaient pas réutilisées (pouvant par exemple être remplacés par un scalaire).
• Recalcul de certaines variables plutôt que de les stocker (gains mémoire, réduction bande passante mémoire, surtout intéressant pour des opérations simples et peu coûteuses).
• Déclaration de la vocation des arguments muets ( intent ) (peu de gains en performance, mais succeptible de faciliter le travail du compilateur).
• Remplacement d'opérations de division par des multiplications.

Après chaque modification, les résultats ont été vérifiés sur un petit cas test (domaine 34x34x30, 30 itérations sur 8 processus et avec les options de débogage et de vérification du compilateur). Les vérifications ont été faites avec la commande h5diff qui permet de comparer le contenu du fichier de sortie. Pour la majorité des transformations, le contenu du fichier est exactement le même (au bit près sauf pour la date du fichier). Pour certaines qui ont impliqués des modifications de l'ordre des opérations ou le précalcul de valeurs utilisées plusieurs fois, des différences au huitième chiffre signficatif peuvent apparaître pour certains points du maillage (au bout de 30 itérations).

Voici le profil obtenu avec la version optimisée dans les mêmes conditions qu'auparavant :

 
  %   cumulative   self              self     total
 time   seconds   seconds    calls   s/call   s/call  name
 16.99     54.74    54.74     1200     0.05     0.10  mod_mardyn__fil_2d_mp_mardyn__fil_2d_
 11.87     92.98    38.24     1200     0.03     0.03  mod_mpi_mp_solve_tri_yzall_
 11.61    130.39    37.41    10800     0.00     0.00  mod_mardyn__hyd_advhb_a_mp_mardyn__hyd_advhb_a_x_
  9.87    162.20    31.81    13200     0.00     0.00  mod_mardyn__hyd_advv_ec_mp_mardyn__hyd_advv_ec_
  9.68    193.40    31.20     1200     0.03     0.03  mod_mpi_mp_solve_tri_xzall_
  8.88    222.00    28.60    10800     0.00     0.00  mod_mardyn__hyd_advhb_a_mp_mardyn__hyd_advhb_a_y_
  7.07    244.77    22.77    13200     0.00     0.01  mod_mardyn__hyd_ff_mp_mardyn__hyd_ff_
  4.59    259.57    14.80    13200     0.00     0.00  mod_mardyn__hyd_ff_rk2_mp_mardyn__hyd_ff_rk2_
  2.91    268.95     9.38     1200     0.01     0.01  mod_mardyn__hyd_ps_mp_mardyn__hyd_ps_
  2.70    277.64     8.69     3600     0.00     0.00  mod_mardyn__hyd_advhl_a_mp_mardyn__hyd_advhl_a_x_
  2.68    286.27     8.64     3600     0.00     0.00  mod_mardyn__hyd_advhl_a_mp_mardyn__hyd_advhl_a_y_
  2.24    293.50     7.23    10800     0.00     0.01  mod_mardyn__hyd_advhb_a_mp_mardyn__hyd_advhb_a_
  2.17    300.49     6.99     1200     0.01     0.01  mod_mardyn__w_mp_mardyn__w_
  1.75    306.12     5.63     3600     0.00     0.00  mod_mardyn__dif_hv_mp_mardyn__dif_hv_
  1.63    311.36     5.24     1200     0.00     0.00  mod_mardyn__hyd_dgzdxy_mp_mardyn__hyd_dgzdxy_
  0.90    314.25     2.89     3600     0.00     0.01  mod_mardyn__hyd_advhl_a_mp_mardyn__hyd_advhl_a_
...

Les optimisations faites sur les sous-programmes qui étaient les plus consommateurs apparaissent bien. Cependant, certains prennent plus de temps qu'initialement (notamment solve_tri_xzall et solve_tri_yzall qui remplacent solve_tri ). Cela est dû au traitement qui ne se fait plus plan par plan mais pour tout le volume en une fois. Le but de cette modification est de réduire fortement le nombre de communications MPI (voir section suivante). La raison de cette perte de performance est que les caches mémoire sont moins bien réutilisés car les données d'un plan occupent souvent un espace suffisamment réduit pour tenir dans les différents niveaux de cache du processeur. Les différentes étapes du solveur tridiagonal peuvent ainsi réutiliser des données qui sont déjà présentes dans les caches sans devoir aller les chercher dans la mémoire centrale qui est beaucoup plus lente. Au contraire, la totalité des données du sous-domaine ne tient plus dans les niveaux de cache et doit être rechargée depuis la mémoire centrale par le processeur. Il s'agit ici d'un compromis pour lequel le choix a été fait de privilégier une forte réduction des surcoûts MPI dans le but d'améliorer l'extensibilité de MAR. Il est à noter qu'il s'agit ici d'une opération assez neutre du point de vue des performances car les gains MPI compensent les pertes de performance séquentielle (pour le cas test utilisé et pour un nombre de processus autour de 32).

Il est aussi à noter que les résultats présentés ci-dessus intègrent les optimisations des communications MPI (bien qu'elles n'apparaissent pas de manière explicite). Les optimisations MPI et séquentielles sont indissociables (sauf pour les exécutions mono-processus).

MAR dans sa version originale fait, un très grand nombre de fois, appel aux sous-programmes de communications MPI (voir tableau TAB.2). Or, une communication MPI est très coûteuse en terme de performance pour plusieurs raisons. La première est que les latences pour accéder au réseau sont de l'ordre de la microseconde sur des machines telles qu'Ada (une microseconde pour un processeur tournant à plusieurs GHz représente plusieurs milliers de cycles d'horloge) et que les débits réseau et mémoire ne sont pas illimités. Ensuite, chaque appel MPI ajoute un surcoût car de nombreuses opérations ont à être effectuées par la bibliothèque (vérification des arguments, préparation du message…). Et enfin, les désynchronisations inévitables entre processus introduisent des temps d'attente dans certaines communications (un processus en attente d'un message devra patienter si l'émetteur n'est pas prêt à lui envoyer ses données). Pour minimiser les surcoûts liés à l'utilisation d'une bibliothèque de communication comme MPI, il est donc essentiel de réduire le nombre de communications effectuées et d'essayer de les regrouper.

TAB.2:Nombre total d'appels MPI de MAR pour un domaine de 578x578x30 points et 1200 pas de temps sur 64 processus (somme sur tous les processus)

Plusieurs stratégies ont été suivies pour réduire le nombre d'appels. La première consiste à regrouper le plus possible les messages. Dans MAR, un grand nombre de communications se faisaient en échangeant des données se trouvant sur les lignes ou les colonnes de plans XY. Les sous-programmes ont été adaptés pour ne plus travailler un plan à la fois, mais pour faire le travail sur l'ensemble des plans simultanément. Cela a ainsi permis de réduire le nombre de ces communications d'un facteur correspondant au nombre de plans. Ces communications échangent maintenant toujous le même volume total de données, mais chaque message est de taille plus importante. Cela permet de réduire les surcoûts liés à MPI car on perd moins de temps dans les latences réseaux, les surcoûts logiciels de la bibliothèque MPI et dans les désynchronisations. Cette approche a néanmoins l'inconvénient de moins bien réutiliser les caches mémoire et donc de dégrader les performances séquentielles de ces sous-programmes (en tout cas, si aucune autre mesure n'est prise).

Certaines communications se sont également avérées non nécessaires. En effet, ces échanges servaient à transmettre des valeurs aux frontières qui avaient déjà été calculées localement (à partir des données précédemmment transmises).

Le tableau TAB.3 montre les gains obtenus en nombre d'appels pour les 3 sous-programmes MPI principaux. Il présente également les gains en temps pour une exécution sur 64 processus. Toutes ces mesures ont été réalisées avec l'outil de mesure de performances SCALASCA (version 2.1).

TAB.3: Temps d'exécution en secondes de MAR pour un domaine de 578x578x30 points et 1200 pas de temps sur 64 processus (somme sur tous les processus)

Attention, dans le cas d'une résolution spatiale d'ordre 4, il a été constaté que la solution dépendait du nombre de processus. Pour ne pas modifier le comportement du code, certaines optimisations n'ont pas pu être appliquées à ce cas particulier. Les performances obtenues dans la version optimisée que nous proposons ne sont donc probablement pas idéales pour cette résolution. Les optimisations faites pour les autres ordres pourront probablement y être (au moins partiellement) appliquées après correction de cette anomalie.

En dehors de l'optimisation de la consommation mémoire et des performances séquentielles et parallèles, d'autres modifications ont été apportées à MAR.

La modification la plus importante a été de rendre l'exécutable indépendant du cas test. Cela s'est fait en ajoutant au fichier MAR_control.dat 4 lignes correspondant aux dimensions du domaine et à la variable n6.

MAR utilise beaucoup de variables servant à délimiter les bornes d'accès aux différents tableaux. Certaines étaient redondantes et ont pu être fusionnées. Des simplifications complémentaires pourraient probablement encore être effectuées.

Les bornes de nombreux tableaux ont également été adaptées (0 remplacés par i0 et j0 par exemple ou réduction de leur taille). Au point de vue du langage Fortran, une partie des boucles de calcul a été réécrite en notation Fortran90 (remplacement de do-loops par des sections régulières de tableaux $a(:,:)=b(:,:)$ par exemple), les write(6,*) ont été remplacés par des write(*,*) plus portables, des if…elseif… ont été transformés en select case plus adaptés.

Certains résultats ont déjà été présentés dans les rubriques précédentes. En dehors de l'aspect gain mémoire, les gains les plus parlants sont ceux concernant les temps d'exécution du code de calcul. Les figures FIG.2 et FIG.3 montrent le temps nécessaire pour effectuer 1200 itérations sur le domaine 578x578x30, ainsi que le facteur d'accélération ( speedup ) obtenu en prenant comme référence le temps d'exécution sur 1 processus de la version originale. Globalement, on y voit que le temps d'exécution a été réduit de 40% et que l'on a une bonne extensibilité du code pour ce cas test jusqu'à 128, voire 256 processus. Au-delà, l'extensibilité n'est plus optimale et l'on voit apparaître du déséquilibrage de charge entre les processus.

FIG.2: Temps d'execution en fonction du nombre de processus du code MAR pour un domaine 578x578x32

FIG.3: Facteur d'accélération (speedup) du code MAR pour un domaine 578x578x32

Dans ce projet de support avancé, de nombreuses modifications ont été apportées au code de calcul MAR dans le but d'améliorer ses performances et réduire son occupation mémoire. Le temps imparti à un tel projet étant limité, certaines opérations n'ont pas été réalisées ou uniquement sur les parties les plus critiques du code. Ici, nous proposons une liste d'actions qui peuvent être entreprises sur le code dans le but d'améliorer encore un peu plus les performances mais aussi de rendre le code plus lisible, moderne et gérable dans le futur :

• Continuer de vérifier les bornes extrêmes des tableaux et les réduire si elles sont trop grandes.
• Supprimer ou fusionner/réutiliser certains tableaux de travail.
• Allouer les tableaux de travail seulement lorsque l'on en a besoin et les désallouer dès qu'ils ne servent plus.

• Optimiser les noyaux de calcul en vérifiant que le compilateur vectorise bien les opérations. Cela n'a pas été fait dans le cadre de ce projet de support avancé bien qu'une bonne vectorisation est essentielle pour tirer parti de l'architecture des processeurs actuels.

• Pour le moment, quel que soit l'ordre de résolution spatiale, les communications échangent une double couche de mailles fantômes. Modifier les types dérivés à l'épaisseur la plus appropriée permettrait de réduire les quantités de données échangées (gain potentiel non négligeable).
• Le nombre de communications peut probablement être encore réduit en modifiant les algorithmes utilisés. Cela implique une remise en question plus forte du code que ce qui a été fait dans ce projet où les interventions ont été relativement ponctuelles. Cette étape serait nécessaire si MAR devait être adapté à des exécutions massivement parallèles (à partir de plusieurs centaines de processus).
• Dans mardyn__fil_2d, les communications sur les différentes variables pourraient être regroupées (en définissant un nouveau type dérivé).
• Une approche hybride ajoutant un niveau de parallélisation OpenMP pourrait être envisagé dans le futur et ce pour des exécutions avec un niveau de parallélisme plus important qu'actuellement.

• Réécrire les boucles en notation Fortran90.
• Déclarer la vocation des arguments muets (intent).
• N'utiliser (use du Fortran) que les modules Fortran nécessaires dans chaque sous-programme (voire uniquement les variables et sous-programmes utilisés de ces modules).
• Simplifier et réduire les variables servant à délimiter les bornes d'accès aux différents tableaux.
• Utiliser une interface implicite pour le passage des tableaux dans les différents sous-programmes.
• Supprimer le module Fortran mod_real qui ne semble plus utilisé.
• Corriger l'indentation des fichiers.
• Autoriser des tailles de maillages non-divisibles par le nombre de processus dans une direction donnée.

Ce projet de support avancé a permis de réduire fortement l'occupation mémoire de MAR permettant à ce code d'exécuter son cas test à la plus grande résolution (578x578x30) sans avoir à passer dans les classes grosses mémoire d'Ada et ainsi à pouvoir utiliser plus de 64 processus (pas de classes grosses mémoire au-delà de 64 processus).

Les optimisations séquentielles et parallèles ont réduit d'environ 40 % le temps d'exécution du code tout en diminuant fortement le nombre d'appels MPI et en améliorant son extensibilité parallèle.

Il est ainsi maintenant possible de réaliser plus de calculs avec les mêmes ressources en heures et d'obtenir des temps de restitution nettement réduits en s'exécutant sur un plus grand nombre de processus qu'auparavant.