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

DOURNAC.ORG
	Français English

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

1.Exécution :

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

Version Langage C : Heat3D_MPI_Solving_C.tar.gz
Version Fortran 90 : Heat3D_MPI_Solving_F90.tar.gz

Avertissement : pour les versions précédentes avant le 23/05/18, bug fixé sur le critère de convergence (le pas d'espace $h$ doit être élevé au carré et non au cube)

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

$ cat param
  size_x
16
  size_y
8
  size_z
32
  x_domains
8
  y_domains
4
  z_domains
2
  MaxStep
1000000
  dt
1.0e-1
  cnv_tol
1.0e-1

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

$ mpirun -np 64 ./explicitPar

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

   Time step =   3.18033787909627519E-006

   Convergence = 0.100000000 after      3333 steps

   Problem size =         4096

   Wall Clock =        2.126691

   Computed solution in "outputPar.dat" file

2.Équation de chaleur instationnaire 3D :

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 à 3 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}}+\dfrac{\partial^{2}\theta}{\partial z^{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 les variable $y$ et $z$. En prenant $h_{x}=size_{x}/N_{x}$, $h_{y}=size_{y}/N_{y}$ et $h_{z}=size_{z}/N_{z}$ avec la convention $\theta(x_{m},y_{m},z_{m})=\theta[i][j][k]$, on obtient :

\begin{align} \dfrac{\partial^{2}\theta}{\partial x^{2}}+\dfrac{\partial^{2}\theta}{\partial y^{2}}+\dfrac{\partial^{2}\theta}{\partial z^{2}} &= \dfrac{\theta[i+1][j][k]-2\theta[i][j][k]+\theta[i-1][j][k]}{h_{x}^{2}}\label{eq6}\\ & +\dfrac{\theta[i][j+1][k]-2\theta[i][j][k]+\theta[i][j-1][k]}{h_{y}^{2}}\notag \\ & +\dfrac{\theta[i][j][k+1]-2\theta[i][j][k]+\theta[i][j][k-1]}{h_{z}^{2}} \notag \end{align}

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},z_{m},t_{n})}{\partial t}=\dfrac{\theta(x_{m},y_{m},z_{m},t_{n+1})-\theta(x_{m},y_{m},z_{m},t_{n})}{\Delta t} \label{eq7} \end{equation} \begin{equation} \Rightarrow \,\,\,\,\dfrac{\partial \theta_{n}(x_{m},y_{m},z_{m})}{\partial t}=\dfrac{\theta_{n+1}[i][j][k]-\theta_{n}[i][j][k]}{\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{align} \theta_{n+1}[i][j][k] &= \theta_{n}[i][j][k]+\kappa \Delta t\bigg[\dfrac{\theta_{n}[i+1][j][k]-2\theta_{n}[i][j][k]+\theta_{n}[i-1][j][k]}{h_{x}^{2}}\label{eq9}\\ &+\dfrac{\theta_{n}[i][j+1][k]-2\theta_{n}[i][j][k]+\theta_{n}[i][j-1][k]}{h_{y}^{2}}\notag \\ &+\dfrac{\theta_{n}[i][j][k+1]-2\theta_{n}[i][j][k]+\theta_{n}[i][j][k-1]}{h_{y}^{2}} \notag \bigg] \end{align}

Cette formule est implémentée de la façon ci-dessous pour la version parallélisée en langage C (fichier explUtilPar.c) ; à noter que les variables xs, xe, ys, ye, zs, ze 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);
   diagz = -2.0 + hz*hz/(2*k0*dt);
   weightx = k0 * dt/(hx*hx);
   weighty = k0 * dt/(hy*hy);
   weightz = k0 * dt/(hz*hz);

   // 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++)
         for (k=zs[me];k<=ze[me];k++)
            x[i][j][k] = weightx*(x0[i-1][j][k] + x0[i+1][j][k] + x0[i][j][k]*diagx)
                       + weighty*(x0[i][j-1][k] + x0[i][j+1][k] + x0[i][j][k]*diagy)
                       + weightz*(x0[i][j][k-1] + x0[i][j][k+1] + x0[i][j][k]*diagz);

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_{z}=h$, le schéma explicite \eqref{eq9} est trouvé stable si :

\begin{equation} \Delta t\leqslant \dfrac{h^{2}}{8\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}{8\kappa}\,\text{min}(h_{x},h_{y},h_{z})^{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 3D (appelé aussi "grille") à calculer en différents sous-domaines cubiques 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, les 6 faces de chaque cube. On qualifie cette décomposition de "Topologie cartésienne de processus" car les sous-domaines sont des cubes ou des parallépipèdes. 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 16 processus avec x_domains=2, y_domains=4 et z_domains=2 :

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

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".

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 faces de chaque cube (voir figure ci-dessous), c'est-à-dire les matrices, aux processus situés aux faces Nord-Sud, Est-Ouest, Avant-Arrière, du processus courant "me".

   // Until convergence
   while(!convergence)
   {
     step = step + 1;
     t = t + dt ;

     // Perform one step of the explicit scheme
     computeNext(x0, x, dt, &resLoc, hx, hy, hz, me, xs, ys, zs, xe, ye, ze, k0);

     // Update the partial solution along the interface
     updateBound(x0, neighBor, comm3d, column_type, me, xs, ys, zs, xe, ye, ze,
                 xcell, ycell, zcell);

     // Sum reduction to get error
     MPI_Allreduce(&resLoc, &result, 1, MPI_DOUBLE, MPI_SUM, comm);

     // Current error
     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 3D du domaine global en 8 sous-domaines, autrement dit avec 8 processus. Un seul plan de 4 processus est affiché, il faut imaginer un second plan parallèle au premier pour représenter l'ensemble des 8 sous-domaines :

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

Comme dans le cas 2D, les cellules phantômes (ghost cells) sont initialisées en début de code à la valeur constante du bord de la grille.

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_heat3d qui va les produire. Nous nous servons ensuite du script Matlab generate_gif_heat3d.m pour générer ce gif. Voici une animation réalisée avec un maillage 32x64x16 :

**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 (64^3, 128^3 et 256^3). Le script run_benchmark_heat3d 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_performances_heat3d.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. Pour un nombre de processus supérieur à 1 (versions parallélisées), les runtimes 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é.

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