1 #include <petsc/private/sectionimpl.h> /*I "petscsection.h" I*/ 2 #include <petscsf.h> 3 #include <petscis.h> 4 #include <petscviewerhdf5.h> 5 #include <petsclayouthdf5.h> 6 7 static PetscErrorCode PetscSectionView_HDF5_SingleField(PetscSection s, PetscViewer viewer) 8 { 9 MPI_Comm comm; 10 PetscInt pStart, pEnd, p, n; 11 PetscBool hasConstraints, includesConstraints; 12 IS dofIS, offIS, cdofIS, coffIS, cindIS; 13 PetscInt *dofs, *offs, *cdofs, *coffs, *cinds, dof, cdof, m, moff, i; 14 15 PetscFunctionBegin; 16 PetscCall(PetscObjectGetComm((PetscObject)s, &comm)); 17 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 18 hasConstraints = (s->bc) ? PETSC_TRUE : PETSC_FALSE; 19 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hasConstraints, 1, MPIU_BOOL, MPI_LOR, comm)); 20 for (p = pStart, n = 0, m = 0; p < pEnd; ++p) { 21 PetscCall(PetscSectionGetDof(s, p, &dof)); 22 if (dof >= 0) { 23 if (hasConstraints) { 24 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 25 m += cdof; 26 } 27 n++; 28 } 29 } 30 PetscCall(PetscMalloc1(n, &dofs)); 31 PetscCall(PetscMalloc1(n, &offs)); 32 if (hasConstraints) { 33 PetscCall(PetscMalloc1(n, &cdofs)); 34 PetscCall(PetscMalloc1(n, &coffs)); 35 PetscCall(PetscMalloc1(m, &cinds)); 36 } 37 for (p = pStart, n = 0, m = 0; p < pEnd; ++p) { 38 PetscCall(PetscSectionGetDof(s, p, &dof)); 39 if (dof >= 0) { 40 dofs[n] = dof; 41 PetscCall(PetscSectionGetOffset(s, p, &offs[n])); 42 if (hasConstraints) { 43 const PetscInt *cpinds; 44 45 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 46 PetscCall(PetscSectionGetConstraintIndices(s, p, &cpinds)); 47 cdofs[n] = cdof; 48 coffs[n] = m; 49 for (i = 0; i < cdof; ++i) cinds[m++] = cpinds[i]; 50 } 51 n++; 52 } 53 } 54 if (hasConstraints) { 55 PetscCallMPI(MPI_Scan(&m, &moff, 1, MPIU_INT, MPI_SUM, comm)); 56 moff -= m; 57 for (p = 0; p < n; ++p) coffs[p] += moff; 58 } 59 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, (void *)&hasConstraints)); 60 PetscCall(PetscSectionGetIncludesConstraints(s, &includesConstraints)); 61 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, (void *)&includesConstraints)); 62 PetscCall(ISCreateGeneral(comm, n, dofs, PETSC_OWN_POINTER, &dofIS)); 63 PetscCall(PetscObjectSetName((PetscObject)dofIS, "atlasDof")); 64 PetscCall(ISView(dofIS, viewer)); 65 PetscCall(ISDestroy(&dofIS)); 66 PetscCall(ISCreateGeneral(comm, n, offs, PETSC_OWN_POINTER, &offIS)); 67 PetscCall(PetscObjectSetName((PetscObject)offIS, "atlasOff")); 68 PetscCall(ISView(offIS, viewer)); 69 PetscCall(ISDestroy(&offIS)); 70 if (hasConstraints) { 71 PetscCall(PetscViewerHDF5PushGroup(viewer, "bc")); 72 PetscCall(ISCreateGeneral(comm, n, cdofs, PETSC_OWN_POINTER, &cdofIS)); 73 PetscCall(PetscObjectSetName((PetscObject)cdofIS, "atlasDof")); 74 PetscCall(ISView(cdofIS, viewer)); 75 PetscCall(ISDestroy(&cdofIS)); 76 PetscCall(ISCreateGeneral(comm, n, coffs, PETSC_OWN_POINTER, &coffIS)); 77 PetscCall(PetscObjectSetName((PetscObject)coffIS, "atlasOff")); 78 PetscCall(ISView(coffIS, viewer)); 79 PetscCall(ISDestroy(&coffIS)); 80 PetscCall(PetscViewerHDF5PopGroup(viewer)); 81 PetscCall(ISCreateGeneral(comm, m, cinds, PETSC_OWN_POINTER, &cindIS)); 82 PetscCall(PetscObjectSetName((PetscObject)cindIS, "bcIndices")); 83 PetscCall(ISView(cindIS, viewer)); 84 PetscCall(ISDestroy(&cindIS)); 85 } 86 PetscFunctionReturn(PETSC_SUCCESS); 87 } 88 89 PetscErrorCode PetscSectionView_HDF5_Internal(PetscSection s, PetscViewer viewer) 90 { 91 PetscInt numFields, f; 92 93 PetscFunctionBegin; 94 PetscCall(PetscViewerHDF5PushGroup(viewer, "section")); 95 PetscCall(PetscSectionGetNumFields(s, &numFields)); 96 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numFields", PETSC_INT, (void *)&numFields)); 97 PetscCall(PetscSectionView_HDF5_SingleField(s, viewer)); 98 for (f = 0; f < numFields; ++f) { 99 char fname[PETSC_MAX_PATH_LEN]; 100 const char *fieldName; 101 PetscInt fieldComponents, c; 102 103 PetscCall(PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f)); 104 PetscCall(PetscViewerHDF5PushGroup(viewer, fname)); 105 PetscCall(PetscSectionGetFieldName(s, f, &fieldName)); 106 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldName", PETSC_STRING, fieldName)); 107 PetscCall(PetscSectionGetFieldComponents(s, f, &fieldComponents)); 108 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldComponents", PETSC_INT, (void *)&fieldComponents)); 109 for (c = 0; c < fieldComponents; ++c) { 110 char cname[PETSC_MAX_PATH_LEN]; 111 const char *componentName; 112 113 PetscCall(PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c)); 114 PetscCall(PetscViewerHDF5PushGroup(viewer, cname)); 115 PetscCall(PetscSectionGetComponentName(s, f, c, &componentName)); 116 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "componentName", PETSC_STRING, componentName)); 117 PetscCall(PetscViewerHDF5PopGroup(viewer)); 118 } 119 PetscCall(PetscSectionView_HDF5_SingleField(s->field[f], viewer)); 120 PetscCall(PetscViewerHDF5PopGroup(viewer)); 121 } 122 PetscCall(PetscViewerHDF5PopGroup(viewer)); 123 PetscFunctionReturn(PETSC_SUCCESS); 124 } 125 126 static PetscErrorCode PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(PetscSection s, IS cindIS, IS coffIS) 127 { 128 MPI_Comm comm; 129 PetscInt pStart, pEnd, p, M, m, i, cdof; 130 const PetscInt *data; 131 PetscInt *cinds; 132 const PetscInt *coffs; 133 PetscInt *coffsets; 134 PetscSF sf; 135 PetscLayout layout; 136 137 PetscFunctionBegin; 138 PetscCall(PetscObjectGetComm((PetscObject)s, &comm)); 139 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 140 PetscCall(ISGetSize(cindIS, &M)); 141 PetscCall(ISGetLocalSize(cindIS, &m)); 142 PetscCall(PetscMalloc1(m, &coffsets)); 143 PetscCall(ISGetIndices(coffIS, &coffs)); 144 for (p = pStart, m = 0; p < pEnd; ++p) { 145 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 146 for (i = 0; i < cdof; ++i) coffsets[m++] = coffs[p - pStart] + i; 147 } 148 PetscCall(ISRestoreIndices(coffIS, &coffs)); 149 PetscCall(PetscSFCreate(comm, &sf)); 150 PetscCall(PetscLayoutCreate(comm, &layout)); 151 PetscCall(PetscLayoutSetSize(layout, M)); 152 PetscCall(PetscLayoutSetLocalSize(layout, m)); 153 PetscCall(PetscLayoutSetBlockSize(layout, 1)); 154 PetscCall(PetscLayoutSetUp(layout)); 155 PetscCall(PetscSFSetGraphLayout(sf, layout, m, NULL, PETSC_OWN_POINTER, coffsets)); 156 PetscCall(PetscLayoutDestroy(&layout)); 157 PetscCall(PetscFree(coffsets)); 158 PetscCall(PetscMalloc1(m, &cinds)); 159 PetscCall(ISGetIndices(cindIS, &data)); 160 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, data, cinds, MPI_REPLACE)); 161 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, data, cinds, MPI_REPLACE)); 162 PetscCall(ISRestoreIndices(cindIS, &data)); 163 PetscCall(PetscSFDestroy(&sf)); 164 PetscCall(PetscSectionSetUpBC(s)); 165 for (p = pStart, m = 0; p < pEnd; ++p) { 166 PetscCall(PetscSectionGetConstraintDof(s, p, &cdof)); 167 PetscCall(PetscSectionSetConstraintIndices(s, p, PetscSafePointerPlusOffset(cinds, m))); 168 m += cdof; 169 } 170 PetscCall(PetscFree(cinds)); 171 PetscFunctionReturn(PETSC_SUCCESS); 172 } 173 174 static PetscErrorCode PetscSectionLoad_HDF5_SingleField(PetscSection s, PetscViewer viewer) 175 { 176 MPI_Comm comm; 177 PetscInt pStart, pEnd, p, N, n, M, m; 178 #if defined(PETSC_USE_DEBUG) 179 PetscInt N1, M1; 180 #endif 181 PetscBool hasConstraints, includesConstraints; 182 IS dofIS, offIS, cdofIS, coffIS, cindIS; 183 const PetscInt *dofs, *offs, *cdofs; 184 PetscLayout map; 185 186 PetscFunctionBegin; 187 PetscCall(PetscObjectGetComm((PetscObject)s, &comm)); 188 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, NULL, (void *)&includesConstraints)); 189 PetscCall(PetscSectionSetIncludesConstraints(s, includesConstraints)); 190 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 191 n = pEnd - pStart; 192 #if defined(PETSC_USE_DEBUG) 193 PetscCallMPI(MPIU_Allreduce(&n, &N1, 1, MPIU_INT, MPI_SUM, comm)); 194 #endif 195 PetscCall(ISCreate(comm, &dofIS)); 196 PetscCall(PetscObjectSetName((PetscObject)dofIS, "atlasDof")); 197 PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N)); 198 #if defined(PETSC_USE_DEBUG) 199 PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->atlasDof: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n); 200 #endif 201 PetscCall(ISGetLayout(dofIS, &map)); 202 PetscCall(PetscLayoutSetSize(map, N)); 203 PetscCall(PetscLayoutSetLocalSize(map, n)); 204 PetscCall(ISLoad(dofIS, viewer)); 205 PetscCall(ISCreate(comm, &offIS)); 206 PetscCall(PetscObjectSetName((PetscObject)offIS, "atlasOff")); 207 PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N)); 208 #if defined(PETSC_USE_DEBUG) 209 PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->atlasOff: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n); 210 #endif 211 PetscCall(ISGetLayout(offIS, &map)); 212 PetscCall(PetscLayoutSetSize(map, N)); 213 PetscCall(PetscLayoutSetLocalSize(map, n)); 214 PetscCall(ISLoad(offIS, viewer)); 215 PetscCall(ISGetIndices(dofIS, &dofs)); 216 PetscCall(ISGetIndices(offIS, &offs)); 217 for (p = pStart, n = 0; p < pEnd; ++p, ++n) { 218 PetscCall(PetscSectionSetDof(s, p, dofs[n])); 219 PetscCall(PetscSectionSetOffset(s, p, offs[n])); 220 } 221 PetscCall(ISRestoreIndices(dofIS, &dofs)); 222 PetscCall(ISRestoreIndices(offIS, &offs)); 223 PetscCall(ISDestroy(&dofIS)); 224 PetscCall(ISDestroy(&offIS)); 225 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, NULL, (void *)&hasConstraints)); 226 if (hasConstraints) { 227 PetscCall(PetscViewerHDF5PushGroup(viewer, "bc")); 228 PetscCall(ISCreate(comm, &cdofIS)); 229 PetscCall(PetscObjectSetName((PetscObject)cdofIS, "atlasDof")); 230 PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N)); 231 #if defined(PETSC_USE_DEBUG) 232 PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->bc->atlasDof: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n); 233 #endif 234 PetscCall(ISGetLayout(cdofIS, &map)); 235 PetscCall(PetscLayoutSetSize(map, N)); 236 PetscCall(PetscLayoutSetLocalSize(map, n)); 237 PetscCall(ISLoad(cdofIS, viewer)); 238 PetscCall(ISGetIndices(cdofIS, &cdofs)); 239 for (p = pStart, n = 0; p < pEnd; ++p, ++n) PetscCall(PetscSectionSetConstraintDof(s, p, cdofs[n])); 240 PetscCall(ISRestoreIndices(cdofIS, &cdofs)); 241 PetscCall(ISDestroy(&cdofIS)); 242 PetscCall(ISCreate(comm, &coffIS)); 243 PetscCall(PetscObjectSetName((PetscObject)coffIS, "atlasOff")); 244 PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N)); 245 #if defined(PETSC_USE_DEBUG) 246 PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->bc->atlasOff: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, N1, N, n); 247 #endif 248 PetscCall(ISGetLayout(coffIS, &map)); 249 PetscCall(PetscLayoutSetSize(map, N)); 250 PetscCall(PetscLayoutSetLocalSize(map, n)); 251 PetscCall(ISLoad(coffIS, viewer)); 252 PetscCall(PetscViewerHDF5PopGroup(viewer)); 253 PetscCall(ISCreate(comm, &cindIS)); 254 PetscCall(PetscObjectSetName((PetscObject)cindIS, "bcIndices")); 255 PetscCall(PetscViewerHDF5ReadSizes(viewer, "bcIndices", NULL, &M)); 256 if (!s->bc) m = 0; 257 else PetscCall(PetscSectionGetStorageSize(s->bc, &m)); 258 #if defined(PETSC_USE_DEBUG) 259 PetscCallMPI(MPIU_Allreduce(&m, &M1, 1, MPIU_INT, MPI_SUM, comm)); 260 PetscCheck(M1 == M, comm, PETSC_ERR_ARG_SIZ, "Unable to load s->bcIndices: sum of local sizes (%" PetscInt_FMT ") != global size (%" PetscInt_FMT "): local size on this process is %" PetscInt_FMT, M1, M, m); 261 #endif 262 PetscCall(ISGetLayout(cindIS, &map)); 263 PetscCall(PetscLayoutSetSize(map, M)); 264 PetscCall(PetscLayoutSetLocalSize(map, m)); 265 PetscCall(ISLoad(cindIS, viewer)); 266 PetscCall(PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(s, cindIS, coffIS)); 267 PetscCall(ISDestroy(&coffIS)); 268 PetscCall(ISDestroy(&cindIS)); 269 } 270 PetscFunctionReturn(PETSC_SUCCESS); 271 } 272 273 PetscErrorCode PetscSectionLoad_HDF5_Internal(PetscSection s, PetscViewer viewer) 274 { 275 MPI_Comm comm; 276 PetscInt N, n, numFields, f; 277 278 PetscFunctionBegin; 279 PetscCall(PetscObjectGetComm((PetscObject)s, &comm)); 280 PetscCall(PetscViewerHDF5PushGroup(viewer, "section")); 281 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numFields", PETSC_INT, NULL, (void *)&numFields)); 282 if (s->pStart < 0 && s->pEnd < 0) n = PETSC_DECIDE; 283 else { 284 PetscCheck(s->pStart == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pStart must be 0 (got %" PetscInt_FMT ")", s->pStart); 285 PetscCheck(s->pEnd >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pEnd must be >= 0, (got %" PetscInt_FMT ")", s->pEnd); 286 n = s->pEnd; 287 } 288 if (numFields > 0) PetscCall(PetscSectionSetNumFields(s, numFields)); 289 PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N)); 290 if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnership(comm, &n, &N)); 291 PetscCall(PetscSectionSetChart(s, 0, n)); 292 PetscCall(PetscSectionLoad_HDF5_SingleField(s, viewer)); 293 for (f = 0; f < numFields; ++f) { 294 char fname[PETSC_MAX_PATH_LEN]; 295 char *fieldName; 296 PetscInt fieldComponents, c; 297 298 PetscCall(PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f)); 299 PetscCall(PetscViewerHDF5PushGroup(viewer, fname)); 300 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldName", PETSC_STRING, NULL, &fieldName)); 301 PetscCall(PetscSectionSetFieldName(s, f, fieldName)); 302 PetscCall(PetscFree(fieldName)); 303 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldComponents", PETSC_INT, NULL, (void *)&fieldComponents)); 304 PetscCall(PetscSectionSetFieldComponents(s, f, fieldComponents)); 305 for (c = 0; c < fieldComponents; ++c) { 306 char cname[PETSC_MAX_PATH_LEN]; 307 char *componentName; 308 309 PetscCall(PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c)); 310 PetscCall(PetscViewerHDF5PushGroup(viewer, cname)); 311 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "componentName", PETSC_STRING, NULL, &componentName)); 312 PetscCall(PetscSectionSetComponentName(s, f, c, componentName)); 313 PetscCall(PetscFree(componentName)); 314 PetscCall(PetscViewerHDF5PopGroup(viewer)); 315 } 316 PetscCall(PetscSectionLoad_HDF5_SingleField(s->field[f], viewer)); 317 PetscCall(PetscViewerHDF5PopGroup(viewer)); 318 } 319 PetscCall(PetscViewerHDF5PopGroup(viewer)); 320 PetscFunctionReturn(PETSC_SUCCESS); 321 } 322