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 { 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(0); 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 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 158 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 159 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 160 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 161 PetscCheck(size >= 3,PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes"); 162 163 /* Save */ 164 mycolor = (PetscMPIInt)(rank >= 2); 165 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm)); 166 if (mycolor == 0) { 167 DM dm; 168 PetscViewer viewer; 169 170 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer)); 171 /* Save exampleDMPlex */ 172 { 173 DM pdm; 174 const PetscInt faces[2] = {6, 1}; 175 PetscSF sf; 176 PetscInt overlap = 1; 177 178 PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, &dm)); 179 PetscCall(DMPlexDistribute(dm, overlap, &sf, &pdm)); 180 if (pdm) { 181 PetscCall(DMDestroy(&dm)); 182 dm = pdm; 183 } 184 PetscCall(PetscSFDestroy(&sf)); 185 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 186 PetscCall(PetscViewerPushFormat(viewer, format)); 187 PetscCall(DMPlexTopologyView(dm, viewer)); 188 PetscCall(DMPlexLabelsView(dm, viewer)); 189 PetscCall(PetscViewerPopFormat(viewer)); 190 } 191 /* Save coordinates */ 192 PetscCall(PetscViewerPushFormat(viewer, format)); 193 PetscCall(DMPlexCoordinatesView(dm, viewer)); 194 PetscCall(PetscViewerPopFormat(viewer)); 195 /* Save exampleVec */ 196 { 197 PetscInt pStart = -1, pEnd = -1; 198 DM sdm; 199 PetscSection section, gsection; 200 PetscBool includesConstraints = PETSC_FALSE; 201 Vec vec; 202 PetscScalar *array = NULL; 203 204 /* Create section */ 205 PetscCall(PetscSectionCreate(comm, §ion)); 206 PetscCall(PetscSectionSetNumFields(section, 2)); 207 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 208 PetscCall(PetscSectionSetChart(section, pStart, pEnd)); 209 switch (rank) { 210 case 0: 211 PetscCall(PetscSectionSetDof(section, 3, 2)); 212 PetscCall(PetscSectionSetDof(section, 12, 3)); 213 PetscCall(PetscSectionSetDof(section, 25, 2)); 214 PetscCall(PetscSectionSetConstraintDof(section, 12, 1)); 215 PetscCall(PetscSectionSetFieldDof(section, 3, 0, 2)); 216 PetscCall(PetscSectionSetFieldDof(section, 12, 0, 2)); 217 PetscCall(PetscSectionSetFieldDof(section, 12, 1, 1)); 218 PetscCall(PetscSectionSetFieldDof(section, 25, 1, 2)); 219 PetscCall(PetscSectionSetFieldConstraintDof(section, 12, 1, 1)); 220 break; 221 case 1: 222 PetscCall(PetscSectionSetDof(section, 0, 2)); 223 PetscCall(PetscSectionSetDof(section, 1, 1)); 224 PetscCall(PetscSectionSetDof(section, 8, 3)); 225 PetscCall(PetscSectionSetDof(section, 20, 2)); 226 PetscCall(PetscSectionSetConstraintDof(section, 8, 1)); 227 PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2)); 228 PetscCall(PetscSectionSetFieldDof(section, 8, 0, 2)); 229 PetscCall(PetscSectionSetFieldDof(section, 1, 1, 1)); 230 PetscCall(PetscSectionSetFieldDof(section, 8, 1, 1)); 231 PetscCall(PetscSectionSetFieldDof(section, 20, 1, 2)); 232 PetscCall(PetscSectionSetFieldConstraintDof(section, 8, 1, 1)); 233 break; 234 } 235 PetscCall(PetscSectionSetUp(section)); 236 { 237 const PetscInt indices[] = {2}; 238 const PetscInt indices1[] = {0}; 239 240 switch (rank) { 241 case 0: 242 PetscCall(PetscSectionSetConstraintIndices(section, 12, indices)); 243 PetscCall(PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1)); 244 break; 245 case 1: 246 PetscCall(PetscSectionSetConstraintIndices(section, 8, indices)); 247 PetscCall(PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1)); 248 break; 249 } 250 } 251 if (user.shell) { 252 PetscSF sf; 253 254 PetscCall(DMShellCreate(comm, &sdm)); 255 PetscCall(DMGetPointSF(dm, &sf)); 256 PetscCall(DMSetPointSF(sdm, sf)); 257 } 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 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 347 { 348 DM distributedDM; 349 PetscInt overlap = 1; 350 PetscPartitioner part; 351 352 PetscCall(DMPlexGetPartitioner(dm, &part)); 353 PetscCall(PetscPartitionerSetFromOptions(part)); 354 /* sfBC: B -> C */ 355 /* B: loaded naive in-memory plex */ 356 /* C: redistributed good in-memory */ 357 PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM)); 358 if (distributedDM) { 359 PetscCall(DMDestroy(&dm)); 360 dm = distributedDM; 361 } 362 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 363 } 364 /* sfXC: X -> C */ 365 PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC)); 366 PetscCall(PetscSFDestroy(&sfXB)); 367 PetscCall(PetscSFDestroy(&sfBC)); 368 } 369 /* Load labels */ 370 PetscCall(PetscViewerPushFormat(viewer, format)); 371 PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC)); 372 PetscCall(PetscViewerPopFormat(viewer)); 373 /* Load coordinates */ 374 PetscCall(PetscViewerPushFormat(viewer, format)); 375 PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC)); 376 PetscCall(PetscViewerPopFormat(viewer)); 377 PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)")); 378 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 379 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 380 /* Load exampleVec */ 381 { 382 DM sdm; 383 PetscSection section, gsection; 384 IS perm; 385 PetscBool includesConstraints = PETSC_FALSE; 386 Vec vec; 387 PetscSF lsf, gsf; 388 389 if (user.shell) { 390 PetscSF sf; 391 392 PetscCall(DMShellCreate(comm, &sdm)); 393 PetscCall(DMGetPointSF(dm, &sf)); 394 PetscCall(DMSetPointSF(sdm, sf)); 395 } else { 396 PetscCall(DMClone(dm, &sdm)); 397 } 398 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName)); 399 PetscCall(PetscSectionCreate(comm, §ion)); 400 { 401 PetscInt pStart = -1, pEnd = -1, p = -1; 402 PetscInt *pinds = NULL; 403 404 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 405 PetscCall(PetscMalloc1(pEnd - pStart, &pinds)); 406 for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p; 407 if (rank == 2) {pinds[10] = 20; pinds[20] = 10;} 408 PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm)); 409 } 410 PetscCall(PetscSectionSetPermutation(section, perm)); 411 PetscCall(ISDestroy(&perm)); 412 PetscCall(DMSetLocalSection(sdm, section)); 413 PetscCall(PetscSectionDestroy(§ion)); 414 PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf)); 415 /* Load as local vector */ 416 PetscCall(DMGetLocalSection(sdm, §ion)); 417 PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section")); 418 PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm))); 419 PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints)); 420 if (user.shell) { 421 PetscInt m = -1; 422 423 PetscCall(VecCreate(comm, &vec)); 424 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m)); 425 else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m)); 426 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 427 PetscCall(VecSetUp(vec)); 428 } else { 429 PetscCall(DMGetLocalVector(sdm, &vec)); 430 } 431 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 432 PetscCall(VecSet(vec, constraintValue)); 433 PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec)); 434 PetscCall(PetscSFDestroy(&lsf)); 435 if (user.shell) { 436 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 437 PetscCall(VecDestroy(&vec)); 438 } else { 439 PetscCall(DMRestoreLocalVector(sdm, &vec)); 440 } 441 /* Load as global vector */ 442 PetscCall(DMGetGlobalSection(sdm, &gsection)); 443 PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section")); 444 PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm))); 445 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 446 if (user.shell) { 447 PetscInt m = -1; 448 449 PetscCall(VecCreate(comm, &vec)); 450 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m)); 451 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m)); 452 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 453 PetscCall(VecSetUp(vec)); 454 } else { 455 PetscCall(DMGetGlobalVector(sdm, &vec)); 456 } 457 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 458 PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec)); 459 PetscCall(PetscSFDestroy(&gsf)); 460 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 461 if (user.shell) { 462 PetscCall(VecDestroy(&vec)); 463 } else { 464 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 465 } 466 PetscCall(DMDestroy(&sdm)); 467 } 468 PetscCall(PetscViewerDestroy(&viewer)); 469 PetscCall(PetscSFDestroy(&sfXC)); 470 PetscCall(DMDestroy(&dm)); 471 } 472 PetscCallMPI(MPI_Comm_free(&comm)); 473 474 /* Finalize */ 475 PetscCall(PetscFinalize()); 476 return 0; 477 } 478 479 /*TEST 480 481 build: 482 requires: hdf5 483 testset: 484 suffix: 0 485 requires: !complex 486 nsize: 4 487 args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail 488 args: -dm_plex_view_hdf5_storage_version 2.0.0 489 test: 490 suffix: parmetis 491 requires: parmetis 492 args: -petscpartitioner_type parmetis 493 test: 494 suffix: ptscotch 495 requires: ptscotch 496 args: -petscpartitioner_type ptscotch 497 498 TEST*/ 499