1 static char help[] = "Performance tests for DMPlex query operations\n\n"; 2 3 #include <petscdmplex.h> 4 5 typedef struct { 6 PetscInt dim; /* The topological mesh dimension */ 7 PetscBool cellSimplex; /* Flag for simplices */ 8 PetscBool spectral; /* Flag for spectral element layout */ 9 PetscBool interpolate; /* Flag for mesh interpolation */ 10 PetscReal refinementLimit; /* Maximum volume of a refined cell */ 11 PetscInt numFields; /* The number of section fields */ 12 PetscInt *numComponents; /* The number of field components */ 13 PetscInt *numDof; /* The dof signature for the section */ 14 PetscBool reuseArray; /* Pass in user allocated array to VecGetClosure() */ 15 /* Test data */ 16 PetscBool errors; /* Treat failures as errors */ 17 PetscInt iterations; /* The number of iterations for a query */ 18 PetscReal maxConeTime; /* Max time per run for DMPlexGetCone() */ 19 PetscReal maxClosureTime; /* Max time per run for DMPlexGetTransitiveClosure() */ 20 PetscReal maxVecClosureTime; /* Max time per run for DMPlexVecGetClosure() */ 21 PetscBool printTimes; /* Print total times, do not check limits */ 22 } AppCtx; 23 24 static PetscErrorCode ProcessOptions(AppCtx *options) 25 { 26 PetscInt len; 27 PetscBool flg; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 options->dim = 2; 32 options->cellSimplex = PETSC_TRUE; 33 options->spectral = PETSC_FALSE; 34 options->interpolate = PETSC_FALSE; 35 options->refinementLimit = 0.0; 36 options->numFields = 0; 37 options->numComponents = NULL; 38 options->numDof = NULL; 39 options->reuseArray = PETSC_FALSE; 40 options->errors = PETSC_FALSE; 41 options->iterations = 1; 42 options->maxConeTime = 0.0; 43 options->maxClosureTime = 0.0; 44 options->maxVecClosureTime = 0.0; 45 options->printTimes = PETSC_FALSE; 46 47 ierr = PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr); 48 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex9.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 49 ierr = PetscOptionsBool("-cellSimplex", "Flag for simplices", "ex9.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr); 50 ierr = PetscOptionsBool("-spectral", "Flag for spectral element layout", "ex9.c", options->spectral, &options->spectral, NULL);CHKERRQ(ierr); 51 ierr = PetscOptionsBool("-interpolate", "Flag for mesh interpolation", "ex9.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr); 52 ierr = PetscOptionsReal("-refinement_limit", "The maximum volume of a refined cell", "ex9.c", options->refinementLimit, &options->refinementLimit, NULL);CHKERRQ(ierr); 53 ierr = PetscOptionsBoundedInt("-num_fields", "The number of section fields", "ex9.c", options->numFields, &options->numFields, NULL, 0);CHKERRQ(ierr); 54 if (options->numFields) { 55 len = options->numFields; 56 ierr = PetscMalloc1(len, &options->numComponents);CHKERRQ(ierr); 57 ierr = PetscOptionsIntArray("-num_components", "The number of components per field", "ex9.c", options->numComponents, &len, &flg);CHKERRQ(ierr); 58 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); 59 } 60 len = (options->dim+1) * PetscMax(1, options->numFields); 61 ierr = PetscMalloc1(len, &options->numDof);CHKERRQ(ierr); 62 ierr = PetscOptionsIntArray("-num_dof", "The dof signature for the section", "ex9.c", options->numDof, &len, &flg);CHKERRQ(ierr); 63 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)); 64 65 /* We are specifying the scalar dof, so augment it for multiple components */ 66 { 67 PetscInt f, d; 68 69 for (f = 0; f < options->numFields; ++f) { 70 for (d = 0; d <= options->dim; ++d) options->numDof[f*(options->dim+1)+d] *= options->numComponents[f]; 71 } 72 } 73 74 ierr = PetscOptionsBool("-reuse_array", "Pass in user allocated array to VecGetClosure()", "ex9.c", options->reuseArray, &options->reuseArray, NULL);CHKERRQ(ierr); 75 ierr = PetscOptionsBool("-errors", "Treat failures as errors", "ex9.c", options->errors, &options->errors, NULL);CHKERRQ(ierr); 76 ierr = PetscOptionsBoundedInt("-iterations", "The number of iterations for a query", "ex9.c", options->iterations, &options->iterations, NULL,0);CHKERRQ(ierr); 77 ierr = PetscOptionsReal("-max_cone_time", "The maximum time per run for DMPlexGetCone()", "ex9.c", options->maxConeTime, &options->maxConeTime, NULL);CHKERRQ(ierr); 78 ierr = PetscOptionsReal("-max_closure_time", "The maximum time per run for DMPlexGetTransitiveClosure()", "ex9.c", options->maxClosureTime, &options->maxClosureTime, NULL);CHKERRQ(ierr); 79 ierr = PetscOptionsReal("-max_vec_closure_time", "The maximum time per run for DMPlexVecGetClosure()", "ex9.c", options->maxVecClosureTime, &options->maxVecClosureTime, NULL);CHKERRQ(ierr); 80 ierr = PetscOptionsBool("-print_times", "Print total times, do not check limits", "ex9.c", options->printTimes, &options->printTimes, NULL);CHKERRQ(ierr); 81 ierr = PetscOptionsEnd();CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 static PetscErrorCode CreateSimplex_2D(MPI_Comm comm, DM *newdm) 86 { 87 DM dm; 88 PetscInt numPoints[2] = {4, 2}; 89 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0}; 90 PetscInt cones[6] = {2, 3, 4, 5, 4, 3}; 91 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0}; 92 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5}; 93 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1}; 94 PetscInt dim = 2, depth = 1, p; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 99 ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr); 100 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 101 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 102 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 103 for (p = 0; p < 4; ++p) { 104 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 105 } 106 *newdm = dm; 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode CreateSimplex_3D(MPI_Comm comm, DM *newdm) 111 { 112 DM dm; 113 PetscInt numPoints[2] = {5, 2}; 114 PetscInt coneSize[23] = {4, 4, 0, 0, 0, 0, 0}; 115 PetscInt cones[8] = {2, 4, 3, 5, 3, 4, 6, 5}; 116 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 117 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}; 118 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1}; 119 PetscInt dim = 3, depth = 1, p; 120 PetscErrorCode ierr; 121 122 PetscFunctionBegin; 123 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 124 ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr); 125 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 126 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 127 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 128 for (p = 0; p < 5; ++p) { 129 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 130 } 131 *newdm = dm; 132 PetscFunctionReturn(0); 133 } 134 135 static PetscErrorCode CreateQuad_2D(MPI_Comm comm, DM *newdm) 136 { 137 DM dm; 138 PetscInt numPoints[2] = {6, 2}; 139 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0}; 140 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4}; 141 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0}; 142 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}; 143 PetscInt markerPoints[12] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1, 7, 1}; 144 PetscInt dim = 2, depth = 1, p; 145 PetscErrorCode ierr; 146 147 PetscFunctionBegin; 148 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 149 ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr); 150 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 151 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 152 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 153 for (p = 0; p < 6; ++p) { 154 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 155 } 156 *newdm = dm; 157 PetscFunctionReturn(0); 158 } 159 160 static PetscErrorCode CreateHex_3D(MPI_Comm comm, DM *newdm) 161 { 162 DM dm; 163 PetscInt numPoints[2] = {12, 2}; 164 PetscInt coneSize[14] = {8, 8, 0,0,0,0,0,0,0,0,0,0,0,0}; 165 PetscInt cones[16] = {2,5,4,3,6,7,8,9, 3,4,11,10,7,12,13,8}; 166 PetscInt coneOrientations[16] = {0,0,0,0,0,0,0,0, 0,0, 0, 0,0, 0, 0,0}; 167 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, 168 -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0, 169 0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0, 0.5,1.0,1.0}; 170 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}; 171 PetscInt dim = 3, depth = 1, p; 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 ierr = DMCreate(comm, &dm);CHKERRQ(ierr); 176 ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr); 177 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 178 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); 179 ierr = DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 180 for (p = 0; p < 12; ++p) { 181 ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr); 182 } 183 *newdm = dm; 184 PetscFunctionReturn(0); 185 } 186 187 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *newdm) 188 { 189 PetscInt dim = user->dim; 190 PetscBool cellSimplex = user->cellSimplex; 191 PetscErrorCode ierr; 192 193 PetscFunctionBegin; 194 switch (dim) { 195 case 2: 196 if (cellSimplex) { 197 ierr = CreateSimplex_2D(comm, newdm);CHKERRQ(ierr); 198 } else { 199 ierr = CreateQuad_2D(comm, newdm);CHKERRQ(ierr); 200 } 201 break; 202 case 3: 203 if (cellSimplex) { 204 ierr = CreateSimplex_3D(comm, newdm);CHKERRQ(ierr); 205 } else { 206 ierr = CreateHex_3D(comm, newdm);CHKERRQ(ierr); 207 } 208 break; 209 default: 210 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim); 211 } 212 if (user->refinementLimit > 0.0) { 213 DM rdm; 214 const char *name; 215 216 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr); 217 ierr = DMPlexSetRefinementLimit(*newdm, user->refinementLimit);CHKERRQ(ierr); 218 ierr = DMRefine(*newdm, PETSC_COMM_SELF, &rdm);CHKERRQ(ierr); 219 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr); 220 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr); 221 ierr = DMDestroy(newdm);CHKERRQ(ierr); 222 *newdm = rdm; 223 } 224 if (user->interpolate) { 225 DM idm; 226 227 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr); 228 ierr = DMDestroy(newdm);CHKERRQ(ierr); 229 *newdm = idm; 230 } 231 ierr = DMSetFromOptions(*newdm);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 static PetscErrorCode TestCone(DM dm, AppCtx *user) 236 { 237 PetscInt numRuns, cStart, cEnd, c, i; 238 PetscReal maxTimePerRun = user->maxConeTime; 239 PetscLogStage stage; 240 PetscLogEvent event; 241 PetscEventPerfInfo eventInfo; 242 MPI_Comm comm; 243 PetscMPIInt rank; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 248 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 249 ierr = PetscLogStageRegister("DMPlex Cone Test", &stage);CHKERRQ(ierr); 250 ierr = PetscLogEventRegister("Cone", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 251 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 252 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 253 ierr = PetscLogEventBegin(event,0,0,0,0);CHKERRQ(ierr); 254 for (i = 0; i < user->iterations; ++i) { 255 for (c = cStart; c < cEnd; ++c) { 256 const PetscInt *cone; 257 258 ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr); 259 } 260 } 261 ierr = PetscLogEventEnd(event,0,0,0,0);CHKERRQ(ierr); 262 ierr = PetscLogStagePop();CHKERRQ(ierr); 263 264 ierr = PetscLogEventGetPerfInfo(stage, event, &eventInfo);CHKERRQ(ierr); 265 numRuns = (cEnd-cStart) * user->iterations; 266 if (eventInfo.count != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be %d", eventInfo.count, 1); 267 if ((PetscInt) eventInfo.flops != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %d should be %d", (PetscInt) eventInfo.flops, 0); 268 if (user->printTimes) { 269 ierr = PetscSynchronizedPrintf(comm, "[%d] Cones: %d Total time: %.3es Average time per cone: %.3es\n", rank, numRuns, eventInfo.time, eventInfo.time/numRuns);CHKERRQ(ierr); 270 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 271 } else if (eventInfo.time > maxTimePerRun * numRuns) { 272 ierr = PetscSynchronizedPrintf(comm, "[%d] Cones: %d Average time per cone: %gs standard: %gs\n", rank, numRuns, eventInfo.time/numRuns, maxTimePerRun);CHKERRQ(ierr); 273 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 274 if (user->errors) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for cone %g > standard %g", eventInfo.time/numRuns, maxTimePerRun); 275 } 276 PetscFunctionReturn(0); 277 } 278 279 static PetscErrorCode TestTransitiveClosure(DM dm, AppCtx *user) 280 { 281 PetscInt numRuns, cStart, cEnd, c, i; 282 PetscReal maxTimePerRun = user->maxClosureTime; 283 PetscLogStage stage; 284 PetscLogEvent event; 285 PetscEventPerfInfo eventInfo; 286 MPI_Comm comm; 287 PetscMPIInt rank; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 292 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 293 ierr = PetscLogStageRegister("DMPlex Transitive Closure Test", &stage);CHKERRQ(ierr); 294 ierr = PetscLogEventRegister("TransitiveClosure", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 295 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 296 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 297 ierr = PetscLogEventBegin(event,0,0,0,0);CHKERRQ(ierr); 298 for (i = 0; i < user->iterations; ++i) { 299 for (c = cStart; c < cEnd; ++c) { 300 PetscInt *closure = NULL; 301 PetscInt closureSize; 302 303 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 304 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 305 } 306 } 307 ierr = PetscLogEventEnd(event,0,0,0,0);CHKERRQ(ierr); 308 ierr = PetscLogStagePop();CHKERRQ(ierr); 309 310 ierr = PetscLogEventGetPerfInfo(stage, event, &eventInfo);CHKERRQ(ierr); 311 numRuns = (cEnd-cStart) * user->iterations; 312 if (eventInfo.count != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be %d", eventInfo.count, 1); 313 if ((PetscInt) eventInfo.flops != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %d should be %d", (PetscInt) eventInfo.flops, 0); 314 if (user->printTimes) { 315 ierr = PetscSynchronizedPrintf(comm, "[%d] Closures: %d Total time: %.3es Average time per cone: %.3es\n", rank, numRuns, eventInfo.time, eventInfo.time/numRuns);CHKERRQ(ierr); 316 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 317 } else if (eventInfo.time > maxTimePerRun * numRuns) { 318 ierr = PetscSynchronizedPrintf(comm, "[%d] Closures: %d Average time per cone: %gs standard: %gs\n", rank, numRuns, eventInfo.time/numRuns, maxTimePerRun);CHKERRQ(ierr); 319 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 320 if (user->errors) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for closure %g > standard %g", eventInfo.time/numRuns, maxTimePerRun); 321 } 322 PetscFunctionReturn(0); 323 } 324 325 static PetscErrorCode TestVecClosure(DM dm, PetscBool useIndex, PetscBool useSpectral, AppCtx *user) 326 { 327 PetscSection s; 328 Vec v; 329 PetscInt numRuns, cStart, cEnd, c, i; 330 PetscScalar tmpArray[64]; 331 PetscScalar *userArray = user->reuseArray ? tmpArray : NULL; 332 PetscReal maxTimePerRun = user->maxVecClosureTime; 333 PetscLogStage stage; 334 PetscLogEvent event; 335 PetscEventPerfInfo eventInfo; 336 MPI_Comm comm; 337 PetscMPIInt rank; 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 342 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 343 if (useIndex) { 344 if (useSpectral) { 345 ierr = PetscLogStageRegister("DMPlex Vector Closure with Index Test", &stage);CHKERRQ(ierr); 346 ierr = PetscLogEventRegister("VecClosureInd", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 347 } else { 348 ierr = PetscLogStageRegister("DMPlex Vector Spectral Closure with Index Test", &stage);CHKERRQ(ierr); 349 ierr = PetscLogEventRegister("VecClosureSpecInd", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 350 } 351 } else { 352 if (useSpectral) { 353 ierr = PetscLogStageRegister("DMPlex Vector Spectral Closure Test", &stage);CHKERRQ(ierr); 354 ierr = PetscLogEventRegister("VecClosureSpec", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 355 } else { 356 ierr = PetscLogStageRegister("DMPlex Vector Closure Test", &stage);CHKERRQ(ierr); 357 ierr = PetscLogEventRegister("VecClosure", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr); 358 } 359 } 360 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 361 ierr = DMSetNumFields(dm, user->numFields);CHKERRQ(ierr); 362 ierr = DMPlexCreateSection(dm, NULL, user->numComponents, user->numDof, 0, NULL, NULL, NULL, NULL, &s);CHKERRQ(ierr); 363 ierr = DMSetLocalSection(dm, s);CHKERRQ(ierr); 364 if (useIndex) {ierr = DMPlexCreateClosureIndex(dm, s);CHKERRQ(ierr);} 365 if (useSpectral) {ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, s);CHKERRQ(ierr);} 366 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); 367 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 368 ierr = DMGetLocalVector(dm, &v);CHKERRQ(ierr); 369 ierr = PetscLogEventBegin(event,0,0,0,0);CHKERRQ(ierr); 370 for (i = 0; i < user->iterations; ++i) { 371 for (c = cStart; c < cEnd; ++c) { 372 PetscScalar *closure = userArray; 373 PetscInt closureSize = 64; 374 375 ierr = DMPlexVecGetClosure(dm, s, v, c, &closureSize, &closure);CHKERRQ(ierr); 376 if (!user->reuseArray) {ierr = DMPlexVecRestoreClosure(dm, s, v, c, &closureSize, &closure);CHKERRQ(ierr);} 377 } 378 } 379 ierr = PetscLogEventEnd(event,0,0,0,0);CHKERRQ(ierr); 380 ierr = DMRestoreLocalVector(dm, &v);CHKERRQ(ierr); 381 ierr = PetscLogStagePop();CHKERRQ(ierr); 382 383 ierr = PetscLogEventGetPerfInfo(stage, event, &eventInfo);CHKERRQ(ierr); 384 numRuns = (cEnd-cStart) * user->iterations; 385 if (eventInfo.count != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event calls %d should be %d", eventInfo.count, 1); 386 if ((PetscInt) eventInfo.flops != 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of event flops %d should be %d", (PetscInt) eventInfo.flops, 0); 387 if (user->printTimes || eventInfo.time > maxTimePerRun * numRuns) { 388 const char *title = "VecClosures"; 389 const char *titleIndex = "VecClosures with Index"; 390 const char *titleSpec = "VecClosures Spectral"; 391 const char *titleSpecIndex = "VecClosures Spectral with Index"; 392 393 if (user->printTimes) { 394 ierr = PetscSynchronizedPrintf(comm, "[%d] %s: %d Total time: %.3es Average time per vector closure: %.3es\n", rank, useIndex ? (useSpectral ? titleSpecIndex : titleIndex) : (useSpectral ? titleSpec : title), numRuns, eventInfo.time, eventInfo.time/numRuns);CHKERRQ(ierr); 395 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 396 } else { 397 ierr = PetscSynchronizedPrintf(comm, "[%d] %s: %d Average time per vector closure: %gs standard: %gs\n", rank, useIndex ? (useSpectral ? titleSpecIndex : titleIndex) : (useSpectral ? titleSpec : title), numRuns, eventInfo.time/numRuns, maxTimePerRun);CHKERRQ(ierr); 398 ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT);CHKERRQ(ierr); 399 if (user->errors) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Average time for vector closure %g > standard %g", eventInfo.time/numRuns, maxTimePerRun); 400 } 401 } 402 PetscFunctionReturn(0); 403 } 404 405 static PetscErrorCode CleanupContext(AppCtx *user) 406 { 407 PetscErrorCode ierr; 408 409 PetscFunctionBegin; 410 ierr = PetscFree(user->numComponents);CHKERRQ(ierr); 411 ierr = PetscFree(user->numDof);CHKERRQ(ierr); 412 PetscFunctionReturn(0); 413 } 414 415 int main(int argc, char **argv) 416 { 417 DM dm; 418 AppCtx user; 419 PetscErrorCode ierr; 420 421 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 422 ierr = ProcessOptions(&user);CHKERRQ(ierr); 423 ierr = PetscLogDefaultBegin();CHKERRQ(ierr); 424 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 425 ierr = TestCone(dm, &user);CHKERRQ(ierr); 426 ierr = TestTransitiveClosure(dm, &user);CHKERRQ(ierr); 427 ierr = TestVecClosure(dm, PETSC_FALSE, PETSC_FALSE, &user);CHKERRQ(ierr); 428 ierr = TestVecClosure(dm, PETSC_TRUE, PETSC_FALSE, &user);CHKERRQ(ierr); 429 if (!user.cellSimplex && user.spectral) { 430 ierr = TestVecClosure(dm, PETSC_FALSE, PETSC_TRUE, &user);CHKERRQ(ierr); 431 ierr = TestVecClosure(dm, PETSC_TRUE, PETSC_TRUE, &user);CHKERRQ(ierr); 432 } 433 ierr = DMDestroy(&dm);CHKERRQ(ierr); 434 ierr = CleanupContext(&user);CHKERRQ(ierr); 435 ierr = PetscFinalize(); 436 return ierr; 437 } 438 439 /*TEST 440 441 build: 442 requires: define(PETSC_USE_LOG) 443 444 # 2D Simplex P_1 scalar tests 445 testset: 446 args: -num_dof 1,0,0 -iterations 2 -print_times 447 test: 448 suffix: correctness_0 449 test: 450 suffix: correctness_1 451 args: -interpolate -dm_refine 2 452 test: 453 suffix: correctness_2 454 requires: triangle 455 args: -interpolate -refinement_limit 1.0e-5 456 test: 457 suffix: 0 458 TODO: Only for performance testing 459 args: -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 1.3e-7 -max_vec_closure_time 3.6e-7 460 test: 461 suffix: 1 462 requires: triangle 463 TODO: Only for performance testing 464 args: -refinement_limit 1.0e-5 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 1.5e-7 -max_vec_closure_time 3.6e-7 465 test: 466 suffix: 2 467 TODO: Only for performance testing 468 args: -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 1.3e-7 -max_vec_closure_time 4.5e-7 469 test: 470 suffix: 3 471 requires: triangle 472 TODO: Only for performance testing 473 args: -refinement_limit 1.0e-5 -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 1.5e-7 -max_vec_closure_time 4.7e-7 474 test: 475 suffix: 4 476 TODO: Only for performance testing 477 args: -interpolate -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.0e-6 478 test: 479 suffix: 5 480 requires: triangle 481 TODO: Only for performance testing 482 args: -interpolate -refinement_limit 1.0e-4 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.0e-6 483 test: 484 suffix: 6 485 TODO: Only for performance testing 486 args: -interpolate -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 10000 -max_cone_time 1.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.1e-6 487 test: 488 suffix: 7 489 requires: triangle 490 TODO: Only for performance testing 491 args: -interpolate -refinement_limit 1.0e-4 -num_fields 1 -num_components 1 -num_dof 1,0,0 -iterations 2 -max_cone_time 2.1e-8 -max_closure_time 6.5e-7 -max_vec_closure_time 1.2e-6 492 493 # 2D Simplex P_1 vector tests 494 # 2D Simplex P_2 scalar tests 495 # 2D Simplex P_2 vector tests 496 # 2D Simplex P_2/P_1 vector/scalar tests 497 # 2D Quad P_1 scalar tests 498 # 2D Quad P_1 vector tests 499 # 2D Quad P_2 scalar tests 500 # 2D Quad P_2 vector tests 501 # 3D Simplex P_1 scalar tests 502 # 3D Simplex P_1 vector tests 503 # 3D Simplex P_2 scalar tests 504 # 3D Simplex P_2 vector tests 505 # 3D Hex P_1 scalar tests 506 # 3D Hex P_1 vector tests 507 # 3D Hex P_2 scalar tests 508 # 3D Hex P_2 vector tests 509 510 TEST*/ 511