Actual source code: ex12.c
1: static char help[] = "Tests save/load of plex/section/vec 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: rakn 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;
136: PetscErrorCode ierr;
139: options->fname[0] = '\0';
140: PetscOptionsBegin(comm, "", "DMPlex View/Load Test Options", "DMPLEX");
141: PetscOptionsString("-fname", "The output mesh file", "ex12.c", options->fname, options->fname, sizeof(options->fname), &flg);
142: PetscOptionsBool("-shell", "Use DMShell to wrap sections", "ex12.c", options->shell, &options->shell, NULL);
143: PetscOptionsEnd();
144: return(0);
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;
157: PetscErrorCode ierr;
159: PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr;
160: ProcessOptions(PETSC_COMM_WORLD, &user);
161: MPI_Comm_size(PETSC_COMM_WORLD, &size);
162: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
163: if (size < 3) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes");
165: /* Save */
166: mycolor = (PetscMPIInt)(rank >= 2);
167: MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm);
168: if (mycolor == 0) {
169: DM dm;
170: PetscViewer viewer;
172: 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;
180: DMPlexCreateBoxMesh(comm, 2, PETSC_FALSE, faces, NULL, NULL, NULL, PETSC_TRUE, &dm);
181: DMPlexDistribute(dm, overlap, &sf, &pdm);
182: if (pdm) {
183: DMDestroy(&dm);
184: dm = pdm;
185: }
186: PetscSFDestroy(&sf);
187: PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
188: PetscViewerPushFormat(viewer, format);
189: DMPlexTopologyView(dm, viewer);
190: DMPlexLabelsView(dm, viewer);
191: PetscViewerPopFormat(viewer);
192: }
193: /* Save coordinates */
194: /* The following block is to replace DMPlexCoordinatesView(). */
195: {
196: DM cdm;
197: Vec coords, newcoords;
198: PetscInt m = -1, M = -1, bs = -1;
199: PetscReal lengthScale = -1;
201: DMGetCoordinateDM(dm, &cdm);
202: PetscObjectSetName((PetscObject)cdm, "coordinateDM");
203: DMPlexSectionView(dm, viewer, cdm);
204: DMGetCoordinates(dm, &coords);
205: VecCreate(PetscObjectComm((PetscObject)coords), &newcoords);
206: PetscObjectSetName((PetscObject)newcoords, "coordinates");
207: VecGetSize(coords, &M);
208: VecGetLocalSize(coords, &m);
209: VecSetSizes(newcoords, m, M);
210: VecGetBlockSize(coords, &bs);
211: VecSetBlockSize(newcoords, bs);
212: VecSetType(newcoords,VECSTANDARD);
213: VecCopy(coords, newcoords);
214: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
215: VecScale(newcoords, lengthScale);
216: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
217: DMPlexGlobalVectorView(dm, viewer, cdm, newcoords);
218: PetscViewerPopFormat(viewer);
219: VecDestroy(&newcoords);
220: }
221: /* Save exampleVec */
222: {
223: PetscInt pStart = -1, pEnd = -1;
224: DM sdm;
225: PetscSection section, gsection;
226: PetscBool includesConstraints = PETSC_FALSE;
227: Vec vec;
228: PetscScalar *array = NULL;
230: /* Create section */
231: PetscSectionCreate(comm, §ion);
232: PetscSectionSetNumFields(section, 2);
233: DMPlexGetChart(dm, &pStart, &pEnd);
234: PetscSectionSetChart(section, pStart, pEnd);
235: switch (rank) {
236: case 0:
237: PetscSectionSetDof(section, 3, 2);
238: PetscSectionSetDof(section, 12, 3);
239: PetscSectionSetDof(section, 25, 2);
240: PetscSectionSetConstraintDof(section, 12, 1);
241: PetscSectionSetFieldDof(section, 3, 0, 2);
242: PetscSectionSetFieldDof(section, 12, 0, 2);
243: PetscSectionSetFieldDof(section, 12, 1, 1);
244: PetscSectionSetFieldDof(section, 25, 1, 2);
245: PetscSectionSetFieldConstraintDof(section, 12, 1, 1);
246: break;
247: case 1:
248: PetscSectionSetDof(section, 0, 2);
249: PetscSectionSetDof(section, 1, 1);
250: PetscSectionSetDof(section, 8, 3);
251: PetscSectionSetDof(section, 20, 2);
252: PetscSectionSetConstraintDof(section, 8, 1);
253: PetscSectionSetFieldDof(section, 0, 0, 2);
254: PetscSectionSetFieldDof(section, 8, 0, 2);
255: PetscSectionSetFieldDof(section, 1, 1, 1);
256: PetscSectionSetFieldDof(section, 8, 1, 1);
257: PetscSectionSetFieldDof(section, 20, 1, 2);
258: PetscSectionSetFieldConstraintDof(section, 8, 1, 1);
259: break;
260: }
261: PetscSectionSetUp(section);
262: {
263: const PetscInt indices[] = {2};
264: const PetscInt indices1[] = {0};
266: switch (rank) {
267: case 0:
268: PetscSectionSetConstraintIndices(section, 12, indices);
269: PetscSectionSetFieldConstraintIndices(section, 12, 1, indices1);
270: break;
271: case 1:
272: PetscSectionSetConstraintIndices(section, 8, indices);
273: PetscSectionSetFieldConstraintIndices(section, 8, 1, indices1);
274: break;
275: }
276: }
277: if (user.shell) {
278: PetscSF sf;
280: DMShellCreate(comm, &sdm);
281: DMGetPointSF(dm, &sf);
282: DMSetPointSF(sdm, sf);
283: }
284: else {
285: DMClone(dm, &sdm);
286: }
287: PetscObjectSetName((PetscObject)sdm, exampleSectionDMName);
288: DMSetLocalSection(sdm, section);
289: PetscSectionDestroy(§ion);
290: DMPlexSectionView(dm, viewer, sdm);
291: /* Create global vector */
292: DMGetGlobalSection(sdm, &gsection);
293: PetscSectionGetIncludesConstraints(gsection, &includesConstraints);
294: if (user.shell) {
295: PetscInt n = -1;
297: VecCreate(comm, &vec);
298: if (includesConstraints) {PetscSectionGetStorageSize(gsection, &n);}
299: else {PetscSectionGetConstrainedStorageSize(gsection, &n);}
300: VecSetSizes(vec, n, PETSC_DECIDE);
301: VecSetUp(vec);
302: } else {
303: DMGetGlobalVector(sdm, &vec);
304: }
305: PetscObjectSetName((PetscObject)vec, exampleVecName);
306: VecGetArrayWrite(vec, &array);
307: if (includesConstraints) {
308: switch (rank) {
309: case 0:
310: break;
311: case 1:
312: array[0] = 1.0;
313: array[1] = 1.1;
314: array[2] = 1.2;
315: array[3] = 1.3;
316: array[4] = 1.4;
317: array[5] = 1.5;
318: array[6] = 1.6;
319: array[7] = 1.7;
320: break;
321: }
322: } else {
323: switch (rank) {
324: case 0:
325: break;
326: case 1:
327: array[0] = 1.0;
328: array[1] = 1.1;
329: array[2] = 1.2;
330: array[3] = 1.3;
331: array[4] = 1.4;
332: array[5] = 1.6;
333: array[6] = 1.7;
334: break;
335: }
336: }
337: VecRestoreArrayWrite(vec, &array);
338: DMPlexGlobalVectorView(dm, viewer, sdm, vec);
339: if (user.shell) {
340: VecDestroy(&vec);
341: } else {
342: DMRestoreGlobalVector(sdm, &vec);
343: }
344: DMDestroy(&sdm);
345: }
346: PetscViewerDestroy(&viewer);
347: DMDestroy(&dm);
348: }
349: MPI_Comm_free(&comm);
350: /* Load */
351: mycolor = (PetscMPIInt)(rank >= 3);
352: MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm);
353: if (mycolor == 0) {
354: DM dm;
355: PetscSF sfXC;
356: PetscViewer viewer;
358: PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer);
359: /* Load exampleDMPlex */
360: {
361: PetscSF sfXB, sfBC;
363: DMCreate(comm, &dm);
364: DMSetType(dm, DMPLEX);
365: PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
366: /* sfXB: X -> B */
367: /* X: set of globalPointNumbers, [0, N) */
368: /* B: loaded naive in-memory plex */
369: PetscViewerPushFormat(viewer, format);
370: DMPlexTopologyLoad(dm, viewer, &sfXB);
371: DMPlexLabelsLoad(dm, viewer);
372: PetscViewerPopFormat(viewer);
373: PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
374: {
375: DM distributedDM;
376: PetscInt overlap = 1;
377: PetscPartitioner part;
379: DMPlexGetPartitioner(dm, &part);
380: PetscPartitionerSetFromOptions(part);
381: /* sfBC: B -> C */
382: /* B: loaded naive in-memory plex */
383: /* C: redistributed good in-memory */
384: DMPlexDistribute(dm, overlap, &sfBC, &distributedDM);
385: if (distributedDM) {
386: DMDestroy(&dm);
387: dm = distributedDM;
388: }
389: PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
390: }
391: /* sfXC: X -> C */
392: PetscSFCompose(sfXB, sfBC, &sfXC);
393: PetscSFDestroy(&sfXB);
394: PetscSFDestroy(&sfBC);
395: }
396: /* Load coordinates */
397: /* The following block is to replace DMPlexCoordinatesLoad() */
398: {
399: DM cdm;
400: PetscSection coordSection;
401: Vec coords;
402: PetscInt m = -1;
403: PetscReal lengthScale = -1;
404: PetscSF lsf, gsf;
406: DMGetCoordinateDM(dm, &cdm);
407: PetscObjectSetName((PetscObject)cdm, "coordinateDM");
408: /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
409: /* gsf: on-disk data -> in-memory global vector associated with cdm's global section */
410: DMPlexSectionLoad(dm, viewer, cdm, sfXC, &gsf, &lsf);
411: VecCreate(comm, &coords);
412: PetscObjectSetName((PetscObject)coords, "coordinates");
413: DMGetLocalSection(cdm, &coordSection);
414: PetscSectionGetStorageSize(coordSection, &m);
415: VecSetSizes(coords, m, PETSC_DECIDE);
416: VecSetUp(coords);
417: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
418: DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords);
419: PetscViewerPopFormat(viewer);
420: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
421: VecScale(coords, 1.0/lengthScale);
422: DMSetCoordinatesLocal(dm, coords);
423: VecDestroy(&coords);
424: PetscSFDestroy(&lsf);
425: PetscSFDestroy(&gsf);
426: PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)");
427: DMViewFromOptions(dm, NULL, "-dm_view");
428: PetscObjectSetName((PetscObject)dm, exampleDMPlexName);
429: }
430: /* Load exampleVec */
431: {
432: DM sdm;
433: PetscSection section, gsection;
434: IS perm;
435: PetscBool includesConstraints = PETSC_FALSE;
436: Vec vec;
437: PetscSF lsf, gsf;
439: if (user.shell) {
440: PetscSF sf;
442: DMShellCreate(comm, &sdm);
443: DMGetPointSF(dm, &sf);
444: DMSetPointSF(sdm, sf);
445: } else {
446: DMClone(dm, &sdm);
447: }
448: PetscObjectSetName((PetscObject)sdm, exampleSectionDMName);
449: PetscSectionCreate(comm, §ion);
450: {
451: PetscInt pStart = -1, pEnd = -1, p = -1;
452: PetscInt *pinds = NULL;
454: DMPlexGetChart(dm, &pStart, &pEnd);
455: PetscMalloc1(pEnd - pStart, &pinds);
456: for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p;
457: if (rank == 2) {pinds[10] = 20; pinds[20] = 10;}
458: ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm);
459: }
460: PetscSectionSetPermutation(section, perm);
461: ISDestroy(&perm);
462: DMSetLocalSection(sdm, section);
463: PetscSectionDestroy(§ion);
464: DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf);
465: /* Load as local vector */
466: DMGetLocalSection(sdm, §ion);
467: PetscObjectSetName((PetscObject)section, "Load: local section");
468: PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm));
469: PetscSectionGetIncludesConstraints(section, &includesConstraints);
470: if (user.shell) {
471: PetscInt m = -1;
473: VecCreate(comm, &vec);
474: if (includesConstraints) {PetscSectionGetStorageSize(section, &m);}
475: else {PetscSectionGetConstrainedStorageSize(section, &m);}
476: VecSetSizes(vec, m, PETSC_DECIDE);
477: VecSetUp(vec);
478: } else {
479: DMGetLocalVector(sdm, &vec);
480: }
481: PetscObjectSetName((PetscObject)vec, exampleVecName);
482: VecSet(vec, constraintValue);
483: DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec);
484: PetscSFDestroy(&lsf);
485: if (user.shell) {
486: VecView(vec, PETSC_VIEWER_STDOUT_(comm));
487: VecDestroy(&vec);
488: } else {
489: DMRestoreLocalVector(sdm, &vec);
490: }
491: /* Load as global vector */
492: DMGetGlobalSection(sdm, &gsection);
493: PetscObjectSetName((PetscObject)gsection, "Load: global section");
494: PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm));
495: PetscSectionGetIncludesConstraints(gsection, &includesConstraints);
496: if (user.shell) {
497: PetscInt m = -1;
499: VecCreate(comm, &vec);
500: if (includesConstraints) {PetscSectionGetStorageSize(gsection, &m);}
501: else {PetscSectionGetConstrainedStorageSize(gsection, &m);}
502: VecSetSizes(vec, m, PETSC_DECIDE);
503: VecSetUp(vec);
504: } else {
505: DMGetGlobalVector(sdm, &vec);
506: }
507: PetscObjectSetName((PetscObject)vec, exampleVecName);
508: DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec);
509: PetscSFDestroy(&gsf);
510: VecView(vec, PETSC_VIEWER_STDOUT_(comm));
511: if (user.shell) {
512: VecDestroy(&vec);
513: } else {
514: DMRestoreGlobalVector(sdm, &vec);
515: }
516: DMDestroy(&sdm);
517: }
518: PetscViewerDestroy(&viewer);
519: PetscSFDestroy(&sfXC);
520: DMDestroy(&dm);
521: }
522: MPI_Comm_free(&comm);
524: /* Finalize */
525: PetscFinalize();
526: return ierr;
527: }
529: /*TEST
531: build:
532: requires: hdf5
533: testset:
534: suffix: 0
535: requires: !complex
536: nsize: 4
537: args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail
538: test:
539: suffix: parmetis
540: requires: parmetis
541: args: -petscpartitioner_type parmetis
542: test:
543: suffix: ptscotch
544: requires: ptscotch
545: args: -petscpartitioner_type ptscotch
547: TEST*/