TP3 : Ajout de la diffusion

Prérequis : Connaissance de base de la grammaire de star-phor et de la programmation en C.

Nous garderons la géométrie complexe et le fichier d'entrée du TP2 afin de pouvoir comparer les résultats après modifications.

L'objectif de ce TP est d'intégrer le phénomène de diffusion dans l'échantillonnage d'un chemin optique en modifiant la fonction compute_MVREA_realization du fichier src/sphor_compute_mvrea.c.

Jusqu'à présent, le milieu était considéré comme purement absorbant (avec réflexions et réfractions sur les interfaces). Nous allons désormais prendre en compte le coefficient de diffusion, défini comme le produit de la section efficace de diffusion par la concentration volumique des espèces.

Partie 1 : Modèle physique

L'objectif de cette partie est de décrire le changement dans le modèle physique afin d'aller vers l'implémentation numérique. Pour plus de détails théoriques, vous pouvez vous référer à la thèse de Nyffenegger-Péré Yaniss : 2.3.3 Formulation statistique de la luminance (https://utheme.utoulouse.fr/s/fr/item/6174).

Partie 1.1 : Coefficient d'extinction

On définit le coefficient d'extinction comme la somme des coefficients d'absorption et de diffusion (k_ext = k_a + k_s). Dans un milieu purement absorbant, seul le coefficient d'absorption intervient dans l'échantillonnage du libre parcours moyen. L'ajout de la diffusion augmente le coefficient d'extinction total, ce qui réduit le libre parcours moyen.

Lors d'une interaction, la nature de l'événement (absorption ou diffusion) est tirée au sort selon le rapport entre le coefficient associé à cet événement et le coefficient d'extinction total.

Remarque : Le code effectuant cette opération se trouve dans la fonction sample_medium_ray_interaction du fichier src/sphor_ran_scattering.c.

Algorithme : Choix du type d'interaction (sample_medium_ray_interaction)

Calcul des coefficients k_a et k_s globaux du volume à la longueur d'onde courante.
Tirage d'un nombre aléatoire uniforme s dans l'intervalle [0, 1[.
Comparaison :
- Si s <= k_a / (k_a + k_s) : Événement d'absorption (MEDIUM_RAY_INTERACTION_ABSORPTION).
- Si s > k_a / (k_a + k_s) : Événement de diffusion (MEDIUM_RAY_INTERACTION_SCATTERING).

Partie 1.2 : Événement de diffusion

Une fois l'événement de diffusion sélectionné, l'espèce responsable de l'interaction est déterminée selon sa contribution relative au coefficient de diffusion total (dépendante de sa concentration et de sa section efficace de diffusion) et une direction est échantillonnée selon une fonction de phase définie dans le fichier d'entrée.

Remarque : Voir les fonctions sample_scatterer_in_volume et sample_scattering_direction dans le fichier src/sphor_ran_scattering.c.

Algorithme : Sélection de l'espèce diffusante et changement de direction

Si l'événement retenu est une diffusion au sein d'un milieu multi-espèces, le traitement s'effectue en deux étapes : l'identification de l'espèce responsable, puis le calcul de la nouvelle trajectoire.

1. Identification de l'espèce diffusante

La probabilité p_i qu'une espèce i interagisse avec le rayon dépend de sa contribution au coefficient de diffusion total. D'un point de vue algorithmique, ce choix s'effectue par un tirage aléatoire discret :

Préparation des données : Une boucle parcourt les propriétés radiatives (prop_rad) du volume. Pour chaque espèce, son coefficient de diffusion k_s,i est calculé selon la longueur d'onde (prop_rad_compute_ks) et stocké dans un tableau dynamique (ks_list).
Tirage de Monte-Carlo : La bibliothèque star-sp configure une distribution de probabilités (ssp_ranst_discrete_setup) proportionnelle selon la formule p_i = k_s,i / somme(k_s). La fonction ssp_ranst_discrete_get tire ensuite au sort l'indice de l'espèce retenue.

2. Échantillonnage de la nouvelle direction Une fois l'espèce identifiée, la direction du rayon est déviée selon la fonction de phase propre à cette espèce.

Délégation : Cette opération passe par la bibliothèque externe star-sf. La bibliothèque se charge d'échantillonner la direction selon le modèle défini dans le fichier d'entrée, supportant actuellement 4 types de fonctions de phase : Henyey-Greenstein, Rayleigh, RDGFA et Discrete.

Partie 2 : Réalisation et nouvel événement de diffusion

Le but de cette partie est de comprendre la fonction compute_MVREA_realisation qui met à jours le poids et le poids au carré d'une réalisation d'un chemin optique et d'aller y ajouter la partie diffusion.

Remarque : Le poids au carré sert à l'estimation de la variance.

Partie 2.1 : La fonction compute_MVREA_realization

Allez lire dans le fichier sphor_compute_mvrea.c la fonction compute_MVREA_realization.

Aide à la lecture : Architecture et conventions du code

La gestion des erreurs (Macro CALL) : Cette macro exécute une fonction, vérifie si elle retourne RES_OK, et si ce n'est pas le cas, elle redirige à error:.
Les paramètres d'Entrée / Sortie (/* in */, /* out */) : Les fonctions utilisent le passage par adresse (pointeurs &) pour renvoyer des résultats. Les commentaires /* in */ et /* out */ vous indiquent quelles variables sont lues et lesquelles sont modifiées.
Le générateur de nombres aléatoires (rng) : Le pointeur rng (Random Number Generator) transporte l'état de la séquence aléatoire. Vous remarquerez qu'il est passé en paramètre à absolument toutes les fonctions impliquant un choix statistique (émission, libre parcours, choix de l'interaction, déviation). Il garantit la reproductibilité de la simulation.

Partie 2.2 : Ajouter la partie correspondant à l'événement de diffusion

Question 2.1: En vous basant sur les algorithmes décrits dans la Partie 1 et l'architecture du code, complétez la fonction compute_MVREA_realization pour gérer le cas d'une diffusion (MEDIUM_RAY_INTERACTION_SCATTERING). Vous devrez appeler les fonctions permettant d'échantillonner l'espèce puis la nouvelle direction.

**Guide :*Consultez le fichier src/sphor_ran_scattering.h pour identifier les fonctions disponibles ainsi que leurs paramètres d'entrée (/ in /) et de sortie (/ out */).

Remarque : La solution se trouve à la toute fin de ce document.

Partie 2.3 : Fichier d'entrée

Dans cette partie, il faut enrichir le fichier d'entrée pour donner la fonction de phase, la section efficace de diffusion ainsi que le facteur d'asymétrie de la fonction de phase.

Pour cela allez lire la page manuel et chercher la grammaire associée :

man star-phor-input

Lire les nouveaux fichiers de données qui sont :

more data/sigma_s.txt
more data/g_colorant.txt

Vous pouvez ouvrir le fichier d'entrée et localiser le volume réactionnel :

vim reacteur.sphin

Question 2.2: À l'aide de la grammaire décrite dans la page de manuel de star-phor-input (man star-phor-input), complétez la propriété radiative "colorant" du fichier reacteur.sphin afin de prendre en compte la diffusion. Vous devrez déclarer : la section efficace de diffusion de l'espèce, fournie dans le fichier data/sigma_s.txt ; sa fonction de phase de Henyey–Greenstein, dont le facteur d'asymétrie g est fourni dans le fichier data/g_asym_colorant.txt. Veillez à la cohérence des unités : l'unité des sections efficaces doit correspondre à celle de la concentration.

Faites un make pour compiler les modifications puis lancez le calcul comme pour les précédents TP :

cd star-phor/
make
star-phor reacteur.sphin

Question 2.2: Quels sont les différences entre ces résultats et ceux du TP2 sans diffusion ? Retrouvez la nouvelle concentration qui minimise les pertes tout en garantissant l'absence de zones d'ombre.

Solution 2.1 :

/* Sample which scattering species in the medium was responsible for the
 * scattering event */
CALL(sample_scatterer_in_volume(sphor,
      /* in */ rng, &ray, volume, wavelength,
      /* out */ &prop_rad));

/* Sample the new direction according to the species' phase function */
CALL(sample_scattering_direction(sphor,
      /* in */ rng, &ray, prop_rad, wavelength,
      /* out */ dir));

Solution 2.2 :

  prop_rad: "colorant"
    scatterer:
      concentration: 5 mol/m^3
      cross_sections:
        abs_cross_sec: sigma_a.txt nm m^2/mol
        scat_cross_sec: sigma_s.txt nm m^2/mol
    phase_fn:
      henyey_greenstein: g_colorant.txt nm
  refractive_index:
    n_real: n_etoh.txt
  sensor:
    response_function: 1

GRAMMAR :

``` ::= 'cross_sections:' [] [] ::= 'abs_cross_sec:' \\ \\

::= 'sca_cross_sec:' \\ \\

::= 'nm' | 'cm' | 'm' | 'cm^-1' # wavenumber (1/wavelength) | '1/cm'

<phase_fn> ::= 'phase_fn:' | ::= 'henyey_greenstein:' # spec-prop-file contains g ::= 'tabulated:' \ \ \ # Polar angle, # azimuthal symmetry # is considered

::= 'rad' | 'deg' | 'cos' # The cosinus of the angle, often refered as # 'mu' in radiative transfer ::= 'sr^-1'

</details>