Actual source code: ex9.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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:     DMPlexSetLabelValue(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:     DMPlexSetLabelValue(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:     DMPlexSetLabelValue(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:     DMPlexSetLabelValue(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:   PetscInt           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:   PetscInt           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:   PetscInt           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:   PetscLogBegin();
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: }