Équation de diffusion à deux dimensions

4.a. Définition

Le schéma implicite à directions alternées consiste à décomposer le pas de temps (n,n+1) en deux sous-pas. Pour le premier, on écrit un schéma implicite pour la dérivée par rapport à x mais explicite pour la dérivée par rapport à y :

$$\frac{U_{i,j}^{n+1/2}-U_{i,j}^n}{(\Delta t/2)}=D\frac{U_{i+1,j}^{n+1/2}-2U_{i,j}^{n+1/2}+U_{i-1,j}^{n+1/2}}{{\left(\Delta x\right)}^2}+D\frac{U_{i,j+1}^n-2U_{i,j}^n+U_{i,j-1}^n}{{\left(\Delta y\right)}^2}+S_{i,j}^n\text{(21)}$$

Le second est implicite dans la direction y :

$$\frac{U_{i,j}^{n+1}-U_{i,j}^{n+1/2}}{(\Delta t/2)}=D\frac{U_{i+1,j}^{n+1/2}-2U_{i,j}^{n+1/2}+U_{i-1,j}^{n+1/2}}{{\left(\Delta x\right)}^2}+D\frac{U_{i,j+1}^{n+1}-2U_{i,j}^{n+1}+U_{i,j-1}^{n+1}}{{\left(\Delta y\right)}^2}+S_{i,j}^n\text{(22)}$$

Chacune de ces deux étapes consiste à résoudre un système tridiagonal, ce qui se fait par élimination directe en utilisant l'algorithme décrit dans équation de diffusion à une dimension.

On écrira les deux systèmes à résoudre sous la forme générale suivante :

$$\begin{array}{ccc}\hfill & A_{i,j}{U}_{i-1,j}^{n+1/2}+B_{i,j}{U}_{i,j}^{n+1/2}+C_{i,j}{U}_{i+1,j}^{n+1/2}=D_{i,j}{U}_{i,j-1}^n+E_{i,j}{U}_{i,j}^n+F_{i,j}{U}_{i,j+1}^n+G_{i,j}\hfill & \hfill \text{(23)}\\ \hfill & A{\text{'}}_{i,j}{U}_{i,j-1}^{n+1}+B{\text{'}}_{i,j}{U}_{i,j}^{n+1}+C{\text{'}}_{i,j}{U}_{i,j+1}^{n+1}=D{\text{'}}_{i,j}{U}_{i-1,j}^{n+1/2}+E{\text{'}}_{i,j}{U}_{i,j}^{n+1/2}+F{\text{'}}_{i,j}{U}_{i+1,j}^{n+1/2}+G{\text{'}}_{i,j}\hfill & \hfill \text{(24)}\end{array}$$

Le premier système tridiagonal est à résoudre pour chaque ligne. Le second est à résoudre pour chaque colonne.

Pour le point A à l'intérieur du domaine, on a donc pour le premier système :

$$\begin{array}{ccc}\hfill & A_{i,j}=\frac{D}{\Delta x^2}\hfill & \hfill \text{(25)}\\ \hfill & B_{i,j}=-\frac{2D}{\Delta x^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(26)}\\ \hfill & C_{i,j}=\frac{D}{\Delta x^2}\hfill & \hfill \text{(27)}\\ \hfill & D_{i,j}=-\frac{D}{\Delta y^2}\hfill & \hfill \text{(28)}\\ \hfill & E_{i,j}=\frac{2D}{\Delta y^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(29)}\\ \hfill & F_{i,j}=-\frac{D}{\Delta y^2}\hfill & \hfill \text{(30)}\\ \hfill & G_{i,j}=-S_{i,j}^n\hfill & \hfill \text{(31)}\end{array}$$

et pour le second :

$$\begin{array}{ccc}\hfill & A{\text{'}}_{i,j}=\frac{D}{\Delta y^2}\hfill & \hfill \text{(32)}\\ \hfill & B{\text{'}}_{i,j}=-\frac{2D}{\Delta y^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(33)}\\ \hfill & C{\text{'}}_{i,j}=\frac{D}{\Delta y^2}\hfill & \hfill \text{(34)}\\ \hfill & D{\text{'}}_{i,j}=-\frac{D}{\Delta x^2}\hfill & \hfill \text{(35)}\\ \hfill & E{\text{'}}_{i,j}=\frac{2D}{\Delta x^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(36)}\\ \hfill & F{\text{'}}_{i,j}=-\frac{D}{\Delta x^2}\hfill & \hfill \text{(37)}\\ \hfill & G{\text{'}}_{i,j}=-S_{i,j}^{n+1/2}\hfill & \hfill \text{(38)}\end{array}$$

En pratique, on multipliera tous ces coefficients par Δt et on simplifiera en posant :

$$\begin{array}{ccc}\hfill & {\alpha }_x=\frac{D\Delta t}{\Delta x^2}\hfill & \hfill \text{(39)}\\ \hfill & {\alpha }_y=\frac{D\Delta t}{\Delta y^2}\hfill & \hfill \text{(40)}\end{array}$$

Dans la plupart des cas, on prendra Δx=Δy.

4.b. Conditions limites

Une condition limite de Dirichlet sur le point (i,j) s'écrit :

$$\begin{array}{ccc}\hfill & A_{i,j}=0\hfill & \hfill \text{(41)}\\ \hfill & B_{i,j}=1\hfill & \hfill \text{(42)}\\ \hfill & C_{i,j}=0\hfill & \hfill \text{(43)}\\ \hfill & D_{i,j}=0\hfill & \hfill \text{(44)}\\ \hfill & E_{i,j}=0\hfill & \hfill \text{(45)}\\ \hfill & F_{i,j}=0\hfill & \hfill \text{(46)}\\ \hfill & G_{i,j}=U_0\hfill & \hfill \text{(47)}\end{array}$$

Considérons une condition limite de Neumann en un point du bord gauche (point D). Cette condition consiste à fixer la dérivée suivante :

$$\frac{\partial u}{\partial x}=w_x\text{(48)}$$

Pour discrétiser cette condition limite, on suit le même principe que pour l'équation de diffusion à une dimension. On commence par écrire une discrétisation d'ordre 2 en faisant intervenir un point fictif situé à gauche de la frontière :

$$\frac{U_{1,j}^n-U_{-1,j}^n}{2\Delta x}=w_x\text{(49)}$$

La discrétisation de la dérivée seconde par rapport à x s'écrit :

$$\frac{{\partial }^{2}u}{\partial x^2}\to \frac{U_{1,j}^n-2U_{0,j}^n+U_{-1,j}^n}{{\left(\Delta x\right)}^2}\text{(50)}$$

Le point fictif i=-1 est éliminé avec la relation $\mathrm{(49)}$ :

$$\frac{{\partial }^{2}u}{\partial x^2}\to \frac{2U_{1,j}^n-2U_{0,j}^n-2\Delta x{w}_x}{{\left(\Delta x\right)}^2}\text{(51)}$$

Pour le schéma implicite à directions alternées, le premier demi-pas, correspondant à la direction x, est défini par :

$$\begin{array}{ccc}\hfill & A_{i,j}=0\hfill & \hfill \text{(52)}\\ \hfill & B_{i,j}=-\frac{2D}{\Delta x^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(53)}\\ \hfill & C_{i,j}=\frac{2D}{\Delta x^2}\hfill & \hfill \text{(54)}\\ \hfill & D_{i,j}=-\frac{D}{\Delta y^2}\hfill & \hfill \text{(55)}\\ \hfill & E_{i,j}=\frac{2D}{\Delta y^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(56)}\\ \hfill & F_{i,j}=-\frac{D}{\Delta y^2}\hfill & \hfill \text{(57)}\\ \hfill & G_{i,j}=\frac{2D{w}_x}{\Delta x}-S_{i,j}^n\hfill & \hfill \text{(58)}\end{array}$$

Le seconde demi-pas, correspondant à la direction y, est défini par :

$$\begin{array}{ccc}\hfill & A{\text{'}}_{i,j}=\frac{D}{\Delta y^2}\hfill & \hfill \text{(59)}\\ \hfill & B{\text{'}}_{i,j}=-\frac{2D}{\Delta y^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(60)}\\ \hfill & C{\text{'}}_{i,j}=\frac{D}{\Delta y^2}\hfill & \hfill \text{(61)}\\ \hfill & D{\text{'}}_{i,j}=0\hfill & \hfill \text{(62)}\\ \hfill & E{\text{'}}_{i,j}=\frac{2D}{\Delta x^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(63)}\\ \hfill & F{\text{'}}_{i,j}=-\frac{2D}{\Delta x^2}\hfill & \hfill \text{(64)}\\ \hfill & G{\text{'}}_{i,j}=\frac{2D{w}_x}{\Delta x}-S_{i,j}^{n+1/2}\hfill & \hfill \text{(65)}\end{array}$$

Voyons comment retrouver ce résultat avec la méthode des volumes finis, qui sera plus simple pour traiter les cas plus complexes. Le contour (o,p,q,r) est défini sur la figure. Voici comment sont discrétisées les différentes intégrales :

$$\begin{array}{ccc}\hfill & {\iint }_{\Sigma }\frac{\partial u}{\partial t}d\Sigma \to \frac{U_{i,j}^{n+1}-U_{i,j}^n}{\Delta t}\frac{1}{2}\Delta x\Delta y\hfill & \hfill \text{(66)}\\ \hfill & {\int }_{op}D\left(\overrightarrow{grad}u\right)\cdot \overrightarrow{n}dl\to D\frac{U_{i,j+1}^n-U_{i,j}^n}{\Delta y}\frac{1}{2}\Delta x\hfill & \hfill \text{(67)}\\ \hfill & {\int }_{pq}D\left(\overrightarrow{grad}u\right)\cdot \overrightarrow{n}dl\to -D{w}_x\Delta y\hfill & \hfill \text{(68)}\\ \hfill & {\int }_{qr}D\left(\overrightarrow{grad}u\right)\cdot \overrightarrow{n}dl\to D\frac{U_{i,j-1}^n-U_{i,j}^n}{\Delta y}\frac{1}{2}\Delta x\hfill & \hfill \text{(69)}\\ \hfill & {\int }_{ro}D\left(\overrightarrow{grad}u\right)\cdot \overrightarrow{n}dl\to D\frac{U_{i+1,j}^n-U_{i,j}^n}{\Delta x}\Delta y\hfill & \hfill \text{(70)}\\ \hfill & {\iint }_{\Sigma }sd\Sigma \to S_{i,j}^n\frac{1}{2}\Delta x\Delta y\hfill & \hfill \text{(71)}\end{array}$$

En écrivant le schéma implicite à directions alternées à partir de ces équations, on retouve bien les coefficients précédents (attention à la division de Δt par deux dans cette opération).

Pour les bords du domaine rectangulaire, il y a 8 types de conditions limites de Neumann (4 côtés et 4 coins). Pour un objet situé à l'intérieur du domaine (contour GHIJKL), il a 4 autres types de coins. Il y a donc en tout 12 types de condition limite de Neumann.

Pour une condition limite de Neumann, il faut une convention pour fixer l'intérieur du domaine de calcul par rapport au sens de parcours d'un contour. On adopte la convention suivante : l'intérieur du domaine est situé à gauche si l'on s'oriente dans le sens du parcours. Par exemple, le contour rectangulaire délimitant le domaine doit être parcouru dans le sens trigonométrique. Le contour de l'objet situé à l'intérieur du domaine devra être parcouru dans le sens horaire. Compte tenu de cette convention, on peut donner une nomenclature simple pour désigner les différentes arêtes et coins. Une arête ou un coin sont définis par deux segments, chacun d'eux défini par deux incréments di et dj égaux à -1, 0 ou +1. Par exemple, le coin inférieur gauche du domaine est défini par (0,-1)-(1,0). Une arête situé sur le bord inférieur est définie par (1,0)-(1,0). Le coin G de l'onjet situé dans le domaine est défini par (-1,0)-(0,1).

À titre d'exemple, voici le résultat pour le coin inférieur gauche :

$$\begin{array}{ccc}\hfill & A_{i,j}=0\hfill & \hfill \text{(72)}\\ \hfill & B_{i,j}=-\frac{4D}{\Delta x^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(73)}\\ \hfill & C_{i,j}=\frac{4D}{\Delta x^2}\hfill & \hfill \text{(74)}\\ \hfill & D_{i,j}=0\hfill & \hfill \text{(75)}\\ \hfill & E_{i,j}=\frac{2D}{\Delta y^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(76)}\\ \hfill & F_{i,j}=-\frac{2D}{\Delta y^2}\hfill & \hfill \text{(77)}\\ \hfill & G_{i,j}=-S_{i,j}^n+\frac{2D{w}_x}{\Delta x}+\frac{2D{w}_y}{\Delta y}\hfill & \hfill \text{(78)}\\ \hfill & A{\text{'}}_{i,j}=0\hfill & \hfill \text{(79)}\\ \hfill & B{\text{'}}_{i,j}=-\frac{2D}{\Delta y^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(80)}\\ \hfill & C{\text{'}}_{i,j}=\frac{2D}{\Delta y^2}\hfill & \hfill \text{(81)}\\ \hfill & D{\text{'}}_{i,j}=0\hfill & \hfill \text{(82)}\\ \hfill & E{\text{'}}_{i,j}=\frac{2D}{\Delta x^2}-\frac{2}{\Delta t}\hfill & \hfill \text{(83)}\\ \hfill & F{\text{'}}_{i,j}=-\frac{4D}{\Delta x^2}\hfill & \hfill \text{(84)}\\ \hfill & G{\text{'}}_{i,j}=-S_{i,j}^{n+1/2}+\frac{2D{w}_x}{\Delta x}+\frac{2D{w}_y}{\Delta y}\hfill & \hfill \text{(85)}\end{array}$$

4.c. Algorithme

La maillage est constitué de matrices à N_x colonnes et N_y lignes. La définition des pas de temps et d'espace, des conditions limites et des sources permet de construire les matrices A,B,C,D,E,F,G,A',B',C',D',E',F',G'.

Il faut aussi construire les matrices β,γ pour la résolution des systèmes sur les lignes, et les matrices β',γ' pour la résolution des systèmes sur les colonnes. Pour chaque ligne j, la récurrence à appliquer pour i variant de 1 à N_x :

$$\begin{array}{ccc}\hfill & {\beta }_{1,j}=B_{1,j}\hfill & \hfill \text{(86)}\\ \hfill & {\gamma }_{1,j}=\frac{C_{1,j}}{{\beta }_{1,j}}\hfill & \hfill \text{(87)}\\ \hfill & {\beta }_{i,j}=B_{i,j}-A_{i,j}{\gamma }_{i-1,j}\hfill & \hfill \text{(88)}\\ \hfill & {\gamma }_{i,j}=\frac{C_{i,j}}{{\beta }_{i,j}}\hfill & \hfill \text{(89)}\end{array}$$

Pour chaque colonne i, la récurrence à appliquer pour j variant de 1 à N_y :

$$\begin{array}{ccc}\hfill & \beta {\text{'}}_{i,1}=B{\text{'}}_{i,1}\hfill & \hfill \text{(90)}\\ \hfill & \gamma {\text{'}}_{i,1}=\frac{C{\text{'}}_{i,1}}{\beta {\text{'}}_{i,1}}\hfill & \hfill \text{(91)}\\ \hfill & \beta {\text{'}}_{i,j}=B{\text{'}}_{i,j}-A{\text{'}}_{i,j}\gamma {\text{'}}_{i,j-1}\hfill & \hfill \text{(92)}\\ \hfill & \gamma {\text{'}}_{i,j}=\frac{C{\text{'}}_{i,j}}{\beta {\text{'}}_{i,j}}\hfill & \hfill \text{(93)}\end{array}$$

Pour chaque intervalle de temps Δt, on doit effectuer une résolution tridiagonale pour chaque ligne, puis pour chaque colonne.

Pour chaque ligne j on commence par calculer la suite suivante pour i variant de 1 à N_x:

$$\begin{array}{ccc}\hfill & x_1=\frac{1}{{\beta }_{1,j}}(D_{1,j}{U}_{1,j-1}^n+E_{1,j}{U}_{1,j}^n+F_{1,j}{U}_{1,j+1}^n+G_{1,j})\hfill & \hfill \text{(94)}\\ \hfill & x_i=\frac{1}{{\beta }_{i,j}}(D_{i,j}{U}_{i,j-1}^n+E_{i,j}{U}_{i,j}^n+F_{i,j}{U}_{i,j+1}^n+G_{i,j}-A_{i,j}{x}_{i-1})\hfill & \hfill \text{(95)}\end{array}$$

Pour la même ligne j, on calcule les nouvelles valeurs à l'instant n+1/2 par récurrence, en faisant varier de i de N_x à 1 :

$$\begin{array}{ccc}\hfill & U_{N_x,j}^{n+\frac{1}{2}}=x_{N_x}\hfill & \hfill \text{(96)}\\ \hfill & U_{i,j}^{n+\frac{1}{2}}=x_i-{\gamma }_{i,j}{U}_{i+1,j}^{n+\frac{1}{2}}\hfill & \hfill \text{(97)}\end{array}$$

Les valeurs à l'instant n+1/2 sont stockées dans la matrice U. On procède alors au parcours des colonnes. Pour chaque colonne i, on calcule la suite suivante pour j variant de 1 à N_y :

$$\begin{array}{ccc}\hfill & x_1=\frac{1}{\beta {\text{'}}_{i,1}}(D{\text{'}}_{i,1}{U}_{i-1,1}^{n+1/2}+E{\text{'}}_{i,1}{U}_{i,1}^{n+1/2}+F{\text{'}}_{i,1}{U}_{i+1,1}^{n+1/2}+G{\text{'}}_{i,1})\hfill & \hfill \text{(98)}\\ \hfill & x_i=\frac{1}{\beta {\text{'}}_{i,j}}(D{\text{'}}_{i,j}{U}_{i-1,j}^{n+1/2}+E{\text{'}}_{i,j}{U}_{i,j}^{n+1/2}+F{\text{'}}_{i,j}{U}_{i+1,j}^{n+1/2}+G{\text{'}}_{i,j}-A{\text{'}}_{i,j}{x}_{i-1})\hfill & \hfill \text{(99)}\end{array}$$

puis pour j variant de N_y à 1 :

$$\begin{array}{ccc}\hfill & U_{i,N_y}^n=x_{N_y}\hfill & \hfill \text{(100)}\\ \hfill & U_{i,j}^n=x_i-\gamma {\text{'}}_{i,j}{U}_{i,j+1}^n\hfill & \hfill \text{(101)}\end{array}$$

Lorsqu'on traite les lignes, il faut écrire dans une seconde matrice U, pour ne pas écraser les valeurs de l'instant n. Il en est de même lorsqu'on traite les colonnes.

Le nombre de matrices utilisées est 20. Pour un maillage 1024x1024 en flottants 32 bits, la mémoire utilisée est de 80 méga-octets.

Les lignes peuvent être parcourues en parallèle (idem pour les colonnes). On peut donc avec profit utiliser un processeur graphique pour faire ces calculs (via l'API CUDA ou OpenCL). Pour le traitement des lignes, chaque unité de calcul traite une ligne. Les opérations effectuées sur les différentes lignes sont rigoureusement identiques, ce qui est le cas idéal pour une implémentation sur GPU. Bien entendu, toutes les matrices doivent être transférées sur la mémoire de la carte graphique avant le calcul.

On remarque que le traitement des bords (qui doit être identique au traitement des autres points) nécessite un accès à des éléments de U hors domaine. Pour résoudre ce problème, on utilise des matrices avec des rangées fictives sur les bords. Le nombre d'éléments dans chaque direction est une puissance de 2 (plus efficace pour l'implémentation GPU). La figure suivante montre un exemple pour N_x=2³ et N_y=2³ (les points sont les nœuds de la grille).

Figure pleine page

Sur cet exemple, les matrices sont 10x10 mais le maillage effectif est 8x8.