Gaétan BAHL (Encadrant : Emmanuel Maître) : Solutions de viscosité d'EDP et schémas numériques associés

De Ensiwiki
Aller à : navigation, rechercher

Cette page présente les résultats du sujet "Solutions de viscosité d'EDP et schémas numériques associés". Étudiant : Gaétan BAHL (2A MMIS 2016)


Cornues.png
Titre du projet Solutions de viscosité d'EDP et schémas numériques associé
Cadre IRL

Labo LJK-IMAG
Équipe EDP
Encadrants Emmanuel Maître

Introduction

Cette Introduction à la Recherche en Laboratoire s’intéresse à une nouvelle notion de solution d’équations aux dérivées partielles, plus générale que les définitions classiques. Ceci permet de considérer des cas que les définitions classiques ne permettent pas de résoudre. Nous nous sommes plus spécifiquement intéressés à ce que cette nouvelle notion peut nous apporter dans le cadre du transport optimal, en implémentant et optimisant un schéma numérique résolvant l’équation de Monge-Ampère.

Équipe d'accueil

J'ai réalisé mon IRL au Laboratoire Jean Kuntzmann (LJK-IMAG), au sein de l’équipe Équations aux Dérivées Partielles (EDP) et sous la tutelle d’Emmanuel Maître, professeur à Grenoble INP et membre de cette équipe.


Travail réalisé

Notion de solution de viscosité

La première partie du travail a été de comprendre la notion de solution de viscosité, car celle-ci n'est pas étudiée en 2A à l'Ensimag. Le module Modèles EDP et Schémas Numériques en Sciences de l’Ingénieur, enseigné par Emmanuel Maître, mon tuteur, permet néanmoins de comprendre la notion classique de solution d'EDP, ce qui constitue une bonne introduction au sujet.

J'ai commencé par lire des articles de recherche à ce sujet (cf. Bibliographie), certains présentant la notion de solution de viscosité [2], d'autres présentant des exemples et des méthodes de résolution [1][4], et d'autres enfin présentant des applications de cette notion [1][5][6]. Celle-ci permet de trouver des solutions d'équations aux dérivées partielles non-linéaires, ce qui a notamment des applications en finance, et en physique (contrôle optimal, théorie des jeux,...).


Un exemple introductif est l'équation eikonale avec  f \equiv 1 , qui s'écrit  - |Du| +1 =0 . Cette équation n'admet pas de solution classique  C^1 , mais elle admet une solution de viscosité (représentée ci-dessous), ce qui est dû à la nature plus faible de ce type de solution.

Solution de viscosité de l'équation eikonale

Je ne vais pas présenter sur cette page la définition exacte d'une solution de viscosité, car celle-ci requiert de définir d'autres notions auparavant (sous-différentielle, sur-différentielle,...). Il faudra donc se référer à mon rapport, à [2], ou encore au module Optimisation Numérique, enseigné par Jérôme Malik au premier semestre de 2A.


Implémentation d'une méthode de résolution

Description du problème

J'ai choisi de m'intéresser à ce que peuvent nous apporter les solutions de viscosité dans le cadre du transport optimal, et d'implémenter une méthode de résolution présentée dans une publication [6]. il s'agissait de résoudre numériquement l'équation de Monge-Ampère, afin d'obtenir une application transformant une densité de probabilité en une autre avec un coût minimal.

Nous disposons de deux mesures de probabilité, \rho_X sur X et \rho_Y sur Y, avec X,Y \in \mathbb{R}^n. Il s'agit de trouver une application \mathbb{M} appartenant à l'ensemble des applications transformant la densité \rho_Y en la densité \rho_X,  \lbrace T : X \to Y, \rho_Y(T) \det(\nabla T) = \rho_X \rbrace .

Cette application doit être optimale au sens de la fonction de coût quadratique :

 \frac{1}{2} \int_X ||x-T(x)||^2\rho_X(x) dx

Cela nous donne donc l'équation de Monge-Ampère suivante :

 \det(D^2 u) = \frac{\rho_X}{\rho_Y(\nabla u)}, \; x \in \Omega

Pour résoudre cette équation, on utilise trois schémas numériques, un schéma monotone, un schéma précis, et un schéma pour la condition au bord. J'ai détaillé ces schémas dans mon rapport.

Implémentation sous Scilab

J'ai commencé par implémenter ces schémas sous Scilab, en utilisant une itération explicite (type Euler). Le schéma précis fonctionne mais il subsiste un soucis que je n'ai pas eu le temps de résoudre dans le schéma monotone. On arrive néanmoins à calculer des solutions. Le problème de l'itération explicite est que la convergence se fait très lentement, la condition CFL étant très restrictive.

Sur la figure ci-dessous, on voit un exemple de déplacement d'une gaussienne, calculé par mon implémentation.

Calcul du déplacement optimal d'une gaussienne

L'application de transport calculée est représentée par un champ de vecteurs sur la figure ci-dessous.

Application de déplacement optimal d'une gaussienne

Pour d'autres exemples, tels que le transport d'une gaussienne vers la densité uniforme, ou le déplacement d'un disque, voir mon rapport.


Implémentation OpenCL

OpenCL (pour Open Compute Language) et un langage proche du C permettant de faire facilement du calcul parallèle pour différentes plate-formes d’execution[8]. En particulier, les processeurs classiques (CPU) et les processeurs graphiques (GPU). Un processeur graphique est une puce constituée de centaines, voire milliers de petits processeurs de flux exécutant tous la même tâche.

Les schémas numériques que nous avons utilisés bénéficient fortement d'une implémentation parallèle, car chaque point de la grille peut être calculé indépendamment. C'est pourquoi j'ai choisi d'utiliser OpenCL afin d'exécuter les itérations sur un GPU.

Les résultats obtenus sont les mêmes qu'avec l'implémentation Scilab, mais le temps d'exécution est grandement diminué. En effet, pour l'exemple du déplacement de la gaussienne vu plus haut, on passe de 1min10 sous Scilab à 0.15s sous OpenCL (grille de 30x30). La puce utilisée pour les mesures est un Intel Core i5-3230m, dont on utilise la partie CPU sous Scilab, et la partie graphique Intel HD 4000 sous OpenCL.

La rapidité de mon implémentation OpenCL permet d'augmenter la taille de la grille. La figure ci-dessous montre la gaussienne après son déplacement sur une grille de 256x256, calculé en 45 secondes. Ce calcul aurait pris un peu moins de 8 semaines sous Scilab.

Gaussienne déplacée sur une grille de 256x256


Conclusion

Durant cette IRL, j'ai pu découvrir une nouvelle notion de solution d'EDP et implémenter une méthode de résolution d'un problème très connu. J'ai également pu mettre en pratique mes connaissances en calcul parallèle sur GPU pour calculer très rapidement des solutions de ce problème.


Il subsiste des pistes d'amélioration, on pourrait notamment terminer l'implémentation du schéma monotone, puis le porter également sous OpenCL.

Bibliographie

[1] M. G. Crandall and P.-L. Lions, “Viscosity solutions of hamilton-jacobi equations,” Transactions of the American Mathematical Society, vol. 277, no. 1, pp. 1–42, 1983.

[2] F. Dragoni, “Introduction to viscosity solutions for nonlinear pdes.”

[3] A. Bressan, “Viscosity solutions of hamilton-jacobi equations and optimal control pro- blems,” Lecture notes, 2010.

[4] M. G. Crandall, H. Ishii, and P.-L. Lions, “User’s guide to viscosity solutions of second order partial differential equations,” Bulletin of the American Mathematical Society, vol. 27, no. 1, pp. 1–67, 1992.

[5] M. Bardi and I. Capuzzo-Dolcetta, Optimal control and viscosity solutions of Hamilton- Jacobi-Bellman equations. Springer Science & Business Media, 2008.

[6] J.-D. Benamou, B. D. Froese, and A. M. Oberman, “Numerical solution of the optimal transportation problem using the monge–ampere equation,” Journal of Computational Physics, vol. 260, pp. 107–126, 2014.

[7] C. Villani, Optimal transport : old and new, vol. 338. Springer Science & Business Media, 2008.

[8] J. E. Stone, D. Gohara, and G. Shi, “Opencl : A parallel programming standard for heterogeneous computing systems,” Computing in science & engineering, vol. 12, no. 1- 3, pp. 66–73, 2010.

[9] K. GROUP et al., “The khronos group,” Beaverton, 2011b. Disponível em, 2009.

Documents