Version MPI du code de résolution numérique de l'équation de chaleur 2D

DOURNAC.ORG
	Français English

Informatique > Version MPI du code de résolution numérique de l'équation de chaleur 2D

1.Exécution :

Les codes séquentiels et parallélisés se trouvent dans les archives suivantes :

Version Langage C : Heat2D_MPI_Solving_C.tar.gz

Version Fortran 90 : Heat2D_MPI_Solving_F90.tar.gz

Veuillez vous assurer à l'exécution du code parallèle (avec la commande "mpirun") que le nombre de processus soit égal au produit du nombre de sous-domaines selon x par le nombre de sous-domaines selon y ( paramètres "x_domains" et "y_domains" situés dans le fichier "param" ) :

$ cat param
  size_x
1024
  size_y
1024
  x_domains
16
  y_domains
16
  MaxStep
100000000
  dt
1.0e-1
  cnv_tol
1.0e-1

Nous allons donc utiliser dans ce cas 256 processus. La ligne de commande est alors :

$ mpirun -np 256 ./explicitPar

  Time step too large in 'param' file - Taking convergence criterion

  Time step = 2.374899779E-07

  Convergence = 0.100000000 after 40595 steps

  Problem size = 1048576

  Wall Clock = 487.038250208

  Computed solution in outputPar.dat file

2.Équation de chaleur instationnaire 2D :

La formule générale de l'équation de Fourier est :

\begin{equation} \dfrac{\partial T}{\partial t}=\kappa\,\Delta T\,\,\,\,\,\,\text{avec}\,\,\,\,\,\Delta=\sum_{i=1}^{n}\dfrac{\partial^{2}}{\partial x_{i}^{2}}\,\,\,\,\text{le Laplacien en dimension $n$} \label{eq1} \end{equation}

Le coefficient $\kappa$ est la conductivité thermique. L'équation de chaleur à 2 dimensions est donc :

\begin{equation} \dfrac{\partial\theta}{\partial t}=\kappa\,\bigg(\dfrac{\partial^{2}\theta}{\partial x^{2}}+\dfrac{\partial^{2}\theta}{\partial y^{2}}\bigg) \label{eq2} \end{equation}

3.Discrétisation :

La discrétisation pour une résolution numérique est de 2 types : spatiale et temporelle

3.1 Discrétisation spatiale :

Nous allons commencer par discrétiser le membre de droite de l'équation ci-dessus, c'est-à-dire dans le domaine spatial. Pour cela, nous écrivons les développements de Taylor suivants :

\begin{equation} \theta_{m+1}=\theta(x_{m}+h)=\theta(x_{m})+h\,\theta'(x_{m})+\dfrac{h^{2}}{2}\,\theta''(x_{m})+\dfrac{h^{3}}{6}\,\theta^{'''}(x_{m})+O(h^{4}) \label{eq3} \end{equation} \begin{equation} \theta_{m-1}=\theta(x_{m}-h)=\theta(x_{m})-h\,\theta'(x_{m})+\dfrac{h^{2}}{2}\,\theta''(x_{m})-\dfrac{h^{3}}{6}\,\theta^{'''}(x_{m})+O(h^{4}) \label{eq4} \end{equation}

En additionnant les relations \eqref{eq3} et \eqref{eq4}, nous obtenons un développement de la dérivée seconde au point $x_{m}$ :

\begin{equation} \theta''(x_{m})=\dfrac{\theta_{m+1}-2\theta_{m}+\theta_{m-1}}{h^{2}}+O(h^{2}) \label{eq5} \end{equation}

On obtient les mêmes relations pour la variable $y$. En prenant $h_{x}=size_{x}/N_{x}$ et $h_{y}=size_{y}/N_{y}$ et avec la convention $\theta(x_{m},y_{m})=\theta[i][j]$, on obtient :

\begin{equation} \dfrac{\partial^{2}\theta}{\partial x^{2}}+\dfrac{\partial^{2}\theta}{\partial y^{2}}=\dfrac{\theta[i+1][j]-2\theta[i][j]+\theta[i-1][j]}{h_{x}^{2}}+\dfrac{\theta[i][j+1]-2\theta[i][j]+\theta[i][j-1]}{h_{y}^{2}} \label{eq6} \end{equation}

3.2 Discrétisation temporelle :

Le membre de gauche de la relation \eqref{eq2} peut s'exprimer ainsi :

\begin{equation} \dfrac{\partial \theta(x_{m},y_{m},t_{n})}{\partial t}=\dfrac{\theta(x_{m},y_{m},t_{n+1})-\theta(x_{m},y_{m},t_{n})}{\Delta t} \label{eq7} \end{equation} \begin{equation} \Rightarrow \,\,\,\,\dfrac{\partial \theta_{n}(x_{m},y_{m})}{\partial t}=\dfrac{\theta_{n+1}[i][j]-\theta_{n}[i][j]}{\Delta t} \label{eq8} \end{equation}

La formule de récurrence s'écrit donc en combinant les relations \eqref{eq2}, \eqref{eq6} et \eqref{eq8} :

\begin{equation} \theta_{n+1}[i][j]=\theta_{n}[i][j]+\kappa \Delta t\bigg[\dfrac{\theta_{n}[i+1][j]-2\theta_{n}[i][j]+ \theta_{n}[i-1][j]}{h_{x}^{2}}+\dfrac{\theta_{n}[i][j+1]-2\theta_{n}[i][j]+\theta_{n}[i][j-1]}{h_{y}^{2}}\bigg] \label{eq9} \end{equation}

Cette formule est implémentée de la façon ci-dessous pour la version parallélisée en langage C (voir fichier explUtilPar.c) ; à noter que les variables xs, xe, ys, ye (avec "s" pour "start" et "e" pour "end") représentent les coordonnées qui délimitent le sous-domaine correspondant au process de rank "me" :

   // Factors for explicit scheme
   diagx = -2.0 + hx*hx/(2*k0*dt);
   diagy = -2.0 + hy*hy/(2*k0*dt);
   weightx = k0 * dt/(hx*hx);
   weighty = k0 * dt/(hy*hy);

   // Perform an explicit update on the points within the domain
   for (i=xs[me];i<=xe[me];i++)
     for (j=ys[me];j<=ye[me];j++)
       x[i][j] = weightx*(x0[i-1][j] + x0[i+1][j] + x0[i][j]*diagx)
               + weighty*(x0[i][j-1] + x0[i][j+1] + x0[i][j]*diagy)

4.Critère de convergence :

La convergence relie la solution exacte à la solution calculée numériquement. Si l'on a $h_{x}=h_{y}=h$, le schéma explicite \eqref{eq9} est trouvé stable si :

\begin{equation} \Delta t\leqslant \dfrac{h^{2}}{4\kappa} \label{eq10} \end{equation}

Cette condition de stabilité nous donne une estimation pour le choix du pas de temps $\Delta t$. Dans l'algorithme, nous choisissons :

\begin{equation} \Delta t=\dfrac{1}{4\kappa}\,\text{min}(h_{x},h_{y})^{2} \label{11} \end{equation}

5.Implémentation MPI :

Nous abordons dans cette partie la manière dont est structurée et implémentée la résolution numérique grâce à l'API "Message Passing Interface" (MPI).

5.1 Topologie cartésienne de processus :

Nous allons décomposer le domaine (on dit aussi grille) à calculer en différents sous-domaines en affectant pour chacun d'entre eux un processus. Le gain sur le temps d'exécution est obtenu en opérant parallèlement les calculs sur chaque sous-domaine et en mettant à jour, pour chaque itération, leurs 4 cotés. On qualifie cette décomposition de "Topologie cartésienne de processus" car les sous-domaines sont carrés ou rectangulaires. Nous prenons pour convention, dans la version Fortran90 et C, les indices (i,j) des tableaux utilisés tels que (i,j)=(ligne,colonne). Ci-dessous une figure montrant une décomposition avec 8 processus selon x_domains=2 et y_domains=4 :

Figure 1 : Rangs des processus avec la convention (ligne,colonne)$\,\equiv\,$(x_domains,y_domains)

Pour créer cette topologie, nous utilisons la fonction MPI_Cart_create. La convention choisie par MPI pour cette fonction est le contraire de celle prise pour nos tableaux 2D, c'est-à-dire (i,j)=(row,column)$\,\equiv\,$(y_domains,x_domains) devient (row,column)$\,\equiv\,$(x_domains,y_domains) comme ceci est illustré sur la figure ci-dessus. Nous devons donc intervertir les axes (Ox) et (Oy) de la manière suivante :

Langage C :

// Invert (Ox,Oy) classic convention
dims[0] = y_domains;
dims[1] = x_domains;
MPI_Cart_create(comm, ndims, dims, periods, reorganisation, &comm2d);

Fortran90 :

! Invert (Ox,Oy) classic convention
dims(1) = y_domains
dims(2) = x_domains
call MPI_Cart_create(MPI_COMM_WORLD, ndims, dims, periods, reorganisation, comm2d, infompi)

5.2 Remarques sur la contiguïté :

En Fortran, les éléments d'un tableau 2D sont contigus en mémoire selon les colonnes : on dit qu'il est "column major". En langage C, c'est le contraire, les éléments sont contigus selon les lignes : il est qualifié de "row major".

Comme nous le verrons ci-dessous dans la partie 5.3, nous sommes amenés à échanger entre les processus des lignes et des colonnes de donnés. Naturellement, le langage C permet de manipuler les data sous forme de lignes et le langage Fortran sous forme de colonnes. Pour obtenir ces 2 fonctionnalités avec chacun de ces langages, on définit, dans la fonction updateBound, le type "colonne" (column_type) et "ligne" (row_type) respectivement pour le langage C et Fortran90.

Langage C :

// Create column data type to communicate with East and West neighBors
MPI_Type_vector(xcell, 1, size_total_y, MPI_DOUBLE, &column_type);
MPI_Type_commit(&column_type);

Fortran90 :

! Create row data type to communicate with South and North neighBors
call MPI_Type_vector(ycell, 1, size_tot_x, MPI_DOUBLE_PRECISION, row_type, infompi)
call MPI_Type_commit(row_type, infompi)

D'un point de vue optimisation, il faut veiller dans les boucles à itérer selon le premier indice ou le deuxième : la boucle la plus intérieure doit être exécutée sur le premier indice en Fortran90 et le deuxième en langage C.

5.3 Communications entre processus :

La parallélisation implique de faire communiquer entre eux les processus de chaque sous-domaine. En effet, comme nous pouvons le voir dans le bout de code ci-dessous, pour une itération donnée avec un processus de rank "me", nous calculons les prochaines valeurs du sous-domaine correspondant grâce à la fonction "computeNext", puis nous envoyons, avec la fonction "updateBound" les lignes et colonnes respectivement aux processus situés au Nord-Sud et Est-Ouest du processus courant "me".

   // Main loop : until convergence
   while (!convergence) {
     // Increment step and time
     step = step+1;
     t = t+dt;
     // Perform one step of the explicit scheme
     computeNext(x0, x, dt, hx, hy, &localDiff, me, xs, ys, xe, ye, k0);
     // Update the partial solution along the interface
     updateBound(x0, neighBor, comm2d, column_type, me, xs, ys, xe, ye, ycell);
     // Sum reduction to get global difference
     MPI_Allreduce(&localDiff, &result, 1, MPI_DOUBLE, MPI_SUM, comm);
     // Current global difference with convergence
     result= sqrt(result);
     // Break if convergence reached or step greater than maxStep
     if ((result<epsilon) || (step>maxStep)) break;
   }

La figure suivante est un exemple de communication pour une décomposition du domaine global en 4 sous-domaines, autrement dit avec 4 processus :

**Figure 2 : Exemple de communications entre 4 processus**

A noter que les cellules phantômes (ghost cells) sont initialisées en début de code à la valeur constante du bord de la grille.

Voici le contenu de la fonction "updateBound" qui permet aux processus de communiquer entre eux :

Langage C :

// North/South communication
// Send boundary to North and receive from South
MPI_Sendrecv(&x[xe[me]][ys[me]], ycell, MPI_DOUBLE, neighBor[N], flag,
             &x[xs[me]-1][ys[me]], ycell, MPI_DOUBLE, neighBor[S], flag,
             comm2d, &status);
// Send boundary to South and receive from North
MPI_Sendrecv(&x[xs[me]][ys[me]], ycell, MPI_DOUBLE, neighBor[S], flag,
             &x[xe[me]+1][ys[me]], ycell, MPI_DOUBLE, neighBor[N], flag,
             comm2d, &status);

// East/West communication
// Send boundary to East and receive from West
MPI_Sendrecv(&x[xs[me]][ye[me]], 1, column_type, neighBor[E], flag,
             &x[xs[me]][ys[me]-1], 1, column_type, neighBor[W], flag,
             comm2d, &status);
// Send boundary to West and receive from East
MPI_Sendrecv(&x[xs[me]][ys[me]], 1, column_type, neighBor[W], flag,
             &x[xs[me]][ye[me]+1], 1, column_type, neighBor[E], flag,
             comm2d, &status);

Fortran90 :

! North/South communication
! Send my boundary to North and receive from South
call MPI_Sendrecv(x(xe(me), ys(me)), 1, row_type, neighBor(N), flag, &
                  x(xs(me)-1, ys(me)), 1, row_type, neighBor(S), flag, &
                  comm2d, status, infompi)
! Send my boundary to South and receive from North
call MPI_Sendrecv(x(xs(me), ys(me)), 1, row_type, neighBor(S), flag, &
                  x(xe(me)+1, ys(me)), 1, row_type, neighBor(N), flag, &
                  comm2d, status, infompi)

! East/West communication 
! Send my boundary to East and receive from West
call MPI_Sendrecv(x(xs(me), ye(me)), xcell, MPI_DOUBLE_PRECISION, neighBor(E), flag, &
                  x(xs(me), ys(me)-1), xcell, MPI_DOUBLE_PRECISION, neighBor(W), flag, &
                  comm2d, status, infompi)
! Send my boundary to West and receive from East
call MPI_Sendrecv(x(xs(me), ys(me)), xcell, MPI_DOUBLE_PRECISION, neighBor(W), flag, &
                  x(xs(me), ye(me)+1), xcell, MPI_DOUBLE_PRECISION, neighBor(E), flag, &
                  comm2d, status, infompi)

5.4 Utilisation de MPI_Gather :

Pour collecter tous les sous-tableaux associés à chaque sous-domaine, nous utilisons la fonction MPI_Gather. Le mode de fonctionnement est différent entre les langages C et F90. Voyons 2 exemples montrant ces différences. Nous prenons les paramètres suivants (situés dans le fichier "param") :

  size_x
8
  size_y
4
  x_domains
4
  y_domains
2
  MaxStep
100000000
  dt
1.0e-1
  cnv_tol
1.0e-1

Une fois la convergence atteinte, on enregistre tous les sous-tableaux 2D en tableaux 1D en respectant les règles d'optimisation liés à la contiguïté spécifique à chaque langage. La fonction MPI_Gather permet alors de collecter tous ces sous-domaines en un seul tableau final 1D qui représente tout le domaine :

Langage C :

   // Gather all subdomains :
   // inner loop on columns index (second index)
   // to optimize since C is row major
   j=1;
   for (i=xs[me];i<=xe[me];i++) {
     for (k=0;k<ycell;k++)
       xtemp[(j-1)*ycell+k] = x0[i][ys[me]+k];
     j=j+1;
   }

   // Perform gathering
   MPI_Gather(xtemp, xcell*ycell, MPI_DOUBLE, xfinal, xcell*ycell, MPI_DOUBLE, 0, comm);

Comme on le voit ci-dessus, la boucle la plus intérieure porte sur le second indice, c'est-à-dire celui des colonnes, respectant ainsi le caractère "row major" du langage C dans l'accès des données.

Fortran90 :

   ! Gather all subdomains :
   ! inner loop on rows index (first index)
   ! to optimize since F90 is column major
   i = 1
   do j=ys(me),ye(me)
     xtemp((i-1)*xcell+1:i*xcell) = x0(xs(me):xe(me),j)
     i = i+1
   end do

   ! Perform gathering
   call MPI_Gather(xtemp, xcell*ycell, MPI_DOUBLE_PRECISION, xfinal, xcell*ycell, &
                   MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, infompi)

En Fortran90, la boucle la plus intérieure porte sur le premier indice, en relation avec le caractère "column major" du langage F90 dans l'accès des données.

6.Résultats :

Pour obtenir un gif animé représentant l'évolution temporelle de la diffusion, nous avons besoin des fichiers de données correspondant à des intervalles de temps différents. Nous utilisons pour cela le script generate_outputs_heat2d qui va les produire. Nous nous servons ensuite du script Matlab generate_gif_heat2d.m pour générer ce gif.

Voici une animation réalisée avec un maillage 256x96 :

**Figure 3 : Gif animé montrant l'évolution de la température sur le domaine**

7.Benchmark :

Pour évaluer les performances du code, nous accomplissons un benchmark en faisant varier le nombre de processus pour 3 tailles de grilles différentes (512^2, 1024^2 et 2048^2). Le script run_benchmark_heat2d permet de générer les temps d'exécution pour chacun de ces deux paramètres.

Le graphique suivant, obtenu avec le script Matlab plot_benchmark_heat2d.m, synthétise ce benchmark :

**Figure 4 : SpeedUp en fonction du nombre de processus et de la taille du domaine**

Les performances pour 1 seul processus correspondent à l'exécution du code séquentiel. Excepté pour N=2048, les runtimes pour un nombre de processus compris entre 2 et 128 sont inférieurs à la version séquentielle. Nous remarquons que les meilleures performances sont obtenues avec 8 processus, indépendamment de la taille du maillage. Ceci peut s'expliquer par le fait que ce benchmark ait été réalisé avec un processeur i7 8 coeurs. En effet, le processeur i7 implémente "l'Hyper-Threading" ; dans notre cas, cette technique multitâche permet d'avoir une optimisation physique pour 8 processus et une optimisation logique pour un nombre plus élevé. A noter que pour un nombre de processus égal à 256, nous retrouvons sensiblement les mêmes temps que ceux obtenus avec le code séquentiel, ce qui traduit le surcoût engendré par le nombre élevé de communications sur le temps de calcul effectif dans chaque sous-domaine de la grille.

ps : contribuez comme moi au projet Cosmology@Home dont le but est d'affiner le modèle décrivant le mieux notre Univers.

Haut de Page