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 PetscFunctionBeginUser; 17 options->dim = 2; 18 options->simplex = PETSC_TRUE; 19 options->its = 1; 20 options->cbs = 8; 21 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 } 30 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 } 38 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 } 44 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 } 50 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 } 56 57 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 58 { 59 PetscDS prob; 60 DMLabel label; 61 const PetscInt id = 1; 62 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, (PetscVoidFn *)trig_u, NULL, user, NULL)); 70 PetscFunctionReturn(PETSC_SUCCESS); 71 } 72 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]; 78 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 } 96 97 /* PetscObjectContainerCompose() compose requires void ** signature on destructor */ 98 static PetscErrorCode PetscFEGeomDestroy_Void(void **ctx) 99 { 100 return PetscFEGeomDestroy((PetscFEGeom **)ctx); 101 } 102 103 PetscErrorCode CellRangeGetFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscFEGeomMode mode, PetscFEGeom **geom) 104 { 105 char composeStr[33] = {0}; 106 PetscObjectId id; 107 PetscContainer container; 108 109 PetscFunctionBegin; 110 PetscCall(PetscObjectGetId((PetscObject)quad, &id)); 111 PetscCall(PetscSNPrintf(composeStr, 32, "CellRangeGetFEGeom_%" PetscInt64_FMT "\n", id)); 112 PetscCall(PetscObjectQuery((PetscObject)cellIS, composeStr, (PetscObject *)&container)); 113 if (container) { 114 PetscCall(PetscContainerGetPointer(container, (void **)geom)); 115 } else { 116 PetscCall(DMFieldCreateFEGeom(coordField, cellIS, quad, mode, geom)); 117 PetscCall(PetscObjectContainerCompose((PetscObject)cellIS, composeStr, *geom, PetscFEGeomDestroy_Void)); 118 } 119 PetscFunctionReturn(PETSC_SUCCESS); 120 } 121 122 PetscErrorCode CellRangeRestoreFEGeom(IS cellIS, DMField coordField, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom) 123 { 124 PetscFunctionBegin; 125 *geom = NULL; 126 PetscFunctionReturn(PETSC_SUCCESS); 127 } 128 129 static PetscErrorCode CreateFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms) 130 { 131 DMField coordField; 132 PetscInt Nf, f, maxDegree; 133 134 PetscFunctionBeginUser; 135 *affineQuad = NULL; 136 *affineGeom = NULL; 137 *quads = NULL; 138 *geoms = NULL; 139 PetscCall(PetscDSGetNumFields(ds, &Nf)); 140 PetscCall(DMGetCoordinateField(dm, &coordField)); 141 PetscCall(DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree)); 142 if (maxDegree <= 1) { 143 PetscCall(DMFieldCreateDefaultQuadrature(coordField, cellIS, affineQuad)); 144 if (*affineQuad) PetscCall(CellRangeGetFEGeom(cellIS, coordField, *affineQuad, PETSC_FEGEOM_BASIC, affineGeom)); 145 } else { 146 PetscCall(PetscCalloc2(Nf, quads, Nf, geoms)); 147 for (f = 0; f < Nf; ++f) { 148 PetscFE fe; 149 150 PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe)); 151 PetscCall(PetscFEGetQuadrature(fe, &(*quads)[f])); 152 PetscCall(PetscObjectReference((PetscObject)(*quads)[f])); 153 PetscCall(CellRangeGetFEGeom(cellIS, coordField, (*quads)[f], PETSC_FEGEOM_BASIC, &(*geoms)[f])); 154 } 155 } 156 PetscFunctionReturn(PETSC_SUCCESS); 157 } 158 159 static PetscErrorCode DestroyFEGeometry(DM dm, PetscDS ds, IS cellIS, PetscQuadrature *affineQuad, PetscFEGeom **affineGeom, PetscQuadrature **quads, PetscFEGeom ***geoms) 160 { 161 DMField coordField; 162 PetscInt Nf, f; 163 164 PetscFunctionBeginUser; 165 PetscCall(PetscDSGetNumFields(ds, &Nf)); 166 PetscCall(DMGetCoordinateField(dm, &coordField)); 167 if (*affineQuad) { 168 PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, *affineQuad, PETSC_FALSE, affineGeom)); 169 PetscCall(PetscQuadratureDestroy(affineQuad)); 170 } else { 171 for (f = 0; f < Nf; ++f) { 172 PetscCall(CellRangeRestoreFEGeom(cellIS, coordField, (*quads)[f], PETSC_FALSE, &(*geoms)[f])); 173 PetscCall(PetscQuadratureDestroy(&(*quads)[f])); 174 } 175 PetscCall(PetscFree2(*quads, *geoms)); 176 } 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 180 static PetscErrorCode TestIntegration(DM dm, PetscInt cbs, PetscInt its) 181 { 182 PetscDS ds; 183 PetscFEGeom *chunkGeom = NULL; 184 PetscQuadrature affineQuad, *quads = NULL; 185 PetscFEGeom *affineGeom, **geoms = NULL; 186 PetscScalar *u, *elemVec; 187 IS cellIS; 188 PetscInt depth, cStart, cEnd, cell, chunkSize = cbs, Nch = 0, Nf, f, totDim, i, k; 189 PetscLogStage stage; 190 PetscLogEvent event; 191 192 PetscFunctionBeginUser; 193 PetscCall(PetscLogStageRegister("PetscFE Residual Integration Test", &stage)); 194 PetscCall(PetscLogEventRegister("FEIntegRes", PETSCFE_CLASSID, &event)); 195 PetscCall(PetscLogStagePush(stage)); 196 PetscCall(DMPlexGetDepth(dm, &depth)); 197 PetscCall(DMGetStratumIS(dm, "depth", depth, &cellIS)); 198 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 199 PetscCall(DMGetCellDS(dm, cStart, &ds, NULL)); 200 PetscCall(PetscDSGetNumFields(ds, &Nf)); 201 PetscCall(PetscDSGetTotalDimension(ds, &totDim)); 202 PetscCall(CreateFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms)); 203 PetscCall(PetscMalloc2(chunkSize * totDim, &u, chunkSize * totDim, &elemVec)); 204 /* Assumptions: 205 - Single field 206 - No input data 207 - No auxiliary data 208 - No time-dependence 209 */ 210 for (i = 0; i < its; ++i) { 211 for (cell = cStart; cell < cEnd; cell += chunkSize, ++Nch) { 212 const PetscInt cS = cell, cE = PetscMin(cS + chunkSize, cEnd), Ne = cE - cS; 213 214 PetscCall(PetscArrayzero(elemVec, chunkSize * totDim)); 215 /* TODO Replace with DMPlexGetCellFields() */ 216 for (k = 0; k < chunkSize * totDim; ++k) u[k] = 1.0; 217 for (f = 0; f < Nf; ++f) { 218 PetscFormKey key; 219 PetscFEGeom *geom = affineGeom ? affineGeom : geoms[f]; 220 /* PetscQuadrature quad = affineQuad ? affineQuad : quads[f]; */ 221 222 key.label = NULL; 223 key.value = 0; 224 key.field = f; 225 key.part = 0; 226 PetscCall(PetscFEGeomGetChunk(geom, cS, cE, &chunkGeom)); 227 PetscCall(PetscLogEventBegin(event, 0, 0, 0, 0)); 228 PetscCall(PetscFEIntegrateResidual(ds, key, Ne, chunkGeom, u, NULL, NULL, NULL, 0.0, elemVec)); 229 PetscCall(PetscLogEventEnd(event, 0, 0, 0, 0)); 230 } 231 } 232 } 233 PetscCall(PetscFEGeomRestoreChunk(affineGeom, cStart, cEnd, &chunkGeom)); 234 PetscCall(DestroyFEGeometry(dm, ds, cellIS, &affineQuad, &affineGeom, &quads, &geoms)); 235 PetscCall(ISDestroy(&cellIS)); 236 PetscCall(PetscFree2(u, elemVec)); 237 PetscCall(PetscLogStagePop()); 238 if (PetscDefined(USE_LOG)) { 239 const char *title = "PETSc FE Residual Integration"; 240 PetscEventPerfInfo eventInfo; 241 PetscInt N = (cEnd - cStart) * Nf * its; 242 PetscReal flopRate, cellRate; 243 244 PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo)); 245 flopRate = eventInfo.time != 0.0 ? eventInfo.flops / eventInfo.time : 0.0; 246 cellRate = eventInfo.time != 0.0 ? N / eventInfo.time : 0.0; 247 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))); 248 } 249 PetscFunctionReturn(PETSC_SUCCESS); 250 } 251 252 static PetscErrorCode TestIntegration2(DM dm, PetscInt cbs, PetscInt its) 253 { 254 Vec X, F; 255 PetscLogStage stage; 256 PetscInt i; 257 258 PetscFunctionBeginUser; 259 PetscCall(PetscLogStageRegister("DMPlex Residual Integration Test", &stage)); 260 PetscCall(PetscLogStagePush(stage)); 261 PetscCall(DMGetLocalVector(dm, &X)); 262 PetscCall(DMGetLocalVector(dm, &F)); 263 for (i = 0; i < its; ++i) PetscCall(DMPlexSNESComputeResidualFEM(dm, X, F, NULL)); 264 PetscCall(DMRestoreLocalVector(dm, &X)); 265 PetscCall(DMRestoreLocalVector(dm, &F)); 266 PetscCall(PetscLogStagePop()); 267 if (PetscDefined(USE_LOG)) { 268 const char *title = "DMPlex Residual Integration"; 269 PetscEventPerfInfo eventInfo; 270 PetscReal flopRate, cellRate; 271 PetscInt cStart, cEnd, Nf, N; 272 PetscLogEvent event; 273 274 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 275 PetscCall(DMGetNumFields(dm, &Nf)); 276 PetscCall(PetscLogEventGetId("DMPlexResidualFE", &event)); 277 PetscCall(PetscLogEventGetPerfInfo(stage, event, &eventInfo)); 278 N = (cEnd - cStart) * Nf * eventInfo.count; 279 flopRate = eventInfo.time != 0.0 ? eventInfo.flops / eventInfo.time : 0.0; 280 cellRate = eventInfo.time != 0.0 ? N / eventInfo.time : 0.0; 281 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))); 282 } 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285 286 int main(int argc, char **argv) 287 { 288 DM dm; 289 AppCtx ctx; 290 PetscMPIInt size; 291 292 PetscFunctionBeginUser; 293 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 294 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 295 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only."); 296 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); 297 PetscCall(PetscLogDefaultBegin()); 298 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 299 PetscCall(DMSetType(dm, DMPLEX)); 300 PetscCall(DMSetFromOptions(dm)); 301 PetscCall(PetscObjectSetName((PetscObject)dm, "Mesh")); 302 PetscCall(PetscObjectViewFromOptions((PetscObject)dm, NULL, "-dm_view")); 303 PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &ctx)); 304 PetscCall(TestIntegration(dm, ctx.cbs, ctx.its)); 305 PetscCall(TestIntegration2(dm, ctx.cbs, ctx.its)); 306 PetscCall(DMDestroy(&dm)); 307 PetscCall(PetscFinalize()); 308 return 0; 309 } 310 311 /*TEST 312 test: 313 suffix: 0 314 requires: triangle 315 args: -dm_view 316 317 test: 318 suffix: 1 319 requires: triangle 320 args: -dm_view -potential_petscspace_degree 1 321 322 test: 323 suffix: 2 324 requires: triangle 325 args: -dm_view -potential_petscspace_degree 2 326 327 test: 328 suffix: 3 329 requires: triangle 330 args: -dm_view -potential_petscspace_degree 3 331 TEST*/ 332