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: {
 16:   PetscFunctionBeginUser;
 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:   PetscCall(PetscOptionsInt("-dim", "The topological dimension", "ex1.c", options->dim, &options->dim, NULL));
 24:   PetscCall(PetscOptionsBool("-simplex", "Simplex or hex cells", "ex1.c", options->simplex, &options->simplex, NULL));
 25:   PetscCall(PetscOptionsInt("-its", "The number of replications for timing", "ex1.c", options->its, &options->its, NULL));
 26:   PetscCall(PetscOptionsInt("-cbs", "The number of cells in an integration block", "ex1.c", options->cbs, &options->cbs, NULL));
 27:   PetscOptionsEnd();
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 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 PETSC_SUCCESS;
 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;

 63:   PetscFunctionBeginUser;
 64:   PetscCall(DMGetDS(dm, &prob));
 65:   PetscCall(PetscDSSetResidual(prob, 0, f0_trig_u, f1_u));
 66:   PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu));
 67:   PetscCall(PetscDSSetExactSolution(prob, 0, trig_u, user));
 68:   PetscCall(DMGetLabel(dm, "marker", &label));
 69:   PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))trig_u, NULL, user, NULL));
 70:   PetscFunctionReturn(PETSC_SUCCESS);
 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];

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

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

101:   PetscFunctionBegin;
102:   PetscCall(PetscFEGeomDestroy(&geom));
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

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

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

129: PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
130: {
131:   PetscFunctionBegin;
132:   *geom = NULL;
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

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

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

157:       PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
158:       PetscCall(PetscFEGetQuadrature(fe, &(*quads)[f]));
159:       PetscCall(PetscObjectReference((PetscObject)(*quads)[f]));
160:       PetscCall(CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]));
161:     }
162:   }
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

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

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

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

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

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

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

251:     PetscCall(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:     PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

259: static PetscErrorCode TestIntegration2(DM dm, PetscInt cbs, PetscInt its)
260: {
261:   Vec           X, F;
262:   PetscLogStage stage;
263:   PetscInt      i;

265:   PetscFunctionBeginUser;
266:   PetscCall(PetscLogStageRegister("DMPlex Residual Integration Test", &stage));
267:   PetscCall(PetscLogStagePush(stage));
268:   PetscCall(DMGetLocalVector(dm, &X));
269:   PetscCall(DMGetLocalVector(dm, &F));
270:   for (i = 0; i < its; ++i) PetscCall(DMPlexSNESComputeResidualFEM(dm, X, F, NULL));
271:   PetscCall(DMRestoreLocalVector(dm, &X));
272:   PetscCall(DMRestoreLocalVector(dm, &F));
273:   PetscCall(PetscLogStagePop());
274:   if (PetscDefined(USE_LOG)) {
275:     const char        *title = "DMPlex Residual Integration";
276:     PetscEventPerfInfo eventInfo;
277:     PetscReal          flopRate, cellRate;
278:     PetscInt           cStart, cEnd, Nf, N;
279:     PetscLogEvent      event;

281:     PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
282:     PetscCall(DMGetNumFields(dm, &Nf));
283:     PetscCall(PetscLogEventGetId("DMPlexResidualFE", &event));
284:     PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo));
285:     N        = (cEnd - cStart) * Nf * eventInfo.count;
286:     flopRate = eventInfo.time != 0.0 ? eventInfo.flops / eventInfo.time : 0.0;
287:     cellRate = eventInfo.time != 0.0 ? N / eventInfo.time : 0.0;
288:     PetscCall(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)));
289:   }
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

293: int main(int argc, char **argv)
294: {
295:   DM          dm;
296:   AppCtx      ctx;
297:   PetscMPIInt size;

299:   PetscFunctionBeginUser;
300:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
301:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
302:   PetscCheck(size <= 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only.");
303:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
304:   PetscCall(PetscLogDefaultBegin());
305:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
306:   PetscCall(DMSetType(dm, DMPLEX));
307:   PetscCall(DMSetFromOptions(dm));
308:   PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh"));
309:   PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_view"));
310:   PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &ctx));
311:   PetscCall(TestIntegration(dm, ctx.cbs, ctx.its));
312:   PetscCall(TestIntegration2(dm, ctx.cbs, ctx.its));
313:   PetscCall(DMDestroy(&dm));
314:   PetscCall(PetscFinalize());
315:   return 0;
316: }

318: /*TEST
319:   test:
320:     suffix: 0
321:     requires: triangle
322:     args: -dm_view

324:   test:
325:     suffix: 1
326:     requires: triangle
327:     args: -dm_view -potential_petscspace_degree 1

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

334:   test:
335:     suffix: 3
336:     requires: triangle
337:     args: -dm_view -potential_petscspace_degree 3
338: TEST*/