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