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