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: #if defined(PETSC_USE_LOG)
197: PetscLogStage stage;
198: PetscLogEvent event;
199: #endif
201: PetscFunctionBeginUser;
202: PetscCall(PetscLogStageRegister("PetscFE Residual Integration Test", &stage));
203: PetscCall(PetscLogEventRegister("FEIntegRes", PETSCFE_CLASSID, &event));
204: PetscCall(PetscLogStagePush(stage));
205: PetscCall(DMPlexGetDepth(dm, &depth));
206: PetscCall(DMGetStratumIS(dm, "depth", depth, &cellIS));
207: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
208: PetscCall(DMGetCellDS(dm, cStart, &ds, NULL));
209: PetscCall(PetscDSGetNumFields(ds, &Nf));
210: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
211: PetscCall(CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms));
212: PetscCall(PetscMalloc2(chunkSize * totDim, &u, chunkSize * totDim, &elemVec));
213: /* Assumptions:
214: - Single field
215: - No input data
216: - No auxiliary data
217: - No time-dependence
218: */
219: for (i = 0; i < its; ++i) {
220: for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) {
221: const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS;
223: PetscCall(PetscArrayzero(elemVec, chunkSize * totDim));
224: /* TODO Replace with DMPlexGetCellFields() */
225: for (k = 0; k < chunkSize * totDim; ++k) u[k] = 1.0;
226: for (f = 0; f < Nf; ++f) {
227: PetscFormKey key;
228: PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f];
229: /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */
231: key.label = NULL;
232: key.value = 0;
233: key.field = f;
234: key.part = 0;
235: PetscCall(PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom));
236: PetscCall(PetscLogEventBegin(event, 0, 0, 0, 0));
237: PetscCall(PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec));
238: PetscCall(PetscLogEventEnd(event, 0, 0, 0, 0));
239: }
240: }
241: }
242: PetscCall(PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom));
243: PetscCall(DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms));
244: PetscCall(ISDestroy(&cellIS));
245: PetscCall(PetscFree2(u, elemVec));
246: PetscCall(PetscLogStagePop());
247: #if defined(PETSC_USE_LOG)
248: {
249: const char *title = "Petsc FE Residual Integration";
250: PetscEventPerfInfo eventInfo;
251: PetscInt N = (cEnd - cStart) * Nf * its;
252: PetscReal flopRate, cellRate;
254: PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo));
255: flopRate = eventInfo.time != 0.0 ? eventInfo.flops / eventInfo.time : 0.0;
256: cellRate = eventInfo.time != 0.0 ? N / eventInfo.time : 0.0;
257: 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)));
258: }
259: #endif
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode TestIntegration2(DM dm, PetscInt cbs, PetscInt its)
264: {
265: Vec X, F;
266: #if defined(PETSC_USE_LOG)
267: PetscLogStage stage;
268: #endif
269: PetscInt i;
271: PetscFunctionBeginUser;
272: PetscCall(PetscLogStageRegister("DMPlex Residual Integration Test", &stage));
273: PetscCall(PetscLogStagePush(stage));
274: PetscCall(DMGetLocalVector(dm, &X));
275: PetscCall(DMGetLocalVector(dm, &F));
276: for (i = 0; i < its; ++i) PetscCall(DMPlexSNESComputeResidualFEM(dm, X, F, NULL));
277: PetscCall(DMRestoreLocalVector(dm, &X));
278: PetscCall(DMRestoreLocalVector(dm, &F));
279: PetscCall(PetscLogStagePop());
280: #if defined(PETSC_USE_LOG)
281: {
282: const char *title = "DMPlex Residual Integration";
283: PetscEventPerfInfo eventInfo;
284: PetscReal flopRate, cellRate;
285: PetscInt cStart, cEnd, Nf, N;
286: PetscLogEvent event;
288: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
289: PetscCall(DMGetNumFields(dm, &Nf));
290: PetscCall(PetscLogEventGetId("DMPlexResidualFE", &event));
291: PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo));
292: N = (cEnd - cStart) * Nf * eventInfo.count;
293: flopRate = eventInfo.time != 0.0 ? eventInfo.flops / eventInfo.time : 0.0;
294: cellRate = eventInfo.time != 0.0 ? N / eventInfo.time : 0.0;
295: 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)));
296: }
297: #endif
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: int main(int argc, char **argv)
302: {
303: DM dm;
304: AppCtx ctx;
305: PetscMPIInt size;
307: PetscFunctionBeginUser;
308: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
309: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
310: PetscCheck(size <= 1, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only.");
311: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx));
312: PetscCall(PetscLogDefaultBegin());
313: PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
314: PetscCall(DMSetType(dm, DMPLEX));
315: PetscCall(DMSetFromOptions(dm));
316: PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh"));
317: PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_view"));
318: PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &ctx));
319: PetscCall(TestIntegration(dm, ctx.cbs, ctx.its));
320: PetscCall(TestIntegration2(dm, ctx.cbs, ctx.its));
321: PetscCall(DMDestroy(&dm));
322: PetscCall(PetscFinalize());
323: return 0;
324: }
326: /*TEST
327: test:
328: suffix: 0
329: requires: triangle
330: args: -dm_view
332: test:
333: suffix: 1
334: requires: triangle
335: args: -dm_view -potential_petscspace_degree 1
337: test:
338: suffix: 2
339: requires: triangle
340: args: -dm_view -potential_petscspace_degree 2
342: test:
343: suffix: 3
344: requires: triangle
345: args: -dm_view -potential_petscspace_degree 3
346: TEST*/