Salome HOME
One more exercise in tuto
authorAnthony Geay <anthony.geay@edf.fr>
Wed, 9 Nov 2016 16:45:52 +0000 (17:45 +0100)
committerAnthony Geay <anthony.geay@edf.fr>
Wed, 9 Nov 2016 16:45:52 +0000 (17:45 +0100)
doc/tutorial/index.rst
doc/tutorial/medcouplingmultiproc_fr.rst [new file with mode: 0644]
doc/tutorial/medcouplingnumpy_fr.rst

index d93c549501600cf42fcb0b205e90d6bfd7cff735..7cba76275b52b8853a5bf58291222ae9bd96923b 100644 (file)
@@ -26,6 +26,7 @@ Table des matières:
    medloader_fr
    medcouplingremapper_fr
    medcouplingnumpy_fr
+   medcouplingmultiproc_fr
    medcouplingcorba_fr
    medcouplingloaderex1_fr
    medcouplingloaderex2_fr
diff --git a/doc/tutorial/medcouplingmultiproc_fr.rst b/doc/tutorial/medcouplingmultiproc_fr.rst
new file mode 100644 (file)
index 0000000..7e0cd3d
--- /dev/null
@@ -0,0 +1,117 @@
+
+MEDCoupling,  multiprocessing
+----------------------------
+
+Cet exercice fait la supposition de Numpy Scipy sont correctement maîtrisés, sinon voir :ref:`medcouplingnumpyptr`.
+On va faire simuler un traitement un peu long (ici l'interpolation d'un maillage de 64000 cells avec un autre de 64000 cells).
+On va faire le traitement d'abord en séquentiel puis en parallèle pour exploiter les coeurs de notre CPU.
+Nous allons utiliser le module ``multiprocessing`` pour cela.
+
+Début de l'implémentation
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Pour commencer l'exercice importer le module Python ``MEDCoupling``, ``MEDCouplingRemapper``, ``numpy``, ``scipy``, ``multiprocessing``
+et ``datetime`` pour chronométrer : ::
+
+       import MEDCoupling as mc
+       import MEDCouplingRemapper as mr
+       import multiprocessing as mp
+       from datetime import datetime
+       from scipy.sparse import csr_matrix
+
+Créer un maillage cartésien régulier 3D avec 40 cells en X, Y et Z
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Créons un maillage cartésien 3D m de pas régulier entre 0. et 1. en X, Y et Z : ::
+
+       nbCells=40
+       arr=mc.DataArrayDouble(nbCells+1) ; arr.iota() ; arr/=nbCells
+       m=mc.MEDCouplingCMesh() ; m.setCoords(arr,arr,arr)
+
+Traduisons m en non structuré pour que le calcul soit plus long : ::
+
+       m=m.buildUnstructured()
+
+Créer une copie m2 de m translatée de la moitié du pas en X, Y et Z ::
+
+       m2=m.deepCopy()
+       t=mc.DataArrayDouble(3)
+       t[:]=1/(2*float(nbCells))
+       m2.translate(t.getValues())
+
+Calculer séquentiellement la matrice d'interpolation M de la projection entre m et m2 en P0P0
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+m sera considéré comme le maillage source et m2 sera considéré comme le maillage cible.
+Profitons en pour chronométrer le temps necessaire pour le traitement séquentiel.
+Utilisons ``MEDCouplingRemapper`` pour cela. ::
+
+       remap=mr.MEDCouplingRemapper()
+       strt=datetime.now()
+       assert(remap.prepare(m,m2,"P0P0")==1)
+       print "time in sequential : %s"%(str(datetime.now()-strt))
+
+Stockons la sparse matrix scipy dans ``matSeq``. ::
+
+       matSeq=remap.getCrudeCSRMatrix()
+
+Calculer cette même matrice M en parallèle avec multiprocessing.
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Commencons par récupérer le nombre de coeur de notre machine. ::
+
+       nbProc=mp.cpu_count()
+
+L'idée est de faire une méthode ``work`` prenant un tuple de longueur 2.
+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``.
+
+L'idée est d'interpoler ``m`` avec ``m2Part``.
+
+On récupèrera ensuite la matrice sparse ``myMat`` issue de ``m`` avec ``m2Part``.
+Ensuite l'idée et de générer une matrice sparse ``mat2`` à partir de ``myMat`` avec les ids globaux de ``m2``. ::
+
+       def work(inp):
+            m2Part,partToGlob=inp
+           myRemap=mr.MEDCouplingRemapper()
+           assert(myRemap.prepare(m,m2Part,"P0P0")==1)
+           myMat=myRemap.getCrudeCSRMatrix()
+           a=mc.DataArrayInt.Range(s.start,s.stop,s.step)
+           indptrnew=mc.DataArrayInt(m2.getNumberOfCells())
+           indptrnew.fillWithZero()
+           d=mc.DataArrayInt(myMat.indptr).deltaShiftIndex()
+           indptrnew[partToGlob]=d
+           indptrnew.computeOffsetsFull()
+           mat2=csr_matrix( (myMat.data,myMat.indices,indptrnew.toNumPyArray()), shape=(m2.getNumberOfCells(),m.getNumberOfCells()))
+           return mat2
+
+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.
+Appelons cette liste python ``workToDo`` qui sera de longueur ``nbProc``.
+On peut se faire aider de ``mc.DataArray.GetSlice``. ::
+
+        workToDo=[]
+        for i in xrange(nbProc):
+              s=mc.DataArray.GetSlice(slice(0,m2.getNumberOfCells(),1),i,nbProc)
+              part=m2[s]
+              partToGlob=mc.DataArrayInt.Range(s.start,s.stop,s.step)
+              workToDo.append((part,partToGlob))
+              pass
+
+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 ! ::
+
+
+       strt=datetime.now()
+       pool = mp.Pool()
+       asyncResult = pool.map_async(work,workToDo)
+       subMatrices = asyncResult.get()
+       print "time in parallel (x%d) : %s"%(nbProc,str(datetime.now()-strt))
+
+.. 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`` !
+         Merci à la pickelisation des objets MEDCoupling :)
+
+Vérfication
+~~~~~~~~~~~
+
+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``. ::
+
+       matPar = sum(subMatrices)
+       matDelta=matSeq-matPar
+       assert(matDelta.nnz==0)
index ccc2911b4cfe46f9bac1ddc661f48d60bfa5e9ff..d41055ef6dbe05be678da869a0accd4464467c09 100644 (file)
@@ -1,3 +1,4 @@
+.. _medcouplingnumpyptr:
 
 MEDCoupling,  NumPy et SciPy
 ----------------------------