Actual source code: ex6.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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*/