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) 191 assert np.allclose(metric.array, metric1.array) 192 self.plex.metricNormalize(metric, metric1, det, 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 pass 236 237class TestPlex_3D(BaseTestPlex_3D, unittest.TestCase): 238 pass 239 240class TestPlex_2D_P3(BaseTestPlex_2D, unittest.TestCase): 241 DOFS = [1, 2, 1] 242 243class TestPlex_3D_P3(BaseTestPlex_3D, unittest.TestCase): 244 DOFS = [1, 2, 1, 0] 245 246class TestPlex_3D_P4(BaseTestPlex_3D, unittest.TestCase): 247 DOFS = [1, 3, 3, 1] 248 249class TestPlex_2D_BoxTensor(BaseTestPlex_2D, unittest.TestCase): 250 CELLS = None 251 COORDS = None 252 def setUp(self): 253 self.plex = PETSc.DMPlex().createBoxMesh([3,3], simplex=False) 254 255class TestPlex_3D_BoxTensor(BaseTestPlex_3D, unittest.TestCase): 256 CELLS = None 257 COORDS = None 258 def setUp(self): 259 self.plex = PETSc.DMPlex().createBoxMesh([3,3,3], simplex=False) 260 261try: 262 raise PETSc.Error 263 PETSc.DMPlex().createBoxMesh([2,2], simplex=True, comm=PETSc.COMM_SELF).destroy() 264except PETSc.Error: 265 pass 266else: 267 class TestPlex_2D_Box(BaseTestPlex_2D, unittest.TestCase): 268 CELLS = None 269 COORDS = None 270 def setUp(self): 271 self.plex = PETSc.DMPlex().createBoxMesh([1,1], simplex=True) 272 273 class TestPlex_2D_Boundary(BaseTestPlex_2D, unittest.TestCase): 274 CELLS = None 275 COORDS = None 276 def setUp(self): 277 boundary = PETSc.DMPlex().create(self.COMM) 278 boundary.createSquareBoundary([0., 0.], [1., 1.], [2, 2]) 279 boundary.setDimension(self.DIM-1) 280 self.plex = PETSc.DMPlex().generate(boundary) 281 282 class TestPlex_3D_Box(BaseTestPlex_3D, unittest.TestCase): 283 CELLS = None 284 COORDS = None 285 def setUp(self): 286 self.plex = PETSc.DMPlex().createBoxMesh([1,1,1], simplex=True) 287 288 class TestPlex_3D_Boundary(BaseTestPlex_3D, unittest.TestCase): 289 CELLS = None 290 COORDS = None 291 def setUp(self): 292 boundary = PETSc.DMPlex().create(self.COMM) 293 boundary.createCubeBoundary([0., 0., 0.], [1., 1., 1.], [1, 1, 1]) 294 boundary.setDimension(self.DIM-1) 295 self.plex = PETSc.DMPlex().generate(boundary) 296 297# -------------------------------------------------------------------- 298 299PETSC_DIR = petsc4py.get_config()['PETSC_DIR'] 300 301def check_dtype(method): 302 def wrapper(self, *args, **kwargs): 303 if PETSc.ScalarType is PETSc.ComplexType: 304 return 305 else: 306 return method(self, *args, **kwargs) 307 return wrapper 308 309def check_package(method): 310 def wrapper(self, *args, **kwargs): 311 if not PETSc.Sys.hasExternalPackage("hdf5"): 312 return 313 elif self.PARTITIONERTYPE != "simple" and \ 314 not PETSc.Sys.hasExternalPackage(self.PARTITIONERTYPE): 315 return 316 else: 317 return method(self, *args, **kwargs) 318 return wrapper 319 320def check_nsize(method): 321 def wrapper(self, *args, **kwargs): 322 if PETSc.COMM_WORLD.size != self.NSIZE: 323 return 324 else: 325 return method(self, *args, **kwargs) 326 return wrapper 327 328class BaseTestPlexHDF5(object): 329 NSIZE = 4 330 NTIMES = 3 331 332 def setUp(self): 333 self.txtvwr = PETSc.Viewer() 334 335 def tearDown(self): 336 if not PETSc.COMM_WORLD.rank: 337 if os.path.exists(self.outfile()): 338 os.remove(self.outfile()) 339 if os.path.exists(self.tmp_output_file()): 340 os.remove(self.tmp_output_file()) 341 self.txtvwr = None 342 343 def _name(self): 344 return "%s_outformat-%s_%s" % (self.SUFFIX, 345 self.OUTFORMAT, 346 self.PARTITIONERTYPE) 347 348 def infile(self): 349 return os.path.join(PETSC_DIR, "share/petsc/datafiles/", 350 "meshes/blockcylinder-50.h5") 351 352 def outfile(self): 353 return os.path.join("./temp_test_dmplex_%s.h5" % self._name()) 354 355 def informat(self): 356 return PETSc.Viewer.Format.HDF5_XDMF 357 358 def outformat(self): 359 d = {"hdf5_petsc": PETSc.Viewer.Format.HDF5_PETSC, 360 "hdf5_xdmf": PETSc.Viewer.Format.HDF5_XDMF} 361 return d[self.OUTFORMAT] 362 363 def partitionerType(self): 364 d = {"simple": PETSc.Partitioner.Type.SIMPLE, 365 "ptscotch": PETSc.Partitioner.Type.PTSCOTCH, 366 "parmetis": PETSc.Partitioner.Type.PARMETIS} 367 return d[self.PARTITIONERTYPE] 368 369 def ref_output_file(self): 370 return os.path.join(PETSC_DIR, "src/dm/impls/plex/tutorials/", 371 "output/ex5_%s.out" % self._name()) 372 373 def tmp_output_file(self): 374 return os.path.join("./temp_test_dmplex_%s.out" % self._name()) 375 376 def outputText(self, msg, comm): 377 if not comm.rank: 378 with open(self.tmp_output_file(), 'a') as f: 379 f.write(msg) 380 381 def outputPlex(self, plex): 382 self.txtvwr.createASCII(self.tmp_output_file(), 383 mode='a', comm=plex.comm) 384 plex.view(viewer=self.txtvwr) 385 self.txtvwr.destroy() 386 387 @check_dtype 388 @check_package 389 @check_nsize 390 def testViewLoadCycle(self): 391 grank = PETSc.COMM_WORLD.rank 392 for i in range(self.NTIMES): 393 if i == 0: 394 infname = self.infile() 395 informt = self.informat() 396 else: 397 infname = self.outfile() 398 informt = self.outformat() 399 if self.HETEROGENEOUS: 400 mycolor = (grank > self.NTIMES - i) 401 else: 402 mycolor = 0 403 try: 404 import mpi4py 405 except ImportError: 406 self.skipTest('mpi4py') # throws special exception to signal test skip 407 mpicomm = PETSc.COMM_WORLD.tompi4py() 408 comm = PETSc.Comm(comm=mpicomm.Split(color=mycolor, key=grank)) 409 if mycolor == 0: 410 self.outputText("Begin cycle %d\n" % i, comm) 411 plex = PETSc.DMPlex() 412 vwr = PETSc.ViewerHDF5() 413 # Create plex 414 plex.create(comm=comm) 415 plex.setName("DMPlex Object") 416 # Load data from XDMF into dm in parallel 417 vwr.create(infname, mode='r', comm=comm) 418 vwr.pushFormat(format=informt) 419 plex.load(viewer=vwr) 420 plex.setOptionsPrefix("loaded_") 421 plex.distributeSetDefault(False) 422 plex.setFromOptions() 423 vwr.popFormat() 424 vwr.destroy() 425 self.outputPlex(plex) 426 # Test DM is indeed distributed 427 flg = plex.isDistributed() 428 self.outputText("Loaded mesh distributed? %s\n" % 429 str(flg).upper(), comm) 430 # Interpolate 431 plex.interpolate() 432 plex.setOptionsPrefix("interpolated_") 433 plex.setFromOptions() 434 self.outputPlex(plex) 435 # Redistribute 436 part = plex.getPartitioner() 437 part.setType(self.partitionerType()) 438 _ = plex.distribute(overlap=0) 439 plex.setName("DMPlex Object") 440 plex.setOptionsPrefix("redistributed_") 441 plex.setFromOptions() 442 self.outputPlex(plex) 443 # Save redistributed dm to XDMF in parallel 444 vwr.create(self.outfile(), mode='w', comm=comm) 445 vwr.pushFormat(format=self.outformat()) 446 plex.setName("DMPlex Object") 447 plex.view(viewer=vwr) 448 vwr.popFormat() 449 vwr.destroy() 450 # Destroy plex 451 plex.destroy() 452 self.outputText("End cycle %d\n--------\n" % i, comm) 453 PETSc.COMM_WORLD.Barrier() 454 # Check that the output is identical to that of plex/tutorial/ex5.c. 455 self.assertTrue(filecmp.cmp(self.tmp_output_file(), 456 self.ref_output_file(), shallow=False), 457 'Contents of the files not the same.') 458 PETSc.COMM_WORLD.Barrier() 459 460class BaseTestPlexHDF5Homogeneous(BaseTestPlexHDF5): 461 """Test save on N / load on N.""" 462 SUFFIX = 0 463 HETEROGENEOUS = False 464 465class BaseTestPlexHDF5Heterogeneous(BaseTestPlexHDF5): 466 """Test save on N / load on M.""" 467 SUFFIX = 1 468 HETEROGENEOUS = True 469 470class TestPlexHDF5PETSCSimpleHomogeneous(BaseTestPlexHDF5Homogeneous, 471 unittest.TestCase): 472 OUTFORMAT = "hdf5_petsc" 473 PARTITIONERTYPE = "simple" 474 475""" 476Skipping. PTScotch produces different distributions when run 477in a sequence in a single session. 478 479class TestPlexHDF5PETSCPTScotchHomogeneous(BaseTestPlexHDF5Homogeneous, 480 unittest.TestCase): 481 OUTFORMAT = "hdf5_petsc" 482 PARTITIONERTYPE = "ptscotch" 483""" 484 485class TestPlexHDF5PETSCParmetisHomogeneous(BaseTestPlexHDF5Homogeneous, 486 unittest.TestCase): 487 OUTFORMAT = "hdf5_petsc" 488 PARTITIONERTYPE = "parmetis" 489 490class TestPlexHDF5XDMFSimpleHomogeneous(BaseTestPlexHDF5Homogeneous, 491 unittest.TestCase): 492 OUTFORMAT = "hdf5_xdmf" 493 PARTITIONERTYPE = "simple" 494 495""" 496Skipping. PTScotch produces different distributions when run 497in a sequence in a single session. 498 499class TestPlexHDF5XDMFPTScotchHomogeneous(BaseTestPlexHDF5Homogeneous, 500 unittest.TestCase): 501 OUTFORMAT = "hdf5_xdmf" 502 PARTITIONERTYPE = "ptscotch" 503""" 504 505class TestPlexHDF5XDMFParmetisHomogeneous(BaseTestPlexHDF5Homogeneous, 506 unittest.TestCase): 507 OUTFORMAT = "hdf5_xdmf" 508 PARTITIONERTYPE = "parmetis" 509 510class TestPlexHDF5PETSCSimpleHeterogeneous(BaseTestPlexHDF5Heterogeneous, 511 unittest.TestCase): 512 OUTFORMAT = "hdf5_petsc" 513 PARTITIONERTYPE = "simple" 514 515""" 516Skipping. PTScotch produces different distributions when run 517in a sequence in a single session. 518 519class TestPlexHDF5PETSCPTScotchHeterogeneous(BaseTestPlexHDF5Heterogeneous, 520 unittest.TestCase): 521 OUTFORMAT = "hdf5_petsc" 522 PARTITIONERTYPE = "ptscotch" 523""" 524 525class TestPlexHDF5PETSCParmetisHeterogeneous(BaseTestPlexHDF5Heterogeneous, 526 unittest.TestCase): 527 OUTFORMAT = "hdf5_petsc" 528 PARTITIONERTYPE = "parmetis" 529 530class TestPlexHDF5XDMFSimpleHeterogeneous(BaseTestPlexHDF5Heterogeneous, 531 unittest.TestCase): 532 OUTFORMAT = "hdf5_xdmf" 533 PARTITIONERTYPE = "simple" 534 535class TestPlexHDF5XDMFPTScotchHeterogeneous(BaseTestPlexHDF5Heterogeneous, 536 unittest.TestCase): 537 OUTFORMAT = "hdf5_xdmf" 538 PARTITIONERTYPE = "ptscotch" 539 540class TestPlexHDF5XDMFParmetisHeterogeneous(BaseTestPlexHDF5Heterogeneous, 541 unittest.TestCase): 542 OUTFORMAT = "hdf5_xdmf" 543 PARTITIONERTYPE = "parmetis" 544 545# -------------------------------------------------------------------- 546 547if __name__ == '__main__': 548 unittest.main() 549