xref: /petsc/src/vec/is/tests/ex5.c (revision 76d901e46dda72c1afe96306c7cb4731c47d4e87)
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   PetscFunctionBegin;
99   options->fname[0] = '\0';
100   options->includes_constraints = PETSC_TRUE;
101   PetscOptionsBegin(comm, "", "PetscSectionView()/Load() in HDF5 Test Options", "DMPLEX");
102   PetscCall(PetscOptionsString("-fname", "The output file", "ex5.c", options->fname, options->fname, sizeof(options->fname), NULL));
103   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));
104   PetscOptionsEnd();
105   PetscFunctionReturn(0);
106 }
107 
108 int main(int argc, char **argv)
109 {
110   MPI_Comm        comm;
111   PetscMPIInt     size, rank, mycolor;
112   AppCtx          user;
113 
114   PetscFunctionBeginUser;
115   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
116   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
117   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
118   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
119   PetscCheck(size >= 3,PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example only works with three or more processes");
120 
121   /* Save */
122   mycolor = (PetscMPIInt)(rank >= 2);
123   PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
124   if (mycolor == 0) {
125     PetscSection  section, gsection;
126     PetscSF       sf;
127     PetscInt      nroots = -1, nleaves = -1, *ilocal;
128     PetscSFNode  *iremote;
129     PetscViewer   viewer;
130 
131     /* Create section */
132     PetscCall(PetscSectionCreate(comm, &section));
133     PetscCall(PetscSectionSetNumFields(section, 2));
134     switch (rank) {
135     case 0:
136       PetscCall(PetscSectionSetChart(section, 0, 4));
137       PetscCall(PetscSectionSetDof(section, 0, 3));
138       PetscCall(PetscSectionSetDof(section, 1, 3));
139       PetscCall(PetscSectionSetDof(section, 2, 5));
140       PetscCall(PetscSectionSetDof(section, 3, 7));
141       PetscCall(PetscSectionSetFieldDof(section, 0, 0, 2));
142       PetscCall(PetscSectionSetFieldDof(section, 1, 0, 3));
143       PetscCall(PetscSectionSetFieldDof(section, 2, 0, 5));
144       PetscCall(PetscSectionSetFieldDof(section, 3, 0, 7));
145       PetscCall(PetscSectionSetFieldDof(section, 0, 1, 1));
146       break;
147     case 1:
148       PetscCall(PetscSectionSetChart(section, 0, 3));
149       PetscCall(PetscSectionSetDof(section, 0, 7));
150       PetscCall(PetscSectionSetDof(section, 1, 5));
151       PetscCall(PetscSectionSetDof(section, 2, 13));
152       PetscCall(PetscSectionSetConstraintDof(section, 2, 1));
153       PetscCall(PetscSectionSetFieldDof(section, 0, 0, 7));
154       PetscCall(PetscSectionSetFieldDof(section, 1, 0, 5));
155       PetscCall(PetscSectionSetFieldDof(section, 2, 0, 11));
156       PetscCall(PetscSectionSetFieldDof(section, 2, 1, 2));
157       PetscCall(PetscSectionSetFieldConstraintDof(section, 2, 0, 1));
158       break;
159     }
160     PetscCall(PetscSectionSetUp(section));
161     if (rank == 1)
162     {
163       const PetscInt indices[] = {7};
164       const PetscInt indices0[] = {7};
165 
166       PetscCall(PetscSectionSetConstraintIndices(section, 2, indices));
167       PetscCall(PetscSectionSetFieldConstraintIndices(section, 2, 0, indices0));
168     }
169     /* Create sf */
170     switch (rank) {
171     case 0:
172       nroots = 4;
173       nleaves = 1;
174       PetscCall(PetscMalloc1(nleaves, &ilocal));
175       PetscCall(PetscMalloc1(nleaves, &iremote));
176       ilocal[0] = 3;
177       iremote[0].rank = 1;
178       iremote[0].index = 0;
179       break;
180     case 1:
181       nroots = 3;
182       nleaves = 1;
183       PetscCall(PetscMalloc1(nleaves, &ilocal));
184       PetscCall(PetscMalloc1(nleaves, &iremote));
185       ilocal[0] = 1;
186       iremote[0].rank = 0;
187       iremote[0].index = 2;
188       break;
189     }
190     PetscCall(PetscSFCreate(comm, &sf));
191     PetscCall(PetscSFSetGraph(sf, nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
192     /* Create global section*/
193     PetscCall(PetscSectionCreateGlobalSection(section, sf, user.includes_constraints, PETSC_FALSE, &gsection));
194     PetscCall(PetscSFDestroy(&sf));
195     /* View */
196     PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer));
197     PetscCall(PetscSectionView(gsection, viewer));
198     PetscCall(PetscViewerDestroy(&viewer));
199     PetscCall(PetscObjectSetName((PetscObject)section, "Save: local section"));
200     PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
201     PetscCall(PetscObjectSetName((PetscObject)gsection, "Save: global section"));
202     PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm)));
203     PetscCall(PetscSectionDestroy(&gsection));
204     PetscCall(PetscSectionDestroy(&section));
205   }
206   PetscCallMPI(MPI_Comm_free(&comm));
207 
208   /* Load */
209   mycolor = (PetscMPIInt)(rank >= 3);
210   PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
211   if (mycolor == 0) {
212     PetscSection  section;
213     PetscInt      chartSize = -1;
214     PetscViewer   viewer;
215 
216     PetscCall(PetscSectionCreate(comm, &section));
217     switch (rank) {
218     case 0:
219       chartSize = 4;
220       break;
221     case 1:
222       chartSize = 0;
223       break;
224     case 2:
225       chartSize = 1;
226       break;
227     }
228     PetscCall(PetscSectionSetChart(section, 0, chartSize));
229     PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer));
230     PetscCall(PetscSectionLoad(section, viewer));
231     PetscCall(PetscViewerDestroy(&viewer));
232     PetscCall(PetscObjectSetName((PetscObject)section, "Load: section"));
233     PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
234     PetscCall(PetscSectionDestroy(&section));
235   }
236   PetscCallMPI(MPI_Comm_free(&comm));
237 
238   /* Finalize */
239   PetscCall(PetscFinalize());
240   return 0;
241 }
242 
243 /*TEST
244 
245   build:
246     requires: hdf5
247     requires: !complex
248   testset:
249     nsize: 4
250     test:
251       suffix: 0
252       args: -fname ex5_dump.h5 -includes_constraints 0
253     test:
254       suffix: 1
255       args: -fname ex5_dump.h5 -includes_constraints 1
256 
257 TEST*/
258