Projets de spécialité 2014 : Modele Electrophysiologique Sur GPU

De Ensiwiki
Aller à : navigation, rechercher
Project schedule.png
Titre du projet Modele Electrophysiologique Sur GPU
Cadre Projets de spécialité
Page principale Projets de spécialité Calcul Scientifique

Encadrants Picard Christophe


L'équipe

Keck Jean-Baptiste (2A - MMIS)

Zirnhelt Gauthier (2A - MMIS)

Contexte

Introduction

La plupart des simulations numériques en ingénierie et en sciences naturelles requièrent des moyens de calcul d'une grande efficacité, aussi bien d'un point de vue algorithmique que d'un point de vue matériel. Ces dernières années ont vu l'émergence de matériels jusque là réservés aux applications graphiques : Les GPU (Graphic Processing Units). Comparativement aux technologies existantes, les cartes graphiques (GPU) sont bon marché, accessibles et peu gourmandes en énergie. Le fait d'effectuer ces types de calculs non graphiques sur GPU est plus communément appelé GPGPU (General Purpose Computing on Graphic Processing Units).

Dotées de milliers de cœurs de calcul et de par leur nature massivement parallèle, les cartes graphiques permettent la démocratisation du calcul intensif et ouvrent cette discipline à de nombreuses applications scientifiques. Dans les stratégies diagnostiques et thérapeutiques en médecine, l'outil informatique est devenu indispensable mais reste fortement inexploité. Il permet d'obtenir et d'analyser précisément l'état d'un patient à un instant donné afin de lui fournir les soins les plus adaptés.

Setup
Illustration du GPGPU : Le processeur épaulé des milliers de cœurs d'un GPU.

Le projet européen Virtual Physiological Human (VPH) a pour but de simuler un corps humain au niveau cellulaire en tant qu'un unique système complexe, d'ici 2050. Cet ambitieux projet combine les compétences des médecins, informaticiens, et mathématiciens. De par le coût prohibitif d'une telle simulation, le GPGPU va jouer un rôle crucial dans l'aboutissement de ce projet. Ce projet permettra à terme de virtualiser l'intégralité du corps humain afin de par exemple simuler l'effet de traitements sur un patient et de détecter leurs conséquences potentiellement néfastes.

Setup
Les différents laboratoires qui travaillent sur le projet VPH.

Dans le cadre de ce projet de spécialité nous nous intéressons à un sous ensemble du corps humain, le cœur, et plus précisément à la simulation des courants trans-membranaires dans les cellules cardiaques.

Objectifs

Notre premier objectif était de simuler un modèle bidimensionnel appartenant au domaine de l'électrophysiologie (modèle cardiaque). Ces simulations devaient pouvoir être calculées et affichées en temps réel dans une interface utilisateur lorsque cela était possible. Parallèlement, nous devions permettre à l'utilisateur de stocker les données de simulation pour pouvoir les analyser par après, de choisir les paramètres de la simulation et ses conditions initiales dans une interface graphique accessible. Nous souhaitions également avoir la possibilité de simuler dans un environnement multi-GPU. Enfin, le dernier objectif, si le temps le permettait, consistait à passer nos simulations en 3D, et de pouvoir visualiser les résultats sous forme de coupes 2D.


Simulation

Modèle

Nous avons choisi comme premier modèle d'expérimentation un modèle générique et simple simulant le comportement des cellules cardiaques (il peut aussi servir pour simuler des neurones).

Le milieu représenté est un milieu excitable et comprend seulement deux variables d'état.

Ce modèle est décrit par le modèle mathématique suivant :


\begin{cases}
\frac{\partial e}{\partial t} = F(e,r) = \nabla {}^t \! D \nabla e - Ke(e-\alpha_1)(e-1) - er)\\ 
\frac{\partial r}{\partial t} = G(e,r) = \left ( \varepsilon + \mu_1\frac{r}{1+\mu_2} \right ) (-r-Ke(e-\alpha_2-1))
\end{cases}

D est un tenseur représentant les caractéristiques du milieu qui nous intéresse (la conductivité dans cet exemple). Dans le cas général, le milieu est anisotrope, c’est-à-dire que D dépend de la direction des tissus.

e représente le potentiel d'action trans-membranaire, r représente la conductance du courant repolarisant et t représente le temps. Les autres grandeurs sont les paramètres du modèle et n'ont pas d'interprétation physique.


En simplifiant l'opérateur \nabla^t D \nabla du modèle en considérant la conductivité constante dans le domaine, le tenseur D devient diagonal :


\begin{cases}
\frac{\partial e}{\partial t} = D\Delta e - Ke(e-\alpha_1)(e-1) - er)\\ 
\frac{\partial r}{\partial t} = \left ( \varepsilon + \mu_1\frac{r}{1+\mu_2} \right ) (-r-Ke(e-\alpha_2-1))
\end{cases}


Ce modèle simplifié fait apparaître un laplacien D \Delta que nous discrétisons plus bas dans cet article.

Discrétisation

Discrétisation temporelle

Notre modèle repose sur la méthode d’Euler explicite appliquée à chacune des équations.

En considérant l'équation différentielle suivante :

\dfrac{du(t)}{dt} = f(t, u(t)) avec u(t=0)=u_0

où f est une fonction continue sur un ouvert de [0,T]\times\R à valeur dans \R et u_0 \in \R.


Avec la formule de Taylor à l'ordre 1 on exprime la dérivée temporelle en t_k : \dfrac{du(t_k)}{dt} = \frac{1}{dt}\left( u(t_{k+1})-u(t_k) \right) + O(dt)


Et donc en substituant la première équation on obtient : u(t_{k+1}) = u(t_k) + dt\;f(t_k, u(t_k)) + O(dt^2)

Si l'on note u^k l'approximation de u en t_k, on obtient l'équation différentielle approchée:

u^{k+1} = u^k + dt\;f(t_k, u^k)


Pour notre modèle on a donc :


\begin{cases}
e^{k+1} = e^k + dt\;F(e^k,r^k) = e^k + dt\;\left( D \Delta e^k-Ke^k(e^k-\alpha_1)(e^k-1)-e^kr^k)) \right)\\ 
r^{k+1} = r^k + dt\;G(e^k,r^k) = r^k + dt\;\left ( \varepsilon + \mu_1\frac{r^k}{1+\mu_2} \right ) \left(-r^k-Ke(e^k-\alpha_2)\right)
\end{cases}

Approximation par différences finies

La première étape consiste à diviser chacune des directions de l’espace en un nombre fini d’intervalles. Pour simplifier, on se place sur \Omega = [0,1]^2. Cela permet d’obtenir une grille régulière de l’espace.

Les noeuds de la grille ont pour coordonnées (x_i,y_j) = (i\;h,j\;h) avec h = \frac{1}{n-1} et 0\le i,j \le n-1 où n est un entier.

Cette simplification du domaine ne permetpas de prendre en compte toute la complexité du coeur, mais suffit pour visualiser certains phénomènes physiologiques.

On note  e_{i,j}^k l'approximation de e(x_i,y_j) en t_k et on fait de même pour la variable r.

En utilisant la formule de Taylor à l'ordre 2, on exprime les dérivées secondes et on les somme :  
\Delta u = \frac{1}{h^2}\left( u(x_{i+1},y_j) + u(x_{i-1},y_j) + u(x_{i},y_{j+1}) + u(x_{i},y_{j-1}) -4u(x_i,y_j) \right)


Ceci définit un laplacien approché à l'ordre deux :

 
\Delta_h u = \frac{1}{h^2}\left( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} -4u_{i,j} \right)

Un autre laplacien approché à l'ordre deux est :

 
\Delta_h u = \frac{1}{h^2}\left( 
u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1}
+u_{i+1,j+1} + u_{i-1,j+1} + u_{i+1,j+1} + u_{i+1,j-1} 
-8u_{i,j} \right)


Ceci peut être résumé dans un stencil comme représenté ci-dessous.

Exemple de stencil de laplacien discret 2D.
Un autre stencil de laplacien 2D.

Pour des raisons pratiques on choisit le premier laplacien. On peut bien voir son application sur la grille régulière suivante :

Setup
Le domaine de simulation est discrétisé sur une grille cartésienne. Le premier stencil de laplacien a été superposé sur la grille.

L'étape de résolution avec ce modèle explicite s'écrit donc :


\forall (i,j) \in [[0,n-1]] \;\;
\begin{cases}
e_{i,j}^{k+1} = e_{i,j}^k + dt\;
\left( 
\frac{D}{h^2} \left( e^k_{i+1,j} + e^k_{i-1,j} + e^k_{i,j+1} + e^k_{i,j-1} -4e^k_{i,j} \right)
-Ke_{i,j}^k(e_{i,j}^k-\alpha_1)(e_{i,j}^k-1)-e_{i,j}^kr_{i,j}^k)) 
\right) \\ 
r_{i,j}^{k+1} = r_{i,j}^k + dt\;\left ( \varepsilon + \mu_1\frac{r_{i,j}^k}{1+\mu_2} \right ) \left(-r_{i,j}^k-Ke_{i,j}^k(e_{i,j}^k-\alpha_2)\right)
\end{cases}

Pour les bords du domaine on prend des conditions aux limites de Neumann (on prolonge les valeurs au bord selon la normale orientée vers l’extérieur du domaine).

Implémentation

Technologies utilisées

Nous avons choisi OpenCL pour le calcul sur GPU et OpenGL pour le rendu, ainsi que Qt4 pour ce qui est de la partie interface graphique utilisateur. Le code est écrit en C++11 avec une toute petite partie en python 2.7. Les systèmes cibles sont ceux de type GNU/Linux, possédant une ou plusieurs cartes graphiques compatibles OpenCL 1.1 ou plus et OpenGL 3.1 ou plus. Le compilateur de référence est gcc 4.8.


On utilise également d'autres librairies :

  • Glx pour gérer les contextes OpenGL
  • Glew pour gérer les extensions OpenGL
  • Glut pour les routines OpenGL (pour afficher du texte à l'écran)
  • Log4cpp pour avoir des logs robustes
  • La librairie python pour interfacer python 2.7 avec le C++ (choix de condition initiale sous forme d'expression ou de fonction python)
  • Boost pour les barrières de synchronisation entre threads

Premiers essais

Sur CPU

Nous avons d'abord effectué des tests sur CPU. Les données étaient exportées au format matriciel pour réaliser des gif à l'aide de gnupot afin d'avoir une sortie visuelle. La génération de gif est plus longue que la simulation elle même, mais le tout reste significativement plus rapide que scilab, et ce sans génération de gif !

Setup
Vidéo de la première simulation. Si cela ne marche pas dans votre navigateur le plus simple est d'utiliser wget.

Sur GPU

Il fallait tout d'abord réaliser des fonctions utilitaires pour OpenCL (macro de gestion d'erreurs, callbacks, parsing de sources OpenCL, ...). Le plus gros problème venait des headers fournis par Nvidia pour le wrapper C++ d'OpenCL. En effet Nvidia fournit le header C++ pour OpenCL 1.2 mais ne fournit que la librairie dynamique OpenCL 1.1 (le matériel Nvidia ne supporte actuellement pas plus). Nous avons donc dû manuellement inclure un header version 1.1.

Après avoir testé quelques exemples de kernels OpenCL, nous avons implémenté le modèle sur GPU, toujours avec une sortie en gif. Nous avons comparé les résultats jusqu'à qu'ils soient identiques aux résultats sur CPU.

Implémentation multi-GPU

OpenGL et Qt

Le problème avec Qt4, c'est qu'on ne sait pas comment est géré le contexte OpenGL derrière l'interface graphique. De plus la documentation de Qt sur OpenGL est quasi inexistante, et il est donc très difficile de s'interfacer avec Qt.

Dans une slide de présentation de Qt5, on peut trouver: "The Dark Times: Achieving windowed, threaded OpenGL support in Qt v4.4" (Qt5 propose apparemment une interface beaucoup plus limpide pour la gestion interne d'OpenGL...).

Comme si cela ne suffisait pas, la documentation sur la création de contextes multiples OpenGL est également très rare et dépend de chaque plateforme. Même le dernier Red Book n'y consacre que 2 pages (sur 935).

Il faut savoir qu'un contexte OpenGL est attaché à un thread. Lorsque l'on réalise une application multithreadée où plusieurs threads sont susceptibles d'appeler des commandes OpenGL, plusieurs possibilités s'offrent à nous :

  • Soit passer un seul contexte à différents threads (dans ce cas il faut manuellement garantir la mutuelle exclusion du changement de contexte entre les threads avec des mutex ou des sémaphores, et bien gérer la machine à états OpenGL)
  • Soit créer plusieurs contextes, un par thread, et laisser le driver ordonnancer les commandes OpenGL lui même.

Après plusieurs tentatives infructueuses, et trois drivers Nvidia plus tard, nous avons enfin réussi à faire cohabiter un contexte OpenGL 3.1 avec Qt4 avec la deuxième possibilité. Différents contextes OpenGL peuvent partager des données entre eux, et nous tirons parti de ceci afin de faire passer des textures entre le solveur et qt.

Sous Linux, la carte graphique qui effectue le rendu OpenGL est celle liée au X Screen de l'écran où se situe la fenêtre liée au contexte (c'est un peu technique, il faut comprendre le fonctionnement de X11). Il est donc quasiment impossible de choisir la carte graphique qui effectue le rendu (à moins d'avoir un X Screen par carte et de créer de faux écrans, ce que nous avons tenté, en vain).

Nous voulions initialement tirer parti de l'interaction entre OpenCL et OpenGL afin de ne plus passer par la RAM pour effectuer le rendu. L'idée initiale était donc de créer un thread par carte graphique, avec son contexte OpenGL et son contexte OpenCL, partagés, afin d'effectuer des rendus 'offscreen' des sous-grilles, puis de laisser le thread de Qt dessiner tous ces bouts de texture à l'écran. Cela s'est finalement avéré impossible car la création d'un contexte partagé OpenGL/OpenCL n'a jamais pu se faire pour des raisons qui restent encore inconnues.

Découpe du domaine

Pour pouvoir charger plusieurs cartes graphiques il faut pouvoir découper le problème en sous-problèmes. On découpe donc la grille initiale en sous-grilles (qui n'ont pas forcément la même taille). Le laplacien présent dans le modèle nous impose d'envoyer en plus de ces sous-grilles les bords de ces sous-grilles pour pouvoir appliquer le stencil et calculer les valeurs aux bords. Celà se voit bien dans la figure ci-dessous :

Découpe d'une grille 2D
Découpe d'une grille 2D. Les coupes sont en rouge, les bords internes de la dernière sous-grille sont en vert et ses bords externes sont en jaune.

Déroulement de la simulation

Après la découpe (sur CPU), on recherche les périphériques OpenCL de type GPU (OpenCL est thread-safe contrairement à OpenGL, sauf pour les paramètres d'un même kernel). Pour chacun de ceux-ci, on crée un 'worker thread' qui va prendre les sous-grilles sur un modèle classique de producteur-consommateur. Chacun de ces threads tente de s'approprier un sous-domaine, l'initialise avec sa condition initiale, et initialise ses bords externes. Tant que la carte graphique dispose d'encore assez de mémoire et qu'il y a encore des sous-domaines disponibles cela continue jusqu'à que tout soit initialisé.

La simulation à proprement parler se déroule sur le même modèle. À chaque itération i, on prend soin de récupérer les bords internes de chaque sous-grilles, qui deviendront les bords externes d'une des 4 grilles voisines (6 en 3D, voir l'image) à l’itération i+1, sur une autre carte graphique.

Tous ces threads sont synchronisés entre eux via de multiples barrières. Il faut également les synchroniser avec Qt. De plus il faut synchroniser OpenCL et OpenGL entre chaque itération pour que le dessin à l'écran se fasse bien (toutes les commandes sont asynchrones).


Pour vous donner une idée de ce que représente cette implémentation avec notre modèle de base voici les arguments du kernel utilisé pour la simulation multi-gpu :


Aperçu du kernel utilisé pour la simulation multi-gpu avec notre modèle de base.

Architecture logicielle

Rapide aperçu

L'application se divise en deux parties principales.

  • La partie interface graphique utilisateur qui comprend surtout du code utilisant Qt et OpenGL.
  • La partie simulation/modèles qui concentre les classes servant à définir et organiser le travail entre les différentes cartes graphiques.

Les diagrammes ci-dessous montrent de façon relativement détaillée les architectures de ces deux parties et le rôle des différentes classes dans celles-ci.


Diagramme UML de la partie GUI
Diagramme UML de la partie interface graphique utilisateur.
Diagramme UML de la partie modèles
Diagramme UML de la partie simulation/modèles.

Répartition du code final

Voici un aperçu de la répartition du code à la fin de ce projet :

Répartition des 10.6k lignes de code
Répartition des 10.6k lignes de code.

Ajout de nouveaux modules

L'architecture de notre application est volontairement modulaire afin de pouvoir ajouter facilement de nouveaux modèles que nous n'avons pas eu le temps d'ajouter nous mêmes. Ceux-ci dérivent tous de la classe Model qui constitue aussi l'interface entre les modèles et la partie graphique.

Malheureusement, nous n'avons pas eu le temps de refactoriser notre code et il faudra aussi compléter deux switch se trouvant dans les classes MainWindow et SidePanel.

L'ajout de nouvelles conditions initiales se fait facilement en prenant pour exemple celles se trouvant déjà dans la classe InitialCondFactory.

Affichage

Les champs scalaires 2D du modèle sont récupérés des cartes graphiques environ 60 fois par seconde lorsque cela est possible (écran 60Hz), et on reconstruit pour chaque variable un tableau contigu dans la RAM. De ces tableaux, on crée des textures 2D dans le contexte OpenGL du solveur. Comme les deux contextes OpenGL sont partagés, le contexte OpenGL de Qt peut lire ces textures et envoyer les commandes nécessaires à l'affichage (génération des rectangles et mapping des textures). Un shader a été écrit pour passer de la représentation scalaire normalisée [0,1] aux couleurs finales.

Les palettes des couleurs (une vingtaine de listes d'une dizaine de teintes) sont stockées dans un Uniform Buffer Object (UBO) ce qui permet d'une part de pouvoir accéder à ces listes dans le shader, mais aussi d'autre part de ne les envoyer qu'une seule fois sur la carte graphique qui effectue le rendu, ce qui évite des transferts inutiles à chaque rendu.

On peut lisser un peu l'affichage des grilles 2D en changeant un paramètre du sampler de l'unité de texture. Avec l'option GL_NEAREST, on affiche la vraie grille de simulation. En passant en GL_LINEAR on obtient une moyenne pondérée avec les 4 valeurs adjacentes ce qui donne un aspect plus lisse, et ce gratuitement (interpolation matérielle).

Quelques résultats

Voici quelques résultats obtenus avec diverses conditions initiales sur notre modèle de base :

Spirales1.png Spirales2.png Spirales3.png
Forme1.png Forme2.png Forme3.png
Cercle.png Forme4.png Forme6.png

Fonctionnement de l'interface utilisateur

Utilisation

L'interface utilisateur se compose d'une zone d'affichage (à gauche) et d'une zone de paramétrage (à droite). Cette dernière est découpée en trois parties:

  • Sélection du modèle à simuler, initialisation des variables et modification des paramètres de celui-ci, ainsi que sélection du nombre d'itérations à effectuer. Ces options ne sont pas modifiables en cours de simulation.
  • Boutons permettant de lancer, pauser, reprendre, ou stopper une simulation. Options permettant de sauvegarder les données de simulation. Ces dernières ne sont pas non plus modifiables en cours de simulation.
  • Sélection des variables à afficher dans la zone de gauche. Changement de palette de couleurs. Ces options sont toujours accessibles, même en cours de simulation comme elles n'affectent que le rendu.

Le bouton "Initialization" ouvre une fenêtre supplémentaire laissant à l'utilisateur le choix de la taille des grilles utilisées ainsi que celui de l'initialisation des grilles correspondant à chaque variable du modèle.

Le bouton "Parameters" ouvre une fenêtre permettant de modifier les valeurs des différents paramètre du modèle.

Le nom des variables rendues à l'écran sont affichés au dessus de celles-ci.

En bas de la fenêtre se situe une barre qui permet de suivre l'avancement de la simulation.

Les raccourcis disponibles sont :

  • 'Echap' pour fermer la fenêtre
  • 's' pour changer le mode d'affichage des textures entre GL_NEAREST et GL_LINEAR (GL_LINEAR par défaut)


Interface utilisateur
L'interface utilisateur

Les différents éléments de l'interface

Interface utilisateur
Fenêtre d'initialisation des variables
Interface utilisateur
Fenêtre de modification des paramètres d'un modèle
Interface utilisateur
L'éditeur de code python pour les condition initiales

La zone d'affichage des variables

Affichage des deux variables du modèle de base à un instant t.
Affichage de ces mêmes variables quelques secondes plus tard, avec une autre palette de couleurs.

Perspectives d'amélioration de l'application

  • Mesures de performance sur une machine avec plus de cartes graphiques et optimisations si nécessaire.
  • Implémentation d’autres modèles en électrophysiologie.
  • Implémentation d’autres méthodes de résolution plus efficaces et plus stables.
  • Affichage de coupes avec une simulation sur grilles en 3 dimensions.
  • Refactorisation du code pour permettre un ajout plus simples de nouveaux modèles ou solveurs.
  • Affichage de l'échelle des couleurs.


Résultats

Ces résultats ont été mesurés dans un environnement multi-GPU (deux GTX 770 comprenant 4Go de mémoire embarquée et 1536 cœurs de calcul chacune) et sur un processeur Intel Core i7-3770K @ 3.50GHz. W est la largeur de la grille considérée, H sa hauteur (on se place en 2D). La simulation effectue 1000 pas.

W H CPU -O2 multithread (ms) Mono GPU (ms) gains/CPU Multi GPU (ms) gains/CPU
128 128 3916 67 5 840% 1372 290%
256 256 12 360 136 9 090% 1 489 831%
512 512 71 531 452 15 825% 1 404 5 095%
1024 1024 274 602 1 213 22 638% 2 436 11 273%
2048 2048 1 103 000 4 154 26 553% 8 585 12 848%
4096 4096 3 454 700 14 329 24 110% 30 891 11 183%
8192 8192 - 55 844 - 102 907 -

On voit qu'on l'on a des gains significatifs par rapport à l'implémentation CPU multithread et avec un impressionnant gain de 265x avec l'implémentation mono-GPU. L'implémentation multi-GPU souffre de la synchronisation entre threads, et de la découpe des grilles, ce qui résulte en une performance deux fois moindre par rapport à l'implémentation mono-GPU. En effet cette implémentation impose le transferts des bords des sous-grilles entre GPU, qui passent d'abord par la RAM tandis que l'implémentation mono-GPU réside uniquement sur une seule carte et les données ne passent jamais par le CPU (à part à l'initialisation). Cependant, l'implémentation multi-GPU est actuellement la seule qui permet de simuler des grosses grilles (au delà de 8192x8192, ou 512x512x256 en 3D) et c'est une bonne chose que l'impact sur la performance n'est que d'un facteur 2 par rapport à une seule carte malgré tous les transferts engendrés !

Bilan et conclusion

  • Ce projet nous a permis de découvrir et maîtriser de nouvelles technologies (OpenCL, Qt), ainsi que de consolider nos connaissances sur celles que nous connaissions déjà (C++, OpenGL).
  • C'est une application concrète du cours de systèmes d'exploitation et programmation concurrente (synchronisation, utilisation de mutex, ...).
  • Il est possible de réaliser des simulations de taille suffisamment grande pour la pratique et avec affichage avec du matériel relativement récent.
  • Nous aurions aimé avoir plus de temps ou une personne en plus afin de pouvoir implémenter plus de modèles et de méthodes de résolution et refactoriser notre code en fin de projet.
  • Plus jamais Qt4 avec du multi-contexte OpenGL.

Fichiers

Sujet de Méthodes Numériques dont notre projet est la continuité logique : Sujet

Liste de modèles de cellules cardiaques : Article

Slides de la soutenance : Fichier:Présentation Projet de spécialité GpGpu.pdf