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;
 15:   PetscErrorCode  ierr;

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

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

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

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

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

116:       PetscSNPrintf(cname, sizeof(cname), "component%D", c);
117:       PetscViewerHDF5PushGroup(viewer, cname);
118:       PetscSectionGetComponentName(s, f, c, &componentName);
119:       PetscViewerHDF5WriteAttribute(viewer, NULL, "componentName", PETSC_STRING, componentName);
120:       PetscViewerHDF5PopGroup(viewer);
121:     }
122:     PetscSectionView_HDF5_SingleField(s->field[f], viewer);
123:     PetscViewerHDF5PopGroup(viewer);
124:   }
125:   PetscViewerHDF5PopGroup(viewer);
126:   return(0);
127: }

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

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

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

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

280: PetscErrorCode PetscSectionLoad_HDF5_Internal(PetscSection s, PetscViewer viewer)
281: {
282:   MPI_Comm        comm;
283:   PetscInt        N, n, numFields, f;
284:   PetscErrorCode  ierr;

287:   PetscObjectGetComm((PetscObject)s, &comm);
288:   PetscViewerHDF5PushGroup(viewer, "section");
289:   PetscViewerHDF5ReadAttribute(viewer, NULL, "numFields", PETSC_INT, NULL, (void *)&numFields);
290:   if (s->pStart < 0 && s->pEnd < 0) n = PETSC_DECIDE;
291:   else {
292:     if (s->pStart != 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pStart must be 0 (got %D)", s->pStart);
293:     if (s->pEnd < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "s->pEnd must be >= 0, (got %D)", s->pEnd);
294:     n = s->pEnd;
295:   }
296:   if (numFields > 0) {PetscSectionSetNumFields(s, numFields);}
297:   PetscViewerHDF5ReadSizes(viewer, "atlasDof", NULL, &N);
298:   if (n == PETSC_DECIDE) {PetscSplitOwnership(comm, &n, &N);}
299:   PetscSectionSetChart(s, 0, n);
300:   PetscSectionLoad_HDF5_SingleField(s, viewer);
301:   for (f = 0; f < numFields; ++f) {
302:     char      fname[PETSC_MAX_PATH_LEN];
303:     char     *fieldName;
304:     PetscInt  fieldComponents, c;

306:     PetscSNPrintf(fname, sizeof(fname), "field%D", f);
307:     PetscViewerHDF5PushGroup(viewer, fname);
308:     PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldName", PETSC_STRING, NULL, &fieldName);
309:     PetscSectionSetFieldName(s, f, fieldName);
310:     PetscFree(fieldName);
311:     PetscViewerHDF5ReadAttribute(viewer, NULL, "fieldComponents", PETSC_INT, NULL, (void *) &fieldComponents);
312:     PetscSectionSetFieldComponents(s, f, fieldComponents);
313:     for (c = 0; c < fieldComponents; ++c) {
314:       char  cname[PETSC_MAX_PATH_LEN];
315:       char *componentName;

317:       PetscSNPrintf(cname, sizeof(cname), "component%D", c);
318:       PetscViewerHDF5PushGroup(viewer, cname);
319:       PetscViewerHDF5ReadAttribute(viewer, NULL, "componentName", PETSC_STRING, NULL, &componentName);
320:       PetscSectionSetComponentName(s, f, c, componentName);
321:       PetscFree(componentName);
322:       PetscViewerHDF5PopGroup(viewer);
323:     }
324:     PetscSectionLoad_HDF5_SingleField(s->field[f], viewer);
325:     PetscViewerHDF5PopGroup(viewer);
326:   }
327:   PetscViewerHDF5PopGroup(viewer);
328:   return(0);
329: }

331: #endif