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*/