Gestion des projets 2019

L’objectif du cours MF210 est de vous faire découvrir la « dynamique des fluides géophysiques » (DFG) (Geophysical Fluid Dynamics (GFD) en anglais) par des cours magistraux et grâce à un projet conduit par groupe sur toute la durée du cours à l’aide d’un code numérique très simple. Le cours met l'accent plutôt sur la circulation des océans et présente ainsi quelques grands processus et phénomènes de la dynamique océanique de grande échelle, par exemple ceux vus pendant les cours MF204 et MF205, ainsi que les approximations des équations de la mécanique des fluides qui permettent de les modéliser de façon simple et les théories qui en découlent.

De façon générale, la DFG concerne l'étude des enveloppes fluides des planètes – leur atmosphère, leurs océans, leur cryosphère, voire leur magma – avec les outils de la mécanique des fluides et de la thermodynamique. Elle utilise donc des équations qui vous sont familières, mais avec des simplifications qui permettent de mieux comprendre les équilibres entre forces et les bilans d'énergie qui expliquent les écoulements observés (décrits par exemple dans le cours MF205 pour les océans de la planète Terre) d'un fluide stratifié dans un repère tournant.

Chaque groupe proposera un sujet que votre culture scientifique ou les transparents des exposés magistraux (disponibles sur le dépôt du cours) ou bien encore les livres Introduction to Geophysical Fluid Dynamics de B. Cushman-Roisin ou Atmosphere-Ocean Dynamics de Adrian E. Gill, disponibles sur le même site ou à la bibliothèque.

Les groupes disposent pour conduire leur projet d’un code en « eau peu profonde » ou « shallow water » à une seule couche écrit en Fortran. Quelques options pour activer des configurations adaptées à l’objectif du cours sont déjà codées (voir les « clefs cpp ») mais elles pourront être adaptées en fonction des sujets traités. Ici, le strict minimum à connaître en Fortran et cpp pour ce projet.

Les posters des projets réalisés les années précédentes se trouvent ici.

Par groupe de 2 ou 3, vous

• définissez ensemble un problème à traiter en vous inspirant des sujets proposés ici, ou tout autre sujet que vous pourrez proposer (nous vérifierons ensemble que le code permet bien de simuler le problème choisi),
• installez le code dans votre répertoire sur le réseau de l’école,
• réalisez quelques simulations de façon à analyser ce problème en faisant varier la configuration du modèle ou les paramètres pour comprendre les phénomènes modélisés et l'impact des erreurs numériques (c'est le cœur du projet),
• et faites une analyse des résultats.

Une synthèse du travail réalisé est faite sous la forme d’un poster présenté lors de la dernière séance. Le poster expose les résultats comme s’il s’agissait d’un poster pour un congrès scientifique. L’alternative est de rédiger un rapport de 4-5 pages qui aura la structure habituelle d’un article scientifique ou du rapport de PRe (Introduction, outils/modèles, résultats, discussion, conclusion). Ce rapport ou le poster et sa présentation servent de base à l’évaluation du cours.

En fonction du problème que vous aurez choisi de traiter, vous devrez définir une configuration de modèle, ce qui consiste à définir la géométrie du domaine, l'état initial, les conditions aux limites et un certain nombre de paramètres numériques.

La taille de la matrice qui représente de domaine est définie dans parameter.h. C’est la ligne

PARAMETER (jpi=xxx,jpj=yyy,jpk=1)

jpi est la dimension en x, jpj en y. Les changer en fonction du modèle à construire. LAISSER TOUJOURS jpk = 1.

Les coordonnées, la bathymétrie et la forme du domaine sont définies dans inigrid.h. Le modèle est écrit en coordonnées curvilinéaires orthogonales très générales, et peut donc être utilisé sur un plan en coordonnées cartésiennes ou polaires, sur une sphère, … Si vous êtes ambitieux, pourquoi ne pas utiliser cette possibilité ! La grille horizontale est définie par les tableaux glamt(jpi,jpj) (pour latitude) et gphit(jpj,jpj) (pour longitude). La grille par défaut est régulière, de pas dx en i et dy en j. dx et dy peuvent être différents. Par défaut, le point d’indice (1,1) est l'origine du repère de coordonnées (0,0). Si vous voulez faire par exemple un bassin équatorial symétrique de part et d’autre de l’équateur avec l’origine sur l’équateur à l’ouest, il faut modifier la définition de gphit en conséquence dans inigrid.h. Dans ce cas, il peut être judicieux de choisir dy plus petit que dx car on sait que certains phénomènes sont piégés au voisinage de l’équateur.

Les frontières du domaine sont définies par le tableau mbathy(jpi,jpj). Par défaut, c’est un bassin rectangulaire fermé avec mbathy(1,:) = mbathy(jpi,:) = mbathy(:,1) = mbathy(:,jpj) = 0, et 1 ailleurs. mbathy(i,j) = 0 signifie donc que le point (i,j) est un point en terre, mbathy(i,j) = 1 un point en mer. Le code contient des exemples commentés pour faire un mur, une île, …

La profondeur du domaine est définie par le tableau proft(jpi,jpj). Par défaut elle est constante. Le code contient des exemples comme un bassin avec un seuil ou une variation linéaire avec y qui peut permettre de vérifier par exemple l’équivalence avec la variation selon y du paramètre de Coriolis (voir 2ème cours). Le paramètre de Coriolis est défini par le tableau f(jpi,jpj). Par défaut il est constant. En activant, la clef cpp betaPlane, on se retrouve dans le plan « beta » (prise en compte des variations du taux de rotation de la planète avec la latitude).

L’état est défini dans inidyn.h. Il y a deux cas possibles.

On part de l’équilibre hydrostatique complet (U=0, h=0) et l’énergie des courants est apportée par le vent (condition de surface dans l’équation de quantité de mouvement) ou les conditions aux limites latérales (cas de la clef cpp obc).

Ou bien, dans l’état initial, h est variable ce qui implique que de l’énergie potentielle est disponible dans le système. La vitesse pour être nulle ou en équilibre géostrophique avec h. Les clef cpp dynWarmPool ou dynGill donnent des exemples où la hauteur d’eau h représentée par le tableau hb(jpi,jpj) est initialisée à une valeur non nulle. Si vous n'utilisez pas un des états initiaux prédéfinis, codez votre propre état initial dans le bloc dynFreeConfig et activer la clef cpp dynFreeConfig.

Si les équations ont été linéarisées (c'est le cas par défaut, la valeur numérique de hb a peu d’importance pourvue qu’elle soit nettement plus petite que la profondeur d’eau. Si vous faites des tests en non linéaire (activer la clef cpp nonlinear), il faut sans doute que hb soit suffisamment grand pour que les effets non linéaires soient visibles. La clef cpp enstrophy change le schéma numérique d'advection qui conserve alors l'enstrophie (carré du rotationnel) totale (par défaut c'est l'énergie qui est conservée).

Les conditions aux limites sont de 3 types.

• Latérales : En général, sur les frontières du domaines, on prend U.n=0 où n est la normale à la frontière (frontière imperméable et glissement), soit U=0 (frontière imperméable et condition de frottement). Le paramètre islip permet d'imposer la condition aux limites U.n=0 (islip = 1) ou U=0 (islip = 0) (par défaut la condition est le glissement U.n=0). On peut également imposer une vitesse, représentant par exemple une entrée d'eau (fleuve, détroit, …). C'est plus délicat à faire en raison de la nature des équations, mais c'est faisable avec quelques précautions (voir la clef obc).
• Au fond : Ce sont les mêmes conditions de frottement ( xlamda != 0., voir dans inipar.h) ou de glissement (xlamda = 0.).
• A la surface : Le vent, variable dans l'espace et dans le temps le cas échéant doit être spécifié. Par défaut, il est nul en permanence. La clef windStress permet d'imposer un vent. Par défaut, il s'agit d'un vent zonal (c'est-à-dire Ouest-Est) constant dans le temps qui varie selon la latitude, idéalisant les alizées au sud du bassin et les vents d'ouest au Nord.

Ils sont définis dans le fichier namelist.txt mais ont tous des valeurs par défaut. Par défaut, le code est sans dimension : la gravité vaut 1, le paramètre de Coriolis vaut 1 ainsi que la profondeur d’eau moyenne. Dans ce cas dx=dy=0.2, ce qui veut dire qu’il faut 5 pas d’espace pour résoudre le « Rayon de déformation de Rossby » qui vaut 1 lui aussi. Avec la clef cpp realWorld, on se situe sur le plan f à 45° nord (soit f=10^-4s^-1), g=10 ms^-2 et la profondeur d’eau 1000m tandis que la résolution horizontale est de 50km. Si vous changez ces valeurs, pensez à changer aussi dx et dy pour avoir toujours quelques pas d’espace pour résoudre le « rayon de déformation ».

Le pas de temps doit satisfaire le critère de Courant-Friedrichs-Lewy (CFL) soit $$U \frac{dt}{dx} < C_{max}$$ où $U$ est une vitesse caractéristique de l’écoulement (y compris la vitesse de phase des ondes) et $C_{max}$ une valeur qui dépend des schémas numériques utilisées. Dans la pratique, on prendra $C_{max}/U$ de l’ordre de 10, soit $dx / dt > 0.1$. Si avec cette valeur le code est instable, il faut diminuer dt ou augmenter légèrement la diffusion numérique xnu. Un critère de stabilité s'applique aussi à xnu qui ne doit pas être trop grand (par défaut à xnu=dxdy/800dt m2s-1, valeur qui satisfait largement ce critère).

La durée de la simulation est définie par le nombre de pas de temps nitend et la fréquence d’écriture des résultats par nwrite pas de temps.

Ce code fait appel au compilateur gfortran et pour l’écriture des résultats à la bibliothèque netcdf. Le format autodescriptif netcdf est reconnu par beaucoup de logiciels comme matlab, python, ferret …

Le fichier makefile permet de compiler le code. C'est dans le fichier makefile que vous définissez les clefs cpp qui permettent d'activer telle ou telle option du code. Il suffit pour cela d'éditer la ligne CPPFLAGS = … (attention aux majuscules/minuscules).

CPPFLAGS      = -DzonalFlow  -DbetaPlane -DrealWorld
#
# ensemble des options cpp disponibles
#
#   -DdynRestart : Pour continuer une simulation à partir de nitend
#   -DdiaTrends  : Pour calculer dis diagnostics integraux (energie, enstrophie, ...)
#   -DwindStress : Active le forçage par le vent (à définir dans inidyn.h)
#   -DrealWorld  : Pour choisir des dimensions typiques de l'Atlantique
#   -Dsill       : Pour choisir la configuration de l'article "Herbaut, 1998"
#   -DbetaPlane  : Pour être dans le plan beta (plan f par défaut)
#   -DdynWarmPool   : Pour choisir la configuration 'El Nino'
#   -DdynKelvinWave : Pour initialiser une onde de Kelvin (iperio = 1)
#   -DdynRossbyWave : Pour initialiser une onde de Rossby (iperio = 1)
#   -DdynGill       : Pour choisir la configuration de l'article "Gill, 1976" 
#   -DdynOBC        : Pour introduire un courant en équilibre géostrophique en bas à gauche du domaine
#   -DdynZonalFlow  : Pour un courant franchissant une montagne (iperio = 1)
#   -Dnonlinear  : Pour conserver les termes non-linéaires dans les équations SW
#   -Denstrophy  : Pour un schéma numérique qui conserve l'enstrophie totale plutôt que l'énergie

Le code peut s’exécuter simplement avec le script ./main.sh qui positionne les bonnes variables d’environnement et fait la liste des simulations réalisées avec les clefs cpp et paramètres de chaque simulation. Deux fichiers sont créés :

• xxx.nc au format netcdf contient l'ensemble des variables du modèle tous les nwrite pas de temps.
• xxx.dg au format ascii contient des quantités intégrales (énergie cinétique, …) calculées à chaque pas de temps permettant d'analyser les sources, puits et échanges d'énergie potentielle et cinétique.

Une visualisation des résultats au format netcdf pour dégrossir peut être faite avec l'utilitaire clicodrome ncview. Il permet de réaliser facilement des gif animés. Lancer ncview avec l'option -frames. Générer une animation avec le menu de ncview, puis quitter. Une suite de fichiers frame.001.png, frame.002.png, … a été créée. Il faut la concaténer, par exemple avec convert (package imagemagick nécessaire).

[mortier@cassis mfe21_F90]$ ncview test5.nc 
[mortier@cassis mfe21_F90]$ ncview -frames test5.nc
[mortier@cassis mfe21_F90]$ convert frame.*.png test5.gif

Panoply est un autre clicodrome qui permet d'aller un peu plus loin et surtout de faire et mettre en page des figures ou des vidéos de bonne qualité. Panoply permet aussi de faire des différences de fichiers, des moyennes, … Homepage de panoply. Panoply est installé dans /home/m/mortier/panoply.

[mortier@cassis mfe21_F90]$ ~mortier/panoply/panoply.sh

Le logiciel ferret permet de faire des visualisations de bonne qualité sans gros effort de codage. Ferret permet aussi de faire des calcul sur les résultats (moyenne, variance, …). Voir par exemple ce tutoriel en français.

Une analyse plus approfondie pourra être faite avec matlab ou python. Le script matlab hoevmoeler.m donne un exemple pour visualiser les résultats avec des diagrammes dit “Hoevmoeler” (longitude-temps ou temps-latitude). Le script matlab movie.m donne un exemple pour visualiser les résultats grâce à une animation.