1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/isimpl.h> 3 #include <petsc/private/vecimpl.h> 4 #include <petsc/private/viewerhdf5impl.h> 5 #include <petsclayouthdf5.h> 6 7 /* Logging support */ 8 PetscLogEvent DMPLEX_DistributionView, DMPLEX_DistributionLoad; 9 10 static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *); 11 static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer, DMPlexStorageVersion); 12 static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer, const char[], DMPlexStorageVersion); 13 static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *); 14 15 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer); 16 17 static PetscErrorCode PetscViewerPrintVersion_Private(PetscViewer viewer, DMPlexStorageVersion version, char str[], size_t len) 18 { 19 PetscFunctionBegin; 20 PetscCall(PetscViewerCheckVersion_Private(viewer, version)); 21 PetscCall(PetscSNPrintf(str, len, "%d.%d.%d", version->major, version->minor, version->subminor)); 22 PetscFunctionReturn(PETSC_SUCCESS); 23 } 24 25 static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer viewer, const char str[], DMPlexStorageVersion *version) 26 { 27 PetscToken t; 28 char *ts; 29 PetscInt i; 30 PetscInt ti[3]; 31 DMPlexStorageVersion v; 32 33 PetscFunctionBegin; 34 PetscCall(PetscTokenCreate(str, '.', &t)); 35 for (i = 0; i < 3; i++) { 36 PetscCall(PetscTokenFind(t, &ts)); 37 PetscCheck(ts, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Malformed version string %s", str); 38 PetscCall(PetscOptionsStringToInt(ts, &ti[i])); 39 } 40 PetscCall(PetscTokenFind(t, &ts)); 41 PetscCheck(!ts, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Malformed version string %s", str); 42 PetscCall(PetscTokenDestroy(&t)); 43 PetscCall(PetscNew(&v)); 44 v->major = (int)ti[0]; 45 v->minor = (int)ti[1]; 46 v->subminor = (int)ti[2]; 47 PetscCall(PetscViewerCheckVersion_Private(viewer, v)); 48 *version = v; 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion v) 53 { 54 PetscFunctionBegin; 55 PetscCall(PetscObjectContainerCompose((PetscObject)viewer, key, v, PetscContainerUserDestroyDefault)); 56 PetscFunctionReturn(PETSC_SUCCESS); 57 } 58 59 static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion *v) 60 { 61 PetscContainer cont; 62 63 PetscFunctionBegin; 64 PetscCall(PetscObjectQuery((PetscObject)viewer, key, (PetscObject *)&cont)); 65 *v = NULL; 66 if (cont) PetscCall(PetscContainerGetPointer(cont, (void **)v)); 67 PetscFunctionReturn(PETSC_SUCCESS); 68 } 69 70 /* 71 Version log: 72 1.0.0 legacy version (default if no "dmplex_storage_version" attribute found in file) 73 1.1.0 legacy version, but output VIZ by default 74 2.0.0 introduce versioning and multiple topologies 75 2.1.0 introduce distributions 76 3.0.0 new checkpointing format in Firedrake paper 77 3.1.0 new format with IS compression 78 */ 79 static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer viewer, DMPlexStorageVersion version) 80 { 81 PetscBool valid = PETSC_FALSE; 82 83 PetscFunctionBegin; 84 switch (version->major) { 85 case 1: 86 switch (version->minor) { 87 case 0: 88 switch (version->subminor) { 89 case 0: 90 valid = PETSC_TRUE; 91 break; 92 } 93 break; 94 case 1: 95 switch (version->subminor) { 96 case 0: 97 valid = PETSC_TRUE; 98 break; 99 } 100 break; 101 } 102 break; 103 case 2: 104 switch (version->minor) { 105 case 0: 106 switch (version->subminor) { 107 case 0: 108 valid = PETSC_TRUE; 109 break; 110 } 111 break; 112 case 1: 113 switch (version->subminor) { 114 case 0: 115 valid = PETSC_TRUE; 116 break; 117 } 118 break; 119 } 120 break; 121 case 3: 122 switch (version->minor) { 123 case 0: 124 switch (version->subminor) { 125 case 0: 126 valid = PETSC_TRUE; 127 break; 128 } 129 break; 130 case 1: 131 switch (version->subminor) { 132 case 0: 133 valid = PETSC_TRUE; 134 break; 135 } 136 break; 137 } 138 break; 139 } 140 PetscCheck(valid, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "DMPlexStorageVersion %d.%d.%d not supported", version->major, version->minor, version->subminor); 141 PetscFunctionReturn(PETSC_SUCCESS); 142 } 143 144 static inline PetscBool DMPlexStorageVersionEQ(DMPlexStorageVersion version, int major, int minor, int subminor) 145 { 146 return (PetscBool)(version->major == major && version->minor == minor && version->subminor == subminor); 147 } 148 149 static inline PetscBool DMPlexStorageVersionGE(DMPlexStorageVersion version, int major, int minor, int subminor) 150 { 151 return (PetscBool)((version->major == major && version->minor == minor && version->subminor >= subminor) || (version->major == major && version->minor > minor) || (version->major > major)); 152 } 153 154 /*@C 155 PetscViewerHDF5SetDMPlexStorageVersionWriting - Set the storage version for writing 156 157 Logically collective 158 159 Input Parameters: 160 + viewer - The `PetscViewer` 161 - version - The storage format version 162 163 Level: advanced 164 165 Note: 166 The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0. 167 168 .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5GetDMPlexStorageVersionWriting()`, `PetscViewerHDF5GetDMPlexStorageVersionReading()`, `PetscViewerHDF5SetDMPlexStorageVersionReading()` 169 @*/ 170 PetscErrorCode PetscViewerHDF5SetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion version) 171 { 172 const char ATTR_NAME[] = "dmplex_storage_version"; 173 DMPlexStorageVersion viewerVersion; 174 PetscBool fileHasVersion; 175 char fileVersion[16], versionStr[16], viewerVersionStr[16]; 176 177 PetscFunctionBegin; 178 PetscValidHeaderSpecificType(viewer, PETSC_VIEWER_CLASSID, 1, PETSCVIEWERHDF5); 179 PetscAssertPointer(version, 2); 180 PetscCall(PetscViewerPrintVersion_Private(viewer, version, versionStr, 16)); 181 PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, &viewerVersion)); 182 if (viewerVersion) { 183 PetscBool flg; 184 185 PetscCall(PetscViewerPrintVersion_Private(viewer, viewerVersion, viewerVersionStr, 16)); 186 PetscCall(PetscStrcmp(versionStr, viewerVersionStr, &flg)); 187 PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but viewer already has version %s - cannot mix versions", versionStr, viewerVersionStr); 188 } 189 190 PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion)); 191 if (fileHasVersion) { 192 PetscBool flg; 193 char *tmp; 194 195 PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp)); 196 PetscCall(PetscStrncpy(fileVersion, tmp, sizeof(fileVersion))); 197 PetscCall(PetscFree(tmp)); 198 PetscCall(PetscStrcmp(fileVersion, versionStr, &flg)); 199 PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but file already has version %s - cannot mix versions", versionStr, fileVersion); 200 } else { 201 PetscCall(PetscViewerHDF5WriteAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, versionStr)); 202 } 203 PetscCall(PetscNew(&viewerVersion)); 204 viewerVersion->major = version->major; 205 viewerVersion->minor = version->minor; 206 viewerVersion->subminor = version->subminor; 207 PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, viewerVersion)); 208 PetscFunctionReturn(PETSC_SUCCESS); 209 } 210 211 /*@C 212 PetscViewerHDF5GetDMPlexStorageVersionWriting - Get the storage version for writing 213 214 Logically collective 215 216 Input Parameter: 217 . viewer - The `PetscViewer` 218 219 Output Parameter: 220 . version - The storage format version 221 222 Options Database Keys: 223 . -dm_plex_view_hdf5_storage_version <num> - Overrides the storage format version 224 225 Level: advanced 226 227 Note: 228 The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0. 229 230 .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5SetDMPlexStorageVersionWriting()`, `PetscViewerHDF5GetDMPlexStorageVersionReading()`, `PetscViewerHDF5SetDMPlexStorageVersionReading()` 231 @*/ 232 PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion *version) 233 { 234 const char ATTR_NAME[] = "dmplex_storage_version"; 235 PetscBool fileHasVersion; 236 char optVersion[16], fileVersion[16]; 237 238 PetscFunctionBegin; 239 PetscValidHeaderSpecificType(viewer, PETSC_VIEWER_CLASSID, 1, PETSCVIEWERHDF5); 240 PetscAssertPointer(version, 2); 241 PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, version)); 242 if (*version) PetscFunctionReturn(PETSC_SUCCESS); 243 244 PetscCall(PetscStrncpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE, sizeof(fileVersion))); 245 PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion)); 246 if (fileHasVersion) { 247 char *tmp; 248 249 PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp)); 250 PetscCall(PetscStrncpy(fileVersion, tmp, sizeof(fileVersion))); 251 PetscCall(PetscFree(tmp)); 252 } 253 PetscCall(PetscStrncpy(optVersion, fileVersion, sizeof(optVersion))); 254 PetscObjectOptionsBegin((PetscObject)viewer); 255 PetscCall(PetscOptionsString("-dm_plex_view_hdf5_storage_version", "DMPlex HDF5 viewer storage version", NULL, optVersion, optVersion, sizeof(optVersion), NULL)); 256 PetscOptionsEnd(); 257 if (!fileHasVersion) { 258 PetscCall(PetscViewerHDF5WriteAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, optVersion)); 259 } else { 260 PetscBool flg; 261 262 PetscCall(PetscStrcmp(fileVersion, optVersion, &flg)); 263 PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but file already has version %s - cannot mix versions", optVersion, fileVersion); 264 } 265 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT)); 266 PetscCall(PetscViewerParseVersion_Private(viewer, optVersion, version)); 267 PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, *version)); 268 PetscFunctionReturn(PETSC_SUCCESS); 269 } 270 271 /*@C 272 PetscViewerHDF5SetDMPlexStorageVersionReading - Set the storage version for reading 273 274 Logically collective 275 276 Input Parameters: 277 + viewer - The `PetscViewer` 278 - version - The storage format version 279 280 Level: advanced 281 282 Note: 283 The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0. 284 285 .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5GetDMPlexStorageVersionReading()`, `PetscViewerHDF5GetDMPlexStorageVersionWriting()`, `PetscViewerHDF5SetDMPlexStorageVersionWriting()` 286 @*/ 287 PetscErrorCode PetscViewerHDF5SetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion version) 288 { 289 const char ATTR_NAME[] = "dmplex_storage_version"; 290 DMPlexStorageVersion viewerVersion; 291 PetscBool fileHasVersion; 292 char versionStr[16], viewerVersionStr[16]; 293 294 PetscFunctionBegin; 295 PetscValidHeaderSpecificType(viewer, PETSC_VIEWER_CLASSID, 1, PETSCVIEWERHDF5); 296 PetscAssertPointer(version, 2); 297 PetscCall(PetscViewerPrintVersion_Private(viewer, version, versionStr, 16)); 298 PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, &viewerVersion)); 299 if (viewerVersion) { 300 PetscBool flg; 301 302 PetscCall(PetscViewerPrintVersion_Private(viewer, viewerVersion, viewerVersionStr, 16)); 303 PetscCall(PetscStrcmp(versionStr, viewerVersionStr, &flg)); 304 PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but viewer already has version %s - cannot mix versions", versionStr, viewerVersionStr); 305 } 306 307 PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion)); 308 if (fileHasVersion) { 309 char *fileVersion; 310 PetscBool flg; 311 312 PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &fileVersion)); 313 PetscCall(PetscStrcmp(fileVersion, versionStr, &flg)); 314 PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but file already has version %s - cannot mix versions", versionStr, fileVersion); 315 PetscCall(PetscFree(fileVersion)); 316 } 317 PetscCall(PetscNew(&viewerVersion)); 318 viewerVersion->major = version->major; 319 viewerVersion->minor = version->minor; 320 viewerVersion->subminor = version->subminor; 321 PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, viewerVersion)); 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 /*@C 326 PetscViewerHDF5GetDMPlexStorageVersionReading - Get the storage version for reading 327 328 Logically collective 329 330 Input Parameter: 331 . viewer - The `PetscViewer` 332 333 Output Parameter: 334 . version - The storage format version 335 336 Options Database Keys: 337 . -dm_plex_view_hdf5_storage_version <num> - Overrides the storage format version 338 339 Level: advanced 340 341 Note: 342 The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0. 343 344 .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5SetDMPlexStorageVersionReading()`, `PetscViewerHDF5GetDMPlexStorageVersionWriting()`, `PetscViewerHDF5SetDMPlexStorageVersionWriting()` 345 @*/ 346 PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion *version) 347 { 348 const char ATTR_NAME[] = "dmplex_storage_version"; 349 char *defaultVersion; 350 char *versionString; 351 352 PetscFunctionBegin; 353 PetscValidHeaderSpecificType(viewer, PETSC_VIEWER_CLASSID, 1, PETSCVIEWERHDF5); 354 PetscAssertPointer(version, 2); 355 PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, version)); 356 if (*version) PetscFunctionReturn(PETSC_SUCCESS); 357 358 //TODO string HDF5 attribute handling is terrible and should be redesigned 359 PetscCall(PetscStrallocpy(DMPLEX_STORAGE_VERSION_FIRST, &defaultVersion)); 360 PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString)); 361 PetscCall(PetscViewerParseVersion_Private(viewer, versionString, version)); 362 PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, *version)); 363 PetscCall(PetscFree(versionString)); 364 PetscCall(PetscFree(defaultVersion)); 365 PetscFunctionReturn(PETSC_SUCCESS); 366 } 367 368 static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[]) 369 { 370 PetscFunctionBegin; 371 if (((PetscObject)dm)->name) { 372 PetscCall(PetscObjectGetName((PetscObject)dm, name)); 373 } else { 374 *name = "plex"; 375 } 376 PetscFunctionReturn(PETSC_SUCCESS); 377 } 378 379 PetscErrorCode DMSequenceGetLength_HDF5_Internal(DM dm, const char seqname[], PetscInt *seqlen, PetscViewer viewer) 380 { 381 hid_t file, group, dset, dspace; 382 hsize_t rdim, *dims; 383 char *groupname; 384 PetscBool has; 385 386 PetscFunctionBegin; 387 PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &groupname)); 388 PetscCall(PetscViewerHDF5HasDataset(viewer, seqname, &has)); 389 PetscCheck(has, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Object (dataset) \"%s\" not stored in group %s", seqname, groupname); 390 391 PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file, &group)); 392 PetscCallHDF5Return(dset, H5Dopen2, (group, seqname, H5P_DEFAULT)); 393 PetscCallHDF5Return(dspace, H5Dget_space, (dset)); 394 PetscCallHDF5ReturnNoCheck(rdim, H5Sget_simple_extent_dims, (dspace, NULL, NULL)); 395 PetscCall(PetscMalloc1(rdim, &dims)); 396 PetscCallHDF5ReturnNoCheck(rdim, H5Sget_simple_extent_dims, (dspace, dims, NULL)); 397 *seqlen = (PetscInt)dims[0]; 398 PetscCall(PetscFree(dims)); 399 PetscCallHDF5(H5Dclose, (dset)); 400 PetscCallHDF5(H5Gclose, (group)); 401 PetscFunctionReturn(PETSC_SUCCESS); 402 } 403 404 static PetscErrorCode DMSequenceView_HDF5(DM dm, const char seqname[], PetscInt seqnum, PetscScalar value, PetscViewer viewer) 405 { 406 Vec stamp; 407 PetscMPIInt rank; 408 409 PetscFunctionBegin; 410 if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS); 411 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank)); 412 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp)); 413 PetscCall(VecSetBlockSize(stamp, 1)); 414 PetscCall(PetscObjectSetName((PetscObject)stamp, seqname)); 415 if (rank == 0) { 416 PetscReal timeScale; 417 PetscBool istime; 418 419 PetscCall(PetscStrncmp(seqname, "time", 5, &istime)); 420 if (istime) { 421 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale)); 422 value *= timeScale; 423 } 424 PetscCall(VecSetValue(stamp, 0, value, INSERT_VALUES)); 425 } 426 PetscCall(VecAssemblyBegin(stamp)); 427 PetscCall(VecAssemblyEnd(stamp)); 428 PetscCall(PetscViewerHDF5PushGroup(viewer, "/")); 429 PetscCall(PetscViewerHDF5PushTimestepping(viewer)); 430 PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */ 431 PetscCall(VecView(stamp, viewer)); 432 PetscCall(PetscViewerHDF5PopTimestepping(viewer)); 433 PetscCall(PetscViewerHDF5PopGroup(viewer)); 434 PetscCall(VecDestroy(&stamp)); 435 PetscFunctionReturn(PETSC_SUCCESS); 436 } 437 438 PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char seqname[], PetscInt seqnum, PetscScalar *value, PetscViewer viewer) 439 { 440 Vec stamp; 441 PetscMPIInt rank; 442 443 PetscFunctionBegin; 444 if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS); 445 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank)); 446 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp)); 447 PetscCall(VecSetBlockSize(stamp, 1)); 448 PetscCall(PetscObjectSetName((PetscObject)stamp, seqname)); 449 PetscCall(PetscViewerHDF5PushGroup(viewer, "/")); 450 PetscCall(PetscViewerHDF5PushTimestepping(viewer)); 451 PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */ 452 PetscCall(VecLoad(stamp, viewer)); 453 PetscCall(PetscViewerHDF5PopTimestepping(viewer)); 454 PetscCall(PetscViewerHDF5PopGroup(viewer)); 455 if (rank == 0) { 456 const PetscScalar *a; 457 PetscReal timeScale; 458 PetscBool istime; 459 460 PetscCall(VecGetArrayRead(stamp, &a)); 461 *value = a[0]; 462 PetscCall(VecRestoreArrayRead(stamp, &a)); 463 PetscCall(PetscStrncmp(seqname, "time", 5, &istime)); 464 if (istime) { 465 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale)); 466 *value /= timeScale; 467 } 468 } 469 PetscCall(VecDestroy(&stamp)); 470 PetscFunctionReturn(PETSC_SUCCESS); 471 } 472 473 static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel) 474 { 475 IS cutcells = NULL; 476 const PetscInt *cutc; 477 PetscInt cellHeight, vStart, vEnd, cStart, cEnd, c; 478 479 PetscFunctionBegin; 480 if (!cutLabel) PetscFunctionReturn(PETSC_SUCCESS); 481 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 482 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 483 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 484 /* Label vertices that should be duplicated */ 485 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel)); 486 PetscCall(DMLabelGetStratumIS(cutLabel, 2, &cutcells)); 487 if (cutcells) { 488 PetscInt n; 489 490 PetscCall(ISGetIndices(cutcells, &cutc)); 491 PetscCall(ISGetLocalSize(cutcells, &n)); 492 for (c = 0; c < n; ++c) { 493 if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) { 494 PetscInt *closure = NULL; 495 PetscInt closureSize, cl, value; 496 497 PetscCall(DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure)); 498 for (cl = 0; cl < closureSize * 2; cl += 2) { 499 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) { 500 PetscCall(DMLabelGetValue(cutLabel, closure[cl], &value)); 501 if (value == 1) PetscCall(DMLabelSetValue(*cutVertexLabel, closure[cl], 1)); 502 } 503 } 504 PetscCall(DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure)); 505 } 506 } 507 PetscCall(ISRestoreIndices(cutcells, &cutc)); 508 } 509 PetscCall(ISDestroy(&cutcells)); 510 PetscFunctionReturn(PETSC_SUCCESS); 511 } 512 513 PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer) 514 { 515 DM dm; 516 DM dmBC; 517 PetscSection section, sectionGlobal; 518 Vec gv; 519 const char *name; 520 PetscViewerFormat format; 521 PetscInt seqnum; 522 PetscReal seqval; 523 PetscBool isseq; 524 525 PetscFunctionBegin; 526 PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq)); 527 PetscCall(VecGetDM(v, &dm)); 528 PetscCall(DMGetLocalSection(dm, §ion)); 529 PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval)); 530 PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar)seqval, viewer)); 531 if (seqnum >= 0) { 532 PetscCall(PetscViewerHDF5PushTimestepping(viewer)); 533 PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); 534 } 535 PetscCall(PetscViewerGetFormat(viewer, &format)); 536 PetscCall(DMGetOutputDM(dm, &dmBC)); 537 PetscCall(DMGetGlobalSection(dmBC, §ionGlobal)); 538 PetscCall(DMGetGlobalVector(dmBC, &gv)); 539 PetscCall(PetscObjectGetName((PetscObject)v, &name)); 540 PetscCall(PetscObjectSetName((PetscObject)gv, name)); 541 PetscCall(DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv)); 542 PetscCall(DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv)); 543 PetscCall(PetscObjectTypeCompare((PetscObject)gv, VECSEQ, &isseq)); 544 if (format == PETSC_VIEWER_HDF5_VIZ) { 545 /* Output visualization representation */ 546 PetscInt numFields, f; 547 DMLabel cutLabel, cutVertexLabel = NULL; 548 549 PetscCall(PetscSectionGetNumFields(section, &numFields)); 550 PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel)); 551 for (f = 0; f < numFields; ++f) { 552 Vec subv; 553 IS is; 554 const char *fname, *fgroup, *componentName, *fname_def = "unnamed"; 555 char subname[PETSC_MAX_PATH_LEN]; 556 PetscInt Nc, c; 557 PetscInt pStart, pEnd; 558 PetscViewerVTKFieldType ft; 559 560 PetscCall(DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft)); 561 if (ft == PETSC_VTK_INVALID) continue; 562 fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields"; 563 PetscCall(PetscSectionGetFieldName(section, f, &fname)); 564 if (!fname) fname = fname_def; 565 566 PetscCall(PetscViewerHDF5PushGroup(viewer, fgroup)); 567 568 if (cutLabel) { 569 const PetscScalar *ga; 570 PetscScalar *suba; 571 PetscInt gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p; 572 573 PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel)); 574 PetscCall(PetscSectionGetFieldComponents(section, f, &Nc)); 575 for (p = pStart; p < pEnd; ++p) { 576 PetscInt gdof, fdof = 0, val; 577 578 PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof)); 579 if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof)); 580 subSize += fdof; 581 PetscCall(DMLabelGetValue(cutVertexLabel, p, &val)); 582 if (val == 1) extSize += fdof; 583 } 584 PetscCall(VecCreate(PetscObjectComm((PetscObject)gv), &subv)); 585 PetscCall(VecSetSizes(subv, subSize + extSize, PETSC_DETERMINE)); 586 PetscCall(VecSetBlockSize(subv, Nc)); 587 PetscCall(VecSetType(subv, VECSTANDARD)); 588 PetscCall(VecGetOwnershipRange(gv, &gstart, NULL)); 589 PetscCall(VecGetArrayRead(gv, &ga)); 590 PetscCall(VecGetArray(subv, &suba)); 591 for (p = pStart; p < pEnd; ++p) { 592 PetscInt gdof, goff, val; 593 594 PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof)); 595 if (gdof > 0) { 596 PetscInt fdof, fc, f2, poff = 0; 597 598 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 599 /* Can get rid of this loop by storing field information in the global section */ 600 for (f2 = 0; f2 < f; ++f2) { 601 PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof)); 602 poff += fdof; 603 } 604 PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof)); 605 for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff + poff + fc - gstart]; 606 PetscCall(DMLabelGetValue(cutVertexLabel, p, &val)); 607 if (val == 1) { 608 for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize + newOff] = ga[goff + poff + fc - gstart]; 609 } 610 } 611 } 612 PetscCall(VecRestoreArrayRead(gv, &ga)); 613 PetscCall(VecRestoreArray(subv, &suba)); 614 PetscCall(DMLabelDestroy(&cutVertexLabel)); 615 } else { 616 PetscCall(PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv)); 617 } 618 PetscCall(PetscStrncpy(subname, name, sizeof(subname))); 619 PetscCall(PetscStrlcat(subname, "_", sizeof(subname))); 620 PetscCall(PetscStrlcat(subname, fname, sizeof(subname))); 621 PetscCall(PetscObjectSetName((PetscObject)subv, subname)); 622 if (isseq) PetscCall(VecView_Seq(subv, viewer)); 623 else PetscCall(VecView_MPI(subv, viewer)); 624 if (ft == PETSC_VTK_POINT_VECTOR_FIELD || ft == PETSC_VTK_CELL_VECTOR_FIELD) { 625 PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "vector")); 626 } else { 627 PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "scalar")); 628 } 629 630 /* Output the component names in the field if available */ 631 PetscCall(PetscSectionGetFieldComponents(section, f, &Nc)); 632 for (c = 0; c < Nc; ++c) { 633 char componentNameLabel[PETSC_MAX_PATH_LEN]; 634 PetscCall(PetscSectionGetComponentName(section, f, c, &componentName)); 635 PetscCall(PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%" PetscInt_FMT, c)); 636 PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, componentNameLabel, PETSC_STRING, componentName)); 637 } 638 639 if (cutLabel) PetscCall(VecDestroy(&subv)); 640 else PetscCall(PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv)); 641 PetscCall(PetscViewerHDF5PopGroup(viewer)); 642 } 643 } else { 644 /* Output full vector */ 645 PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields")); 646 if (isseq) PetscCall(VecView_Seq(gv, viewer)); 647 else PetscCall(VecView_MPI(gv, viewer)); 648 PetscCall(PetscViewerHDF5PopGroup(viewer)); 649 } 650 if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer)); 651 PetscCall(DMRestoreGlobalVector(dmBC, &gv)); 652 PetscFunctionReturn(PETSC_SUCCESS); 653 } 654 655 PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer) 656 { 657 DMPlexStorageVersion version; 658 DM dm; 659 Vec locv; 660 PetscObject isZero; 661 const char *name; 662 PetscReal time; 663 664 PetscFunctionBegin; 665 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version)); 666 PetscCall(PetscInfo(v, "Writing Vec %s storage version %d.%d.%d\n", v->hdr.name, version->major, version->minor, version->subminor)); 667 668 PetscCall(VecGetDM(v, &dm)); 669 PetscCall(DMGetLocalVector(dm, &locv)); 670 PetscCall(PetscObjectGetName((PetscObject)v, &name)); 671 PetscCall(PetscObjectSetName((PetscObject)locv, name)); 672 PetscCall(PetscObjectQuery((PetscObject)v, "__Vec_bc_zero__", &isZero)); 673 PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", isZero)); 674 PetscCall(DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv)); 675 PetscCall(DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv)); 676 PetscCall(DMGetOutputSequenceNumber(dm, NULL, &time)); 677 PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL)); 678 PetscCall(VecView_Plex_Local_HDF5_Internal(locv, viewer)); 679 if (DMPlexStorageVersionEQ(version, 1, 1, 0)) { 680 PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields")); 681 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ)); 682 PetscCall(VecView_Plex_Local_HDF5_Internal(locv, viewer)); 683 PetscCall(PetscViewerPopFormat(viewer)); 684 PetscCall(PetscViewerHDF5PopGroup(viewer)); 685 } 686 PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", NULL)); 687 PetscCall(DMRestoreLocalVector(dm, &locv)); 688 PetscFunctionReturn(PETSC_SUCCESS); 689 } 690 691 PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer) 692 { 693 PetscBool isseq; 694 695 PetscFunctionBegin; 696 PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq)); 697 PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields")); 698 if (isseq) PetscCall(VecView_Seq(v, viewer)); 699 else PetscCall(VecView_MPI(v, viewer)); 700 PetscCall(PetscViewerHDF5PopGroup(viewer)); 701 PetscFunctionReturn(PETSC_SUCCESS); 702 } 703 704 PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer) 705 { 706 DM dm; 707 Vec locv; 708 const char *name; 709 PetscInt seqnum; 710 711 PetscFunctionBegin; 712 PetscCall(VecGetDM(v, &dm)); 713 PetscCall(DMGetLocalVector(dm, &locv)); 714 PetscCall(PetscObjectGetName((PetscObject)v, &name)); 715 PetscCall(PetscObjectSetName((PetscObject)locv, name)); 716 PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL)); 717 PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields")); 718 if (seqnum >= 0) { 719 PetscCall(PetscViewerHDF5PushTimestepping(viewer)); 720 PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); 721 } 722 PetscCall(VecLoad_Plex_Local(locv, viewer)); 723 if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer)); 724 PetscCall(PetscViewerHDF5PopGroup(viewer)); 725 PetscCall(DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v)); 726 PetscCall(DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v)); 727 PetscCall(DMRestoreLocalVector(dm, &locv)); 728 PetscFunctionReturn(PETSC_SUCCESS); 729 } 730 731 PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer) 732 { 733 DM dm; 734 PetscInt seqnum; 735 736 PetscFunctionBegin; 737 PetscCall(VecGetDM(v, &dm)); 738 PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL)); 739 PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields")); 740 if (seqnum >= 0) { 741 PetscCall(PetscViewerHDF5PushTimestepping(viewer)); 742 PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); 743 } 744 PetscCall(VecLoad_Default(v, viewer)); 745 if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer)); 746 PetscCall(PetscViewerHDF5PopGroup(viewer)); 747 PetscFunctionReturn(PETSC_SUCCESS); 748 } 749 750 static PetscErrorCode DMPlexDistributionView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer) 751 { 752 MPI_Comm comm; 753 PetscMPIInt size, rank; 754 PetscInt size_petsc_int; 755 const char *topologydm_name, *distribution_name; 756 const PetscInt *gpoint; 757 PetscInt pStart, pEnd, p; 758 PetscSF pointSF; 759 PetscInt nroots, nleaves; 760 const PetscInt *ilocal; 761 const PetscSFNode *iremote; 762 IS chartSizesIS, ownersIS, gpointsIS; 763 PetscInt *chartSize, *owners, *gpoints; 764 765 PetscFunctionBegin; 766 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 767 PetscCallMPI(MPI_Comm_size(comm, &size)); 768 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 769 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 770 PetscCall(DMPlexDistributionGetName(dm, &distribution_name)); 771 if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS); 772 PetscCall(PetscLogEventBegin(DMPLEX_DistributionView, viewer, 0, 0, 0)); 773 size_petsc_int = (PetscInt)size; 774 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "comm_size", PETSC_INT, (void *)&size_petsc_int)); 775 PetscCall(ISGetIndices(globalPointNumbers, &gpoint)); 776 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 777 PetscCall(PetscMalloc1(1, &chartSize)); 778 *chartSize = pEnd - pStart; 779 PetscCall(PetscMalloc1(*chartSize, &owners)); 780 PetscCall(PetscMalloc1(*chartSize, &gpoints)); 781 PetscCall(DMGetPointSF(dm, &pointSF)); 782 PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote)); 783 for (p = pStart; p < pEnd; ++p) { 784 PetscInt gp = gpoint[p - pStart]; 785 786 owners[p - pStart] = rank; 787 gpoints[p - pStart] = (gp < 0 ? -(gp + 1) : gp); 788 } 789 for (p = 0; p < nleaves; ++p) { 790 PetscInt ilocalp = (ilocal ? ilocal[p] : p); 791 792 owners[ilocalp] = iremote[p].rank; 793 } 794 PetscCall(ISCreateGeneral(comm, 1, chartSize, PETSC_OWN_POINTER, &chartSizesIS)); 795 PetscCall(ISCreateGeneral(comm, *chartSize, owners, PETSC_OWN_POINTER, &ownersIS)); 796 PetscCall(ISCreateGeneral(comm, *chartSize, gpoints, PETSC_OWN_POINTER, &gpointsIS)); 797 PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes")); 798 PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners")); 799 PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers")); 800 PetscCall(ISView(chartSizesIS, viewer)); 801 PetscCall(ISView(ownersIS, viewer)); 802 PetscCall(ISView(gpointsIS, viewer)); 803 PetscCall(ISDestroy(&chartSizesIS)); 804 PetscCall(ISDestroy(&ownersIS)); 805 PetscCall(ISDestroy(&gpointsIS)); 806 PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint)); 807 PetscCall(PetscLogEventEnd(DMPLEX_DistributionView, viewer, 0, 0, 0)); 808 PetscFunctionReturn(PETSC_SUCCESS); 809 } 810 811 static PetscErrorCode DMPlexTopologyView_HDF5_Inner_Private(DM dm, IS globalPointNumbers, PetscViewer viewer, PetscInt pStart, PetscInt pEnd, const char pointsName[], const char coneSizesName[], const char conesName[], const char orientationsName[]) 812 { 813 IS coneSizesIS, conesIS, orientationsIS; 814 PetscInt *coneSizes, *cones, *orientations; 815 const PetscInt *gpoint; 816 PetscInt nPoints = 0, conesSize = 0; 817 PetscInt p, c, s; 818 MPI_Comm comm; 819 820 PetscFunctionBegin; 821 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 822 PetscCall(ISGetIndices(globalPointNumbers, &gpoint)); 823 for (p = pStart; p < pEnd; ++p) { 824 if (gpoint[p] >= 0) { 825 PetscInt coneSize; 826 827 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 828 nPoints += 1; 829 conesSize += coneSize; 830 } 831 } 832 PetscCall(PetscMalloc1(nPoints, &coneSizes)); 833 PetscCall(PetscMalloc1(conesSize, &cones)); 834 PetscCall(PetscMalloc1(conesSize, &orientations)); 835 for (p = pStart, c = 0, s = 0; p < pEnd; ++p) { 836 if (gpoint[p] >= 0) { 837 const PetscInt *cone, *ornt; 838 PetscInt coneSize, cp; 839 840 PetscCall(DMPlexGetConeSize(dm, p, &coneSize)); 841 PetscCall(DMPlexGetCone(dm, p, &cone)); 842 PetscCall(DMPlexGetConeOrientation(dm, p, &ornt)); 843 coneSizes[s] = coneSize; 844 for (cp = 0; cp < coneSize; ++cp, ++c) { 845 cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]] + 1) : gpoint[cone[cp]]; 846 orientations[c] = ornt[cp]; 847 } 848 ++s; 849 } 850 } 851 PetscCheck(s == nPoints, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %" PetscInt_FMT " != %" PetscInt_FMT, s, nPoints); 852 PetscCheck(c == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %" PetscInt_FMT " != %" PetscInt_FMT, c, conesSize); 853 PetscCall(ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS)); 854 PetscCall(ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS)); 855 PetscCall(ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS)); 856 PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName)); 857 PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName)); 858 PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName)); 859 PetscCall(ISView(coneSizesIS, viewer)); 860 PetscCall(ISView(conesIS, viewer)); 861 PetscCall(ISView(orientationsIS, viewer)); 862 PetscCall(ISDestroy(&coneSizesIS)); 863 PetscCall(ISDestroy(&conesIS)); 864 PetscCall(ISDestroy(&orientationsIS)); 865 if (pointsName) { 866 IS pointsIS; 867 PetscInt *points; 868 869 PetscCall(PetscMalloc1(nPoints, &points)); 870 for (p = pStart, c = 0, s = 0; p < pEnd; ++p) { 871 if (gpoint[p] >= 0) { 872 points[s] = gpoint[p]; 873 ++s; 874 } 875 } 876 PetscCall(ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS)); 877 PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName)); 878 PetscCall(ISView(pointsIS, viewer)); 879 PetscCall(ISDestroy(&pointsIS)); 880 } 881 PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint)); 882 PetscFunctionReturn(PETSC_SUCCESS); 883 } 884 885 static PetscErrorCode DMPlexTopologyView_HDF5_Legacy_Private(DM dm, IS globalPointNumbers, PetscViewer viewer) 886 { 887 const char *pointsName, *coneSizesName, *conesName, *orientationsName; 888 PetscInt pStart, pEnd, dim; 889 890 PetscFunctionBegin; 891 pointsName = "order"; 892 coneSizesName = "cones"; 893 conesName = "cells"; 894 orientationsName = "orientation"; 895 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 896 PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers, viewer, pStart, pEnd, pointsName, coneSizesName, conesName, orientationsName)); 897 PetscCall(DMGetDimension(dm, &dim)); 898 PetscCall(PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *)&dim)); 899 PetscFunctionReturn(PETSC_SUCCESS); 900 } 901 902 //TODO get this numbering right away without needing this function 903 /* Renumber global point numbers so that they are 0-based per stratum */ 904 static PetscErrorCode RenumberGlobalPointNumbersPerStratum_Private(DM dm, IS globalPointNumbers, IS *newGlobalPointNumbers, IS *strataPermutation) 905 { 906 PetscInt d, depth, p, n; 907 PetscInt *offsets; 908 const PetscInt *gpn; 909 PetscInt *ngpn; 910 MPI_Comm comm; 911 912 PetscFunctionBegin; 913 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 914 PetscCall(ISGetLocalSize(globalPointNumbers, &n)); 915 PetscCall(ISGetIndices(globalPointNumbers, &gpn)); 916 PetscCall(PetscMalloc1(n, &ngpn)); 917 PetscCall(DMPlexGetDepth(dm, &depth)); 918 PetscCall(PetscMalloc1(depth + 1, &offsets)); 919 for (d = 0; d <= depth; d++) { 920 PetscInt pStart, pEnd; 921 922 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 923 offsets[d] = PETSC_INT_MAX; 924 for (p = pStart; p < pEnd; p++) { 925 if (gpn[p] >= 0 && gpn[p] < offsets[d]) offsets[d] = gpn[p]; 926 } 927 } 928 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, offsets, depth + 1, MPIU_INT, MPI_MIN, comm)); 929 for (d = 0; d <= depth; d++) { 930 PetscInt pStart, pEnd; 931 932 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 933 for (p = pStart; p < pEnd; p++) ngpn[p] = gpn[p] - PetscSign(gpn[p]) * offsets[d]; 934 } 935 PetscCall(ISRestoreIndices(globalPointNumbers, &gpn)); 936 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)globalPointNumbers), n, ngpn, PETSC_OWN_POINTER, newGlobalPointNumbers)); 937 { 938 PetscInt *perm; 939 940 PetscCall(PetscMalloc1(depth + 1, &perm)); 941 for (d = 0; d <= depth; d++) perm[d] = d; 942 PetscCall(PetscSortIntWithPermutation(depth + 1, offsets, perm)); 943 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, depth + 1, perm, PETSC_OWN_POINTER, strataPermutation)); 944 } 945 PetscCall(PetscFree(offsets)); 946 PetscFunctionReturn(PETSC_SUCCESS); 947 } 948 949 static PetscErrorCode DMPlexTopologyView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer) 950 { 951 IS globalPointNumbers0, strataPermutation; 952 const char *coneSizesName, *conesName, *orientationsName; 953 PetscInt depth, d; 954 MPI_Comm comm; 955 956 PetscFunctionBegin; 957 coneSizesName = "cone_sizes"; 958 conesName = "cones"; 959 orientationsName = "orientations"; 960 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 961 PetscCall(DMPlexGetDepth(dm, &depth)); 962 { 963 PetscInt dim; 964 965 PetscCall(DMGetDimension(dm, &dim)); 966 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "cell_dim", PETSC_INT, &dim)); 967 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "depth", PETSC_INT, &depth)); 968 } 969 970 PetscCall(RenumberGlobalPointNumbersPerStratum_Private(dm, globalPointNumbers, &globalPointNumbers0, &strataPermutation)); 971 /* TODO dirty trick to save serial IS using the same parallel viewer */ 972 { 973 IS spOnComm; 974 PetscInt n = 0, N; 975 const PetscInt *idx = NULL; 976 const PetscInt *old; 977 PetscMPIInt rank; 978 979 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 980 PetscCall(ISGetLocalSize(strataPermutation, &N)); 981 PetscCall(ISGetIndices(strataPermutation, &old)); 982 if (!rank) { 983 n = N; 984 idx = old; 985 } 986 PetscCall(ISCreateGeneral(comm, n, idx, PETSC_COPY_VALUES, &spOnComm)); 987 PetscCall(ISRestoreIndices(strataPermutation, &old)); 988 PetscCall(ISDestroy(&strataPermutation)); 989 strataPermutation = spOnComm; 990 } 991 PetscCall(PetscObjectSetName((PetscObject)strataPermutation, "permutation")); 992 PetscCall(ISView(strataPermutation, viewer)); 993 PetscCall(PetscViewerHDF5PushGroup(viewer, "strata")); 994 for (d = 0; d <= depth; d++) { 995 PetscInt pStart, pEnd; 996 char group[128]; 997 998 PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, d)); 999 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 1000 PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 1001 PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers0, viewer, pStart, pEnd, NULL, coneSizesName, conesName, orientationsName)); 1002 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1003 } 1004 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */ 1005 PetscCall(ISDestroy(&globalPointNumbers0)); 1006 PetscCall(ISDestroy(&strataPermutation)); 1007 PetscFunctionReturn(PETSC_SUCCESS); 1008 } 1009 1010 PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer) 1011 { 1012 DMPlexStorageVersion version; 1013 const char *topologydm_name; 1014 char group[PETSC_MAX_PATH_LEN]; 1015 1016 PetscFunctionBegin; 1017 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version)); 1018 PetscCall(PetscInfo(dm, "Writing DM %s storage version %d.%d.%d\n", dm->hdr.name, version->major, version->minor, version->subminor)); 1019 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 1020 if (DMPlexStorageVersionGE(version, 2, 0, 0)) { 1021 PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name)); 1022 } else { 1023 PetscCall(PetscStrncpy(group, "/", sizeof(group))); 1024 } 1025 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 1026 1027 PetscCall(PetscViewerHDF5PushGroup(viewer, "topology")); 1028 if (version->major < 3) { 1029 PetscCall(DMPlexTopologyView_HDF5_Legacy_Private(dm, globalPointNumbers, viewer)); 1030 } else { 1031 /* since DMPlexStorageVersion 3.0.0 */ 1032 PetscCall(DMPlexTopologyView_HDF5_Private(dm, globalPointNumbers, viewer)); 1033 } 1034 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */ 1035 1036 if (DMPlexStorageVersionGE(version, 2, 1, 0)) { 1037 const char *distribution_name; 1038 1039 PetscCall(DMPlexDistributionGetName(dm, &distribution_name)); 1040 PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions")); 1041 PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL)); 1042 PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name)); 1043 PetscCall(DMPlexDistributionView_HDF5_Private(dm, globalPointNumbers, viewer)); 1044 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */ 1045 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */ 1046 } 1047 1048 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1049 PetscFunctionReturn(PETSC_SUCCESS); 1050 } 1051 1052 static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS) 1053 { 1054 PetscSF sfPoint; 1055 DMLabel cutLabel, cutVertexLabel = NULL; 1056 IS globalVertexNumbers, cutvertices = NULL; 1057 const PetscInt *gcell, *gvertex, *cutverts = NULL; 1058 PetscInt *vertices; 1059 PetscInt conesSize = 0; 1060 PetscInt dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v; 1061 1062 PetscFunctionBegin; 1063 *numCorners = 0; 1064 PetscCall(DMGetDimension(dm, &dim)); 1065 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1066 PetscCall(ISGetIndices(globalCellNumbers, &gcell)); 1067 1068 for (cell = cStart; cell < cEnd; ++cell) { 1069 PetscInt *closure = NULL; 1070 PetscInt closureSize, v, Nc = 0; 1071 1072 if (gcell[cell] < 0) continue; 1073 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 1074 for (v = 0; v < closureSize * 2; v += 2) { 1075 if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc; 1076 } 1077 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 1078 conesSize += Nc; 1079 if (!numCornersLocal) numCornersLocal = Nc; 1080 else if (numCornersLocal != Nc) numCornersLocal = 1; 1081 } 1082 PetscCallMPI(MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm))); 1083 PetscCheck(!numCornersLocal || !(numCornersLocal != *numCorners || *numCorners == 0), PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes"); 1084 /* Handle periodic cuts by identifying vertices which should be duplicated */ 1085 PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel)); 1086 PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel)); 1087 if (cutVertexLabel) PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices)); 1088 if (cutvertices) { 1089 PetscCall(ISGetIndices(cutvertices, &cutverts)); 1090 PetscCall(ISGetLocalSize(cutvertices, &vExtra)); 1091 } 1092 PetscCall(DMGetPointSF(dm, &sfPoint)); 1093 if (cutLabel) { 1094 const PetscInt *ilocal; 1095 const PetscSFNode *iremote; 1096 PetscInt nroots, nleaves; 1097 1098 PetscCall(PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote)); 1099 if (nleaves < 0) { 1100 PetscCall(PetscObjectReference((PetscObject)sfPoint)); 1101 } else { 1102 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), &sfPoint)); 1103 PetscCall(PetscSFSetGraph(sfPoint, nroots + vExtra, nleaves, (PetscInt *)ilocal, PETSC_COPY_VALUES, (PetscSFNode *)iremote, PETSC_COPY_VALUES)); 1104 } 1105 } else { 1106 PetscCall(PetscObjectReference((PetscObject)sfPoint)); 1107 } 1108 /* Number all vertices */ 1109 PetscCall(DMPlexCreateNumbering_Plex(dm, vStart, vEnd + vExtra, 0, NULL, sfPoint, &globalVertexNumbers)); 1110 PetscCall(PetscSFDestroy(&sfPoint)); 1111 /* Create cones */ 1112 PetscCall(ISGetIndices(globalVertexNumbers, &gvertex)); 1113 PetscCall(PetscMalloc1(conesSize, &vertices)); 1114 for (cell = cStart, v = 0; cell < cEnd; ++cell) { 1115 PetscInt *closure = NULL; 1116 PetscInt closureSize, Nc = 0, p, value = -1; 1117 PetscBool replace; 1118 1119 if (gcell[cell] < 0) continue; 1120 if (cutLabel) PetscCall(DMLabelGetValue(cutLabel, cell, &value)); 1121 replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE; 1122 PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 1123 for (p = 0; p < closureSize * 2; p += 2) { 1124 if ((closure[p] >= vStart) && (closure[p] < vEnd)) closure[Nc++] = closure[p]; 1125 } 1126 PetscCall(DMPlexReorderCell(dm, cell, closure)); 1127 for (p = 0; p < Nc; ++p) { 1128 PetscInt nv, gv = gvertex[closure[p] - vStart]; 1129 1130 if (replace) { 1131 PetscCall(PetscFindInt(closure[p], vExtra, cutverts, &nv)); 1132 if (nv >= 0) gv = gvertex[vEnd - vStart + nv]; 1133 } 1134 vertices[v++] = gv < 0 ? -(gv + 1) : gv; 1135 } 1136 PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure)); 1137 } 1138 PetscCall(ISRestoreIndices(globalVertexNumbers, &gvertex)); 1139 PetscCall(ISDestroy(&globalVertexNumbers)); 1140 PetscCall(ISRestoreIndices(globalCellNumbers, &gcell)); 1141 if (cutvertices) PetscCall(ISRestoreIndices(cutvertices, &cutverts)); 1142 PetscCall(ISDestroy(&cutvertices)); 1143 PetscCall(DMLabelDestroy(&cutVertexLabel)); 1144 PetscCheck(v == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %" PetscInt_FMT " != %" PetscInt_FMT, v, conesSize); 1145 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS)); 1146 PetscCall(PetscLayoutSetBlockSize((*cellIS)->map, *numCorners)); 1147 PetscCall(PetscObjectSetName((PetscObject)*cellIS, "cells")); 1148 PetscFunctionReturn(PETSC_SUCCESS); 1149 } 1150 1151 static PetscErrorCode DMPlexTopologyView_HDF5_XDMF_Private(DM dm, IS globalCellNumbers, PetscViewer viewer) 1152 { 1153 DM cdm; 1154 DMLabel depthLabel, ctLabel; 1155 IS cellIS; 1156 PetscInt dim, depth, cellHeight, c, n = 0; 1157 1158 PetscFunctionBegin; 1159 PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz")); 1160 PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL)); 1161 1162 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1163 PetscCall(DMGetDimension(dm, &dim)); 1164 PetscCall(DMPlexGetDepth(dm, &depth)); 1165 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1166 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 1167 PetscCall(DMPlexGetDepthLabel(dm, &depthLabel)); 1168 PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel)); 1169 for (c = 0; c < DM_NUM_POLYTOPES; ++c) { 1170 const DMPolytopeType ict = (DMPolytopeType)c; 1171 PetscInt pStart, pEnd, dep, numCorners; 1172 PetscBool output = PETSC_FALSE, doOutput; 1173 1174 if (ict == DM_POLYTOPE_FV_GHOST) continue; 1175 PetscCall(DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd)); 1176 if (pStart >= 0) { 1177 PetscCall(DMLabelGetValue(depthLabel, pStart, &dep)); 1178 if (dep == depth - cellHeight) output = PETSC_TRUE; 1179 } 1180 PetscCallMPI(MPIU_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm))); 1181 if (!doOutput) continue; 1182 PetscCall(CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS)); 1183 if (!n) { 1184 PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/topology")); 1185 } else { 1186 char group[PETSC_MAX_PATH_LEN]; 1187 1188 PetscCall(PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%" PetscInt_FMT, n)); 1189 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 1190 } 1191 PetscCall(ISView(cellIS, viewer)); 1192 PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_corners", PETSC_INT, (void *)&numCorners)); 1193 PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_dim", PETSC_INT, (void *)&dim)); 1194 PetscCall(ISDestroy(&cellIS)); 1195 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1196 ++n; 1197 } 1198 PetscFunctionReturn(PETSC_SUCCESS); 1199 } 1200 1201 static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer) 1202 { 1203 DM cdm; 1204 Vec coordinates, newcoords; 1205 PetscReal lengthScale; 1206 PetscInt m, M, bs; 1207 1208 PetscFunctionBegin; 1209 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 1210 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1211 PetscCall(DMGetCoordinates(dm, &coordinates)); 1212 PetscCall(VecCreate(PetscObjectComm((PetscObject)coordinates), &newcoords)); 1213 PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices")); 1214 PetscCall(VecGetSize(coordinates, &M)); 1215 PetscCall(VecGetLocalSize(coordinates, &m)); 1216 PetscCall(VecSetSizes(newcoords, m, M)); 1217 PetscCall(VecGetBlockSize(coordinates, &bs)); 1218 PetscCall(VecSetBlockSize(newcoords, bs)); 1219 PetscCall(VecSetType(newcoords, VECSTANDARD)); 1220 PetscCall(VecCopy(coordinates, newcoords)); 1221 PetscCall(VecScale(newcoords, lengthScale)); 1222 /* Did not use DMGetGlobalVector() in order to bypass default group assignment */ 1223 PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry")); 1224 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE)); 1225 PetscCall(VecView(newcoords, viewer)); 1226 PetscCall(PetscViewerPopFormat(viewer)); 1227 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1228 PetscCall(VecDestroy(&newcoords)); 1229 PetscFunctionReturn(PETSC_SUCCESS); 1230 } 1231 1232 PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer) 1233 { 1234 DM cdm; 1235 Vec coords, newcoords; 1236 PetscInt m, M, bs; 1237 PetscReal lengthScale; 1238 PetscBool viewSection = PETSC_TRUE; 1239 const char *topologydm_name, *coordinatedm_name, *coordinates_name; 1240 1241 PetscFunctionBegin; 1242 { 1243 PetscViewerFormat format; 1244 DMPlexStorageVersion version; 1245 1246 PetscCall(PetscViewerGetFormat(viewer, &format)); 1247 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version)); 1248 if (!DMPlexStorageVersionGE(version, 2, 0, 0) || format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ) { 1249 PetscCall(DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer)); 1250 PetscFunctionReturn(PETSC_SUCCESS); 1251 } 1252 } 1253 /* since 2.0.0 */ 1254 PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_view_coordinate_section", &viewSection, NULL)); 1255 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1256 PetscCall(DMGetCoordinates(dm, &coords)); 1257 PetscCall(PetscObjectGetName((PetscObject)cdm, &coordinatedm_name)); 1258 PetscCall(PetscObjectGetName((PetscObject)coords, &coordinates_name)); 1259 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 1260 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 1261 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 1262 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name)); 1263 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name)); 1264 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1265 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1266 if (viewSection) PetscCall(DMPlexSectionView(dm, viewer, cdm)); 1267 PetscCall(VecCreate(PetscObjectComm((PetscObject)coords), &newcoords)); 1268 PetscCall(PetscObjectSetName((PetscObject)newcoords, coordinates_name)); 1269 PetscCall(VecGetSize(coords, &M)); 1270 PetscCall(VecGetLocalSize(coords, &m)); 1271 PetscCall(VecSetSizes(newcoords, m, M)); 1272 PetscCall(VecGetBlockSize(coords, &bs)); 1273 PetscCall(VecSetBlockSize(newcoords, bs)); 1274 PetscCall(VecSetType(newcoords, VECSTANDARD)); 1275 PetscCall(VecCopy(coords, newcoords)); 1276 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 1277 PetscCall(VecScale(newcoords, lengthScale)); 1278 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE)); 1279 PetscCall(DMPlexGlobalVectorView(dm, viewer, cdm, newcoords)); 1280 PetscCall(PetscViewerPopFormat(viewer)); 1281 PetscCall(VecDestroy(&newcoords)); 1282 PetscFunctionReturn(PETSC_SUCCESS); 1283 } 1284 1285 static PetscErrorCode DMPlexCoordinatesView_HDF5_XDMF_Private(DM dm, PetscViewer viewer) 1286 { 1287 DM cdm; 1288 Vec coordinatesLocal, newcoords; 1289 PetscSection cSection, cGlobalSection; 1290 PetscScalar *coords, *ncoords; 1291 DMLabel cutLabel, cutVertexLabel = NULL; 1292 const PetscReal *L; 1293 PetscReal lengthScale; 1294 PetscInt vStart, vEnd, v, bs, N, coordSize, dof, off, d; 1295 PetscBool localized, embedded; 1296 1297 PetscFunctionBegin; 1298 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1299 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 1300 PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 1301 PetscCall(VecGetBlockSize(coordinatesLocal, &bs)); 1302 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 1303 if (localized == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS); 1304 PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L)); 1305 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1306 PetscCall(DMGetLocalSection(cdm, &cSection)); 1307 PetscCall(DMGetGlobalSection(cdm, &cGlobalSection)); 1308 PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel)); 1309 N = 0; 1310 1311 PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel)); 1312 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &newcoords)); 1313 PetscCall(PetscSectionGetDof(cSection, vStart, &dof)); 1314 PetscCall(PetscPrintf(PETSC_COMM_SELF, "DOF: %" PetscInt_FMT "\n", dof)); 1315 embedded = (PetscBool)(L && dof == 2 && !cutLabel); 1316 if (cutVertexLabel) { 1317 PetscCall(DMLabelGetStratumSize(cutVertexLabel, 1, &v)); 1318 N += dof * v; 1319 } 1320 for (v = vStart; v < vEnd; ++v) { 1321 PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof)); 1322 if (dof < 0) continue; 1323 if (embedded) N += dof + 1; 1324 else N += dof; 1325 } 1326 if (embedded) PetscCall(VecSetBlockSize(newcoords, bs + 1)); 1327 else PetscCall(VecSetBlockSize(newcoords, bs)); 1328 PetscCall(VecSetSizes(newcoords, N, PETSC_DETERMINE)); 1329 PetscCall(VecSetType(newcoords, VECSTANDARD)); 1330 PetscCall(VecGetArray(coordinatesLocal, &coords)); 1331 PetscCall(VecGetArray(newcoords, &ncoords)); 1332 coordSize = 0; 1333 for (v = vStart; v < vEnd; ++v) { 1334 PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof)); 1335 PetscCall(PetscSectionGetOffset(cSection, v, &off)); 1336 if (dof < 0) continue; 1337 if (embedded) { 1338 if (L && (L[0] > 0.0) && (L[1] > 0.0)) { 1339 PetscReal theta, phi, r, R; 1340 /* XY-periodic */ 1341 /* Suppose its an y-z circle, then 1342 \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0) 1343 and the circle in that plane is 1344 \hat r cos(phi) + \hat x sin(phi) */ 1345 theta = 2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]; 1346 phi = 2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]; 1347 r = L[0] / (2.0 * PETSC_PI * 2.0 * L[1]); 1348 R = L[1] / (2.0 * PETSC_PI); 1349 ncoords[coordSize++] = PetscSinReal(phi) * r; 1350 ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi)); 1351 ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi)); 1352 } else if (L && (L[0] > 0.0)) { 1353 /* X-periodic */ 1354 ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI)); 1355 ncoords[coordSize++] = coords[off + 1]; 1356 ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI)); 1357 } else if (L && (L[1] > 0.0)) { 1358 /* Y-periodic */ 1359 ncoords[coordSize++] = coords[off + 0]; 1360 ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI)); 1361 ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI)); 1362 #if 0 1363 } else if ((bd[0] == DM_BOUNDARY_TWIST)) { 1364 PetscReal phi, r, R; 1365 /* Mobius strip */ 1366 /* Suppose its an x-z circle, then 1367 \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0) 1368 and in that plane we rotate by pi as we go around the circle 1369 \hat r cos(phi/2) + \hat y sin(phi/2) */ 1370 phi = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0]; 1371 R = L[0]; 1372 r = PetscRealPart(coords[off+1]) - L[1]/2.0; 1373 ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0)); 1374 ncoords[coordSize++] = PetscSinReal(phi/2.0) * r; 1375 ncoords[coordSize++] = PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0)); 1376 #endif 1377 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain"); 1378 } else { 1379 for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off + d]; 1380 } 1381 } 1382 if (cutVertexLabel) { 1383 IS vertices; 1384 const PetscInt *verts; 1385 PetscInt n; 1386 1387 PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &vertices)); 1388 if (vertices) { 1389 PetscCall(ISGetIndices(vertices, &verts)); 1390 PetscCall(ISGetLocalSize(vertices, &n)); 1391 for (v = 0; v < n; ++v) { 1392 PetscCall(PetscSectionGetDof(cSection, verts[v], &dof)); 1393 PetscCall(PetscSectionGetOffset(cSection, verts[v], &off)); 1394 for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off + d] + ((L[d] > 0.) ? L[d] : 0.0); 1395 } 1396 PetscCall(ISRestoreIndices(vertices, &verts)); 1397 PetscCall(ISDestroy(&vertices)); 1398 } 1399 } 1400 PetscCheck(coordSize == N, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %" PetscInt_FMT " != %" PetscInt_FMT, coordSize, N); 1401 PetscCall(DMLabelDestroy(&cutVertexLabel)); 1402 PetscCall(VecRestoreArray(coordinatesLocal, &coords)); 1403 PetscCall(VecRestoreArray(newcoords, &ncoords)); 1404 PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices")); 1405 PetscCall(VecScale(newcoords, lengthScale)); 1406 PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz")); 1407 PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL)); 1408 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1409 PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/geometry")); 1410 PetscCall(VecView(newcoords, viewer)); 1411 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1412 PetscCall(VecDestroy(&newcoords)); 1413 PetscFunctionReturn(PETSC_SUCCESS); 1414 } 1415 1416 PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer) 1417 { 1418 const char *topologydm_name; 1419 const PetscInt *gpoint; 1420 PetscInt numLabels; 1421 PetscBool omitCelltypes = PETSC_FALSE; 1422 DMPlexStorageVersion version; 1423 char group[PETSC_MAX_PATH_LEN]; 1424 1425 PetscFunctionBegin; 1426 PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_omit_celltypes", &omitCelltypes, NULL)); 1427 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version)); 1428 PetscCall(ISGetIndices(globalPointNumbers, &gpoint)); 1429 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 1430 if (DMPlexStorageVersionGE(version, 2, 0, 0)) { 1431 PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name)); 1432 } else { 1433 PetscCall(PetscStrncpy(group, "/labels", sizeof(group))); 1434 } 1435 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 1436 PetscCall(DMGetNumLabels(dm, &numLabels)); 1437 for (PetscInt l = 0; l < numLabels; ++l) { 1438 DMLabel label; 1439 const char *name; 1440 IS valueIS, pvalueIS, globalValueIS; 1441 const PetscInt *values; 1442 PetscInt numValues, v; 1443 PetscBool isDepth, isCelltype, output; 1444 1445 PetscCall(DMGetLabelByNum(dm, l, &label)); 1446 PetscCall(PetscObjectGetName((PetscObject)label, &name)); 1447 PetscCall(DMGetLabelOutput(dm, name, &output)); 1448 PetscCall(PetscStrncmp(name, "depth", 10, &isDepth)); 1449 PetscCall(PetscStrncmp(name, "celltype", 10, &isCelltype)); 1450 // TODO Should only filter out celltype if it can be calculated 1451 if (isDepth || (isCelltype && omitCelltypes) || !output) continue; 1452 PetscCall(PetscViewerHDF5PushGroup(viewer, name)); 1453 PetscCall(DMLabelGetValueIS(label, &valueIS)); 1454 /* Must copy to a new IS on the global comm */ 1455 PetscCall(ISGetLocalSize(valueIS, &numValues)); 1456 PetscCall(ISGetIndices(valueIS, &values)); 1457 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS)); 1458 PetscCall(ISRestoreIndices(valueIS, &values)); 1459 PetscCall(ISAllGather(pvalueIS, &globalValueIS)); 1460 PetscCall(ISDestroy(&pvalueIS)); 1461 PetscCall(ISSortRemoveDups(globalValueIS)); 1462 PetscCall(ISGetLocalSize(globalValueIS, &numValues)); 1463 PetscCall(ISGetIndices(globalValueIS, &values)); 1464 for (v = 0; v < numValues; ++v) { 1465 IS stratumIS, globalStratumIS; 1466 const PetscInt *spoints = NULL; 1467 PetscInt *gspoints, n = 0, gn, p; 1468 const char *iname = "indices"; 1469 char group[PETSC_MAX_PATH_LEN]; 1470 1471 PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, values[v])); 1472 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 1473 PetscCall(DMLabelGetStratumIS(label, values[v], &stratumIS)); 1474 1475 if (stratumIS) PetscCall(ISGetLocalSize(stratumIS, &n)); 1476 if (stratumIS) PetscCall(ISGetIndices(stratumIS, &spoints)); 1477 for (gn = 0, p = 0; p < n; ++p) 1478 if (gpoint[spoints[p]] >= 0) ++gn; 1479 PetscCall(PetscMalloc1(gn, &gspoints)); 1480 for (gn = 0, p = 0; p < n; ++p) 1481 if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]]; 1482 if (stratumIS) PetscCall(ISRestoreIndices(stratumIS, &spoints)); 1483 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS)); 1484 PetscCall(PetscObjectSetName((PetscObject)globalStratumIS, iname)); 1485 1486 PetscCall(ISView(globalStratumIS, viewer)); 1487 PetscCall(ISDestroy(&globalStratumIS)); 1488 PetscCall(ISDestroy(&stratumIS)); 1489 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1490 } 1491 PetscCall(ISRestoreIndices(globalValueIS, &values)); 1492 PetscCall(ISDestroy(&globalValueIS)); 1493 PetscCall(ISDestroy(&valueIS)); 1494 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1495 } 1496 PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint)); 1497 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1498 PetscFunctionReturn(PETSC_SUCCESS); 1499 } 1500 1501 /* We only write cells and vertices. Does this screw up parallel reading? */ 1502 PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer) 1503 { 1504 IS globalPointNumbers; 1505 PetscViewerFormat format; 1506 PetscBool viz_geom = PETSC_FALSE, xdmf_topo = PETSC_FALSE, petsc_topo = PETSC_FALSE; 1507 1508 PetscFunctionBegin; 1509 PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers)); 1510 PetscCall(DMPlexCoordinatesView_HDF5_Internal(dm, viewer)); 1511 1512 PetscCall(PetscViewerGetFormat(viewer, &format)); 1513 switch (format) { 1514 case PETSC_VIEWER_HDF5_VIZ: 1515 viz_geom = PETSC_TRUE; 1516 xdmf_topo = PETSC_TRUE; 1517 break; 1518 case PETSC_VIEWER_HDF5_XDMF: 1519 xdmf_topo = PETSC_TRUE; 1520 break; 1521 case PETSC_VIEWER_HDF5_PETSC: 1522 petsc_topo = PETSC_TRUE; 1523 break; 1524 case PETSC_VIEWER_DEFAULT: 1525 case PETSC_VIEWER_NATIVE: 1526 viz_geom = PETSC_TRUE; 1527 xdmf_topo = PETSC_TRUE; 1528 petsc_topo = PETSC_TRUE; 1529 break; 1530 default: 1531 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]); 1532 } 1533 1534 if (viz_geom) PetscCall(DMPlexCoordinatesView_HDF5_XDMF_Private(dm, viewer)); 1535 if (xdmf_topo) PetscCall(DMPlexTopologyView_HDF5_XDMF_Private(dm, globalPointNumbers, viewer)); 1536 if (petsc_topo) { 1537 PetscBool viewLabels = PETSC_TRUE; 1538 1539 PetscCall(DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer)); 1540 PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_view_labels", &viewLabels, NULL)); 1541 if (viewLabels) PetscCall(DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer)); 1542 } 1543 1544 PetscCall(ISDestroy(&globalPointNumbers)); 1545 PetscFunctionReturn(PETSC_SUCCESS); 1546 } 1547 1548 PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm) 1549 { 1550 MPI_Comm comm; 1551 const char *topologydm_name; 1552 const char *sectiondm_name; 1553 PetscSection gsection; 1554 1555 PetscFunctionBegin; 1556 PetscCall(PetscObjectGetComm((PetscObject)sectiondm, &comm)); 1557 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 1558 PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name)); 1559 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 1560 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 1561 PetscCall(PetscViewerHDF5PushGroup(viewer, "dms")); 1562 PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name)); 1563 PetscCall(DMGetGlobalSection(sectiondm, &gsection)); 1564 /* Save raw section */ 1565 PetscCall(PetscSectionView(gsection, viewer)); 1566 /* Save plex wrapper */ 1567 { 1568 PetscInt pStart, pEnd, p, n; 1569 IS globalPointNumbers; 1570 const PetscInt *gpoints; 1571 IS orderIS; 1572 PetscInt *order; 1573 1574 PetscCall(PetscSectionGetChart(gsection, &pStart, &pEnd)); 1575 PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers)); 1576 PetscCall(ISGetIndices(globalPointNumbers, &gpoints)); 1577 for (p = pStart, n = 0; p < pEnd; ++p) 1578 if (gpoints[p] >= 0) n++; 1579 /* "order" is an array of global point numbers. 1580 When loading, it is used with topology/order array 1581 to match section points with plex topology points. */ 1582 PetscCall(PetscMalloc1(n, &order)); 1583 for (p = pStart, n = 0; p < pEnd; ++p) 1584 if (gpoints[p] >= 0) order[n++] = gpoints[p]; 1585 PetscCall(ISRestoreIndices(globalPointNumbers, &gpoints)); 1586 PetscCall(ISDestroy(&globalPointNumbers)); 1587 PetscCall(ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS)); 1588 PetscCall(PetscObjectSetName((PetscObject)orderIS, "order")); 1589 PetscCall(ISView(orderIS, viewer)); 1590 PetscCall(ISDestroy(&orderIS)); 1591 } 1592 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1593 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1594 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1595 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1596 PetscFunctionReturn(PETSC_SUCCESS); 1597 } 1598 1599 PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec) 1600 { 1601 const char *topologydm_name; 1602 const char *sectiondm_name; 1603 const char *vec_name; 1604 PetscInt bs; 1605 1606 PetscFunctionBegin; 1607 /* Check consistency */ 1608 { 1609 PetscSF pointsf, pointsf1; 1610 1611 PetscCall(DMGetPointSF(dm, &pointsf)); 1612 PetscCall(DMGetPointSF(sectiondm, &pointsf1)); 1613 PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm"); 1614 } 1615 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 1616 PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name)); 1617 PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name)); 1618 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 1619 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 1620 PetscCall(PetscViewerHDF5PushGroup(viewer, "dms")); 1621 PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name)); 1622 PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs")); 1623 PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name)); 1624 PetscCall(VecGetBlockSize(vec, &bs)); 1625 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs)); 1626 PetscCall(VecSetBlockSize(vec, 1)); 1627 /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but, */ 1628 /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view */ 1629 /* is set to VecView_Plex, which would save vec in a predefined location. */ 1630 /* To save vec in where we want, we create a new Vec (temp) with */ 1631 /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */ 1632 { 1633 Vec temp; 1634 const PetscScalar *array; 1635 PetscLayout map; 1636 1637 PetscCall(VecCreate(PetscObjectComm((PetscObject)vec), &temp)); 1638 PetscCall(PetscObjectSetName((PetscObject)temp, vec_name)); 1639 PetscCall(VecGetLayout(vec, &map)); 1640 PetscCall(VecSetLayout(temp, map)); 1641 PetscCall(VecSetUp(temp)); 1642 PetscCall(VecGetArrayRead(vec, &array)); 1643 PetscCall(VecPlaceArray(temp, array)); 1644 PetscCall(VecView(temp, viewer)); 1645 PetscCall(VecResetArray(temp)); 1646 PetscCall(VecRestoreArrayRead(vec, &array)); 1647 PetscCall(VecDestroy(&temp)); 1648 } 1649 PetscCall(VecSetBlockSize(vec, bs)); 1650 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1651 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1652 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1653 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1654 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1655 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1656 PetscFunctionReturn(PETSC_SUCCESS); 1657 } 1658 1659 PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec) 1660 { 1661 MPI_Comm comm; 1662 const char *topologydm_name; 1663 const char *sectiondm_name; 1664 const char *vec_name; 1665 PetscSection section; 1666 PetscBool includesConstraints; 1667 Vec gvec; 1668 PetscInt m, bs; 1669 1670 PetscFunctionBegin; 1671 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1672 /* Check consistency */ 1673 { 1674 PetscSF pointsf, pointsf1; 1675 1676 PetscCall(DMGetPointSF(dm, &pointsf)); 1677 PetscCall(DMGetPointSF(sectiondm, &pointsf1)); 1678 PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm"); 1679 } 1680 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 1681 PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name)); 1682 PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name)); 1683 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 1684 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 1685 PetscCall(PetscViewerHDF5PushGroup(viewer, "dms")); 1686 PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name)); 1687 PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs")); 1688 PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name)); 1689 PetscCall(VecGetBlockSize(vec, &bs)); 1690 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs)); 1691 PetscCall(VecCreate(comm, &gvec)); 1692 PetscCall(PetscObjectSetName((PetscObject)gvec, vec_name)); 1693 PetscCall(DMGetGlobalSection(sectiondm, §ion)); 1694 PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints)); 1695 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m)); 1696 else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m)); 1697 PetscCall(VecSetSizes(gvec, m, PETSC_DECIDE)); 1698 PetscCall(VecSetUp(gvec)); 1699 PetscCall(DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec)); 1700 PetscCall(DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec)); 1701 PetscCall(VecView(gvec, viewer)); 1702 PetscCall(VecDestroy(&gvec)); 1703 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1704 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1705 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1706 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1707 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1708 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1709 PetscFunctionReturn(PETSC_SUCCESS); 1710 } 1711 1712 struct _n_LoadLabelsCtx { 1713 MPI_Comm comm; 1714 PetscMPIInt rank; 1715 DM dm; 1716 PetscViewer viewer; 1717 DMLabel label; 1718 PetscSF sfXC; 1719 PetscLayout layoutX; 1720 }; 1721 typedef struct _n_LoadLabelsCtx *LoadLabelsCtx; 1722 1723 static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx) 1724 { 1725 PetscFunctionBegin; 1726 PetscCall(PetscNew(ctx)); 1727 PetscCall(PetscObjectReference((PetscObject)((*ctx)->dm = dm))); 1728 PetscCall(PetscObjectReference((PetscObject)((*ctx)->viewer = viewer))); 1729 PetscCall(PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm)); 1730 PetscCallMPI(MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank)); 1731 (*ctx)->sfXC = sfXC; 1732 if (sfXC) { 1733 PetscInt nX; 1734 1735 PetscCall(PetscObjectReference((PetscObject)sfXC)); 1736 PetscCall(PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL)); 1737 PetscCall(PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX)); 1738 } 1739 PetscFunctionReturn(PETSC_SUCCESS); 1740 } 1741 1742 static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx) 1743 { 1744 PetscFunctionBegin; 1745 if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS); 1746 PetscCall(DMDestroy(&(*ctx)->dm)); 1747 PetscCall(PetscViewerDestroy(&(*ctx)->viewer)); 1748 PetscCall(PetscSFDestroy(&(*ctx)->sfXC)); 1749 PetscCall(PetscLayoutDestroy(&(*ctx)->layoutX)); 1750 PetscCall(PetscFree(*ctx)); 1751 PetscFunctionReturn(PETSC_SUCCESS); 1752 } 1753 1754 /* 1755 A: on-disk points 1756 X: global points [0, NX) 1757 C: distributed plex points 1758 */ 1759 static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS) 1760 { 1761 MPI_Comm comm = ctx->comm; 1762 PetscSF sfXC = ctx->sfXC; 1763 PetscLayout layoutX = ctx->layoutX; 1764 PetscSF sfXA; 1765 const PetscInt *A_points; 1766 PetscInt nX, nC; 1767 PetscInt n; 1768 1769 PetscFunctionBegin; 1770 PetscCall(PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL)); 1771 PetscCall(ISGetLocalSize(stratumIS, &n)); 1772 PetscCall(ISGetIndices(stratumIS, &A_points)); 1773 PetscCall(PetscSFCreate(comm, &sfXA)); 1774 PetscCall(PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points)); 1775 PetscCall(ISCreate(comm, newStratumIS)); 1776 PetscCall(ISSetType(*newStratumIS, ISGENERAL)); 1777 { 1778 PetscInt i; 1779 PetscBool *A_mask, *X_mask, *C_mask; 1780 1781 PetscCall(PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask)); 1782 for (i = 0; i < n; i++) A_mask[i] = PETSC_TRUE; 1783 PetscCall(PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE)); 1784 PetscCall(PetscSFReduceEnd(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE)); 1785 PetscCall(PetscSFBcastBegin(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR)); 1786 PetscCall(PetscSFBcastEnd(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR)); 1787 PetscCall(ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask)); 1788 PetscCall(PetscFree3(A_mask, X_mask, C_mask)); 1789 } 1790 PetscCall(PetscSFDestroy(&sfXA)); 1791 PetscCall(ISRestoreIndices(stratumIS, &A_points)); 1792 PetscFunctionReturn(PETSC_SUCCESS); 1793 } 1794 1795 static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data) 1796 { 1797 LoadLabelsCtx ctx = (LoadLabelsCtx)op_data; 1798 PetscViewer viewer = ctx->viewer; 1799 DMLabel label = ctx->label; 1800 MPI_Comm comm = ctx->comm; 1801 IS stratumIS; 1802 const PetscInt *ind; 1803 PetscInt value, N, i; 1804 1805 PetscCall(PetscOptionsStringToInt(vname, &value)); 1806 PetscCall(ISCreate(comm, &stratumIS)); 1807 PetscCall(PetscObjectSetName((PetscObject)stratumIS, "indices")); 1808 PetscCall(PetscViewerHDF5PushGroup(viewer, vname)); /* labels/<lname>/<vname> */ 1809 1810 if (!ctx->sfXC) { 1811 /* Force serial load */ 1812 PetscCall(PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N)); 1813 PetscCall(PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0)); 1814 PetscCall(PetscLayoutSetSize(stratumIS->map, N)); 1815 } 1816 PetscCall(ISLoad(stratumIS, viewer)); 1817 1818 if (ctx->sfXC) { 1819 IS newStratumIS; 1820 1821 PetscCallHDF5(ReadLabelStratumHDF5_Distribute_Private, (stratumIS, ctx, &newStratumIS)); 1822 PetscCall(ISDestroy(&stratumIS)); 1823 stratumIS = newStratumIS; 1824 } 1825 1826 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1827 PetscCall(ISGetLocalSize(stratumIS, &N)); 1828 PetscCall(ISGetIndices(stratumIS, &ind)); 1829 for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, ind[i], value)); 1830 PetscCall(ISRestoreIndices(stratumIS, &ind)); 1831 PetscCall(ISDestroy(&stratumIS)); 1832 return 0; 1833 } 1834 1835 /* TODO: Fix this code, it is returning PETSc error codes when it should be translating them to herr_t codes */ 1836 static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data) 1837 { 1838 LoadLabelsCtx ctx = (LoadLabelsCtx)op_data; 1839 DM dm = ctx->dm; 1840 hsize_t idx = 0; 1841 PetscErrorCode ierr; 1842 PetscBool flg; 1843 herr_t err; 1844 1845 PetscCall(DMHasLabel(dm, lname, &flg)); 1846 if (flg) PetscCall(DMRemoveLabel(dm, lname, NULL)); 1847 ierr = DMCreateLabel(dm, lname); 1848 if (ierr) return (herr_t)ierr; 1849 ierr = DMGetLabel(dm, lname, &ctx->label); 1850 if (ierr) return (herr_t)ierr; 1851 ierr = PetscViewerHDF5PushGroup(ctx->viewer, lname); 1852 if (ierr) return (herr_t)ierr; 1853 /* Iterate over the label's strata */ 1854 PetscCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0)); 1855 ierr = PetscViewerHDF5PopGroup(ctx->viewer); 1856 if (ierr) return (herr_t)ierr; 1857 return err; 1858 } 1859 1860 PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC) 1861 { 1862 const char *topologydm_name; 1863 LoadLabelsCtx ctx; 1864 hsize_t idx = 0; 1865 char group[PETSC_MAX_PATH_LEN]; 1866 DMPlexStorageVersion version; 1867 PetscBool distributed, hasGroup; 1868 1869 PetscFunctionBegin; 1870 PetscCall(DMPlexIsDistributed(dm, &distributed)); 1871 if (distributed) PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load"); 1872 PetscCall(LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx)); 1873 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 1874 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version)); 1875 if (DMPlexStorageVersionGE(version, 2, 0, 0)) { 1876 PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name)); 1877 } else { 1878 PetscCall(PetscStrncpy(group, "labels", sizeof(group))); 1879 } 1880 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 1881 PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup)); 1882 if (hasGroup) { 1883 hid_t fileId, groupId; 1884 1885 PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &fileId, &groupId)); 1886 /* Iterate over labels */ 1887 PetscCallHDF5(H5Literate, (groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx)); 1888 PetscCallHDF5(H5Gclose, (groupId)); 1889 } 1890 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1891 PetscCall(LoadLabelsCtxDestroy(&ctx)); 1892 PetscFunctionReturn(PETSC_SUCCESS); 1893 } 1894 1895 static PetscErrorCode DMPlexDistributionLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF sf, PetscSF *distsf, DM *distdm) 1896 { 1897 MPI_Comm comm; 1898 PetscMPIInt size, rank; 1899 PetscInt dist_size; 1900 const char *distribution_name; 1901 PetscInt p, lsize; 1902 IS chartSizesIS, ownersIS, gpointsIS; 1903 const PetscInt *chartSize, *owners, *gpoints; 1904 PetscLayout layout; 1905 PetscBool has; 1906 1907 PetscFunctionBegin; 1908 *distsf = NULL; 1909 *distdm = NULL; 1910 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 1911 PetscCallMPI(MPI_Comm_size(comm, &size)); 1912 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1913 PetscCall(DMPlexDistributionGetName(dm, &distribution_name)); 1914 if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS); 1915 PetscCall(PetscLogEventBegin(DMPLEX_DistributionLoad, viewer, 0, 0, 0)); 1916 PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &has)); 1917 if (!has) { 1918 char *full_group; 1919 1920 PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &full_group)); 1921 PetscCheck(has, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Distribution %s cannot be found: HDF5 group %s not found in file", distribution_name, full_group); 1922 } 1923 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "comm_size", PETSC_INT, NULL, (void *)&dist_size)); 1924 PetscCheck(dist_size == (PetscInt)size, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mismatching comm sizes: comm size of this session (%d) != comm size used for %s (%" PetscInt_FMT ")", size, distribution_name, dist_size); 1925 PetscCall(ISCreate(comm, &chartSizesIS)); 1926 PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes")); 1927 PetscCall(ISCreate(comm, &ownersIS)); 1928 PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners")); 1929 PetscCall(ISCreate(comm, &gpointsIS)); 1930 PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers")); 1931 PetscCall(PetscLayoutSetLocalSize(chartSizesIS->map, 1)); 1932 PetscCall(ISLoad(chartSizesIS, viewer)); 1933 PetscCall(ISGetIndices(chartSizesIS, &chartSize)); 1934 PetscCall(PetscLayoutSetLocalSize(ownersIS->map, *chartSize)); 1935 PetscCall(PetscLayoutSetLocalSize(gpointsIS->map, *chartSize)); 1936 PetscCall(ISLoad(ownersIS, viewer)); 1937 PetscCall(ISLoad(gpointsIS, viewer)); 1938 PetscCall(ISGetIndices(ownersIS, &owners)); 1939 PetscCall(ISGetIndices(gpointsIS, &gpoints)); 1940 PetscCall(PetscSFCreate(comm, distsf)); 1941 PetscCall(PetscSFSetFromOptions(*distsf)); 1942 PetscCall(PetscLayoutCreate(comm, &layout)); 1943 PetscCall(PetscSFGetGraph(sf, &lsize, NULL, NULL, NULL)); 1944 PetscCall(PetscLayoutSetLocalSize(layout, lsize)); 1945 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 1946 PetscCall(PetscLayoutSetUp(layout)); 1947 PetscCall(PetscSFSetGraphLayout(*distsf, layout, *chartSize, NULL, PETSC_OWN_POINTER, gpoints)); 1948 PetscCall(PetscLayoutDestroy(&layout)); 1949 /* Migrate DM */ 1950 { 1951 PetscInt pStart, pEnd; 1952 PetscSFNode *buffer0, *buffer1, *buffer2; 1953 1954 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1955 PetscCall(PetscMalloc2(pEnd - pStart, &buffer0, lsize, &buffer1)); 1956 PetscCall(PetscMalloc1(*chartSize, &buffer2)); 1957 { 1958 PetscSF workPointSF; 1959 PetscInt workNroots, workNleaves; 1960 const PetscInt *workIlocal; 1961 const PetscSFNode *workIremote; 1962 1963 for (p = pStart; p < pEnd; ++p) { 1964 buffer0[p - pStart].rank = rank; 1965 buffer0[p - pStart].index = p - pStart; 1966 } 1967 PetscCall(DMGetPointSF(dm, &workPointSF)); 1968 PetscCall(PetscSFGetGraph(workPointSF, &workNroots, &workNleaves, &workIlocal, &workIremote)); 1969 for (p = 0; p < workNleaves; ++p) { 1970 PetscInt workIlocalp = (workIlocal ? workIlocal[p] : p); 1971 1972 buffer0[workIlocalp].rank = -1; 1973 } 1974 } 1975 for (p = 0; p < lsize; ++p) buffer1[p].rank = -1; 1976 for (p = 0; p < *chartSize; ++p) buffer2[p].rank = -1; 1977 PetscCall(PetscSFReduceBegin(sf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC)); 1978 PetscCall(PetscSFReduceEnd(sf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC)); 1979 PetscCall(PetscSFBcastBegin(*distsf, MPIU_SF_NODE, buffer1, buffer2, MPI_REPLACE)); 1980 PetscCall(PetscSFBcastEnd(*distsf, MPIU_SF_NODE, buffer1, buffer2, MPI_REPLACE)); 1981 if (PetscDefined(USE_DEBUG)) { 1982 for (p = 0; p < *chartSize; ++p) { 1983 PetscCheck(buffer2[p].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Found negative root rank %" PetscInt_FMT " at local point %" PetscInt_FMT " on rank %d when making migrationSF", buffer2[p].rank, p, rank); 1984 } 1985 } 1986 PetscCall(PetscFree2(buffer0, buffer1)); 1987 PetscCall(DMCreate(comm, distdm)); 1988 PetscCall(DMSetType(*distdm, DMPLEX)); 1989 PetscCall(PetscObjectSetName((PetscObject)*distdm, ((PetscObject)dm)->name)); 1990 PetscCall(DMPlexDistributionSetName(*distdm, distribution_name)); 1991 { 1992 PetscSF migrationSF; 1993 1994 PetscCall(PetscSFCreate(comm, &migrationSF)); 1995 PetscCall(PetscSFSetFromOptions(migrationSF)); 1996 PetscCall(PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER)); 1997 PetscCall(PetscSFSetUp(migrationSF)); 1998 PetscCall(DMPlexMigrate(dm, migrationSF, *distdm)); 1999 PetscCall(PetscSFDestroy(&migrationSF)); 2000 } 2001 } 2002 /* Set pointSF */ 2003 { 2004 PetscSF pointSF; 2005 PetscInt *ilocal, nleaves, q; 2006 PetscSFNode *iremote, *buffer0, *buffer1; 2007 2008 PetscCall(PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1)); 2009 for (p = 0, nleaves = 0; p < *chartSize; ++p) { 2010 if (owners[p] == rank) { 2011 buffer0[p].rank = rank; 2012 } else { 2013 buffer0[p].rank = -1; 2014 nleaves++; 2015 } 2016 buffer0[p].index = p; 2017 } 2018 for (p = 0; p < lsize; ++p) buffer1[p].rank = -1; 2019 PetscCall(PetscSFReduceBegin(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC)); 2020 PetscCall(PetscSFReduceEnd(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC)); 2021 for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1; 2022 PetscCall(PetscSFBcastBegin(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE)); 2023 PetscCall(PetscSFBcastEnd(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE)); 2024 if (PetscDefined(USE_DEBUG)) { 2025 for (p = 0; p < *chartSize; ++p) PetscCheck(buffer0[p].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Found negative root rank %" PetscInt_FMT " at local point %" PetscInt_FMT " on rank %d when making pointSF", buffer0[p].rank, p, rank); 2026 } 2027 PetscCall(PetscMalloc1(nleaves, &ilocal)); 2028 PetscCall(PetscMalloc1(nleaves, &iremote)); 2029 for (p = 0, q = 0; p < *chartSize; ++p) { 2030 if (buffer0[p].rank != rank) { 2031 ilocal[q] = p; 2032 iremote[q].rank = buffer0[p].rank; 2033 iremote[q].index = buffer0[p].index; 2034 q++; 2035 } 2036 } 2037 PetscCheck(q == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching leaf sizes: %" PetscInt_FMT " != %" PetscInt_FMT, q, nleaves); 2038 PetscCall(PetscFree2(buffer0, buffer1)); 2039 PetscCall(PetscSFCreate(comm, &pointSF)); 2040 PetscCall(PetscSFSetFromOptions(pointSF)); 2041 PetscCall(PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 2042 PetscCall(DMSetPointSF(*distdm, pointSF)); 2043 { 2044 DM cdm; 2045 2046 PetscCall(DMGetCoordinateDM(*distdm, &cdm)); 2047 PetscCall(DMSetPointSF(cdm, pointSF)); 2048 } 2049 PetscCall(PetscSFDestroy(&pointSF)); 2050 } 2051 PetscCall(ISRestoreIndices(chartSizesIS, &chartSize)); 2052 PetscCall(ISRestoreIndices(ownersIS, &owners)); 2053 PetscCall(ISRestoreIndices(gpointsIS, &gpoints)); 2054 PetscCall(ISDestroy(&chartSizesIS)); 2055 PetscCall(ISDestroy(&ownersIS)); 2056 PetscCall(ISDestroy(&gpointsIS)); 2057 /* Record that overlap has been manually created. */ 2058 /* This is to pass `DMPlexCheckPointSF()`, which checks that */ 2059 /* pointSF does not contain cells in the leaves if overlap = 0. */ 2060 PetscCall(DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL)); 2061 PetscCall(DMPlexDistributeSetDefault(*distdm, PETSC_FALSE)); 2062 PetscCall(DMPlexReorderSetDefault(*distdm, DM_REORDER_DEFAULT_FALSE)); 2063 PetscCall(PetscLogEventEnd(DMPLEX_DistributionLoad, viewer, 0, 0, 0)); 2064 PetscFunctionReturn(PETSC_SUCCESS); 2065 } 2066 2067 // Serial load of topology 2068 static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf) 2069 { 2070 MPI_Comm comm; 2071 const char *pointsName, *coneSizesName, *conesName, *orientationsName; 2072 IS pointsIS, coneSizesIS, conesIS, orientationsIS; 2073 const PetscInt *points, *coneSizes, *cones, *orientations; 2074 PetscInt *cone, *ornt; 2075 PetscInt dim, N, Np, pEnd, p, q, maxConeSize = 0, c; 2076 PetscMPIInt size, rank; 2077 2078 PetscFunctionBegin; 2079 pointsName = "order"; 2080 coneSizesName = "cones"; 2081 conesName = "cells"; 2082 orientationsName = "orientation"; 2083 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2084 PetscCallMPI(MPI_Comm_size(comm, &size)); 2085 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2086 PetscCall(ISCreate(comm, &pointsIS)); 2087 PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName)); 2088 PetscCall(ISCreate(comm, &coneSizesIS)); 2089 PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName)); 2090 PetscCall(ISCreate(comm, &conesIS)); 2091 PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName)); 2092 PetscCall(ISCreate(comm, &orientationsIS)); 2093 PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName)); 2094 PetscCall(PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim)); 2095 PetscCall(DMSetDimension(dm, dim)); 2096 { 2097 /* Force serial load */ 2098 PetscCall(PetscInfo(dm, "Loading DM %s in serial\n", dm->hdr.name)); 2099 PetscCall(PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np)); 2100 PetscCall(PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0)); 2101 PetscCall(PetscLayoutSetSize(pointsIS->map, Np)); 2102 pEnd = rank == 0 ? Np : 0; 2103 PetscCall(PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np)); 2104 PetscCall(PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0)); 2105 PetscCall(PetscLayoutSetSize(coneSizesIS->map, Np)); 2106 PetscCall(PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N)); 2107 PetscCall(PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0)); 2108 PetscCall(PetscLayoutSetSize(conesIS->map, N)); 2109 PetscCall(PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N)); 2110 PetscCall(PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0)); 2111 PetscCall(PetscLayoutSetSize(orientationsIS->map, N)); 2112 } 2113 PetscCall(ISLoad(pointsIS, viewer)); 2114 PetscCall(ISLoad(coneSizesIS, viewer)); 2115 PetscCall(ISLoad(conesIS, viewer)); 2116 PetscCall(ISLoad(orientationsIS, viewer)); 2117 /* Create Plex */ 2118 PetscCall(DMPlexSetChart(dm, 0, pEnd)); 2119 PetscCall(ISGetIndices(pointsIS, &points)); 2120 PetscCall(ISGetIndices(coneSizesIS, &coneSizes)); 2121 for (p = 0; p < pEnd; ++p) { 2122 PetscCall(DMPlexSetConeSize(dm, points[p], coneSizes[p])); 2123 maxConeSize = PetscMax(maxConeSize, coneSizes[p]); 2124 } 2125 PetscCall(DMSetUp(dm)); 2126 PetscCall(ISGetIndices(conesIS, &cones)); 2127 PetscCall(ISGetIndices(orientationsIS, &orientations)); 2128 PetscCall(PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt)); 2129 for (p = 0, q = 0; p < pEnd; ++p) { 2130 for (c = 0; c < coneSizes[p]; ++c, ++q) { 2131 cone[c] = cones[q]; 2132 ornt[c] = orientations[q]; 2133 } 2134 PetscCall(DMPlexSetCone(dm, points[p], cone)); 2135 PetscCall(DMPlexSetConeOrientation(dm, points[p], ornt)); 2136 } 2137 PetscCall(PetscFree2(cone, ornt)); 2138 /* Create global section migration SF */ 2139 if (sf) { 2140 PetscLayout layout; 2141 PetscInt *globalIndices; 2142 2143 PetscCall(PetscMalloc1(pEnd, &globalIndices)); 2144 /* plex point == globalPointNumber in this case */ 2145 for (p = 0; p < pEnd; ++p) globalIndices[p] = p; 2146 PetscCall(PetscLayoutCreate(comm, &layout)); 2147 PetscCall(PetscLayoutSetSize(layout, Np)); 2148 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 2149 PetscCall(PetscLayoutSetUp(layout)); 2150 PetscCall(PetscSFCreate(comm, sf)); 2151 PetscCall(PetscSFSetFromOptions(*sf)); 2152 PetscCall(PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices)); 2153 PetscCall(PetscLayoutDestroy(&layout)); 2154 PetscCall(PetscFree(globalIndices)); 2155 } 2156 /* Clean-up */ 2157 PetscCall(ISRestoreIndices(pointsIS, &points)); 2158 PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes)); 2159 PetscCall(ISRestoreIndices(conesIS, &cones)); 2160 PetscCall(ISRestoreIndices(orientationsIS, &orientations)); 2161 PetscCall(ISDestroy(&pointsIS)); 2162 PetscCall(ISDestroy(&coneSizesIS)); 2163 PetscCall(ISDestroy(&conesIS)); 2164 PetscCall(ISDestroy(&orientationsIS)); 2165 /* Fill in the rest of the topology structure */ 2166 PetscCall(DMPlexSymmetrize(dm)); 2167 PetscCall(DMPlexStratify(dm)); 2168 PetscFunctionReturn(PETSC_SUCCESS); 2169 } 2170 2171 /* Representation of two DMPlex strata in 0-based global numbering */ 2172 struct _n_PlexLayer { 2173 PetscInt d; 2174 IS conesIS, orientationsIS; 2175 PetscSection coneSizesSection; 2176 PetscLayout vertexLayout; 2177 PetscSF overlapSF, l2gSF; //TODO maybe confusing names (in DMPlex in general) 2178 PetscInt offset, conesOffset, leafOffset; 2179 }; 2180 typedef struct _n_PlexLayer *PlexLayer; 2181 2182 static PetscErrorCode PlexLayerDestroy(PlexLayer *layer) 2183 { 2184 PetscFunctionBegin; 2185 if (!*layer) PetscFunctionReturn(PETSC_SUCCESS); 2186 PetscCall(PetscSectionDestroy(&(*layer)->coneSizesSection)); 2187 PetscCall(ISDestroy(&(*layer)->conesIS)); 2188 PetscCall(ISDestroy(&(*layer)->orientationsIS)); 2189 PetscCall(PetscSFDestroy(&(*layer)->overlapSF)); 2190 PetscCall(PetscSFDestroy(&(*layer)->l2gSF)); 2191 PetscCall(PetscLayoutDestroy(&(*layer)->vertexLayout)); 2192 PetscCall(PetscFree(*layer)); 2193 PetscFunctionReturn(PETSC_SUCCESS); 2194 } 2195 2196 static PetscErrorCode PlexLayerCreate_Private(PlexLayer *layer) 2197 { 2198 PetscFunctionBegin; 2199 PetscCall(PetscNew(layer)); 2200 (*layer)->d = -1; 2201 (*layer)->offset = -1; 2202 (*layer)->conesOffset = -1; 2203 (*layer)->leafOffset = -1; 2204 PetscFunctionReturn(PETSC_SUCCESS); 2205 } 2206 2207 // Parallel load of a depth stratum 2208 static PetscErrorCode PlexLayerLoad_Private(PlexLayer layer, PetscViewer viewer, PetscInt d, PetscLayout pointsLayout) 2209 { 2210 char path[128]; 2211 MPI_Comm comm; 2212 const char *coneSizesName, *conesName, *orientationsName; 2213 IS coneSizesIS, conesIS, orientationsIS; 2214 PetscSection coneSizesSection; 2215 PetscLayout vertexLayout = NULL; 2216 PetscInt s; 2217 2218 PetscFunctionBegin; 2219 coneSizesName = "cone_sizes"; 2220 conesName = "cones"; 2221 orientationsName = "orientations"; 2222 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm)); 2223 2224 /* query size of next lower depth stratum (next lower dimension) */ 2225 if (d > 0) { 2226 PetscInt NVertices; 2227 2228 PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT "/%s", d - 1, coneSizesName)); 2229 PetscCall(PetscViewerHDF5ReadSizes(viewer, path, NULL, &NVertices)); 2230 PetscCall(PetscLayoutCreate(comm, &vertexLayout)); 2231 PetscCall(PetscLayoutSetSize(vertexLayout, NVertices)); 2232 PetscCall(PetscLayoutSetUp(vertexLayout)); 2233 } 2234 2235 PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT, d)); 2236 PetscCall(PetscViewerHDF5PushGroup(viewer, path)); 2237 2238 /* create coneSizesSection from stored IS coneSizes */ 2239 { 2240 const PetscInt *coneSizes; 2241 2242 PetscCall(ISCreate(comm, &coneSizesIS)); 2243 PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName)); 2244 if (pointsLayout) PetscCall(ISSetLayout(coneSizesIS, pointsLayout)); 2245 PetscCall(ISLoad(coneSizesIS, viewer)); 2246 if (!pointsLayout) PetscCall(ISGetLayout(coneSizesIS, &pointsLayout)); 2247 PetscCall(ISGetIndices(coneSizesIS, &coneSizes)); 2248 PetscCall(PetscSectionCreate(comm, &coneSizesSection)); 2249 //TODO different start ? 2250 PetscCall(PetscSectionSetChart(coneSizesSection, 0, pointsLayout->n)); 2251 for (s = 0; s < pointsLayout->n; ++s) PetscCall(PetscSectionSetDof(coneSizesSection, s, coneSizes[s])); 2252 PetscCall(PetscSectionSetUp(coneSizesSection)); 2253 PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes)); 2254 { 2255 PetscLayout tmp = NULL; 2256 2257 /* We need to keep the layout until the end of function */ 2258 PetscCall(PetscLayoutReference((PetscLayout)pointsLayout, &tmp)); 2259 } 2260 PetscCall(ISDestroy(&coneSizesIS)); 2261 } 2262 2263 /* use value layout of coneSizesSection as layout of cones and orientations */ 2264 { 2265 PetscLayout conesLayout; 2266 2267 PetscCall(PetscSectionGetValueLayout(comm, coneSizesSection, &conesLayout)); 2268 PetscCall(ISCreate(comm, &conesIS)); 2269 PetscCall(ISCreate(comm, &orientationsIS)); 2270 PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName)); 2271 PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName)); 2272 PetscCall(PetscLayoutDuplicate(conesLayout, &conesIS->map)); 2273 PetscCall(PetscLayoutDuplicate(conesLayout, &orientationsIS->map)); 2274 PetscCall(ISLoad(conesIS, viewer)); 2275 PetscCall(ISLoad(orientationsIS, viewer)); 2276 PetscCall(PetscLayoutDestroy(&conesLayout)); 2277 } 2278 2279 /* check assertion that layout of points is the same as point layout of coneSizesSection */ 2280 { 2281 PetscLayout pointsLayout0; 2282 PetscBool flg; 2283 2284 PetscCall(PetscSectionGetPointLayout(comm, coneSizesSection, &pointsLayout0)); 2285 PetscCall(PetscLayoutCompare(pointsLayout, pointsLayout0, &flg)); 2286 PetscCheck(flg, comm, PETSC_ERR_PLIB, "points layout != coneSizesSection point layout"); 2287 PetscCall(PetscLayoutDestroy(&pointsLayout0)); 2288 } 2289 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2290 PetscCall(PetscLayoutDestroy(&pointsLayout)); 2291 2292 layer->d = d; 2293 layer->conesIS = conesIS; 2294 layer->coneSizesSection = coneSizesSection; 2295 layer->orientationsIS = orientationsIS; 2296 layer->vertexLayout = vertexLayout; 2297 PetscFunctionReturn(PETSC_SUCCESS); 2298 } 2299 2300 static PetscErrorCode PlexLayerDistribute_Private(PlexLayer layer, PetscSF cellLocalToGlobalSF) 2301 { 2302 IS newConesIS, newOrientationsIS; 2303 PetscSection newConeSizesSection; 2304 MPI_Comm comm; 2305 2306 PetscFunctionBegin; 2307 PetscCall(PetscObjectGetComm((PetscObject)cellLocalToGlobalSF, &comm)); 2308 PetscCall(PetscSectionCreate(comm, &newConeSizesSection)); 2309 //TODO rename to something like ISDistribute() with optional PetscSection argument 2310 PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->conesIS, newConeSizesSection, &newConesIS)); 2311 PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->orientationsIS, newConeSizesSection, &newOrientationsIS)); 2312 2313 PetscCall(PetscObjectSetName((PetscObject)newConeSizesSection, ((PetscObject)layer->coneSizesSection)->name)); 2314 PetscCall(PetscObjectSetName((PetscObject)newConesIS, ((PetscObject)layer->conesIS)->name)); 2315 PetscCall(PetscObjectSetName((PetscObject)newOrientationsIS, ((PetscObject)layer->orientationsIS)->name)); 2316 PetscCall(PetscSectionDestroy(&layer->coneSizesSection)); 2317 PetscCall(ISDestroy(&layer->conesIS)); 2318 PetscCall(ISDestroy(&layer->orientationsIS)); 2319 layer->coneSizesSection = newConeSizesSection; 2320 layer->conesIS = newConesIS; 2321 layer->orientationsIS = newOrientationsIS; 2322 PetscFunctionReturn(PETSC_SUCCESS); 2323 } 2324 2325 //TODO share code with DMPlexBuildFromCellListParallel() 2326 #include <petsc/private/hashseti.h> 2327 static PetscErrorCode PlexLayerCreateSFs_Private(PlexLayer layer, PetscSF *vertexOverlapSF, PetscSF *sfXC) 2328 { 2329 PetscLayout vertexLayout = layer->vertexLayout; 2330 PetscSection coneSection = layer->coneSizesSection; 2331 IS cellVertexData = layer->conesIS; 2332 IS coneOrientations = layer->orientationsIS; 2333 PetscSF vl2gSF, vOverlapSF; 2334 PetscInt *verticesAdj; 2335 PetscInt i, n, numVerticesAdj; 2336 const PetscInt *cvd, *co = NULL; 2337 MPI_Comm comm; 2338 2339 PetscFunctionBegin; 2340 PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm)); 2341 PetscCall(PetscSectionGetStorageSize(coneSection, &n)); 2342 { 2343 PetscInt n0; 2344 2345 PetscCall(ISGetLocalSize(cellVertexData, &n0)); 2346 PetscCheck(n == n0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local size of IS cellVertexData = %" PetscInt_FMT " != %" PetscInt_FMT " = storage size of PetscSection coneSection", n0, n); 2347 PetscCall(ISGetIndices(cellVertexData, &cvd)); 2348 } 2349 if (coneOrientations) { 2350 PetscInt n0; 2351 2352 PetscCall(ISGetLocalSize(coneOrientations, &n0)); 2353 PetscCheck(n == n0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local size of IS coneOrientations = %" PetscInt_FMT " != %" PetscInt_FMT " = storage size of PetscSection coneSection", n0, n); 2354 PetscCall(ISGetIndices(coneOrientations, &co)); 2355 } 2356 /* Get/check global number of vertices */ 2357 { 2358 PetscInt NVerticesInCells = PETSC_INT_MIN; 2359 2360 /* NVerticesInCells = max(cellVertexData) + 1 */ 2361 for (i = 0; i < n; i++) 2362 if (cvd[i] > NVerticesInCells) NVerticesInCells = cvd[i]; 2363 ++NVerticesInCells; 2364 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, comm)); 2365 2366 if (vertexLayout->n == PETSC_DECIDE && vertexLayout->N == PETSC_DECIDE) vertexLayout->N = NVerticesInCells; 2367 else 2368 PetscCheck(vertexLayout->N == PETSC_DECIDE || vertexLayout->N >= NVerticesInCells, comm, PETSC_ERR_ARG_SIZ, "Specified global number of vertices %" PetscInt_FMT " must be greater than or equal to the number of unique vertices in the cell-vertex dataset %" PetscInt_FMT, 2369 vertexLayout->N, NVerticesInCells); 2370 PetscCall(PetscLayoutSetUp(vertexLayout)); 2371 } 2372 /* Find locally unique vertices in cellVertexData */ 2373 { 2374 PetscHSetI vhash; 2375 PetscInt off = 0; 2376 2377 PetscCall(PetscHSetICreate(&vhash)); 2378 for (i = 0; i < n; ++i) PetscCall(PetscHSetIAdd(vhash, cvd[i])); 2379 PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj)); 2380 PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj)); 2381 PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj)); 2382 PetscAssert(off == numVerticesAdj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assertion failed: off == numVerticesAdj (%" PetscInt_FMT " != %" PetscInt_FMT ")", off, numVerticesAdj); 2383 PetscCall(PetscHSetIDestroy(&vhash)); 2384 } 2385 /* We must sort vertices to preserve numbering */ 2386 PetscCall(PetscSortInt(numVerticesAdj, verticesAdj)); 2387 /* Connect locally unique vertices with star forest */ 2388 PetscCall(PetscSFCreateByMatchingIndices(vertexLayout, numVerticesAdj, verticesAdj, NULL, 0, numVerticesAdj, verticesAdj, NULL, 0, &vl2gSF, &vOverlapSF)); 2389 PetscCall(PetscObjectSetName((PetscObject)vOverlapSF, "overlapSF")); 2390 PetscCall(PetscObjectSetName((PetscObject)vl2gSF, "localToGlobalSF")); 2391 2392 PetscCall(PetscFree(verticesAdj)); 2393 *vertexOverlapSF = vOverlapSF; 2394 *sfXC = vl2gSF; 2395 PetscFunctionReturn(PETSC_SUCCESS); 2396 } 2397 2398 static PetscErrorCode PlexLayerCreateCellSFs_Private(PlexLayer layer, PetscSF *cellOverlapSF, PetscSF *cellLocalToGlobalSF) 2399 { 2400 PetscSection coneSection = layer->coneSizesSection; 2401 PetscInt nCells; 2402 MPI_Comm comm; 2403 2404 PetscFunctionBegin; 2405 PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm)); 2406 { 2407 PetscInt cStart; 2408 2409 PetscCall(PetscSectionGetChart(coneSection, &cStart, &nCells)); 2410 PetscCheck(cStart == 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "coneSection must start at 0"); 2411 } 2412 /* Create overlapSF as empty SF with the right number of roots */ 2413 PetscCall(PetscSFCreate(comm, cellOverlapSF)); 2414 PetscCall(PetscSFSetGraph(*cellOverlapSF, nCells, 0, NULL, PETSC_USE_POINTER, NULL, PETSC_USE_POINTER)); 2415 PetscCall(PetscSFSetUp(*cellOverlapSF)); 2416 /* Create localToGlobalSF as identity mapping */ 2417 { 2418 PetscLayout map; 2419 2420 PetscCall(PetscLayoutCreateFromSizes(comm, nCells, PETSC_DECIDE, 1, &map)); 2421 PetscCall(PetscSFCreateFromLayouts(map, map, cellLocalToGlobalSF)); 2422 PetscCall(PetscSFSetUp(*cellLocalToGlobalSF)); 2423 PetscCall(PetscLayoutDestroy(&map)); 2424 } 2425 PetscFunctionReturn(PETSC_SUCCESS); 2426 } 2427 2428 static PetscErrorCode DMPlexTopologyBuildFromLayers_Private(DM dm, PetscInt depth, PlexLayer *layers, IS strataPermutation) 2429 { 2430 const PetscInt *permArr; 2431 PetscInt d, nPoints; 2432 MPI_Comm comm; 2433 2434 PetscFunctionBegin; 2435 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2436 PetscCall(ISGetIndices(strataPermutation, &permArr)); 2437 2438 /* Count points, strata offsets and cones offsets (taking strataPermutation into account) */ 2439 { 2440 PetscInt stratumOffset = 0; 2441 PetscInt conesOffset = 0; 2442 2443 for (d = 0; d <= depth; d++) { 2444 const PetscInt e = permArr[d]; 2445 const PlexLayer l = layers[e]; 2446 PetscInt lo, n, size; 2447 2448 PetscCall(PetscSectionGetChart(l->coneSizesSection, &lo, &n)); 2449 PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &size)); 2450 PetscCheck(lo == 0, comm, PETSC_ERR_PLIB, "starting point should be 0 in coneSizesSection %" PetscInt_FMT, d); 2451 l->offset = stratumOffset; 2452 l->conesOffset = conesOffset; 2453 stratumOffset += n; 2454 conesOffset += size; 2455 } 2456 nPoints = stratumOffset; 2457 } 2458 2459 /* Set interval for all plex points */ 2460 //TODO we should store starting point of plex 2461 PetscCall(DMPlexSetChart(dm, 0, nPoints)); 2462 2463 /* Set up plex coneSection from layer coneSections */ 2464 { 2465 PetscSection coneSection; 2466 2467 PetscCall(DMPlexGetConeSection(dm, &coneSection)); 2468 for (d = 0; d <= depth; d++) { 2469 const PlexLayer l = layers[d]; 2470 PetscInt n, q; 2471 2472 PetscCall(PetscSectionGetChart(l->coneSizesSection, NULL, &n)); 2473 for (q = 0; q < n; q++) { 2474 const PetscInt p = l->offset + q; 2475 PetscInt coneSize; 2476 2477 PetscCall(PetscSectionGetDof(l->coneSizesSection, q, &coneSize)); 2478 PetscCall(PetscSectionSetDof(coneSection, p, coneSize)); 2479 } 2480 } 2481 } 2482 //TODO this is terrible, DMSetUp_Plex() should be DMPlexSetUpSections() or so 2483 PetscCall(DMSetUp(dm)); 2484 2485 /* Renumber cones points from layer-global numbering to plex-local numbering */ 2486 { 2487 PetscInt *cones, *ornts; 2488 2489 PetscCall(DMPlexGetCones(dm, &cones)); 2490 PetscCall(DMPlexGetConeOrientations(dm, &ornts)); 2491 for (d = 1; d <= depth; d++) { 2492 const PlexLayer l = layers[d]; 2493 PetscInt i, lConesSize; 2494 PetscInt *lCones; 2495 const PetscInt *lOrnts; 2496 PetscInt *pCones = &cones[l->conesOffset]; 2497 PetscInt *pOrnts = &ornts[l->conesOffset]; 2498 2499 PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &lConesSize)); 2500 /* Get cones in local plex numbering */ 2501 { 2502 ISLocalToGlobalMapping l2g; 2503 PetscLayout vertexLayout = l->vertexLayout; 2504 PetscSF vertexSF = layers[d - 1]->l2gSF; /* vertices of this layer are cells of previous layer */ 2505 const PetscInt *gCones; 2506 PetscInt lConesSize0; 2507 2508 PetscCall(ISGetLocalSize(l->conesIS, &lConesSize0)); 2509 PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(conesIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize); 2510 PetscCall(ISGetLocalSize(l->orientationsIS, &lConesSize0)); 2511 PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(orientationsIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize); 2512 2513 PetscCall(PetscMalloc1(lConesSize, &lCones)); 2514 PetscCall(ISGetIndices(l->conesIS, &gCones)); 2515 PetscCall(ISLocalToGlobalMappingCreateSF(vertexSF, vertexLayout->rstart, &l2g)); 2516 PetscCall(ISGlobalToLocalMappingApply(l2g, IS_GTOLM_MASK, lConesSize, gCones, &lConesSize0, lCones)); 2517 PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "global to local does not cover all indices (%" PetscInt_FMT " of %" PetscInt_FMT ")", lConesSize0, lConesSize); 2518 PetscCall(ISLocalToGlobalMappingDestroy(&l2g)); 2519 PetscCall(ISRestoreIndices(l->conesIS, &gCones)); 2520 } 2521 PetscCall(ISGetIndices(l->orientationsIS, &lOrnts)); 2522 /* Set cones, need to add stratum offset */ 2523 for (i = 0; i < lConesSize; i++) { 2524 pCones[i] = lCones[i] + layers[d - 1]->offset; /* cone points of current layer are points of previous layer */ 2525 pOrnts[i] = lOrnts[i]; 2526 } 2527 PetscCall(PetscFree(lCones)); 2528 PetscCall(ISRestoreIndices(l->orientationsIS, &lOrnts)); 2529 } 2530 } 2531 PetscCall(DMPlexSymmetrize(dm)); 2532 PetscCall(DMPlexStratify(dm)); 2533 PetscCall(ISRestoreIndices(strataPermutation, &permArr)); 2534 PetscFunctionReturn(PETSC_SUCCESS); 2535 } 2536 2537 static PetscErrorCode PlexLayerConcatenateSFs_Private(MPI_Comm comm, PetscInt depth, PlexLayer layers[], IS strataPermutation, PetscSF *overlapSF, PetscSF *l2gSF) 2538 { 2539 PetscInt d; 2540 PetscSF *osfs, *lsfs; 2541 PetscInt *leafOffsets; 2542 const PetscInt *permArr; 2543 2544 PetscFunctionBegin; 2545 PetscCall(ISGetIndices(strataPermutation, &permArr)); 2546 PetscCall(PetscCalloc3(depth + 1, &osfs, depth + 1, &lsfs, depth + 1, &leafOffsets)); 2547 for (d = 0; d <= depth; d++) { 2548 const PetscInt e = permArr[d]; 2549 2550 PetscAssert(e == layers[e]->d, PETSC_COMM_SELF, PETSC_ERR_PLIB, "assertion: e == layers[e]->d"); 2551 osfs[d] = layers[e]->overlapSF; 2552 lsfs[d] = layers[e]->l2gSF; 2553 leafOffsets[d] = layers[e]->offset; 2554 } 2555 PetscCall(PetscSFConcatenate(comm, depth + 1, osfs, PETSCSF_CONCATENATE_ROOTMODE_LOCAL, leafOffsets, overlapSF)); 2556 PetscCall(PetscSFConcatenate(comm, depth + 1, lsfs, PETSCSF_CONCATENATE_ROOTMODE_GLOBAL, leafOffsets, l2gSF)); 2557 PetscCall(PetscFree3(osfs, lsfs, leafOffsets)); 2558 PetscCall(ISRestoreIndices(strataPermutation, &permArr)); 2559 PetscFunctionReturn(PETSC_SUCCESS); 2560 } 2561 2562 // Parallel load of topology 2563 static PetscErrorCode DMPlexTopologyLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF *sfXC) 2564 { 2565 PlexLayer *layers; 2566 IS strataPermutation; 2567 PetscLayout pointsLayout; 2568 PetscInt depth; 2569 PetscInt d; 2570 MPI_Comm comm; 2571 2572 PetscFunctionBegin; 2573 { 2574 PetscInt dim; 2575 2576 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "depth", PETSC_INT, NULL, &depth)); 2577 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "cell_dim", PETSC_INT, NULL, &dim)); 2578 PetscCall(DMSetDimension(dm, dim)); 2579 } 2580 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2581 2582 PetscCall(PetscInfo(dm, "Loading DM %s in parallel\n", dm->hdr.name)); 2583 { 2584 IS spOnComm; 2585 2586 PetscCall(ISCreate(comm, &spOnComm)); 2587 PetscCall(PetscObjectSetName((PetscObject)spOnComm, "permutation")); 2588 PetscCall(ISLoad(spOnComm, viewer)); 2589 /* have the same serial IS on every rank */ 2590 PetscCall(ISAllGather(spOnComm, &strataPermutation)); 2591 PetscCall(PetscObjectSetName((PetscObject)strataPermutation, ((PetscObject)spOnComm)->name)); 2592 PetscCall(ISDestroy(&spOnComm)); 2593 } 2594 2595 /* Create layers, load raw data for each layer */ 2596 PetscCall(PetscViewerHDF5PushGroup(viewer, "strata")); 2597 PetscCall(PetscMalloc1(depth + 1, &layers)); 2598 for (d = depth, pointsLayout = NULL; d >= 0; pointsLayout = layers[d]->vertexLayout, d--) { 2599 PetscCall(PlexLayerCreate_Private(&layers[d])); 2600 PetscCall(PlexLayerLoad_Private(layers[d], viewer, d, pointsLayout)); 2601 } 2602 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */ 2603 2604 for (d = depth; d >= 0; d--) { 2605 /* Redistribute cells and vertices for each applicable layer */ 2606 if (d < depth) PetscCall(PlexLayerDistribute_Private(layers[d], layers[d]->l2gSF)); 2607 /* Create vertex overlap SF and vertex localToGlobal SF for each applicable layer */ 2608 if (d > 0) PetscCall(PlexLayerCreateSFs_Private(layers[d], &layers[d - 1]->overlapSF, &layers[d - 1]->l2gSF)); 2609 } 2610 /* Build trivial SFs for the cell layer as well */ 2611 PetscCall(PlexLayerCreateCellSFs_Private(layers[depth], &layers[depth]->overlapSF, &layers[depth]->l2gSF)); 2612 2613 /* Build DMPlex topology from the layers */ 2614 PetscCall(DMPlexTopologyBuildFromLayers_Private(dm, depth, layers, strataPermutation)); 2615 2616 /* Build overall point SF alias overlap SF */ 2617 { 2618 PetscSF overlapSF; 2619 2620 PetscCall(PlexLayerConcatenateSFs_Private(comm, depth, layers, strataPermutation, &overlapSF, sfXC)); 2621 PetscCall(DMSetPointSF(dm, overlapSF)); 2622 PetscCall(PetscSFDestroy(&overlapSF)); 2623 } 2624 2625 for (d = depth; d >= 0; d--) PetscCall(PlexLayerDestroy(&layers[d])); 2626 PetscCall(PetscFree(layers)); 2627 PetscCall(ISDestroy(&strataPermutation)); 2628 PetscFunctionReturn(PETSC_SUCCESS); 2629 } 2630 2631 PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC) 2632 { 2633 DMPlexStorageVersion version; 2634 const char *topologydm_name; 2635 char group[PETSC_MAX_PATH_LEN]; 2636 PetscSF sfwork = NULL; 2637 2638 PetscFunctionBegin; 2639 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 2640 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version)); 2641 if (DMPlexStorageVersionGE(version, 2, 0, 0)) { 2642 PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name)); 2643 } else { 2644 PetscCall(PetscStrncpy(group, "/", sizeof(group))); 2645 } 2646 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 2647 2648 PetscCall(PetscViewerHDF5PushGroup(viewer, "topology")); 2649 if (version->major < 3) { 2650 PetscCall(DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork)); 2651 } else { 2652 /* since DMPlexStorageVersion 3.0.0 */ 2653 PetscCall(DMPlexTopologyLoad_HDF5_Private(dm, viewer, &sfwork)); 2654 } 2655 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */ 2656 2657 if (DMPlexStorageVersionGE(version, 2, 0, 0)) { 2658 DM distdm; 2659 PetscSF distsf; 2660 const char *distribution_name; 2661 PetscBool exists; 2662 2663 PetscCall(DMPlexDistributionGetName(dm, &distribution_name)); 2664 /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */ 2665 PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions")); 2666 PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &exists)); 2667 if (exists) { 2668 PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name)); 2669 PetscCall(DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm)); 2670 if (distdm) { 2671 PetscCall(DMPlexReplace_Internal(dm, &distdm)); 2672 PetscCall(PetscSFDestroy(&sfwork)); 2673 sfwork = distsf; 2674 } 2675 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */ 2676 } 2677 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */ 2678 } 2679 if (sfXC) { 2680 *sfXC = sfwork; 2681 } else { 2682 PetscCall(PetscSFDestroy(&sfwork)); 2683 } 2684 2685 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2686 PetscFunctionReturn(PETSC_SUCCESS); 2687 } 2688 2689 /* If the file is old, it not only has different path to the coordinates, but */ 2690 /* does not contain coordinateDMs, so must fall back to the old implementation. */ 2691 static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer) 2692 { 2693 PetscSection coordSection; 2694 Vec coordinates; 2695 PetscReal lengthScale; 2696 PetscInt spatialDim, N, numVertices, vStart, vEnd, v; 2697 PetscMPIInt rank; 2698 2699 PetscFunctionBegin; 2700 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 2701 /* Read geometry */ 2702 PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry")); 2703 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates)); 2704 PetscCall(PetscObjectSetName((PetscObject)coordinates, "vertices")); 2705 { 2706 /* Force serial load */ 2707 PetscCall(PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N)); 2708 PetscCall(VecSetSizes(coordinates, rank == 0 ? N : 0, N)); 2709 PetscCall(VecSetBlockSize(coordinates, spatialDim)); 2710 } 2711 PetscCall(VecLoad(coordinates, viewer)); 2712 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2713 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 2714 PetscCall(VecScale(coordinates, 1.0 / lengthScale)); 2715 PetscCall(VecGetLocalSize(coordinates, &numVertices)); 2716 PetscCall(VecGetBlockSize(coordinates, &spatialDim)); 2717 numVertices /= spatialDim; 2718 /* Create coordinates */ 2719 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 2720 PetscCheck(numVertices == vEnd - vStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of coordinates loaded %" PetscInt_FMT " does not match number of vertices %" PetscInt_FMT, numVertices, vEnd - vStart); 2721 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2722 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 2723 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spatialDim)); 2724 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 2725 for (v = vStart; v < vEnd; ++v) { 2726 PetscCall(PetscSectionSetDof(coordSection, v, spatialDim)); 2727 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spatialDim)); 2728 } 2729 PetscCall(PetscSectionSetUp(coordSection)); 2730 PetscCall(DMSetCoordinates(dm, coordinates)); 2731 PetscCall(VecDestroy(&coordinates)); 2732 PetscFunctionReturn(PETSC_SUCCESS); 2733 } 2734 2735 PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC) 2736 { 2737 DMPlexStorageVersion version; 2738 DM cdm; 2739 Vec coords; 2740 PetscInt blockSize; 2741 PetscReal lengthScale; 2742 PetscSF lsf; 2743 const char *topologydm_name; 2744 char *coordinatedm_name, *coordinates_name; 2745 2746 PetscFunctionBegin; 2747 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version)); 2748 if (!DMPlexStorageVersionGE(version, 2, 0, 0)) { 2749 PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer)); 2750 PetscFunctionReturn(PETSC_SUCCESS); 2751 } 2752 /* else: since DMPlexStorageVersion 2.0.0 */ 2753 PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load"); 2754 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 2755 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 2756 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 2757 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name)); 2758 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name)); 2759 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2760 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2761 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2762 PetscCall(PetscObjectSetName((PetscObject)cdm, coordinatedm_name)); 2763 PetscCall(PetscFree(coordinatedm_name)); 2764 /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */ 2765 PetscCall(DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf)); 2766 PetscCall(DMCreateLocalVector(cdm, &coords)); 2767 PetscCall(PetscObjectSetName((PetscObject)coords, coordinates_name)); 2768 PetscCall(PetscFree(coordinates_name)); 2769 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE)); 2770 PetscCall(DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords)); 2771 PetscCall(PetscViewerPopFormat(viewer)); 2772 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 2773 PetscCall(VecScale(coords, 1.0 / lengthScale)); 2774 PetscCall(DMSetCoordinatesLocal(dm, coords)); 2775 PetscCall(VecGetBlockSize(coords, &blockSize)); 2776 PetscCall(DMSetCoordinateDim(dm, blockSize)); 2777 PetscCall(VecDestroy(&coords)); 2778 PetscCall(PetscSFDestroy(&lsf)); 2779 PetscFunctionReturn(PETSC_SUCCESS); 2780 } 2781 2782 PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer) 2783 { 2784 DMPlexStorageVersion version; 2785 2786 PetscFunctionBegin; 2787 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version)); 2788 PetscCall(PetscInfo(dm, "Loading DM %s storage version %d.%d.%d\n", dm->hdr.name, version->major, version->minor, version->subminor)); 2789 if (!DMPlexStorageVersionGE(version, 2, 0, 0)) { 2790 PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL)); 2791 PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL)); 2792 PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer)); 2793 } else { 2794 PetscSF sfXC; 2795 2796 /* since DMPlexStorageVersion 2.0.0 */ 2797 PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC)); 2798 PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC)); 2799 PetscCall(DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC)); 2800 PetscCall(PetscSFDestroy(&sfXC)); 2801 } 2802 PetscFunctionReturn(PETSC_SUCCESS); 2803 } 2804 2805 static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF) 2806 { 2807 MPI_Comm comm; 2808 PetscInt pStart, pEnd, p, m; 2809 PetscInt *goffs, *ilocal; 2810 PetscBool rootIncludeConstraints, leafIncludeConstraints; 2811 2812 PetscFunctionBegin; 2813 PetscCall(PetscObjectGetComm((PetscObject)leafSection, &comm)); 2814 PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd)); 2815 PetscCall(PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints)); 2816 PetscCall(PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints)); 2817 if (rootIncludeConstraints && leafIncludeConstraints) PetscCall(PetscSectionGetStorageSize(leafSection, &m)); 2818 else PetscCall(PetscSectionGetConstrainedStorageSize(leafSection, &m)); 2819 PetscCall(PetscMalloc1(m, &ilocal)); 2820 PetscCall(PetscMalloc1(m, &goffs)); 2821 /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */ 2822 /* for the top-level section (not for each field), so one must have */ 2823 /* rootSection->pointMajor == PETSC_TRUE. */ 2824 PetscCheck(rootSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering"); 2825 /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */ 2826 PetscCheck(leafSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering"); 2827 for (p = pStart, m = 0; p < pEnd; ++p) { 2828 PetscInt dof, cdof, i, j, off, goff; 2829 const PetscInt *cinds; 2830 2831 PetscCall(PetscSectionGetDof(leafSection, p, &dof)); 2832 if (dof < 0) continue; 2833 goff = globalOffsets[p - pStart]; 2834 PetscCall(PetscSectionGetOffset(leafSection, p, &off)); 2835 PetscCall(PetscSectionGetConstraintDof(leafSection, p, &cdof)); 2836 PetscCall(PetscSectionGetConstraintIndices(leafSection, p, &cinds)); 2837 for (i = 0, j = 0; i < dof; ++i) { 2838 PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]); 2839 2840 if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) { 2841 ilocal[m] = off++; 2842 goffs[m++] = goff++; 2843 } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off; 2844 else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff; 2845 if (constrained) ++j; 2846 } 2847 } 2848 PetscCall(PetscSFCreate(comm, sectionSF)); 2849 PetscCall(PetscSFSetFromOptions(*sectionSF)); 2850 PetscCall(PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs)); 2851 PetscCall(PetscFree(goffs)); 2852 PetscFunctionReturn(PETSC_SUCCESS); 2853 } 2854 2855 PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf) 2856 { 2857 MPI_Comm comm; 2858 PetscMPIInt size, rank; 2859 const char *topologydm_name; 2860 const char *sectiondm_name; 2861 PetscSection sectionA, sectionB; 2862 PetscBool has; 2863 PetscInt nX, n, i; 2864 PetscSF sfAB; 2865 2866 PetscFunctionBegin; 2867 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2868 PetscCallMPI(MPI_Comm_size(comm, &size)); 2869 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2870 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 2871 PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name)); 2872 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 2873 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 2874 PetscCall(PetscViewerHDF5PushGroup(viewer, "dms")); 2875 PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name)); 2876 /* A: on-disk points */ 2877 /* X: list of global point numbers, [0, NX) */ 2878 /* B: plex points */ 2879 /* Load raw section (sectionA) */ 2880 PetscCall(PetscSectionCreate(comm, §ionA)); 2881 PetscCall(PetscViewerHDF5HasGroup(viewer, "section", &has)); 2882 if (has) PetscCall(PetscSectionLoad(sectionA, viewer)); 2883 else { 2884 // TODO If section is missing, create the default affine section with dim dofs on each vertex. Use PetscSplitOwnership() to split vertices. 2885 // How do I know the total number of vertices? 2886 PetscInt dim, Nf = 1, Nv, nv = PETSC_DECIDE; 2887 2888 PetscCall(DMGetDimension(dm, &dim)); 2889 PetscCall(DMPlexGetDepthStratumGlobalSize(dm, 0, &Nv)); 2890 PetscCall(PetscSectionSetNumFields(sectionA, Nf)); 2891 PetscCall(PetscSectionSetFieldName(sectionA, 0, "Cartesian")); 2892 PetscCall(PetscSectionSetFieldComponents(sectionA, 0, dim)); 2893 for (PetscInt c = 0; c < dim; ++c) { 2894 char axis = 'X' + (char)c; 2895 2896 PetscCall(PetscSectionSetComponentName(sectionA, 0, c, &axis)); 2897 } 2898 PetscCall(PetscSplitOwnership(comm, &nv, &Nv)); 2899 PetscCall(PetscSectionSetChart(sectionA, 0, nv)); 2900 for (PetscInt p = 0; p < nv; ++p) { 2901 PetscCall(PetscSectionSetDof(sectionA, p, dim)); 2902 PetscCall(PetscSectionSetFieldDof(sectionA, p, 0, dim)); 2903 } 2904 PetscCall(PetscSectionSetUp(sectionA)); 2905 } 2906 PetscCall(PetscSectionGetChart(sectionA, NULL, &n)); 2907 /* Create sfAB: A -> B */ 2908 #if defined(PETSC_USE_DEBUG) 2909 { 2910 PetscInt N, N1; 2911 2912 PetscCall(PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1)); 2913 PetscCallMPI(MPIU_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm)); 2914 PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Mismatching sizes: on-disk order array size (%" PetscInt_FMT ") != number of loaded section points (%" PetscInt_FMT ")", N1, N); 2915 } 2916 #endif 2917 { 2918 IS orderIS; 2919 const PetscInt *gpoints; 2920 PetscSF sfXA, sfAX; 2921 PetscLayout layout; 2922 PetscSFNode *owners, *buffer; 2923 PetscInt nleaves; 2924 PetscInt *ilocal; 2925 PetscSFNode *iremote; 2926 2927 /* Create sfAX: A -> X */ 2928 PetscCall(ISCreate(comm, &orderIS)); 2929 PetscCall(PetscObjectSetName((PetscObject)orderIS, "order")); 2930 PetscCall(PetscLayoutSetLocalSize(orderIS->map, n)); 2931 PetscCall(ISLoad(orderIS, viewer)); 2932 PetscCall(PetscLayoutCreate(comm, &layout)); 2933 PetscCall(PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL)); 2934 PetscCall(PetscLayoutSetLocalSize(layout, nX)); 2935 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 2936 PetscCall(PetscLayoutSetUp(layout)); 2937 PetscCall(PetscSFCreate(comm, &sfXA)); 2938 PetscCall(ISGetIndices(orderIS, &gpoints)); 2939 PetscCall(PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints)); 2940 PetscCall(ISRestoreIndices(orderIS, &gpoints)); 2941 PetscCall(ISDestroy(&orderIS)); 2942 PetscCall(PetscLayoutDestroy(&layout)); 2943 PetscCall(PetscMalloc1(n, &owners)); 2944 PetscCall(PetscMalloc1(nX, &buffer)); 2945 for (i = 0; i < n; ++i) { 2946 owners[i].rank = rank; 2947 owners[i].index = i; 2948 } 2949 for (i = 0; i < nX; ++i) { 2950 buffer[i].rank = -1; 2951 buffer[i].index = -1; 2952 } 2953 PetscCall(PetscSFReduceBegin(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC)); 2954 PetscCall(PetscSFReduceEnd(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC)); 2955 PetscCall(PetscSFDestroy(&sfXA)); 2956 PetscCall(PetscFree(owners)); 2957 for (i = 0, nleaves = 0; i < nX; ++i) 2958 if (buffer[i].rank >= 0) nleaves++; 2959 PetscCall(PetscMalloc1(nleaves, &ilocal)); 2960 PetscCall(PetscMalloc1(nleaves, &iremote)); 2961 for (i = 0, nleaves = 0; i < nX; ++i) { 2962 if (buffer[i].rank >= 0) { 2963 ilocal[nleaves] = i; 2964 iremote[nleaves].rank = buffer[i].rank; 2965 iremote[nleaves].index = buffer[i].index; 2966 nleaves++; 2967 } 2968 } 2969 PetscCall(PetscFree(buffer)); 2970 PetscCall(PetscSFCreate(comm, &sfAX)); 2971 PetscCall(PetscSFSetFromOptions(sfAX)); 2972 PetscCall(PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 2973 PetscCall(PetscSFCompose(sfAX, sfXB, &sfAB)); 2974 PetscCall(PetscSFDestroy(&sfAX)); 2975 } 2976 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2977 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2978 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2979 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2980 /* Create plex section (sectionB) */ 2981 PetscCall(DMGetLocalSection(sectiondm, §ionB)); 2982 if (lsf || gsf) { 2983 PetscLayout layout; 2984 PetscInt M, m; 2985 PetscInt *offsetsA; 2986 PetscBool includesConstraintsA; 2987 2988 PetscCall(PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB)); 2989 PetscCall(PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA)); 2990 if (includesConstraintsA) PetscCall(PetscSectionGetStorageSize(sectionA, &m)); 2991 else PetscCall(PetscSectionGetConstrainedStorageSize(sectionA, &m)); 2992 PetscCallMPI(MPIU_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm)); 2993 PetscCall(PetscLayoutCreate(comm, &layout)); 2994 PetscCall(PetscLayoutSetSize(layout, M)); 2995 PetscCall(PetscLayoutSetUp(layout)); 2996 if (lsf) { 2997 PetscSF lsfABdata; 2998 2999 PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata)); 3000 *lsf = lsfABdata; 3001 } 3002 if (gsf) { 3003 PetscSection gsectionB, gsectionB1; 3004 PetscBool includesConstraintsB; 3005 PetscSF gsfABdata, pointsf; 3006 3007 PetscCall(DMGetGlobalSection(sectiondm, &gsectionB1)); 3008 PetscCall(PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB)); 3009 PetscCall(DMGetPointSF(sectiondm, &pointsf)); 3010 PetscCall(PetscSectionCreateGlobalSection(sectionB, pointsf, PETSC_TRUE, includesConstraintsB, PETSC_TRUE, &gsectionB)); 3011 PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata)); 3012 PetscCall(PetscSectionDestroy(&gsectionB)); 3013 *gsf = gsfABdata; 3014 } 3015 PetscCall(PetscLayoutDestroy(&layout)); 3016 PetscCall(PetscFree(offsetsA)); 3017 } else { 3018 PetscCall(PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB)); 3019 } 3020 PetscCall(PetscSFDestroy(&sfAB)); 3021 PetscCall(PetscSectionDestroy(§ionA)); 3022 PetscFunctionReturn(PETSC_SUCCESS); 3023 } 3024 3025 PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec) 3026 { 3027 MPI_Comm comm; 3028 const char *topologydm_name; 3029 const char *sectiondm_name; 3030 const char *vec_name; 3031 Vec vecA; 3032 PetscInt mA, m, bs; 3033 const PetscInt *ilocal; 3034 const PetscScalar *src; 3035 PetscScalar *dest; 3036 3037 PetscFunctionBegin; 3038 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3039 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 3040 PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name)); 3041 PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name)); 3042 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 3043 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 3044 PetscCall(PetscViewerHDF5PushGroup(viewer, "dms")); 3045 PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name)); 3046 PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs")); 3047 PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name)); 3048 PetscCall(VecCreate(comm, &vecA)); 3049 PetscCall(PetscObjectSetName((PetscObject)vecA, vec_name)); 3050 PetscCall(PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL)); 3051 /* Check consistency */ 3052 { 3053 PetscSF pointsf, pointsf1; 3054 PetscInt m1, i, j; 3055 3056 PetscCall(DMGetPointSF(dm, &pointsf)); 3057 PetscCall(DMGetPointSF(sectiondm, &pointsf1)); 3058 PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm"); 3059 #if defined(PETSC_USE_DEBUG) 3060 { 3061 PetscInt MA, MA1; 3062 3063 PetscCallMPI(MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm)); 3064 PetscCall(PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1)); 3065 PetscCheck(MA1 == MA, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%" PetscInt_FMT ") != On-disk vector data size (%" PetscInt_FMT ")", MA, MA1); 3066 } 3067 #endif 3068 PetscCall(VecGetLocalSize(vec, &m1)); 3069 PetscCheck(m1 >= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%" PetscInt_FMT ") < SF leaf size (%" PetscInt_FMT ")", m1, m); 3070 for (i = 0; i < m; ++i) { 3071 j = ilocal ? ilocal[i] : i; 3072 PetscCheck(j >= 0 && j < m1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Leaf's %" PetscInt_FMT "-th index, %" PetscInt_FMT ", not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", i, j, (PetscInt)0, m1); 3073 } 3074 } 3075 PetscCall(VecSetSizes(vecA, mA, PETSC_DECIDE)); 3076 PetscCall(VecLoad(vecA, viewer)); 3077 PetscCall(VecGetArrayRead(vecA, &src)); 3078 PetscCall(VecGetArray(vec, &dest)); 3079 PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE)); 3080 PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE)); 3081 PetscCall(VecRestoreArray(vec, &dest)); 3082 PetscCall(VecRestoreArrayRead(vecA, &src)); 3083 PetscCall(VecDestroy(&vecA)); 3084 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs)); 3085 PetscCall(VecSetBlockSize(vec, bs)); 3086 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3087 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3088 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3089 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3090 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3091 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3092 PetscFunctionReturn(PETSC_SUCCESS); 3093 } 3094