1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 // libCEED + PETSc Example: CEED BPs 18 // 19 // This example demonstrates a simple usage of libCEED with PETSc to solve the 20 // CEED BP benchmark problems, see http://ceed.exascaleproject.org/bps. 21 // 22 // The code uses higher level communication protocols in DMPlex. 23 // 24 // Build with: 25 // 26 // make bps [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 27 // 28 // Sample runs: 29 // 30 // ./bps -problem bp1 -degree 3 31 // ./bps -problem bp2 -degree 3 32 // ./bps -problem bp3 -degree 3 33 // ./bps -problem bp4 -degree 3 34 // ./bps -problem bp5 -degree 3 -ceed /cpu/self 35 // ./bps -problem bp6 -degree 3 -ceed /gpu/cuda 36 // 37 //TESTARGS -ceed {ceed_resource} -test -problem bp5 -degree 3 -ksp_max_it_clip 15,15 38 39 /// @file 40 /// CEED BPs example using PETSc with DMPlex 41 /// See bpsraw.c for a "raw" implementation using a structured grid. 42 const char help[] = "Solve CEED BPs using PETSc with DMPlex\n"; 43 44 #include <ceed.h> 45 #include <petscdmplex.h> 46 #include <petscksp.h> 47 #include <stdbool.h> 48 #include <string.h> 49 #include "setup.h" 50 51 // ----------------------------------------------------------------------------- 52 // Utilities 53 // ----------------------------------------------------------------------------- 54 55 // Utility function, compute three factors of an integer 56 static void Split3(PetscInt size, PetscInt m[3], bool reverse) { 57 for (PetscInt d=0,sizeleft=size; d<3; d++) { 58 PetscInt try = (PetscInt)PetscCeilReal(PetscPowReal(sizeleft, 1./(3 - d))); 59 while (try * (sizeleft / try) != sizeleft) try++; 60 m[reverse ? 2-d : d] = try; 61 sizeleft /= try; 62 } 63 } 64 65 static int Max3(const PetscInt a[3]) { 66 return PetscMax(a[0], PetscMax(a[1], a[2])); 67 } 68 69 static int Min3(const PetscInt a[3]) { 70 return PetscMin(a[0], PetscMin(a[1], a[2])); 71 } 72 73 // ----------------------------------------------------------------------------- 74 // Parameter structure for running problems 75 // ----------------------------------------------------------------------------- 76 typedef struct RunParams_ *RunParams; 77 struct RunParams_ { 78 MPI_Comm comm; 79 PetscBool test_mode, read_mesh, userlnodes, 80 petschavecuda, write_solution; 81 char *filename, *hostname; 82 PetscInt localnodes, degree, qextra, dim, ncompu, *melem; 83 PetscInt ksp_max_it_clip[2]; 84 PetscMPIInt rankspernode; 85 bpType bpchoice; 86 PetscLogStage solvestage; 87 }; 88 89 // ----------------------------------------------------------------------------- 90 // Main body of program, called in a loop for performance benchmarking purposes 91 // ----------------------------------------------------------------------------- 92 static PetscErrorCode RunWithDM(RunParams rp, DM dm, 93 const char *ceedresource) { 94 PetscErrorCode ierr; 95 double my_rt_start, my_rt, rt_min, rt_max; 96 PetscInt xlsize, lsize, gsize; 97 PetscScalar *r; 98 Vec X, Xloc, rhs, rhsloc; 99 Mat matO; 100 KSP ksp; 101 UserO userO; 102 Ceed ceed; 103 CeedData ceeddata; 104 CeedQFunction qferror; 105 CeedOperator operror; 106 CeedVector rhsceed, target; 107 VecType vectype; 108 PetscMemType memtype; 109 110 PetscFunctionBeginUser; 111 // Set up libCEED 112 CeedInit(ceedresource, &ceed); 113 CeedMemType memtypebackend; 114 CeedGetPreferredMemType(ceed, &memtypebackend); 115 116 ierr = DMGetVecType(dm, &vectype); CHKERRQ(ierr); 117 if (!vectype) { // Not yet set by user -dm_vec_type 118 switch (memtypebackend) { 119 case CEED_MEM_HOST: vectype = VECSTANDARD; break; 120 case CEED_MEM_DEVICE: { 121 const char *resolved; 122 CeedGetResource(ceed, &resolved); 123 if (strstr(resolved, "/gpu/cuda")) vectype = VECCUDA; 124 else if (strstr(resolved, "/gpu/hip/occa")) 125 vectype = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 126 else if (strstr(resolved, "/gpu/hip")) vectype = VECHIP; 127 else vectype = VECSTANDARD; 128 } 129 } 130 ierr = DMSetVecType(dm, vectype); CHKERRQ(ierr); 131 } 132 133 // Create global and local solution vectors 134 ierr = DMCreateGlobalVector(dm, &X); CHKERRQ(ierr); 135 ierr = VecGetLocalSize(X, &lsize); CHKERRQ(ierr); 136 ierr = VecGetSize(X, &gsize); CHKERRQ(ierr); 137 ierr = DMCreateLocalVector(dm, &Xloc); CHKERRQ(ierr); 138 ierr = VecGetSize(Xloc, &xlsize); CHKERRQ(ierr); 139 ierr = VecDuplicate(X, &rhs); CHKERRQ(ierr); 140 141 // Operator 142 ierr = PetscMalloc1(1, &userO); CHKERRQ(ierr); 143 ierr = MatCreateShell(rp->comm, lsize, lsize, gsize, gsize, 144 userO, &matO); CHKERRQ(ierr); 145 ierr = MatShellSetOperation(matO, MATOP_MULT, 146 (void(*)(void))MatMult_Ceed); CHKERRQ(ierr); 147 ierr = MatShellSetOperation(matO, MATOP_GET_DIAGONAL, 148 (void(*)(void))MatGetDiag); CHKERRQ(ierr); 149 ierr = MatShellSetVecType(matO, vectype); CHKERRQ(ierr); 150 151 // Print summary 152 if (!rp->test_mode) { 153 PetscInt P = rp->degree + 1, Q = P + rp->qextra; 154 155 const char *usedresource; 156 CeedGetResource(ceed, &usedresource); 157 158 VecType vectype; 159 ierr = VecGetType(X, &vectype); CHKERRQ(ierr); 160 161 PetscInt cStart, cEnd; 162 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr); 163 PetscMPIInt comm_size; 164 ierr = MPI_Comm_size(rp->comm, &comm_size); CHKERRQ(ierr); 165 ierr = PetscPrintf(rp->comm, 166 "\n-- CEED Benchmark Problem %d -- libCEED + PETSc --\n" 167 " MPI:\n" 168 " Hostname : %s\n" 169 " Total ranks : %d\n" 170 " Ranks per compute node : %d\n" 171 " PETSc:\n" 172 " PETSc Vec Type : %s\n" 173 " libCEED:\n" 174 " libCEED Backend : %s\n" 175 " libCEED Backend MemType : %s\n" 176 " Mesh:\n" 177 " Number of 1D Basis Nodes (P) : %d\n" 178 " Number of 1D Quadrature Points (Q) : %d\n" 179 " Global nodes : %D\n" 180 " Local Elements : %D\n" 181 " Owned nodes : %D\n" 182 " DoF per node : %D\n", 183 rp->bpchoice+1, 184 rp->hostname, 185 comm_size, rp->rankspernode, 186 vectype, usedresource, 187 CeedMemTypes[memtypebackend], 188 P, Q, gsize/rp->ncompu, cEnd - cStart, lsize/rp->ncompu, 189 rp->ncompu); 190 CHKERRQ(ierr); 191 } 192 193 // Create RHS vector 194 ierr = VecDuplicate(Xloc, &rhsloc); CHKERRQ(ierr); 195 ierr = VecZeroEntries(rhsloc); CHKERRQ(ierr); 196 ierr = VecGetArrayAndMemType(rhsloc, &r, &memtype); CHKERRQ(ierr); 197 CeedVectorCreate(ceed, xlsize, &rhsceed); 198 CeedVectorSetArray(rhsceed, MemTypeP2C(memtype), CEED_USE_POINTER, r); 199 200 ierr = PetscMalloc1(1, &ceeddata); CHKERRQ(ierr); 201 ierr = SetupLibceedByDegree(dm, ceed, rp->degree, rp->dim, rp->qextra, 202 rp->ncompu, gsize, xlsize, rp->bpchoice, ceeddata, 203 true, rhsceed, &target); CHKERRQ(ierr); 204 205 // Gather RHS 206 CeedVectorTakeArray(rhsceed, MemTypeP2C(memtype), NULL); 207 ierr = VecRestoreArrayAndMemType(rhsloc, &r); CHKERRQ(ierr); 208 ierr = VecZeroEntries(rhs); CHKERRQ(ierr); 209 ierr = DMLocalToGlobal(dm, rhsloc, ADD_VALUES, rhs); CHKERRQ(ierr); 210 CeedVectorDestroy(&rhsceed); 211 212 // Create the error QFunction 213 CeedQFunctionCreateInterior(ceed, 1, bpOptions[rp->bpchoice].error, 214 bpOptions[rp->bpchoice].errorfname, &qferror); 215 CeedQFunctionAddInput(qferror, "u", rp->ncompu, CEED_EVAL_INTERP); 216 CeedQFunctionAddInput(qferror, "true_soln", rp->ncompu, CEED_EVAL_NONE); 217 CeedQFunctionAddOutput(qferror, "error", rp->ncompu, CEED_EVAL_NONE); 218 219 // Create the error operator 220 CeedOperatorCreate(ceed, qferror, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 221 &operror); 222 CeedOperatorSetField(operror, "u", ceeddata->Erestrictu, 223 ceeddata->basisu, CEED_VECTOR_ACTIVE); 224 CeedOperatorSetField(operror, "true_soln", ceeddata->Erestrictui, 225 CEED_BASIS_COLLOCATED, target); 226 CeedOperatorSetField(operror, "error", ceeddata->Erestrictui, 227 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 228 229 // Set up Mat 230 userO->comm = rp->comm; 231 userO->dm = dm; 232 userO->Xloc = Xloc; 233 ierr = VecDuplicate(Xloc, &userO->Yloc); CHKERRQ(ierr); 234 userO->xceed = ceeddata->xceed; 235 userO->yceed = ceeddata->yceed; 236 userO->op = ceeddata->opapply; 237 userO->ceed = ceed; 238 239 ierr = KSPCreate(rp->comm, &ksp); CHKERRQ(ierr); 240 { 241 PC pc; 242 ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); 243 if (rp->bpchoice == CEED_BP1 || rp->bpchoice == CEED_BP2) { 244 ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr); 245 ierr = PCJacobiSetType(pc, PC_JACOBI_ROWSUM); CHKERRQ(ierr); 246 } else { 247 ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr); 248 } 249 ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr); 250 ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr); 251 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 252 PETSC_DEFAULT); CHKERRQ(ierr); 253 } 254 ierr = KSPSetOperators(ksp, matO, matO); CHKERRQ(ierr); 255 256 // First run's performance log is not considered for benchmarking purposes 257 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1); 258 CHKERRQ(ierr); 259 my_rt_start = MPI_Wtime(); 260 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 261 my_rt = MPI_Wtime() - my_rt_start; 262 ierr = MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, rp->comm); 263 CHKERRQ(ierr); 264 // Set maxits based on first iteration timing 265 if (my_rt > 0.02) { 266 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 267 rp->ksp_max_it_clip[0]); 268 CHKERRQ(ierr); 269 } else { 270 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 271 rp->ksp_max_it_clip[1]); 272 CHKERRQ(ierr); 273 } 274 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); 275 276 // Timed solve 277 ierr = VecZeroEntries(X); CHKERRQ(ierr); 278 ierr = PetscBarrier((PetscObject)ksp); CHKERRQ(ierr); 279 280 // -- Performance logging 281 ierr = PetscLogStagePush(rp->solvestage); CHKERRQ(ierr); 282 283 // -- Solve 284 my_rt_start = MPI_Wtime(); 285 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 286 my_rt = MPI_Wtime() - my_rt_start; 287 288 // -- Performance logging 289 ierr = PetscLogStagePop(); 290 291 // Output results 292 { 293 KSPType ksptype; 294 KSPConvergedReason reason; 295 PetscReal rnorm; 296 PetscInt its; 297 ierr = KSPGetType(ksp, &ksptype); CHKERRQ(ierr); 298 ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr); 299 ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr); 300 ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr); 301 if (!rp->test_mode || reason < 0 || rnorm > 1e-8) { 302 ierr = PetscPrintf(rp->comm, 303 " KSP:\n" 304 " KSP Type : %s\n" 305 " KSP Convergence : %s\n" 306 " Total KSP Iterations : %D\n" 307 " Final rnorm : %e\n", 308 ksptype, KSPConvergedReasons[reason], its, 309 (double)rnorm); CHKERRQ(ierr); 310 } 311 if (!rp->test_mode) { 312 ierr = PetscPrintf(rp->comm," Performance:\n"); CHKERRQ(ierr); 313 } 314 { 315 PetscReal maxerror; 316 ierr = ComputeErrorMax(userO, operror, X, target, &maxerror); 317 CHKERRQ(ierr); 318 PetscReal tol = 5e-2; 319 if (!rp->test_mode || maxerror > tol) { 320 ierr = MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, rp->comm); 321 CHKERRQ(ierr); 322 ierr = MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, rp->comm); 323 CHKERRQ(ierr); 324 ierr = PetscPrintf(rp->comm, 325 " Pointwise Error (max) : %e\n" 326 " CG Solve Time : %g (%g) sec\n", 327 (double)maxerror, rt_max, rt_min); CHKERRQ(ierr); 328 } 329 } 330 if (!rp->test_mode) { 331 ierr = PetscPrintf(rp->comm, 332 " DoFs/Sec in CG : %g (%g) million\n", 333 1e-6*gsize*its/rt_max, 334 1e-6*gsize*its/rt_min); CHKERRQ(ierr); 335 } 336 } 337 338 if (rp->write_solution) { 339 PetscViewer vtkviewersoln; 340 341 ierr = PetscViewerCreate(rp->comm, &vtkviewersoln); CHKERRQ(ierr); 342 ierr = PetscViewerSetType(vtkviewersoln, PETSCVIEWERVTK); CHKERRQ(ierr); 343 ierr = PetscViewerFileSetName(vtkviewersoln, "solution.vtk"); CHKERRQ(ierr); 344 ierr = VecView(X, vtkviewersoln); CHKERRQ(ierr); 345 ierr = PetscViewerDestroy(&vtkviewersoln); CHKERRQ(ierr); 346 } 347 348 // Cleanup 349 ierr = VecDestroy(&X); CHKERRQ(ierr); 350 ierr = VecDestroy(&Xloc); CHKERRQ(ierr); 351 ierr = VecDestroy(&userO->Yloc); CHKERRQ(ierr); 352 ierr = MatDestroy(&matO); CHKERRQ(ierr); 353 ierr = PetscFree(userO); CHKERRQ(ierr); 354 ierr = CeedDataDestroy(0, ceeddata); CHKERRQ(ierr); 355 356 ierr = VecDestroy(&rhs); CHKERRQ(ierr); 357 ierr = VecDestroy(&rhsloc); CHKERRQ(ierr); 358 ierr = KSPDestroy(&ksp); CHKERRQ(ierr); 359 CeedVectorDestroy(&target); 360 CeedQFunctionDestroy(&qferror); 361 CeedOperatorDestroy(&operror); 362 CeedDestroy(&ceed); 363 PetscFunctionReturn(0); 364 } 365 366 static PetscErrorCode Run(RunParams rp, 367 PetscInt num_resources, char *const *ceedresources, 368 PetscInt num_bpchoices, const bpType *bpchoices) { 369 PetscInt ierr; 370 DM dm; 371 372 PetscFunctionBeginUser; 373 // Setup DM 374 if (rp->read_mesh) { 375 ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, rp->filename, PETSC_TRUE, &dm); 376 CHKERRQ(ierr); 377 } else { 378 if (rp->userlnodes) { 379 // Find a nicely composite number of elements no less than global nodes 380 PetscMPIInt size; 381 ierr = MPI_Comm_size(rp->comm, &size); CHKERRQ(ierr); 382 for (PetscInt gelem = 383 PetscMax(1, size * rp->localnodes / PetscPowInt(rp->degree, rp->dim)); 384 ; 385 gelem++) { 386 Split3(gelem, rp->melem, true); 387 if (Max3(rp->melem) / Min3(rp->melem) <= 2) break; 388 } 389 } 390 ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, rp->dim, PETSC_FALSE, rp->melem, 391 NULL, NULL, NULL, PETSC_TRUE, &dm); CHKERRQ(ierr); 392 } 393 394 { 395 DM dmDist = NULL; 396 PetscPartitioner part; 397 398 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 399 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 400 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 401 if (dmDist) { 402 ierr = DMDestroy(&dm); CHKERRQ(ierr); 403 dm = dmDist; 404 } 405 } 406 // Disable default VECSTANDARD *after* distribution (which creates a Vec) 407 ierr = DMSetVecType(dm, NULL); CHKERRQ(ierr); 408 409 for (PetscInt b = 0; b < num_bpchoices; b++) { 410 DM dm_deg; 411 VecType vectype; 412 PetscInt qextra = rp->qextra; 413 rp->bpchoice = bpchoices[b]; 414 rp->ncompu = bpOptions[rp->bpchoice].ncompu; 415 rp->qextra = qextra < 0 ? bpOptions[rp->bpchoice].qextra : qextra; 416 ierr = DMClone(dm, &dm_deg); CHKERRQ(ierr); 417 ierr = DMGetVecType(dm, &vectype); CHKERRQ(ierr); 418 ierr = DMSetVecType(dm_deg, vectype); CHKERRQ(ierr); 419 // Create DM 420 ierr = SetupDMByDegree(dm_deg, rp->degree, rp->ncompu, rp->bpchoice); 421 CHKERRQ(ierr); 422 for (PetscInt r = 0; r < num_resources; r++) { 423 ierr = RunWithDM(rp, dm_deg, ceedresources[r]); CHKERRQ(ierr); 424 } 425 ierr = DMDestroy(&dm_deg); CHKERRQ(ierr); 426 rp->qextra = qextra; 427 } 428 429 ierr = DMDestroy(&dm); CHKERRQ(ierr); 430 PetscFunctionReturn(0); 431 } 432 433 int main(int argc, char **argv) { 434 PetscInt ierr, commsize; 435 RunParams rp; 436 MPI_Comm comm; 437 char filename[PETSC_MAX_PATH_LEN]; 438 char *ceedresources[30]; 439 PetscInt num_ceedresources = 30; 440 char hostname[PETSC_MAX_PATH_LEN]; 441 442 PetscInt dim = 3, melem[3] = {3, 3, 3}; 443 PetscInt num_degrees = 30, degree[30] = {}, num_localnodes = 2, localnodes[2] = {}; 444 PetscMPIInt rankspernode; 445 PetscBool degree_set; 446 bpType bpchoices[10]; 447 PetscInt num_bpchoices = 10; 448 449 // Check PETSc CUDA support 450 PetscBool petschavecuda; 451 // *INDENT-OFF* 452 #ifdef PETSC_HAVE_CUDA 453 petschavecuda = PETSC_TRUE; 454 #else 455 petschavecuda = PETSC_FALSE; 456 #endif 457 // *INDENT-ON* 458 459 // Initialize PETSc 460 ierr = PetscInitialize(&argc, &argv, NULL, help); 461 if (ierr) return ierr; 462 comm = PETSC_COMM_WORLD; 463 ierr = MPI_Comm_size(comm, &commsize); 464 if (ierr != MPI_SUCCESS) return ierr; 465 #if defined(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY) 466 { 467 MPI_Comm splitcomm; 468 ierr = MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, 469 &splitcomm); 470 CHKERRQ(ierr); 471 ierr = MPI_Comm_size(splitcomm, &rankspernode); CHKERRQ(ierr); 472 ierr = MPI_Comm_free(&splitcomm); CHKERRQ(ierr); 473 } 474 #else 475 rankspernode = -1; // Unknown 476 #endif 477 478 // Setup all parameters needed in Run() 479 ierr = PetscMalloc1(1, &rp); CHKERRQ(ierr); 480 rp->comm = comm; 481 482 // Read command line options 483 ierr = PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); 484 CHKERRQ(ierr); 485 { 486 PetscBool set; 487 ierr = PetscOptionsEnumArray("-problem", "CEED benchmark problem to solve", 488 NULL, 489 bpTypes, (PetscEnum *)bpchoices, &num_bpchoices, &set); 490 CHKERRQ(ierr); 491 if (!set) { 492 bpchoices[0] = CEED_BP1; 493 num_bpchoices = 1; 494 } 495 } 496 rp->test_mode = PETSC_FALSE; 497 ierr = PetscOptionsBool("-test", 498 "Testing mode (do not print unless error is large)", 499 NULL, rp->test_mode, &rp->test_mode, NULL); CHKERRQ(ierr); 500 rp->write_solution = PETSC_FALSE; 501 ierr = PetscOptionsBool("-write_solution", "Write solution for visualization", 502 NULL, rp->write_solution, &rp->write_solution, NULL); 503 CHKERRQ(ierr); 504 degree[0] = rp->test_mode ? 3 : 2; 505 ierr = PetscOptionsIntArray("-degree", 506 "Polynomial degree of tensor product basis", NULL, 507 degree, &num_degrees, °ree_set); CHKERRQ(ierr); 508 if (!degree_set) 509 num_degrees = 1; 510 rp->qextra = PETSC_DECIDE; 511 ierr = PetscOptionsInt("-qextra", 512 "Number of extra quadrature points (-1 for auto)", NULL, 513 rp->qextra, &rp->qextra, NULL); CHKERRQ(ierr); 514 { 515 PetscBool set; 516 ierr = PetscOptionsStringArray("-ceed", 517 "CEED resource specifier (comma-separated list)", NULL, 518 ceedresources, &num_ceedresources, &set); CHKERRQ(ierr); 519 if (!set) { 520 ierr = PetscStrallocpy( "/cpu/self", &ceedresources[0]); CHKERRQ(ierr); 521 num_ceedresources = 1; 522 } 523 } 524 ierr = PetscGetHostName(hostname, sizeof hostname); CHKERRQ(ierr); 525 ierr = PetscOptionsString("-hostname", "Hostname for output", NULL, hostname, 526 hostname, sizeof(hostname), NULL); CHKERRQ(ierr); 527 rp->read_mesh = PETSC_FALSE; 528 ierr = PetscOptionsString("-mesh", "Read mesh from file", NULL, filename, 529 filename, sizeof(filename), &rp->read_mesh); 530 CHKERRQ(ierr); 531 rp->filename = filename; 532 if (!rp->read_mesh) { 533 PetscInt tmp = dim; 534 ierr = PetscOptionsIntArray("-cells", "Number of cells per dimension", NULL, 535 melem, &tmp, NULL); CHKERRQ(ierr); 536 } 537 localnodes[0] = 1000; 538 ierr = PetscOptionsIntArray("-local_nodes", 539 "Target number of locally owned nodes per " 540 "process (single value or min,max)", 541 NULL, localnodes, &num_localnodes, &rp->userlnodes); 542 CHKERRQ(ierr); 543 if (num_localnodes < 2) 544 localnodes[1] = 2 * localnodes[0]; 545 { 546 PetscInt two = 2; 547 rp->ksp_max_it_clip[0] = 5; 548 rp->ksp_max_it_clip[1] = 20; 549 ierr = PetscOptionsIntArray("-ksp_max_it_clip", 550 "Min and max number of iterations to use during benchmarking", 551 NULL, rp->ksp_max_it_clip, &two, NULL); CHKERRQ(ierr); 552 } 553 if (!degree_set) { 554 PetscInt maxdegree = 8; 555 ierr = PetscOptionsInt("-max_degree", 556 "Range of degrees [1, maxdegree] to run with", 557 NULL, maxdegree, &maxdegree, NULL); 558 CHKERRQ(ierr); 559 for (PetscInt i = 0; i < maxdegree; i++) 560 degree[i] = i + 1; 561 num_degrees = maxdegree; 562 } 563 { 564 PetscBool flg; 565 PetscInt p = rankspernode; 566 ierr = PetscOptionsInt("-p", "Number of MPI ranks per node", NULL, 567 p, &p, &flg); 568 CHKERRQ(ierr); 569 if (flg) rankspernode = p; 570 } 571 572 ierr = PetscOptionsEnd(); 573 CHKERRQ(ierr); 574 575 // Register PETSc logging stage 576 ierr = PetscLogStageRegister("Solve Stage", &rp->solvestage); 577 CHKERRQ(ierr); 578 579 rp->petschavecuda = petschavecuda; 580 rp->hostname = hostname; 581 rp->dim = dim; 582 rp->melem = melem; 583 rp->rankspernode = rankspernode; 584 585 for (PetscInt d = 0; d < num_degrees; d++) { 586 PetscInt deg = degree[d]; 587 for (PetscInt n = localnodes[0]; n < localnodes[1]; n *= 2) { 588 rp->degree = deg; 589 rp->localnodes = n; 590 ierr = Run(rp, num_ceedresources, ceedresources, 591 num_bpchoices, bpchoices); CHKERRQ(ierr); 592 } 593 } 594 // Clear memory 595 ierr = PetscFree(rp); CHKERRQ(ierr); 596 for (PetscInt i=0; i<num_ceedresources; i++) { 597 ierr = PetscFree(ceedresources[i]); CHKERRQ(ierr); 598 } 599 return PetscFinalize(); 600 } 601