Actual source code: ex1.c

  1: static const char help[] = "Performance Tests for FE Integration";

  3: #include <petscdmplex.h>
  4: #include <petscfe.h>
  5: #include <petscds.h>

  7: typedef struct {
  8:   PetscInt  dim;     /* The topological dimension */
  9:   PetscBool simplex; /* True for simplices, false for hexes */
 10:   PetscInt  its;     /* Number of replications for timing */
 11:   PetscInt  cbs;     /* Number of cells in an integration block */
 12: } AppCtx;

 14: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 15: {
 17:   options->dim     = 2;
 18:   options->simplex = PETSC_TRUE;
 19:   options->its     = 1;
 20:   options->cbs     = 8;

 22:   PetscOptionsBegin(comm, "", "FE Integration Performance Options", "PETSCFE");
 23:   PetscOptionsInt("-dim", "The topological dimension", "ex1.c", options->dim, &options->dim, NULL);
 24:   PetscOptionsBool("-simplex", "Simplex or hex cells", "ex1.c", options->simplex, &options->simplex, NULL);
 25:   PetscOptionsInt("-its", "The number of replications for timing", "ex1.c", options->its, &options->its, NULL);
 26:   PetscOptionsInt("-cbs", "The number of cells in an integration block", "ex1.c", options->cbs, &options->cbs, NULL);
 27:   PetscOptionsEnd();
 28:   return 0;
 29: }

 31: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
 32: {
 33:   PetscInt d;
 34:   *u = 0.0;
 35:   for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]);
 36:   return 0;
 37: }

 39: static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
 40: {
 41:   PetscInt d;
 42:   for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
 43: }

 45: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
 46: {
 47:   PetscInt d;
 48:   for (d = 0; d < dim; ++d) f1[d] = u_x[d];
 49: }

 51: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
 52: {
 53:   PetscInt d;
 54:   for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
 55: }

 57: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
 58: {
 59:   PetscDS        prob;
 60:   DMLabel        label;
 61:   const PetscInt id = 1;

 64:   DMGetDS(dm, &prob);
 65:   PetscDSSetResidual(prob, 0, f0_trig_u, f1_u);
 66:   PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
 67:   PetscDSSetExactSolution(prob, 0, trig_u, user);
 68:   DMGetLabel(dm, "marker", &label);
 69:   DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL);
 70:   return 0;
 71: }

 73: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
 74: {
 75:   DM      cdm = dm;
 76:   PetscFE fe;
 77:   char    prefix[PETSC_MAX_PATH_LEN];

 80:   /* Create finite element */
 81:   PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
 82:   PetscFECreateDefault(PetscObjectComm((PetscObject)dm), user->dim, 1, user->simplex, name ? prefix : NULL, -1, &fe);
 83:   PetscObjectSetName((PetscObject)fe, name);
 84:   /* Set discretization and boundary conditions for each mesh */
 85:   DMSetField(dm, 0, NULL, (PetscObject)fe);
 86:   DMCreateDS(dm);
 87:   (*setup)(dm, user);
 88:   while (cdm) {
 89:     DMCopyDisc(dm, cdm);
 90:     /* TODO: Check whether the boundary of coarse meshes is marked */
 91:     DMGetCoarseDM(cdm, &cdm);
 92:   }
 93:   PetscFEDestroy(&fe);
 94:   return 0;
 95: }

 97: static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom(void *ctx)
 98: {
 99:   PetscFEGeom *geom = (PetscFEGeom *)ctx;

101:   PetscFEGeomDestroy(&geom);
102:   return 0;
103: }

105: PetscErrorCode CellRangeGetFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
106: {
107:   char           composeStr[33] = {0};
108:   PetscObjectId  id;
109:   PetscContainer container;

111:   PetscObjectGetId((PetscObject)quad, &id);
112:   PetscSNPrintf(composeStr, 32, "CellRangeGetFEGeom_%" PetscInt64_FMT "\n", id);
113:   PetscObjectQuery((PetscObject)cellIS, composeStr, (PetscObject *)&container);
114:   if (container) {
115:     PetscContainerGetPointer(container, (void **)geom);
116:   } else {
117:     DMFieldCreateFEGeom(coordField, cellIS, quad, faceData, geom);
118:     PetscContainerCreate(PETSC_COMM_SELF, &container);
119:     PetscContainerSetPointer(container, (void *)*geom);
120:     PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
121:     PetscObjectCompose((PetscObject)cellIS, composeStr, (PetscObject)container);
122:     PetscContainerDestroy(&container);
123:   }
124:   return 0;
125: }

127: PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
128: {
129:   *geom = NULL;
130:   return 0;
131: }

133: static PetscErrorCode CreateFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
134: {
135:   DMField  coordField;
136:   PetscInt Nf, f, maxDegree;

139:   *affineQuad = NULL;
140:   *affineGeom = NULL;
141:   *quads      = NULL;
142:   *geoms      = NULL;
143:   PetscDSGetNumFields(ds, &Nf);
144:   DMGetCoordinateField(dm, &coordField);
145:   DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
146:   if (maxDegree <= 1) {
147:     DMFieldCreateDefaultQuadrature(coordField, cellIS, affineQuad);
148:     if (*affineQuad) CellRangeGetFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom);
149:   } else {
150:     PetscCalloc2(Nf, quads, Nf, geoms);
151:     for (f = 0; f < Nf; ++f) {
152:       PetscFE fe;

154:       PetscDSGetDiscretization(ds, f, (PetscObject *)&fe);
155:       PetscFEGetQuadrature(fe, &(*quads)[f]);
156:       PetscObjectReference((PetscObject)(*quads)[f]);
157:       CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]);
158:     }
159:   }
160:   return 0;
161: }

163: static PetscErrorCode DestroyFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms)
164: {
165:   DMField  coordField;
166:   PetscInt Nf, f;

169:   PetscDSGetNumFields(ds, &Nf);
170:   DMGetCoordinateField(dm, &coordField);
171:   if (*affineQuad) {
172:     CellRangeRestoreFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom);
173:     PetscQuadratureDestroy(affineQuad);
174:   } else {
175:     for (f = 0; f < Nf; ++f) {
176:       CellRangeRestoreFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]);
177:       PetscQuadratureDestroy(&(*quads)[f]);
178:     }
179:     PetscFree2(*quads, *geoms);
180:   }
181:   return 0;
182: }

184: static PetscErrorCode TestIntegration(DM dm, PetscInt cbs, PetscInt its)
185: {
186:   PetscDS         ds;
187:   PetscFEGeom    *chunkGeom = NULL;
188:   PetscQuadrature affineQuad, *quads  = NULL;
189:   PetscFEGeom    *affineGeom, **geoms = NULL;
190:   PetscScalar    *u, *elemVec;
191:   IS              cellIS;
192:   PetscInt        depth, cStart, cEnd, cell, chunkSize = cbs, Nch = 0, Nf, f, totDim, i, k;
193: #if defined(PETSC_USE_LOG)
194:   PetscLogStage stage;
195:   PetscLogEvent event;
196: #endif

199:   PetscLogStageRegister("PetscFE Residual Integration Test", &stage);
200:   PetscLogEventRegister("FEIntegRes", PETSCFE_CLASSID, &event);
201:   PetscLogStagePush(stage);
202:   DMPlexGetDepth(dm, &depth);
203:   DMGetStratumIS(dm, "depth", depth, &cellIS);
204:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
205:   DMGetCellDS(dm, cStart, &ds);
206:   PetscDSGetNumFields(ds, &Nf);
207:   PetscDSGetTotalDimension(ds, &totDim);
208:   CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms);
209:   PetscMalloc2(chunkSize * totDim, &u, chunkSize * totDim, &elemVec);
210:   /* Assumptions:
211:     - Single field
212:     - No input data
213:     - No auxiliary data
214:     - No time-dependence
215:   */
216:   for (i = 0; i < its; ++i) {
217:     for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) {
218:       const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS;

220:       PetscArrayzero(elemVec, chunkSize * totDim);
221:       /* TODO Replace with DMPlexGetCellFields() */
222:       for (k = 0; k < chunkSize * totDim; ++k) u[k] = 1.0;
223:       for (f = 0; f < Nf; ++f) {
224:         PetscFormKey key;
225:         PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
226:         /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */

228:         key.label = NULL;
229:         key.value = 0;
230:         key.field = f;
231:         key.part  = 0;
232:         PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom);
233:         PetscLogEventBegin(event, 0, 0, 0, 0);
234:         PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec);
235:         PetscLogEventEnd(event, 0, 0, 0, 0);
236:       }
237:     }
238:   }
239:   PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom);
240:   DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms);
241:   ISDestroy(&cellIS);
242:   PetscFree2(u, elemVec);
243:   PetscLogStagePop();
244: #if defined(PETSC_USE_LOG)
245:   {
246:     const char        *title = "Petsc FE Residual Integration";
247:     PetscEventPerfInfo eventInfo;
248:     PetscInt           N = (cEnd - cStart) * Nf * its;
249:     PetscReal          flopRate, cellRate;

251:     PetscLogEventGetPerfInfo(stage, event, &eventInfo);
252:     flopRate = eventInfo.time != 0.0 ? eventInfo.flops / eventInfo.time : 0.0;
253:     cellRate = eventInfo.time != 0.0 ? N / eventInfo.time : 0.0;
254:     PetscPrintf(PetscObjectComm((PetscObject)dm), "%s: %" PetscInt_FMT " integrals %" PetscInt_FMT " chunks %" PetscInt_FMT " reps\n  Cell rate: %.2f/s flop rate: %.2f MF/s\n", title, N, Nch, its, (double)cellRate, (double)(flopRate / 1.e6));
255:   }
256: #endif
257:   return 0;
258: }

260: static PetscErrorCode TestIntegration2(DM dm, PetscInt cbs, PetscInt its)
261: {
262:   Vec X, F;
263: #if defined(PETSC_USE_LOG)
264:   PetscLogStage stage;
265: #endif
266:   PetscInt i;

269:   PetscLogStageRegister("DMPlex Residual Integration Test", &stage);
270:   PetscLogStagePush(stage);
271:   DMGetLocalVector(dm, &X);
272:   DMGetLocalVector(dm, &F);
273:   for (i = 0; i < its; ++i) DMPlexSNESComputeResidualFEM(dm, X, F, NULL);
274:   DMRestoreLocalVector(dm, &X);
275:   DMRestoreLocalVector(dm, &F);
276:   PetscLogStagePop();
277: #if defined(PETSC_USE_LOG)
278:   {
279:     const char        *title = "DMPlex Residual Integration";
280:     PetscEventPerfInfo eventInfo;
281:     PetscReal          flopRate, cellRate;
282:     PetscInt           cStart, cEnd, Nf, N;
283:     PetscLogEvent      event;

285:     DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
286:     DMGetNumFields(dm, &Nf);
287:     PetscLogEventGetId("DMPlexResidualFE", &event);
288:     PetscLogEventGetPerfInfo(stage, event, &eventInfo);
289:     N        = (cEnd - cStart) * Nf * eventInfo.count;
290:     flopRate = eventInfo.time != 0.0 ? eventInfo.flops / eventInfo.time : 0.0;
291:     cellRate = eventInfo.time != 0.0 ? N / eventInfo.time : 0.0;
292:     PetscPrintf(PetscObjectComm((PetscObject)dm), "%s: %" PetscInt_FMT " integrals %d reps\n  Cell rate: %.2f/s flop rate: %.2f MF/s\n", title, N, eventInfo.count, (double)cellRate, (double)(flopRate / 1.e6));
293:   }
294: #endif
295:   return 0;
296: }

298: int main(int argc, char **argv)
299: {
300:   DM          dm;
301:   AppCtx      ctx;
302:   PetscMPIInt size;

305:   PetscInitialize(&argc, &argv, NULL, help);
306:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
308:   ProcessOptions(PETSC_COMM_WORLD, &ctx);
309:   PetscLogDefaultBegin();
310:   DMCreate(PETSC_COMM_WORLD, &dm);
311:   DMSetType(dm, DMPLEX);
312:   DMSetFromOptions(dm);
313:   PetscObjectSetName((PetscObject)dm, "Mesh");
314:   PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_view");
315:   SetupDiscretization(dm, "potential", SetupPrimalProblem, &ctx);
316:   TestIntegration(dm, ctx.cbs, ctx.its);
317:   TestIntegration2(dm, ctx.cbs, ctx.its);
318:   DMDestroy(&dm);
319:   PetscFinalize();
320:   return 0;
321: }

323: /*TEST
324:   test:
325:     suffix: 0
326:     requires: triangle
327:     args: -dm_view

329:   test:
330:     suffix: 1
331:     requires: triangle
332:     args: -dm_view -potential_petscspace_degree 1

334:   test:
335:     suffix: 2
336:     requires: triangle
337:     args: -dm_view -potential_petscspace_degree 2

339:   test:
340:     suffix: 3
341:     requires: triangle
342:     args: -dm_view -potential_petscspace_degree 3
343: TEST*/