1 static const char help[] = "Performance Tests for FE Integration"; 2 3 #include <petscdmplex.h> 4 #include <petscfe.h> 5 #include <petscds.h> 6 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; 13 14 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 15 { 16 PetscErrorCode ierr; 17 18 PetscFunctionBeginUser; 19 options->dim = 2; 20 options->simplex = PETSC_TRUE; 21 options->its = 1; 22 options->cbs = 8; 23 24 ierr = PetscOptionsBegin(comm, "", "FE Integration Performance Options", "PETSCFE");CHKERRQ(ierr); 25 ierr = PetscOptionsInt("-dim", "The topological dimension", "ex1.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 26 ierr = PetscOptionsBool("-simplex", "Simplex or hex cells", "ex1.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 27 ierr = PetscOptionsInt("-its", "The number of replications for timing", "ex1.c", options->its, &options->its, NULL);CHKERRQ(ierr); 28 ierr = PetscOptionsInt("-cbs", "The number of cells in an integration block", "ex1.c", options->cbs, &options->cbs, NULL);CHKERRQ(ierr); 29 ierr = PetscOptionsEnd();CHKERRQ(ierr); 30 PetscFunctionReturn(0); 31 } 32 33 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 34 { 35 PetscInt d; 36 *u = 0.0; 37 for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0*PETSC_PI*x[d]); 38 return 0; 39 } 40 41 static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 42 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 43 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 44 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 45 { 46 PetscInt d; 47 for (d = 0; d < dim; ++d) f0[0] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]); 48 } 49 50 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 51 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 52 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 53 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 54 { 55 PetscInt d; 56 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 57 } 58 59 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 60 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 61 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 62 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 63 { 64 PetscInt d; 65 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 66 } 67 68 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 69 { 70 PetscDS prob; 71 DMLabel label; 72 const PetscInt id = 1; 73 PetscErrorCode ierr; 74 75 PetscFunctionBeginUser; 76 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 77 ierr = PetscDSSetResidual(prob, 0, f0_trig_u, f1_u);CHKERRQ(ierr); 78 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 79 ierr = PetscDSSetExactSolution(prob, 0, trig_u, user);CHKERRQ(ierr); 80 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 81 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) trig_u, NULL, user, NULL);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 86 { 87 DM cdm = dm; 88 PetscFE fe; 89 char prefix[PETSC_MAX_PATH_LEN]; 90 PetscErrorCode ierr; 91 92 PetscFunctionBeginUser; 93 /* Create finite element */ 94 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 95 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), user->dim, 1, user->simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 96 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 97 /* Set discretization and boundary conditions for each mesh */ 98 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 99 ierr = DMCreateDS(dm);CHKERRQ(ierr); 100 ierr = (*setup)(dm, user);CHKERRQ(ierr); 101 while (cdm) { 102 ierr = DMCopyDisc(dm,cdm);CHKERRQ(ierr); 103 /* TODO: Check whether the boundary of coarse meshes is marked */ 104 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 105 } 106 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom(void *ctx) 111 { 112 PetscFEGeom *geom = (PetscFEGeom *) ctx; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = PetscFEGeomDestroy(&geom);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 PetscErrorCode CellRangeGetFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 121 { 122 char composeStr[33] = {0}; 123 PetscObjectId id; 124 PetscContainer container; 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ierr = PetscObjectGetId((PetscObject) quad, &id);CHKERRQ(ierr); 129 ierr = PetscSNPrintf(composeStr, 32, "CellRangeGetFEGeom_%x\n", id);CHKERRQ(ierr); 130 ierr = PetscObjectQuery((PetscObject) cellIS, composeStr, (PetscObject *) &container);CHKERRQ(ierr); 131 if (container) { 132 ierr = PetscContainerGetPointer(container, (void **) geom);CHKERRQ(ierr); 133 } else { 134 ierr = DMFieldCreateFEGeom(coordField, cellIS, quad, faceData, geom);CHKERRQ(ierr); 135 ierr = PetscContainerCreate(PETSC_COMM_SELF, &container);CHKERRQ(ierr); 136 ierr = PetscContainerSetPointer(container, (void *) *geom);CHKERRQ(ierr); 137 ierr = PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);CHKERRQ(ierr); 138 ierr = PetscObjectCompose((PetscObject) cellIS, composeStr, (PetscObject) container);CHKERRQ(ierr); 139 ierr = PetscContainerDestroy(&container);CHKERRQ(ierr); 140 } 141 PetscFunctionReturn(0); 142 } 143 144 PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 145 { 146 PetscFunctionBegin; 147 *geom = NULL; 148 PetscFunctionReturn(0); 149 } 150 151 static PetscErrorCode CreateFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms) 152 { 153 DMField coordField; 154 PetscInt Nf, f, maxDegree; 155 PetscErrorCode ierr; 156 157 PetscFunctionBeginUser; 158 *affineQuad = NULL; 159 *affineGeom = NULL; 160 *quads = NULL; 161 *geoms = NULL; 162 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 163 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 164 ierr = DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);CHKERRQ(ierr); 165 if (maxDegree <= 1) { 166 ierr = DMFieldCreateDefaultQuadrature(coordField, cellIS, affineQuad);CHKERRQ(ierr); 167 if (*affineQuad) {ierr = CellRangeGetFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom);CHKERRQ(ierr);} 168 } else { 169 ierr = PetscCalloc2(Nf, quads, Nf, geoms);CHKERRQ(ierr); 170 for (f = 0; f < Nf; ++f) { 171 PetscFE fe; 172 173 ierr = PetscDSGetDiscretization(ds, f, (PetscObject *) &fe);CHKERRQ(ierr); 174 ierr = PetscFEGetQuadrature(fe, &(*quads)[f]);CHKERRQ(ierr); 175 ierr = PetscObjectReference((PetscObject) (*quads)[f]);CHKERRQ(ierr); 176 ierr = CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]);CHKERRQ(ierr); 177 } 178 } 179 PetscFunctionReturn(0); 180 } 181 182 static PetscErrorCode DestroyFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms) 183 { 184 DMField coordField; 185 PetscInt Nf, f; 186 PetscErrorCode ierr; 187 188 PetscFunctionBeginUser; 189 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 190 ierr = DMGetCoordinateField(dm, &coordField);CHKERRQ(ierr); 191 if (*affineQuad) { 192 ierr = CellRangeRestoreFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom);CHKERRQ(ierr); 193 ierr = PetscQuadratureDestroy(affineQuad);CHKERRQ(ierr); 194 } else { 195 for (f = 0; f < Nf; ++f) { 196 ierr = CellRangeRestoreFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f]);CHKERRQ(ierr); 197 ierr = PetscQuadratureDestroy(&(*quads)[f]);CHKERRQ(ierr); 198 } 199 ierr = PetscFree2(*quads, *geoms);CHKERRQ(ierr); 200 } 201 PetscFunctionReturn(0); 202 } 203 204 static PetscErrorCode TestIntegration(DM dm, PetscInt cbs, PetscInt its) 205 { 206 PetscDS ds; 207 PetscFEGeom *chunkGeom = NULL; 208 PetscQuadrature affineQuad, *quads = NULL; 209 PetscFEGeom *affineGeom, **geoms = NULL; 210 PetscScalar *u, *elemVec; 211 IS cellIS; 212 PetscInt depth, cStart, cEnd, cell, chunkSize = cbs, Nch = 0, Nf, f, totDim, i, k; 213 #if defined(PETSC_USE_LOG) 214 PetscLogStage stage; 215 PetscLogEvent event; 216 #endif 217 PetscErrorCode ierr; 218 219 PetscFunctionBeginUser; 220 ierr = PetscLogStageRegister("PetscFE Residual Integration Test", &stage);CHKERRQ(ierr); 221 ierr = PetscLogEventRegister("FEIntegRes", PETSCFE_CLASSID, &event);CHKERRQ(ierr); 222 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 223 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 224 ierr = DMGetStratumIS(dm, "depth", depth, &cellIS);CHKERRQ(ierr); 225 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 226 ierr = DMGetCellDS(dm, cStart, &ds);CHKERRQ(ierr); 227 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr); 228 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr); 229 ierr = CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms);CHKERRQ(ierr); 230 ierr = PetscMalloc2(chunkSize*totDim, &u, chunkSize*totDim, &elemVec);CHKERRQ(ierr); 231 /* Assumptions: 232 - Single field 233 - No input data 234 - No auxiliary data 235 - No time-dependence 236 */ 237 for (i = 0; i < its; ++i) { 238 for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) { 239 const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS; 240 241 ierr = PetscArrayzero(elemVec, chunkSize*totDim);CHKERRQ(ierr); 242 /* TODO Replace with DMPlexGetCellFields() */ 243 for (k = 0; k < chunkSize*totDim; ++k) u[k] = 1.0; 244 for (f = 0; f < Nf; ++f) { 245 PetscFormKey key; 246 PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f]; 247 /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */ 248 249 key.label = NULL; key.value = 0; key.field = f; 250 ierr = PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom);CHKERRQ(ierr); 251 ierr = PetscLogEventBegin(event,0,0,0,0);CHKERRQ(ierr); 252 ierr = PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec);CHKERRQ(ierr); 253 ierr = PetscLogEventEnd(event,0,0,0,0);CHKERRQ(ierr); 254 } 255 } 256 } 257 ierr = PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom);CHKERRQ(ierr); 258 ierr = DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms);CHKERRQ(ierr); 259 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 260 ierr = PetscFree2(u, elemVec);CHKERRQ(ierr); 261 ierr = PetscLogStagePop();CHKERRQ(ierr); 262 #if defined(PETSC_USE_LOG) 263 { 264 const char *title = "Petsc FE Residual Integration"; 265 PetscEventPerfInfo eventInfo; 266 PetscInt N = (cEnd - cStart)*Nf*its; 267 PetscReal flopRate, cellRate; 268 269 ierr = PetscLogEventGetPerfInfo(stage, event, &eventInfo);CHKERRQ(ierr); 270 flopRate = eventInfo.time != 0.0 ? eventInfo.flops/eventInfo.time : 0.0; 271 cellRate = eventInfo.time != 0.0 ? N/eventInfo.time : 0.0; 272 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "%s: %D integrals %D chunks %D reps\n Cell rate: %.2f/s flop rate: %.2f MF/s\n", title, N, Nch, its, (double)cellRate, (double)(flopRate/1.e6));CHKERRQ(ierr); 273 } 274 #endif 275 PetscFunctionReturn(0); 276 } 277 278 static PetscErrorCode TestIntegration2(DM dm, PetscInt cbs, PetscInt its) 279 { 280 Vec X, F; 281 #if defined(PETSC_USE_LOG) 282 PetscLogStage stage; 283 #endif 284 PetscInt i; 285 PetscErrorCode ierr; 286 287 PetscFunctionBeginUser; 288 ierr = PetscLogStageRegister("DMPlex Residual Integration Test", &stage);CHKERRQ(ierr); 289 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 290 ierr = DMGetLocalVector(dm, &X);CHKERRQ(ierr); 291 ierr = DMGetLocalVector(dm, &F);CHKERRQ(ierr); 292 for (i = 0; i < its; ++i) { 293 ierr = DMPlexSNESComputeResidualFEM(dm, X, F, NULL);CHKERRQ(ierr); 294 } 295 ierr = DMRestoreLocalVector(dm, &X);CHKERRQ(ierr); 296 ierr = DMRestoreLocalVector(dm, &F);CHKERRQ(ierr); 297 ierr = PetscLogStagePop();CHKERRQ(ierr); 298 #if defined(PETSC_USE_LOG) 299 { 300 const char *title = "DMPlex Residual Integration"; 301 PetscEventPerfInfo eventInfo; 302 PetscReal flopRate, cellRate; 303 PetscInt cStart, cEnd, Nf, N; 304 PetscLogEvent event; 305 306 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 307 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 308 ierr = PetscLogEventGetId("DMPlexResidualFE", &event);CHKERRQ(ierr); 309 ierr = PetscLogEventGetPerfInfo(stage, event, &eventInfo);CHKERRQ(ierr); 310 N = (cEnd - cStart)*Nf*eventInfo.count; 311 flopRate = eventInfo.time != 0.0 ? eventInfo.flops/eventInfo.time : 0.0; 312 cellRate = eventInfo.time != 0.0 ? N/eventInfo.time : 0.0; 313 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "%s: %D integrals %D reps\n Cell rate: %.2f/s flop rate: %.2f MF/s\n", title, N, eventInfo.count, (double)cellRate, (double)(flopRate/1.e6));CHKERRQ(ierr); 314 } 315 #endif 316 PetscFunctionReturn(0); 317 } 318 319 int main(int argc, char **argv) 320 { 321 DM dm; 322 AppCtx ctx; 323 PetscMPIInt size; 324 PetscErrorCode ierr; 325 326 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 327 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); 328 if (size > 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a uniprocessor example only."); 329 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 330 ierr = PetscLogDefaultBegin();CHKERRQ(ierr); 331 ierr = DMCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr); 332 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); 333 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 334 ierr = PetscObjectSetName((PetscObject) dm, "Mesh");CHKERRQ(ierr); 335 ierr = PetscObjectViewFromOptions((PetscObject) dm, NULL, "-dm_view");CHKERRQ(ierr); 336 ierr = SetupDiscretization(dm, "potential", SetupPrimalProblem, &ctx);CHKERRQ(ierr); 337 ierr = TestIntegration(dm, ctx.cbs, ctx.its);CHKERRQ(ierr); 338 ierr = TestIntegration2(dm, ctx.cbs, ctx.its);CHKERRQ(ierr); 339 ierr = DMDestroy(&dm);CHKERRQ(ierr); 340 ierr = PetscFinalize(); 341 return ierr; 342 } 343 344 /*TEST 345 test: 346 suffix: 0 347 requires: triangle 348 args: -dm_view 349 350 test: 351 suffix: 1 352 requires: triangle 353 args: -dm_view -potential_petscspace_degree 1 354 355 test: 356 suffix: 2 357 requires: triangle 358 args: -dm_view -potential_petscspace_degree 2 359 360 test: 361 suffix: 3 362 requires: triangle 363 args: -dm_view -potential_petscspace_degree 3 364 TEST*/ 365