Équation de Poisson : module Python

1. Introduction

Ce document présente le module poisson, qui permet d'effectuer la résolution numérique de l'équation de Poisson 2D (applications en électromagnétisme et en thermodynamique) par la méthode des différences finies :

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} = s (x, y)

où u(x,y) est la fonction inconnue et s(x,y) la fonction source, éventuellement nulle (équation de Laplace). L'équation ci-dessus s'applique à un problème en géométrie XY.

Le module permet également de résoudre l'équation de Poisson en géométrie axiale :

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial u}{\partial r}) + \frac{\partial^{2} u}{\partial z^{2}} = s (z, r) (1)

avec u(z,r) la fonction inconnue en coordonnées cylindriques, et l'équation de Poisson vectorielle en géométrie axiale :

\frac{\partial}{\partial r} (\frac{1}{r} \frac{\partial (r A_{θ})}{\partial r}) + \frac{\partial^{2} A_{θ}}{\partial z^{2}} = s_{θ} (z, r) (2)

Deux méthodes itératives de résolution sont possibles :

Méthode de Gauss-Seidel avec sur-relaxation
Méthodes multigrilles (cycle en V et multigrille complet)

Pour la méthode de Gauss-Seidel avec sur-relaxation, le module comporte un programme tournant sur le processeur central (CPU), et un programme utilisant les capacités de calculs parallèles des cartes graphiques (GPU), via l'interface de programmation OpenCL.

2. Installation du module

Fichiers wheel pour installation par pip install fichier.whl :

Source : poisson-0.3.zip
Binaire Python2.7 Win32 : poisson-0.2.win32-py2.7.exe
Binaire Python3.7 Win32 : poisson-0.3-cp37-cp37m-win32.whl
Binaire Python3.8 Win32 : poisson-0.3-cp38-cp38-win32.whl
Binaire Python3.9 Win32 : poisson-0.3-cp39-cp39-win32.whl
Binaire Python3.7 Win64 : poisson-0.3-cp37-cp37m-win_amd64.whl
Binaire Python3.8 Win64 : poisson-0.3-cp38-cp38-win_amd64.whl
Binaire Python3.9 Win64 : poisson-0.3-cp39-cp39-win_amd64.whl
Binaire Python3.11 Win64 : poisson-0.3-cp311-cp311-win_amd64.whl

Les binaires sous forme de paquets wheel sont aussi disponibles sur https://pypi.org/project/poisson/0.3. L'installation se fait avec pip. Pour cela, ouvrir une console puis aller dans le dossier Scripts de la distribution, par exemple C:\Anaconda3\Scripts, et exécuter :

pip install poisson

La bibliothèque pthreads est utilisée pour le multi-threading sur le processeur central (CPU). Cette bibliothèque est disponible en standard sur les distributions Linux. Pour la compilation sur Windows, il faut installer la version pthreads-win32. La DLL correspondante est incluse dans la distribution binaire pour win32.

Pour compiler et installer sous linux : sudo python setup.py install. Pour compiler sans OpenCL, ouvrir le fichier setup.py et enlever opencl_include et define_macros.

3. Documentation

Exemples d'utilisation :

Scripts de test :

Méthodes numériques

4. Description de l'interface

4.a. Définition du maillage

On définit un maillage réduit, qui servira à définir les objets, et un nombre de niveaux de grilles.

La classe principale est Poisson. Elle se trouve dans le module poisson.main. Son constructeur effectue l'allocation de l'espace mémoire pour le maillage :

object = Poisson.Poisson(pw_min,ph_min,levels,width,height,x0=0,y0=0)

Constructeur : définition d'un maillage 2D

pw_min, ph_min : puissances définissant les dimensions du maillage réduit. Les nombres de mailles sont 2^pw_min dans la direction X et 2^ph_min dans la direction Y.
levels : nombre de niveaux de grilles. Le maillage le plus fin est défini par les puissances pw_min+levels-1 et ph_min+levels-1.
width : largeur du domaine.
height : hauteur du domaine.
x0, y0 : coordonnées du point origine.

Le nombre de niveaux de grilles est le nombre maximal de grilles utilisables dans la méthode multigrille. On aura donc intérêt à définir le problème sur une petite grille et à fair varier le nombre de niveaux de grilles en fonction de la finesse souhaitée pour le calcul.

La figure suivante montre un exemple de maillage avec pw_min=ph_min=2 et levels=3. Il y a trois grilles disponibles pour l'algorithme multigrilles. Les objets seront définis sur la grille 4x4.

Figure pleine page

La position du point origine ne joue aucun rôle dans le calcul numérique. Elle est utilisée pour la présentation graphique des résultats. Pour le calcul numérique, seul importe le rapport de la largeur sur la hauteur. Par exemple, si la largeur est deux fois plus grande que la hauteur, et si le nombre de mailles est identique sur les deux dimensions, alors le pas d'espace est deux fois plus grand dans la direction X. En revanche, ces dimensions interviennent dans le calcul des dérivées.

Une seule instance de cette classe est utilisable, car le programme C sous-jacent ne comporte qu'une instance des données.

Poisson.close()

Libération de l'espace mémoire réservé pour le maillage. À appeler avant d'ouvrir une nouvelle instance de la classe.

Pour un problème en coordonnées cylindriques et à symétrie axiale, on utilise la classe PoissonAxial.

object = Poisson.PoissonAxial(pw_min,ph_min,levels,width,height,x0=0,y0=0)

Constructeur : définition d'un maillage 2D

pw_min, ph_min : puissances définissant les dimensions du maillage réduit. Les nombres de mailles sont 2^pw_min dans la direction X et 2^ph_min dans la direction Y.
levels : nombre de niveaux de grilles. Le maillage le plus fin est défini par les puissances pw_min+levels-1 et ph_min+levels-1.
width : largeur du domaine.
height : hauteur du domaine. L'intervalle de r est [0,height].
x0, y0 : coordonnées du point origine.

La première dimension du maillage (axe Ox en coordonnées cartésiennes) correspond à l'axe de symétrie Oz. Les indices des points sur cet axe sont donc (i,0).

Les fonctions de la classe Poisson décrites ci-dessous sont disponibles aussi pour la classe PoissonAxial.

Pour l'équation de Poisson vectorielle en coordonnées cylindriques et à symétrie axiale, on utilise la classe PoissonAxialVect.

object = Poisson.PoissonAxialVect(pw_min,ph_min,levels,width,height,x0=0,y0=0)

Constructeur : définition d'un maillage 2D

pw_min, ph_min : puissances définissant les dimensions du maillage réduit. Les nombres de mailles sont 2^pw_min dans la direction X et 2^ph_min dans la direction Y.
levels : nombre de niveaux de grilles. Le maillage le plus fin est défini par les puissances pw_min+levels-1 et ph_min+levels-1.
width : largeur du domaine.
height : hauteur du domaine. L'intervalle de r est [0,height].
x0, y0 : coordonnées du point origine.

4.b. Équation de Poisson

La fonction suivante permet d'obtenir la discrétisation du laplacien sur tout le domaine :

Poisson.laplacien()

Discrétisation du laplacien

Dans le cas de la géométrie axiale, le laplacien est aussi discrétisé sur l'axe de symétrie Oz.

Une source peut être ajoutée en un point du maillage au moyen de la fonction suivante :

Poisson.source(i,j,s)

Source ponctuelle

i,j : indices du point sur le maillage réduit
s : valeur de s(x,y)

On peut aussi ajouter une source étendue sous forme d'un bloc rectangulaire :

Poisson.source_rect(ci,cj,width,height,s)

Source rectangulaire

ci,cj : indices du centre sur le maillage réduit.
width,height : demi-largeur et demi-hauteur du rectangle sur le maillage réduit.
s : valeur de s(x,y) dans le rectangle.

La fonction suivante permet de définir une source sous forme d'un polygone (ouvert ou fermé). Si le polygone est fermé, seul le bord est affecté.

Poisson.source_polygon(istart,jstart,lilist,ljlist,value)

Définit un polygone avec une source.

istart,jstart : indices du point de départ sur le maillage réduit.
lilist : liste des déplacements parallèles à X.
ljlist : liste des déplacements parallèles à Y.
value : valeur à imposer sur le polygone.

4.c. Conditions limites sur les bords

Les deux fonctions suivantes permettent de définir les conditions limites sur les bords du domaine rectangulaire.

Poisson.dirichlet_borders(value)

Condition limite de Dirichlet sur les bords du domaine rectangulaire.

value : valeur sur les bords de u

Dans le cas de la géométrie axiale, la condition ne s'applique pas sur l'axe de symétrie Oz.

Poisson.neumann_borders(source,derivX1,derivX2,derivY1,derivY2)

Condition limite de Neumann sur les bords du domaine, c'est-à-dire valeur de la dérivée selon la normale de chaque bord.

source : valeur de s(x,y) au voisinage du bord (0 pour l'équation de Laplace)
derivX1 : valeur de la dérivée par rapport à x à imposer sur le bord gauche
derivX2 : valeur de la dérivée par rapport à x à imposer sur le bord droit
derivY1 : valeur de la dérivée par rapport à y à imposer sur le bord inférieur
derivY2 : valeur de la dérivée par rapport à y à imposer sur le bord supérieur

Dans le cas de la géométrie axiale, la condition ne s'applique pas sur l'axe de symétrie Oz (bord inférieur).

Pour combiner des conditions de Neumann et Dirichlet sur les bords du domaine, il faut définir explicitement les bords avec un polygone, comme expliqué ci-dessous.

4.d. Polygone

On peut définir un polygone dont les côtés sont parallèles aux axes x et y (ou z et r en géométrie axiale).

Poisson.dirichlet_polygon(istart,jstart,lilist,ljlist,value)

Définit un polygone avec une condition limite de dirichlet.

istart,jstart : indices du point de départ sur le maillage réduit.
lilist : liste des déplacements parallèles à X.
ljlist : liste des déplacements parallèles à Y.
value : valeur à imposer sur le polygone.

La figure suivante montre le polygone ouvert obtenu par dirichlet_polygon(1,2,[2,0],[0,1],1).

Figure pleine page

Pour obtenir un polygone fermé, il faut se déplacer jusqu'au premier point. On peut aussi définir un polygone sur les bords du domaine, pour modifier la condition limite ou même pour la définir entièrement. Tous les points des grilles jusqu'à la grille la plus fine (ici 32x32) sont bien sûr affectés par cette condition limite.

Poisson.neumann_polygon(istart,jstart,lilist,ljlist,derivXlist,derivYlist,source)

Définit un polygone fermé avec une condition limite de Neumann.

istart,jstart : indices du point de départ sur le maillage réduit.
lilist : liste des déplacements parallèles à X.
ljlist : liste des déplacements parallèles à Y.
derivXlist : liste des dérivées par rapport à X.
derivYlist : liste des dérivées par rapport à Y.
source : source à appliquer sur le polygone.

Les dérivées par rapport à X n'affectent que les côtés parallèles à Y. En principe, un polygone avec une condition de Neumann doit être fermé : il faut donc se déplacer jusqu'au point initial. Les angles recoivent un traitement spécial, qui combine les dérivées selon X et selon Y. Il est important de savoir où se trouve le domaine de calcul par rapport au polygone. Celui-ci est défini par rapport au sens du parcours : si le sens de parcours est le sens horaire alors le domaine se trouve à l'extérieur du polygone.

Poisson.interface_polygon(istart,jstart,lilist,ljlist,source,a1,a2)

Définit un polygone qui délimite deux milieux différents (conductivités thermiques ou permittivités électriques différentes).

istart,jstart : indices du point de départ sur le maillage réduit.
lilist : liste des déplacements parallèles à X.
ljlist : liste des déplacements parallèles à Y.
source : source à appliquer sur le polygone.
a1,a2 : conductivités à l'intérieur et à l'extérieur (pour le sens horaire).

Ce polygone peut être ouvert. Lorsqu'on parcourt un côté, le milieu de conductivité a1 se trouve à gauche.

4.e. Itérations de Gauss-Seidel sur CPU

Le système linéaire obtenu par discrétisation des équations (équation de Poisson et conditions limites) est résolu par la méthode d'itération de Gauss-Seidel sur le maillage le plus fin. Ce paragraphe présente les fonctions de la classe Poisson qui permettent d'effectuer ce calcul sur le processeur central (CPU).

Poisson.iterations(niter,threads=1,omega=1.0)

Itérations de Gauss-Seidel.

niter : nombre d'itérations
threads : nombre de threads utilisés dans le CPU (1, 2 ou 4)
omega : paramètre de sur-relaxation

Pour contrôler la convergence, la fonctions suivante effectue des blocs d'itérations et calcule la norme de la matrice U à chaque bloc.

[ni,norm]=Poisson.iterations_norm(niter,nblock,threads=1,omega=1.0)

Itérations avec calcul de la norme de la matrice des valeurs de u.

niter : nombre d'itérations dans un bloc
nblock : nombre de blocs
threads : nombre de threads utilisés dans le CPU (1, 2 ou 4)
omega : paramètre de sur-relaxation

ni : liste des nombres d'itérations
norm : liste des normes

Sur les processeurs comportant au moins 4 cœurs, l'exécution avec 4 threads est plus rapide.

La fonction suivante fournit le paramètre de sur-relaxation optimale pour le problème de Poisson avec des conditions limites de Dirichlet sur les bords.

omega=Poisson.get_optimal_omega()

Renvoit le paramètre optimal, qui dépend de la taille du maillage.

omega : paramètre de sur-relaxation optimal

Pour les autres types de conditions limites, on peut partir de cette valeur mais la valeur optimale peut être notablement différente.

4.f. Itérations de Gauss-Seidel sur GPU

Ces fonctions ne sont disponibles que sur la version avec OpenCL.

Les itérations peuvent être effectuées sur processeur graphique (GPU). Une plateforme openCL doit être installée.

La première fonction permet d'afficher les plateformes openCL présentes sur le système, et pour chaque plateforme les périphériques associés. Les plateformes et les périphériques sont numérotés à partir de 0.

Si par exemple on dispose d'une carte graphique NVIDIA et d'un contrôleur graphique Intel, il est probable que la première soit accessible par la plateforme 0 et la seconde par la platforme 1 (avec un seul périphérique par plateforme).

Poisson.opencl_platforms()

Affiche sur la console les plateformes openCL et leurs périphériques associés.

Par défaut, la plateforme 0 et le périphérique 0 sont sélectionnés. La fonction suivante permet de sélectionner une plateforme openCL et un périphérique :

Poisson.set_opencl_platform_device(platform,device)

Sélection d'une plateforme et d'un périphérique pour effectuer les itérations.

platform : numéro de la plateforme, 0 pour la première
device : numéro du périphérique, 0 pour le premier

Les deux fonctions suivantes effectuent les itérations, et sont analogues aux fonctions définies plus haut pour le CPU :

Poisson.opencl_iterations(niter,omega=1.0)

Itérations de Gauss-Seidel sur plateforme openCL.

niter : nombre d'itérations
omega : paramètre de sur-relaxation, compris entre 1 et 2

[ni,norm]=Poisson.opencl_iterations_norm(niter,nblock,omega=1.0)

Itérations avec calcul de la norme de la matrice des valeurs de u.

niter : nombre d'itérations dans un bloc
nblock : nombre de blocs
omega : paramètre de sur-relaxation, compris entre 1 et 2

ni : liste des nombres d'itérations
norm : liste des normes

4.g. Méthodes multigrilles

Les méthodes multigrilles sont implémentées pour le processeur central (CPU).

[ni,norm]=Poisson.multigrid_Vcycles_norm(nsteps,nlevels,niter,nu1,nu2,nblocks,tolerance=1.0e-8)

Effectue des cycles en V.

nsteps : nombre de cycles par bloc
nlevels : nombre de niveaux de grilles utilisés (au maximum le nombre défini dans le constructeur)
niter : nombre d'itérations sur la grille la plus grosse, pour obtenir la solution sur cette grille
nu1 : nombre d'itérations pour les phases descendantes (restrictions)
nu2 : nombre d'itérations pour les phases montantes (prolongations)
nblocks : nombre de blocs maximal. Pour chaque bloc, la norme de U est calculée.
tolerance : tolérance pour la norme du résidu. Si celle-ci est inférieure à cette tolérance, le calcul est stoppé.

ni : liste des nombres d'itérations
norm : liste des normes

Au sujet de ces paramètres, voir Méthode itérative multigrilles

Les nombres d'itérations sont les nombres équivalents pour la grille fine, c'est-à-dire les nombres calculés en tenant compte de la réduction du nombre de points des grilles (4,16, etc). On peut donc comparer la courbe obtenue avec celle donnée par la méthode de Gauss-Seidel. Néanmoins, les transferts intergrilles (prolongation et restriction) ne sont pas pris en compte dans le calcul. On peut estimer que la vitesse de convergence donnée par cette courbe est sur-estimée d'environ 10 à 20 pour cents.

[ni,norm]=Poisson.multigrid_full_norm(mcycle,nlevels,niter,nu1,nu2)

Effectue un cycle multigrille complet

mcycles : nombre de cycles en V pour la correction à chaque niveau de grille
nlevels : nombre de niveaux de grilles utilisés (au maximum le nombre défini dans le constructeur)
niter : nombre d'itérations sur la grille la plus grosse, pour obtenir la solution sur cette grille
nu1 : nombre d'itérations pour les phases descendantes (restrictions)
nu2 : nombre d'itérations pour les phases montantes (prolongations)

ni : liste des nombres d'itérations
norm : liste des normes

Pour le multigrille complet, la liste des normes ne contient que deux valeurs, la valeur initiale et la valeur finale.

Si les paramètres sont bien choisis, une méthode multigrille nécessite un nombre d'itérations proportionnel au nombre de points du maillage, ce qui en fait une méthode optimale.

4.h. Récupération des données

La matrice U contient les valeurs de u(x,y) aux points de la maille. Elle est fournie sous forme d'un tableau de type numpy.ndarray par la fonction suivante :

U=Poisson.get_array(symetry=0)

symetry : si égal à 1, un tableau symétrique par rapport à l'axe Ox est renvoyé, un tableau antisymétrique pour -1.

U : tableau numpy contenant les valeurs de u

Les dérivées par rapport à x et y sont obtenues avec les deux fonctions suivantes :

DX=Poisson.get_derivX(symetry=0)

symetry : si égal à 1, un tableau symétrique par rapport à l'axe Ox est renvoyé, un tableau antisymétrique pour -1.

DX : tableau numpy contenant les valeurs de $\frac{\partial u}{\partial x}$

DY=Poisson.get_derivY(symetry=0)

symetry : si égal à 1, un tableau symétrique par rapport à l'axe Ox est renvoyé, un tableau antisymétrique pour -1.

DY : tableau numpy contenant les valeurs de $\frac{\partial u}{\partial y}$

En géométrie axiale, les fonction PoissonAxial.get_derivZ, PoissonAxial.get_derivR, PoissonAxialVect.get_derivZ et PoissonAxialVect.get_derivZ sont analogues au précédentes. Dans ce cas, il peut être nécessaire de calculer la dérivée suivante (par exemple pour calculer le rotationnel du potentiel vecteur en magnétostatique) :

\frac{1}{r} \frac{\partial (r u)}{\partial r}

La fonction suivante effectue ce calcul :

DY=PoissonAxialVect.get_derivRUR(symetry=0)

symetry : si égal à 1, un tableau symétrique par rapport à l'axe Ox est renvoyé, un tableau antisymétrique pour -1.

DY : tableau numpy contenant les valeurs de $\frac{1}{r} \frac{\partial (r u)}{\partial r}$

Pour faciliter le tracé des résultats, il peut être utile de calculer les valeurs des abscisses et ordonnées des points du maillage. Les fonctions suivantes permettent d'obtenir ces listes. Le point origine du domaine est celui spécifié par les arguments x0 et y0 du constructeur.

x = Poisson.get_x()

x : tableau numpy contenant les valeurs de x

y = Poisson.get_y()

y : tableau numpy contenant les valeurs de y

z = PoissonAxial.get_z()

z : tableau numpy contenant les valeurs de z

r = PoissonAxial.get_R(symetry=0)

r : tableau numpy contenant les valeurs de r
symetry : si non nul, renvoit un tableau symétrique comportant aussi des valeurs négatives de r

z = PoissonAxialVect.get_z()

z : tableau numpy contenant les valeurs de z

r = PoissonAxialVect.get_R(symetry=0)

r : tableau numpy contenant les valeurs de r
symetry : si non nul, renvoit un tableau symétrique comportant aussi des valeurs négatives de r

extent = Poisson.get_extent()

Renvoit les bornes du domaines

extent : (xmin,xmax,ymin,ymax)

mask=Poisson.get_mask()

Renvoit une image comportant un tracé des objets définis sur le maillage.

mask : image ayant la taille du maillage, à utiliser avec matplotlib.pyplot.imshow.

4.i. Initialisation de U

Poisson.init_array()

Initialise à zéro la matrice des valeurs de u

Poisson.sinus_mode(mx,my)

Ajoute un mode sinusoïdal à la matrice des valeurs de u (pour les tests de convergences)

mx,my : indices des modes