src/ParaMEDMEM_Swig/test_InterpKernelDEC.py

   1 #!/usr/bin/env python
   2 #  -*- coding: iso-8859-1 -*-
   3 # Copyright (C) 2007-2023  CEA, EDF
   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
  30 class ParaMEDMEM_IK_DEC_Tests(unittest.TestCase):
  31     """ See test_StructuredCoincidentDEC_py_1() for a quick start.
  32     """
  33     def generateFullSource(self):
  34         """ The complete source mesh: 4 squares each divided in 2 diagonaly (so 8 cells in total) """
  35         msh  = self.generateFullTarget()
  36         msh.simplexize(0)
  37         msh.setName("src_mesh")
  38         fld = MEDCouplingFieldDouble(ON_CELLS, ONE_TIME)
  39         fld.setMesh(msh); fld.setName("source_F");
  40         da = DataArrayDouble(msh.getNumberOfCells())
  41         da.iota()
  42         da *= 2
  43         fld.setArray(da)
  44         return msh, fld
  45
  46     def generateFullTarget(self):
  47         """ The complete target mesh: 4 squares """
  48         m1 = MEDCouplingCMesh("tgt_msh")
  49         da = DataArrayDouble([0,1,2])
  50         m1.setCoords(da, da)
  51         msh = m1.buildUnstructured()
  52         return msh
  53
  54     #
  55     # Below, the two functions emulating the set up of a piece of the source and target mesh
  56     # on each proc. Obviously in real world problems, this comes from your code and is certainly
  57     # not computed by cuting again from scratch the full-size mesh!!
  58     #
  59     def getPartialSource(self, rank):
  60         """ Will return an empty mesh piece for rank=2 and 3 """
  61         msh, f = self.generateFullSource()
  62         if rank == 0:
  63             sub_m, sub_f = msh[0:4], f[0:4]
  64         elif rank == 1:
  65             sub_m, sub_f = msh[4:8], f[4:8]
  66         sub_m.zipCoords()
  67         return sub_m, sub_f
  68
  69     def getPartialTarget(self, rank):
  70         """ One square for each rank """
  71         msh = self.generateFullTarget()
  72         if rank == 2:
  73             sub_m = msh[[0,2]]
  74         elif rank == 3:
  75             sub_m = msh[[1,3]]
  76         sub_m.zipCoords()
  77         # Receiving side must prepare an empty field that will be filled by DEC:
  78         fld = MEDCouplingFieldDouble(ON_CELLS, ONE_TIME)
  79         da = DataArrayDouble(sub_m.getNumberOfCells())
  80         fld.setArray(da)
  81         fld.setName("tgt_F")
  82         fld.setMesh(sub_m)
  83         return sub_m, fld
  84
  85     def testInterpKernelDEC_ctor(self):
  86         """ Test the various Python ctors """
  87         size = MPI.COMM_WORLD.size
  88         if size != 4:
  89             print("Should be run on 4 procs!")
  90             return
  91         # Define two processor groups
  92         nproc_source = 2
  93         l1, l2 = range(nproc_source), range(size - nproc_source, size)
  94         # With 2 iterables:
  95         i1 = InterpKernelDEC.New(l1, l2)
  96         # Should also work directly:
  97         i2 = InterpKernelDEC(l1, l2)
  98         # With 2 proc groups:
  99         interface = CommInterface()
 100         source_group = MPIProcessorGroup(interface, list(l1))
 101         target_group = MPIProcessorGroup(interface, list(l2))
 102         i3 = InterpKernelDEC.New(source_group, target_group)
 103         # Should also work directly:
 104         i4 = InterpKernelDEC(source_group, target_group)
 105         # With 2 iterables and a custom comm:
 106         i5 = InterpKernelDEC.New(l1, l2, MPI.COMM_WORLD)
 107         # Work directly with the **hack**
 108         i6 = InterpKernelDEC(l1, l2, MPI._addressof(MPI.COMM_WORLD))
 109         # Should fail with 2 proc groups **and** a communicator
 110         self.assertRaises(InterpKernelException, InterpKernelDEC.New, source_group, target_group, MPI.COMM_WORLD)
 111         self.assertRaises(Exception, InterpKernelDEC, source_group, target_group, MPI.COMM_WORLD)
 112         i6.release(); i5.release(); i4.release(); i3.release(); i2.release(); i1.release()
 113         source_group.release()
 114         target_group.release()
 115
 116     @WriteInTmpDir
 117     def testInterpKernelDEC_2D_py_1(self):
 118         """ This test illustrates a basic use of the InterpKernelDEC.
 119         Look at the C++ documentation of the class for more informations.
 120         """
 121         size = MPI.COMM_WORLD.size
 122         rank = MPI.COMM_WORLD.rank
 123         if size != 4:
 124             print("Should be run on 4 procs!")
 125             return
 126
 127         # Define two processor groups
 128         nproc_source = 2
 129         procs_source = list(range(nproc_source))
 130         procs_target = list(range(size - nproc_source, size))
 131
 132         interface = CommInterface()
 133         source_group = MPIProcessorGroup(interface, procs_source)
 134         target_group = MPIProcessorGroup(interface, procs_target)
 135         idec = InterpKernelDEC(source_group, target_group)
 136
 137         # Write out full size meshes/fields for inspection
 138         if rank == 0:
 139             _, fld = self.generateFullSource()
 140             mshT = self.generateFullTarget()
 141             WriteField("./source_field_FULL.med", fld, True)
 142             WriteUMesh("./target_mesh_FULL.med", mshT, True)
 143
 144         MPI.COMM_WORLD.Barrier()  # really necessary??
 145
 146         #
 147         # OK, let's go DEC !!
 148         #
 149         if source_group.containsMyRank():
 150             _, fieldS = self.getPartialSource(rank)
 151             fieldS.setNature(IntensiveMaximum)   # The only policy supported for now ...
 152             WriteField("./source_field_part_%d.med" % rank, fieldS, True)
 153             idec.attachLocalField(fieldS)
 154             idec.synchronize()
 155             idec.sendData()
 156
 157         if target_group.containsMyRank():
 158             mshT, fieldT = self.getPartialTarget(rank)
 159             fieldT.setNature(IntensiveMaximum)
 160             WriteUMesh("./target_mesh_part_%d.med" % rank, mshT, True)
 161             idec.attachLocalField(fieldT)
 162             idec.synchronize()
 163             idec.recvData()
 164             # Now the actual checks:
 165             if rank == 2:
 166                 self.assertEqual(fieldT.getArray().getValues(), [1.0, 9.0])
 167             elif rank == 3:
 168                 self.assertEqual(fieldT.getArray().getValues(), [5.0, 13.0])
 169
 170         # Release DEC (this involves MPI exchanges -- notably the release of the communicator -- so better be done before MPI.Finalize()
 171         idec.release()
 172         source_group.release()
 173         target_group.release()
 174         MPI.COMM_WORLD.Barrier()
 175
 176     @WriteInTmpDir
 177     def test_InterpKernelDEC_2D_py_2(self):
 178         """ More involved test using Para* objects.
 179         """
 180         size = MPI.COMM_WORLD.size
 181         rank = MPI.COMM_WORLD.rank
 182         if size != 5:
 183             print("Should be run on 5 procs!")
 184             return
 185
 186         print(rank)
 187         nproc_source = 3
 188         procs_source = list(range(nproc_source))
 189         procs_target = list(range(size - nproc_source + 1, size))
 190
 191         interface = CommInterface()
 192         target_group = MPIProcessorGroup(interface, procs_target)
 193         source_group = MPIProcessorGroup(interface, procs_source)
 194         dec = InterpKernelDEC(source_group, target_group)
 195
 196         data_dir = os.path.join(os.environ['MEDCOUPLING_ROOT_DIR'], "share", "resources", "med")
 197         if not os.path.isdir(data_dir):
 198             data_dir = os.environ.get('MED_RESOURCES_DIR',"::").split(":")[1]
 199
 200         filename_xml1 = os.path.join(data_dir, "square1_split")
 201         filename_xml2 = os.path.join(data_dir, "square2_split")
 202
 203         MPI.COMM_WORLD.Barrier()
 204         if source_group.containsMyRank():
 205             filename = filename_xml1 + str(rank+1) + ".med"
 206             meshname = "Mesh_2_" + str(rank+1)
 207             mesh=ReadUMeshFromFile(filename,meshname,0)
 208             paramesh=ParaMESH(mesh,source_group,"source mesh")
 209             comptopo = ComponentTopology()
 210             parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo)
 211             parafield.getField().setNature(IntensiveMaximum)
 212             nb_local=mesh.getNumberOfCells()
 213             value = [1.0]*nb_local
 214             parafield.getField().setValues(value)
 215             icocofield = ICoCoMEDDoubleField(parafield.getField())
 216             dec.attachLocalField(icocofield)
 217             pass
 218         else:
 219             filename = filename_xml2 + str(rank - nproc_source + 1) + ".med"
 220             meshname = "Mesh_3_" + str(rank - nproc_source + 1)
 221             mesh=ReadUMeshFromFile(filename,meshname,0)
 222             paramesh=ParaMESH(mesh,target_group,"target mesh")
 223             comptopo = ComponentTopology()
 224             parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo)
 225             parafield.getField().setNature(IntensiveMaximum)
 226             nb_local=mesh.getNumberOfCells()
 227             value = [0.0]*nb_local
 228             parafield.getField().setValues(value)
 229             icocofield = ICoCoMEDDoubleField(parafield.getField())
 230             dec.attachLocalField(icocofield)
 231             pass
 232
 233         if source_group.containsMyRank():
 234             field_before_int = parafield.getVolumeIntegral(0,True)
 235             dec.synchronize()
 236             dec.setForcedRenormalization(False)
 237             dec.sendData()
 238             dec.recvData()
 239             field_after_int=parafield.getVolumeIntegral(0,True);
 240             self.assertTrue(math.fabs(field_after_int-field_before_int)<1e-8)
 241             pass
 242         else:
 243             dec.synchronize()
 244             dec.setForcedRenormalization(False)
 245             dec.recvData()
 246             dec.sendData()
 247             pass
 248         ## end
 249
 250         # Some clean up that still needs MPI communication, so to be done before MPI_Finalize()
 251         parafield.release()
 252         paramesh.release()
 253         dec.release()
 254         target_group.release()
 255         source_group.release()
 256         MPI.COMM_WORLD.Barrier()
 257
 258     def test_InterpKernelDEC_default(self):
 259         """
 260         [EDF27375] : Put a default value when non intersecting case
 261         """
 262         size = MPI.COMM_WORLD.size
 263         rank = MPI.COMM_WORLD.rank
 264         if size != 4:
 265             print("Should be run on 4 procs!")
 266             return
 267         nproc_source = 2
 268         procs_source = list(range(nproc_source))
 269         procs_target = list(range(size - nproc_source, size))
 270
 271         interface = CommInterface()
 272         target_group = MPIProcessorGroup(interface, procs_target)
 273         source_group = MPIProcessorGroup(interface, procs_source)
 274         dec = InterpKernelDEC(source_group, target_group)
 275         #
 276         MPI.COMM_WORLD.Barrier()
 277         if source_group.containsMyRank():
 278             mesh = eval("Source_Proc_{}".format(rank))()
 279             nb_local=mesh.getNumberOfCells()
 280             field = MEDCouplingFieldDouble(ON_CELLS)
 281             field.setNature(IntensiveMaximum)
 282             field.setMesh(mesh)
 283             arr = DataArrayDouble(nb_local) ; arr.iota() ; arr += rank
 284             field.setArray(arr)
 285             dec.attachLocalField(field)
 286             dec.synchronizeWithDefaultValue(-2000.0)
 287             dec.sendData()
 288             # target -> source
 289             dec.recvData()
 290             if rank == 0:
 291                 self.assertTrue(field.getArray().isEqual(DataArrayDouble([0.6,0.6,-2000]),1e-12))
 292                 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([2])) )
 293             if rank == 1:
 294                 self.assertTrue(field.getArray().isEqual(DataArrayDouble([1.0,-2000]),1e-12))
 295                 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([1])) )
 296         else:
 297             mesh = eval("Target_Proc_{}".format(rank))()
 298             nb_local=mesh.getNumberOfCells()
 299             field = MEDCouplingFieldDouble(ON_CELLS)
 300             field.setNature(IntensiveMaximum)
 301             field.setMesh(mesh)
 302             arr = DataArrayDouble(nb_local) ; arr[:] = -20
 303             field.setArray(arr)
 304             dec.attachLocalField(field)
 305             dec.synchronizeWithDefaultValue(-1000.0)
 306             dec.recvData()
 307             if rank == 2:
 308                 # matrix S0 / T2 = [[(0,S0,1),(1,S0,1.5)]
 309                 # IntensiveMaximum => [[(0,S0,1/2.5),(1,S0,1.5/2.5)]
 310                 #
 311                 self.assertTrue(field.getArray().isEqual(DataArrayDouble([0.6]),1e-12))
 312                 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([])) )
 313             if rank == 3:
 314                 # matrix S1 / T3 = [[],[(0,S1,1.0)],[(0,S1,2.0)],[]]
 315                 # IntensiveMaximum => [[],[(0,S1,1.0/1.0)],[(0,S1,2.0/2.0)],[]]
 316                 self.assertTrue(field.getArray().isEqual(DataArrayDouble([-1000.0, 1.0, 1.0, -1000.0]),1e-8))
 317                 self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([0,3])) )
 318             # target -> source
 319             dec.sendData()
 320
 321         # Some clean up that still needs MPI communication, so to be done before MPI_Finalize()
 322         dec.release()
 323         target_group.release()
 324         source_group.release()
 325         MPI.COMM_WORLD.Barrier()
 326
 327 def Source_Proc_0():
 328     coo = DataArrayDouble([(0,2),(2,2),(4,2),(0,4),(2,4),(4,4),(0,6),(2,6)])
 329     m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
 330     m.insertNextCell(NORM_QUAD4,[0,1,4,3])
 331     m.insertNextCell(NORM_QUAD4,[1,2,5,4])
 332     m.insertNextCell(NORM_QUAD4,[3,4,7,6])
 333     return m
 334
 335 def Source_Proc_1():
 336     coo = DataArrayDouble([(6,2),(8,2),(10,2),(6,4),(8,4),(10,4)])
 337     m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
 338     m.insertNextCell(NORM_QUAD4,[0,1,4,3])
 339     m.insertNextCell(NORM_QUAD4,[1,2,5,4])
 340     return m
 341
 342 def Target_Proc_2():
 343     coo = DataArrayDouble([(1,0),(3.5,0),(1,3),(3.5,3)])
 344     m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
 345     m.insertNextCell(NORM_QUAD4,[0,1,3,2])
 346     return m
 347
 348 def Target_Proc_3():
 349     coo = DataArrayDouble([(6,0),(7,0),(8,0),(9,0),(10,0),
 350                            (6,1),(7,1),(9,1),(10,1),
 351                            (7,3),(8,3),
 352                            (6,4),(7,4)])
 353     m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
 354     m.insertNextCell(NORM_QUAD4,[0,1,6,5])
 355     m.insertNextCell(NORM_QUAD4,[1,2,10,9])
 356     m.insertNextCell(NORM_QUAD4,[5,6,12,11])
 357     m.insertNextCell(NORM_QUAD4,[3,4,8,7])
 358     return m
 359
 360 if __name__ == "__main__":
 361     unittest.main()
 362     MPI.Finalize()
 363