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 PetscErrorCode ierr; 137 138 PetscFunctionBegin; 139 options->fname[0] = '\0'; 140 ierr = PetscOptionsBegin(comm, "", "DMPlex View/Load Test Options", "DMPLEX");PetscCall(ierr); 141 PetscCall(PetscOptionsString("-fname", "The output mesh file", "ex12.c", options->fname, options->fname, sizeof(options->fname), &flg)); 142 PetscCall(PetscOptionsBool("-shell", "Use DMShell to wrap sections", "ex12.c", options->shell, &options->shell, NULL)); 143 ierr = PetscOptionsEnd();PetscCall(ierr); 144 PetscFunctionReturn(0); 145 } 146 147 int main(int argc, char **argv) 148 { 149 MPI_Comm comm; 150 PetscMPIInt size, rank, mycolor; 151 const char exampleDMPlexName[] = "exampleDMPlex"; 152 const char exampleSectionDMName[] = "exampleSectionDM"; 153 const char exampleVecName[] = "exampleVec"; 154 PetscScalar constraintValue = 1.5; 155 PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC; 156 AppCtx user; 157 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 PetscCheckFalse(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, &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 } 259 else { 260 PetscCall(DMClone(dm, &sdm)); 261 } 262 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName)); 263 PetscCall(DMSetLocalSection(sdm, section)); 264 PetscCall(PetscSectionDestroy(§ion)); 265 PetscCall(DMPlexSectionView(dm, viewer, sdm)); 266 /* Create global vector */ 267 PetscCall(DMGetGlobalSection(sdm, &gsection)); 268 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 269 if (user.shell) { 270 PetscInt n = -1; 271 272 PetscCall(VecCreate(comm, &vec)); 273 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &n)); 274 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &n)); 275 PetscCall(VecSetSizes(vec, n, PETSC_DECIDE)); 276 PetscCall(VecSetUp(vec)); 277 } else { 278 PetscCall(DMGetGlobalVector(sdm, &vec)); 279 } 280 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 281 PetscCall(VecGetArrayWrite(vec, &array)); 282 if (includesConstraints) { 283 switch (rank) { 284 case 0: 285 break; 286 case 1: 287 array[0] = 1.0; 288 array[1] = 1.1; 289 array[2] = 1.2; 290 array[3] = 1.3; 291 array[4] = 1.4; 292 array[5] = 1.5; 293 array[6] = 1.6; 294 array[7] = 1.7; 295 break; 296 } 297 } else { 298 switch (rank) { 299 case 0: 300 break; 301 case 1: 302 array[0] = 1.0; 303 array[1] = 1.1; 304 array[2] = 1.2; 305 array[3] = 1.3; 306 array[4] = 1.4; 307 array[5] = 1.6; 308 array[6] = 1.7; 309 break; 310 } 311 } 312 PetscCall(VecRestoreArrayWrite(vec, &array)); 313 PetscCall(DMPlexGlobalVectorView(dm, viewer, sdm, vec)); 314 if (user.shell) { 315 PetscCall(VecDestroy(&vec)); 316 } else { 317 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 318 } 319 PetscCall(DMDestroy(&sdm)); 320 } 321 PetscCall(PetscViewerDestroy(&viewer)); 322 PetscCall(DMDestroy(&dm)); 323 } 324 PetscCallMPI(MPI_Comm_free(&comm)); 325 /* Load */ 326 mycolor = (PetscMPIInt)(rank >= 3); 327 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm)); 328 if (mycolor == 0) { 329 DM dm; 330 PetscSF sfXC; 331 PetscViewer viewer; 332 333 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer)); 334 /* Load exampleDMPlex */ 335 { 336 PetscSF sfXB, sfBC; 337 338 PetscCall(DMCreate(comm, &dm)); 339 PetscCall(DMSetType(dm, DMPLEX)); 340 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 341 /* sfXB: X -> B */ 342 /* X: set of globalPointNumbers, [0, N) */ 343 /* B: loaded naive in-memory plex */ 344 PetscCall(PetscViewerPushFormat(viewer, format)); 345 PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB)); 346 PetscCall(PetscViewerPopFormat(viewer)); 347 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 348 { 349 DM distributedDM; 350 PetscInt overlap = 1; 351 PetscPartitioner part; 352 353 PetscCall(DMPlexGetPartitioner(dm, &part)); 354 PetscCall(PetscPartitionerSetFromOptions(part)); 355 /* sfBC: B -> C */ 356 /* B: loaded naive in-memory plex */ 357 /* C: redistributed good in-memory */ 358 PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM)); 359 if (distributedDM) { 360 PetscCall(DMDestroy(&dm)); 361 dm = distributedDM; 362 } 363 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 364 } 365 /* sfXC: X -> C */ 366 PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC)); 367 PetscCall(PetscSFDestroy(&sfXB)); 368 PetscCall(PetscSFDestroy(&sfBC)); 369 } 370 /* Load labels */ 371 PetscCall(PetscViewerPushFormat(viewer, format)); 372 PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC)); 373 PetscCall(PetscViewerPopFormat(viewer)); 374 /* Load coordinates */ 375 PetscCall(PetscViewerPushFormat(viewer, format)); 376 PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC)); 377 PetscCall(PetscViewerPopFormat(viewer)); 378 PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)")); 379 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 380 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName)); 381 /* Load exampleVec */ 382 { 383 DM sdm; 384 PetscSection section, gsection; 385 IS perm; 386 PetscBool includesConstraints = PETSC_FALSE; 387 Vec vec; 388 PetscSF lsf, gsf; 389 390 if (user.shell) { 391 PetscSF sf; 392 393 PetscCall(DMShellCreate(comm, &sdm)); 394 PetscCall(DMGetPointSF(dm, &sf)); 395 PetscCall(DMSetPointSF(sdm, sf)); 396 } else { 397 PetscCall(DMClone(dm, &sdm)); 398 } 399 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName)); 400 PetscCall(PetscSectionCreate(comm, §ion)); 401 { 402 PetscInt pStart = -1, pEnd = -1, p = -1; 403 PetscInt *pinds = NULL; 404 405 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 406 PetscCall(PetscMalloc1(pEnd - pStart, &pinds)); 407 for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p; 408 if (rank == 2) {pinds[10] = 20; pinds[20] = 10;} 409 PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm)); 410 } 411 PetscCall(PetscSectionSetPermutation(section, perm)); 412 PetscCall(ISDestroy(&perm)); 413 PetscCall(DMSetLocalSection(sdm, section)); 414 PetscCall(PetscSectionDestroy(§ion)); 415 PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf)); 416 /* Load as local vector */ 417 PetscCall(DMGetLocalSection(sdm, §ion)); 418 PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section")); 419 PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm))); 420 PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints)); 421 if (user.shell) { 422 PetscInt m = -1; 423 424 PetscCall(VecCreate(comm, &vec)); 425 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m)); 426 else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m)); 427 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 428 PetscCall(VecSetUp(vec)); 429 } else { 430 PetscCall(DMGetLocalVector(sdm, &vec)); 431 } 432 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 433 PetscCall(VecSet(vec, constraintValue)); 434 PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec)); 435 PetscCall(PetscSFDestroy(&lsf)); 436 if (user.shell) { 437 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 438 PetscCall(VecDestroy(&vec)); 439 } else { 440 PetscCall(DMRestoreLocalVector(sdm, &vec)); 441 } 442 /* Load as global vector */ 443 PetscCall(DMGetGlobalSection(sdm, &gsection)); 444 PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section")); 445 PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm))); 446 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 447 if (user.shell) { 448 PetscInt m = -1; 449 450 PetscCall(VecCreate(comm, &vec)); 451 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m)); 452 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m)); 453 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 454 PetscCall(VecSetUp(vec)); 455 } else { 456 PetscCall(DMGetGlobalVector(sdm, &vec)); 457 } 458 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 459 PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec)); 460 PetscCall(PetscSFDestroy(&gsf)); 461 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 462 if (user.shell) { 463 PetscCall(VecDestroy(&vec)); 464 } else { 465 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 466 } 467 PetscCall(DMDestroy(&sdm)); 468 } 469 PetscCall(PetscViewerDestroy(&viewer)); 470 PetscCall(PetscSFDestroy(&sfXC)); 471 PetscCall(DMDestroy(&dm)); 472 } 473 PetscCallMPI(MPI_Comm_free(&comm)); 474 475 /* Finalize */ 476 PetscCall(PetscFinalize()); 477 return 0; 478 } 479 480 /*TEST 481 482 build: 483 requires: hdf5 484 testset: 485 suffix: 0 486 requires: !complex 487 nsize: 4 488 args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail 489 args: -dm_plex_view_hdf5_storage_version 2.0.0 490 test: 491 suffix: parmetis 492 requires: parmetis 493 args: -petscpartitioner_type parmetis 494 test: 495 suffix: ptscotch 496 requires: ptscotch 497 args: -petscpartitioner_type ptscotch 498 499 TEST*/ 500