xref: /petsc/src/vec/is/tests/ex5.c (revision dd39110ba674af879e763a5939b18b2b35b4c87f)
1 static char help[] = "Tests PetscSectionView()/Load() with HDF5.\n\n";
2 
3 #include <petscdmshell.h>
4 #include <petscdmplex.h>
5 #include <petscsection.h>
6 #include <petscsf.h>
7 #include <petsclayouthdf5.h>
8 
9 /* Save/Load abstract sections
10 
11 =====================
12  Save on 2 processes
13 =====================
14 
15 section:
16                          0   1   2   3
17   rank 0: Dof (Field 0)  2   3   5   7
18           Dof (Field 1)  1   0   0   0
19 
20                          0   1   2
21   rank 1: Dof (Field 0)  7   5  11 <- DoF 7 is constrained
22           Dof (Field 1)  0   0   2
23 
24 sf:
25   [0] 3 <- (1, 0)
26   [1] 1 <- (0, 2)
27 
28 global section (includesConstraints = PETSC_FALSE):
29                          0   1   2   3
30   rank 0: Dof (Field 0)  2   3   5  -8
31           Off (Field 0)  0   3   6  -12
32           Dof (Field 1)  1   0   0  -1
33           Off (Field 1)  2   6  11  -19
34 
35                          0   1   2
36   rank 1: Dof (Field 0)  7  -6  11
37           Off (Field 0) 11  -7  18
38           Dof (Field 1)  0  -1   2
39           Off (Field 1) 18 -12  28
40 
41 global section (includesConstraints = PETSC_TRUE):
42                          0   1   2   3
43   rank 0: Dof (Field 0)  2   3   5  -8
44           Off (Field 0)  0   3   6  -12
45           Dof (Field 1)  1   0   0  -1
46           Off (Field 1)  2   6  11  -19
47 
48                          0   1   2
49   rank 1: Dof (Field 0)  7  -6  11
50           Off (Field 0) 11  -7  18
51           Dof (Field 1)  0  -1   2
52           Off (Field 1) 18 -12  29
53 
54 =====================
55  Load on 3 Processes
56 =====================
57 
58 (Set chartSize = 4, 0, 1 for rank 0, 1, 2, respectively)
59 
60 global section (includesConstraints = PETSC_FALSE):
61 
62   rank 0: Dof (Field 0)  2   3   5   7
63           Off (Field 0)  0   3   6  11
64           Dof (Field 1)  1   0   0   0
65           Off (Field 1)  2   6  11  18
66 
67   rank 1: Dof (Field 0)
68           Dof (Field 1)
69 
70   rank 2: Dof (Field 0) 11
71           Off (Field 0) 18
72           Dof (Field 1)  2
73           Off (Field 1) 28
74 
75 global section (includesConstraints = PETSC_TRUE):
76 
77   rank 0: Dof (Field 0)  2   3   5   7
78           Off (Field 0)  0   3   6  11
79           Dof (Field 1)  1   0   0   0
80           Off (Field 1)  2   6  11  18
81 
82   rank 1: Dof (Field 0)
83           Dof (Field 1)
84 
85   rank 2: Dof (Field 0) 11
86           Off (Field 0) 18
87           Dof (Field 1)  2
88           Off (Field 1) 29
89 */
90 
91 typedef struct {
92   char      fname[PETSC_MAX_PATH_LEN]; /* Output mesh filename */
93   PetscBool includes_constraints;      /* Flag for if global section is to include constrained DoFs or not */
94 } AppCtx;
95 
96 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
97 {
98   PetscErrorCode  ierr;
99 
100   PetscFunctionBegin;
101   options->fname[0] = '\0';
102   options->includes_constraints = PETSC_TRUE;
103   ierr = PetscOptionsBegin(comm, "", "PetscSectionView()/Load() in HDF5 Test Options", "DMPLEX");PetscCall(ierr);
104   PetscCall(PetscOptionsString("-fname", "The output file", "ex5.c", options->fname, options->fname, sizeof(options->fname), NULL));
105   PetscCall(PetscOptionsBool("-includes_constraints", "Flag for if global section is to include constrained DoFs or not", "ex5.c", options->includes_constraints, &options->includes_constraints, NULL));
106   ierr = PetscOptionsEnd();PetscCall(ierr);
107   PetscFunctionReturn(0);
108 }
109 
110 int main(int argc, char **argv)
111 {
112   MPI_Comm        comm;
113   PetscMPIInt     size, rank, mycolor;
114   AppCtx          user;
115 
116   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
117   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
118   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
119   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
120   PetscCheck(size >= 3,PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes");
121 
122   /* Save */
123   mycolor = (PetscMPIInt)(rank >= 2);
124   PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
125   if (mycolor == 0) {
126     PetscSection  section, gsection;
127     PetscSF       sf;
128     PetscInt      nroots = -1, nleaves = -1, *ilocal;
129     PetscSFNode  *iremote;
130     PetscViewer   viewer;
131 
132     /* Create section */
133     PetscCall(PetscSectionCreate(comm, &section));
134     PetscCall(PetscSectionSetNumFields(section, 2));
135     switch (rank) {
136     case 0:
137       PetscCall(PetscSectionSetChart(section, 0, 4));
138       PetscCall(PetscSectionSetDof(section, 0, 3));
139       PetscCall(PetscSectionSetDof(section, 1, 3));
140       PetscCall(PetscSectionSetDof(section, 2, 5));
141       PetscCall(PetscSectionSetDof(section, 3, 7));
142       PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2));
143       PetscCall(PetscSectionSetFieldDof(section, 1, 0, 3));
144       PetscCall(PetscSectionSetFieldDof(section, 2, 0, 5));
145       PetscCall(PetscSectionSetFieldDof(section, 3, 0, 7));
146       PetscCall(PetscSectionSetFieldDof(section, 0, 1, 1));
147       break;
148     case 1:
149       PetscCall(PetscSectionSetChart(section, 0, 3));
150       PetscCall(PetscSectionSetDof(section, 0, 7));
151       PetscCall(PetscSectionSetDof(section, 1, 5));
152       PetscCall(PetscSectionSetDof(section, 2, 13));
153       PetscCall(PetscSectionSetConstraintDof(section, 2, 1));
154       PetscCall(PetscSectionSetFieldDof(section, 0, 0, 7));
155       PetscCall(PetscSectionSetFieldDof(section, 1, 0, 5));
156       PetscCall(PetscSectionSetFieldDof(section, 2, 0, 11));
157       PetscCall(PetscSectionSetFieldDof(section, 2, 1, 2));
158       PetscCall(PetscSectionSetFieldConstraintDof(section, 2, 0, 1));
159       break;
160     }
161     PetscCall(PetscSectionSetUp(section));
162     if (rank == 1)
163     {
164       const PetscInt indices[] = {7};
165       const PetscInt indices0[] = {7};
166 
167       PetscCall(PetscSectionSetConstraintIndices(section, 2, indices));
168       PetscCall(PetscSectionSetFieldConstraintIndices(section, 2, 0, indices0));
169     }
170     /* Create sf */
171     switch (rank) {
172     case 0:
173       nroots = 4;
174       nleaves = 1;
175       PetscCall(PetscMalloc1(nleaves, &ilocal));
176       PetscCall(PetscMalloc1(nleaves, &iremote));
177       ilocal[0] = 3;
178       iremote[0].rank = 1;
179       iremote[0].index = 0;
180       break;
181     case 1:
182       nroots = 3;
183       nleaves = 1;
184       PetscCall(PetscMalloc1(nleaves, &ilocal));
185       PetscCall(PetscMalloc1(nleaves, &iremote));
186       ilocal[0] = 1;
187       iremote[0].rank = 0;
188       iremote[0].index = 2;
189       break;
190     }
191     PetscCall(PetscSFCreate(comm, &sf));
192     PetscCall(PetscSFSetGraph(sf, nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
193     /* Create global section*/
194     PetscCall(PetscSectionCreateGlobalSection(section, sf, user.includes_constraints, PETSC_FALSE, &gsection));
195     PetscCall(PetscSFDestroy(&sf));
196     /* View */
197     PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer));
198     PetscCall(PetscSectionView(gsection, viewer));
199     PetscCall(PetscViewerDestroy(&viewer));
200     PetscCall(PetscObjectSetName((PetscObject)section, "Save: local section"));
201     PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
202     PetscCall(PetscObjectSetName((PetscObject)gsection, "Save: global section"));
203     PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm)));
204     PetscCall(PetscSectionDestroy(&gsection));
205     PetscCall(PetscSectionDestroy(&section));
206   }
207   PetscCallMPI(MPI_Comm_free(&comm));
208 
209   /* Load */
210   mycolor = (PetscMPIInt)(rank >= 3);
211   PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
212   if (mycolor == 0) {
213     PetscSection  section;
214     PetscInt      chartSize = -1;
215     PetscViewer   viewer;
216 
217     PetscCall(PetscSectionCreate(comm, &section));
218     switch (rank) {
219     case 0:
220       chartSize = 4;
221       break;
222     case 1:
223       chartSize = 0;
224       break;
225     case 2:
226       chartSize = 1;
227       break;
228     }
229     PetscCall(PetscSectionSetChart(section, 0, chartSize));
230     PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer));
231     PetscCall(PetscSectionLoad(section, viewer));
232     PetscCall(PetscViewerDestroy(&viewer));
233     PetscCall(PetscObjectSetName((PetscObject)section, "Load: section"));
234     PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
235     PetscCall(PetscSectionDestroy(&section));
236   }
237   PetscCallMPI(MPI_Comm_free(&comm));
238 
239   /* Finalize */
240   PetscCall(PetscFinalize());
241   return 0;
242 }
243 
244 /*TEST
245 
246   build:
247     requires: hdf5
248     requires: !complex
249   testset:
250     nsize: 4
251     test:
252       suffix: 0
253       args: -fname ex5_dump.h5 -includes_constraints 0
254     test:
255       suffix: 1
256       args: -fname ex5_dump.h5 -includes_constraints 1
257 
258 TEST*/
259