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