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, MPI_C_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, MPI_C_BOOL, A_mask, X_mask, MPI_REPLACE)); 1827 PetscCall(PetscSFReduceEnd(sfXA, MPI_C_BOOL, A_mask, X_mask, MPI_REPLACE)); 1828 PetscCall(PetscSFBcastBegin(sfXC, MPI_C_BOOL, X_mask, C_mask, MPI_LOR)); 1829 PetscCall(PetscSFBcastEnd(sfXC, MPI_C_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) 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); 2026 } 2027 PetscCall(PetscFree2(buffer0, buffer1)); 2028 PetscCall(DMCreate(comm, distdm)); 2029 PetscCall(DMSetType(*distdm, DMPLEX)); 2030 PetscCall(PetscObjectSetName((PetscObject)*distdm, ((PetscObject)dm)->name)); 2031 PetscCall(DMPlexDistributionSetName(*distdm, distribution_name)); 2032 { 2033 PetscSF migrationSF; 2034 2035 PetscCall(PetscSFCreate(comm, &migrationSF)); 2036 PetscCall(PetscSFSetFromOptions(migrationSF)); 2037 PetscCall(PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER)); 2038 PetscCall(PetscSFSetUp(migrationSF)); 2039 PetscCall(DMPlexMigrate(dm, migrationSF, *distdm)); 2040 PetscCall(PetscSFDestroy(&migrationSF)); 2041 } 2042 } 2043 /* Set pointSF */ 2044 { 2045 PetscSF pointSF; 2046 PetscInt *ilocal, nleaves, q; 2047 PetscSFNode *iremote, *buffer0, *buffer1; 2048 2049 PetscCall(PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1)); 2050 for (p = 0, nleaves = 0; p < *chartSize; ++p) { 2051 if (owners[p] == rank) { 2052 buffer0[p].rank = rank; 2053 } else { 2054 buffer0[p].rank = -1; 2055 nleaves++; 2056 } 2057 buffer0[p].index = p; 2058 } 2059 for (p = 0; p < lsize; ++p) buffer1[p].rank = -1; 2060 PetscCall(PetscSFReduceBegin(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC)); 2061 PetscCall(PetscSFReduceEnd(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC)); 2062 for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1; 2063 PetscCall(PetscSFBcastBegin(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE)); 2064 PetscCall(PetscSFBcastEnd(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE)); 2065 if (PetscDefined(USE_DEBUG)) { 2066 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); 2067 } 2068 PetscCall(PetscMalloc1(nleaves, &ilocal)); 2069 PetscCall(PetscMalloc1(nleaves, &iremote)); 2070 for (p = 0, q = 0; p < *chartSize; ++p) { 2071 if (buffer0[p].rank != rank) { 2072 ilocal[q] = p; 2073 iremote[q].rank = buffer0[p].rank; 2074 iremote[q].index = buffer0[p].index; 2075 q++; 2076 } 2077 } 2078 PetscCheck(q == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching leaf sizes: %" PetscInt_FMT " != %" PetscInt_FMT, q, nleaves); 2079 PetscCall(PetscFree2(buffer0, buffer1)); 2080 PetscCall(PetscSFCreate(comm, &pointSF)); 2081 PetscCall(PetscSFSetFromOptions(pointSF)); 2082 PetscCall(PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 2083 PetscCall(DMSetPointSF(*distdm, pointSF)); 2084 { 2085 DM cdm; 2086 2087 PetscCall(DMGetCoordinateDM(*distdm, &cdm)); 2088 PetscCall(DMSetPointSF(cdm, pointSF)); 2089 } 2090 PetscCall(PetscSFDestroy(&pointSF)); 2091 } 2092 PetscCall(ISRestoreIndices(chartSizesIS, &chartSize)); 2093 PetscCall(ISRestoreIndices(ownersIS, &owners)); 2094 PetscCall(ISRestoreIndices(gpointsIS, &gpoints)); 2095 PetscCall(ISDestroy(&chartSizesIS)); 2096 PetscCall(ISDestroy(&ownersIS)); 2097 PetscCall(ISDestroy(&gpointsIS)); 2098 /* Record that overlap has been manually created. */ 2099 /* This is to pass `DMPlexCheckPointSF()`, which checks that */ 2100 /* pointSF does not contain cells in the leaves if overlap = 0. */ 2101 PetscCall(DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL)); 2102 PetscCall(DMPlexDistributeSetDefault(*distdm, PETSC_FALSE)); 2103 PetscCall(DMPlexReorderSetDefault(*distdm, DM_REORDER_DEFAULT_FALSE)); 2104 PetscCall(PetscLogEventEnd(DMPLEX_DistributionLoad, viewer, 0, 0, 0)); 2105 PetscFunctionReturn(PETSC_SUCCESS); 2106 } 2107 2108 // Serial load of topology 2109 static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf) 2110 { 2111 MPI_Comm comm; 2112 const char *pointsName, *coneSizesName, *conesName, *orientationsName; 2113 IS pointsIS, coneSizesIS, conesIS, orientationsIS; 2114 const PetscInt *points, *coneSizes, *cones, *orientations; 2115 PetscInt *cone, *ornt; 2116 PetscInt dim, N, Np, pEnd, p, q, maxConeSize = 0, c; 2117 PetscMPIInt size, rank; 2118 2119 PetscFunctionBegin; 2120 pointsName = "order"; 2121 coneSizesName = "cones"; 2122 conesName = "cells"; 2123 orientationsName = "orientation"; 2124 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2125 PetscCallMPI(MPI_Comm_size(comm, &size)); 2126 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2127 PetscCall(ISCreate(comm, &pointsIS)); 2128 PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName)); 2129 PetscCall(ISCreate(comm, &coneSizesIS)); 2130 PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName)); 2131 PetscCall(ISCreate(comm, &conesIS)); 2132 PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName)); 2133 PetscCall(ISCreate(comm, &orientationsIS)); 2134 PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName)); 2135 PetscCall(PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim)); 2136 PetscCall(DMSetDimension(dm, dim)); 2137 { 2138 /* Force serial load */ 2139 PetscCall(PetscInfo(dm, "Loading DM %s in serial\n", dm->hdr.name)); 2140 PetscCall(PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np)); 2141 PetscCall(PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0)); 2142 PetscCall(PetscLayoutSetSize(pointsIS->map, Np)); 2143 pEnd = rank == 0 ? Np : 0; 2144 PetscCall(PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np)); 2145 PetscCall(PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0)); 2146 PetscCall(PetscLayoutSetSize(coneSizesIS->map, Np)); 2147 PetscCall(PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N)); 2148 PetscCall(PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0)); 2149 PetscCall(PetscLayoutSetSize(conesIS->map, N)); 2150 PetscCall(PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N)); 2151 PetscCall(PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0)); 2152 PetscCall(PetscLayoutSetSize(orientationsIS->map, N)); 2153 } 2154 PetscCall(ISLoad(pointsIS, viewer)); 2155 PetscCall(ISLoad(coneSizesIS, viewer)); 2156 PetscCall(ISLoad(conesIS, viewer)); 2157 PetscCall(ISLoad(orientationsIS, viewer)); 2158 /* Create Plex */ 2159 PetscCall(DMPlexSetChart(dm, 0, pEnd)); 2160 PetscCall(ISGetIndices(pointsIS, &points)); 2161 PetscCall(ISGetIndices(coneSizesIS, &coneSizes)); 2162 for (p = 0; p < pEnd; ++p) { 2163 PetscCall(DMPlexSetConeSize(dm, points[p], coneSizes[p])); 2164 maxConeSize = PetscMax(maxConeSize, coneSizes[p]); 2165 } 2166 PetscCall(DMSetUp(dm)); 2167 PetscCall(ISGetIndices(conesIS, &cones)); 2168 PetscCall(ISGetIndices(orientationsIS, &orientations)); 2169 PetscCall(PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt)); 2170 for (p = 0, q = 0; p < pEnd; ++p) { 2171 for (c = 0; c < coneSizes[p]; ++c, ++q) { 2172 cone[c] = cones[q]; 2173 ornt[c] = orientations[q]; 2174 } 2175 PetscCall(DMPlexSetCone(dm, points[p], cone)); 2176 PetscCall(DMPlexSetConeOrientation(dm, points[p], ornt)); 2177 } 2178 PetscCall(PetscFree2(cone, ornt)); 2179 /* Create global section migration SF */ 2180 if (sf) { 2181 PetscLayout layout; 2182 PetscInt *globalIndices; 2183 2184 PetscCall(PetscMalloc1(pEnd, &globalIndices)); 2185 /* plex point == globalPointNumber in this case */ 2186 for (p = 0; p < pEnd; ++p) globalIndices[p] = p; 2187 PetscCall(PetscLayoutCreate(comm, &layout)); 2188 PetscCall(PetscLayoutSetSize(layout, Np)); 2189 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 2190 PetscCall(PetscLayoutSetUp(layout)); 2191 PetscCall(PetscSFCreate(comm, sf)); 2192 PetscCall(PetscSFSetFromOptions(*sf)); 2193 PetscCall(PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices)); 2194 PetscCall(PetscLayoutDestroy(&layout)); 2195 PetscCall(PetscFree(globalIndices)); 2196 } 2197 /* Clean-up */ 2198 PetscCall(ISRestoreIndices(pointsIS, &points)); 2199 PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes)); 2200 PetscCall(ISRestoreIndices(conesIS, &cones)); 2201 PetscCall(ISRestoreIndices(orientationsIS, &orientations)); 2202 PetscCall(ISDestroy(&pointsIS)); 2203 PetscCall(ISDestroy(&coneSizesIS)); 2204 PetscCall(ISDestroy(&conesIS)); 2205 PetscCall(ISDestroy(&orientationsIS)); 2206 /* Fill in the rest of the topology structure */ 2207 PetscCall(DMPlexSymmetrize(dm)); 2208 PetscCall(DMPlexStratify(dm)); 2209 PetscFunctionReturn(PETSC_SUCCESS); 2210 } 2211 2212 /* Representation of two DMPlex strata in 0-based global numbering */ 2213 struct _n_PlexLayer { 2214 PetscInt d; 2215 IS conesIS, orientationsIS; 2216 PetscSection coneSizesSection; 2217 PetscLayout vertexLayout; 2218 PetscSF overlapSF, l2gSF; //TODO maybe confusing names (in DMPlex in general) 2219 PetscInt offset, conesOffset, leafOffset; 2220 }; 2221 typedef struct _n_PlexLayer *PlexLayer; 2222 2223 static PetscErrorCode PlexLayerDestroy(PlexLayer *layer) 2224 { 2225 PetscFunctionBegin; 2226 if (!*layer) PetscFunctionReturn(PETSC_SUCCESS); 2227 PetscCall(PetscSectionDestroy(&(*layer)->coneSizesSection)); 2228 PetscCall(ISDestroy(&(*layer)->conesIS)); 2229 PetscCall(ISDestroy(&(*layer)->orientationsIS)); 2230 PetscCall(PetscSFDestroy(&(*layer)->overlapSF)); 2231 PetscCall(PetscSFDestroy(&(*layer)->l2gSF)); 2232 PetscCall(PetscLayoutDestroy(&(*layer)->vertexLayout)); 2233 PetscCall(PetscFree(*layer)); 2234 PetscFunctionReturn(PETSC_SUCCESS); 2235 } 2236 2237 static PetscErrorCode PlexLayerCreate_Private(PlexLayer *layer) 2238 { 2239 PetscFunctionBegin; 2240 PetscCall(PetscNew(layer)); 2241 (*layer)->d = -1; 2242 (*layer)->offset = -1; 2243 (*layer)->conesOffset = -1; 2244 (*layer)->leafOffset = -1; 2245 PetscFunctionReturn(PETSC_SUCCESS); 2246 } 2247 2248 // Parallel load of a depth stratum 2249 static PetscErrorCode PlexLayerLoad_Private(PlexLayer layer, PetscViewer viewer, PetscInt d, PetscLayout pointsLayout) 2250 { 2251 char path[128]; 2252 MPI_Comm comm; 2253 const char *coneSizesName, *conesName, *orientationsName; 2254 IS coneSizesIS, conesIS, orientationsIS; 2255 PetscSection coneSizesSection; 2256 PetscLayout vertexLayout = NULL; 2257 PetscInt s; 2258 2259 PetscFunctionBegin; 2260 coneSizesName = "cone_sizes"; 2261 conesName = "cones"; 2262 orientationsName = "orientations"; 2263 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm)); 2264 2265 /* query size of next lower depth stratum (next lower dimension) */ 2266 if (d > 0) { 2267 PetscInt NVertices; 2268 2269 PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT "/%s", d - 1, coneSizesName)); 2270 PetscCall(PetscViewerHDF5ReadSizes(viewer, path, NULL, &NVertices)); 2271 PetscCall(PetscLayoutCreate(comm, &vertexLayout)); 2272 PetscCall(PetscLayoutSetSize(vertexLayout, NVertices)); 2273 PetscCall(PetscLayoutSetUp(vertexLayout)); 2274 } 2275 2276 PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT, d)); 2277 PetscCall(PetscViewerHDF5PushGroup(viewer, path)); 2278 2279 /* create coneSizesSection from stored IS coneSizes */ 2280 { 2281 const PetscInt *coneSizes; 2282 2283 PetscCall(ISCreate(comm, &coneSizesIS)); 2284 PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName)); 2285 if (pointsLayout) PetscCall(ISSetLayout(coneSizesIS, pointsLayout)); 2286 PetscCall(ISLoad(coneSizesIS, viewer)); 2287 if (!pointsLayout) PetscCall(ISGetLayout(coneSizesIS, &pointsLayout)); 2288 PetscCall(ISGetIndices(coneSizesIS, &coneSizes)); 2289 PetscCall(PetscSectionCreate(comm, &coneSizesSection)); 2290 //TODO different start ? 2291 PetscCall(PetscSectionSetChart(coneSizesSection, 0, pointsLayout->n)); 2292 for (s = 0; s < pointsLayout->n; ++s) PetscCall(PetscSectionSetDof(coneSizesSection, s, coneSizes[s])); 2293 PetscCall(PetscSectionSetUp(coneSizesSection)); 2294 PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes)); 2295 { 2296 PetscLayout tmp = NULL; 2297 2298 /* We need to keep the layout until the end of function */ 2299 PetscCall(PetscLayoutReference((PetscLayout)pointsLayout, &tmp)); 2300 } 2301 PetscCall(ISDestroy(&coneSizesIS)); 2302 } 2303 2304 /* use value layout of coneSizesSection as layout of cones and orientations */ 2305 { 2306 PetscLayout conesLayout; 2307 2308 PetscCall(PetscSectionGetValueLayout(comm, coneSizesSection, &conesLayout)); 2309 PetscCall(ISCreate(comm, &conesIS)); 2310 PetscCall(ISCreate(comm, &orientationsIS)); 2311 PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName)); 2312 PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName)); 2313 PetscCall(PetscLayoutDuplicate(conesLayout, &conesIS->map)); 2314 PetscCall(PetscLayoutDuplicate(conesLayout, &orientationsIS->map)); 2315 PetscCall(ISLoad(conesIS, viewer)); 2316 PetscCall(ISLoad(orientationsIS, viewer)); 2317 PetscCall(PetscLayoutDestroy(&conesLayout)); 2318 } 2319 2320 /* check assertion that layout of points is the same as point layout of coneSizesSection */ 2321 { 2322 PetscLayout pointsLayout0; 2323 PetscBool flg; 2324 2325 PetscCall(PetscSectionGetPointLayout(comm, coneSizesSection, &pointsLayout0)); 2326 PetscCall(PetscLayoutCompare(pointsLayout, pointsLayout0, &flg)); 2327 PetscCheck(flg, comm, PETSC_ERR_PLIB, "points layout != coneSizesSection point layout"); 2328 PetscCall(PetscLayoutDestroy(&pointsLayout0)); 2329 } 2330 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2331 PetscCall(PetscLayoutDestroy(&pointsLayout)); 2332 2333 layer->d = d; 2334 layer->conesIS = conesIS; 2335 layer->coneSizesSection = coneSizesSection; 2336 layer->orientationsIS = orientationsIS; 2337 layer->vertexLayout = vertexLayout; 2338 PetscFunctionReturn(PETSC_SUCCESS); 2339 } 2340 2341 static PetscErrorCode PlexLayerDistribute_Private(PlexLayer layer, PetscSF cellLocalToGlobalSF) 2342 { 2343 IS newConesIS, newOrientationsIS; 2344 PetscSection newConeSizesSection; 2345 MPI_Comm comm; 2346 2347 PetscFunctionBegin; 2348 PetscCall(PetscObjectGetComm((PetscObject)cellLocalToGlobalSF, &comm)); 2349 PetscCall(PetscSectionCreate(comm, &newConeSizesSection)); 2350 //TODO rename to something like ISDistribute() with optional PetscSection argument 2351 PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->conesIS, newConeSizesSection, &newConesIS)); 2352 PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->orientationsIS, newConeSizesSection, &newOrientationsIS)); 2353 2354 PetscCall(PetscObjectSetName((PetscObject)newConeSizesSection, ((PetscObject)layer->coneSizesSection)->name)); 2355 PetscCall(PetscObjectSetName((PetscObject)newConesIS, ((PetscObject)layer->conesIS)->name)); 2356 PetscCall(PetscObjectSetName((PetscObject)newOrientationsIS, ((PetscObject)layer->orientationsIS)->name)); 2357 PetscCall(PetscSectionDestroy(&layer->coneSizesSection)); 2358 PetscCall(ISDestroy(&layer->conesIS)); 2359 PetscCall(ISDestroy(&layer->orientationsIS)); 2360 layer->coneSizesSection = newConeSizesSection; 2361 layer->conesIS = newConesIS; 2362 layer->orientationsIS = newOrientationsIS; 2363 PetscFunctionReturn(PETSC_SUCCESS); 2364 } 2365 2366 //TODO share code with DMPlexBuildFromCellListParallel() 2367 #include <petsc/private/hashseti.h> 2368 static PetscErrorCode PlexLayerCreateSFs_Private(PlexLayer layer, PetscSF *vertexOverlapSF, PetscSF *sfXC) 2369 { 2370 PetscLayout vertexLayout = layer->vertexLayout; 2371 PetscSection coneSection = layer->coneSizesSection; 2372 IS cellVertexData = layer->conesIS; 2373 IS coneOrientations = layer->orientationsIS; 2374 PetscSF vl2gSF, vOverlapSF; 2375 PetscInt *verticesAdj; 2376 PetscInt i, n, numVerticesAdj; 2377 const PetscInt *cvd, *co = NULL; 2378 MPI_Comm comm; 2379 2380 PetscFunctionBegin; 2381 PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm)); 2382 PetscCall(PetscSectionGetStorageSize(coneSection, &n)); 2383 { 2384 PetscInt n0; 2385 2386 PetscCall(ISGetLocalSize(cellVertexData, &n0)); 2387 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); 2388 PetscCall(ISGetIndices(cellVertexData, &cvd)); 2389 } 2390 if (coneOrientations) { 2391 PetscInt n0; 2392 2393 PetscCall(ISGetLocalSize(coneOrientations, &n0)); 2394 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); 2395 PetscCall(ISGetIndices(coneOrientations, &co)); 2396 } 2397 /* Get/check global number of vertices */ 2398 { 2399 PetscInt NVerticesInCells = PETSC_INT_MIN; 2400 2401 /* NVerticesInCells = max(cellVertexData) + 1 */ 2402 for (i = 0; i < n; i++) 2403 if (cvd[i] > NVerticesInCells) NVerticesInCells = cvd[i]; 2404 ++NVerticesInCells; 2405 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, comm)); 2406 2407 if (vertexLayout->n == PETSC_DECIDE && vertexLayout->N == PETSC_DECIDE) vertexLayout->N = NVerticesInCells; 2408 else 2409 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, 2410 vertexLayout->N, NVerticesInCells); 2411 PetscCall(PetscLayoutSetUp(vertexLayout)); 2412 } 2413 /* Find locally unique vertices in cellVertexData */ 2414 { 2415 PetscHSetI vhash; 2416 PetscInt off = 0; 2417 2418 PetscCall(PetscHSetICreate(&vhash)); 2419 for (i = 0; i < n; ++i) PetscCall(PetscHSetIAdd(vhash, cvd[i])); 2420 PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj)); 2421 PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj)); 2422 PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj)); 2423 PetscAssert(off == numVerticesAdj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assertion failed: off == numVerticesAdj (%" PetscInt_FMT " != %" PetscInt_FMT ")", off, numVerticesAdj); 2424 PetscCall(PetscHSetIDestroy(&vhash)); 2425 } 2426 /* We must sort vertices to preserve numbering */ 2427 PetscCall(PetscSortInt(numVerticesAdj, verticesAdj)); 2428 /* Connect locally unique vertices with star forest */ 2429 PetscCall(PetscSFCreateByMatchingIndices(vertexLayout, numVerticesAdj, verticesAdj, NULL, 0, numVerticesAdj, verticesAdj, NULL, 0, &vl2gSF, &vOverlapSF)); 2430 PetscCall(PetscObjectSetName((PetscObject)vOverlapSF, "overlapSF")); 2431 PetscCall(PetscObjectSetName((PetscObject)vl2gSF, "localToGlobalSF")); 2432 2433 PetscCall(PetscFree(verticesAdj)); 2434 *vertexOverlapSF = vOverlapSF; 2435 *sfXC = vl2gSF; 2436 PetscFunctionReturn(PETSC_SUCCESS); 2437 } 2438 2439 static PetscErrorCode PlexLayerCreateCellSFs_Private(PlexLayer layer, PetscSF *cellOverlapSF, PetscSF *cellLocalToGlobalSF) 2440 { 2441 PetscSection coneSection = layer->coneSizesSection; 2442 PetscInt nCells; 2443 MPI_Comm comm; 2444 2445 PetscFunctionBegin; 2446 PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm)); 2447 { 2448 PetscInt cStart; 2449 2450 PetscCall(PetscSectionGetChart(coneSection, &cStart, &nCells)); 2451 PetscCheck(cStart == 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "coneSection must start at 0"); 2452 } 2453 /* Create overlapSF as empty SF with the right number of roots */ 2454 PetscCall(PetscSFCreate(comm, cellOverlapSF)); 2455 PetscCall(PetscSFSetGraph(*cellOverlapSF, nCells, 0, NULL, PETSC_USE_POINTER, NULL, PETSC_USE_POINTER)); 2456 PetscCall(PetscSFSetUp(*cellOverlapSF)); 2457 /* Create localToGlobalSF as identity mapping */ 2458 { 2459 PetscLayout map; 2460 2461 PetscCall(PetscLayoutCreateFromSizes(comm, nCells, PETSC_DECIDE, 1, &map)); 2462 PetscCall(PetscSFCreateFromLayouts(map, map, cellLocalToGlobalSF)); 2463 PetscCall(PetscSFSetUp(*cellLocalToGlobalSF)); 2464 PetscCall(PetscLayoutDestroy(&map)); 2465 } 2466 PetscFunctionReturn(PETSC_SUCCESS); 2467 } 2468 2469 static PetscErrorCode DMPlexTopologyBuildFromLayers_Private(DM dm, PetscInt depth, PlexLayer *layers, IS strataPermutation) 2470 { 2471 const PetscInt *permArr; 2472 PetscInt d, nPoints; 2473 MPI_Comm comm; 2474 2475 PetscFunctionBegin; 2476 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2477 PetscCall(ISGetIndices(strataPermutation, &permArr)); 2478 2479 /* Count points, strata offsets and cones offsets (taking strataPermutation into account) */ 2480 { 2481 PetscInt stratumOffset = 0; 2482 PetscInt conesOffset = 0; 2483 2484 for (d = 0; d <= depth; d++) { 2485 const PetscInt e = permArr[d]; 2486 const PlexLayer l = layers[e]; 2487 PetscInt lo, n, size; 2488 2489 PetscCall(PetscSectionGetChart(l->coneSizesSection, &lo, &n)); 2490 PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &size)); 2491 PetscCheck(lo == 0, comm, PETSC_ERR_PLIB, "starting point should be 0 in coneSizesSection %" PetscInt_FMT, d); 2492 l->offset = stratumOffset; 2493 l->conesOffset = conesOffset; 2494 stratumOffset += n; 2495 conesOffset += size; 2496 } 2497 nPoints = stratumOffset; 2498 } 2499 2500 /* Set interval for all plex points */ 2501 //TODO we should store starting point of plex 2502 PetscCall(DMPlexSetChart(dm, 0, nPoints)); 2503 2504 /* Set up plex coneSection from layer coneSections */ 2505 { 2506 PetscSection coneSection; 2507 2508 PetscCall(DMPlexGetConeSection(dm, &coneSection)); 2509 for (d = 0; d <= depth; d++) { 2510 const PlexLayer l = layers[d]; 2511 PetscInt n, q; 2512 2513 PetscCall(PetscSectionGetChart(l->coneSizesSection, NULL, &n)); 2514 for (q = 0; q < n; q++) { 2515 const PetscInt p = l->offset + q; 2516 PetscInt coneSize; 2517 2518 PetscCall(PetscSectionGetDof(l->coneSizesSection, q, &coneSize)); 2519 PetscCall(PetscSectionSetDof(coneSection, p, coneSize)); 2520 } 2521 } 2522 } 2523 //TODO this is terrible, DMSetUp_Plex() should be DMPlexSetUpSections() or so 2524 PetscCall(DMSetUp(dm)); 2525 2526 /* Renumber cones points from layer-global numbering to plex-local numbering */ 2527 { 2528 PetscInt *cones, *ornts; 2529 2530 PetscCall(DMPlexGetCones(dm, &cones)); 2531 PetscCall(DMPlexGetConeOrientations(dm, &ornts)); 2532 for (d = 1; d <= depth; d++) { 2533 const PlexLayer l = layers[d]; 2534 PetscInt i, lConesSize; 2535 PetscInt *lCones; 2536 const PetscInt *lOrnts; 2537 PetscInt *pCones = &cones[l->conesOffset]; 2538 PetscInt *pOrnts = &ornts[l->conesOffset]; 2539 2540 PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &lConesSize)); 2541 /* Get cones in local plex numbering */ 2542 { 2543 ISLocalToGlobalMapping l2g; 2544 PetscLayout vertexLayout = l->vertexLayout; 2545 PetscSF vertexSF = layers[d - 1]->l2gSF; /* vertices of this layer are cells of previous layer */ 2546 const PetscInt *gCones; 2547 PetscInt lConesSize0; 2548 2549 PetscCall(ISGetLocalSize(l->conesIS, &lConesSize0)); 2550 PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(conesIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize); 2551 PetscCall(ISGetLocalSize(l->orientationsIS, &lConesSize0)); 2552 PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(orientationsIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize); 2553 2554 PetscCall(PetscMalloc1(lConesSize, &lCones)); 2555 PetscCall(ISGetIndices(l->conesIS, &gCones)); 2556 PetscCall(ISLocalToGlobalMappingCreateSF(vertexSF, vertexLayout->rstart, &l2g)); 2557 PetscCall(ISGlobalToLocalMappingApply(l2g, IS_GTOLM_MASK, lConesSize, gCones, &lConesSize0, lCones)); 2558 PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "global to local does not cover all indices (%" PetscInt_FMT " of %" PetscInt_FMT ")", lConesSize0, lConesSize); 2559 PetscCall(ISLocalToGlobalMappingDestroy(&l2g)); 2560 PetscCall(ISRestoreIndices(l->conesIS, &gCones)); 2561 } 2562 PetscCall(ISGetIndices(l->orientationsIS, &lOrnts)); 2563 /* Set cones, need to add stratum offset */ 2564 for (i = 0; i < lConesSize; i++) { 2565 pCones[i] = lCones[i] + layers[d - 1]->offset; /* cone points of current layer are points of previous layer */ 2566 pOrnts[i] = lOrnts[i]; 2567 } 2568 PetscCall(PetscFree(lCones)); 2569 PetscCall(ISRestoreIndices(l->orientationsIS, &lOrnts)); 2570 } 2571 } 2572 PetscCall(DMPlexSymmetrize(dm)); 2573 PetscCall(DMPlexStratify(dm)); 2574 PetscCall(ISRestoreIndices(strataPermutation, &permArr)); 2575 PetscFunctionReturn(PETSC_SUCCESS); 2576 } 2577 2578 static PetscErrorCode PlexLayerConcatenateSFs_Private(MPI_Comm comm, PetscInt depth, PlexLayer layers[], IS strataPermutation, PetscSF *overlapSF, PetscSF *l2gSF) 2579 { 2580 PetscInt d; 2581 PetscSF *osfs, *lsfs; 2582 PetscInt *leafOffsets; 2583 const PetscInt *permArr; 2584 2585 PetscFunctionBegin; 2586 PetscCall(ISGetIndices(strataPermutation, &permArr)); 2587 PetscCall(PetscCalloc3(depth + 1, &osfs, depth + 1, &lsfs, depth + 1, &leafOffsets)); 2588 for (d = 0; d <= depth; d++) { 2589 const PetscInt e = permArr[d]; 2590 2591 PetscAssert(e == layers[e]->d, PETSC_COMM_SELF, PETSC_ERR_PLIB, "assertion: e == layers[e]->d"); 2592 osfs[d] = layers[e]->overlapSF; 2593 lsfs[d] = layers[e]->l2gSF; 2594 leafOffsets[d] = layers[e]->offset; 2595 } 2596 PetscCall(PetscSFConcatenate(comm, depth + 1, osfs, PETSCSF_CONCATENATE_ROOTMODE_LOCAL, leafOffsets, overlapSF)); 2597 PetscCall(PetscSFConcatenate(comm, depth + 1, lsfs, PETSCSF_CONCATENATE_ROOTMODE_GLOBAL, leafOffsets, l2gSF)); 2598 PetscCall(PetscFree3(osfs, lsfs, leafOffsets)); 2599 PetscCall(ISRestoreIndices(strataPermutation, &permArr)); 2600 PetscFunctionReturn(PETSC_SUCCESS); 2601 } 2602 2603 // Parallel load of topology 2604 static PetscErrorCode DMPlexTopologyLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF *sfXC) 2605 { 2606 PlexLayer *layers; 2607 IS strataPermutation; 2608 PetscLayout pointsLayout; 2609 PetscInt depth; 2610 PetscInt d; 2611 MPI_Comm comm; 2612 2613 PetscFunctionBegin; 2614 { 2615 PetscInt dim; 2616 2617 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "depth", PETSC_INT, NULL, &depth)); 2618 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "cell_dim", PETSC_INT, NULL, &dim)); 2619 PetscCall(DMSetDimension(dm, dim)); 2620 } 2621 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2622 2623 PetscCall(PetscInfo(dm, "Loading DM %s in parallel\n", dm->hdr.name)); 2624 { 2625 IS spOnComm; 2626 2627 PetscCall(ISCreate(comm, &spOnComm)); 2628 PetscCall(PetscObjectSetName((PetscObject)spOnComm, "permutation")); 2629 PetscCall(ISLoad(spOnComm, viewer)); 2630 /* have the same serial IS on every rank */ 2631 PetscCall(ISAllGather(spOnComm, &strataPermutation)); 2632 PetscCall(PetscObjectSetName((PetscObject)strataPermutation, ((PetscObject)spOnComm)->name)); 2633 PetscCall(ISDestroy(&spOnComm)); 2634 } 2635 2636 /* Create layers, load raw data for each layer */ 2637 PetscCall(PetscViewerHDF5PushGroup(viewer, "strata")); 2638 PetscCall(PetscMalloc1(depth + 1, &layers)); 2639 for (d = depth, pointsLayout = NULL; d >= 0; pointsLayout = layers[d]->vertexLayout, d--) { 2640 PetscCall(PlexLayerCreate_Private(&layers[d])); 2641 PetscCall(PlexLayerLoad_Private(layers[d], viewer, d, pointsLayout)); 2642 } 2643 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */ 2644 2645 for (d = depth; d >= 0; d--) { 2646 /* Redistribute cells and vertices for each applicable layer */ 2647 if (d < depth) PetscCall(PlexLayerDistribute_Private(layers[d], layers[d]->l2gSF)); 2648 /* Create vertex overlap SF and vertex localToGlobal SF for each applicable layer */ 2649 if (d > 0) PetscCall(PlexLayerCreateSFs_Private(layers[d], &layers[d - 1]->overlapSF, &layers[d - 1]->l2gSF)); 2650 } 2651 /* Build trivial SFs for the cell layer as well */ 2652 PetscCall(PlexLayerCreateCellSFs_Private(layers[depth], &layers[depth]->overlapSF, &layers[depth]->l2gSF)); 2653 2654 /* Build DMPlex topology from the layers */ 2655 PetscCall(DMPlexTopologyBuildFromLayers_Private(dm, depth, layers, strataPermutation)); 2656 2657 /* Build overall point SF alias overlap SF */ 2658 { 2659 PetscSF overlapSF; 2660 2661 PetscCall(PlexLayerConcatenateSFs_Private(comm, depth, layers, strataPermutation, &overlapSF, sfXC)); 2662 PetscCall(DMSetPointSF(dm, overlapSF)); 2663 PetscCall(PetscSFDestroy(&overlapSF)); 2664 } 2665 2666 for (d = depth; d >= 0; d--) PetscCall(PlexLayerDestroy(&layers[d])); 2667 PetscCall(PetscFree(layers)); 2668 PetscCall(ISDestroy(&strataPermutation)); 2669 PetscFunctionReturn(PETSC_SUCCESS); 2670 } 2671 2672 PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC) 2673 { 2674 DMPlexStorageVersion version; 2675 const char *topologydm_name; 2676 char group[PETSC_MAX_PATH_LEN]; 2677 PetscSF sfwork = NULL; 2678 2679 PetscFunctionBegin; 2680 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 2681 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version)); 2682 if (DMPlexStorageVersionGE(version, 2, 0, 0)) { 2683 PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name)); 2684 } else { 2685 PetscCall(PetscStrncpy(group, "/", sizeof(group))); 2686 } 2687 PetscCall(PetscViewerHDF5PushGroup(viewer, group)); 2688 2689 PetscCall(PetscViewerHDF5PushGroup(viewer, "topology")); 2690 if (version->major < 3) { 2691 PetscCall(DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork)); 2692 } else { 2693 /* since DMPlexStorageVersion 3.0.0 */ 2694 PetscCall(DMPlexTopologyLoad_HDF5_Private(dm, viewer, &sfwork)); 2695 } 2696 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */ 2697 2698 if (DMPlexStorageVersionGE(version, 2, 0, 0)) { 2699 DM distdm; 2700 PetscSF distsf; 2701 const char *distribution_name; 2702 PetscBool exists; 2703 2704 PetscCall(DMPlexDistributionGetName(dm, &distribution_name)); 2705 /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */ 2706 PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions")); 2707 PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &exists)); 2708 if (exists) { 2709 PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name)); 2710 PetscCall(DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm)); 2711 if (distdm) { 2712 PetscCall(DMPlexReplace_Internal(dm, &distdm)); 2713 PetscCall(PetscSFDestroy(&sfwork)); 2714 sfwork = distsf; 2715 } 2716 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */ 2717 } 2718 PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */ 2719 } 2720 if (sfXC) { 2721 *sfXC = sfwork; 2722 } else { 2723 PetscCall(PetscSFDestroy(&sfwork)); 2724 } 2725 2726 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2727 PetscFunctionReturn(PETSC_SUCCESS); 2728 } 2729 2730 /* If the file is old, it not only has different path to the coordinates, but */ 2731 /* does not contain coordinateDMs, so must fall back to the old implementation. */ 2732 static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer) 2733 { 2734 PetscSection coordSection; 2735 Vec coordinates; 2736 PetscReal lengthScale; 2737 PetscInt spatialDim, N, numVertices, vStart, vEnd, v; 2738 PetscMPIInt rank; 2739 2740 PetscFunctionBegin; 2741 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 2742 /* Read geometry */ 2743 PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry")); 2744 PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates)); 2745 PetscCall(PetscObjectSetName((PetscObject)coordinates, "vertices")); 2746 { 2747 /* Force serial load */ 2748 PetscCall(PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N)); 2749 PetscCall(VecSetSizes(coordinates, rank == 0 ? N : 0, N)); 2750 PetscCall(VecSetBlockSize(coordinates, spatialDim)); 2751 } 2752 PetscCall(VecLoad(coordinates, viewer)); 2753 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2754 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 2755 PetscCall(VecScale(coordinates, 1.0 / lengthScale)); 2756 PetscCall(VecGetLocalSize(coordinates, &numVertices)); 2757 PetscCall(VecGetBlockSize(coordinates, &spatialDim)); 2758 numVertices /= spatialDim; 2759 /* Create coordinates */ 2760 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 2761 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); 2762 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 2763 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 2764 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spatialDim)); 2765 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 2766 for (v = vStart; v < vEnd; ++v) { 2767 PetscCall(PetscSectionSetDof(coordSection, v, spatialDim)); 2768 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spatialDim)); 2769 } 2770 PetscCall(PetscSectionSetUp(coordSection)); 2771 PetscCall(DMSetCoordinates(dm, coordinates)); 2772 PetscCall(VecDestroy(&coordinates)); 2773 PetscFunctionReturn(PETSC_SUCCESS); 2774 } 2775 2776 PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC) 2777 { 2778 DMPlexStorageVersion version; 2779 DM cdm; 2780 Vec coords; 2781 PetscInt blockSize; 2782 PetscReal lengthScale; 2783 PetscSF lsf; 2784 const char *topologydm_name; 2785 char *coordinatedm_name, *coordinates_name; 2786 2787 PetscFunctionBegin; 2788 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version)); 2789 if (!DMPlexStorageVersionGE(version, 2, 0, 0)) { 2790 PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer)); 2791 PetscFunctionReturn(PETSC_SUCCESS); 2792 } 2793 /* else: since DMPlexStorageVersion 2.0.0 */ 2794 PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load"); 2795 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 2796 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 2797 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 2798 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name)); 2799 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name)); 2800 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2801 PetscCall(PetscViewerHDF5PopGroup(viewer)); 2802 PetscCall(DMGetCoordinateDM(dm, &cdm)); 2803 PetscCall(PetscObjectSetName((PetscObject)cdm, coordinatedm_name)); 2804 PetscCall(PetscFree(coordinatedm_name)); 2805 /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */ 2806 PetscCall(DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf)); 2807 PetscCall(DMCreateLocalVector(cdm, &coords)); 2808 PetscCall(PetscObjectSetName((PetscObject)coords, coordinates_name)); 2809 PetscCall(PetscFree(coordinates_name)); 2810 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE)); 2811 PetscCall(DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords)); 2812 PetscCall(PetscViewerPopFormat(viewer)); 2813 PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale)); 2814 PetscCall(VecScale(coords, 1.0 / lengthScale)); 2815 PetscCall(DMSetCoordinatesLocal(dm, coords)); 2816 PetscCall(VecGetBlockSize(coords, &blockSize)); 2817 PetscCall(DMSetCoordinateDim(dm, blockSize)); 2818 PetscCall(VecDestroy(&coords)); 2819 PetscCall(PetscSFDestroy(&lsf)); 2820 PetscFunctionReturn(PETSC_SUCCESS); 2821 } 2822 2823 PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer) 2824 { 2825 DMPlexStorageVersion version; 2826 2827 PetscFunctionBegin; 2828 PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version)); 2829 PetscCall(PetscInfo(dm, "Loading DM %s storage version %d.%d.%d\n", dm->hdr.name, version->major, version->minor, version->subminor)); 2830 if (!DMPlexStorageVersionGE(version, 2, 0, 0)) { 2831 PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL)); 2832 PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL)); 2833 PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer)); 2834 } else { 2835 PetscSF sfXC; 2836 2837 /* since DMPlexStorageVersion 2.0.0 */ 2838 PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC)); 2839 PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC)); 2840 PetscCall(DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC)); 2841 PetscCall(PetscSFDestroy(&sfXC)); 2842 } 2843 PetscFunctionReturn(PETSC_SUCCESS); 2844 } 2845 2846 static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF) 2847 { 2848 MPI_Comm comm; 2849 PetscInt pStart, pEnd, p, m; 2850 PetscInt *goffs, *ilocal; 2851 PetscBool rootIncludeConstraints, leafIncludeConstraints; 2852 2853 PetscFunctionBegin; 2854 PetscCall(PetscObjectGetComm((PetscObject)leafSection, &comm)); 2855 PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd)); 2856 PetscCall(PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints)); 2857 PetscCall(PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints)); 2858 if (rootIncludeConstraints && leafIncludeConstraints) PetscCall(PetscSectionGetStorageSize(leafSection, &m)); 2859 else PetscCall(PetscSectionGetConstrainedStorageSize(leafSection, &m)); 2860 PetscCall(PetscMalloc1(m, &ilocal)); 2861 PetscCall(PetscMalloc1(m, &goffs)); 2862 /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */ 2863 /* for the top-level section (not for each field), so one must have */ 2864 /* rootSection->pointMajor == PETSC_TRUE. */ 2865 PetscCheck(rootSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering"); 2866 /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */ 2867 PetscCheck(leafSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering"); 2868 for (p = pStart, m = 0; p < pEnd; ++p) { 2869 PetscInt dof, cdof, i, j, off, goff; 2870 const PetscInt *cinds; 2871 2872 PetscCall(PetscSectionGetDof(leafSection, p, &dof)); 2873 if (dof < 0) continue; 2874 goff = globalOffsets[p - pStart]; 2875 PetscCall(PetscSectionGetOffset(leafSection, p, &off)); 2876 PetscCall(PetscSectionGetConstraintDof(leafSection, p, &cdof)); 2877 PetscCall(PetscSectionGetConstraintIndices(leafSection, p, &cinds)); 2878 for (i = 0, j = 0; i < dof; ++i) { 2879 PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]); 2880 2881 if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) { 2882 ilocal[m] = off++; 2883 goffs[m++] = goff++; 2884 } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off; 2885 else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff; 2886 if (constrained) ++j; 2887 } 2888 } 2889 PetscCall(PetscSFCreate(comm, sectionSF)); 2890 PetscCall(PetscSFSetFromOptions(*sectionSF)); 2891 PetscCall(PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs)); 2892 PetscCall(PetscFree(goffs)); 2893 PetscFunctionReturn(PETSC_SUCCESS); 2894 } 2895 2896 PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf) 2897 { 2898 MPI_Comm comm; 2899 PetscMPIInt size, rank; 2900 const char *topologydm_name; 2901 const char *sectiondm_name; 2902 PetscSection sectionA, sectionB; 2903 PetscBool has; 2904 PetscInt nX, n, i; 2905 PetscSF sfAB; 2906 2907 PetscFunctionBegin; 2908 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 2909 PetscCallMPI(MPI_Comm_size(comm, &size)); 2910 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 2911 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 2912 PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name)); 2913 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 2914 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 2915 PetscCall(PetscViewerHDF5PushGroup(viewer, "dms")); 2916 PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name)); 2917 /* A: on-disk points */ 2918 /* X: list of global point numbers, [0, NX) */ 2919 /* B: plex points */ 2920 /* Load raw section (sectionA) */ 2921 PetscCall(PetscSectionCreate(comm, §ionA)); 2922 PetscCall(PetscViewerHDF5HasGroup(viewer, "section", &has)); 2923 if (has) PetscCall(PetscSectionLoad(sectionA, viewer)); 2924 else { 2925 // TODO If section is missing, create the default affine section with dim dofs on each vertex. Use PetscSplitOwnership() to split vertices. 2926 // How do I know the total number of vertices? 2927 PetscInt dim, Nf = 1, Nv, nv = PETSC_DECIDE; 2928 2929 PetscCall(DMGetDimension(dm, &dim)); 2930 PetscCall(DMPlexGetDepthStratumGlobalSize(dm, 0, &Nv)); 2931 PetscCall(PetscSectionSetNumFields(sectionA, Nf)); 2932 PetscCall(PetscSectionSetFieldName(sectionA, 0, "Cartesian")); 2933 PetscCall(PetscSectionSetFieldComponents(sectionA, 0, dim)); 2934 for (PetscInt c = 0; c < dim; ++c) { 2935 char axis = 'X' + (char)c; 2936 2937 PetscCall(PetscSectionSetComponentName(sectionA, 0, c, &axis)); 2938 } 2939 PetscCall(PetscSplitOwnership(comm, &nv, &Nv)); 2940 PetscCall(PetscSectionSetChart(sectionA, 0, nv)); 2941 for (PetscInt p = 0; p < nv; ++p) { 2942 PetscCall(PetscSectionSetDof(sectionA, p, dim)); 2943 PetscCall(PetscSectionSetFieldDof(sectionA, p, 0, dim)); 2944 } 2945 PetscCall(PetscSectionSetUp(sectionA)); 2946 } 2947 PetscCall(PetscSectionGetChart(sectionA, NULL, &n)); 2948 /* Create sfAB: A -> B */ 2949 #if defined(PETSC_USE_DEBUG) 2950 { 2951 PetscInt N, N1; 2952 2953 PetscCall(PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1)); 2954 PetscCallMPI(MPIU_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm)); 2955 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); 2956 } 2957 #endif 2958 { 2959 IS orderIS; 2960 const PetscInt *gpoints; 2961 PetscSF sfXA, sfAX; 2962 PetscLayout layout; 2963 PetscSFNode *owners, *buffer; 2964 PetscInt nleaves; 2965 PetscInt *ilocal; 2966 PetscSFNode *iremote; 2967 2968 /* Create sfAX: A -> X */ 2969 PetscCall(ISCreate(comm, &orderIS)); 2970 PetscCall(PetscObjectSetName((PetscObject)orderIS, "order")); 2971 PetscCall(PetscLayoutSetLocalSize(orderIS->map, n)); 2972 PetscCall(ISLoad(orderIS, viewer)); 2973 PetscCall(PetscLayoutCreate(comm, &layout)); 2974 PetscCall(PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL)); 2975 PetscCall(PetscLayoutSetLocalSize(layout, nX)); 2976 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 2977 PetscCall(PetscLayoutSetUp(layout)); 2978 PetscCall(PetscSFCreate(comm, &sfXA)); 2979 PetscCall(ISGetIndices(orderIS, &gpoints)); 2980 PetscCall(PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints)); 2981 PetscCall(ISRestoreIndices(orderIS, &gpoints)); 2982 PetscCall(ISDestroy(&orderIS)); 2983 PetscCall(PetscLayoutDestroy(&layout)); 2984 PetscCall(PetscMalloc1(n, &owners)); 2985 PetscCall(PetscMalloc1(nX, &buffer)); 2986 for (i = 0; i < n; ++i) { 2987 owners[i].rank = rank; 2988 owners[i].index = i; 2989 } 2990 for (i = 0; i < nX; ++i) { 2991 buffer[i].rank = -1; 2992 buffer[i].index = -1; 2993 } 2994 PetscCall(PetscSFReduceBegin(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC)); 2995 PetscCall(PetscSFReduceEnd(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC)); 2996 PetscCall(PetscSFDestroy(&sfXA)); 2997 PetscCall(PetscFree(owners)); 2998 for (i = 0, nleaves = 0; i < nX; ++i) 2999 if (buffer[i].rank >= 0) nleaves++; 3000 PetscCall(PetscMalloc1(nleaves, &ilocal)); 3001 PetscCall(PetscMalloc1(nleaves, &iremote)); 3002 for (i = 0, nleaves = 0; i < nX; ++i) { 3003 if (buffer[i].rank >= 0) { 3004 ilocal[nleaves] = i; 3005 iremote[nleaves].rank = buffer[i].rank; 3006 iremote[nleaves].index = buffer[i].index; 3007 nleaves++; 3008 } 3009 } 3010 PetscCall(PetscFree(buffer)); 3011 PetscCall(PetscSFCreate(comm, &sfAX)); 3012 PetscCall(PetscSFSetFromOptions(sfAX)); 3013 PetscCall(PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 3014 PetscCall(PetscSFCompose(sfAX, sfXB, &sfAB)); 3015 PetscCall(PetscSFDestroy(&sfAX)); 3016 } 3017 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3018 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3019 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3020 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3021 /* Create plex section (sectionB) */ 3022 PetscCall(DMGetLocalSection(sectiondm, §ionB)); 3023 if (lsf || gsf) { 3024 PetscLayout layout; 3025 PetscInt M, m; 3026 PetscInt *offsetsA; 3027 PetscBool includesConstraintsA; 3028 3029 PetscCall(PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB)); 3030 PetscCall(PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA)); 3031 if (includesConstraintsA) PetscCall(PetscSectionGetStorageSize(sectionA, &m)); 3032 else PetscCall(PetscSectionGetConstrainedStorageSize(sectionA, &m)); 3033 PetscCallMPI(MPIU_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm)); 3034 PetscCall(PetscLayoutCreate(comm, &layout)); 3035 PetscCall(PetscLayoutSetSize(layout, M)); 3036 PetscCall(PetscLayoutSetUp(layout)); 3037 if (lsf) { 3038 PetscSF lsfABdata; 3039 3040 PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata)); 3041 *lsf = lsfABdata; 3042 } 3043 if (gsf) { 3044 PetscSection gsectionB, gsectionB1; 3045 PetscBool includesConstraintsB; 3046 PetscSF gsfABdata, pointsf; 3047 3048 PetscCall(DMGetGlobalSection(sectiondm, &gsectionB1)); 3049 PetscCall(PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB)); 3050 PetscCall(DMGetPointSF(sectiondm, &pointsf)); 3051 PetscCall(PetscSectionCreateGlobalSection(sectionB, pointsf, PETSC_TRUE, includesConstraintsB, PETSC_TRUE, &gsectionB)); 3052 PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata)); 3053 PetscCall(PetscSectionDestroy(&gsectionB)); 3054 *gsf = gsfABdata; 3055 } 3056 PetscCall(PetscLayoutDestroy(&layout)); 3057 PetscCall(PetscFree(offsetsA)); 3058 } else { 3059 PetscCall(PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB)); 3060 } 3061 PetscCall(PetscSFDestroy(&sfAB)); 3062 PetscCall(PetscSectionDestroy(§ionA)); 3063 PetscFunctionReturn(PETSC_SUCCESS); 3064 } 3065 3066 PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec) 3067 { 3068 MPI_Comm comm; 3069 const char *topologydm_name; 3070 const char *sectiondm_name; 3071 const char *vec_name; 3072 Vec vecA; 3073 PetscInt mA, m, bs; 3074 const PetscInt *ilocal; 3075 const PetscScalar *src; 3076 PetscScalar *dest; 3077 3078 PetscFunctionBegin; 3079 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 3080 PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name)); 3081 PetscCall(PetscObjectGetName((PetscObject)sectiondm, §iondm_name)); 3082 PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name)); 3083 PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies")); 3084 PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name)); 3085 PetscCall(PetscViewerHDF5PushGroup(viewer, "dms")); 3086 PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name)); 3087 PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs")); 3088 PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name)); 3089 PetscCall(VecCreate(comm, &vecA)); 3090 PetscCall(PetscObjectSetName((PetscObject)vecA, vec_name)); 3091 PetscCall(PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL)); 3092 /* Check consistency */ 3093 { 3094 PetscSF pointsf, pointsf1; 3095 PetscInt m1, i, j; 3096 3097 PetscCall(DMGetPointSF(dm, &pointsf)); 3098 PetscCall(DMGetPointSF(sectiondm, &pointsf1)); 3099 PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm"); 3100 #if defined(PETSC_USE_DEBUG) 3101 { 3102 PetscInt MA, MA1; 3103 3104 PetscCallMPI(MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm)); 3105 PetscCall(PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1)); 3106 PetscCheck(MA1 == MA, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%" PetscInt_FMT ") != On-disk vector data size (%" PetscInt_FMT ")", MA, MA1); 3107 } 3108 #endif 3109 PetscCall(VecGetLocalSize(vec, &m1)); 3110 PetscCheck(m1 >= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%" PetscInt_FMT ") < SF leaf size (%" PetscInt_FMT ")", m1, m); 3111 for (i = 0; i < m; ++i) { 3112 j = ilocal ? ilocal[i] : i; 3113 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); 3114 } 3115 } 3116 PetscCall(VecSetSizes(vecA, mA, PETSC_DECIDE)); 3117 PetscCall(VecLoad(vecA, viewer)); 3118 PetscCall(VecGetArrayRead(vecA, &src)); 3119 PetscCall(VecGetArray(vec, &dest)); 3120 PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE)); 3121 PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE)); 3122 PetscCall(VecRestoreArray(vec, &dest)); 3123 PetscCall(VecRestoreArrayRead(vecA, &src)); 3124 PetscCall(VecDestroy(&vecA)); 3125 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs)); 3126 PetscCall(VecSetBlockSize(vec, bs)); 3127 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3128 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3129 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3130 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3131 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3132 PetscCall(PetscViewerHDF5PopGroup(viewer)); 3133 PetscFunctionReturn(PETSC_SUCCESS); 3134 } 3135