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