Actual source code: ex21.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;

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, &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, &section));
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(&section));
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:       PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
347:       {
348:         DM               distributedDM;
349:         PetscInt         overlap = 1;
350:         PetscPartitioner part;

352:         PetscCall(DMPlexGetPartitioner(dm, &part));
353:         PetscCall(PetscPartitionerSetFromOptions(part));
354:         /* sfBC: B -> C                    */
355:         /* B: loaded naive in-memory plex  */
356:         /* C: redistributed good in-memory */
357:         PetscCall(DMPlexDistribute(dm, overlap, &sfBC, &distributedDM));
358:         if (distributedDM) {
359:           PetscCall(DMDestroy(&dm));
360:           dm = distributedDM;
361:         }
362:         PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
363:       }
364:       /* sfXC: X -> C */
365:       PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC));
366:       PetscCall(PetscSFDestroy(&sfXB));
367:       PetscCall(PetscSFDestroy(&sfBC));
368:     }
369:     /* Load labels */
370:     PetscCall(PetscViewerPushFormat(viewer, format));
371:     PetscCall(DMPlexLabelsLoad(dm, viewer, sfXC));
372:     PetscCall(PetscViewerPopFormat(viewer));
373:     /* Load coordinates */
374:     PetscCall(PetscViewerPushFormat(viewer, format));
375:     PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXC));
376:     PetscCall(PetscViewerPopFormat(viewer));
377:     PetscCall(PetscObjectSetName((PetscObject)dm, "Load: DM (with coordinates)"));
378:     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
379:     PetscCall(PetscObjectSetName((PetscObject)dm, exampleDMPlexName));
380:     /* Load exampleVec */
381:     {
382:       DM           sdm;
383:       PetscSection section, gsection;
384:       IS           perm;
385:       PetscBool    includesConstraints = PETSC_FALSE;
386:       Vec          vec;
387:       PetscSF      lsf, gsf;

389:       if (user.shell) {
390:         PetscSF sf;

392:         PetscCall(DMShellCreate(comm, &sdm));
393:         PetscCall(DMGetPointSF(dm, &sf));
394:         PetscCall(DMSetPointSF(sdm, sf));
395:       } else {
396:         PetscCall(DMClone(dm, &sdm));
397:       }
398:       PetscCall(PetscObjectSetName((PetscObject)sdm, exampleSectionDMName));
399:       PetscCall(PetscSectionCreate(comm, &section));
400:       {
401:         PetscInt  pStart = -1, pEnd = -1, p = -1;
402:         PetscInt *pinds = NULL;

404:         PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
405:         PetscCall(PetscMalloc1(pEnd - pStart, &pinds));
406:         for (p = 0; p < pEnd - pStart; ++p) pinds[p] = p;
407:         if (rank == 2) {
408:           pinds[10] = 20;
409:           pinds[20] = 10;
410:         }
411:         PetscCall(ISCreateGeneral(comm, pEnd - pStart, pinds, PETSC_OWN_POINTER, &perm));
412:       }
413:       PetscCall(PetscSectionSetPermutation(section, perm));
414:       PetscCall(ISDestroy(&perm));
415:       PetscCall(DMSetLocalSection(sdm, section));
416:       PetscCall(PetscSectionDestroy(&section));
417:       PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXC, &gsf, &lsf));
418:       /* Load as local vector */
419:       PetscCall(DMGetLocalSection(sdm, &section));
420:       PetscCall(PetscObjectSetName((PetscObject)section, "Load: local section"));
421:       PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
422:       PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
423:       if (user.shell) {
424:         PetscInt m = -1;

426:         PetscCall(VecCreate(comm, &vec));
427:         if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
428:         else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
429:         PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
430:         PetscCall(VecSetUp(vec));
431:       } else {
432:         PetscCall(DMGetLocalVector(sdm, &vec));
433:       }
434:       PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
435:       PetscCall(VecSet(vec, constraintValue));
436:       PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, lsf, vec));
437:       PetscCall(PetscSFDestroy(&lsf));
438:       if (user.shell) {
439:         PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
440:         PetscCall(VecDestroy(&vec));
441:       } else {
442:         PetscCall(DMRestoreLocalVector(sdm, &vec));
443:       }
444:       /* Load as global vector */
445:       PetscCall(DMGetGlobalSection(sdm, &gsection));
446:       PetscCall(PetscObjectSetName((PetscObject)gsection, "Load: global section"));
447:       PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm)));
448:       PetscCall(PetscSectionGetIncludesConstraints(gsection, &includesConstraints));
449:       if (user.shell) {
450:         PetscInt m = -1;

452:         PetscCall(VecCreate(comm, &vec));
453:         if (includesConstraints) PetscCall(PetscSectionGetStorageSize(gsection, &m));
454:         else PetscCall(PetscSectionGetConstrainedStorageSize(gsection, &m));
455:         PetscCall(VecSetSizes(vec, m, PETSC_DECIDE));
456:         PetscCall(VecSetUp(vec));
457:       } else {
458:         PetscCall(DMGetGlobalVector(sdm, &vec));
459:       }
460:       PetscCall(PetscObjectSetName((PetscObject)vec, exampleVecName));
461:       PetscCall(DMPlexGlobalVectorLoad(dm, viewer, sdm, gsf, vec));
462:       PetscCall(PetscSFDestroy(&gsf));
463:       PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_(comm)));
464:       if (user.shell) {
465:         PetscCall(VecDestroy(&vec));
466:       } else {
467:         PetscCall(DMRestoreGlobalVector(sdm, &vec));
468:       }
469:       PetscCall(DMDestroy(&sdm));
470:     }
471:     PetscCall(PetscViewerDestroy(&viewer));
472:     PetscCall(PetscSFDestroy(&sfXC));
473:     PetscCall(DMDestroy(&dm));
474:   }
475:   PetscCallMPI(MPI_Comm_free(&comm));

477:   /* Finalize */
478:   PetscCall(PetscFinalize());
479:   return 0;
480: }

482: /*TEST

484:   build:
485:     requires: hdf5
486:   testset:
487:     suffix: 0
488:     requires: !complex
489:     nsize: 4
490:     args: -fname ex12_dump.h5 -shell {{True False}separate output} -dm_view ascii::ascii_info_detail
491:     args: -dm_plex_view_hdf5_storage_version 2.0.0
492:     test:
493:       suffix: parmetis
494:       requires: parmetis
495:       args: -petscpartitioner_type parmetis
496:     test:
497:       suffix: ptscotch
498:       requires: ptscotch
499:       args: -petscpartitioner_type ptscotch

501: TEST*/