Particules ayant une géométrie complexe

L'objectif est ici de reproduire les calculs présentés dans le tableau 1 de l'article Dauchet-et-al_2023 fourni avec ce TP (dans le répertoire doc/). Par rapport au TP1, dans lequel nous n'avons étudié que des sphères monodisperses, l'enjeu est ici de prendre en main la grammaire utilisée pour décrire la géométrie des diffuseurs. Cette géométrie est spécifiée dans un fichier yaml passé dernière l'option -i de la commande schiff (voir TP1 et man schiff). La grammaire utilisée est détaillée dans la documentation de référence :

man schiff-geometry

Récupération du matériel

Récupérer les fichiers correspondant aux géométries étudiées dans l'article (après avoir vérifié que vous êtes bien dans le répertoire TP2/) :

cp ../data/case*.yaml .

Le numéro de chaque "cas" indique la ligne correspondante du Tableau 1 de l'article Dauchet-et-al_2023.

Le fichier mis à disposition avec le TP pour les indices de réfraction peut directement être utilisé pour reproduire les résultats de l'article (seule la gamme spectrale diffère, mais les indices sont les mêmes) :

cp ../data/index.txt .

L'objectif de ce TP est de suivre les commandes données dans le Supplemental document associé à l'article Dauchet-et-al_2023. Ouvrez l'article et aller à la page 19 du pdf :

mupdf ../doc/Dauchet-et-al_opticsLetters_2023.pdf 19 &

Observer les géométries

Les figures S5, S6, S7 de l'article (page 20 à 22 du pdf) sont un rendu d'image des particules étudiées. Pour produire ces figures, on utilise d'abord l'option -G de la commande schiff qui permet d'exporter la donnée géométrique sous format .obj sans faire de calcul d'électromagnétisme (voir man schiff), puis on produit en suite un rendu d'image de cette donnée géométrique à l'aide du logiciel paraview (par exemple).

Voici une illustration pour le cas 10 (voir la figure S5) :

schiff -i case10.yaml -G 1 -o case10.obj
paraview case10.obj

Remarque : l'objet de ce TP n'est pas de prendre en main paraview ; voici quelques commandes basiques. Pour afficher la géométrie, il suffit de cliquer sur l'oeil à gauche du nom du fichier considéré (dans le cas présent case10.obj), dans le bandeau de gauche de l'interface. Puis il est possible de faire pivoter la géométrie pour l'observer sous différents angles avec la souris : maintenir le bouton gauche de la souris en la déplaçant. Un clic sur la molette permet la translation et un clic droit de zoomer sur la géométrie.

Fermer la fenêtre paraview.

Les géométries de la figure S5 n'incluent pas de distribution de la forme de la particule ; seule l'orientation et éventuellement la taille sont distribuées, par les rapports d'échelle (élongation, etc.) sont eux constants. Nous n'avons donc généré qu'un seul fichier .obj lors de l'appel à schiff. Dans les figures S6 et S7, on représente des formes distribuées : nous choisissons de représenter 35 données géométriques échantillonnées sur la distribution décrite dans case6.yaml et case11.yaml.

Utiliser les commandes données dans la légende des figures S6 et S7 afin de produire les 35 fichier .obj correspondants puis les visualiser avec paraview.

schiff -i case6.yaml -G 35 -o temp_output && csplit temp_output -z /^g\ \
/ {*} -f particle -b %d.obj
schiff -i case11.yaml -G 35 -o temp_output && csplit temp_output -z /^g\ \
/ {*} -f particle -b %d.obj

Supprimer le fichier intermédiaire généré :

rm temp_output

Pilotage de la convergence via le nombre de chemins échantillonnés

On veut analyser le temps de calcul nécessaire pour atteindre 1% de convergence sur le calcul des trois sections efficaces et sur le calcul de la fonction de phase à 1° (c'est souvent cette dernière grandeur qui est la plus délicate à faire converger). On se propose ici de piloter la convergence via l'option -d spécifiant le nombre de chemins échantillonnés par réalisation de l'algorithme de Monte Carlo, et non via le nombre de réalisations de l'algorithme de Monte Carlo que nous fixerons à 10000 (voir man schiff).

Prenons le cas 10 de la publication en exemple. Dans un premier temps, fixer le nombre de chemins à 5 :

schiff -i case10.yaml -g 10000 -d 5 -D -l 150 -a 181 -A 181 -w 0.4 \
-o results.txt index.txt

Analyser la convergence des résultats avec la même démarche que dans le TP 1 :

../schiff-tools/schiff_pretty_results.sh results.txt

Vérifier la bonne convergence des sections efficaces, en regardant si l'écart type relatif est inférieur à 1% :

more xsections.txt

Vérifier la bonne convergence de la fonction de phase à 1° (regarder la 3eme ligne de func_0.4.txt, correspondant à un angle de 0.0349066 rad), et relever l'écart-type relatif :

more func_0.4.txt

Tant que l'écart-type relatif d'une de ces grandeurs est supérieur à 1%, augmenter le nombre de chemins, et vérifier que vous retrouvez un nombre similaire à celui donné en page 19 de l'article pour le cas 10 (à savoir 63).

Un analyse précise de la convergence en fonction du nombre de chemins échantillonnés pour chaque réalisation (option -d) demande de s'appuyer sur l'article Weitz-et-al_2016, mais dans le cas présent, on remarque que l'on retrouve une propriété bien connue des estimateurs de Monte-Carlo : l'écart-type diminue proportionnellement à l'inverse de la racine carrée du nombre de chemins échantillonnés. Ainsi, si on note std5 l'écart-type initialement relevé avec 5 chemins, vous devriez au final avoir besoin de (std5/0.01)^2 chemins pour atteindre 1% de convergence. Vérifier que vous retrouvez cette propriété (à l'incertitude statistique sur la variance près). En pratique, cela signifie que le nombre de chemins à échantillonner pour atteindre 1% de convergence ne nécessite que 2 ou 3 calculs :

Appliquer la démarche précédente aux cas 1, 3, 4, 5, 6, 8... 12 de l'article Dauchet-et-al_2023 pour vérifier que le critère de convergence à 1% est bien respecté avec les commandes de la page 19 du pdf.

Maintenant que nous nous en sommes placé à 1% de convergence, nous allons pouvoir analyser les temps de calcul.

Temps de calcul

A l'exécution, la commande schiff renvoie dans la sortie standard le temps de calcul, subdivisé en 1) le temps de chargement et traitement de la géométrie (''Initialized the geometry distribution in...'') et 2) le temps de résolution du problème par la méthode de Monte Carlo (''Computation performed in...''). Attention : l'analyse des temps de calcul demande d'avoir une machine sur laquelle très peu de processus sont en cours d'exécution, sans quoi les résultats ne sont ni reproductibles, ni représentatifs. Ne pas hésiter à exécuter plusieurs fois la même commande afin de vérifier que le temps de calcul est bien reproductible (à une certaine précision).

Analyser les temps de calcul observés au regard du paragraphe Results, page 3 de l'article. Dans ce paragraphe, nous adressons le temps de calcul total - la somme des temps de chargement et de résolution - mais vous pouvez aussi étudier le poids respectif de ces deux phases. En particulier, vous devriez retrouver les résultats décrits à partir de ''(...) the number of samples and the computation time are only multiplied by 2 when accounting for orientation distribution (compare cases 2 vs. 3 and 7 vs. 8), and by 2 when further adding size distribution (compare cases 3 vs. 4 with 8 vs. 9).'' et jusqu'à la fin du paragraphe.

Analyse des propriétés radiatives obtenues

Les propriétés radiatives évaluées sont analysées dans le TP4.