xref: /petsc/src/vec/is/tests/ex5.c (revision 4ad8454beace47809662cdae21ee081016eaa39a)
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(PETSC_SUCCESS);
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       const PetscInt indices[]  = {7};
163       const PetscInt indices0[] = {7};
164 
165       PetscCall(PetscSectionSetConstraintIndices(section, 2, indices));
166       PetscCall(PetscSectionSetFieldConstraintIndices(section, 2, 0, indices0));
167     }
168     /* Create sf */
169     switch (rank) {
170     case 0:
171       nroots  = 4;
172       nleaves = 1;
173       PetscCall(PetscMalloc1(nleaves, &ilocal));
174       PetscCall(PetscMalloc1(nleaves, &iremote));
175       ilocal[0]        = 3;
176       iremote[0].rank  = 1;
177       iremote[0].index = 0;
178       break;
179     case 1:
180       nroots  = 3;
181       nleaves = 1;
182       PetscCall(PetscMalloc1(nleaves, &ilocal));
183       PetscCall(PetscMalloc1(nleaves, &iremote));
184       ilocal[0]        = 1;
185       iremote[0].rank  = 0;
186       iremote[0].index = 2;
187       break;
188     }
189     PetscCall(PetscSFCreate(comm, &sf));
190     PetscCall(PetscSFSetGraph(sf, nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
191     /* Create global section*/
192     PetscCall(PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, user.includes_constraints, PETSC_FALSE, &gsection));
193     PetscCall(PetscSFDestroy(&sf));
194     /* View */
195     PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_WRITE, &viewer));
196     PetscCall(PetscSectionView(gsection, viewer));
197     PetscCall(PetscViewerDestroy(&viewer));
198     PetscCall(PetscObjectSetName((PetscObject)section, "Save: local section"));
199     PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
200     PetscCall(PetscObjectSetName((PetscObject)gsection, "Save: global section"));
201     PetscCall(PetscSectionView(gsection, PETSC_VIEWER_STDOUT_(comm)));
202     PetscCall(PetscSectionDestroy(&gsection));
203     PetscCall(PetscSectionDestroy(&section));
204   }
205   PetscCallMPI(MPI_Comm_free(&comm));
206 
207   /* Load */
208   mycolor = (PetscMPIInt)(rank >= 3);
209   PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, rank, &comm));
210   if (mycolor == 0) {
211     PetscSection section;
212     PetscInt     chartSize = -1;
213     PetscViewer  viewer;
214 
215     PetscCall(PetscSectionCreate(comm, &section));
216     switch (rank) {
217     case 0:
218       chartSize = 4;
219       break;
220     case 1:
221       chartSize = 0;
222       break;
223     case 2:
224       chartSize = 1;
225       break;
226     }
227     PetscCall(PetscSectionSetChart(section, 0, chartSize));
228     PetscCall(PetscViewerHDF5Open(comm, user.fname, FILE_MODE_READ, &viewer));
229     PetscCall(PetscSectionLoad(section, viewer));
230     PetscCall(PetscViewerDestroy(&viewer));
231     PetscCall(PetscObjectSetName((PetscObject)section, "Load: section"));
232     PetscCall(PetscSectionView(section, PETSC_VIEWER_STDOUT_(comm)));
233     PetscCall(PetscSectionDestroy(&section));
234   }
235   PetscCallMPI(MPI_Comm_free(&comm));
236 
237   /* Finalize */
238   PetscCall(PetscFinalize());
239   return 0;
240 }
241 
242 /*TEST
243 
244   build:
245     requires: hdf5
246     requires: !complex
247   testset:
248     nsize: 4
249     test:
250       suffix: 0
251       args: -fname ex5_dump.h5 -includes_constraints 0
252     test:
253       suffix: 1
254       args: -fname ex5_dump.h5 -includes_constraints 1
255 
256 TEST*/
257