]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/ParaMEDMEM_Swig/test_InterpKernelDEC.py
Salome HOME
Merge branch 'V9_11_BR'
[tools/medcoupling.git] / 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 if __name__ == "__main__":
259     unittest.main()
260     MPI.Finalize()
261