Nous nous intéressons ici à la résolution de l'équation de Poisson 3D dans le cas de la gravitation universelle. Le but est de calculer le potentiel gravitationnel d'une galaxie à disque à partir de la position des particules représentant des étoiles; ce calcul de potentiel permet par exemple dans le code de simulation BmyOGG de générer les vitesses initiales de ces étoiles.
La méthode utilisée est celle basée sur les transformées de Fourier rapides (FFT) telles qu'elles sont implémentées dans l'API GPU/OpenCL. Nous commencerons par exprimer l'équation de Poisson 3D dans l'espace de Fourier avec la fonction de Green, tout d'abord dans le cas continu puis dans le cas discret. Il s'agira ensuite de représenter au mieux la distribution de matière sur le maillage choisi grâce à la méthode "Particle Mesh" (PM) et d'effectuer les FFT permettant de calculer le potentiel gravitationnel dans l'espace de départ.
Le code source de ce projet est contenu dans l'archive suivante :
On démontre que la fonction de Green, solution de l'équation ci-dessus, est :
\begin{equation}
G(\vec{r})=-\dfrac{1}{4\pi|\vec{r}|}
\label{green-start}
\end{equation}
La solution de l'équation de Poisson s'écrit alors sous la forme suivante :
Il reste alors à effectuer une TF inverse pour retrouver $\Phi$ dans l'espace de départ, c'est-à-dire $\Phi(\vec{r})$.
Nous nous intéresssons maintenant à la résolution numérique de l'équation de Poisson : ceci nous impose de travailler dans le cas discret (en remplaçant la TF par la FFT); Notre algorithme va donc se décomposer ainsi :
Lire un fichier de donnés sur les positions initiales des particules (ici des étoiles dans une galaxie)
Définir à partir de ces données un maillage cartésien $(x,y,z)$
Discrétiser la fonction de Green $G$ sur le maillage précédemment défini
Discrétiser la distribution de matière $\rho$ sur le maillage avec la méthode d'interpolation "Particle-Mesh"
Appliquer une FFT sur la fonction de Green $G$
Appliquer une FFT sur la distribution de matière $\rho$
Appliquer une FFT Inverse sur le produit des 2 FFT précédentes pour obtenir le potentiel
La fonction de Green dans l'espace de départ s'écrit selon \eqref{green-start}. Nous remarquons qu'il existe une singularité pour $|\vec{r}|=0$. Ce problème peut être contourné en posant $G(x=0,y=0,z=0)=1$. Ci-dessous sur un maillage de taille $N=256$, c'est-à-dire sur [i][j][k] = [-128,127][-128,127][-128,127], la fonction de Green en valeur absolue pour $z=0$ ($k=128=N/2$):
Figure 1 : Fonction de Green discrétisée (en valeur absolue)
Pour représenter la distribution de matière, nous utilisons un jeu de données sur les coordonnées $(x,y,z)$ d'étoiles dans une galaxie spirale. Ces données ont été générées selon la méthode de rejet utilisée dans le code BmyOGG. Ci-dessous une représentation graphique WebGL (avec rotation souris et zoom molette) de cette galaxie spirale (10240 étoiles). Pour vérifier si votre navigateur supporte WebGL ou si vous voulez activer cette fonctionnalité, veuillez vous rendre sur Enabling WebGL.
Méthode "Particle Mesh" :
La méthode "Particle Mesh" (PM) est une méthode numérique pour déterminer les forces dans un système de particules. Les particules pouvant être des atomes, des étoiles ou les composants d'un fluide, cette méthode s'applique à de nombreux domaines comme la physique des plasmas ou l'astrophysique. Le principe de base est que ce système de particules peut être converti en un maillage (mesh) de valeurs de densité. Le potentiel est alors résolu sur cette grille.
De nombreuses méthodes existent pour convertir un système de particules en un maillage. La méthode la plus simple est d'assigner la masse de chaque particule à son point le plus proche dans la grille : c'est la méthode "Nearest-Grid-Point" (NGP). Une autre façon, nommée "Cloud-In-Cell" (CIC), consiste à modéliser chaque particule comme un cube de densité constante et chaque particule contribue à la masse de plusieurs cellules. Dans le code, les 2 méthodes sont au choix mais l'interpolation la plus réaliste reste celle de la méthode CIC.
Ci-dessous un schéma représentant, dans le cas 2D, la répartition sur les cellules de la masse d'une particule située au point $(x,y)$ (méthode CIC).
Figure 2 : Représentation du schéma "Cloud-In-Cell" (CIC)
Nous allons désormais voir comment la répartition de la masse se discrétise, en notant :
\begin{equation}
\begin{align}
& \Delta x = fmod(x,H_{x}) \\
& \Delta y = fmod(y,H_{y})
\end{align}
\end{equation}
avec $fmod()$ le modulo étendu aux floats. On obtient alors dans chacun des 4 suivants les équations \eqref{area1}, \eqref{area2}, \eqref{area3}, \eqref{area4}. Le cas 1 correspond à la figure précédente et les équations \eqref{area1} sont listées respectivement selon les surfaces
cas 1 - Si $\Delta x > \dfrac{H_{x}}{2}$ et $\Delta y > \dfrac{H_{y}}{2}$ :
Il gaut gérer aussi les cas où nous sommes sur les bords du maillage : par exemple avec $i=1$ et $\Delta x<\dfrac{H_{x}}{2}$ ou $i=N$ avec $\Delta x>\dfrac{H_{x}}{2}$, $N$ étant le nombre de points selon l'axe $Ox$. Il faut alors veiller à ne pas tenir compte des cellules extérieures au maillage.
Résultats :
Nous exposons ici les résultats pour un maillage de taille 256x256x256. La fonction de densité avec le jeu de données Positions_Disk.dat et la méthode d'interpolation "CIC" est représentée sur la figure suivante :
Figure 3 : Fonction de densité obtenue avec la méthode "CIC"
Sur cette figure, la valeur de la densité n'est pas significative car un facteur a été appliquée pour rester dans le "range" du type float.
Le calcul des FFT se fait en utilisant la librairie OpenCL clMathLibraries/clFFT. La taille maximale des tableaux est de 2^24 pour le type float (single precision) et de 2^22 pour le type double (double precision). Nous utilisons ici le type float afin de pouvoir effectuer les FFT sur des tableaux 3D de taille maximale 256x256x256.
Après avoir lu le fichier des positions des étoiles, nous définissons un maillage cartésien 3D (qui s'adapte automatiquement en fonction du min,max des positions). Les positions sont exprimées en kiloparsec (1kpc = 3,086e+21 cm). Il faut ensuite implémenter numériquement la fonction de Green comme indiqué à la partie (3) ainsi que la fonction de densité. Nous appliquons une FFT sur chacun des tableaux 3D de ces fonctions, faisons le produit de ces 2 tableaux puis terminons par une FFT inverse pour obtenir le potentiel gravitationnel sur le maillage 3D. Tous les calculs sont effectués en unités CGS : le potentiel gravitationnel final sera donc exprimé en $cm^{2}.s^{-2}$.
Enfin, il ne faut pas oublier les facteurs intervenant dans la transformée de Fourier inverse d'un produit de 2 TF. En effet, chaque TF fait apparaître un facteur $(space_{x}\,space_{y}\,space_{z})$ et la TF inverse 3D inclut le facteur $1/(N^{3}\,space_{x}\,space_{y}\,space_{z})$. Il reste donc à multiplier la solution par le facteur $(space_{x}\,space_{y}\,space_{z})$.
Etant donné la difficulté pour représenter les valeurs du potentiel gravitationnel dans un domaine 3D, nous faisons le choix de faire un "slice" plot, c'est-à-dire d'afficher les valeurs dans un espace 3D coupé selon les tranches (Oxz) et (Oyz) avec une colormap permettant de voir les valeurs du potentiel sur ces tranches.
Ci-dessous les résultats obtenus pour un maillage 256x256x256 avec le script Matlab plot_potential_slice.m :
Figure 4 : Potentiel gravitationnel avec calcul FFT/GPU sur un maillage 256x256x256
La validation du calcul fait avec la librairie OpenCL/FFT peut être effectuée de 2 manières. La première avec le calcul direct du potentiel à partir de sa définition :
\begin{equation}
\Phi(\vec{r})=\sum_{i=1}^{N}\,-\frac{Gm_{i}}{|\vec{r}-\vec{r_{i}}|}
\label{def-potential}
\end{equation}
La deuxième où l'on utilise les routines FFT multi-dimensionnelles de Matlab. Le script compute_potential_Matlab.m implémente ces 2 options.
Ci-dessous le calcul direct du potentiel sur un maillage 256x256x256 :
Figure 5 : Calcul direct du potentiel gravitationnel sur un maillage 256x256x256
Tout comme pour la version C, on commence par charger le fichier d'entré des positions, puis on discrétise la fonction de Green et la fonction de densité. Leurs FFT sont alors calculées avec la fonction "fftn" de Matlab pour le cas 3D. Sur la figure ci-dessous le potentiel obtenu :
Figure 6 : Potentiel gravitationnel avec FFT/Matlab sur un maillage 256x256x256
Finalement, les résultats de ces 2 tests de validation sont comparables à ceux obtenus avec FFT/GPU. Comparons la valeur moyenne et l'écart type entre la méthode FFT/GPU, le calcul direct et le calcul par FFT/Matlab :
Nous voyons que les valeurs restent dans le même ordre de grandeur pour les 3 méthodes. Bien qu'il puisse exister localement des différences dues à des artefacts numériques, la méthode que nous avons implémentée est validée.
ps : contribuez comme moi au projet Cosmology@Home dont le but est d'affiner le modèle décrivant le mieux notre Univers.
