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: 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;
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;
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: }
89: PetscErrorCode PetscSectionView_HDF5_Internal(PetscSection s, PetscViewer viewer)
90: {
91: PetscInt numFields, f;
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;
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;
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: }
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;
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: }
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;
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: }
273: PetscErrorCode PetscSectionLoad_HDF5_Internal(PetscSection s, PetscViewer viewer)
274: {
275: MPI_Comm comm;
276: PetscInt N, n, numFields, f;
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;
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;
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: }