Actual source code: sectionhdf5.c

  1: #include <petsc/private/sectionimpl.h>
  2: #include <petscsf.h>
  3: #include <petscis.h>
  4: #include <petscviewerhdf5.h>
  5: #include <petsclayouthdf5.h>

  7: #if defined(PETSC_HAVE_HDF5)
  8: static PetscErrorCode PetscSectionView_HDF5_SingleField(PetscSection s, PetscViewer viewer)
  9: {
 10:   MPI_Comm  comm;
 11:   PetscInt  pStart, pEnd, p, n;
 12:   PetscBool hasConstraints, includesConstraints;
 13:   IS        dofIS, offIS, cdofIS, coffIS, cindIS;
 14:   PetscInt *dofs, *offs, *cdofs, *coffs, *cinds, dof, cdof, m, moff, i;

 16:   PetscFunctionBegin;
 17:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
 18:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
 19:   hasConstraints = (s->bc) ? PETSC_TRUE : PETSC_FALSE;
 20:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &hasConstraints, 1, MPIU_BOOL, MPI_LOR, comm));
 21:   for (p = pStart, n = 0, m = 0; p < pEnd; ++p) {
 22:     PetscCall(PetscSectionGetDof(s, p, &dof));
 23:     if (dof >= 0) {
 24:       if (hasConstraints) {
 25:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
 26:         m += cdof;
 27:       }
 28:       n++;
 29:     }
 30:   }
 31:   PetscCall(PetscMalloc1(n, &dofs));
 32:   PetscCall(PetscMalloc1(n, &offs));
 33:   if (hasConstraints) {
 34:     PetscCall(PetscMalloc1(n, &cdofs));
 35:     PetscCall(PetscMalloc1(n, &coffs));
 36:     PetscCall(PetscMalloc1(m, &cinds));
 37:   }
 38:   for (p = pStart, n = 0, m = 0; p < pEnd; ++p) {
 39:     PetscCall(PetscSectionGetDof(s, p, &dof));
 40:     if (dof >= 0) {
 41:       dofs[n] = dof;
 42:       PetscCall(PetscSectionGetOffset(s, p, &offs[n]));
 43:       if (hasConstraints) {
 44:         const PetscInt *cpinds;

 46:         PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
 47:         PetscCall(PetscSectionGetConstraintIndices(s, p, &cpinds));
 48:         cdofs[n] = cdof;
 49:         coffs[n] = m;
 50:         for (i = 0; i < cdof; ++i) cinds[m++] = cpinds[i];
 51:       }
 52:       n++;
 53:     }
 54:   }
 55:   if (hasConstraints) {
 56:     PetscCallMPI(MPI_Scan(&m, &moff, 1, MPIU_INT, MPI_SUM, comm));
 57:     moff -= m;
 58:     for (p = 0; p < n; ++p) coffs[p] += moff;
 59:   }
 60:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, (void *)&hasConstraints));
 61:   PetscCall(PetscSectionGetIncludesConstraints(s, &includesConstraints));
 62:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, (void *)&includesConstraints));
 63:   PetscCall(ISCreateGeneral(comm, n, dofs, PETSC_OWN_POINTER, &dofIS));
 64:   PetscCall(PetscObjectSetName((PetscObject)dofIS, "atlasDof"));
 65:   PetscCall(ISView(dofIS, viewer));
 66:   PetscCall(ISDestroy(&dofIS));
 67:   PetscCall(ISCreateGeneral(comm, n, offs, PETSC_OWN_POINTER, &offIS));
 68:   PetscCall(PetscObjectSetName((PetscObject)offIS, "atlasOff"));
 69:   PetscCall(ISView(offIS, viewer));
 70:   PetscCall(ISDestroy(&offIS));
 71:   if (hasConstraints) {
 72:     PetscCall(PetscViewerHDF5PushGroup(viewer, "bc"));
 73:     PetscCall(ISCreateGeneral(comm, n, cdofs, PETSC_OWN_POINTER, &cdofIS));
 74:     PetscCall(PetscObjectSetName((PetscObject)cdofIS, "atlasDof"));
 75:     PetscCall(ISView(cdofIS, viewer));
 76:     PetscCall(ISDestroy(&cdofIS));
 77:     PetscCall(ISCreateGeneral(comm, n, coffs, PETSC_OWN_POINTER, &coffIS));
 78:     PetscCall(PetscObjectSetName((PetscObject)coffIS, "atlasOff"));
 79:     PetscCall(ISView(coffIS, viewer));
 80:     PetscCall(ISDestroy(&coffIS));
 81:     PetscCall(PetscViewerHDF5PopGroup(viewer));
 82:     PetscCall(ISCreateGeneral(comm, m, cinds, PETSC_OWN_POINTER, &cindIS));
 83:     PetscCall(PetscObjectSetName((PetscObject)cindIS, "bcIndices"));
 84:     PetscCall(ISView(cindIS, viewer));
 85:     PetscCall(ISDestroy(&cindIS));
 86:   }
 87:   PetscFunctionReturn(PETSC_SUCCESS);
 88: }

 90: PetscErrorCode PetscSectionView_HDF5_Internal(PetscSection s, PetscViewer viewer)
 91: {
 92:   PetscInt numFields, f;

 94:   PetscFunctionBegin;
 95:   PetscCall(PetscViewerHDF5PushGroup(viewer, "section"));
 96:   PetscCall(PetscSectionGetNumFields(s, &numFields));
 97:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numFields", PETSC_INT, (void *)&numFields));
 98:   PetscCall(PetscSectionView_HDF5_SingleField(s, viewer));
 99:   for (f = 0; f < numFields; ++f) {
100:     char        fname[PETSC_MAX_PATH_LEN];
101:     const char *fieldName;
102:     PetscInt    fieldComponents, c;

104:     PetscCall(PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f));
105:     PetscCall(PetscViewerHDF5PushGroup(viewer, fname));
106:     PetscCall(PetscSectionGetFieldName(s, f, &fieldName));
107:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldName", PETSC_STRING, fieldName));
108:     PetscCall(PetscSectionGetFieldComponents(s, f, &fieldComponents));
109:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "fieldComponents", PETSC_INT, (void *)&fieldComponents));
110:     for (c = 0; c < fieldComponents; ++c) {
111:       char        cname[PETSC_MAX_PATH_LEN];
112:       const char *componentName;

114:       PetscCall(PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c));
115:       PetscCall(PetscViewerHDF5PushGroup(viewer, cname));
116:       PetscCall(PetscSectionGetComponentName(s, f, c, &componentName));
117:       PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "componentName", PETSC_STRING, componentName));
118:       PetscCall(PetscViewerHDF5PopGroup(viewer));
119:     }
120:     PetscCall(PetscSectionView_HDF5_SingleField(s->field[f], viewer));
121:     PetscCall(PetscViewerHDF5PopGroup(viewer));
122:   }
123:   PetscCall(PetscViewerHDF5PopGroup(viewer));
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: static PetscErrorCode PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(PetscSection s, IS cindIS, IS coffIS)
128: {
129:   MPI_Comm        comm;
130:   PetscInt        pStart, pEnd, p, M, m, i, cdof;
131:   const PetscInt *data;
132:   PetscInt       *cinds;
133:   const PetscInt *coffs;
134:   PetscInt       *coffsets;
135:   PetscSF         sf;
136:   PetscLayout     layout;

138:   PetscFunctionBegin;
139:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
140:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
141:   PetscCall(ISGetSize(cindIS, &M));
142:   PetscCall(ISGetLocalSize(cindIS, &m));
143:   PetscCall(PetscMalloc1(m, &coffsets));
144:   PetscCall(ISGetIndices(coffIS, &coffs));
145:   for (p = pStart, m = 0; p < pEnd; ++p) {
146:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
147:     for (i = 0; i < cdof; ++i) coffsets[m++] = coffs[p - pStart] + i;
148:   }
149:   PetscCall(ISRestoreIndices(coffIS, &coffs));
150:   PetscCall(PetscSFCreate(comm, &sf));
151:   PetscCall(PetscLayoutCreate(comm, &layout));
152:   PetscCall(PetscLayoutSetSize(layout, M));
153:   PetscCall(PetscLayoutSetLocalSize(layout, m));
154:   PetscCall(PetscLayoutSetBlockSize(layout, 1));
155:   PetscCall(PetscLayoutSetUp(layout));
156:   PetscCall(PetscSFSetGraphLayout(sf, layout, m, NULL, PETSC_OWN_POINTER, coffsets));
157:   PetscCall(PetscLayoutDestroy(&layout));
158:   PetscCall(PetscFree(coffsets));
159:   PetscCall(PetscMalloc1(m, &cinds));
160:   PetscCall(ISGetIndices(cindIS, &data));
161:   PetscCall(PetscSFBcastBegin(sf, MPIU_INT, data, cinds, MPI_REPLACE));
162:   PetscCall(PetscSFBcastEnd(sf, MPIU_INT, data, cinds, MPI_REPLACE));
163:   PetscCall(ISRestoreIndices(cindIS, &data));
164:   PetscCall(PetscSFDestroy(&sf));
165:   PetscCall(PetscSectionSetUpBC(s));
166:   for (p = pStart, m = 0; p < pEnd; ++p) {
167:     PetscCall(PetscSectionGetConstraintDof(s, p, &cdof));
168:     PetscCall(PetscSectionSetConstraintIndices(s, p, &cinds[m]));
169:     m += cdof;
170:   }
171:   PetscCall(PetscFree(cinds));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: static PetscErrorCode PetscSectionLoad_HDF5_SingleField(PetscSection s, PetscViewer viewer)
176: {
177:   MPI_Comm comm;
178:   PetscInt pStart, pEnd, p, N, n, M, m;
179:   #if defined(PETSC_USE_DEBUG)
180:   PetscInt N1, M1;
181:   #endif
182:   PetscBool       hasConstraints, includesConstraints;
183:   IS              dofIS, offIS, cdofIS, coffIS, cindIS;
184:   const PetscInt *dofs, *offs, *cdofs;
185:   PetscLayout     map;

187:   PetscFunctionBegin;
188:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
189:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "includesConstraints", PETSC_BOOL, NULL, (void *)&includesConstraints));
190:   PetscCall(PetscSectionSetIncludesConstraints(s, includesConstraints));
191:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
192:   n = pEnd - pStart;
193:   #if defined(PETSC_USE_DEBUG)
194:   PetscCall(MPIU_Allreduce(&n, &N1, 1, MPIU_INT, MPI_SUM, comm));
195:   #endif
196:   PetscCall(ISCreate(comm, &dofIS));
197:   PetscCall(PetscObjectSetName((PetscObject)dofIS, "atlasDof"));
198:   PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N));
199:   #if defined(PETSC_USE_DEBUG)
200:   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);
201:   #endif
202:   PetscCall(ISGetLayout(dofIS, &map));
203:   PetscCall(PetscLayoutSetSize(map, N));
204:   PetscCall(PetscLayoutSetLocalSize(map, n));
205:   PetscCall(ISLoad(dofIS, viewer));
206:   PetscCall(ISCreate(comm, &offIS));
207:   PetscCall(PetscObjectSetName((PetscObject)offIS, "atlasOff"));
208:   PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N));
209:   #if defined(PETSC_USE_DEBUG)
210:   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);
211:   #endif
212:   PetscCall(ISGetLayout(offIS, &map));
213:   PetscCall(PetscLayoutSetSize(map, N));
214:   PetscCall(PetscLayoutSetLocalSize(map, n));
215:   PetscCall(ISLoad(offIS, viewer));
216:   PetscCall(ISGetIndices(dofIS, &dofs));
217:   PetscCall(ISGetIndices(offIS, &offs));
218:   for (p = pStart, n = 0; p < pEnd; ++p, ++n) {
219:     PetscCall(PetscSectionSetDof(s, p, dofs[n]));
220:     PetscCall(PetscSectionSetOffset(s, p, offs[n]));
221:   }
222:   PetscCall(ISRestoreIndices(dofIS, &dofs));
223:   PetscCall(ISRestoreIndices(offIS, &offs));
224:   PetscCall(ISDestroy(&dofIS));
225:   PetscCall(ISDestroy(&offIS));
226:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "hasConstraints", PETSC_BOOL, NULL, (void *)&hasConstraints));
227:   if (hasConstraints) {
228:     PetscCall(PetscViewerHDF5PushGroup(viewer, "bc"));
229:     PetscCall(ISCreate(comm, &cdofIS));
230:     PetscCall(PetscObjectSetName((PetscObject)cdofIS, "atlasDof"));
231:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N));
232:   #if defined(PETSC_USE_DEBUG)
233:     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);
234:   #endif
235:     PetscCall(ISGetLayout(cdofIS, &map));
236:     PetscCall(PetscLayoutSetSize(map, N));
237:     PetscCall(PetscLayoutSetLocalSize(map, n));
238:     PetscCall(ISLoad(cdofIS, viewer));
239:     PetscCall(ISGetIndices(cdofIS, &cdofs));
240:     for (p = pStart, n = 0; p < pEnd; ++p, ++n) PetscCall(PetscSectionSetConstraintDof(s, p, cdofs[n]));
241:     PetscCall(ISRestoreIndices(cdofIS, &cdofs));
242:     PetscCall(ISDestroy(&cdofIS));
243:     PetscCall(ISCreate(comm, &coffIS));
244:     PetscCall(PetscObjectSetName((PetscObject)coffIS, "atlasOff"));
245:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasOff", NULL, &N));
246:   #if defined(PETSC_USE_DEBUG)
247:     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);
248:   #endif
249:     PetscCall(ISGetLayout(coffIS, &map));
250:     PetscCall(PetscLayoutSetSize(map, N));
251:     PetscCall(PetscLayoutSetLocalSize(map, n));
252:     PetscCall(ISLoad(coffIS, viewer));
253:     PetscCall(PetscViewerHDF5PopGroup(viewer));
254:     PetscCall(ISCreate(comm, &cindIS));
255:     PetscCall(PetscObjectSetName((PetscObject)cindIS, "bcIndices"));
256:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "bcIndices", NULL, &M));
257:     if (!s->bc) m = 0;
258:     else PetscCall(PetscSectionGetStorageSize(s->bc, &m));
259:   #if defined(PETSC_USE_DEBUG)
260:     PetscCall(MPIU_Allreduce(&m, &M1, 1, MPIU_INT, MPI_SUM, comm));
261:     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);
262:   #endif
263:     PetscCall(ISGetLayout(cindIS, &map));
264:     PetscCall(PetscLayoutSetSize(map, M));
265:     PetscCall(PetscLayoutSetLocalSize(map, m));
266:     PetscCall(ISLoad(cindIS, viewer));
267:     PetscCall(PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(s, cindIS, coffIS));
268:     PetscCall(ISDestroy(&coffIS));
269:     PetscCall(ISDestroy(&cindIS));
270:   }
271:   PetscFunctionReturn(PETSC_SUCCESS);
272: }

274: PetscErrorCode PetscSectionLoad_HDF5_Internal(PetscSection s, PetscViewer viewer)
275: {
276:   MPI_Comm comm;
277:   PetscInt N, n, numFields, f;

279:   PetscFunctionBegin;
280:   PetscCall(PetscObjectGetComm((PetscObject)s, &comm));
281:   PetscCall(PetscViewerHDF5PushGroup(viewer, "section"));
282:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numFields", PETSC_INT, NULL, (void *)&numFields));
283:   if (s->pStart < 0 && s->pEnd < 0) n = PETSC_DECIDE;
284:   else {
285:     PetscCheck(s->pStart == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pStart must be 0 (got %" PetscInt_FMT ")", s->pStart);
286:     PetscCheck(s->pEnd >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pEnd must be >= 0, (got %" PetscInt_FMT ")", s->pEnd);
287:     n = s->pEnd;
288:   }
289:   if (numFields > 0) PetscCall(PetscSectionSetNumFields(s, numFields));
290:   PetscCall(PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N));
291:   if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnership(comm, &n, &N));
292:   PetscCall(PetscSectionSetChart(s, 0, n));
293:   PetscCall(PetscSectionLoad_HDF5_SingleField(s, viewer));
294:   for (f = 0; f < numFields; ++f) {
295:     char     fname[PETSC_MAX_PATH_LEN];
296:     char    *fieldName;
297:     PetscInt fieldComponents, c;

299:     PetscCall(PetscSNPrintf(fname, sizeof(fname), "field%" PetscInt_FMT, f));
300:     PetscCall(PetscViewerHDF5PushGroup(viewer, fname));
301:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldName", PETSC_STRING, NULL, &fieldName));
302:     PetscCall(PetscSectionSetFieldName(s, f, fieldName));
303:     PetscCall(PetscFree(fieldName));
304:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldComponents", PETSC_INT, NULL, (void *)&fieldComponents));
305:     PetscCall(PetscSectionSetFieldComponents(s, f, fieldComponents));
306:     for (c = 0; c < fieldComponents; ++c) {
307:       char  cname[PETSC_MAX_PATH_LEN];
308:       char *componentName;

310:       PetscCall(PetscSNPrintf(cname, sizeof(cname), "component%" PetscInt_FMT, c));
311:       PetscCall(PetscViewerHDF5PushGroup(viewer, cname));
312:       PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "componentName", PETSC_STRING, NULL, &componentName));
313:       PetscCall(PetscSectionSetComponentName(s, f, c, componentName));
314:       PetscCall(PetscFree(componentName));
315:       PetscCall(PetscViewerHDF5PopGroup(viewer));
316:     }
317:     PetscCall(PetscSectionLoad_HDF5_SingleField(s->field[f], viewer));
318:     PetscCall(PetscViewerHDF5PopGroup(viewer));
319:   }
320:   PetscCall(PetscViewerHDF5PopGroup(viewer));
321:   PetscFunctionReturn(PETSC_SUCCESS);
322: }

324: #endif