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