1 static char help[] = "Tests save/load of plex/section/vec on different numbers of processes in HDF5.\n\n";
2
3 #include <petscdmshell.h>
4 #include <petscdmplex.h>
5 #include <petscsection.h>
6 #include <petscsf.h>
7 #include <petsclayouthdf5.h>
8
9 /* A six-element mesh
10
11 =====================
12 Save on 2 processes
13 =====================
14
15 exampleDMPlex: Local numbering:
16
17 7---17--8---18--9--19--(12)(24)(13)
18 | | | | |
19 rank 0: 20 0 21 1 22 2 (25) (3)(26)
20 | | | | |
21 4---14--5---15--6--16--(10)(23)(11)
22
23 (13)(25)--8--17---9--18--10--19--11
24 | | | | |
25 rank 1: (26) (3) 20 0 21 1 22 2 23
26 | | | | |
27 (12)(24)--4--14---5--15---6--16---7
28
29 exampleDMPlex: globalPointNumbering:
30
31 9--23--10--24--11--25--16--32--17--33--18--34--19
32 | | | | | | |
33 26 0 27 1 28 2 35 3 36 4 37 5 38
34 | | | | | | |
35 6--20---7--21---8--22--12--29--13--30--14--31--15
36
37 exampleSectionDM:
38 - includesConstraints = TRUE for local section (default)
39 - includesConstraints = FALSE for global section (default)
40
41 exampleSectionDM: Dofs (Field 0):
42
43 0---0---0---0---0---0---2---0---0---0---0---0---0
44 | | | | | | |
45 0 0 0 0 0 0 0 2 0 0 0 0 0
46 | | | | | | |
47 0---0---0---0---0---0---0---0---0---0---0---0---0
48
49 exampleSectionDM: Dofs (Field 1): constrained
50 /
51 0---0---0---0---0---0---1---0---0---0---0---0---0
52 | | | | | | |
53 0 0 0 0 0 0 2 0 0 1 0 0 0
54 | | | | | | |
55 0---0---0---0---0---0---0---0---0---0---0---0---0
56
57 exampleSectionDM: Offsets (total) in global section:
58
59 0---0---0---0---0---0---3---5---5---5---5---5---5
60 | | | | | | |
61 0 0 0 0 0 0 5 0 7 2 7 3 7
62 | | | | | | |
63 0---0---0---0---0---0---3---5---3---5---3---5---3
64
65 exampleVec: Values (Field 0): (1.3, 1.4)
66 /
67 +-------+-------+-------*-------+-------+-------+
68 | | | | | | |
69 | | | | * (1.0, 1.1)| |
70 | | | | | | |
71 +-------+-------+-------+-------+-------+-------+
72
73 exampleVec: Values (Field 1): (1.5,) constrained
74 /
75 +-------+-------+-------*-------+-------+-------+
76 | | | | | | |
77 | | (1.6, 1.7) * | * (1.2,) |
78 | | | | | | |
79 +-------+-------+-------+-------+-------+-------+
80
81 exampleVec: as global vector
82
83 rank 0: []
84 rank 1: [1.0, 1.1, 1.2, 1.3, 1.4, 1.6, 1.7]
85
86 =====================
87 Load on 3 Processes
88 =====================
89
90 exampleDMPlex: Loaded/Distributed:
91
92 5--13---6--14--(8)(18)(10)
93 | | | |
94 rank 0: 15 0 16 1 (19)(2)(20)
95 | | | |
96 3--11---4--12--(7)(17)-(9)
97
98 (9)(21)--5--15---7--18-(12)(24)(13)
99 | | | | |
100 rank 1: (22) (2) 16 0 19 1 (25) (3)(26)
101 | | | | |
102 (8)(20)--4--14---6--17-(10)(23)(11)
103
104 +-> (10)(19)--6--13---7--14---8
105 permute | | | | |
106 rank 2: +-> (20) (2) 15 0 16 1 17
107 | | | |
108 (9)(18)--3--11---4--12---5
109
110 exampleSectionDM:
111 - includesConstraints = TRUE for local section (default)
112 - includesConstraints = FALSE for global section (default)
113
114 exampleVec: as local vector:
115
116 rank 0: [1.3, 1.4, 1.5, 1.6, 1.7]
117 rank 1: [1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7]
118 rank 2: [1.2, 1.0, 1.1, 1.6, 1.7, 1.3, 1.4, 1.5]
119
120 exampleVec: as global vector:
121
122 rank 0: []
123 rank 1: [1.0, 1.1, 1.3, 1.4, 1.6, 1.7]
124 rank 2: [1.2]
125
126 */
127
128 typedef struct {
129 char fname[PETSC_MAX_PATH_LEN]; /* Output mesh filename */
130 PetscBool shell; /* Use DMShell to wrap sections */
131 } AppCtx;
132
ProcessOptions(MPI_Comm comm,AppCtx * options)133 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
134 {
135 PetscBool flg;
136
137 PetscFunctionBegin;
138 options->fname[0] = '\0';
139 PetscOptionsBegin(comm, "", "DMPlex View/Load Test Options", "DMPLEX");
140 PetscCall(PetscOptionsString("-fname", "The output mesh file", "ex12.c", options->fname, options->fname, sizeof(options->fname), &flg));
141 PetscCall(PetscOptionsBool("-shell", "Use DMShell to wrap sections", "ex12.c", options->shell, &options->shell, NULL));
142 PetscOptionsEnd();
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
main(int argc,char ** argv)146 int main(int argc, char **argv)
147 {
148 MPI_Comm comm;
149 PetscMPIInt size, rank, mycolor;
150 const char exampleDMPlexName[] = "exampleDMPlex";
151 const char exampleSectionDMName[] = "exampleSectionDM";
152 const char exampleVecName[] = "exampleVec";
153 PetscScalar constraintValue = 1.5;
154 PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
155 AppCtx user;
156
157 PetscFunctionBeginUser;
158 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
159 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
160 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
161 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
162 PetscCheck(size >= 3, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes");
163
164 /* Save */
165 mycolor = (PetscMPIInt)(rank >= 2);
166 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
167 if (mycolor == 0) {
168 DM dm;
169 PetscViewer viewer;
170
171 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer));
172 /* Save exampleDMPlex */
173 {
174 DM pdm;
175 const PetscInt faces[2] = {6, 1};
176 PetscSF sf;
177 PetscInt overlap = 1;
178
179 PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &dm));
180 PetscCall(DMPlexDistribute(dm, overlap, &sf, &pdm));
181 if (pdm) {
182 PetscCall(DMDestroy(&dm));
183 dm = pdm;
184 }
185 PetscCall(PetscSFDestroy(&sf));
186 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
187 PetscCall(PetscViewerPushFormat(viewer, format));
188 PetscCall(DMPlexTopologyView(dm, viewer));
189 PetscCall(DMPlexLabelsView(dm, viewer));
190 PetscCall(PetscViewerPopFormat(viewer));
191 }
192 /* Save coordinates */
193 PetscCall(PetscViewerPushFormat(viewer, format));
194 PetscCall(DMPlexCoordinatesView(dm, viewer));
195 PetscCall(PetscViewerPopFormat(viewer));
196 /* Save exampleVec */
197 {
198 PetscInt pStart = -1, pEnd = -1;
199 DM sdm;
200 PetscSection section, gsection;
201 PetscBool includesConstraints = PETSC_FALSE;
202 Vec vec;
203 PetscScalar *array = NULL;
204
205 /* Create section */
206 PetscCall(PetscSectionCreate(comm, §ion));
207 PetscCall(PetscSectionSetNumFields(section, 2));
208 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
209 PetscCall(PetscSectionSetChart(section, pStart, pEnd));
210 switch (rank) {
211 case 0:
212 PetscCall(PetscSectionSetDof(section, 3, 2));
213 PetscCall(PetscSectionSetDof(section, 12, 3));
214 PetscCall(PetscSectionSetDof(section, 25, 2));
215 PetscCall(PetscSectionSetConstraintDof(section, 12, 1));
216 PetscCall(PetscSectionSetFieldDof(section, 3, 0, 2));
217 PetscCall(PetscSectionSetFieldDof(section, 12, 0, 2));
218 PetscCall(PetscSectionSetFieldDof(section, 12, 1, 1));
219 PetscCall(PetscSectionSetFieldDof(section, 25, 1, 2));
220 PetscCall(PetscSectionSetFieldConstraintDof(section, 12, 1, 1));
221 break;
222 case 1:
223 PetscCall(PetscSectionSetDof(section, 0, 2));
224 PetscCall(PetscSectionSetDof(section, 1, 1));
225 PetscCall(PetscSectionSetDof(section, 8, 3));
226 PetscCall(PetscSectionSetDof(section, 20, 2));
227 PetscCall(PetscSectionSetConstraintDof(section, 8, 1));
228 PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2));
229 PetscCall(PetscSectionSetFieldDof(section, 8, 0, 2));
230 PetscCall(PetscSectionSetFieldDof(section, 1, 1, 1));
231 PetscCall(PetscSectionSetFieldDof(section, 8, 1, 1));
232 PetscCall(PetscSectionSetFieldDof(section, 20, 1, 2));
233 PetscCall(PetscSectionSetFieldConstraintDof(section, 8, 1, 1));
234 break;
235 }
236 PetscCall(PetscSectionSetUp(section));
237 {
238 const PetscInt indices[] = {2};
239 const PetscInt indices1[] = {0};
240
241 switch (rank) {
242 case 0:
243 PetscCall(PetscSectionSetConstraintIndices(section, 12, indices));
244 PetscCall(PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1));
245 break;
246 case 1:
247 PetscCall(PetscSectionSetConstraintIndices(section, 8, indices));
248 PetscCall(PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1));
249 break;
250 }
251 }
252 if (user.shell) {
253 PetscSF sf;
254
255 PetscCall(DMShellCreate(comm, &sdm));
256 PetscCall(DMGetPointSF(dm, &sf));
257 PetscCall(DMSetPointSF(sdm, sf));
258 } else {
259 PetscCall(DMClone(dm, &sdm));
260 }
261 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName));
262 PetscCall(DMSetLocalSection(sdm, section));
263 PetscCall(PetscSectionDestroy(§ion));
264 PetscCall(DMPlexSectionView(dm, viewer, sdm));
265 /* Create global vector */
266 PetscCall(DMGetGlobalSection(sdm, &gsection));
267 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints));
268 if (user.shell) {
269 PetscInt n = -1;
270
271 PetscCall(VecCreate(comm, &vec));
272 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &n));
273 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &n));
274 PetscCall(VecSetSizes(vec, n, PETSC_DECIDE));
275 PetscCall(VecSetUp(vec));
276 } else {
277 PetscCall(DMGetGlobalVector(sdm, &vec));
278 }
279 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
280 PetscCall(VecGetArrayWrite(vec, &array));
281 if (includesConstraints) {
282 switch (rank) {
283 case 0:
284 break;
285 case 1:
286 array[0] = 1.0;
287 array[1] = 1.1;
288 array[2] = 1.2;
289 array[3] = 1.3;
290 array[4] = 1.4;
291 array[5] = 1.5;
292 array[6] = 1.6;
293 array[7] = 1.7;
294 break;
295 }
296 } else {
297 switch (rank) {
298 case 0:
299 break;
300 case 1:
301 array[0] = 1.0;
302 array[1] = 1.1;
303 array[2] = 1.2;
304 array[3] = 1.3;
305 array[4] = 1.4;
306 array[5] = 1.6;
307 array[6] = 1.7;
308 break;
309 }
310 }
311 PetscCall(VecRestoreArrayWrite(vec, &array));
312 PetscCall(DMPlexGlobalVectorView(dm, viewer, sdm, vec));
313 if (user.shell) {
314 PetscCall(VecDestroy(&vec));
315 } else {
316 PetscCall(DMRestoreGlobalVector(sdm, &vec));
317 }
318 PetscCall(DMDestroy(&sdm));
319 }
320 PetscCall(PetscViewerDestroy(&viewer));
321 PetscCall(DMDestroy(&dm));
322 }
323 PetscCallMPI(MPI_Comm_free(&comm));
324 /* Load */
325 mycolor = (PetscMPIInt)(rank >= 3);
326 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
327 if (mycolor == 0) {
328 DM dm;
329 PetscSF sfXC;
330 PetscViewer viewer;
331
332 PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer));
333 /* Load exampleDMPlex */
334 {
335 PetscSF sfXB, sfBC;
336
337 PetscCall(DMCreate(comm, &dm));
338 PetscCall(DMSetType(dm, DMPLEX));
339 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
340 /* sfXB: X -> B */
341 /* X: set of globalPointNumbers, [0, N) */
342 /* B: loaded naive in-memory plex */
343 PetscCall(PetscViewerPushFormat(viewer, format));
344 PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB));
345 PetscCall(PetscViewerPopFormat(viewer));
346 {
347 DM distributedDM;
348 PetscInt overlap = 1;
349 PetscPartitioner part;
350
351 PetscCall(DMPlexGetPartitioner(dm, &part));
352 PetscCall(PetscPartitionerSetFromOptions(part));
353 /* sfBC: B -> C */
354 /* B: loaded naive in-memory plex */
355 /* C: redistributed good in-memory */
356 PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM));
357 if (distributedDM) {
358 PetscCall(DMDestroy(&dm));
359 dm = distributedDM;
360 }
361 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
362 }
363 /* sfXC: X -> C */
364 PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
365 PetscCall(PetscSFDestroy(&sfXB));
366 PetscCall(PetscSFDestroy(&sfBC));
367 }
368 /* Load labels */
369 PetscCall(PetscViewerPushFormat(viewer, format));
370 PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC));
371 PetscCall(PetscViewerPopFormat(viewer));
372 /* Load coordinates */
373 PetscCall(PetscViewerPushFormat(viewer, format));
374 PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC));
375 PetscCall(PetscViewerPopFormat(viewer));
376 PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)"));
377 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
378 PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
379 /* Load exampleVec */
380 {
381 DM sdm;
382 PetscSection section, gsection;
383 IS perm;
384 PetscBool includesConstraints = PETSC_FALSE;
385 Vec vec;
386 PetscSF lsf, gsf;
387
388 if (user.shell) {
389 PetscSF sf;
390
391 PetscCall(DMShellCreate(comm, &sdm));
392 PetscCall(DMGetPointSF(dm, &sf));
393 PetscCall(DMSetPointSF(sdm, sf));
394 } else {
395 PetscCall(DMClone(dm, &sdm));
396 }
397 PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName));
398 PetscCall(PetscSectionCreate(comm, §ion));
399 {
400 PetscInt pStart = -1, pEnd = -1, p = -1;
401 PetscInt *pinds = NULL;
402
403 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
404 PetscCall(PetscMalloc1(pEnd - pStart, &pinds));
405 for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p;
406 if (rank == 2) {
407 pinds[10] = 20;
408 pinds[20] = 10;
409 }
410 PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm));
411 }
412 PetscCall(PetscSectionSetPermutation(section, perm));
413 PetscCall(ISDestroy(&perm));
414 PetscCall(DMSetLocalSection(sdm, section));
415 PetscCall(PetscSectionDestroy(§ion));
416 PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf));
417 /* Load as local vector */
418 PetscCall(DMGetLocalSection(sdm, §ion));
419 PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section"));
420 PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
421 PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
422 if (user.shell) {
423 PetscInt m = -1;
424
425 PetscCall(VecCreate(comm, &vec));
426 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
427 else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
428 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
429 PetscCall(VecSetUp(vec));
430 } else {
431 PetscCall(DMGetLocalVector(sdm, &vec));
432 }
433 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
434 PetscCall(VecSet(vec, constraintValue));
435 PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec));
436 PetscCall(PetscSFDestroy(&lsf));
437 if (user.shell) {
438 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
439 PetscCall(VecDestroy(&vec));
440 } else {
441 PetscCall(DMRestoreLocalVector(sdm, &vec));
442 }
443 /* Load as global vector */
444 PetscCall(DMGetGlobalSection(sdm, &gsection));
445 PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section"));
446 PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm)));
447 PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints));
448 if (user.shell) {
449 PetscInt m = -1;
450
451 PetscCall(VecCreate(comm, &vec));
452 if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m));
453 else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m));
454 PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
455 PetscCall(VecSetUp(vec));
456 } else {
457 PetscCall(DMGetGlobalVector(sdm, &vec));
458 }
459 PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
460 PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec));
461 PetscCall(PetscSFDestroy(&gsf));
462 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
463 if (user.shell) {
464 PetscCall(VecDestroy(&vec));
465 } else {
466 PetscCall(DMRestoreGlobalVector(sdm, &vec));
467 }
468 PetscCall(DMDestroy(&sdm));
469 }
470 PetscCall(PetscViewerDestroy(&viewer));
471 PetscCall(PetscSFDestroy(&sfXC));
472 PetscCall(DMDestroy(&dm));
473 }
474 PetscCallMPI(MPI_Comm_free(&comm));
475
476 /* Finalize */
477 PetscCall(PetscFinalize());
478 return 0;
479 }
480
481 /*TEST
482
483 build:
484 requires: hdf5
485 testset:
486 suffix: 0
487 requires: !complex
488 nsize: 4
489 args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail
490 args: -dm_plex_view_hdf5_storage_version 2.0.0
491 test:
492 suffix: parmetis
493 requires: parmetis
494 args: -petscpartitioner_type parmetis
495 test:
496 suffix: ptscotch
497 requires: ptscotch
498 args: -petscpartitioner_type ptscotch
499
500 TEST*/
501