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