xref: /petsc/src/dm/impls/plex/tutorials/ex5.c (revision b698fc57f0bea7237255b29c1b77df0acc362ffd)
1 static char help[] = "Demonstrate HDF5/XDMF load-save-reload cycle\n\n";
2 
3 #include <petscdmplex.h>
4 #include <petscviewerhdf5.h>
5 #define EX "ex5.c"
6 
7 typedef struct {
8   char              infile[PETSC_MAX_PATH_LEN];  /* Input mesh filename */
9   char              outfile[PETSC_MAX_PATH_LEN]; /* Dump/reload mesh filename */
10   PetscViewerFormat informat;                    /* Input mesh format */
11   PetscViewerFormat outformat;                   /* Dump/reload mesh format */
12   PetscBool         redistribute;                /* Redistribute the mesh */
13   PetscInt          ntimes;                      /* How many times do the cycle */
14 } AppCtx;
15 
16 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
17 {
18   PetscBool      flg;
19   PetscErrorCode ierr;
20 
21   PetscFunctionBeginUser;
22   options->infile[0]    = '\0';
23   options->outfile[0]   = '\0';
24   options->informat     = PETSC_VIEWER_HDF5_XDMF;
25   options->outformat    = PETSC_VIEWER_HDF5_XDMF;
26   options->redistribute = PETSC_TRUE;
27   options->ntimes       = 2;
28   ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
29   ierr = PetscOptionsString("-infile", "The input mesh file", EX, options->infile, options->infile, sizeof(options->infile), &flg);CHKERRQ(ierr);
30   if (!flg) SETERRQ(comm, PETSC_ERR_USER_INPUT, "-infile needs to be specified");
31   ierr = PetscOptionsString("-outfile", "The output mesh file (by default it's the same as infile)", EX, options->outfile, options->outfile, sizeof(options->outfile), &flg);CHKERRQ(ierr);
32   if (!flg) SETERRQ(comm, PETSC_ERR_USER_INPUT, "-outfile needs to be specified");
33   ierr = PetscOptionsEnum("-informat", "Input mesh format", EX, PetscViewerFormats, (PetscEnum)options->informat, (PetscEnum*)&options->informat, NULL);CHKERRQ(ierr);
34   ierr = PetscOptionsEnum("-outformat", "Dump/reload mesh format", EX, PetscViewerFormats, (PetscEnum)options->outformat, (PetscEnum*)&options->outformat, NULL);CHKERRQ(ierr);
35   ierr = PetscOptionsBool("-redistribute", "Redistribute the mesh", EX, options->redistribute, &options->redistribute, NULL);CHKERRQ(ierr);
36   ierr = PetscOptionsInt("-ntimes", "How many times do the cycle", EX, options->ntimes, &options->ntimes, NULL);CHKERRQ(ierr);
37   ierr = PetscOptionsEnd();
38   PetscFunctionReturn(0);
39 };
40 
41 //TODO test DMLabel I/O (not yet working for PETSC_VIEWER_HDF5_XDMF)
42 int main(int argc, char **argv)
43 {
44   AppCtx            user;
45   PetscInt          i;
46   PetscBool         flg;
47   MPI_Comm          comm;
48   PetscMPIInt       size;
49   PetscErrorCode    ierr;
50   const char        *filename;
51   PetscViewerFormat format;
52 
53   ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
54   comm = PETSC_COMM_WORLD;
55   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
56   ierr = ProcessOptions(comm, &user);CHKERRQ(ierr);
57 
58   /* Use infile for the initial load */
59   filename = user.infile;
60   format   = user.informat;
61 
62   for (i=0; i<user.ntimes; i++) {
63     DM                dm;
64     PetscPartitioner  part;
65     PetscViewer       v;
66 
67     ierr = PetscPrintf(comm, "Begin cycle %D\n",i);CHKERRQ(ierr);
68 
69     /* Load data from XDMF into dm in parallel */
70     /* We could also use
71         ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, user.filename, PETSC_TRUE, &dm);CHKERRQ(ierr);
72       This currently support a few more formats than DMLoad().
73     */
74     ierr = PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &v);CHKERRQ(ierr);
75     ierr = PetscViewerPushFormat(v, format);CHKERRQ(ierr);
76     ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
77     ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
78     ierr = DMSetOptionsPrefix(dm,"loaded_");CHKERRQ(ierr);
79     ierr = DMLoad(dm, v);CHKERRQ(ierr);
80     ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
81     ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
82     ierr = PetscViewerPopFormat(v);CHKERRQ(ierr);
83     ierr = PetscViewerDestroy(&v);CHKERRQ(ierr);
84 
85     /* Use outfile for all I/O except the very initial load */
86     filename = user.outfile;
87     format   = user.outformat;
88 
89     /* We just test/demonstrate DM is indeed distributed - unneeded in the application code */
90     ierr = DMPlexIsDistributed(dm, &flg);CHKERRQ(ierr);
91     ierr = PetscPrintf(comm, "Loaded mesh distributed? %s\n", PetscBools[flg]);
92 
93     /* Interpolate */
94     //TODO we want to be able to do this from options in DMSetFromOptions() probably
95     //TODO we want to be able to do this in-place
96     {
97       DM idm;
98 
99       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
100       ierr = DMDestroy(&dm);CHKERRQ(ierr);
101       dm   = idm;
102         ierr = DMSetOptionsPrefix(dm,"interpolated_");CHKERRQ(ierr);
103         ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
104         ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
105     }
106 
107     /* Redistribute */
108     //TODO we want to be able to do this from options in DMSetFromOptions() probably
109     if (user.redistribute) {
110       DM dmdist;
111 
112       ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr);
113       ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
114       ierr = DMPlexDistribute(dm, 0, NULL, &dmdist);CHKERRQ(ierr);
115       //TODO we want to be able to do this in-place
116       if (dmdist) {
117         ierr = DMDestroy(&dm);CHKERRQ(ierr);
118         dm   = dmdist;
119         ierr = DMSetOptionsPrefix(dm,"redistributed_");CHKERRQ(ierr);
120         ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
121         ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
122       }
123     }
124 
125     /* Save redistributed dm to XDMF in parallel and destroy it */
126     ierr = PetscViewerHDF5Open(comm, user.outfile, FILE_MODE_WRITE, &v);CHKERRQ(ierr);
127     ierr = PetscViewerPushFormat(v, format);CHKERRQ(ierr);
128     ierr = DMView(dm, v);CHKERRQ(ierr);
129     ierr = PetscViewerPopFormat(v);CHKERRQ(ierr);
130     ierr = PetscViewerDestroy(&v);CHKERRQ(ierr);
131     ierr = DMDestroy(&dm);CHKERRQ(ierr);
132 
133     ierr = PetscPrintf(comm, "End   cycle %D\n--------\n",i);CHKERRQ(ierr);
134   }
135 
136   /* Final clean-up */
137   ierr = PetscFinalize();
138   return ierr;
139 }
140 
141 /*TEST
142   build:
143     requires: hdf5
144   testset:
145     suffix: 0
146     requires: !complex
147     nsize: {{2 4}}
148     args: -infile ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.h5 -informat hdf5_xdmf
149     args: -outfile ex5_dump.h5 -outformat {{hdf5_xdmf hdf5_petsc}separate output}
150     args: -ntimes 3
151     test:
152       # this partitioner should not shuffle anything, it should yield the same partititioning as the XDMF reader - added just for testing
153       suffix: simple
154       args: -petscpartitioner_type simple
155     test:
156       suffix: parmetis
157       requires: parmetis
158       args: -petscpartitioner_type parmetis
159     test:
160       suffix: ptscotch
161       requires: ptscotch
162       args: -petscpartitioner_type ptscotch
163 TEST*/
164