1 static char help[] = "Tests save/load of plex/section/vec 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 rakn 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 PetscBool flg; 135 136 PetscFunctionBegin; 137 options->fname[0] = '\0'; 138 PetscOptionsBegin(comm, "", "DMPlex View/Load Test Options", "DMPLEX"); 139 PetscCall(PetscOptionsString("-fname", "The output mesh file", "ex12.c", options->fname, options->fname, sizeof(options->fname), &flg)); 140 PetscCall(PetscOptionsBool("-shell", "Use DMShell to wrap sections", "ex12.c", options->shell, &options->shell, NULL)); 141 PetscOptionsEnd(); 142 PetscFunctionReturn(0); 143 } 144 145 int main(int argc, char **argv) { 146 MPI_Comm comm; 147 PetscMPIInt size, rank, mycolor; 148 const char exampleDMPlexName[] = "exampleDMPlex"; 149 const char exampleSectionDMName[] = "exampleSectionDM"; 150 const char exampleVecName[] = "exampleVec"; 151 PetscScalar constraintValue = 1.5; 152 PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC; 153 AppCtx user; 154 155 PetscFunctionBeginUser; 156 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 157 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 158 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 159 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 160 PetscCheck(size >= 3, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes"); 161 162 /* Save */ 163 mycolor = (PetscMPIInt)(rank >= 2); 164 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm)); 165 if (mycolor == 0) { 166 DM dm; 167 PetscViewer viewer; 168 169 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer)); 170 /* Save exampleDMPlex */ 171 { 172 DM pdm; 173 const PetscInt faces[2] = {6, 1}; 174 PetscSF sf; 175 PetscInt overlap = 1; 176 177 PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, &dm)); 178 PetscCall(DMPlexDistribute(dm, overlap, &sf, &pdm)); 179 if (pdm) { 180 PetscCall(DMDestroy(&dm)); 181 dm = pdm; 182 } 183 PetscCall(PetscSFDestroy(&sf)); 184 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 185 PetscCall(PetscViewerPushFormat(viewer, format)); 186 PetscCall(DMPlexTopologyView(dm, viewer)); 187 PetscCall(DMPlexLabelsView(dm, viewer)); 188 PetscCall(PetscViewerPopFormat(viewer)); 189 } 190 /* Save coordinates */ 191 PetscCall(PetscViewerPushFormat(viewer, format)); 192 PetscCall(DMPlexCoordinatesView(dm, viewer)); 193 PetscCall(PetscViewerPopFormat(viewer)); 194 /* Save exampleVec */ 195 { 196 PetscInt pStart = -1, pEnd = -1; 197 DM sdm; 198 PetscSection section, gsection; 199 PetscBool includesConstraints = PETSC_FALSE; 200 Vec vec; 201 PetscScalar *array = NULL; 202 203 /* Create section */ 204 PetscCall(PetscSectionCreate(comm, §ion)); 205 PetscCall(PetscSectionSetNumFields(section, 2)); 206 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 207 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 208 switch (rank) { 209 case 0: 210 PetscCall(PetscSectionSetDof(section, 3, 2)); 211 PetscCall(PetscSectionSetDof(section, 12, 3)); 212 PetscCall(PetscSectionSetDof(section, 25, 2)); 213 PetscCall(PetscSectionSetConstraintDof(section, 12, 1)); 214 PetscCall(PetscSectionSetFieldDof(section, 3, 0, 2)); 215 PetscCall(PetscSectionSetFieldDof(section, 12, 0, 2)); 216 PetscCall(PetscSectionSetFieldDof(section, 12, 1, 1)); 217 PetscCall(PetscSectionSetFieldDof(section, 25, 1, 2)); 218 PetscCall(PetscSectionSetFieldConstraintDof(section, 12, 1, 1)); 219 break; 220 case 1: 221 PetscCall(PetscSectionSetDof(section, 0, 2)); 222 PetscCall(PetscSectionSetDof(section, 1, 1)); 223 PetscCall(PetscSectionSetDof(section, 8, 3)); 224 PetscCall(PetscSectionSetDof(section, 20, 2)); 225 PetscCall(PetscSectionSetConstraintDof(section, 8, 1)); 226 PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2)); 227 PetscCall(PetscSectionSetFieldDof(section, 8, 0, 2)); 228 PetscCall(PetscSectionSetFieldDof(section, 1, 1, 1)); 229 PetscCall(PetscSectionSetFieldDof(section, 8, 1, 1)); 230 PetscCall(PetscSectionSetFieldDof(section, 20, 1, 2)); 231 PetscCall(PetscSectionSetFieldConstraintDof(section, 8, 1, 1)); 232 break; 233 } 234 PetscCall(PetscSectionSetUp(section)); 235 { 236 const PetscInt indices[] = {2}; 237 const PetscInt indices1[] = {0}; 238 239 switch (rank) { 240 case 0: 241 PetscCall(PetscSectionSetConstraintIndices(section, 12, indices)); 242 PetscCall(PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1)); 243 break; 244 case 1: 245 PetscCall(PetscSectionSetConstraintIndices(section, 8, indices)); 246 PetscCall(PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1)); 247 break; 248 } 249 } 250 if (user.shell) { 251 PetscSF sf; 252 253 PetscCall(DMShellCreate(comm, &sdm)); 254 PetscCall(DMGetPointSF(dm, &sf)); 255 PetscCall(DMSetPointSF(sdm, sf)); 256 } else { 257 PetscCall(DMClone(dm, &sdm)); 258 } 259 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName)); 260 PetscCall(DMSetLocalSection(sdm, section)); 261 PetscCall(PetscSectionDestroy(§ion)); 262 PetscCall(DMPlexSectionView(dm, viewer, sdm)); 263 /* Create global vector */ 264 PetscCall(DMGetGlobalSection(sdm, &gsection)); 265 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 266 if (user.shell) { 267 PetscInt n = -1; 268 269 PetscCall(VecCreate(comm, &vec)); 270 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &n)); 271 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &n)); 272 PetscCall(VecSetSizes(vec, n, PETSC_DECIDE)); 273 PetscCall(VecSetUp(vec)); 274 } else { 275 PetscCall(DMGetGlobalVector(sdm, &vec)); 276 } 277 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 278 PetscCall(VecGetArrayWrite(vec, &array)); 279 if (includesConstraints) { 280 switch (rank) { 281 case 0: break; 282 case 1: 283 array[0] = 1.0; 284 array[1] = 1.1; 285 array[2] = 1.2; 286 array[3] = 1.3; 287 array[4] = 1.4; 288 array[5] = 1.5; 289 array[6] = 1.6; 290 array[7] = 1.7; 291 break; 292 } 293 } else { 294 switch (rank) { 295 case 0: break; 296 case 1: 297 array[0] = 1.0; 298 array[1] = 1.1; 299 array[2] = 1.2; 300 array[3] = 1.3; 301 array[4] = 1.4; 302 array[5] = 1.6; 303 array[6] = 1.7; 304 break; 305 } 306 } 307 PetscCall(VecRestoreArrayWrite(vec, &array)); 308 PetscCall(DMPlexGlobalVectorView(dm, viewer, sdm, vec)); 309 if (user.shell) { 310 PetscCall(VecDestroy(&vec)); 311 } else { 312 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 313 } 314 PetscCall(DMDestroy(&sdm)); 315 } 316 PetscCall(PetscViewerDestroy(&viewer)); 317 PetscCall(DMDestroy(&dm)); 318 } 319 PetscCallMPI(MPI_Comm_free(&comm)); 320 /* Load */ 321 mycolor = (PetscMPIInt)(rank >= 3); 322 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm)); 323 if (mycolor == 0) { 324 DM dm; 325 PetscSF sfXC; 326 PetscViewer viewer; 327 328 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer)); 329 /* Load exampleDMPlex */ 330 { 331 PetscSF sfXB, sfBC; 332 333 PetscCall(DMCreate(comm, &dm)); 334 PetscCall(DMSetType(dm, DMPLEX)); 335 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 336 /* sfXB: X -> B */ 337 /* X: set of globalPointNumbers, [0, N) */ 338 /* B: loaded naive in-memory plex */ 339 PetscCall(PetscViewerPushFormat(viewer, format)); 340 PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB)); 341 PetscCall(PetscViewerPopFormat(viewer)); 342 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 343 { 344 DM distributedDM; 345 PetscInt overlap = 1; 346 PetscPartitioner part; 347 348 PetscCall(DMPlexGetPartitioner(dm, &part)); 349 PetscCall(PetscPartitionerSetFromOptions(part)); 350 /* sfBC: B -> C */ 351 /* B: loaded naive in-memory plex */ 352 /* C: redistributed good in-memory */ 353 PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM)); 354 if (distributedDM) { 355 PetscCall(DMDestroy(&dm)); 356 dm = distributedDM; 357 } 358 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 359 } 360 /* sfXC: X -> C */ 361 PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC)); 362 PetscCall(PetscSFDestroy(&sfXB)); 363 PetscCall(PetscSFDestroy(&sfBC)); 364 } 365 /* Load labels */ 366 PetscCall(PetscViewerPushFormat(viewer, format)); 367 PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC)); 368 PetscCall(PetscViewerPopFormat(viewer)); 369 /* Load coordinates */ 370 PetscCall(PetscViewerPushFormat(viewer, format)); 371 PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC)); 372 PetscCall(PetscViewerPopFormat(viewer)); 373 PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)")); 374 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 375 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 376 /* Load exampleVec */ 377 { 378 DM sdm; 379 PetscSection section, gsection; 380 IS perm; 381 PetscBool includesConstraints = PETSC_FALSE; 382 Vec vec; 383 PetscSF lsf, gsf; 384 385 if (user.shell) { 386 PetscSF sf; 387 388 PetscCall(DMShellCreate(comm, &sdm)); 389 PetscCall(DMGetPointSF(dm, &sf)); 390 PetscCall(DMSetPointSF(sdm, sf)); 391 } else { 392 PetscCall(DMClone(dm, &sdm)); 393 } 394 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName)); 395 PetscCall(PetscSectionCreate(comm, §ion)); 396 { 397 PetscInt pStart = -1, pEnd = -1, p = -1; 398 PetscInt *pinds = NULL; 399 400 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 401 PetscCall(PetscMalloc1(pEnd - pStart, &pinds)); 402 for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p; 403 if (rank == 2) { 404 pinds[10] = 20; 405 pinds[20] = 10; 406 } 407 PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm)); 408 } 409 PetscCall(PetscSectionSetPermutation(section, perm)); 410 PetscCall(ISDestroy(&perm)); 411 PetscCall(DMSetLocalSection(sdm, section)); 412 PetscCall(PetscSectionDestroy(§ion)); 413 PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf)); 414 /* Load as local vector */ 415 PetscCall(DMGetLocalSection(sdm, §ion)); 416 PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section")); 417 PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm))); 418 PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints)); 419 if (user.shell) { 420 PetscInt m = -1; 421 422 PetscCall(VecCreate(comm, &vec)); 423 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m)); 424 else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m)); 425 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 426 PetscCall(VecSetUp(vec)); 427 } else { 428 PetscCall(DMGetLocalVector(sdm, &vec)); 429 } 430 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 431 PetscCall(VecSet(vec, constraintValue)); 432 PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec)); 433 PetscCall(PetscSFDestroy(&lsf)); 434 if (user.shell) { 435 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 436 PetscCall(VecDestroy(&vec)); 437 } else { 438 PetscCall(DMRestoreLocalVector(sdm, &vec)); 439 } 440 /* Load as global vector */ 441 PetscCall(DMGetGlobalSection(sdm, &gsection)); 442 PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section")); 443 PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm))); 444 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 445 if (user.shell) { 446 PetscInt m = -1; 447 448 PetscCall(VecCreate(comm, &vec)); 449 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m)); 450 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m)); 451 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 452 PetscCall(VecSetUp(vec)); 453 } else { 454 PetscCall(DMGetGlobalVector(sdm, &vec)); 455 } 456 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 457 PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec)); 458 PetscCall(PetscSFDestroy(&gsf)); 459 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 460 if (user.shell) { 461 PetscCall(VecDestroy(&vec)); 462 } else { 463 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 464 } 465 PetscCall(DMDestroy(&sdm)); 466 } 467 PetscCall(PetscViewerDestroy(&viewer)); 468 PetscCall(PetscSFDestroy(&sfXC)); 469 PetscCall(DMDestroy(&dm)); 470 } 471 PetscCallMPI(MPI_Comm_free(&comm)); 472 473 /* Finalize */ 474 PetscCall(PetscFinalize()); 475 return 0; 476 } 477 478 /*TEST 479 480 build: 481 requires: hdf5 482 testset: 483 suffix: 0 484 requires: !complex 485 nsize: 4 486 args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail 487 args: -dm_plex_view_hdf5_storage_version 2.0.0 488 test: 489 suffix: parmetis 490 requires: parmetis 491 args: -petscpartitioner_type parmetis 492 test: 493 suffix: ptscotch 494 requires: ptscotch 495 args: -petscpartitioner_type ptscotch 496 497 TEST*/ 498