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;
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;
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};
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;
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;
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;
332: PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer));
333: /* Load exampleDMPlex */
334: {
335: PetscSF sfXB, sfBC;
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;
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;
388: if (user.shell) {
389: PetscSF sf;
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;
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;
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;
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));
476: /* Finalize */
477: PetscCall(PetscFinalize());
478: return 0;
479: }
481: /*TEST
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
500: TEST*/