Actual source code: ex21.c
1: static char help[] = "Tests save/load of plex/section/vec on different numbers of processes in HDF5.\n\n";
3: #include <petscdmshell.h>
4: #include <petscdmplex.h>
5: #include <petscsection.h>
6: #include <petscsf.h>
7: #include <petsclayouthdf5.h>
9: /* A six-element mesh
11: =====================
12: Save on 2 processes
13: =====================
15: exampleDMPlex: Local numbering:
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)
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
29: exampleDMPlex: globalPointNumbering:
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
37: exampleSectionDM:
38: - includesConstraints = TRUE for local section (default)
39: - includesConstraints = FALSE for global section (default)
41: exampleSectionDM: Dofs (Field 0):
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
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
57: exampleSectionDM: Offsets (total) in global section:
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
65: exampleVec: Values (Field 0): (1.3, 1.4)
66: /
67: +-------+-------+-------*-------+-------+-------+
68: | | | | | | |
69: | | | | * (1.0, 1.1)| |
70: | | | | | | |
71: +-------+-------+-------+-------+-------+-------+
73: exampleVec: Values (Field 1): (1.5,) constrained
74: /
75: +-------+-------+-------*-------+-------+-------+
76: | | | | | | |
77: | | (1.6, 1.7) * | * (1.2,) |
78: | | | | | | |
79: +-------+-------+-------+-------+-------+-------+
81: exampleVec: as global vector
83: rank 0: []
84: rank 1: [1.0, 1.1, 1.2, 1.3, 1.4, 1.6, 1.7]
86: =====================
87: Load on 3 Processes
88: =====================
90: exampleDMPlex: Loaded/Distributed:
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)
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)
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
110: exampleSectionDM:
111: - includesConstraints = TRUE for local section (default)
112: - includesConstraints = FALSE for global section (default)
114: exampleVec: as local vector:
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]
120: exampleVec: as global vector:
122: rank 0: []
123: rank 1: [1.0, 1.1, 1.3, 1.4, 1.6, 1.7]
124: rank 2: [1.2]
126: */
128: typedef struct {
129: char fname[PETSC_MAX_PATH_LEN]; /* Output mesh filename */
130: PetscBool shell; /* Use DMShell to wrap sections */
131: } AppCtx;
133: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
134: {
135: PetscBool flg;
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: }
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;
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");
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;
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: PetscPartitioner part;
180: PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &dm));
181: PetscCall(DMPlexGetPartitioner(dm, &part));
182: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, "save_"));
183: PetscCall(PetscPartitionerSetFromOptions(part));
184: PetscCall(DMPlexDistribute(dm, overlap, &sf, &pdm));
185: if (pdm) {
186: PetscCall(DMDestroy(&dm));
187: dm = pdm;
188: }
189: PetscCall(PetscSFDestroy(&sf));
190: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
191: PetscCall(PetscViewerPushFormat(viewer, format));
192: PetscCall(DMPlexTopologyView(dm, viewer));
193: PetscCall(DMPlexLabelsView(dm, viewer));
194: PetscCall(PetscViewerPopFormat(viewer));
195: }
196: /* Save coordinates */
197: PetscCall(PetscViewerPushFormat(viewer, format));
198: PetscCall(DMPlexCoordinatesView(dm, viewer));
199: PetscCall(PetscViewerPopFormat(viewer));
200: /* Save exampleVec */
201: {
202: PetscInt pStart = -1, pEnd = -1;
203: DM sdm;
204: PetscSection section, gsection;
205: PetscBool includesConstraints = PETSC_FALSE;
206: Vec vec;
207: PetscScalar *array = NULL;
209: /* Create section */
210: PetscCall(PetscSectionCreate(comm, §ion));
211: PetscCall(PetscSectionSetNumFields(section, 2));
212: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
213: PetscCall(PetscSectionSetChart(section, pStart, pEnd));
214: switch (rank) {
215: case 0:
216: PetscCall(PetscSectionSetDof(section, 3, 2));
217: PetscCall(PetscSectionSetDof(section, 12, 3));
218: PetscCall(PetscSectionSetDof(section, 25, 2));
219: PetscCall(PetscSectionSetConstraintDof(section, 12, 1));
220: PetscCall(PetscSectionSetFieldDof(section, 3, 0, 2));
221: PetscCall(PetscSectionSetFieldDof(section, 12, 0, 2));
222: PetscCall(PetscSectionSetFieldDof(section, 12, 1, 1));
223: PetscCall(PetscSectionSetFieldDof(section, 25, 1, 2));
224: PetscCall(PetscSectionSetFieldConstraintDof(section, 12, 1, 1));
225: break;
226: case 1:
227: PetscCall(PetscSectionSetDof(section, 0, 2));
228: PetscCall(PetscSectionSetDof(section, 1, 1));
229: PetscCall(PetscSectionSetDof(section, 8, 3));
230: PetscCall(PetscSectionSetDof(section, 20, 2));
231: PetscCall(PetscSectionSetConstraintDof(section, 8, 1));
232: PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2));
233: PetscCall(PetscSectionSetFieldDof(section, 8, 0, 2));
234: PetscCall(PetscSectionSetFieldDof(section, 1, 1, 1));
235: PetscCall(PetscSectionSetFieldDof(section, 8, 1, 1));
236: PetscCall(PetscSectionSetFieldDof(section, 20, 1, 2));
237: PetscCall(PetscSectionSetFieldConstraintDof(section, 8, 1, 1));
238: break;
239: }
240: PetscCall(PetscSectionSetUp(section));
241: {
242: const PetscInt indices[] = {2};
243: const PetscInt indices1[] = {0};
245: switch (rank) {
246: case 0:
247: PetscCall(PetscSectionSetConstraintIndices(section, 12, indices));
248: PetscCall(PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1));
249: break;
250: case 1:
251: PetscCall(PetscSectionSetConstraintIndices(section, 8, indices));
252: PetscCall(PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1));
253: break;
254: }
255: }
256: if (user.shell) {
257: PetscSF sf;
259: PetscCall(DMShellCreate(comm, &sdm));
260: PetscCall(DMGetPointSF(dm, &sf));
261: PetscCall(DMSetPointSF(sdm, sf));
262: } else {
263: PetscCall(DMClone(dm, &sdm));
264: }
265: PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName));
266: PetscCall(DMSetLocalSection(sdm, section));
267: PetscCall(PetscSectionDestroy(§ion));
268: PetscCall(DMPlexSectionView(dm, viewer, sdm));
269: /* Create global vector */
270: PetscCall(DMGetGlobalSection(sdm, &gsection));
271: PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints));
272: if (user.shell) {
273: PetscInt n = -1;
275: PetscCall(VecCreate(comm, &vec));
276: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &n));
277: else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &n));
278: PetscCall(VecSetSizes(vec, n, PETSC_DECIDE));
279: PetscCall(VecSetUp(vec));
280: } else {
281: PetscCall(DMGetGlobalVector(sdm, &vec));
282: }
283: PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
284: PetscCall(VecGetArrayWrite(vec, &array));
285: if (includesConstraints) {
286: switch (rank) {
287: case 0:
288: break;
289: case 1:
290: array[0] = 1.0;
291: array[1] = 1.1;
292: array[2] = 1.2;
293: array[3] = 1.3;
294: array[4] = 1.4;
295: array[5] = 1.5;
296: array[6] = 1.6;
297: array[7] = 1.7;
298: break;
299: }
300: } else {
301: switch (rank) {
302: case 0:
303: break;
304: case 1:
305: array[0] = 1.0;
306: array[1] = 1.1;
307: array[2] = 1.2;
308: array[3] = 1.3;
309: array[4] = 1.4;
310: array[5] = 1.6;
311: array[6] = 1.7;
312: break;
313: }
314: }
315: PetscCall(VecRestoreArrayWrite(vec, &array));
316: PetscCall(DMPlexGlobalVectorView(dm, viewer, sdm, vec));
317: if (user.shell) {
318: PetscCall(VecDestroy(&vec));
319: } else {
320: PetscCall(DMRestoreGlobalVector(sdm, &vec));
321: }
322: PetscCall(DMDestroy(&sdm));
323: }
324: PetscCall(PetscViewerDestroy(&viewer));
325: PetscCall(DMDestroy(&dm));
326: }
327: PetscCallMPI(MPI_Comm_free(&comm));
328: /* Load */
329: mycolor = (PetscMPIInt)(rank >= 3);
330: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
331: if (mycolor == 0) {
332: DM dm;
333: PetscSF sfXC;
334: PetscViewer viewer;
336: PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer));
337: /* Load exampleDMPlex */
338: {
339: PetscSF sfXB, sfBC;
341: PetscCall(DMCreate(comm, &dm));
342: PetscCall(DMSetType(dm, DMPLEX));
343: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
344: /* sfXB: X -> B */
345: /* X: set of globalPointNumbers, [0, N) */
346: /* B: loaded naive in-memory plex */
347: PetscCall(PetscViewerPushFormat(viewer, format));
348: PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB));
349: PetscCall(PetscViewerPopFormat(viewer));
350: {
351: DM distributedDM;
352: PetscInt overlap = 1;
353: PetscPartitioner part;
355: PetscCall(DMPlexGetPartitioner(dm, &part));
356: PetscCall(PetscPartitionerSetFromOptions(part));
357: /* sfBC: B -> C */
358: /* B: loaded naive in-memory plex */
359: /* C: redistributed good in-memory */
360: PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM));
361: if (distributedDM) {
362: PetscCall(DMDestroy(&dm));
363: dm = distributedDM;
364: }
365: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
366: }
367: /* sfXC: X -> C */
368: PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
369: PetscCall(PetscSFDestroy(&sfXB));
370: PetscCall(PetscSFDestroy(&sfBC));
371: }
372: /* Load labels */
373: PetscCall(PetscViewerPushFormat(viewer, format));
374: PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC));
375: PetscCall(PetscViewerPopFormat(viewer));
376: /* Load coordinates */
377: PetscCall(PetscViewerPushFormat(viewer, format));
378: PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC));
379: PetscCall(PetscViewerPopFormat(viewer));
380: PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)"));
381: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
382: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
383: /* Load exampleVec */
384: {
385: DM sdm;
386: PetscSection section, gsection;
387: IS perm;
388: PetscBool includesConstraints = PETSC_FALSE;
389: Vec vec;
390: PetscSF lsf, gsf;
392: if (user.shell) {
393: PetscSF sf;
395: PetscCall(DMShellCreate(comm, &sdm));
396: PetscCall(DMGetPointSF(dm, &sf));
397: PetscCall(DMSetPointSF(sdm, sf));
398: } else {
399: PetscCall(DMClone(dm, &sdm));
400: }
401: PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName));
402: PetscCall(PetscSectionCreate(comm, §ion));
403: {
404: PetscInt pStart = -1, pEnd = -1, p = -1;
405: PetscInt *pinds = NULL;
407: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
408: PetscCall(PetscMalloc1(pEnd - pStart, &pinds));
409: for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p;
410: if (rank == 2) {
411: pinds[10] = 20;
412: pinds[20] = 10;
413: }
414: PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm));
415: }
416: PetscCall(PetscSectionSetPermutation(section, perm));
417: PetscCall(ISDestroy(&perm));
418: PetscCall(DMSetLocalSection(sdm, section));
419: PetscCall(PetscSectionDestroy(§ion));
420: PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf));
421: /* Load as local vector */
422: PetscCall(DMGetLocalSection(sdm, §ion));
423: PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section"));
424: PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
425: PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
426: if (user.shell) {
427: PetscInt m = -1;
429: PetscCall(VecCreate(comm, &vec));
430: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
431: else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
432: PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
433: PetscCall(VecSetUp(vec));
434: } else {
435: PetscCall(DMGetLocalVector(sdm, &vec));
436: }
437: PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
438: PetscCall(VecSet(vec, constraintValue));
439: PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec));
440: PetscCall(PetscSFDestroy(&lsf));
441: if (user.shell) {
442: PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
443: PetscCall(VecDestroy(&vec));
444: } else {
445: PetscCall(DMRestoreLocalVector(sdm, &vec));
446: }
447: /* Load as global vector */
448: PetscCall(DMGetGlobalSection(sdm, &gsection));
449: PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section"));
450: PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm)));
451: PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints));
452: if (user.shell) {
453: PetscInt m = -1;
455: PetscCall(VecCreate(comm, &vec));
456: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m));
457: else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m));
458: PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
459: PetscCall(VecSetUp(vec));
460: } else {
461: PetscCall(DMGetGlobalVector(sdm, &vec));
462: }
463: PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
464: PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec));
465: PetscCall(PetscSFDestroy(&gsf));
466: PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
467: if (user.shell) {
468: PetscCall(VecDestroy(&vec));
469: } else {
470: PetscCall(DMRestoreGlobalVector(sdm, &vec));
471: }
472: PetscCall(DMDestroy(&sdm));
473: }
474: PetscCall(PetscViewerDestroy(&viewer));
475: PetscCall(PetscSFDestroy(&sfXC));
476: PetscCall(DMDestroy(&dm));
477: }
478: PetscCallMPI(MPI_Comm_free(&comm));
480: /* Finalize */
481: PetscCall(PetscFinalize());
482: return 0;
483: }
485: /*TEST
487: build:
488: requires: hdf5
489: testset:
490: suffix: 0
491: requires: !complex
492: nsize: 4
493: args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail
494: args: -dm_plex_view_hdf5_storage_version 2.0.0
495: args: -save_petscpartitioner_type simple
496: test:
497: suffix: parmetis
498: requires: parmetis
499: args: -petscpartitioner_type parmetis
500: test:
501: suffix: ptscotch
502: requires: ptscotch
503: args: -petscpartitioner_type ptscotch
505: TEST*/