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
PetscSectionView_HDF5_SingleField(PetscSection s,PetscViewer viewer)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, MPI_C_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
PetscSectionView_HDF5_Internal(PetscSection s,PetscViewer viewer)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
PetscSectionLoad_HDF5_SingleField_SetConstraintIndices(PetscSection s,IS cindIS,IS coffIS)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
PetscSectionLoad_HDF5_SingleField(PetscSection s,PetscViewer viewer)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
PetscSectionLoad_HDF5_Internal(PetscSection s,PetscViewer viewer)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