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: options->shell = PETSC_FALSE;
140: PetscOptionsBegin(comm, "", "DMPlex View/Load Test Options", "DMPLEX");
141: PetscCall(PetscOptionsString("-fname", "The output mesh file", "ex12.c", options->fname, options->fname, sizeof(options->fname), &flg));
142: PetscCall(PetscOptionsBool("-shell", "Use DMShell to wrap sections", "ex12.c", options->shell, &options->shell, NULL));
143: PetscOptionsEnd();
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: int main(int argc, char **argv)
148: {
149: MPI_Comm comm;
150: PetscMPIInt size, rank, mycolor;
151: const char exampleDMPlexName[] = "exampleDMPlex";
152: const char exampleSectionDMName[] = "exampleSectionDM";
153: const char exampleVecName[] = "exampleVec";
154: PetscScalar constraintValue = 1.5;
155: PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
156: AppCtx user;
158: PetscFunctionBeginUser;
159: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
160: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
161: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
162: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
163: PetscCheck(size >= 3, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes");
165: /* Save */
166: mycolor = (PetscMPIInt)(rank >= 2);
167: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
168: if (mycolor == 0) {
169: DM dm;
170: PetscViewer viewer;
172: PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer));
173: /* Save exampleDMPlex */
174: {
175: DM pdm;
176: const PetscInt faces[2] = {6, 1};
177: PetscSF sf;
178: PetscInt overlap = 1;
179: PetscPartitioner part;
181: PetscCall(DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &dm));
182: PetscCall(DMPlexGetPartitioner(dm, &part));
183: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, "save_"));
184: PetscCall(PetscPartitionerSetFromOptions(part));
185: PetscCall(DMPlexDistribute(dm, overlap, &sf, &pdm));
186: if (pdm) {
187: PetscCall(DMDestroy(&dm));
188: dm = pdm;
189: }
190: PetscCall(PetscSFDestroy(&sf));
191: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
192: PetscCall(PetscViewerPushFormat(viewer, format));
193: PetscCall(DMPlexTopologyView(dm, viewer));
194: PetscCall(DMPlexLabelsView(dm, viewer));
195: PetscCall(PetscViewerPopFormat(viewer));
196: }
197: /* Save coordinates */
198: PetscCall(PetscViewerPushFormat(viewer, format));
199: PetscCall(DMPlexCoordinatesView(dm, viewer));
200: PetscCall(PetscViewerPopFormat(viewer));
201: /* Save exampleVec */
202: {
203: PetscInt pStart = -1, pEnd = -1;
204: DM sdm;
205: PetscSection section, gsection;
206: PetscBool includesConstraints = PETSC_FALSE;
207: Vec vec;
208: PetscScalar *array = NULL;
210: /* Create section */
211: PetscCall(PetscSectionCreate(comm, §ion));
212: PetscCall(PetscSectionSetNumFields(section, 2));
213: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
214: PetscCall(PetscSectionSetChart(section, pStart, pEnd));
215: switch (rank) {
216: case 0:
217: PetscCall(PetscSectionSetDof(section, 3, 2));
218: PetscCall(PetscSectionSetDof(section, 12, 3));
219: PetscCall(PetscSectionSetDof(section, 25, 2));
220: PetscCall(PetscSectionSetConstraintDof(section, 12, 1));
221: PetscCall(PetscSectionSetFieldDof(section, 3, 0, 2));
222: PetscCall(PetscSectionSetFieldDof(section, 12, 0, 2));
223: PetscCall(PetscSectionSetFieldDof(section, 12, 1, 1));
224: PetscCall(PetscSectionSetFieldDof(section, 25, 1, 2));
225: PetscCall(PetscSectionSetFieldConstraintDof(section, 12, 1, 1));
226: break;
227: case 1:
228: PetscCall(PetscSectionSetDof(section, 0, 2));
229: PetscCall(PetscSectionSetDof(section, 1, 1));
230: PetscCall(PetscSectionSetDof(section, 8, 3));
231: PetscCall(PetscSectionSetDof(section, 20, 2));
232: PetscCall(PetscSectionSetConstraintDof(section, 8, 1));
233: PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2));
234: PetscCall(PetscSectionSetFieldDof(section, 8, 0, 2));
235: PetscCall(PetscSectionSetFieldDof(section, 1, 1, 1));
236: PetscCall(PetscSectionSetFieldDof(section, 8, 1, 1));
237: PetscCall(PetscSectionSetFieldDof(section, 20, 1, 2));
238: PetscCall(PetscSectionSetFieldConstraintDof(section, 8, 1, 1));
239: break;
240: }
241: PetscCall(PetscSectionSetUp(section));
242: {
243: const PetscInt indices[] = {2};
244: const PetscInt indices1[] = {0};
246: switch (rank) {
247: case 0:
248: PetscCall(PetscSectionSetConstraintIndices(section, 12, indices));
249: PetscCall(PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1));
250: break;
251: case 1:
252: PetscCall(PetscSectionSetConstraintIndices(section, 8, indices));
253: PetscCall(PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1));
254: break;
255: }
256: }
257: if (user.shell) {
258: PetscSF sf;
260: PetscCall(DMShellCreate(comm, &sdm));
261: PetscCall(DMGetPointSF(dm, &sf));
262: PetscCall(DMSetPointSF(sdm, sf));
263: } else {
264: PetscCall(DMClone(dm, &sdm));
265: }
266: PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName));
267: PetscCall(DMSetLocalSection(sdm, section));
268: PetscCall(PetscSectionDestroy(§ion));
269: PetscCall(DMPlexSectionView(dm, viewer, sdm));
270: /* Create global vector */
271: PetscCall(DMGetGlobalSection(sdm, &gsection));
272: PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints));
273: if (user.shell) {
274: PetscInt n = -1;
276: PetscCall(VecCreate(comm, &vec));
277: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &n));
278: else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &n));
279: PetscCall(VecSetSizes(vec, n, PETSC_DECIDE));
280: PetscCall(VecSetUp(vec));
281: } else {
282: PetscCall(DMGetGlobalVector(sdm, &vec));
283: }
284: PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
285: PetscCall(VecGetArrayWrite(vec, &array));
286: if (includesConstraints) {
287: switch (rank) {
288: case 0:
289: break;
290: case 1:
291: array[0] = 1.0;
292: array[1] = 1.1;
293: array[2] = 1.2;
294: array[3] = 1.3;
295: array[4] = 1.4;
296: array[5] = 1.5;
297: array[6] = 1.6;
298: array[7] = 1.7;
299: break;
300: }
301: } else {
302: switch (rank) {
303: case 0:
304: break;
305: case 1:
306: array[0] = 1.0;
307: array[1] = 1.1;
308: array[2] = 1.2;
309: array[3] = 1.3;
310: array[4] = 1.4;
311: array[5] = 1.6;
312: array[6] = 1.7;
313: break;
314: }
315: }
316: PetscCall(VecRestoreArrayWrite(vec, &array));
317: PetscCall(DMPlexGlobalVectorView(dm, viewer, sdm, vec));
318: if (user.shell) {
319: PetscCall(VecDestroy(&vec));
320: } else {
321: PetscCall(DMRestoreGlobalVector(sdm, &vec));
322: }
323: PetscCall(DMDestroy(&sdm));
324: }
325: PetscCall(PetscViewerDestroy(&viewer));
326: PetscCall(DMDestroy(&dm));
327: }
328: PetscCallMPI(MPI_Comm_free(&comm));
329: /* Load */
330: mycolor = (PetscMPIInt)(rank >= 3);
331: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
332: if (mycolor == 0) {
333: DM dm;
334: PetscSF sfXC;
335: PetscViewer viewer;
337: PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer));
338: /* Load exampleDMPlex */
339: {
340: PetscSF sfXB, sfBC;
342: PetscCall(DMCreate(comm, &dm));
343: PetscCall(DMSetType(dm, DMPLEX));
344: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
345: /* sfXB: X -> B */
346: /* X: set of globalPointNumbers, [0, N) */
347: /* B: loaded naive in-memory plex */
348: PetscCall(PetscViewerPushFormat(viewer, format));
349: PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB));
350: PetscCall(PetscViewerPopFormat(viewer));
351: {
352: DM distributedDM;
353: PetscInt overlap = 1;
354: PetscPartitioner part;
356: PetscCall(DMPlexGetPartitioner(dm, &part));
357: PetscCall(PetscPartitionerSetFromOptions(part));
358: /* sfBC: B -> C */
359: /* B: loaded naive in-memory plex */
360: /* C: redistributed good in-memory */
361: PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM));
362: if (distributedDM) {
363: PetscCall(DMDestroy(&dm));
364: dm = distributedDM;
365: }
366: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
367: }
368: /* sfXC: X -> C */
369: PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
370: PetscCall(PetscSFDestroy(&sfXB));
371: PetscCall(PetscSFDestroy(&sfBC));
372: }
373: /* Load labels */
374: PetscCall(PetscViewerPushFormat(viewer, format));
375: PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC));
376: PetscCall(PetscViewerPopFormat(viewer));
377: /* Load coordinates */
378: PetscCall(PetscViewerPushFormat(viewer, format));
379: PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC));
380: PetscCall(PetscViewerPopFormat(viewer));
381: PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)"));
382: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
383: PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
384: /* Load exampleVec */
385: {
386: DM sdm;
387: PetscSection section, gsection;
388: IS perm;
389: PetscBool includesConstraints = PETSC_FALSE;
390: Vec vec;
391: PetscSF lsf, gsf;
393: if (user.shell) {
394: PetscSF sf;
396: PetscCall(DMShellCreate(comm, &sdm));
397: PetscCall(DMGetPointSF(dm, &sf));
398: PetscCall(DMSetPointSF(sdm, sf));
399: } else {
400: PetscCall(DMClone(dm, &sdm));
401: }
402: PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName));
403: PetscCall(PetscSectionCreate(comm, §ion));
404: {
405: PetscInt pStart = -1, pEnd = -1, p = -1;
406: PetscInt *pinds = NULL;
408: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
409: PetscCall(PetscMalloc1(pEnd - pStart, &pinds));
410: for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p;
411: if (rank == 2) {
412: pinds[10] = 20;
413: pinds[20] = 10;
414: }
415: PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm));
416: }
417: PetscCall(PetscSectionSetPermutation(section, perm));
418: PetscCall(ISDestroy(&perm));
419: PetscCall(DMSetLocalSection(sdm, section));
420: PetscCall(PetscSectionDestroy(§ion));
421: PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf));
422: /* Load as local vector */
423: PetscCall(DMGetLocalSection(sdm, §ion));
424: PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section"));
425: PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
426: PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
427: if (user.shell) {
428: PetscInt m = -1;
430: PetscCall(VecCreate(comm, &vec));
431: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
432: else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
433: PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
434: PetscCall(VecSetUp(vec));
435: } else {
436: PetscCall(DMGetLocalVector(sdm, &vec));
437: }
438: PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
439: PetscCall(VecSet(vec, constraintValue));
440: PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec));
441: PetscCall(PetscSFDestroy(&lsf));
442: if (user.shell) {
443: PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
444: PetscCall(VecDestroy(&vec));
445: } else {
446: PetscCall(DMRestoreLocalVector(sdm, &vec));
447: }
448: /* Load as global vector */
449: PetscCall(DMGetGlobalSection(sdm, &gsection));
450: PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section"));
451: PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm)));
452: PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints));
453: if (user.shell) {
454: PetscInt m = -1;
456: PetscCall(VecCreate(comm, &vec));
457: if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m));
458: else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m));
459: PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
460: PetscCall(VecSetUp(vec));
461: } else {
462: PetscCall(DMGetGlobalVector(sdm, &vec));
463: }
464: PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
465: PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec));
466: PetscCall(PetscSFDestroy(&gsf));
467: PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
468: if (user.shell) {
469: PetscCall(VecDestroy(&vec));
470: } else {
471: PetscCall(DMRestoreGlobalVector(sdm, &vec));
472: }
473: PetscCall(DMDestroy(&sdm));
474: }
475: PetscCall(PetscViewerDestroy(&viewer));
476: PetscCall(PetscSFDestroy(&sfXC));
477: PetscCall(DMDestroy(&dm));
478: }
479: PetscCallMPI(MPI_Comm_free(&comm));
481: /* Finalize */
482: PetscCall(PetscFinalize());
483: return 0;
484: }
486: /*TEST
488: build:
489: requires: hdf5
490: testset:
491: suffix: 0
492: requires: !complex
493: nsize: 4
494: args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail
495: args: -dm_plex_view_hdf5_storage_version 2.0.0
496: args: -save_petscpartitioner_type simple
497: test:
498: suffix: parmetis
499: requires: parmetis
500: args: -petscpartitioner_type parmetis
501: test:
502: suffix: ptscotch
503: requires: ptscotch
504: args: -petscpartitioner_type ptscotch
506: TEST*/