1 static char help[] = "Spectral element access patterns with Plex\n\n";
2
3 #include <petscdmplex.h>
4
5 typedef struct {
6 PetscInt Nf; /* Number of fields */
7 PetscInt *Nc; /* Number of components per field */
8 PetscInt *k; /* Spectral order per field */
9 } AppCtx;
10
ProcessOptions(MPI_Comm comm,AppCtx * options)11 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12 {
13 PetscInt len;
14 PetscBool flg;
15
16 PetscFunctionBeginUser;
17 options->Nf = 0;
18 options->Nc = NULL;
19 options->k = NULL;
20
21 PetscOptionsBegin(comm, "", "SEM Problem Options", "DMPLEX");
22 PetscCall(PetscOptionsBoundedInt("-num_fields", "The number of fields", "ex6.c", options->Nf, &options->Nf, NULL, 0));
23 if (options->Nf) {
24 len = options->Nf;
25 PetscCall(PetscMalloc1(len, &options->Nc));
26 PetscCall(PetscOptionsIntArray("-num_components", "The number of components per field", "ex6.c", options->Nc, &len, &flg));
27 PetscCheck(!flg || !(len != options->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of components array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->Nf);
28 len = options->Nf;
29 PetscCall(PetscMalloc1(len, &options->k));
30 PetscCall(PetscOptionsIntArray("-order", "The spectral order per field", "ex6.c", options->k, &len, &flg));
31 PetscCheck(!flg || !(len != options->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of order array is %" PetscInt_FMT " should be %" PetscInt_FMT, len, options->Nf);
32 }
33 PetscOptionsEnd();
34 PetscFunctionReturn(PETSC_SUCCESS);
35 }
36
LoadData2D(DM dm,PetscInt Ni,PetscInt Nj,PetscInt clSize,Vec u,AppCtx * user)37 static PetscErrorCode LoadData2D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt clSize, Vec u, AppCtx *user)
38 {
39 PetscInt i, j, f, c;
40 PetscScalar *closure;
41
42 PetscFunctionBeginUser;
43 PetscCall(PetscMalloc1(clSize, &closure));
44 for (j = 0; j < Nj; ++j) {
45 for (i = 0; i < Ni; ++i) {
46 PetscInt ki, kj, o = 0;
47 PetscCall(PetscArrayzero(closure, clSize));
48
49 for (f = 0; f < user->Nf; ++f) {
50 PetscInt ioff = i * user->k[f], joff = j * user->k[f];
51
52 for (kj = 0; kj <= user->k[f]; ++kj) {
53 for (ki = 0; ki <= user->k[f]; ++ki) {
54 for (c = 0; c < user->Nc[f]; ++c) closure[o++] = ((kj + joff) * (Ni * user->k[f] + 1) + ki + ioff) * user->Nc[f] + c;
55 }
56 }
57 }
58 PetscCall(DMPlexVecSetClosure(dm, NULL, u, j * Ni + i, closure, INSERT_VALUES));
59 }
60 }
61 PetscCall(PetscFree(closure));
62 PetscFunctionReturn(PETSC_SUCCESS);
63 }
64
LoadData3D(DM dm,PetscInt Ni,PetscInt Nj,PetscInt Nk,PetscInt clSize,Vec u,AppCtx * user)65 static PetscErrorCode LoadData3D(DM dm, PetscInt Ni, PetscInt Nj, PetscInt Nk, PetscInt clSize, Vec u, AppCtx *user)
66 {
67 PetscInt i, j, k, f, c;
68 PetscScalar *closure;
69
70 PetscFunctionBeginUser;
71 PetscCall(PetscMalloc1(clSize, &closure));
72 for (k = 0; k < Nk; ++k) {
73 for (j = 0; j < Nj; ++j) {
74 for (i = 0; i < Ni; ++i) {
75 PetscInt ki, kj, kk, o = 0;
76 PetscCall(PetscArrayzero(closure, clSize));
77
78 for (f = 0; f < user->Nf; ++f) {
79 PetscInt ioff = i * user->k[f], joff = j * user->k[f], koff = k * user->k[f];
80
81 for (kk = 0; kk <= user->k[f]; ++kk) {
82 for (kj = 0; kj <= user->k[f]; ++kj) {
83 for (ki = 0; ki <= user->k[f]; ++ki) {
84 for (c = 0; c < user->Nc[f]; ++c) closure[o++] = (((kk + koff) * (Nj * user->k[f] + 1) + kj + joff) * (Ni * user->k[f] + 1) + ki + ioff) * user->Nc[f] + c;
85 }
86 }
87 }
88 }
89 PetscCall(DMPlexVecSetClosure(dm, NULL, u, (k * Nj + j) * Ni + i, closure, INSERT_VALUES));
90 }
91 }
92 }
93 PetscCall(PetscFree(closure));
94 PetscFunctionReturn(PETSC_SUCCESS);
95 }
96
CheckPoint(DM dm,Vec u,PetscInt point,AppCtx * user)97 static PetscErrorCode CheckPoint(DM dm, Vec u, PetscInt point, AppCtx *user)
98 {
99 PetscSection s;
100 PetscScalar *a;
101 const PetscScalar *array;
102 PetscInt dof, d;
103
104 PetscFunctionBeginUser;
105 PetscCall(DMGetLocalSection(dm, &s));
106 PetscCall(VecGetArrayRead(u, &array));
107 PetscCall(DMPlexPointLocalRead(dm, point, array, &a));
108 PetscCall(PetscSectionGetDof(s, point, &dof));
109 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Point %" PetscInt_FMT ": ", point));
110 for (d = 0; d < dof; ++d) {
111 if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
112 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(a[d])));
113 }
114 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
115 PetscCall(VecRestoreArrayRead(u, &array));
116 PetscFunctionReturn(PETSC_SUCCESS);
117 }
118
ReadData2D(DM dm,Vec u,AppCtx * user)119 static PetscErrorCode ReadData2D(DM dm, Vec u, AppCtx *user)
120 {
121 PetscInt cStart, cEnd, cell;
122
123 PetscFunctionBeginUser;
124 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
125 for (cell = cStart; cell < cEnd; ++cell) {
126 PetscScalar *closure = NULL;
127 PetscInt closureSize, ki, kj, f, c, foff = 0;
128
129 PetscCall(DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure));
130 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell));
131 for (f = 0; f < user->Nf; ++f) {
132 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT "\n", f));
133 for (kj = user->k[f]; kj >= 0; --kj) {
134 for (ki = 0; ki <= user->k[f]; ++ki) {
135 if (ki > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
136 for (c = 0; c < user->Nc[f]; ++c) {
137 if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
138 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(closure[(kj * (user->k[f] + 1) + ki) * user->Nc[f] + c + foff])));
139 }
140 }
141 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
142 }
143 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
144 foff += PetscSqr(user->k[f] + 1);
145 }
146 PetscCall(DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure));
147 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
148 }
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151
ReadData3D(DM dm,Vec u,AppCtx * user)152 static PetscErrorCode ReadData3D(DM dm, Vec u, AppCtx *user)
153 {
154 PetscInt cStart, cEnd, cell;
155
156 PetscFunctionBeginUser;
157 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
158 for (cell = cStart; cell < cEnd; ++cell) {
159 PetscScalar *closure = NULL;
160 PetscInt closureSize, ki, kj, kk, f, c, foff = 0;
161
162 PetscCall(DMPlexVecGetClosure(dm, NULL, u, cell, &closureSize, &closure));
163 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT "\n", cell));
164 for (f = 0; f < user->Nf; ++f) {
165 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT "\n", f));
166 for (kk = user->k[f]; kk >= 0; --kk) {
167 for (kj = user->k[f]; kj >= 0; --kj) {
168 for (ki = 0; ki <= user->k[f]; ++ki) {
169 if (ki > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " "));
170 for (c = 0; c < user->Nc[f]; ++c) {
171 if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
172 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%2.0f", (double)PetscRealPart(closure[((kk * (user->k[f] + 1) + kj) * (user->k[f] + 1) + ki) * user->Nc[f] + c + foff])));
173 }
174 }
175 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
176 }
177 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
178 }
179 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
180 foff += PetscSqr(user->k[f] + 1);
181 }
182 PetscCall(DMPlexVecRestoreClosure(dm, NULL, u, cell, &closureSize, &closure));
183 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n\n"));
184 }
185 PetscFunctionReturn(PETSC_SUCCESS);
186 }
187
SetSymmetries(DM dm,PetscSection s,AppCtx * user)188 static PetscErrorCode SetSymmetries(DM dm, PetscSection s, AppCtx *user)
189 {
190 PetscInt dim, f, o, i, j, k, c, d;
191 DMLabel depthLabel;
192
193 PetscFunctionBegin;
194 PetscCall(DMGetDimension(dm, &dim));
195 PetscCall(DMGetLabel(dm, "depth", &depthLabel));
196 for (f = 0; f < user->Nf; f++) {
197 PetscSectionSym sym;
198
199 if (user->k[f] < 3) continue; /* No symmetries needed for order < 3, because no cell, facet, edge or vertex has more than one node */
200 PetscCall(PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)s), depthLabel, &sym));
201
202 for (d = 0; d <= dim; d++) {
203 if (d == 1) {
204 PetscInt numDof = user->k[f] - 1;
205 PetscInt numComp = user->Nc[f];
206 PetscInt minOrnt = -1;
207 PetscInt maxOrnt = 1;
208 PetscInt **perms;
209
210 PetscCall(PetscCalloc1(maxOrnt - minOrnt, &perms));
211 for (o = minOrnt; o < maxOrnt; o++) {
212 PetscInt *perm;
213
214 if (!o) { /* identity */
215 perms[o - minOrnt] = NULL;
216 } else {
217 PetscCall(PetscMalloc1(numDof * numComp, &perm));
218 for (i = numDof - 1, k = 0; i >= 0; i--) {
219 for (j = 0; j < numComp; j++, k++) perm[k] = i * numComp + j;
220 }
221 perms[o - minOrnt] = perm;
222 }
223 }
224 PetscCall(PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL));
225 } else if (d == 2) {
226 PetscInt perEdge = user->k[f] - 1;
227 PetscInt numDof = perEdge * perEdge;
228 PetscInt numComp = user->Nc[f];
229 PetscInt minOrnt = -4;
230 PetscInt maxOrnt = 4;
231 PetscInt **perms;
232
233 PetscCall(PetscCalloc1(maxOrnt - minOrnt, &perms));
234 for (o = minOrnt; o < maxOrnt; o++) {
235 PetscInt *perm;
236
237 if (!o) continue; /* identity */
238 PetscCall(PetscMalloc1(numDof * numComp, &perm));
239 /* We want to perm[k] to list which *localArray* position the *sectionArray* position k should go to for the given orientation*/
240 switch (o) {
241 case 0:
242 break; /* identity */
243 case -2: /* flip along (-1,-1)--( 1, 1), which swaps edges 0 and 3 and edges 1 and 2. This swaps the i and j variables */
244 for (i = 0, k = 0; i < perEdge; i++) {
245 for (j = 0; j < perEdge; j++, k++) {
246 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + i) * numComp + c;
247 }
248 }
249 break;
250 case -1: /* flip along (-1, 0)--( 1, 0), which swaps edges 0 and 2. This reverses the i variable */
251 for (i = 0, k = 0; i < perEdge; i++) {
252 for (j = 0; j < perEdge; j++, k++) {
253 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + j) * numComp + c;
254 }
255 }
256 break;
257 case -4: /* flip along ( 1,-1)--(-1, 1), which swaps edges 0 and 1 and edges 2 and 3. This swaps the i and j variables and reverse both */
258 for (i = 0, k = 0; i < perEdge; i++) {
259 for (j = 0; j < perEdge; j++, k++) {
260 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + (perEdge - 1 - i)) * numComp + c;
261 }
262 }
263 break;
264 case -3: /* flip along ( 0,-1)--( 0, 1), which swaps edges 3 and 1. This reverses the j variable */
265 for (i = 0, k = 0; i < perEdge; i++) {
266 for (j = 0; j < perEdge; j++, k++) {
267 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * i + (perEdge - 1 - j)) * numComp + c;
268 }
269 }
270 break;
271 case 1: /* rotate section edge 1 to local edge 0. This swaps the i and j variables and then reverses the j variable */
272 for (i = 0, k = 0; i < perEdge; i++) {
273 for (j = 0; j < perEdge; j++, k++) {
274 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - j) + i) * numComp + c;
275 }
276 }
277 break;
278 case 2: /* rotate section edge 2 to local edge 0. This reverse both i and j variables */
279 for (i = 0, k = 0; i < perEdge; i++) {
280 for (j = 0; j < perEdge; j++, k++) {
281 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * (perEdge - 1 - i) + (perEdge - 1 - j)) * numComp + c;
282 }
283 }
284 break;
285 case 3: /* rotate section edge 3 to local edge 0. This swaps the i and j variables and then reverses the i variable */
286 for (i = 0, k = 0; i < perEdge; i++) {
287 for (j = 0; j < perEdge; j++, k++) {
288 for (c = 0; c < numComp; c++) perm[k * numComp + c] = (perEdge * j + (perEdge - 1 - i)) * numComp + c;
289 }
290 }
291 break;
292 default:
293 break;
294 }
295 perms[o - minOrnt] = perm;
296 }
297 PetscCall(PetscSectionSymLabelSetStratum(sym, d, numDof * numComp, minOrnt, maxOrnt, PETSC_OWN_POINTER, (const PetscInt **)perms, NULL));
298 }
299 }
300 PetscCall(PetscSectionSetFieldSym(s, f, sym));
301 PetscCall(PetscSectionSymDestroy(&sym));
302 }
303 PetscCall(PetscSectionViewFromOptions(s, NULL, "-section_with_sym_view"));
304 PetscFunctionReturn(PETSC_SUCCESS);
305 }
306
main(int argc,char ** argv)307 int main(int argc, char **argv)
308 {
309 DM dm;
310 PetscSection s;
311 Vec u;
312 AppCtx user;
313 PetscInt dim, size = 0, f;
314
315 PetscFunctionBeginUser;
316 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
317 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
318 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
319 PetscCall(DMSetType(dm, DMPLEX));
320 PetscCall(DMSetFromOptions(dm));
321 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
322 PetscCall(DMGetDimension(dm, &dim));
323 /* Create a section for SEM order k */
324 {
325 PetscInt *numDof, d;
326
327 PetscCall(PetscMalloc1(user.Nf * (dim + 1), &numDof));
328 for (f = 0; f < user.Nf; ++f) {
329 for (d = 0; d <= dim; ++d) numDof[f * (dim + 1) + d] = PetscPowInt(user.k[f] - 1, d) * user.Nc[f];
330 size += PetscPowInt(user.k[f] + 1, d) * user.Nc[f];
331 }
332 PetscCall(DMSetNumFields(dm, user.Nf));
333 PetscCall(DMPlexCreateSection(dm, NULL, user.Nc, numDof, 0, NULL, NULL, NULL, NULL, &s));
334 PetscCall(SetSymmetries(dm, s, &user));
335 PetscCall(PetscFree(numDof));
336 }
337 PetscCall(DMSetLocalSection(dm, s));
338 /* Create spectral ordering and load in data */
339 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
340 PetscCall(DMGetLocalVector(dm, &u));
341 switch (dim) {
342 case 2:
343 PetscCall(LoadData2D(dm, 2, 2, size, u, &user));
344 break;
345 case 3:
346 PetscCall(LoadData3D(dm, 2, 2, 2, size, u, &user));
347 break;
348 }
349 /* Remove ordering and check some values */
350 PetscCall(PetscSectionSetClosurePermutation(s, (PetscObject)dm, dim, NULL));
351 switch (dim) {
352 case 2:
353 PetscCall(CheckPoint(dm, u, 0, &user));
354 PetscCall(CheckPoint(dm, u, 13, &user));
355 PetscCall(CheckPoint(dm, u, 15, &user));
356 PetscCall(CheckPoint(dm, u, 19, &user));
357 break;
358 case 3:
359 PetscCall(CheckPoint(dm, u, 0, &user));
360 PetscCall(CheckPoint(dm, u, 13, &user));
361 PetscCall(CheckPoint(dm, u, 15, &user));
362 PetscCall(CheckPoint(dm, u, 19, &user));
363 break;
364 }
365 /* Recreate spectral ordering and read out data */
366 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, s));
367 switch (dim) {
368 case 2:
369 PetscCall(ReadData2D(dm, u, &user));
370 break;
371 case 3:
372 PetscCall(ReadData3D(dm, u, &user));
373 break;
374 }
375 PetscCall(DMRestoreLocalVector(dm, &u));
376 PetscCall(PetscSectionDestroy(&s));
377 PetscCall(DMDestroy(&dm));
378 PetscCall(PetscFree(user.Nc));
379 PetscCall(PetscFree(user.k));
380 PetscCall(PetscFinalize());
381 return 0;
382 }
383
384 /*TEST
385
386 # Spectral ordering 2D 0-5
387 testset:
388 args: -dm_plex_simplex 0 -dm_plex_box_faces 2,2
389
390 test:
391 suffix: 0
392 args: -num_fields 1 -num_components 1 -order 2
393 test:
394 suffix: 1
395 args: -num_fields 1 -num_components 1 -order 3
396 test:
397 suffix: 2
398 args: -num_fields 1 -num_components 1 -order 5
399 test:
400 suffix: 3
401 args: -num_fields 1 -num_components 2 -order 2
402 test:
403 suffix: 4
404 args: -num_fields 2 -num_components 1,1 -order 2,2
405 test:
406 suffix: 5
407 args: -num_fields 2 -num_components 1,2 -order 2,3
408
409 # Spectral ordering 3D 6-11
410 testset:
411 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,2
412
413 test:
414 suffix: 6
415 args: -num_fields 1 -num_components 1 -order 2
416 test:
417 suffix: 7
418 args: -num_fields 1 -num_components 1 -order 3
419 test:
420 suffix: 8
421 args: -num_fields 1 -num_components 1 -order 5
422 test:
423 suffix: 9
424 args: -num_fields 1 -num_components 2 -order 2
425 test:
426 suffix: 10
427 args: -num_fields 2 -num_components 1,1 -order 2,2
428 test:
429 suffix: 11
430 args: -num_fields 2 -num_components 1,2 -order 2,3
431
432 TEST*/
433