From: Anthony Geay Date: Wed, 9 Nov 2016 16:45:52 +0000 (+0100) Subject: One more exercise in tuto X-Git-Tag: V8_2_0b1~2 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=f2a31df2f6d7e59e6461238614abf3909c2f9466;p=tools%2Fmedcoupling.git One more exercise in tuto --- diff --git a/doc/tutorial/index.rst b/doc/tutorial/index.rst index d93c54950..7cba76275 100644 --- a/doc/tutorial/index.rst +++ b/doc/tutorial/index.rst @@ -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 index 000000000..7e0cd3dfa --- /dev/null +++ b/doc/tutorial/medcouplingmultiproc_fr.rst @@ -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) diff --git a/doc/tutorial/medcouplingnumpy_fr.rst b/doc/tutorial/medcouplingnumpy_fr.rst index ccc2911b4..d41055ef6 100644 --- a/doc/tutorial/medcouplingnumpy_fr.rst +++ b/doc/tutorial/medcouplingnumpy_fr.rst @@ -1,3 +1,4 @@ +.. _medcouplingnumpyptr: MEDCoupling, NumPy et SciPy ----------------------------