Modèle particulaire 2D et 3D sur GPU pour plasma froid magnétisé

Transport des particules 

Dans les simulations PIC, la discrétisation des équations du mouvement est basée sur la résolution des équations différentielles ordinaires de Newton-Lorentz. L’algorithme le plus couramment utilisé et l’un des plus simples est celui connu sous le nom de saute-mouton (en anglais, leap-frog). Ainsi, la position et la vitesse des particules sont intégrées sur un pas en temps en suivant le schéma Fig. 2.2 : Figure 2.2 – Schéma représentatif de l’algorithme saute-mouton (leap-frog). La discrétisation en temps est échelonnée à chaque demi pas en temps avec le calcul de la force et des positions des particules (F,x ) et le calcul des vitesses (v). La position des particules et le champ électrostatique sont définis au pas de temps ti et sont utilisés pour calculer la force de Lorentz (cf. Eq. 2.2) agissant sur les particules. La force électrostatique est interpolée à la position de la particule à partir des valeurs obtenues sur la grille. L’interpolation linéaire requiert les deux points du maillage les plus proches par rapport à la position de la particule, dans chaque direction. Nous verrons dans la section concernant l’interpolation des champs électrostatiques et magnétiques (cf. Section 2.6), que les points nécessaires pour cette interpolation passent de deux points en 1D à quatre points en 2D et enfin à huit points en 3D. L’avancement des particules au pas de temps ti+1 se fait en connaissant les vitesses calculées au pas de temps ti+ 1 2 . Le schéma utilisé pour le calcul des vitesses est appelé algorithme de Boris (Réfs. [32, 42]). C’est un schéma numérique centré en temps, appliqué dans de nombreuses simulations PIC [33]. Il se trouve être particulièrement bien adapté pour les codes PIC. Cette technique sépare les effets du champ électrostatique et du champ 2.3 Transport des particules 46 magnétique en plusieurs étapes, au temps ti− 1 2 . Ces étapes consistent en une demie accélération due au champ électrostatique, une rotation due au champ magnétique et une deuxième demie accélération due au champ électrostatique. En considérant les vitesses calculées au pas de temps ti+ 1 2 , les nouvelles positions des particules peuvent être calculées. Le schéma saute-mouton est donné par les relations suivantes : xti+1 = xti + vt i+ 1 2 δt vt i+ 1 2 = vt i− 1 2 + δt Fti (2.6) O`u Fti est la force interpolée à la position de la particule au début du pas de temps ti . x et v sont respectivement la position et la vitesse de la particule et δt correspond à l’intervalle de temps compris entre ti et ti+1. 

 Implémentation sur GPU

 Les cartes graphiques actuelles permettent de lancer en parallèle un nombre relativement important de threads. La solution la plus simple et la plus évidente au premier abord, pour la parallélisation des particules, est d’affecter un thread par particule. Ceci signifie qu’il est possible d’avoir autant de threads que de particules pour le calcul des nouvelles positions à chaque pas de temps. Cependant, il existe d’autres possibilités, telle que celle utilisée dans l’article de P. Abreu (Réf. [34]), o`u l’on définit des groupements de particules et o`u chaque groupe est administré par un seul thread. Ceci peut avoir un avantage dans le cas d’un nombre de particules élevé et peut permettre de mieux gérer l’occupation mémoire, l’utilisation des registres (mémoire à accès directe) ou encore la bande passante. Enfin, dans le but d’obtenir une occupation maximale du GPU, il a été estimé que le nombre de threads par block est de l’ordre de 64 ou de 128 threads. Ceci est également corroboré par plusieurs travaux de simulations PIC sur GPU . Ce nombre reste néanmoins dépendant de l’usage de la mémoire, du type de carte graphique utilisé et de la variation du nombre total de particule lors de la simulation. Ces nombres ne doivent donc pas être considérés comme invariants et dépendent sensiblement du problème et des outils utilisés. Dans certaines simulations, lorsqu’il n’y a pas de créations ou de pertes de particules chargées (e.g. dans le cas de conditions limites périodiques et o`u les réactions d’ionisation n’interviennent pas), chaque particule possède son propre indice durant toute la simulation et leur suivi, i.e la fa¸con d’ordonner les particules dans le tableau, ne pose pas de problème significatif. Chapitre 2. Méthode PIC et Parallélisation sur GPU 47 Au contraire, lorsqu’il y a création de particules chargées (e.g. lors de collisions ionisantes) ou des pertes (e.g. absorptions par les électrodes ou les parois du diélectrique, processus d’attachement, . . .), les indices des particules dans le tableau doivent être réorganisés à chaque pas de temps. Les pertes de particules induisent une diminution du nombre total de particules et font apparaître des ”trous” dans les tableaux de positions et de vitesses, i.e. des emplacements mémoire du tableau se libèrent. Le but de la méthode présentée ici est de réordonner les particules dans le tableau de telle sorte que les particules toujours présentes dans le système soient rangées de fa¸con consécutive. Dans le cas des créations de particules, la méthode consiste juste à déposer les nouvelles particules à partir du dernier élément de la liste des positions et des vitesses et est présentée dans la section 2.7. Dans notre cas, la perte des particules à lieu seulement au niveau des parois du domaine de simulation. La méthode suivante peut être utilisée de la même fa¸con pour d’autres types de pertes. Lors du réarrangement en parallèle des particules dans les tableaux (de positions et de vitesses), à chaque pas en temps et afin d’éviter l’attribution d’un même indice à deux particules différentes, on utilise des opérations de préfixages couramment appelées ”Prefix Sum” ou encore ”Scan” . Cette méthode est décrite dans ce qui suit. Après que les particules aient été déplacées au cours d’un pas de temps δt donné, les positions de chacune d’entre elles dans le domaine de simulation sont vérifiées puis une valeur indiquant si la particule réside toujours dans le système ou non est inscrite dans un tableau que l’on nomme ici ”pfx”. Une valeur de 1 est assignée à la particule si celle-ci est toujours dans le domaine de simulation, sinon on lui assigne la valeur nulle. pfx[i] = [a0, a1, . . . , an−1] avec ( ai = 0 ou 1 et i = 0, . . . , n n = nombre total de particules (2.7) Les valeurs du tableau de préfixes (pfx ) sont, dans une seconde étape, utilisées afin d’effectuer la somme de ces préfixes. Cette étape consiste à incrémenter la valeur du i eme élément du tableau par la valeur de l’élément i − 1. Ainsi, si la particule réside dans le domaine de simulation, la valeur est incrémentée de 1, sinon elle ne change pas. Cette somme peut être schématisée par la relation donnée ci-dessous : pfxSum[i] = [a0,(a0 ⊕ a1), . . . ,(a0 ⊕ a1 ⊕ . . . ⊕ an−1)] avec i = 0, . . . , n O`u ⊕ correspond à l’opérateur d’addition binaire (opérateur d’addition en Mathématique). NVIDIA fournit un algorithme de scan optimisé qui peut être utilisé pour réaliser 2.3 Transport des particules 48 cette opération appelée ”Prefix Sum” . Comme présenté dans la figure 2.3, par l’opération de préfixage puis la somme des préfixes (e.g. opération de scan), les nouveaux indices des particules se déduisent facilement du tableau pfx et pfxsum. Figure 2.3 – Illustration de l’opération de réarrangement des particules (les 0s dans le tableau de Prefix correspondent aux particules situées à l’extérieur du domaine de simulation). La nouvelle position des particules dans le tableau, i.e. le nouvel indice de la particule i, correspond à la valeur de l’élément i du tableau pfxSum. Ceci n’est vrai que pour une valeur non nulle de l’élément i du tableau pfx. Listing 2.1– Kernel : changement d’indices 1 for ( chaque particule Pi) 3 { int index = pfxSum [i ]; /* Nouvel indice */ 5 if ( pfx [ i] == 1) /* pour les particules dont la valeur */ { /* de pfx[i] = 1 */ 7 NewPos [ index -1] = pos [i ]; NewVel [ index -1] = vel [i ]; 9 } } 11 Npart = pfxSum [n -1] /* Nombre de particules restant */ /* avec (n -1) = dernier element de la liste */ Finalement, toutes les particules quittant le domaine peuvent être retirées du tableau de particules ou indexées à la fin de celui-ci. De plus le nombre de particules restant est également déduit par la lecture du dernier élément du tableau pfxsum. Un algorithme simplifié est présenté List. 2.1. 

 Densité de charge

 Dans les simulations PIC, les particules chargées doivent être assignées sur une grille à chaque itération en temps dans le but de calculer la distribution spatiale de charge. Une grille rectilinéaire est généralement utilisée afin d’obtenir une représentation discrétisée du champ sur le domaine de simulation par le biais de la résolution de l’équation de Poisson. La stratégie couramment employée dans les simulations PIC consiste à interpoler la charge des particules sur les noeuds de la grille. Au pas de temps δt, connaissant la position de chaque particule, on en déduit la cellule dans laquelle elle se trouve. La charge de chaque particule dans une cellule est distribuée sur les quatre noeuds (en 2D) définissant la cellule (schématiquement représenté Fig. 2.4(a)), en tenant compte de la distance entre la particule et les noeuds. a) b) Figure 2.4 – (a) Assignement de charges sur les noeuds. Une particule dans une cellule contribue à la densité aux 4 noeuds définissant la cellule (en 2D), (b) Illustration de l’interpolation linéaire de la densité par la formule des ”aires opposées”. Chaque Aij multiplié par le coefficient de poids wi correspond à la contribution de la particule pi pour la fonction de densité de particule au noeud sur la diagonale opposée. 

 Méthode d’interpolation 

De manière générale, l’opération d’interpolation d’un ensemble de particules {p1, p2, . . . , pNT } sur la grille correspond au calcul de la fonction de distribution de la densité sur une grille rectilinéaire uniforme, dont la valeur de la fonction de distribution de densité f à chaque noeud ns est donnée par : f(ns) = X NT i=1 wiK(ns, pi) (2.8) 2.4 Densité de charge 50 O`u s est un indice multiple 2 , wi sont des coefficients de poids des particules et K correspondent à la fonction d’interpolation, i.e. à un facteur de forme. A l’ordre 0 le facteur de forme ` K à deux dimensions est défini comme : K0 ((0, 0),(x, y)) = ( 1, si |x| < 1 2 , et |y| < 1 2 0, sinon (2.9) o`u l’on assigne la charge de la particule au noeud le plus proche. Cette méthode est plus généralement appelée Nearest Grid Point (NGP). A l’ordre 1 le facteur de forme ` K à deux dimensions correspond à une interpolation linéaire définit comme : K1 ((0, 0),(x, y)) = ( (1 − |x|)(1 − |y|), si |x| < 1, et |y| < 1 0, sinon (2.10) Dans la plupart des applications, K est une fonction d’interpolation linéaire qui, à une dimension, est représentée par une fonction ”chapeau”. L’utilisation de ce type d’interpolation linéaire implique que pour chaque noeud ns, seules les particules contenues dans des cellules incidentes à ns verront leurs contributions sommées dans l’Eq. (2.11) : f(ns) = X pi∈P(ns) wiK(ns, pi) (2.11) O`u P(ns) représente l’ensemble des particules contenues dans la cellule incidente au noeud ns. Cette expression peut être interprétée géométriquement par la mesure Euclidienne du volume opposé au noeud ns par rapport à la particule pi . Soit à deux dimensions, la contribution de la particule pi peut être calculée par la formule des aires opposées Eq. (2.29) et représentée schématiquement Fig. 2.4(b).