doc/tutorial/medcouplingmultiproc_fr.rst

   1
   2 MEDCoupling,  multiprocessing
   3 ----------------------------
   4
   5 Cet exercice fait la supposition de Numpy Scipy sont correctement maîtrisés, sinon voir :ref:`medcouplingnumpyptr`.
   6 On va faire simuler un traitement un peu long (ici l'interpolation d'un maillage de 64000 cells avec un autre de 64000 cells).
   7 On va faire le traitement d'abord en séquentiel puis en parallèle pour exploiter les coeurs de notre CPU.
   8 Nous allons utiliser le module ``multiprocessing`` pour cela.
   9
  10 Début de l'implémentation
  11 ~~~~~~~~~~~~~~~~~~~~~~~~~
  12
  13 Pour commencer l'exercice importer le module Python ``medcoupling``, ``MEDCouplingRemapper``, ``numpy``, ``scipy``, ``multiprocessing``
  14 et ``datetime`` pour chronométrer : ::
  15
  16         import medcoupling as mc
  17         import MEDCouplingRemapper as mr
  18         import multiprocessing as mp
  19         from datetime import datetime
  20         from scipy.sparse import csr_matrix
  21
  22 Créer un maillage cartésien régulier 3D avec 40 cells en X, Y et Z
  23 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  24 Créons un maillage cartésien 3D m de pas régulier entre 0. et 1. en X, Y et Z : ::
  25
  26         nbCells=40
  27         arr=mc.DataArrayDouble(nbCells+1) ; arr.iota() ; arr/=nbCells
  28         m=mc.MEDCouplingCMesh() ; m.setCoords(arr,arr,arr)
  29
  30 Traduisons m en non structuré pour que le calcul soit plus long : ::
  31
  32         m=m.buildUnstructured()
  33
  34 Créer une copie m2 de m translatée de la moitié du pas en X, Y et Z ::
  35
  36         m2=m.deepCopy()
  37         t=mc.DataArrayDouble(3)
  38         t[:]=1/(2*float(nbCells))
  39         m2.translate(t.getValues())
  40
  41 Calculer séquentiellement la matrice d'interpolation M de la projection entre m et m2 en P0P0
  42 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  43
  44 m sera considéré comme le maillage source et m2 sera considéré comme le maillage cible.
  45 Profitons en pour chronométrer le temps necessaire pour le traitement séquentiel.
  46 Utilisons ``MEDCouplingRemapper`` pour cela. ::
  47
  48         remap=mr.MEDCouplingRemapper()
  49         strt=datetime.now()
  50         assert(remap.prepare(m,m2,"P0P0")==1)
  51         print("time in sequential : %s"%(str(datetime.now()-strt)))
  52
  53 Stockons la sparse matrix scipy dans ``matSeq``. ::
  54
  55         matSeq=remap.getCrudeCSRMatrix()
  56
  57 Calculer cette même matrice M en parallèle avec multiprocessing.
  58 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  59
  60 Commencons par récupérer le nombre de coeur de notre machine. ::
  61
  62         nbProc=mp.cpu_count()
  63
  64 L'idée est de faire une méthode ``work`` prenant un tuple de longueur 2.
  65 Le premier élément du tuple est la partie du maillage ``m2`` considérée. Le 2eme élément donne la correspondance entre les cells id de ``m2Part`` les cells id de ``m2``.
  66
  67 L'idée est d'interpoler ``m`` avec ``m2Part``.
  68
  69 On récupèrera ensuite la matrice sparse ``myMat`` issue de ``m`` avec ``m2Part``.
  70 Ensuite l'idée est de générer une matrice sparse ``mat2`` à partir de ``myMat`` avec les ids globaux de ``m2``. ::
  71
  72         def work(inp):
  73             m2Part,partToGlob=inp
  74             myRemap=mr.MEDCouplingRemapper()
  75             assert(myRemap.prepare(m,m2Part,"P0P0")==1)
  76             myMat=myRemap.getCrudeCSRMatrix()
  77             indptrnew=mc.DataArrayInt(m2.getNumberOfCells())
  78             indptrnew.fillWithZero()
  79             d=mc.DataArrayInt(myMat.indptr).deltaShiftIndex()
  80             indptrnew[partToGlob]=d
  81             indptrnew.computeOffsetsFull()
  82             mat2=csr_matrix( (myMat.data,myMat.indices,indptrnew.toNumPyArray()), shape=(m2.getNumberOfCells(),m.getNumberOfCells()))
  83             return mat2
  84
  85 Il s'agit désormais de faire la liste des inputs à donner aux ``nbProc`` instances de ``work`` qui seront exécutés en parallèle.
  86 Appelons cette liste python ``workToDo`` qui sera de longueur ``nbProc``.
  87 On peut se faire aider de ``mc.DataArray.GetSlice``. ::
  88
  89         workToDo=[]
  90         for i in list(range(nbProc)):
  91               s=mc.DataArray.GetSlice(slice(0,m2.getNumberOfCells(),1),i,nbProc)
  92               part=m2[s]
  93               partToGlob=mc.DataArrayInt.Range(s.start,s.stop,s.step)
  94               workToDo.append((part,partToGlob))
  95               pass
  96
  97 On est maintenant prêt pour faire travailler chacun des coeurs indépendamment. Pour ce faire, on crée un ``mp.Pool`` et on assigne à chaque worker le travail ``work`` avec autant de worker que de coeurs. Et chronométrons le tout ! ::
  98
  99
 100         strt=datetime.now()
 101         pool = mp.Pool()
 102         asyncResult = pool.map_async(work,workToDo)
 103         subMatrices = asyncResult.get()
 104         print("time in parallel (x%d) : %s"%(nbProc,str(datetime.now()-strt)))
 105
 106 .. note:: A noter la magie ! On a transféré entre le process maitre et les process esclave sans même s'en rendre compte les maillages et les DataArrayInt contenus dans ``workToDo`` !
 107           Merci à la pickelisation des objets MEDCoupling :)
 108
 109 Vérfication
 110 ~~~~~~~~~~~
 111
 112 Vérifions que les matrices sont les mêmes ! Sommons ``subMatrices`` (``matPar``) et regardons le nombre de non zéros de la différence entre la ``matPar`` et ``matSeq``. ::
 113
 114         matPar = sum(subMatrices)
 115         matDelta=matSeq-matPar
 116         assert(matDelta.nnz==0)