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(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, &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 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) { 408 pinds[10] = 20; 409 pinds[20] = 10; 410 } 411 PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm)); 412 } 413 PetscCall(PetscSectionSetPermutation(section, perm)); 414 PetscCall(ISDestroy(&perm)); 415 PetscCall(DMSetLocalSection(sdm, section)); 416 PetscCall(PetscSectionDestroy(§ion)); 417 PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf)); 418 /* Load as local vector */ 419 PetscCall(DMGetLocalSection(sdm, §ion)); 420 PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section")); 421 PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm))); 422 PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints)); 423 if (user.shell) { 424 PetscInt m = -1; 425 426 PetscCall(VecCreate(comm, &vec)); 427 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m)); 428 else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m)); 429 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 430 PetscCall(VecSetUp(vec)); 431 } else { 432 PetscCall(DMGetLocalVector(sdm, &vec)); 433 } 434 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 435 PetscCall(VecSet(vec, constraintValue)); 436 PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec)); 437 PetscCall(PetscSFDestroy(&lsf)); 438 if (user.shell) { 439 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 440 PetscCall(VecDestroy(&vec)); 441 } else { 442 PetscCall(DMRestoreLocalVector(sdm, &vec)); 443 } 444 /* Load as global vector */ 445 PetscCall(DMGetGlobalSection(sdm, &gsection)); 446 PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section")); 447 PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm))); 448 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints)); 449 if (user.shell) { 450 PetscInt m = -1; 451 452 PetscCall(VecCreate(comm, &vec)); 453 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m)); 454 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m)); 455 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE)); 456 PetscCall(VecSetUp(vec)); 457 } else { 458 PetscCall(DMGetGlobalVector(sdm, &vec)); 459 } 460 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName)); 461 PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec)); 462 PetscCall(PetscSFDestroy(&gsf)); 463 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm))); 464 if (user.shell) { 465 PetscCall(VecDestroy(&vec)); 466 } else { 467 PetscCall(DMRestoreGlobalVector(sdm, &vec)); 468 } 469 PetscCall(DMDestroy(&sdm)); 470 } 471 PetscCall(PetscViewerDestroy(&viewer)); 472 PetscCall(PetscSFDestroy(&sfXC)); 473 PetscCall(DMDestroy(&dm)); 474 } 475 PetscCallMPI(MPI_Comm_free(&comm)); 476 477 /* Finalize */ 478 PetscCall(PetscFinalize()); 479 return 0; 480 } 481 482 /*TEST 483 484 build: 485 requires: hdf5 486 testset: 487 suffix: 0 488 requires: !complex 489 nsize: 4 490 args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail 491 args: -dm_plex_view_hdf5_storage_version 2.0.0 492 test: 493 suffix: parmetis 494 requires: parmetis 495 args: -petscpartitioner_type parmetis 496 test: 497 suffix: ptscotch 498 requires: ptscotch 499 args: -petscpartitioner_type ptscotch 500 501 TEST*/ 502