xref: /petsc/src/binding/petsc4py/test/test_dmplex.py (revision 8af18dd833b18fc2df02377e75d76d273dc343a9) 
1import petsc4py
2from petsc4py import PETSc
3import unittest
4import os
5import filecmp
6import numpy as np
7
8# --------------------------------------------------------------------
9
10ERR_ARG_OUTOFRANGE = 63
11
12class BaseTestPlex(object):
13
14    COMM = PETSc.COMM_WORLD
15    DIM = 1
16    CELLS = [[0, 1], [1, 2]]
17    COORDS = [[0.], [0.5], [1.]]
18    COMP = 1
19    DOFS = [1, 0]
20
21    def setUp(self):
22        self.plex = PETSc.DMPlex().createFromCellList(self.DIM,
23                                                      self.CELLS,
24                                                      self.COORDS,
25                                                      comm=self.COMM)
26
27    def tearDown(self):
28        self.plex.destroy()
29        self.plex = None
30
31    def testTopology(self):
32        rank = self.COMM.rank
33        dim = self.plex.getDimension()
34        pStart, pEnd = self.plex.getChart()
35        cStart, cEnd = self.plex.getHeightStratum(0)
36        vStart, vEnd = self.plex.getDepthStratum(0)
37        numDepths = self.plex.getLabelSize("depth")
38        coords_raw = self.plex.getCoordinates().getArray()
39        coords = np.reshape(coords_raw, (vEnd - vStart, dim))
40        self.assertEqual(dim, self.DIM)
41        self.assertEqual(numDepths, self.DIM+1)
42        if rank == 0 and self.CELLS is not None:
43            self.assertEqual(cEnd-cStart, len(self.CELLS))
44        if rank == 0 and self.COORDS is not None:
45            self.assertEqual(vEnd-vStart, len(self.COORDS))
46            self.assertTrue((coords == self.COORDS).all())
47
48    def testClosure(self):
49        pStart, pEnd = self.plex.getChart()
50        for p in range(pStart, pEnd):
51            closure = self.plex.getTransitiveClosure(p)[0]
52            for c in closure:
53                cone = self.plex.getCone(c)
54                self.assertEqual(self.plex.getConeSize(c), len(cone))
55                for i in cone:
56                    self.assertIn(i, closure)
57            star = self.plex.getTransitiveClosure(p, useCone=False)[0]
58            for s in star:
59                support = self.plex.getSupport(s)
60                self.assertEqual(self.plex.getSupportSize(s), len(support))
61                for i in support:
62                    self.assertIn(i, star)
63
64    def testAdjacency(self):
65        PETSc.DMPlex.setAdjacencyUseAnchors(self.plex, False)
66        flag = PETSc.DMPlex.getAdjacencyUseAnchors(self.plex)
67        self.assertFalse(flag)
68        PETSc.DMPlex.setAdjacencyUseAnchors(self.plex, True)
69        flag = PETSc.DMPlex.getAdjacencyUseAnchors(self.plex)
70        self.assertTrue(flag)
71        PETSc.DMPlex.setBasicAdjacency(self.plex, False, False)
72        flagA, flagB = PETSc.DMPlex.getBasicAdjacency(self.plex)
73        self.assertFalse(flagA)
74        self.assertFalse(flagB)
75        PETSc.DMPlex.setBasicAdjacency(self.plex, True, True)
76        flagA, flagB = PETSc.DMPlex.getBasicAdjacency(self.plex)
77        self.assertTrue(flagA)
78        self.assertTrue(flagB)
79        pStart, pEnd = self.plex.getChart()
80        for p in range(pStart, pEnd):
81            adjacency = self.plex.getAdjacency(p)
82            self.assertTrue(p in adjacency)
83            self.assertTrue(len(adjacency) > 1)
84
85    def testSectionDofs(self):
86        self.plex.setNumFields(1)
87        section = self.plex.createSection([self.COMP], [self.DOFS])
88        size = section.getStorageSize()
89        entity_dofs = [self.plex.getStratumSize("depth", d) *
90                       self.DOFS[d] for d in range(self.DIM+1)]
91        self.assertEqual(sum(entity_dofs), size)
92
93    def testSectionClosure(self):
94        section = self.plex.createSection([self.COMP], [self.DOFS])
95        self.plex.setSection(section)
96        vec = self.plex.createLocalVec()
97        pStart, pEnd = self.plex.getChart()
98        for p in range(pStart, pEnd):
99            for i in range(section.getDof(p)):
100                off = section.getOffset(p)
101                vec.setValue(off+i, p)
102
103        for p in range(pStart, pEnd):
104            point_closure = self.plex.getTransitiveClosure(p)[0]
105            dof_closure = self.plex.vecGetClosure(section, vec, p)
106            for p in dof_closure:
107                self.assertIn(p, point_closure)
108
109    def testBoundaryLabel(self):
110        pStart, pEnd = self.plex.getChart()
111        if (pEnd - pStart == 0): return
112
113        self.assertFalse(self.plex.hasLabel("boundary"))
114        self.plex.markBoundaryFaces("boundary")
115        self.assertTrue(self.plex.hasLabel("boundary"))
116
117        faces = self.plex.getStratumIS("boundary", 1)
118        for f in faces.getIndices():
119            points, orient = self.plex.getTransitiveClosure(f, useCone=True)
120            for p in points:
121                self.plex.setLabelValue("boundary", p, 1)
122
123        for p in range(pStart, pEnd):
124            if self.plex.getLabelValue("boundary", p) != 1:
125                self.plex.setLabelValue("boundary", p, 2)
126
127        numBoundary = self.plex.getStratumSize("boundary", 1)
128        numInterior = self.plex.getStratumSize("boundary", 2)
129        self.assertNotEqual(numBoundary, pEnd - pStart)
130        self.assertNotEqual(numInterior, pEnd - pStart)
131        self.assertEqual(numBoundary + numInterior, pEnd - pStart)
132
133    def testMetric(self):
134        if self.DIM == 1: return
135        self.plex.distribute()
136        if self.CELLS is None and not self.plex.isSimplex(): return
137
138        h_min = 1.0e-30
139        h_max = 1.0e+30
140        a_max = 1.0e+10
141        target = 10.0
142        p = 1.0
143        beta = 1.3
144        self.plex.metricSetIsotropic(False)
145        self.plex.metricSetRestrictAnisotropyFirst(False)
146        self.plex.metricSetNoInsertion(False)
147        self.plex.metricSetNoSwapping(False)
148        self.plex.metricSetNoMovement(False)
149        self.plex.metricSetVerbosity(-1)
150        self.plex.metricSetNumIterations(3)
151        self.plex.metricSetMinimumMagnitude(h_min)
152        self.plex.metricSetMaximumMagnitude(h_max)
153        self.plex.metricSetMaximumAnisotropy(a_max)
154        self.plex.metricSetTargetComplexity(target)
155        self.plex.metricSetNormalizationOrder(p)
156        self.plex.metricSetGradationFactor(beta)
157
158        self.assertFalse(self.plex.metricIsIsotropic())
159        self.assertFalse(self.plex.metricRestrictAnisotropyFirst())
160        self.assertFalse(self.plex.metricNoInsertion())
161        self.assertFalse(self.plex.metricNoSwapping())
162        self.assertFalse(self.plex.metricNoMovement())
163        assert self.plex.metricGetVerbosity() == -1
164        assert self.plex.metricGetNumIterations() == 3
165        assert np.isclose(self.plex.metricGetMinimumMagnitude(), h_min)
166        assert np.isclose(self.plex.metricGetMaximumMagnitude(), h_max)
167        assert np.isclose(self.plex.metricGetMaximumAnisotropy(), a_max)
168        assert np.isclose(self.plex.metricGetTargetComplexity(), target)
169        assert np.isclose(self.plex.metricGetNormalizationOrder(), p)
170        assert np.isclose(self.plex.metricGetGradationFactor(), beta)
171
172        metric1 = self.plex.metricCreateUniform(1.0)
173        metric2 = self.plex.metricCreateUniform(2.0)
174        metric = self.plex.metricAverage2(metric1, metric2)
175        metric2.array[:] *= 0.75
176        assert np.allclose(metric.array, metric2.array)
177        metric = self.plex.metricIntersection2(metric1, metric2)
178        assert np.allclose(metric.array, metric1.array)
179        self.plex.metricEnforceSPD(metric)
180        assert np.allclose(metric.array, metric1.array)
181
182    def testAdapt(self):
183        if self.DIM == 1: return
184        self.plex.distribute()
185        if self.CELLS is None and not self.plex.isSimplex(): return
186        if sum(self.DOFS) > 1: return
187        metric = self.plex.metricCreateUniform(9.0)
188        try:
189            newplex = self.plex.adaptMetric(metric,"")
190        except PETSc.Error as exc:
191            if exc.ierr != ERR_ARG_OUTOFRANGE: raise
192
193
194# --------------------------------------------------------------------
195
196class BaseTestPlex_2D(BaseTestPlex):
197    DIM = 2
198    CELLS = [[0, 1, 3], [1, 3, 4], [1, 2, 4], [2, 4, 5],
199             [3, 4, 6], [4, 6, 7], [4, 5, 7], [5, 7, 8]]
200    COORDS = [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
201              [0.0, 0.5], [0.5, 0.5], [1.0, 0.5],
202              [0.0, 1.0], [0.5, 1.0], [1.0, 1.0]]
203    DOFS = [1, 0, 0]
204
205class BaseTestPlex_3D(BaseTestPlex):
206    DIM = 3
207    CELLS = [[0, 2, 3, 7], [0, 2, 6, 7], [0, 4, 6, 7],
208             [0, 1, 3, 7], [0, 1, 5, 7], [0, 4, 5, 7]]
209    COORDS = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [1., 1., 0.],
210              [0., 0., 1.], [1., 0., 1.], [0., 1., 1.], [1., 1., 1.]]
211    DOFS = [1, 0, 0, 0]
212
213# --------------------------------------------------------------------
214
215class TestPlex_1D(BaseTestPlex, unittest.TestCase):
216    pass
217
218class TestPlex_2D(BaseTestPlex_2D, unittest.TestCase):
219    pass
220
221class TestPlex_3D(BaseTestPlex_3D, unittest.TestCase):
222    pass
223
224class TestPlex_2D_P3(BaseTestPlex_2D, unittest.TestCase):
225    DOFS = [1, 2, 1]
226
227class TestPlex_3D_P3(BaseTestPlex_3D, unittest.TestCase):
228    DOFS = [1, 2, 1, 0]
229
230class TestPlex_3D_P4(BaseTestPlex_3D, unittest.TestCase):
231    DOFS = [1, 3, 3, 1]
232
233class TestPlex_2D_BoxTensor(BaseTestPlex_2D, unittest.TestCase):
234    CELLS = None
235    COORDS = None
236    def setUp(self):
237        self.plex = PETSc.DMPlex().createBoxMesh([3,3], simplex=False)
238
239class TestPlex_3D_BoxTensor(BaseTestPlex_3D, unittest.TestCase):
240    CELLS = None
241    COORDS = None
242    def setUp(self):
243        self.plex = PETSc.DMPlex().createBoxMesh([3,3,3], simplex=False)
244
245try:
246    raise PETSc.Error
247    PETSc.DMPlex().createBoxMesh([2,2], simplex=True, comm=PETSc.COMM_SELF).destroy()
248except PETSc.Error:
249    pass
250else:
251    class TestPlex_2D_Box(BaseTestPlex_2D, unittest.TestCase):
252        CELLS = None
253        COORDS = None
254        def setUp(self):
255            self.plex = PETSc.DMPlex().createBoxMesh([1,1], simplex=True)
256
257    class TestPlex_2D_Boundary(BaseTestPlex_2D, unittest.TestCase):
258        CELLS = None
259        COORDS = None
260        def setUp(self):
261            boundary = PETSc.DMPlex().create(self.COMM)
262            boundary.createSquareBoundary([0., 0.], [1., 1.], [2, 2])
263            boundary.setDimension(self.DIM-1)
264            self.plex = PETSc.DMPlex().generate(boundary)
265
266    class TestPlex_3D_Box(BaseTestPlex_3D, unittest.TestCase):
267        CELLS = None
268        COORDS = None
269        def setUp(self):
270            self.plex = PETSc.DMPlex().createBoxMesh([1,1,1], simplex=True)
271
272    class TestPlex_3D_Boundary(BaseTestPlex_3D, unittest.TestCase):
273        CELLS = None
274        COORDS = None
275        def setUp(self):
276            boundary = PETSc.DMPlex().create(self.COMM)
277            boundary.createCubeBoundary([0., 0., 0.], [1., 1., 1.], [1, 1, 1])
278            boundary.setDimension(self.DIM-1)
279            self.plex = PETSc.DMPlex().generate(boundary)
280
281# --------------------------------------------------------------------
282
283PETSC_DIR = petsc4py.get_config()['PETSC_DIR']
284
285def check_dtype(method):
286    def wrapper(self, *args, **kwargs):
287        if PETSc.ScalarType is PETSc.ComplexType:
288            return
289        else:
290            return method(self, *args, **kwargs)
291    return wrapper
292
293def check_package(method):
294    def wrapper(self, *args, **kwargs):
295        if not PETSc.Sys.hasExternalPackage("hdf5"):
296            return
297        elif self.PARTITIONERTYPE != "simple" and \
298           not PETSc.Sys.hasExternalPackage(self.PARTITIONERTYPE):
299            return
300        else:
301            return method(self, *args, **kwargs)
302    return wrapper
303
304def check_nsize(method):
305    def wrapper(self, *args, **kwargs):
306        if PETSc.COMM_WORLD.size != self.NSIZE:
307            return
308        else:
309            return method(self, *args, **kwargs)
310    return wrapper
311
312class BaseTestPlexHDF5(object):
313    NSIZE = 4
314    NTIMES = 3
315
316    def setUp(self):
317        self.txtvwr = PETSc.Viewer()
318
319    def tearDown(self):
320        if not PETSc.COMM_WORLD.rank:
321            if os.path.exists(self.outfile()):
322                os.remove(self.outfile())
323            if os.path.exists(self.tmp_output_file()):
324                os.remove(self.tmp_output_file())
325        self.txtvwr = None
326
327    def _name(self):
328        return "%s_outformat-%s_%s" % (self.SUFFIX,
329                                       self.OUTFORMAT,
330                                       self.PARTITIONERTYPE)
331
332    def infile(self):
333        return os.path.join(PETSC_DIR, "share/petsc/datafiles/",
334                            "meshes/blockcylinder-50.h5")
335
336    def outfile(self):
337        return os.path.join("./temp_test_dmplex_%s.h5" % self._name())
338
339    def informat(self):
340        return PETSc.Viewer.Format.HDF5_XDMF
341
342    def outformat(self):
343        d = {"hdf5_petsc": PETSc.Viewer.Format.HDF5_PETSC,
344             "hdf5_xdmf": PETSc.Viewer.Format.HDF5_XDMF}
345        return d[self.OUTFORMAT]
346
347    def partitionerType(self):
348        d = {"simple": PETSc.Partitioner.Type.SIMPLE,
349             "ptscotch": PETSc.Partitioner.Type.PTSCOTCH,
350             "parmetis": PETSc.Partitioner.Type.PARMETIS}
351        return d[self.PARTITIONERTYPE]
352
353    def ref_output_file(self):
354        return os.path.join(PETSC_DIR, "src/dm/impls/plex/tutorials/",
355                            "output/ex5_%s.out" % self._name())
356
357    def tmp_output_file(self):
358        return os.path.join("./temp_test_dmplex_%s.out" % self._name())
359
360    def outputText(self, msg, comm):
361        if not comm.rank:
362            with open(self.tmp_output_file(), 'a') as f:
363                f.write(msg)
364
365    def outputPlex(self, plex):
366        self.txtvwr.createASCII(self.tmp_output_file(),
367                                mode='a', comm=plex.comm)
368        plex.view(viewer=self.txtvwr)
369        self.txtvwr.destroy()
370
371    @check_dtype
372    @check_package
373    @check_nsize
374    def testViewLoadCycle(self):
375        grank = PETSc.COMM_WORLD.rank
376        for i in range(self.NTIMES):
377            if i == 0:
378                infname = self.infile()
379                informt = self.informat()
380            else:
381                infname = self.outfile()
382                informt = self.outformat()
383            if self.HETEROGENEOUS:
384                mycolor = (grank > self.NTIMES - i)
385            else:
386                mycolor = 0
387            try:
388                import mpi4py
389            except ImportError:
390                self.skipTest('mpi4py') # throws special exception to signal test skip
391            mpicomm = PETSc.COMM_WORLD.tompi4py()
392            comm = PETSc.Comm(comm=mpicomm.Split(color=mycolor, key=grank))
393            if mycolor == 0:
394                self.outputText("Begin cycle %d\n" % i, comm)
395                plex = PETSc.DMPlex()
396                vwr = PETSc.ViewerHDF5()
397                # Create plex
398                plex.create(comm=comm)
399                plex.setName("DMPlex Object")
400                # Load data from XDMF into dm in parallel
401                vwr.create(infname, mode='r', comm=comm)
402                vwr.pushFormat(format=informt)
403                plex.load(viewer=vwr)
404                plex.setOptionsPrefix("loaded_")
405                plex.setFromOptions()
406                vwr.popFormat()
407                vwr.destroy()
408                self.outputPlex(plex)
409                # Test DM is indeed distributed
410                flg = plex.isDistributed()
411                self.outputText("Loaded mesh distributed? %s\n" %
412                                str(flg).upper(), comm)
413                # Interpolate
414                plex.interpolate()
415                plex.setOptionsPrefix("interpolated_")
416                plex.setFromOptions()
417                self.outputPlex(plex)
418                # Redistribute
419                part = plex.getPartitioner()
420                part.setType(self.partitionerType())
421                _ = plex.distribute(overlap=0)
422                plex.setOptionsPrefix("redistributed_")
423                plex.setFromOptions()
424                self.outputPlex(plex)
425                # Save redistributed dm to XDMF in parallel
426                vwr.create(self.outfile(), mode='w', comm=comm)
427                vwr.pushFormat(format=self.outformat())
428                plex.setName("DMPlex Object")
429                plex.view(viewer=vwr)
430                vwr.popFormat()
431                vwr.destroy()
432                # Destroy plex
433                plex.destroy()
434                self.outputText("End   cycle %d\n--------\n" % i, comm)
435            PETSc.COMM_WORLD.Barrier()
436        # Check that the output is identical to that of plex/tutorial/ex5.c.
437        self.assertTrue(filecmp.cmp(self.tmp_output_file(),
438                                    self.ref_output_file(), shallow=False),
439                        'Contents of the files not the same.')
440        PETSc.COMM_WORLD.Barrier()
441
442class BaseTestPlexHDF5Homogeneous(BaseTestPlexHDF5):
443    """Test save on N / load on N."""
444    SUFFIX = 0
445    HETEROGENEOUS = False
446
447class BaseTestPlexHDF5Heterogeneous(BaseTestPlexHDF5):
448    """Test save on N / load on M."""
449    SUFFIX = 1
450    HETEROGENEOUS = True
451
452class TestPlexHDF5PETSCSimpleHomogeneous(BaseTestPlexHDF5Homogeneous,
453                                         unittest.TestCase):
454    OUTFORMAT = "hdf5_petsc"
455    PARTITIONERTYPE = "simple"
456
457"""
458Skipping. PTScotch produces different distributions when run
459in a sequence in a single session.
460
461class TestPlexHDF5PETSCPTScotchHomogeneous(BaseTestPlexHDF5Homogeneous,
462                                           unittest.TestCase):
463    OUTFORMAT = "hdf5_petsc"
464    PARTITIONERTYPE = "ptscotch"
465"""
466
467class TestPlexHDF5PETSCParmetisHomogeneous(BaseTestPlexHDF5Homogeneous,
468                                           unittest.TestCase):
469    OUTFORMAT = "hdf5_petsc"
470    PARTITIONERTYPE = "parmetis"
471
472class TestPlexHDF5XDMFSimpleHomogeneous(BaseTestPlexHDF5Homogeneous,
473                                        unittest.TestCase):
474    OUTFORMAT = "hdf5_xdmf"
475    PARTITIONERTYPE = "simple"
476
477"""
478Skipping. PTScotch produces different distributions when run
479in a sequence in a single session.
480
481class TestPlexHDF5XDMFPTScotchHomogeneous(BaseTestPlexHDF5Homogeneous,
482                                          unittest.TestCase):
483    OUTFORMAT = "hdf5_xdmf"
484    PARTITIONERTYPE = "ptscotch"
485"""
486
487class TestPlexHDF5XDMFParmetisHomogeneous(BaseTestPlexHDF5Homogeneous,
488                                          unittest.TestCase):
489    OUTFORMAT = "hdf5_xdmf"
490    PARTITIONERTYPE = "parmetis"
491
492class TestPlexHDF5PETSCSimpleHeterogeneous(BaseTestPlexHDF5Heterogeneous,
493                                           unittest.TestCase):
494    OUTFORMAT = "hdf5_petsc"
495    PARTITIONERTYPE = "simple"
496
497"""
498Skipping. PTScotch produces different distributions when run
499in a sequence in a single session.
500
501class TestPlexHDF5PETSCPTScotchHeterogeneous(BaseTestPlexHDF5Heterogeneous,
502                                             unittest.TestCase):
503    OUTFORMAT = "hdf5_petsc"
504    PARTITIONERTYPE = "ptscotch"
505"""
506
507class TestPlexHDF5PETSCParmetisHeterogeneous(BaseTestPlexHDF5Heterogeneous,
508                                             unittest.TestCase):
509    OUTFORMAT = "hdf5_petsc"
510    PARTITIONERTYPE = "parmetis"
511
512class TestPlexHDF5XDMFSimpleHeterogeneous(BaseTestPlexHDF5Heterogeneous,
513                                          unittest.TestCase):
514    OUTFORMAT = "hdf5_xdmf"
515    PARTITIONERTYPE = "simple"
516
517class TestPlexHDF5XDMFPTScotchHeterogeneous(BaseTestPlexHDF5Heterogeneous,
518                                            unittest.TestCase):
519    OUTFORMAT = "hdf5_xdmf"
520    PARTITIONERTYPE = "ptscotch"
521
522class TestPlexHDF5XDMFParmetisHeterogeneous(BaseTestPlexHDF5Heterogeneous,
523                                            unittest.TestCase):
524    OUTFORMAT = "hdf5_xdmf"
525    PARTITIONERTYPE = "parmetis"
526
527# --------------------------------------------------------------------
528
529if __name__ == '__main__':
530    unittest.main()
531