src/ParaMEDMEM_Swig/test_OverlapDEC.py

   1 #!/usr/bin/env python
   2 #  -*- coding: iso-8859-1 -*-
   3 # Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License, or (at your option) any later version.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21
  22 from medcoupling import *
  23 from ParaMEDMEMTestTools import WriteInTmpDir
  24 import sys, os
  25 import unittest
  26 import math
  27 from mpi4py import MPI
  28
  29 class ParaMEDMEM_O_DEC_Tests(unittest.TestCase):
  30     """ This test illustrates a basic use of the OverlapDEC and shows notably that not all
  31     processors must possess a piece of the source and/or target mesh.
  32     Look at the C++ documentation of the class for more informations.
  33     In this case, the source mesh is only stored on 2 procs, whereas the target is on 4.
  34     Since only a single group of processor is defined in the setup, the 2 idle procs on the source side are just providing an empty mesh,
  35     thus indicating that they don't participate in the source definition.
  36
  37     Main method is testOverlapDEC_2D_py_1()
  38     """
  39
  40     def generateFullSource(self):
  41         """ The complete source mesh: 4 squares each divided in 2 diagonaly (so 8 cells in total) """
  42         msh  = self.generateFullTarget()
  43         msh.simplexize(0)
  44         msh.setName("src_mesh")
  45         fld = MEDCouplingFieldDouble(ON_CELLS, ONE_TIME)
  46         fld.setMesh(msh); fld.setName("source_F");
  47         da = DataArrayDouble(msh.getNumberOfCells())
  48         da.iota()
  49         da *= 2
  50         fld.setArray(da)
  51         return msh, fld
  52
  53     def generateFullTarget(self):
  54         """ The complete target mesh: 4 squares """
  55         m1 = MEDCouplingCMesh("tgt_msh")
  56         da = DataArrayDouble([0,1,2])
  57         m1.setCoords(da, da)
  58         msh = m1.buildUnstructured()
  59         return msh
  60
  61     #
  62     # Below, the two functions emulating the set up of a piece of the source and target mesh
  63     # on each proc. Obviously in real world problems, this comes from your code and is certainly
  64     # not computed by cuting again from scratch the full-size mesh!!
  65     #
  66     def getPartialSource(self, rank):
  67         """ Will return an empty mesh piece for rank=2 and 3 """
  68         msh, f = self.generateFullSource()
  69         if rank in [2,3]:
  70             sub_m, sub_f = msh[[]], f[[]]  # Little trick to select nothing in the mesh, thus producing an empty mesh
  71         elif rank == 0:
  72             sub_m, sub_f = msh[0:4], f[0:4]
  73         elif rank == 1:
  74             sub_m, sub_f = msh[4:8], f[4:8]
  75         sub_m.zipCoords()
  76         return sub_m, sub_f
  77
  78     def getPartialTarget(self, rank):
  79         """ One square for each rank """
  80         msh = self.generateFullTarget()
  81         sub_m = msh[rank]
  82         sub_m.zipCoords()
  83         # Receiving side must prepare an empty field that will be filled by DEC:
  84         fld = MEDCouplingFieldDouble(ON_CELLS, ONE_TIME)
  85         da = DataArrayDouble(sub_m.getNumberOfCells())
  86         fld.setArray(da)
  87         fld.setName("tgt_F")
  88         fld.setMesh(sub_m)
  89         return sub_m, fld
  90
  91     @WriteInTmpDir
  92     def testOverlapDEC_2D_py_1(self):
  93         """ The main method of the test """
  94         size = MPI.COMM_WORLD.size
  95         rank = MPI.COMM_WORLD.rank
  96         if size != 4:
  97             raise RuntimeError("Should be run on 4 procs!")
  98
  99         # Define (single) processor group - note the difference with InterpKernelDEC which needs two groups.
 100         proc_group = list(range(size))   # No need for ProcessorGroup object here.
 101         odec = OverlapDEC(proc_group)
 102
 103         # Write out full size meshes/fields for inspection
 104         if rank == 0:
 105             _, fld = self.generateFullSource()
 106             mshT = self.generateFullTarget()
 107             WriteField("./source_field_FULL.med", fld, True)
 108             WriteUMesh("./target_mesh_FULL.med", mshT, True)
 109
 110         MPI.COMM_WORLD.Barrier()  # really necessary??
 111
 112         #
 113         # OK, let's go DEC !!
 114         #
 115         _, fieldS = self.getPartialSource(rank)
 116         fieldS.setNature(IntensiveMaximum)   # The only policy supported for now ...
 117         mshT, fieldT = self.getPartialTarget(rank)
 118         fieldT.setNature(IntensiveMaximum)
 119         if rank not in [2,3]:
 120             WriteField("./source_field_part_%d.med" % rank, fieldS, True)
 121         WriteUMesh("./target_mesh_part_%d.med" % rank, mshT, True)
 122
 123         odec.attachSourceLocalField(fieldS)
 124         odec.attachTargetLocalField(fieldT)
 125         odec.synchronize()
 126         odec.sendRecvData()
 127
 128         # Now the actual checks:
 129         if rank == 0:
 130             self.assertEqual(fieldT.getArray().getValues(), [1.0])
 131         elif rank == 1:
 132             self.assertEqual(fieldT.getArray().getValues(), [5.0])
 133         elif rank == 2:
 134             self.assertEqual(fieldT.getArray().getValues(), [9.0])
 135         elif rank == 3:
 136             self.assertEqual(fieldT.getArray().getValues(), [13.0])
 137
 138         # Release DEC (this involves MPI exchanges -- notably the release of the communicator -- so better be done before MPI.Finalize()
 139         odec.release()
 140
 141         MPI.COMM_WORLD.Barrier()
 142
 143 if __name__ == "__main__":
 144     unittest.main()
 145     MPI.Finalize()