1 static char help[] = "Tests save/load of plex/section/vec on different numbers of processes in HDF5.\n\n"; 2 3 #include <petscdmshell.h> 4 #include <petscdmplex.h> 5 #include <petscsection.h> 6 #include <petscsf.h> 7 #include <petsclayouthdf5.h> 8 9 /* A six-element mesh 10 11 ===================== 12 Save on 2 processes 13 ===================== 14 15 exampleDMPlex: Local numbering: 16 17 7---17--8---18--9--19--(12)(24)(13) 18 | | | | | 19 rank 0: 20 0 21 1 22 2 (25) (3)(26) 20 | | | | | 21 4---14--5---15--6--16--(10)(23)(11) 22 23 (13)(25)--8--17---9--18--10--19--11 24 | | | | | 25 rank 1: (26) (3) 20 0 21 1 22 2 23 26 | | | | | 27 (12)(24)--4--14---5--15---6--16---7 28 29 exampleDMPlex: globalPointNumbering: 30 31 9--23--10--24--11--25--16--32--17--33--18--34--19 32 | | | | | | | 33 26 0 27 1 28 2 35 3 36 4 37 5 38 34 | | | | | | | 35 6--20---7--21---8--22--12--29--13--30--14--31--15 36 37 exampleSectionDM: 38 - includesConstraints = TRUE for local section (default) 39 - includesConstraints = FALSE for global section (default) 40 41 exampleSectionDM: Dofs (Field 0): 42 43 0---0---0---0---0---0---2---0---0---0---0---0---0 44 | | | | | | | 45 0 0 0 0 0 0 0 2 0 0 0 0 0 46 | | | | | | | 47 0---0---0---0---0---0---0---0---0---0---0---0---0 48 49 exampleSectionDM: Dofs (Field 1): constrained 50 / 51 0---0---0---0---0---0---1---0---0---0---0---0---0 52 | | | | | | | 53 0 0 0 0 0 0 2 0 0 1 0 0 0 54 | | | | | | | 55 0---0---0---0---0---0---0---0---0---0---0---0---0 56 57 exampleSectionDM: Offsets (total) in global section: 58 59 0---0---0---0---0---0---3---5---5---5---5---5---5 60 | | | | | | | 61 0 0 0 0 0 0 5 0 7 2 7 3 7 62 | | | | | | | 63 0---0---0---0---0---0---3---5---3---5---3---5---3 64 65 exampleVec: Values (Field 0): (1.3, 1.4) 66 / 67 +-------+-------+-------*-------+-------+-------+ 68 | | | | | | | 69 | | | | * (1.0, 1.1)| | 70 | | | | | | | 71 +-------+-------+-------+-------+-------+-------+ 72 73 exampleVec: Values (Field 1): (1.5,) constrained 74 / 75 +-------+-------+-------*-------+-------+-------+ 76 | | | | | | | 77 | | (1.6, 1.7) * | * (1.2,) | 78 | | | | | | | 79 +-------+-------+-------+-------+-------+-------+ 80 81 exampleVec: as global vector 82 83 rank 0: [] 84 rank 1: [1.0, 1.1, 1.2, 1.3, 1.4, 1.6, 1.7] 85 86 ===================== 87 Load on 3 Processes 88 ===================== 89 90 exampleDMPlex: Loaded/Distributed: 91 92 5--13---6--14--(8)(18)(10) 93 | | | | 94 rank 0: 15 0 16 1 (19)(2)(20) 95 | | | | 96 3--11---4--12--(7)(17)-(9) 97 98 (9)(21)--5--15---7--18-(12)(24)(13) 99 | | | | | 100 rank 1: (22) (2) 16 0 19 1 (25) (3)(26) 101 | | | | | 102 (8)(20)--4--14---6--17-(10)(23)(11) 103 104 +-> (10)(19)--6--13---7--14---8 105 permute | | | | | 106 rank 2: +-> (20) (2) 15 0 16 1 17 107 | | | | 108 (9)(18)--3--11---4--12---5 109 110 exampleSectionDM: 111 - includesConstraints = TRUE for local section (default) 112 - includesConstraints = FALSE for global section (default) 113 114 exampleVec: as local vector: 115 116 rank 0: [1.3, 1.4, 1.5, 1.6, 1.7] 117 rank 1: [1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7] 118 rank 2: [1.2, 1.0, 1.1, 1.6, 1.7, 1.3, 1.4, 1.5] 119 120 exampleVec: as global vector: 121 122 rank 0: [] 123 rank 1: [1.0, 1.1, 1.3, 1.4, 1.6, 1.7] 124 rank 2: [1.2] 125 126 */ 127 128 typedef struct { 129 char fname[PETSC_MAX_PATH_LEN]; /* Output mesh filename */ 130 PetscBool shell; /* Use DMShell to wrap sections */ 131 } AppCtx; 132 133 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 134 { 135 PetscBool flg; 136 137 PetscFunctionBegin; 138 options->fname[0] = '\0'; 139 PetscOptionsBegin(comm, "", "DMPlex View/Load Test Options", "DMPLEX"); 140 PetscCall(PetscOptionsString("-fname", "The output mesh file", "ex12.c", options->fname, options->fname, sizeof(options->fname), &flg)); 141 PetscCall(PetscOptionsBool("-shell", "Use DMShell to wrap sections", "ex12.c", options->shell, &options->shell, NULL)); 142 PetscOptionsEnd(); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 int main(int argc, char **argv) 147 { 148 MPI_Comm comm; 149 PetscMPIInt size, rank, mycolor; 150 const char exampleDMPlexName[] = "exampleDMPlex"; 151 const char exampleSectionDMName[] = "exampleSectionDM"; 152 const char exampleVecName[] = "exampleVec"; 153 PetscScalar constraintValue = 1.5; 154 PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC; 155 AppCtx user; 156 157 PetscFunctionBeginUser; 158 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 159 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 160 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 161 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 162 PetscCheck(size >= 3, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes"); 163 164 /* Save */ 165 mycolor = (PetscMPIInt)(rank >= 2); 166 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm)); 167 if (mycolor == 0) { 168 DM dm; 169 PetscViewer viewer; 170 171 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer)); 172 /* Save exampleDMPlex */ 173 { 174 DM pdm; 175 const PetscInt faces[2] = {6, 1}; 176 PetscSF sf; 177 PetscInt overlap = 1; 178 179 PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &dm)); 180 PetscCall(DMPlexDistribute(dm, overlap, &sf, &pdm)); 181 if (pdm) { 182 PetscCall(DMDestroy(&dm)); 183 dm = pdm; 184 } 185 PetscCall(PetscSFDestroy(&sf)); 186 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 187 PetscCall(PetscViewerPushFormat(viewer, format)); 188 PetscCall(DMPlexTopologyView(dm, viewer)); 189 PetscCall(DMPlexLabelsView(dm, viewer)); 190 PetscCall(PetscViewerPopFormat(viewer)); 191 } 192 /* Save coordinates */ 193 PetscCall(PetscViewerPushFormat(viewer, format)); 194 PetscCall(DMPlexCoordinatesView(dm, viewer)); 195 PetscCall(PetscViewerPopFormat(viewer)); 196 /* Save exampleVec */ 197 { 198 PetscInt pStart = -1, pEnd = -1; 199 DM sdm; 200 PetscSection section, gsection; 201 PetscBool includesConstraints = PETSC_FALSE; 202 Vec vec; 203 PetscScalar *array = NULL; 204 205 /* Create section */ 206 PetscCall(PetscSectionCreate(comm, §ion)); 207 PetscCall(PetscSectionSetNumFields(section, 2)); 208 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 209 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 210 switch (rank) { 211 case 0: 212 PetscCall(PetscSectionSetDof(section, 3, 2)); 213 PetscCall(PetscSectionSetDof(section, 12, 3)); 214 PetscCall(PetscSectionSetDof(section, 25, 2)); 215 PetscCall(PetscSectionSetConstraintDof(section, 12, 1)); 216 PetscCall(PetscSectionSetFieldDof(section, 3, 0, 2)); 217 PetscCall(PetscSectionSetFieldDof(section, 12, 0, 2)); 218 PetscCall(PetscSectionSetFieldDof(section, 12, 1, 1)); 219 PetscCall(PetscSectionSetFieldDof(section, 25, 1, 2)); 220 PetscCall(PetscSectionSetFieldConstraintDof(section, 12, 1, 1)); 221 break; 222 case 1: 223 PetscCall(PetscSectionSetDof(section, 0, 2)); 224 PetscCall(PetscSectionSetDof(section, 1, 1)); 225 PetscCall(PetscSectionSetDof(section, 8, 3)); 226 PetscCall(PetscSectionSetDof(section, 20, 2)); 227 PetscCall(PetscSectionSetConstraintDof(section, 8, 1)); 228 PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2)); 229 PetscCall(PetscSectionSetFieldDof(section, 8, 0, 2)); 230 PetscCall(PetscSectionSetFieldDof(section, 1, 1, 1)); 231 PetscCall(PetscSectionSetFieldDof(section, 8, 1, 1)); 232 PetscCall(PetscSectionSetFieldDof(section, 20, 1, 2)); 233 PetscCall(PetscSectionSetFieldConstraintDof(section, 8, 1, 1)); 234 break; 235 } 236 PetscCall(PetscSectionSetUp(section)); 237 { 238 const PetscInt indices[] = {2}; 239 const PetscInt indices1[] = {0}; 240 241 switch (rank) { 242 case 0: 243 PetscCall(PetscSectionSetConstraintIndices(section, 12, indices)); 244 PetscCall(PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1)); 245 break; 246 case 1: 247 PetscCall(PetscSectionSetConstraintIndices(section, 8, indices)); 248 PetscCall(PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1)); 249 break; 250 } 251 } 252 if (user.shell) { 253 PetscSF sf; 254 255 PetscCall(DMShellCreate(comm, &sdm)); 256 PetscCall(DMGetPointSF(dm, &sf)); 257 PetscCall(DMSetPointSF(sdm, sf)); 258 } else { 259 PetscCall(DMClone(dm, &sdm)); 260 } 261 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName)); 262 PetscCall(DMSetLocalSection(sdm, section)); 263 PetscCall(PetscSectionDestroy(§ion)); 264 PetscCall(DMPlexSectionView(dm, viewer, sdm)); 265 /* Create global vector */ 266 PetscCall(DMGetGlobalSection(sdm, &gsection)); 267 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 268 if (user.shell) { 269 PetscInt n = -1; 270 271 PetscCall(VecCreate(comm, &vec)); 272 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &n)); 273 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &n)); 274 PetscCall(VecSetSizes(vec, n, PETSC_DECIDE)); 275 PetscCall(VecSetUp(vec)); 276 } else { 277 PetscCall(DMGetGlobalVector(sdm, &vec)); 278 } 279 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 280 PetscCall(VecGetArrayWrite(vec, &array)); 281 if (includesConstraints) { 282 switch (rank) { 283 case 0: 284 break; 285 case 1: 286 array[0] = 1.0; 287 array[1] = 1.1; 288 array[2] = 1.2; 289 array[3] = 1.3; 290 array[4] = 1.4; 291 array[5] = 1.5; 292 array[6] = 1.6; 293 array[7] = 1.7; 294 break; 295 } 296 } else { 297 switch (rank) { 298 case 0: 299 break; 300 case 1: 301 array[0] = 1.0; 302 array[1] = 1.1; 303 array[2] = 1.2; 304 array[3] = 1.3; 305 array[4] = 1.4; 306 array[5] = 1.6; 307 array[6] = 1.7; 308 break; 309 } 310 } 311 PetscCall(VecRestoreArrayWrite(vec, &array)); 312 PetscCall(DMPlexGlobalVectorView(dm, viewer, sdm, vec)); 313 if (user.shell) { 314 PetscCall(VecDestroy(&vec)); 315 } else { 316 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 317 } 318 PetscCall(DMDestroy(&sdm)); 319 } 320 PetscCall(PetscViewerDestroy(&viewer)); 321 PetscCall(DMDestroy(&dm)); 322 } 323 PetscCallMPI(MPI_Comm_free(&comm)); 324 /* Load */ 325 mycolor = (PetscMPIInt)(rank >= 3); 326 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm)); 327 if (mycolor == 0) { 328 DM dm; 329 PetscSF sfXC; 330 PetscViewer viewer; 331 332 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer)); 333 /* Load exampleDMPlex */ 334 { 335 PetscSF sfXB, sfBC; 336 337 PetscCall(DMCreate(comm, &dm)); 338 PetscCall(DMSetType(dm, DMPLEX)); 339 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 340 /* sfXB: X -> B */ 341 /* X: set of globalPointNumbers, [0, N) */ 342 /* B: loaded naive in-memory plex */ 343 PetscCall(PetscViewerPushFormat(viewer, format)); 344 PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB)); 345 PetscCall(PetscViewerPopFormat(viewer)); 346 { 347 DM distributedDM; 348 PetscInt overlap = 1; 349 PetscPartitioner part; 350 351 PetscCall(DMPlexGetPartitioner(dm, &part)); 352 PetscCall(PetscPartitionerSetFromOptions(part)); 353 /* sfBC: B -> C */ 354 /* B: loaded naive in-memory plex */ 355 /* C: redistributed good in-memory */ 356 PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM)); 357 if (distributedDM) { 358 PetscCall(DMDestroy(&dm)); 359 dm = distributedDM; 360 } 361 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 362 } 363 /* sfXC: X -> C */ 364 PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC)); 365 PetscCall(PetscSFDestroy(&sfXB)); 366 PetscCall(PetscSFDestroy(&sfBC)); 367 } 368 /* Load labels */ 369 PetscCall(PetscViewerPushFormat(viewer, format)); 370 PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC)); 371 PetscCall(PetscViewerPopFormat(viewer)); 372 /* Load coordinates */ 373 PetscCall(PetscViewerPushFormat(viewer, format)); 374 PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC)); 375 PetscCall(PetscViewerPopFormat(viewer)); 376 PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)")); 377 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 378 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 379 /* Load exampleVec */ 380 { 381 DM sdm; 382 PetscSection section, gsection; 383 IS perm; 384 PetscBool includesConstraints = PETSC_FALSE; 385 Vec vec; 386 PetscSF lsf, gsf; 387 388 if (user.shell) { 389 PetscSF sf; 390 391 PetscCall(DMShellCreate(comm, &sdm)); 392 PetscCall(DMGetPointSF(dm, &sf)); 393 PetscCall(DMSetPointSF(sdm, sf)); 394 } else { 395 PetscCall(DMClone(dm, &sdm)); 396 } 397 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName)); 398 PetscCall(PetscSectionCreate(comm, §ion)); 399 { 400 PetscInt pStart = -1, pEnd = -1, p = -1; 401 PetscInt *pinds = NULL; 402 403 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 404 PetscCall(PetscMalloc1(pEnd - pStart, &pinds)); 405 for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p; 406 if (rank == 2) { 407 pinds[10] = 20; 408 pinds[20] = 10; 409 } 410 PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm)); 411 } 412 PetscCall(PetscSectionSetPermutation(section, perm)); 413 PetscCall(ISDestroy(&perm)); 414 PetscCall(DMSetLocalSection(sdm, section)); 415 PetscCall(PetscSectionDestroy(§ion)); 416 PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf)); 417 /* Load as local vector */ 418 PetscCall(DMGetLocalSection(sdm, §ion)); 419 PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section")); 420 PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm))); 421 PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints)); 422 if (user.shell) { 423 PetscInt m = -1; 424 425 PetscCall(VecCreate(comm, &vec)); 426 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m)); 427 else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m)); 428 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 429 PetscCall(VecSetUp(vec)); 430 } else { 431 PetscCall(DMGetLocalVector(sdm, &vec)); 432 } 433 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 434 PetscCall(VecSet(vec, constraintValue)); 435 PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec)); 436 PetscCall(PetscSFDestroy(&lsf)); 437 if (user.shell) { 438 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 439 PetscCall(VecDestroy(&vec)); 440 } else { 441 PetscCall(DMRestoreLocalVector(sdm, &vec)); 442 } 443 /* Load as global vector */ 444 PetscCall(DMGetGlobalSection(sdm, &gsection)); 445 PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section")); 446 PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm))); 447 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 448 if (user.shell) { 449 PetscInt m = -1; 450 451 PetscCall(VecCreate(comm, &vec)); 452 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m)); 453 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m)); 454 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 455 PetscCall(VecSetUp(vec)); 456 } else { 457 PetscCall(DMGetGlobalVector(sdm, &vec)); 458 } 459 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 460 PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec)); 461 PetscCall(PetscSFDestroy(&gsf)); 462 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 463 if (user.shell) { 464 PetscCall(VecDestroy(&vec)); 465 } else { 466 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 467 } 468 PetscCall(DMDestroy(&sdm)); 469 } 470 PetscCall(PetscViewerDestroy(&viewer)); 471 PetscCall(PetscSFDestroy(&sfXC)); 472 PetscCall(DMDestroy(&dm)); 473 } 474 PetscCallMPI(MPI_Comm_free(&comm)); 475 476 /* Finalize */ 477 PetscCall(PetscFinalize()); 478 return 0; 479 } 480 481 /*TEST 482 483 build: 484 requires: hdf5 485 testset: 486 suffix: 0 487 requires: !complex 488 nsize: 4 489 args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail 490 args: -dm_plex_view_hdf5_storage_version 2.0.0 491 test: 492 suffix: parmetis 493 requires: parmetis 494 args: -petscpartitioner_type parmetis 495 test: 496 suffix: ptscotch 497 requires: ptscotch 498 args: -petscpartitioner_type ptscotch 499 500 TEST*/ 501