Actual source code: ex9.c
petsc-3.7.7 2017-09-25
1: static char help[] = "Performance tests for DMPlex query operations\n\n";
3: #include <petscdmplex.h>
5: typedef struct {
6: PetscInt dim; /* The topological mesh dimension */
7: PetscBool cellSimplex; /* Flag for simplices */
8: PetscBool interpolate; /* Flag for mesh interpolation */
9: PetscReal refinementLimit; /* Maximum volume of a refined cell */
10: PetscInt numFields; /* The number of section fields */
11: PetscInt *numComponents; /* The number of field components */
12: PetscInt *numDof; /* The dof signature for the section */
13: PetscBool reuseArray; /* Pass in user allocated array to VecGetClosure() */
14: /* Test data */
15: PetscBool errors; /* Treat failures as errors */
16: PetscInt iterations; /* The number of iterations for a query */
17: PetscReal maxConeTime; /* Max time per run for DMPlexGetCone() */
18: PetscReal maxClosureTime; /* Max time per run for DMPlexGetTransitiveClosure() */
19: PetscReal maxVecClosureTime; /* Max time per run for DMPlexVecGetClosure() */
20: } AppCtx;
24: static PetscErrorCode ProcessOptions(AppCtx *options)
25: {
26: PetscInt len;
27: PetscBool flg;
31: options->dim = 2;
32: options->cellSimplex = PETSC_TRUE;
33: options->interpolate = PETSC_FALSE;
34: options->refinementLimit = 0.0;
35: options->numFields = 0;
36: options->numComponents = NULL;
37: options->numDof = NULL;
38: options->reuseArray = PETSC_FALSE;
39: options->errors = PETSC_FALSE;
40: options->iterations = 1;
41: options->maxConeTime = 0.0;
42: options->maxClosureTime = 0.0;
43: options->maxVecClosureTime = 0.0;
45: PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");
46: PetscOptionsInt("-dim", "The topological mesh dimension", "ex9.c", options->dim, &options->dim, NULL);
47: PetscOptionsBool("-cellSimplex", "Flag for simplices", "ex9.c", options->cellSimplex, &options->cellSimplex, NULL);
48: PetscOptionsBool("-interpolate", "Flag for mesh interpolation", "ex9.c", options->interpolate, &options->interpolate, NULL);
49: PetscOptionsReal("-refinement_limit", "The maximum volume of a refined cell", "ex9.c", options->refinementLimit, &options->refinementLimit, NULL);
50: PetscOptionsInt("-num_fields", "The number of section fields", "ex9.c", options->numFields, &options->numFields, NULL);
51: if (options->numFields) {
52: len = options->numFields;
53: PetscMalloc1(len, &options->numComponents);
54: PetscOptionsIntArray("-num_components", "The number of components per field", "ex9.c", options->numComponents, &len, &flg);
55: if (flg && (len != options->numFields)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of components array is %d should be %d", len, options->numFields);
56: }
57: len = (options->dim+1) * PetscMax(1, options->numFields);
58: PetscMalloc1(len, &options->numDof);
59: PetscOptionsIntArray("-num_dof", "The dof signature for the section", "ex9.c", options->numDof, &len, &flg);
60: if (flg && (len != (options->dim+1) * PetscMax(1, options->numFields))) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of dof array is %d should be %d", len, (options->dim+1) * PetscMax(1, options->numFields));
62: PetscOptionsBool("-reuse_array", "Pass in user allocated array to VecGetClosure()", "ex9.c", options->reuseArray, &options->reuseArray, NULL);
63: PetscOptionsBool("-errors", "Treat failures as errors", "ex9.c", options->errors, &options->errors, NULL);
64: PetscOptionsInt("-iterations", "The number of iterations for a query", "ex9.c", options->iterations, &options->iterations, NULL);
65: PetscOptionsReal("-max_cone_time", "The maximum time per run for DMPlexGetCone()", "ex9.c", options->maxConeTime, &options->maxConeTime, NULL);
66: PetscOptionsReal("-max_closure_time", "The maximum time per run for DMPlexGetTransitiveClosure()", "ex9.c", options->maxClosureTime, &options->maxClosureTime, NULL);
67: PetscOptionsReal("-max_vec_closure_time", "The maximum time per run for DMPlexVecGetClosure()", "ex9.c", options->maxVecClosureTime, &options->maxVecClosureTime, NULL);
68: PetscOptionsEnd();
69: return(0);
70: };
74: static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, DM *newdm)
75: {
76: DM dm;
77: PetscInt numPoints[2] = {4, 2};
78: PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0};
79: PetscInt cones[6] = {2, 3, 4, 5, 4, 3};
80: PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
81: PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
82: PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1};
83: PetscInt dim = 2, depth = 1, p;
87: DMCreate(comm, &dm);
88: PetscObjectSetName((PetscObject) dm, "triangular");
89: DMSetType(dm, DMPLEX);
90: DMSetDimension(dm, dim);
91: DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);
92: for (p = 0; p < 4; ++p) {
93: DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);
94: }
95: *newdm = dm;
96: return(0);
97: }
101: static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, DM *newdm)
102: {
103: DM dm;
104: PetscInt numPoints[2] = {5, 2};
105: PetscInt coneSize[23] = {4, 4, 0, 0, 0, 0, 0};
106: PetscInt cones[8] = {2, 4, 3, 5, 3, 4, 6, 5};
107: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
108: PetscScalar vertexCoords[15] = {0.0, 0.0, -0.5, 0.0, -0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5};
109: PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
110: PetscInt dim = 3, depth = 1, p;
114: DMCreate(comm, &dm);
115: PetscObjectSetName((PetscObject) dm, "tetrahedral");
116: DMSetType(dm, DMPLEX);
117: DMSetDimension(dm, dim);
118: DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);
119: for (p = 0; p < 5; ++p) {
120: DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);
121: }
122: *newdm = dm;
123: return(0);
124: }
128: static PetscErrorCode CreateQuad_2D(MPI_Comm comm, DM *newdm)
129: {
130: DM dm;
131: PetscInt numPoints[2] = {6, 2};
132: PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0};
133: PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4};
134: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
135: PetscScalar vertexCoords[12] = {-0.5, 0.0, 0.0, 0.0, 0.0, 1.0, -0.5, 1.0, 0.5, 0.0, 0.5, 1.0};
136: PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1};
137: PetscInt dim = 2, depth = 1, p;
141: DMCreate(comm, &dm);
142: PetscObjectSetName((PetscObject) dm, "quadrilateral");
143: DMSetType(dm, DMPLEX);
144: DMSetDimension(dm, dim);
145: DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);
146: for (p = 0; p < 6; ++p) {
147: DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);
148: }
149: *newdm = dm;
150: return(0);
151: }
155: static PetscErrorCode CreateHex_3D(MPI_Comm comm, DM *newdm)
156: {
157: DM dm;
158: PetscInt numPoints[2] = {12, 2};
159: PetscInt coneSize[14] = {8, 8, 0,0,0,0,0,0,0,0,0,0,0,0};
160: PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8};
161: PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0, 0,0, 0, 0,0};
162: PetscScalar vertexCoords[36] = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0,
163: -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0,
164: 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0};
165: PetscInt markerPoints[24] = {2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1,10,1,11,1,12,1,13,1};
166: PetscInt dim = 3, depth = 1, p;
170: DMCreate(comm, &dm);
171: PetscObjectSetName((PetscObject) dm, "hexahedral");
172: DMSetType(dm, DMPLEX);
173: DMSetDimension(dm, dim);
174: DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);
175: for(p = 0; p < 12; ++p) {
176: DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);
177: }
178: *newdm = dm;
179: return(0);
180: }
184: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *newdm)
185: {
186: PetscInt dim = user->dim;
187: PetscBool cellSimplex = user->cellSimplex;
191: switch (dim) {
192: case 2:
193: if (cellSimplex) {
194: CreateSimplex_2D(comm, newdm);
195: } else {
196: CreateQuad_2D(comm, newdm);
197: }
198: break;
199: case 3:
200: if (cellSimplex) {
201: CreateSimplex_3D(comm, newdm);
202: } else {
203: CreateHex_3D(comm, newdm);
204: }
205: break;
206: default:
207: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
208: }
209: if (user->refinementLimit > 0.0) {
210: DM rdm;
211: const char *name;
213: DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);
214: DMPlexSetRefinementLimit(*newdm, user->refinementLimit);
215: DMRefine(*newdm, PETSC_COMM_SELF, &rdm);
216: PetscObjectGetName((PetscObject) *newdm, &name);
217: PetscObjectSetName((PetscObject) rdm, name);
218: DMDestroy(newdm);
219: *newdm = rdm;
220: }
221: if (user->interpolate) {
222: DM idm = NULL;
223: const char *name;
225: DMPlexInterpolate(*newdm, &idm);
226: PetscObjectGetName((PetscObject) *newdm, &name);
227: PetscObjectSetName((PetscObject) idm, name);
228: DMPlexCopyCoordinates(*newdm, idm);
229: DMDestroy(newdm);
230: *newdm = idm;
231: }
232: return(0);
233: }
237: static PetscErrorCode TestCone(DM dm, AppCtx *user)
238: {
239: PetscInt numRuns, cStart, cEnd, c, i;
240: PetscReal maxTimePerRun = user->maxConeTime;
241: PetscLogStage stage;
242: PetscLogEvent event;
243: PetscEventPerfInfo eventInfo;
244: PetscErrorCode ierr;
247: PetscLogStageRegister("DMPlex Cone Test", &stage);
248: PetscLogEventRegister("Cone", PETSC_OBJECT_CLASSID, &event);
249: PetscLogStagePush(stage);
250: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
251: PetscLogEventBegin(event,0,0,0,0);
252: for (i = 0; i < user->iterations; ++i) {
253: for (c = cStart; c < cEnd; ++c) {
254: const PetscInt *cone;
256: DMPlexGetCone(dm, c, &cone);
257: }
258: }
259: PetscLogEventEnd(event,0,0,0,0);
260: PetscLogStagePop();
262: PetscLogEventGetPerfInfo(stage, event, &eventInfo);
263: numRuns = (cEnd-cStart) * user->iterations;
264: if (eventInfo.count != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be %d", eventInfo.count, 1);
265: if ((PetscInt) eventInfo.flops != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %d should be %d", (PetscInt) eventInfo.flops, 0);
266: if (eventInfo.time > maxTimePerRun * numRuns) {
267: PetscPrintf(PETSC_COMM_SELF, "Cones: %d Average time per cone: %gs standard: %gs\n", numRuns, eventInfo.time/numRuns, maxTimePerRun);
268: if (user->errors) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for cone %g > standard %g", eventInfo.time/numRuns, maxTimePerRun);
269: }
270: return(0);
271: }
275: static PetscErrorCode TestTransitiveClosure(DM dm, AppCtx *user)
276: {
277: PetscInt numRuns, cStart, cEnd, c, i;
278: PetscReal maxTimePerRun = user->maxClosureTime;
279: PetscLogStage stage;
280: PetscLogEvent event;
281: PetscEventPerfInfo eventInfo;
282: PetscErrorCode ierr;
285: PetscLogStageRegister("DMPlex Transitive Closure Test", &stage);
286: PetscLogEventRegister("TransitiveClosure", PETSC_OBJECT_CLASSID, &event);
287: PetscLogStagePush(stage);
288: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
289: PetscLogEventBegin(event,0,0,0,0);
290: for (i = 0; i < user->iterations; ++i) {
291: for (c = cStart; c < cEnd; ++c) {
292: PetscInt *closure = NULL;
293: PetscInt closureSize;
295: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
296: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
297: }
298: }
299: PetscLogEventEnd(event,0,0,0,0);
300: PetscLogStagePop();
302: PetscLogEventGetPerfInfo(stage, event, &eventInfo);
303: numRuns = (cEnd-cStart) * user->iterations;
304: if (eventInfo.count != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be %d", eventInfo.count, 1);
305: if ((PetscInt) eventInfo.flops != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %d should be %d", (PetscInt) eventInfo.flops, 0);
306: if (eventInfo.time > maxTimePerRun * numRuns) {
307: PetscPrintf(PETSC_COMM_SELF, "Closures: %d Average time per cone: %gs standard: %gs\n", numRuns, eventInfo.time/numRuns, maxTimePerRun);
308: if (user->errors) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for closure %g > standard %g", eventInfo.time/numRuns, maxTimePerRun);
309: }
310: return(0);
311: }
315: static PetscErrorCode TestVecClosure(DM dm, PetscBool useIndex, AppCtx *user)
316: {
317: PetscSection s;
318: Vec v;
319: PetscInt numRuns, cStart, cEnd, c, i;
320: PetscScalar tmpArray[64];
321: PetscScalar *userArray = user->reuseArray ? tmpArray : NULL;
322: PetscReal maxTimePerRun = user->maxVecClosureTime;
323: PetscLogStage stage;
324: PetscLogEvent event;
325: PetscEventPerfInfo eventInfo;
326: PetscErrorCode ierr;
329: if (useIndex) {
330: PetscLogStageRegister("DMPlex Vector Closure with Index Test", &stage);
331: PetscLogEventRegister("VecClosureInd", PETSC_OBJECT_CLASSID, &event);
332: } else {
333: PetscLogStageRegister("DMPlex Vector Closure Test", &stage);
334: PetscLogEventRegister("VecClosure", PETSC_OBJECT_CLASSID, &event);
335: }
336: PetscLogStagePush(stage);
337: DMPlexCreateSection(dm, user->dim, user->numFields, user->numComponents, user->numDof, 0, NULL, NULL, NULL, NULL, &s);
338: DMSetDefaultSection(dm, s);
339: if (useIndex) {DMPlexCreateClosureIndex(dm, s);}
340: PetscSectionDestroy(&s);
341: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
342: DMGetLocalVector(dm, &v);
343: PetscLogEventBegin(event,0,0,0,0);
344: for (i = 0; i < user->iterations; ++i) {
345: for (c = cStart; c < cEnd; ++c) {
346: PetscScalar *closure = userArray;
347: PetscInt closureSize = 64;;
349: DMPlexVecGetClosure(dm, s, v, c, &closureSize, &closure);
350: if (!user->reuseArray) {DMPlexVecRestoreClosure(dm, s, v, c, &closureSize, &closure);}
351: }
352: }
353: PetscLogEventEnd(event,0,0,0,0);
354: DMRestoreLocalVector(dm, &v);
355: PetscLogStagePop();
357: PetscLogEventGetPerfInfo(stage, event, &eventInfo);
358: numRuns = (cEnd-cStart) * user->iterations;
359: if (eventInfo.count != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be %d", eventInfo.count, 1);
360: if ((PetscInt) eventInfo.flops != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %d should be %d", (PetscInt) eventInfo.flops, 0);
361: if (eventInfo.time > maxTimePerRun * numRuns) {
362: PetscPrintf(PETSC_COMM_SELF, "VecClosures: %d Average time per cone: %gs standard: %gs\n", numRuns, eventInfo.time/numRuns, maxTimePerRun);
363: if (user->errors) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for vector closure %g > standard %g", eventInfo.time/numRuns, maxTimePerRun);
364: }
365: return(0);
366: }
370: static PetscErrorCode CleanupContext(AppCtx *user)
371: {
375: PetscFree(user->numComponents);
376: PetscFree(user->numDof);
377: return(0);
378: }
382: int main(int argc, char **argv)
383: {
384: DM dm;
385: AppCtx user;
388: PetscInitialize(&argc, &argv, NULL, help);
389: ProcessOptions(&user);
390: PetscLogDefaultBegin();
391: CreateMesh(PETSC_COMM_SELF, &user, &dm);
392: TestCone(dm, &user);
393: TestTransitiveClosure(dm, &user);
394: TestVecClosure(dm, PETSC_FALSE, &user);
395: TestVecClosure(dm, PETSC_TRUE, &user);
396: DMDestroy(&dm);
397: CleanupContext(&user);
398: PetscFinalize();
399: return 0;
400: }