TP2 : insensibilité du temps de calcul aux rapports d'échelle

Conditions de réalisation

Nous supposerons que le TP est réalisé par un utilisateur travaillant depuis le système d'exploitation EDIX.

TP2 - Partie 1 : environnement de travail

Partie 1.1 : les débuts dans EDIX

Référez vous à la partie 1.1 du TP1.

Partie 1.2 : préparation de l'environnement de TP

S'il n'existe pas encore, on crée un sous-répertoire dédié au tp htrdr-atmosphere et on va dedans

mkdir tp_htrdr-atmosphere
cd tp_htrdr-atmosphere

On crée un lien symbolique (un raccourci) vers le fichier profile définissant l'environnement du TP :

ln -sf /home/lafrier/environnement_tp_htrdr-atmosphere/etc/profile profile

Pour pouvoir exécuter htrdr-atmosphere, il faudra réaliser l'opération suivante à chaque fois qu'un nouveau terminal sera ouvert :

source ~/mon_nom/tp_htrdr-atmosphere/profile

Elements techniques de compréhension

Afficher le fichier pour voir que des variables vont être modifiées.

Ces variables sont enrichies avec les chemins vers le dossier d'installation de htrdr-atmosphere :

etudiant:~/mon_nom/tp_htrdr-atmosphere$ cat \
~/mon_nom/tp_htrdr-atmosphere/profile

ou, de manière équivalente :

etudiant:~/mon_nom/tp_htrdr-atmosphere$ cat \
/home/lafrier/star-build/local/etc/htrdr-atmosphere.profile

Vérifier que les variables d'environnement ont bien été mises à jour :
etudiant:~/mon_nom/tp_htrdr-atmosphere$ echo ${CPATH}

On copie le répertoire dédié au TP2

cp -r "${TP_HTRDR_ATMOSPHERE}"/TP2/ .

TP2 - Partie 2 : insensibilité à la complexité de la géométrie du sol

Partie 2.1 : produire une image avec htrdr-atmosphere

Aller dans le répertoire de travail dédié au TP2

cd TP2

Afficher les fichiers de données :

ls data

On voit qu'on a différents types de fichiers, d'extensions :

txt : fichier décrivant les propriétés d'une atmosphère sans nuages
obj : fichier décrivant la géométrie du sol
mtls : liste des matériaux constituant le sol (un type de matériaux et une température)
mrumtl : propriétés radiatives des matériaux en fonction de la longueur d'onde
nc : champs de concentration en eau nuageuse au format netCDF
htcp : champs de concentration en eau nuageuse au format htcp

On va ici se concentrer sur les géométries de surface (voir TP3 pour les autres types de données). On affiche dans le terminal la taille mémoire occupée par le fichier data/mountain.obj, ainsi que le nombre de facettes constituant la géométrie :

ls -lh data/mountain.obj
grep "^f " data/mountain.obj | wc -l

Puis on lance le calcul d'une image avec htrdr-atmosphere :

htrdr-atmosphere -vf \
    -a data/clear_sky_opt_prop_SMLSATM.txt \
    -c data/cubic_cloud.htcp \
    -m data/Mie_LUT_Clouds.nc \
    -D 0,60 \
    -i spp=1024\
    -C pos=500,0,750:tgt=1000,1000,700:fov=60 \
    -R -g data/mountain.obj \
    -M data/mountain.mtls \
    -o image.ht

On décrit brièvement les options utilisées avec la commande :

On exécute le code htrdr-atmosphere sur l'atmosphère clear_sky_opt_prop_SMLSATM.txt (option -a) avec un nuage cubique (option -c) associée à des données de Mie (option -m)
Avec l'option -v on a choisi un mode "verbeux" avec beaucoup d'informations, indiquant notamment les géométries chargées et les temps de calcul.
Avec l'option -D on a réglé la position du soleil
Avec l'option -i on a réglé le nombre de réalisations Monte Carlo par pixels
Avec l'option -C on a choisi que la grandeur radiative simulée soit une image, pour un capteur positionné en 500,0,750 (coordonnées en mètres) et visant un point en 1000,1000,700 avec un angle d'ouverture de 60 degrés.
Avec les options -R -g et -M on décrit la géométrie de sol
Avec l'option -o on a choisi que la sortie du calcul soit écrite dans un fichier image.ht (plutôt que directement dans le terminal)

Pour avoir plus d'information sur les différentes options disponibles, on se réfèrera au manuel de la commande :

man htrdr-atmosphere

Une fois le calcul terminé, on convertit le fichier de sortie du format text au format image avec la commande htpp

htpp -i exposure=0.1 -o image.ppm image.ht

On peut alors visualiser l'image

mupdf image.ppm

Partie 2.2 : temps de calcul en fonction de la complexité du sol

On a réalisé un calcul avec de nombreuses réalisations par pixel (1024) afin d'avoir une image de bonne qualité, suffisamment convergée. Dans la suite, on fera des calculs avec environ dix fois moins de réalisations (128) par pixel pour accélérer les calculs. En fonction de la machine sur laquelle vous travaillez, les calculs peuvent rester trop longs pour pouvoir effectuer les exercices, donc n'hésitez pas à baisser encore le nombre d'échantillons par pixel.

Comme pour le TP1, on va exécuter htrdr-atmosphere en boucle, cette fois-ci pour des géométries de sol de plus en plus détaillées. Etudier le script loop_htrdr.sh et l'exécuter

bash loop_htrdr.sh

Le fichier resultats.dat contient alors cinq colonnes : nombre de facettes de la géométrie, temps de chargement et structuration de la donnée géométrique, temps de chargement et structuration de la donnée nuages, estimation et écart type du temps de calcul par réalisation.

cat resultats.dat

Et on peut regarder la figure qui a été produite par le script

mupdf figure.png

On voit que le temps de chargement de la donnée géométrique augmente au fur et à mesure que le nombre de facettes augmente, tandis que le temps de calcul n'évolue quasiment pas. On a ainsi reproduit la figure 1 de Villefranque et al. (2019), également affichée ci dessous :

figure1_villefranque_et_al_2019

TP2 - Partie 3 : insensibilité à la complexité du champ nuageux

Partie 3.1 : temps de calcul en fonction de la complexité du champs nuageux

Etudier la sensibilité du temps de calcul à la complexité des données volumiques

=> Sur le même principe qu'à la partie précédente, exécuter htrdr-atmosphere en boucle pour les différents champs nuageux qui sont dans le répertoire data/ et examiner l'évolution du temps de chargement et de calcul. Pour cela, on peut copier et adapter le script loop_htrdr.sh

cp loop_htrdr.sh loop_htrdr_clouds.sh

Aide TP

Il faut récupérer le nombre de voxels du champ nuageux à la place du nombre de facettes de la géométrie du sol. Le nom des fichiers contient cette information e.g. le fichier data/cubic_cloud_003.nc est constitué de 3x3 voxels sur l'horizontale. Tous les champs sont constitués de 160 couches sur la verticale.

Aide TP

Au sein de la boucle, il faudra utiliser la commande les2htcp pour convertir les champs nuageux fournis au format netCDF (data/cubic_cloud_*nc) au format htcp attendu par htrdr-atmosphere.

Voir la syntaxe de la commande avec les2htcp -h

Aide shell

Par défaut, les scripts écrivent toujours dans les mêmes fichiers :
resultats.dat et figure.png. Si vous voulez sauvegarder votre travail au fur
et à mesure, vous devrez renommer ces fichiers après chaque exercice afin
qu'ils ne soient pas écrasés à la prochaine exécution, ou bien les déplacer dans un répertoire spécifique. Pour renommer un fichier :

etudiant:~/mon_nom/tp_htrdr-atmosphere/TP2$ mv resultats.dat Q1_resultats.dat
etudiant:~/mon_nom/tp_htrdr-atmosphere/TP2$ mv figure.png Q1_figure.png

Pour créer un répertoire et déplacer des fichiers dedans :

etudiant:~/mon_nom/tp_htrdr-atmosphere/TP2$ mkdir sensibilite_geom
etudiant:~/mon_nom/tp_htrdr-atmosphere/TP2$ mv *.png *.dat *.ht *.out sensibilite_geom

Aide TP

Le script reponse_part3.1.sh fournit une solution à la question précédente

On a ainsi reproduit les résultats montrés sur la figure 5 de Villefranque et al. (2019, également reproduite ci-dessous :

figure5_villefranque_et_al_2019

Partie 3.2 : collisions nulles et structures accélératrices

On vient de montrer que le temps de calcul est insensible à la complexité de la donné volumique (nombre de mailles pour décrire le champs nuageux). Cela est rendu possible par l'utilisation combinée d'algorithmes à collisions nulles et de grilles accélératrices, construites au lancement du calcul. Les algorithmes à collisions nulles reviennent à échantillonner des collisions dans un champ d'extinction majorant et à rejeter le surplus de collisions (les collisions "nulles"), une reformulation qui n'introduit aucun biais dans la résolution du transfert mais permet de se ramener à un champ virtuellement homogène, et ainsi de se débarrasser de la non-linéarité de la loi de Beer, comme décrit dans cet article donnant trois points de vue sur les collisions nulles ou cet article décrivant la méthode dans le contexte du rayonnement atmosphérique. Les deux articles sont également fournis dans le répertoire ${TP_HTRDR_ATMOSPHERE}/doc/ et peuvent être ouverts sous edix :

mupdf ${TP_HTRDR_ATMOSPHERE}/doc/elhafi_three_2020.pdf 
mupdf ${TP_HTRDR_ATMOSPHERE}/doc/villefranque_path_2019.pdf

figure2_villefranque_et_al_2019

La figure 2 de Villefranque et al. (2019), reproduite ci-dessus, illustre le concept des algorithmes à collisions nulles sur un champ maillé. En l'absence de collisions nulles (fig a), l'opérateur de lancer de rayon accède à la donnée de chaque maille traversée, afin de déterminer le lieu de la prochaine collision. Cet opérateur est fondamentalement sensible à la complexité de la donnée volumique, c'est à dire que plus le champ est finement maillé, plus l'opérateur sera lent. Avec un algorithme à collisions nulles (fig b) un majorant du coefficient d'extinction est utilisé pour échantillonner le lieu de prochaine collision, et l'accès à la donnée brute se fait uniquement en ces lieux échantillonnés. La probabilité de rejeter une collision est le ratio du coefficient d'extinction local sur le majorant. A condition que le coût d'un accès à la donnée brute ne soit pas plus élevé pour un champ maillé plus finement (e.g. si le champ est régulièrement maillé, l'opérateur d'accès à la donnée pour une position donnée est insensible au nombre de mailles), alors l'opérateur de lancer de rayon sera insensible à la complexité de la donnée volumique. En effet, le nombre moyen d'accès à la donnée le long du rayon ne dépend pas du maillage original mais de la fréquence de collisions échantillonnées à partir du coefficient d'extinction majorant. Par contre, si le champ est hautement hétérogène et qu'un majorant unique est utilisé pour tout le domaine, un temps non négligeable peut être passé à échantillonner et rejeter des collisions dans les régions optiquement minces du domaine (par exemple dans la partie de l'atmosphère au dessus des nuages). D'où l'idée de construire un champ de majorant adapté aux hétérogénéités du milieu (fig c). C'est ce champ de majorant que l'on appelle structure accélératrice, au sens où elle accélère l'opérateur de lancer de rayon par rapport à un opérateur sans collisions nulles, et par rapport à un opérateur avec collisions nulles mais utilisant un unique majorant.

Dans htrdr-atmosphere, la construction des grilles accélératrices est controlée par un paramètre qui fixe le seuil en épaisseur optique à partir duquel les voxels du champ majorant ne sont plus fusionnés. Une valeur faible de ce paramètre revient à utiliser un champ de majorant qui coincide avec la discrétisation de la donnée initiale, nous ramenant alors à la situation sans collisions nulles (fig a). Une valeur élevée de ce paramètre revient à fusionner tous les voxels, nous ramenant à la situation du majorant unique (fig b). Il y a donc un compromis à trouver entre ces deux situations extrêmes.

Etudier la sensibilité du temps de calcul au paramètre de construction des grilles accélératrices

=> reprendre et adapter le script loop_htrdr.sh pour boucler sur la valeur du paramètre (option -T) au lieu des géométries. Utiliser un champ nuageux suffisamment détaillé (par exemple data/cubic_cloud_300.htcp) pour voir le temps de calcul évoluer.

Aide TP
Le script reponse_part3.2.sh fournit une solution à la question précédente

On reproduit ici les résultats montrés sur la figure 8 de Villefranque et al. (2019), reproduite ci dessous :

figure8_villefranque_et_al_2019

Etudier la sensibilité de l'occupation mémoire au paramètre de construction des grilles accélératrices

=> Utiliser le script print_memory_usage_svx.sh sur les fichiers log*out produits par la boucle précédente, par exemple

echo "# threshold [-] memoryuse [MB]" > resultats_mem.dat 
for log in log*out ; do 
  threshold=$(echo $log|sed -e "s/.*_T\(.*\).out/\1/g")
  memoryuse=$(./print_memory_usage_svx.sh $log)
  echo "$threshold $memoryuse" >> resultats_mem.dat
done