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 is intentionally "raw", using only low-level communication 23 // primitives. 24 // 25 // Build with: 26 // 27 // make bps [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 28 // 29 // Sample runs: 30 // 31 // bps -problem bp1 32 // bps -problem bp2 -ceed /cpu/self 33 // bps -problem bp3 -ceed /gpu/occa 34 // bps -problem bp4 -ceed /cpu/occa 35 // bps -problem bp5 -ceed /omp/occa 36 // bps -problem bp6 -ceed /ocl/occa 37 // 38 //TESTARGS -ceed {ceed_resource} -test -problem bp3 39 40 /// @file 41 /// CEED BPs example using PETSc 42 /// See bpsdmplex.c for an implementation using DMPlex unstructured grids. 43 const char help[] = "Solve CEED BPs using PETSc\n"; 44 45 #include <stdbool.h> 46 #include <string.h> 47 #include "common.h" 48 #include "bp1.h" 49 #include "bp2.h" 50 #include "bp3.h" 51 #include "bp4.h" 52 53 #define PATH(BASE) __DIR__ #BASE 54 55 static void Split3(PetscInt size, PetscInt m[3], bool reverse) { 56 for (PetscInt d=0,sizeleft=size; d<3; d++) { 57 PetscInt try = (PetscInt)PetscCeilReal(PetscPowReal(sizeleft, 1./(3 - d))); 58 while (try * (sizeleft / try) != sizeleft) try++; 59 m[reverse ? 2-d : d] = try; 60 sizeleft /= try; 61 } 62 } 63 64 static PetscInt Max3(const PetscInt a[3]) { 65 return PetscMax(a[0], PetscMax(a[1], a[2])); 66 } 67 static PetscInt Min3(const PetscInt a[3]) { 68 return PetscMin(a[0], PetscMin(a[1], a[2])); 69 } 70 static void GlobalDof(const PetscInt p[3], const PetscInt irank[3], 71 PetscInt degree, const PetscInt melem[3], 72 PetscInt mdof[3]) { 73 for (int d=0; d<3; d++) 74 mdof[d] = degree*melem[d] + (irank[d] == p[d]-1); 75 } 76 static PetscInt GlobalStart(const PetscInt p[3], const PetscInt irank[3], 77 PetscInt degree, const PetscInt melem[3]) { 78 PetscInt start = 0; 79 // Dumb brute-force is easier to read 80 for (PetscInt i=0; i<p[0]; i++) { 81 for (PetscInt j=0; j<p[1]; j++) { 82 for (PetscInt k=0; k<p[2]; k++) { 83 PetscInt mdof[3], ijkrank[] = {i,j,k}; 84 if (i == irank[0] && j == irank[1] && k == irank[2]) return start; 85 GlobalDof(p, ijkrank, degree, melem, mdof); 86 start += mdof[0] * mdof[1] * mdof[2]; 87 } 88 } 89 } 90 return -1; 91 } 92 static int CreateRestriction(Ceed ceed, const CeedInt melem[3], 93 CeedInt P, CeedInt ncomp, 94 CeedElemRestriction *Erestrict) { 95 const PetscInt nelem = melem[0]*melem[1]*melem[2]; 96 PetscInt mdof[3], *idx, *idxp; 97 98 // Get indicies 99 for (int d=0; d<3; d++) mdof[d] = melem[d]*(P-1) + 1; 100 idxp = idx = malloc(nelem*P*P*P*sizeof idx[0]); 101 for (CeedInt i=0; i<melem[0]; i++) { 102 for (CeedInt j=0; j<melem[1]; j++) { 103 for (CeedInt k=0; k<melem[2]; k++,idxp += P*P*P) { 104 for (CeedInt ii=0; ii<P; ii++) { 105 for (CeedInt jj=0; jj<P; jj++) { 106 for (CeedInt kk=0; kk<P; kk++) { 107 if (0) { // This is the C-style (i,j,k) ordering that I prefer 108 idxp[(ii*P+jj)*P+kk] = (((i*(P-1)+ii)*mdof[1] 109 + (j*(P-1)+jj))*mdof[2] 110 + (k*(P-1)+kk)); 111 } else { // (k,j,i) ordering for consistency with MFEM example 112 idxp[ii+P*(jj+P*kk)] = (((i*(P-1)+ii)*mdof[1] 113 + (j*(P-1)+jj))*mdof[2] 114 + (k*(P-1)+kk)); 115 } 116 } 117 } 118 } 119 } 120 } 121 } 122 123 // Setup CEED restriction 124 CeedElemRestrictionCreate(ceed, nelem, P*P*P, mdof[0]*mdof[1]*mdof[2], ncomp, 125 CEED_MEM_HOST, CEED_OWN_POINTER, idx, Erestrict); 126 127 PetscFunctionReturn(0); 128 } 129 130 // Data for PETSc 131 typedef struct User_ *User; 132 struct User_ { 133 MPI_Comm comm; 134 VecScatter ltog; // Scatter for all entries 135 VecScatter ltog0; // Skip Dirichlet values 136 VecScatter gtogD; // global-to-global; only Dirichlet values 137 Vec Xloc, Yloc; 138 CeedVector xceed, yceed; 139 CeedOperator op; 140 CeedVector rho; 141 Ceed ceed; 142 }; 143 144 // BP Options 145 typedef enum { 146 CEED_BP1 = 0, CEED_BP2 = 1, CEED_BP3 = 2, 147 CEED_BP4 = 3, CEED_BP5 = 4, CEED_BP6 = 5 148 } bpType; 149 static const char *const bpTypes[] = {"bp1","bp2","bp3","bp4","bp5","bp6", 150 "bpType","CEED_BP",0}; 151 152 // BP specific data 153 typedef struct { 154 CeedInt vscale, qdatasize, qextra; 155 CeedQFunctionUser setup, apply, error; 156 const char setupfname[PETSC_MAX_PATH_LEN], applyfname[PETSC_MAX_PATH_LEN], 157 errorfname[PETSC_MAX_PATH_LEN]; 158 CeedEvalMode inmode, outmode; 159 CeedQuadMode qmode; 160 } bpData; 161 162 bpData bpOptions[6] = { 163 [CEED_BP1] = { 164 .vscale = 1, 165 .qdatasize = 1, 166 .qextra = 1, 167 .setup = SetupMass, 168 .apply = Mass, 169 .error = Error, 170 .setupfname = PATH(bp1.h:SetupMass), 171 .applyfname = PATH(bp1.h:Mass), 172 .errorfname = PATH(common.h:Error), 173 .inmode = CEED_EVAL_INTERP, 174 .outmode = CEED_EVAL_INTERP, 175 .qmode = CEED_GAUSS}, 176 [CEED_BP2] = { 177 .vscale = 3, 178 .qdatasize = 1, 179 .qextra = 1, 180 .setup = SetupMass3, 181 .apply = Mass3, 182 .error = Error3, 183 .setupfname = PATH(bp2.h:SetupMass3), 184 .applyfname = PATH(bp2.h:Mass3), 185 .errorfname = PATH(common.h:Error3), 186 .inmode = CEED_EVAL_INTERP, 187 .outmode = CEED_EVAL_INTERP, 188 .qmode = CEED_GAUSS}, 189 [CEED_BP3] = { 190 .vscale = 1, 191 .qdatasize = 6, 192 .qextra = 1, 193 .setup = SetupDiff, 194 .apply = Diff, 195 .error = Error, 196 .setupfname = PATH(bp3.h:SetupDiff), 197 .applyfname = PATH(bp3.h:Diff), 198 .errorfname = PATH(common.h:Error), 199 .inmode = CEED_EVAL_GRAD, 200 .outmode = CEED_EVAL_GRAD, 201 .qmode = CEED_GAUSS}, 202 [CEED_BP4] = { 203 .vscale = 3, 204 .qdatasize = 6, 205 .qextra = 1, 206 .setup = SetupDiff3, 207 .apply = Diff3, 208 .error = Error3, 209 .setupfname = PATH(bp4.h:SetupDiff3), 210 .applyfname = PATH(bp4.h:Diff3), 211 .errorfname = PATH(common.h:Error3), 212 .inmode = CEED_EVAL_GRAD, 213 .outmode = CEED_EVAL_GRAD, 214 .qmode = CEED_GAUSS}, 215 [CEED_BP5] = { 216 .vscale = 1, 217 .qdatasize = 6, 218 .qextra = 0, 219 .setup = SetupDiff, 220 .apply = Diff, 221 .error = Error, 222 .setupfname = PATH(bp3.h:SetupDiff), 223 .applyfname = PATH(bp3.h:Diff), 224 .errorfname = PATH(common.h:Error), 225 .inmode = CEED_EVAL_GRAD, 226 .outmode = CEED_EVAL_GRAD, 227 .qmode = CEED_GAUSS_LOBATTO}, 228 [CEED_BP6] = { 229 .vscale = 3, 230 .qdatasize = 6, 231 .qextra = 0, 232 .setup = SetupDiff3, 233 .apply = Diff3, 234 .error = Error3, 235 .setupfname = PATH(bp4.h:SetupDiff3), 236 .applyfname = PATH(bp4.h:Diff3), 237 .errorfname = PATH(common.h:Error3), 238 .inmode = CEED_EVAL_GRAD, 239 .outmode = CEED_EVAL_GRAD, 240 .qmode = CEED_GAUSS_LOBATTO} 241 }; 242 243 // This function uses libCEED to compute the action of the mass matrix 244 static PetscErrorCode MatMult_Mass(Mat A, Vec X, Vec Y) { 245 PetscErrorCode ierr; 246 User user; 247 PetscScalar *x, *y; 248 249 PetscFunctionBeginUser; 250 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 251 ierr = VecScatterBegin(user->ltog, X, user->Xloc, INSERT_VALUES, 252 SCATTER_REVERSE); CHKERRQ(ierr); 253 ierr = VecScatterEnd(user->ltog, X, user->Xloc, INSERT_VALUES, SCATTER_REVERSE); 254 CHKERRQ(ierr); 255 ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr); 256 257 ierr = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr); 258 ierr = VecGetArray(user->Yloc, &y); CHKERRQ(ierr); 259 CeedVectorSetArray(user->xceed, CEED_MEM_HOST, CEED_USE_POINTER, x); 260 CeedVectorSetArray(user->yceed, CEED_MEM_HOST, CEED_USE_POINTER, y); 261 262 CeedOperatorApply(user->op, user->xceed, user->yceed, 263 CEED_REQUEST_IMMEDIATE); 264 ierr = CeedVectorSyncArray(user->yceed, CEED_MEM_HOST); CHKERRQ(ierr); 265 266 ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr); 267 ierr = VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr); 268 269 if (Y) { 270 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 271 ierr = VecScatterBegin(user->ltog, user->Yloc, Y, ADD_VALUES, SCATTER_FORWARD); 272 CHKERRQ(ierr); 273 ierr = VecScatterEnd(user->ltog, user->Yloc, Y, ADD_VALUES, SCATTER_FORWARD); 274 CHKERRQ(ierr); 275 } 276 PetscFunctionReturn(0); 277 } 278 279 // This function uses libCEED to compute the action of the Laplacian with 280 // Dirichlet boundary conditions 281 static PetscErrorCode MatMult_Diff(Mat A, Vec X, Vec Y) { 282 PetscErrorCode ierr; 283 User user; 284 PetscScalar *x, *y; 285 286 PetscFunctionBeginUser; 287 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 288 289 // Global-to-local 290 ierr = VecScatterBegin(user->ltog0, X, user->Xloc, INSERT_VALUES, 291 SCATTER_REVERSE); CHKERRQ(ierr); 292 ierr = VecScatterEnd(user->ltog0, X, user->Xloc, INSERT_VALUES, 293 SCATTER_REVERSE); 294 CHKERRQ(ierr); 295 ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr); 296 297 // Setup CEED vectors 298 ierr = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr); 299 ierr = VecGetArray(user->Yloc, &y); CHKERRQ(ierr); 300 CeedVectorSetArray(user->xceed, CEED_MEM_HOST, CEED_USE_POINTER, x); 301 CeedVectorSetArray(user->yceed, CEED_MEM_HOST, CEED_USE_POINTER, y); 302 303 // Apply CEED operator 304 CeedOperatorApply(user->op, user->xceed, user->yceed, 305 CEED_REQUEST_IMMEDIATE); 306 ierr = CeedVectorSyncArray(user->yceed, CEED_MEM_HOST); CHKERRQ(ierr); 307 308 // Restore PETSc vectors 309 ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr); 310 ierr = VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr); 311 312 // Local-to-global 313 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 314 ierr = VecScatterBegin(user->gtogD, X, Y, INSERT_VALUES, SCATTER_FORWARD); 315 CHKERRQ(ierr); 316 ierr = VecScatterEnd(user->gtogD, X, Y, INSERT_VALUES, SCATTER_FORWARD); 317 CHKERRQ(ierr); 318 ierr = VecScatterBegin(user->ltog0, user->Yloc, Y, ADD_VALUES, SCATTER_FORWARD); 319 CHKERRQ(ierr); 320 ierr = VecScatterEnd(user->ltog0, user->Yloc, Y, ADD_VALUES, SCATTER_FORWARD); 321 CHKERRQ(ierr); 322 323 PetscFunctionReturn(0); 324 } 325 326 // This function calculates the error in the final solution 327 static PetscErrorCode ComputeErrorMax(User user, CeedOperator op_error, Vec X, 328 CeedVector target, PetscReal *maxerror) { 329 PetscErrorCode ierr; 330 PetscScalar *x; 331 CeedVector collocated_error; 332 CeedInt length; 333 334 PetscFunctionBeginUser; 335 CeedVectorGetLength(target, &length); 336 CeedVectorCreate(user->ceed, length, &collocated_error); 337 338 // Global-to-local 339 ierr = VecScatterBegin(user->ltog, X, user->Xloc, INSERT_VALUES, 340 SCATTER_REVERSE); CHKERRQ(ierr); 341 ierr = VecScatterEnd(user->ltog, X, user->Xloc, INSERT_VALUES, SCATTER_REVERSE); 342 CHKERRQ(ierr); 343 344 // Setup CEED vector 345 ierr = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr); 346 CeedVectorSetArray(user->xceed, CEED_MEM_HOST, CEED_USE_POINTER, x); 347 348 // Apply CEED operator 349 CeedOperatorApply(op_error, user->xceed, collocated_error, 350 CEED_REQUEST_IMMEDIATE); 351 352 // Restore PETSc vector 353 VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr); 354 355 // Reduce max error 356 *maxerror = 0; 357 const CeedScalar *e; 358 CeedVectorGetArrayRead(collocated_error, CEED_MEM_HOST, &e); 359 for (CeedInt i=0; i<length; i++) { 360 *maxerror = PetscMax(*maxerror, PetscAbsScalar(e[i])); 361 } 362 CeedVectorRestoreArrayRead(collocated_error, &e); 363 ierr = MPI_Allreduce(MPI_IN_PLACE, maxerror, 364 1, MPIU_REAL, MPIU_MAX, user->comm); CHKERRQ(ierr); 365 366 // Cleanup 367 CeedVectorDestroy(&collocated_error); 368 369 PetscFunctionReturn(0); 370 } 371 372 int main(int argc, char **argv) { 373 PetscInt ierr; 374 MPI_Comm comm; 375 char ceedresource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 376 double my_rt_start, my_rt, rt_min, rt_max; 377 PetscInt degree, qextra, localdof, localelem, melem[3], mdof[3], p[3], 378 irank[3], ldof[3], lsize, vscale = 1; 379 PetscScalar *r; 380 PetscBool test_mode, benchmark_mode; 381 PetscMPIInt size, rank; 382 Vec X, Xloc, rhs, rhsloc; 383 Mat mat; 384 KSP ksp; 385 VecScatter ltog, ltog0, gtogD; 386 User user; 387 Ceed ceed; 388 CeedBasis basisx, basisu; 389 CeedElemRestriction Erestrictx, Erestrictu, Erestrictxi, Erestrictui, 390 Erestrictqdi; 391 CeedQFunction qf_setup, qf_apply, qf_error; 392 CeedOperator op_setup, op_apply, op_error; 393 CeedVector xcoord, rho, rhsceed, target; 394 CeedInt P, Q; 395 bpType bpChoice; 396 397 ierr = PetscInitialize(&argc, &argv, NULL, help); 398 if (ierr) return ierr; 399 comm = PETSC_COMM_WORLD; 400 ierr = PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); CHKERRQ(ierr); 401 bpChoice = CEED_BP1; 402 ierr = PetscOptionsEnum("-problem", 403 "CEED benchmark problem to solve", NULL, 404 bpTypes, (PetscEnum)bpChoice, (PetscEnum*)&bpChoice, 405 NULL); CHKERRQ(ierr); 406 vscale = bpOptions[bpChoice].vscale; 407 test_mode = PETSC_FALSE; 408 ierr = PetscOptionsBool("-test", 409 "Testing mode (do not print unless error is large)", 410 NULL, test_mode, &test_mode, NULL); CHKERRQ(ierr); 411 benchmark_mode = PETSC_FALSE; 412 ierr = PetscOptionsBool("-benchmark", 413 "Benchmarking mode (prints benchmark statistics)", 414 NULL, benchmark_mode, &benchmark_mode, NULL); 415 CHKERRQ(ierr); 416 degree = test_mode ? 3 : 1; 417 ierr = PetscOptionsInt("-degree", "Polynomial degree of tensor product basis", 418 NULL, degree, °ree, NULL); CHKERRQ(ierr); 419 qextra = bpOptions[bpChoice].qextra; 420 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 421 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 422 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 423 NULL, ceedresource, ceedresource, 424 sizeof(ceedresource), NULL); CHKERRQ(ierr); 425 localdof = 1000; 426 ierr = PetscOptionsInt("-local", 427 "Target number of locally owned degrees of freedom per process", 428 NULL, localdof, &localdof, NULL); CHKERRQ(ierr); 429 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 430 P = degree + 1; 431 Q = P + qextra; 432 433 // Determine size of process grid 434 ierr = MPI_Comm_size(comm, &size); CHKERRQ(ierr); 435 Split3(size, p, false); 436 437 // Find a nicely composite number of elements no less than localdof 438 for (localelem = PetscMax(1, localdof / (degree*degree*degree)); ; 439 localelem++) { 440 Split3(localelem, melem, true); 441 if (Max3(melem) / Min3(melem) <= 2) break; 442 } 443 444 // Find my location in the process grid 445 ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr); 446 for (int d=0,rankleft=rank; d<3; d++) { 447 const int pstride[3] = {p[1] *p[2], p[2], 1}; 448 irank[d] = rankleft / pstride[d]; 449 rankleft -= irank[d] * pstride[d]; 450 } 451 452 GlobalDof(p, irank, degree, melem, mdof); 453 454 // Setup global vector 455 ierr = VecCreate(comm, &X); CHKERRQ(ierr); 456 ierr = VecSetSizes(X, mdof[0]*mdof[1]*mdof[2]*vscale, PETSC_DECIDE); 457 CHKERRQ(ierr); 458 ierr = VecSetUp(X); CHKERRQ(ierr); 459 460 // Print summary 461 if (!test_mode) { 462 CeedInt gsize; 463 ierr = VecGetSize(X, &gsize); CHKERRQ(ierr); 464 ierr = PetscPrintf(comm, 465 "\n-- CEED Benchmark Problem %d -- libCEED + PETSc --\n" 466 " libCEED:\n" 467 " libCEED Backend : %s\n" 468 " Mesh:\n" 469 " Number of 1D Basis Nodes (p) : %d\n" 470 " Number of 1D Quadrature Points (q) : %d\n" 471 " Global DOFs : %D\n" 472 " Process Decomposition : %D %D %D\n" 473 " Local Elements : %D = %D %D %D\n" 474 " Owned DOFs : %D = %D %D %D\n", 475 bpChoice+1, ceedresource, P, Q, gsize/vscale, p[0], 476 p[1], p[2], localelem, melem[0], melem[1], melem[2], 477 mdof[0]*mdof[1]*mdof[2], mdof[0], mdof[1], mdof[2]); 478 CHKERRQ(ierr); 479 } 480 481 { 482 lsize = 1; 483 for (int d=0; d<3; d++) { 484 ldof[d] = melem[d]*degree + 1; 485 lsize *= ldof[d]; 486 } 487 ierr = VecCreate(PETSC_COMM_SELF, &Xloc); CHKERRQ(ierr); 488 ierr = VecSetSizes(Xloc, lsize*vscale, PETSC_DECIDE); CHKERRQ(ierr); 489 ierr = VecSetUp(Xloc); CHKERRQ(ierr); 490 491 // Create local-to-global scatter 492 PetscInt *ltogind, *ltogind0, *locind, l0count; 493 IS ltogis, ltogis0, locis; 494 PetscInt gstart[2][2][2], gmdof[2][2][2][3]; 495 496 for (int i=0; i<2; i++) { 497 for (int j=0; j<2; j++) { 498 for (int k=0; k<2; k++) { 499 PetscInt ijkrank[3] = {irank[0]+i, irank[1]+j, irank[2]+k}; 500 gstart[i][j][k] = GlobalStart(p, ijkrank, degree, melem); 501 GlobalDof(p, ijkrank, degree, melem, gmdof[i][j][k]); 502 } 503 } 504 } 505 506 ierr = PetscMalloc1(lsize, <ogind); CHKERRQ(ierr); 507 ierr = PetscMalloc1(lsize, <ogind0); CHKERRQ(ierr); 508 ierr = PetscMalloc1(lsize, &locind); CHKERRQ(ierr); 509 l0count = 0; 510 for (PetscInt i=0,ir,ii; ir=i>=mdof[0], ii=i-ir*mdof[0], i<ldof[0]; i++) { 511 for (PetscInt j=0,jr,jj; jr=j>=mdof[1], jj=j-jr*mdof[1], j<ldof[1]; j++) { 512 for (PetscInt k=0,kr,kk; kr=k>=mdof[2], kk=k-kr*mdof[2], k<ldof[2]; k++) { 513 PetscInt here = (i*ldof[1]+j)*ldof[2]+k; 514 ltogind[here] = 515 gstart[ir][jr][kr] + (ii*gmdof[ir][jr][kr][1]+jj)*gmdof[ir][jr][kr][2]+kk; 516 if ((irank[0] == 0 && i == 0) 517 || (irank[1] == 0 && j == 0) 518 || (irank[2] == 0 && k == 0) 519 || (irank[0]+1 == p[0] && i+1 == ldof[0]) 520 || (irank[1]+1 == p[1] && j+1 == ldof[1]) 521 || (irank[2]+1 == p[2] && k+1 == ldof[2])) 522 continue; 523 ltogind0[l0count] = ltogind[here]; 524 locind[l0count++] = here; 525 } 526 } 527 } 528 ierr = ISCreateBlock(comm, vscale, lsize, ltogind, PETSC_OWN_POINTER, 529 <ogis); CHKERRQ(ierr); 530 ierr = VecScatterCreate(Xloc, NULL, X, ltogis, <og); CHKERRQ(ierr); 531 CHKERRQ(ierr); 532 ierr = ISCreateBlock(comm, vscale, l0count, ltogind0, PETSC_OWN_POINTER, 533 <ogis0); CHKERRQ(ierr); 534 ierr = ISCreateBlock(comm, vscale, l0count, locind, PETSC_OWN_POINTER, 535 &locis); CHKERRQ(ierr); 536 ierr = VecScatterCreate(Xloc, locis, X, ltogis0, <og0); CHKERRQ(ierr); 537 { 538 // Create global-to-global scatter for Dirichlet values (everything not in 539 // ltogis0, which is the range of ltog0) 540 PetscInt xstart, xend, *indD, countD = 0; 541 IS isD; 542 const PetscScalar *x; 543 ierr = VecZeroEntries(Xloc); CHKERRQ(ierr); 544 ierr = VecSet(X, 1.0); CHKERRQ(ierr); 545 ierr = VecScatterBegin(ltog0, Xloc, X, INSERT_VALUES, SCATTER_FORWARD); 546 CHKERRQ(ierr); 547 ierr = VecScatterEnd(ltog0, Xloc, X, INSERT_VALUES, SCATTER_FORWARD); 548 CHKERRQ(ierr); 549 ierr = VecGetOwnershipRange(X, &xstart, &xend); CHKERRQ(ierr); 550 ierr = PetscMalloc1(xend-xstart, &indD); CHKERRQ(ierr); 551 ierr = VecGetArrayRead(X, &x); CHKERRQ(ierr); 552 for (PetscInt i=0; i<xend-xstart; i++) { 553 if (x[i] == 1.) indD[countD++] = xstart + i; 554 } 555 ierr = VecRestoreArrayRead(X, &x); CHKERRQ(ierr); 556 ierr = ISCreateGeneral(comm, countD, indD, PETSC_COPY_VALUES, &isD); 557 CHKERRQ(ierr); 558 ierr = PetscFree(indD); CHKERRQ(ierr); 559 ierr = VecScatterCreate(X, isD, X, isD, >ogD); CHKERRQ(ierr); 560 ierr = ISDestroy(&isD); CHKERRQ(ierr); 561 } 562 ierr = ISDestroy(<ogis); CHKERRQ(ierr); 563 ierr = ISDestroy(<ogis0); CHKERRQ(ierr); 564 ierr = ISDestroy(&locis); CHKERRQ(ierr); 565 } 566 567 // Set up libCEED 568 CeedInit(ceedresource, &ceed); 569 CeedBasisCreateTensorH1Lagrange(ceed, 3, vscale, P, Q, 570 bpOptions[bpChoice].qmode, &basisu); 571 CeedBasisCreateTensorH1Lagrange(ceed, 3, 3, 2, Q, 572 bpOptions[bpChoice].qmode, &basisx); 573 574 CreateRestriction(ceed, melem, P, vscale, &Erestrictu); 575 CreateRestriction(ceed, melem, 2, 3, &Erestrictx); 576 CeedInt nelem = melem[0]*melem[1]*melem[2]; 577 CeedElemRestrictionCreateIdentity(ceed, nelem, Q*Q*Q, nelem*Q*Q*Q, vscale, 578 &Erestrictui); 579 CeedElemRestrictionCreateIdentity(ceed, nelem, 580 bpOptions[bpChoice].qdatasize*Q*Q*Q, 581 bpOptions[bpChoice].qdatasize*nelem*Q*Q*Q, 582 1, &Erestrictqdi); 583 CeedElemRestrictionCreateIdentity(ceed, nelem, Q*Q*Q, nelem*Q*Q*Q, 1, 584 &Erestrictxi); 585 { 586 CeedScalar *xloc; 587 CeedInt shape[3] = {melem[0]+1, melem[1]+1, melem[2]+1}, len = 588 shape[0]*shape[1]*shape[2]; 589 xloc = malloc(len*3*sizeof xloc[0]); 590 for (CeedInt i=0; i<shape[0]; i++) { 591 for (CeedInt j=0; j<shape[1]; j++) { 592 for (CeedInt k=0; k<shape[2]; k++) { 593 xloc[((i*shape[1]+j)*shape[2]+k) + 0*len] = 1.*(irank[0]*melem[0]+i) / 594 (p[0]*melem[0]); 595 xloc[((i*shape[1]+j)*shape[2]+k) + 1*len] = 1.*(irank[1]*melem[1]+j) / 596 (p[1]*melem[1]); 597 xloc[((i*shape[1]+j)*shape[2]+k) + 2*len] = 1.*(irank[2]*melem[2]+k) / 598 (p[2]*melem[2]); 599 } 600 } 601 } 602 CeedVectorCreate(ceed, len*3, &xcoord); 603 CeedVectorSetArray(xcoord, CEED_MEM_HOST, CEED_OWN_POINTER, xloc); 604 } 605 606 // Create the Q-function that builds the operator (i.e. computes its 607 // quadrature data) and set its context data 608 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].setup, 609 bpOptions[bpChoice].setupfname, &qf_setup); 610 CeedQFunctionAddInput(qf_setup, "x", 3, CEED_EVAL_INTERP); 611 CeedQFunctionAddInput(qf_setup, "dx", 3, CEED_EVAL_GRAD); 612 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 613 CeedQFunctionAddOutput(qf_setup, "rho", bpOptions[bpChoice].qdatasize, 614 CEED_EVAL_NONE); 615 CeedQFunctionAddOutput(qf_setup, "true_soln", vscale, CEED_EVAL_NONE); 616 CeedQFunctionAddOutput(qf_setup, "rhs", vscale, CEED_EVAL_INTERP); 617 618 // Set up PDE operator 619 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].apply, 620 bpOptions[bpChoice].applyfname, &qf_apply); 621 // Add inputs and outputs 622 CeedQFunctionAddInput(qf_apply, "u", vscale, bpOptions[bpChoice].inmode); 623 CeedQFunctionAddInput(qf_apply, "rho", bpOptions[bpChoice].qdatasize, 624 CEED_EVAL_NONE); 625 CeedQFunctionAddOutput(qf_apply, "v", vscale, bpOptions[bpChoice].outmode); 626 627 // Create the error qfunction 628 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].error, 629 bpOptions[bpChoice].errorfname, &qf_error); 630 CeedQFunctionAddInput(qf_error, "u", vscale, CEED_EVAL_INTERP); 631 CeedQFunctionAddInput(qf_error, "true_soln", vscale, CEED_EVAL_NONE); 632 CeedQFunctionAddOutput(qf_error, "error", vscale, CEED_EVAL_NONE); 633 634 // Create the persistent vectors that will be needed in setup 635 CeedInt nqpts; 636 CeedBasisGetNumQuadraturePoints(basisu, &nqpts); 637 CeedVectorCreate(ceed, bpOptions[bpChoice].qdatasize*nelem*nqpts, &rho); 638 CeedVectorCreate(ceed, nelem*nqpts*vscale, &target); 639 CeedVectorCreate(ceed, lsize*vscale, &rhsceed); 640 641 // Create the operator that builds the quadrature data for the ceed operator 642 CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup); 643 CeedOperatorSetField(op_setup, "x", Erestrictx, CEED_NOTRANSPOSE, 644 basisx, CEED_VECTOR_ACTIVE); 645 CeedOperatorSetField(op_setup, "dx", Erestrictx, CEED_NOTRANSPOSE, 646 basisx, CEED_VECTOR_ACTIVE); 647 CeedOperatorSetField(op_setup, "weight", Erestrictxi, CEED_NOTRANSPOSE, 648 basisx, CEED_VECTOR_NONE); 649 CeedOperatorSetField(op_setup, "rho", Erestrictqdi, CEED_NOTRANSPOSE, 650 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 651 CeedOperatorSetField(op_setup, "true_soln", Erestrictui, CEED_NOTRANSPOSE, 652 CEED_BASIS_COLLOCATED, target); 653 CeedOperatorSetField(op_setup, "rhs", Erestrictu, CEED_TRANSPOSE, 654 basisu, rhsceed); 655 656 // Create the mass or diff operator 657 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 658 CeedOperatorSetField(op_apply, "u", Erestrictu, CEED_TRANSPOSE, 659 basisu, CEED_VECTOR_ACTIVE); 660 CeedOperatorSetField(op_apply, "rho", Erestrictqdi, CEED_NOTRANSPOSE, 661 CEED_BASIS_COLLOCATED, rho); 662 CeedOperatorSetField(op_apply, "v", Erestrictu, CEED_TRANSPOSE, 663 basisu, CEED_VECTOR_ACTIVE); 664 665 // Create the error operator 666 CeedOperatorCreate(ceed, qf_error, NULL, NULL, &op_error); 667 CeedOperatorSetField(op_error, "u", Erestrictu, CEED_TRANSPOSE, 668 basisu, CEED_VECTOR_ACTIVE); 669 CeedOperatorSetField(op_error, "true_soln", Erestrictui, CEED_NOTRANSPOSE, 670 CEED_BASIS_COLLOCATED, target); 671 CeedOperatorSetField(op_error, "error", Erestrictui, CEED_NOTRANSPOSE, 672 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 673 674 675 // Set up Mat 676 ierr = PetscMalloc1(1, &user); CHKERRQ(ierr); 677 user->comm = comm; 678 user->ltog = ltog; 679 if (bpChoice != CEED_BP1 && bpChoice != CEED_BP2) { 680 user->ltog0 = ltog0; 681 user->gtogD = gtogD; 682 } 683 user->Xloc = Xloc; 684 ierr = VecDuplicate(Xloc, &user->Yloc); CHKERRQ(ierr); 685 CeedVectorCreate(ceed, lsize*vscale, &user->xceed); 686 CeedVectorCreate(ceed, lsize*vscale, &user->yceed); 687 user->op = op_apply; 688 user->rho = rho; 689 user->ceed = ceed; 690 691 ierr = MatCreateShell(comm, mdof[0]*mdof[1]*mdof[2]*vscale, 692 mdof[0]*mdof[1]*mdof[2]*vscale, 693 PETSC_DECIDE, PETSC_DECIDE, user, &mat); CHKERRQ(ierr); 694 if (bpChoice == CEED_BP1 || bpChoice == CEED_BP2) { 695 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Mass); 696 CHKERRQ(ierr); 697 } else { 698 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Diff); 699 CHKERRQ(ierr); 700 } 701 ierr = MatCreateVecs(mat, &rhs, NULL); CHKERRQ(ierr); 702 703 // Get RHS vector 704 ierr = VecDuplicate(Xloc, &rhsloc); CHKERRQ(ierr); 705 ierr = VecZeroEntries(rhsloc); CHKERRQ(ierr); 706 ierr = VecGetArray(rhsloc, &r); CHKERRQ(ierr); 707 CeedVectorSetArray(rhsceed, CEED_MEM_HOST, CEED_USE_POINTER, r); 708 709 // Setup rho, rhs, and target 710 CeedOperatorApply(op_setup, xcoord, rho, CEED_REQUEST_IMMEDIATE); 711 ierr = CeedVectorSyncArray(rhsceed, CEED_MEM_HOST); CHKERRQ(ierr); 712 CeedVectorDestroy(&xcoord); 713 714 // Gather RHS 715 ierr = VecRestoreArray(rhsloc, &r); CHKERRQ(ierr); 716 ierr = VecZeroEntries(rhs); CHKERRQ(ierr); 717 ierr = VecScatterBegin(ltog, rhsloc, rhs, ADD_VALUES, SCATTER_FORWARD); 718 CHKERRQ(ierr); 719 ierr = VecScatterEnd(ltog, rhsloc, rhs, ADD_VALUES, SCATTER_FORWARD); 720 CHKERRQ(ierr); 721 CeedVectorDestroy(&rhsceed); 722 723 ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr); 724 { 725 PC pc; 726 ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); 727 if (bpChoice == CEED_BP1 || bpChoice == CEED_BP2) { 728 ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr); 729 ierr = PCJacobiSetType(pc, PC_JACOBI_ROWSUM); CHKERRQ(ierr); 730 } else { 731 ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr); 732 } 733 ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr); 734 ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr); 735 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 736 PETSC_DEFAULT); CHKERRQ(ierr); 737 } 738 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); 739 ierr = KSPSetOperators(ksp, mat, mat); CHKERRQ(ierr); 740 // First run, if benchmarking 741 if (benchmark_mode) { 742 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1); 743 CHKERRQ(ierr); 744 my_rt_start = MPI_Wtime(); 745 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 746 my_rt = MPI_Wtime() - my_rt_start; 747 // Set maxits based on first iteration timing 748 if (my_rt > 0.02) { 749 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 5); 750 CHKERRQ(ierr); 751 } else { 752 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 20); 753 CHKERRQ(ierr); 754 } 755 } 756 // Timed solve 757 my_rt_start = MPI_Wtime(); 758 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 759 my_rt = MPI_Wtime() - my_rt_start; 760 { 761 KSPType ksptype; 762 KSPConvergedReason reason; 763 PetscReal rnorm; 764 PetscInt its; 765 ierr = KSPGetType(ksp, &ksptype); CHKERRQ(ierr); 766 ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr); 767 ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr); 768 ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr); 769 if (!test_mode || reason < 0 || rnorm > 1e-8) { 770 ierr = PetscPrintf(comm, 771 " KSP:\n" 772 " KSP Type : %s\n" 773 " KSP Convergence : %s\n" 774 " Total KSP Iterations : %D\n" 775 " Final rnorm : %e\n", 776 ksptype, KSPConvergedReasons[reason], its, 777 (double)rnorm); CHKERRQ(ierr); 778 } 779 if (benchmark_mode && (!test_mode)) { 780 CeedInt gsize; 781 ierr = VecGetSize(X, &gsize); CHKERRQ(ierr); 782 MPI_Reduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, 0, comm); 783 MPI_Reduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, 0, comm); 784 ierr = PetscPrintf(comm, 785 " Performance:\n" 786 " CG Solve Time : %g (%g) sec\n" 787 " DOFs/Sec in CG : %g (%g) million\n", 788 rt_max, rt_min, 1e-6*gsize*its/rt_max, 789 1e-6*gsize*its/rt_min); CHKERRQ(ierr); 790 } 791 } 792 793 { 794 PetscReal maxerror; 795 ierr = ComputeErrorMax(user, op_error, X, target, &maxerror); CHKERRQ(ierr); 796 PetscReal tol = (bpChoice == CEED_BP1 || bpChoice == CEED_BP2) ? 5e-3 : 5e-2; 797 if (!test_mode || maxerror > tol) { 798 ierr = PetscPrintf(comm, 799 " Pointwise Error (max) : %e\n", 800 (double)maxerror); CHKERRQ(ierr); 801 } 802 } 803 804 ierr = VecDestroy(&rhs); CHKERRQ(ierr); 805 ierr = VecDestroy(&rhsloc); CHKERRQ(ierr); 806 ierr = VecDestroy(&X); CHKERRQ(ierr); 807 ierr = VecDestroy(&user->Xloc); CHKERRQ(ierr); 808 ierr = VecDestroy(&user->Yloc); CHKERRQ(ierr); 809 ierr = VecScatterDestroy(<og); CHKERRQ(ierr); 810 ierr = VecScatterDestroy(<og0); CHKERRQ(ierr); 811 ierr = VecScatterDestroy(>ogD); CHKERRQ(ierr); 812 ierr = MatDestroy(&mat); CHKERRQ(ierr); 813 ierr = KSPDestroy(&ksp); CHKERRQ(ierr); 814 815 CeedVectorDestroy(&user->xceed); 816 CeedVectorDestroy(&user->yceed); 817 CeedVectorDestroy(&user->rho); 818 CeedVectorDestroy(&target); 819 CeedOperatorDestroy(&op_setup); 820 CeedOperatorDestroy(&op_apply); 821 CeedOperatorDestroy(&op_error); 822 CeedElemRestrictionDestroy(&Erestrictu); 823 CeedElemRestrictionDestroy(&Erestrictx); 824 CeedElemRestrictionDestroy(&Erestrictui); 825 CeedElemRestrictionDestroy(&Erestrictxi); 826 CeedElemRestrictionDestroy(&Erestrictqdi); 827 CeedQFunctionDestroy(&qf_setup); 828 CeedQFunctionDestroy(&qf_apply); 829 CeedQFunctionDestroy(&qf_error); 830 CeedBasisDestroy(&basisu); 831 CeedBasisDestroy(&basisx); 832 CeedDestroy(&ceed); 833 ierr = PetscFree(user); CHKERRQ(ierr); 834 return PetscFinalize(); 835 } 836