1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 // libCEED + PETSc Example: CEED BPs 9 // 10 // This example demonstrates a simple usage of libCEED with PETSc to solve the 11 // CEED BP benchmark problems, see http://ceed.exascaleproject.org/bps. 12 // 13 // The code is intentionally "raw", using only low-level communication 14 // primitives. 15 // 16 // Build with: 17 // 18 // make bpsraw [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 19 // 20 // Sample runs: 21 // 22 // ./bpsraw -problem bp1 23 // ./bpsraw -problem bp2 24 // ./bpsraw -problem bp3 25 // ./bpsraw -problem bp4 26 // ./bpsraw -problem bp5 -ceed /cpu/self 27 // ./bpsraw -problem bp6 -ceed /gpu/cuda 28 // 29 //TESTARGS -ceed {ceed_resource} -test -problem bp2 -degree 5 -q_extra 1 -ksp_max_it_clip 15,15 30 31 /// @file 32 /// CEED BPs example using PETSc 33 /// See bps.c for an implementation using DMPlex unstructured grids. 34 const char help[] = "Solve CEED BPs using PETSc\n"; 35 36 #include <ceed.h> 37 #include <petscksp.h> 38 #include <petscsys.h> 39 #include <stdbool.h> 40 #include <string.h> 41 #include "qfunctions/bps/bp1.h" 42 #include "qfunctions/bps/bp2.h" 43 #include "qfunctions/bps/bp3.h" 44 #include "qfunctions/bps/bp4.h" 45 #include "qfunctions/bps/common.h" 46 47 #if PETSC_VERSION_LT(3,12,0) 48 #ifdef PETSC_HAVE_CUDA 49 #include <petsccuda.h> 50 // Note: With PETSc prior to version 3.12.0, providing the source path to 51 // include 'cublas_v2.h' will be needed to use 'petsccuda.h'. 52 #endif 53 #endif 54 55 static CeedMemType MemTypeP2C(PetscMemType mem_type) { 56 return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 57 } 58 59 static void Split3(PetscInt size, PetscInt m[3], bool reverse) { 60 for (PetscInt d=0, size_left=size; d<3; d++) { 61 PetscInt try = (PetscInt)PetscCeilReal(PetscPowReal(size_left, 1./(3 - d))); 62 while (try * (size_left / try) != size_left) try++; 63 m[reverse ? 2-d : d] = try; 64 size_left /= try; 65 } 66 } 67 68 static PetscInt Max3(const PetscInt a[3]) { 69 return PetscMax(a[0], PetscMax(a[1], a[2])); 70 } 71 static PetscInt Min3(const PetscInt a[3]) { 72 return PetscMin(a[0], PetscMin(a[1], a[2])); 73 } 74 static void GlobalNodes(const PetscInt p[3], const PetscInt i_rank[3], 75 PetscInt degree, const PetscInt mesh_elem[3], 76 PetscInt m_nodes[3]) { 77 for (int d=0; d<3; d++) 78 m_nodes[d] = degree*mesh_elem[d] + (i_rank[d] == p[d]-1); 79 } 80 static PetscInt GlobalStart(const PetscInt p[3], const PetscInt i_rank[3], 81 PetscInt degree, const PetscInt mesh_elem[3]) { 82 PetscInt start = 0; 83 // Dumb brute-force is easier to read 84 for (PetscInt i=0; i<p[0]; i++) { 85 for (PetscInt j=0; j<p[1]; j++) { 86 for (PetscInt k=0; k<p[2]; k++) { 87 PetscInt m_nodes[3], ijk_rank[] = {i,j,k}; 88 if (i == i_rank[0] && j == i_rank[1] && k == i_rank[2]) return start; 89 GlobalNodes(p, ijk_rank, degree, mesh_elem, m_nodes); 90 start += m_nodes[0] * m_nodes[1] * m_nodes[2]; 91 } 92 } 93 } 94 return -1; 95 } 96 static PetscErrorCode CreateRestriction(Ceed ceed, const CeedInt mesh_elem[3], 97 CeedInt P, 98 CeedInt num_comp, CeedElemRestriction *elem_restr) { 99 const PetscInt num_elem = mesh_elem[0]*mesh_elem[1]*mesh_elem[2]; 100 PetscInt m_nodes[3], *idx, *idx_p; 101 102 PetscFunctionBeginUser; 103 // Get indicies 104 for (int d=0; d<3; d++) m_nodes[d] = mesh_elem[d]*(P-1) + 1; 105 idx_p = idx = malloc(num_elem*P*P*P*sizeof idx[0]); 106 for (CeedInt i=0; i<mesh_elem[0]; i++) 107 for (CeedInt j=0; j<mesh_elem[1]; j++) 108 for (CeedInt k=0; k<mesh_elem[2]; k++,idx_p += P*P*P) 109 for (CeedInt ii=0; ii<P; ii++) 110 for (CeedInt jj=0; jj<P; jj++) 111 for (CeedInt kk=0; kk<P; kk++) { 112 if (0) { // This is the C-style (i,j,k) ordering that I prefer 113 idx_p[(ii*P+jj)*P+kk] = num_comp*(((i*(P-1)+ii)*m_nodes[1] 114 + (j*(P-1)+jj))*m_nodes[2] 115 + (k*(P-1)+kk)); 116 } else { // (k,j,i) ordering for consistency with MFEM example 117 idx_p[ii+P*(jj+P*kk)] = num_comp*(((i*(P-1)+ii)*m_nodes[1] 118 + (j*(P-1)+jj))*m_nodes[2] 119 + (k*(P-1)+kk)); 120 } 121 } 122 123 // Setup CEED restriction 124 CeedElemRestrictionCreate(ceed, num_elem, P*P*P, num_comp, 1, 125 m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp, 126 CEED_MEM_HOST, CEED_OWN_POINTER, idx, elem_restr); 127 128 PetscFunctionReturn(0); 129 } 130 131 // Data for PETSc 132 typedef struct OperatorApplyContext_ *OperatorApplyContext; 133 struct OperatorApplyContext_ { 134 MPI_Comm comm; 135 VecScatter l_to_g; // Scatter for all entries 136 VecScatter l_to_g_0; // Skip Dirichlet values 137 VecScatter g_to_g_D; // global-to-global; only Dirichlet values 138 Vec X_loc, Y_loc; 139 CeedVector x_ceed, y_ceed; 140 CeedOperator op; 141 CeedVector q_data; 142 Ceed ceed; 143 }; 144 145 // BP Options 146 typedef enum { 147 CEED_BP1 = 0, CEED_BP2 = 1, CEED_BP3 = 2, 148 CEED_BP4 = 3, CEED_BP5 = 4, CEED_BP6 = 5 149 } BPType; 150 static const char *const bp_types[] = {"bp1","bp2","bp3","bp4","bp5","bp6", 151 "BPType","CEED_BP",0 152 }; 153 154 // BP specific data 155 typedef struct { 156 CeedInt num_comp_u, q_data_size, q_extra; 157 CeedQFunctionUser setup_geo, setup_rhs, apply, error; 158 const char *setup_geo_loc, *setup_rhs_loc, *apply_loc, *error_loc; 159 CeedEvalMode in_mode, out_mode; 160 CeedQuadMode q_mode; 161 } BPData; 162 163 BPData bp_options[6] = { 164 [CEED_BP1] = { 165 .num_comp_u = 1, 166 .q_data_size = 1, 167 .q_extra = 1, 168 .setup_geo = SetupMassGeo, 169 .setup_rhs = SetupMassRhs, 170 .apply = Mass, 171 .error = Error, 172 .setup_geo_loc = SetupMassGeo_loc, 173 .setup_rhs_loc = SetupMassRhs_loc, 174 .apply_loc = Mass_loc, 175 .error_loc = Error_loc, 176 .in_mode = CEED_EVAL_INTERP, 177 .out_mode = CEED_EVAL_INTERP, 178 .q_mode = CEED_GAUSS 179 }, 180 [CEED_BP2] = { 181 .num_comp_u = 3, 182 .q_data_size = 1, 183 .q_extra = 1, 184 .setup_geo = SetupMassGeo, 185 .setup_rhs = SetupMassRhs3, 186 .apply = Mass3, 187 .error = Error3, 188 .setup_geo_loc = SetupMassGeo_loc, 189 .setup_rhs_loc = SetupMassRhs3_loc, 190 .apply_loc = Mass3_loc, 191 .error_loc = Error3_loc, 192 .in_mode = CEED_EVAL_INTERP, 193 .out_mode = CEED_EVAL_INTERP, 194 .q_mode = CEED_GAUSS 195 }, 196 [CEED_BP3] = { 197 .num_comp_u = 1, 198 .q_data_size = 7, 199 .q_extra = 1, 200 .setup_geo = SetupDiffGeo, 201 .setup_rhs = SetupDiffRhs, 202 .apply = Diff, 203 .error = Error, 204 .setup_geo_loc = SetupDiffGeo_loc, 205 .setup_rhs_loc = SetupDiffRhs_loc, 206 .apply_loc = Diff_loc, 207 .error_loc = Error_loc, 208 .in_mode = CEED_EVAL_GRAD, 209 .out_mode = CEED_EVAL_GRAD, 210 .q_mode = CEED_GAUSS 211 }, 212 [CEED_BP4] = { 213 .num_comp_u = 3, 214 .q_data_size = 7, 215 .q_extra = 1, 216 .setup_geo = SetupDiffGeo, 217 .setup_rhs = SetupDiffRhs3, 218 .apply = Diff3, 219 .error = Error3, 220 .setup_geo_loc = SetupDiffGeo_loc, 221 .setup_rhs_loc = SetupDiffRhs3_loc, 222 .apply_loc = Diff3_loc, 223 .error_loc = Error3_loc, 224 .in_mode = CEED_EVAL_GRAD, 225 .out_mode = CEED_EVAL_GRAD, 226 .q_mode = CEED_GAUSS 227 }, 228 [CEED_BP5] = { 229 .num_comp_u = 1, 230 .q_data_size = 7, 231 .q_extra = 0, 232 .setup_geo = SetupDiffGeo, 233 .setup_rhs = SetupDiffRhs, 234 .apply = Diff, 235 .error = Error, 236 .setup_geo_loc = SetupDiffGeo_loc, 237 .setup_rhs_loc = SetupDiffRhs_loc, 238 .apply_loc = Diff_loc, 239 .error_loc = Error_loc, 240 .in_mode = CEED_EVAL_GRAD, 241 .out_mode = CEED_EVAL_GRAD, 242 .q_mode = CEED_GAUSS_LOBATTO 243 }, 244 [CEED_BP6] = { 245 .num_comp_u = 3, 246 .q_data_size = 7, 247 .q_extra = 0, 248 .setup_geo = SetupDiffGeo, 249 .setup_rhs = SetupDiffRhs3, 250 .apply = Diff3, 251 .error = Error3, 252 .setup_geo_loc = SetupDiffGeo_loc, 253 .setup_rhs_loc = SetupDiffRhs3_loc, 254 .apply_loc = Diff3_loc, 255 .error_loc = Error3_loc, 256 .in_mode = CEED_EVAL_GRAD, 257 .out_mode = CEED_EVAL_GRAD, 258 .q_mode = CEED_GAUSS_LOBATTO 259 } 260 }; 261 262 // This function uses libCEED to compute the action of the mass matrix 263 static PetscErrorCode MatMult_Mass(Mat A, Vec X, Vec Y) { 264 PetscErrorCode ierr; 265 OperatorApplyContext op_apply_ctx; 266 PetscScalar *x, *y; 267 PetscMemType x_mem_type, y_mem_type; 268 269 PetscFunctionBeginUser; 270 271 ierr = MatShellGetContext(A, &op_apply_ctx); CHKERRQ(ierr); 272 273 // Global-to-local 274 ierr = VecScatterBegin(op_apply_ctx->l_to_g, X, op_apply_ctx->X_loc, 275 INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 276 ierr = VecScatterEnd(op_apply_ctx->l_to_g, X, op_apply_ctx->X_loc, 277 INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 278 279 // Setup libCEED vectors 280 ierr = VecGetArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x, 281 &x_mem_type); CHKERRQ(ierr); 282 ierr = VecGetArrayAndMemType(op_apply_ctx->Y_loc, &y, &y_mem_type); 283 CHKERRQ(ierr); 284 CeedVectorSetArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), 285 CEED_USE_POINTER, x); 286 CeedVectorSetArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), 287 CEED_USE_POINTER, y); 288 289 // Apply libCEED operator 290 CeedOperatorApply(op_apply_ctx->op, op_apply_ctx->x_ceed, op_apply_ctx->y_ceed, 291 CEED_REQUEST_IMMEDIATE); 292 293 // Restore PETSc vectors 294 CeedVectorTakeArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), NULL); 295 CeedVectorTakeArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), NULL); 296 ierr = VecRestoreArrayReadAndMemType(op_apply_ctx->X_loc, 297 (const PetscScalar **)&x); CHKERRQ(ierr); 298 ierr = VecRestoreArrayAndMemType(op_apply_ctx->Y_loc, &y); CHKERRQ(ierr); 299 300 // Local-to-global 301 if (Y) { 302 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 303 ierr = VecScatterBegin(op_apply_ctx->l_to_g, op_apply_ctx->Y_loc, Y, ADD_VALUES, 304 SCATTER_FORWARD); CHKERRQ(ierr); 305 ierr = VecScatterEnd(op_apply_ctx->l_to_g, op_apply_ctx->Y_loc, Y, ADD_VALUES, 306 SCATTER_FORWARD); CHKERRQ(ierr); 307 } 308 PetscFunctionReturn(0); 309 } 310 311 // This function uses libCEED to compute the action of the Laplacian with 312 // Dirichlet boundary conditions 313 static PetscErrorCode MatMult_Diff(Mat A, Vec X, Vec Y) { 314 PetscErrorCode ierr; 315 OperatorApplyContext op_apply_ctx; 316 PetscScalar *x, *y; 317 PetscMemType x_mem_type, y_mem_type; 318 319 PetscFunctionBeginUser; 320 321 ierr = MatShellGetContext(A, &op_apply_ctx); CHKERRQ(ierr); 322 323 // Global-to-local 324 ierr = VecScatterBegin(op_apply_ctx->l_to_g_0, X, op_apply_ctx->X_loc, 325 INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 326 ierr = VecScatterEnd(op_apply_ctx->l_to_g_0, X, op_apply_ctx->X_loc, 327 INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 328 329 // Setup libCEED vectors 330 ierr = VecGetArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x, 331 &x_mem_type); CHKERRQ(ierr); 332 ierr = VecGetArrayAndMemType(op_apply_ctx->Y_loc, &y, &y_mem_type); 333 CHKERRQ(ierr); 334 CeedVectorSetArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), 335 CEED_USE_POINTER, x); 336 CeedVectorSetArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), 337 CEED_USE_POINTER, y); 338 339 // Apply libCEED operator 340 CeedOperatorApply(op_apply_ctx->op, op_apply_ctx->x_ceed, op_apply_ctx->y_ceed, 341 CEED_REQUEST_IMMEDIATE); 342 343 // Restore PETSc vectors 344 CeedVectorTakeArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), NULL); 345 CeedVectorTakeArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), NULL); 346 ierr = VecRestoreArrayReadAndMemType(op_apply_ctx->X_loc, 347 (const PetscScalar **)&x); 348 CHKERRQ(ierr); 349 ierr = VecRestoreArrayAndMemType(op_apply_ctx->Y_loc, &y); CHKERRQ(ierr); 350 351 // Local-to-global 352 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 353 ierr = VecScatterBegin(op_apply_ctx->g_to_g_D, X, Y, INSERT_VALUES, 354 SCATTER_FORWARD); 355 CHKERRQ(ierr); 356 ierr = VecScatterEnd(op_apply_ctx->g_to_g_D, X, Y, INSERT_VALUES, 357 SCATTER_FORWARD); CHKERRQ(ierr); 358 ierr = VecScatterBegin(op_apply_ctx->l_to_g_0, op_apply_ctx->Y_loc, Y, 359 ADD_VALUES, SCATTER_FORWARD); CHKERRQ(ierr); 360 ierr = VecScatterEnd(op_apply_ctx->l_to_g_0, op_apply_ctx->Y_loc, Y, ADD_VALUES, 361 SCATTER_FORWARD); 362 CHKERRQ(ierr); 363 364 PetscFunctionReturn(0); 365 } 366 367 // This function calculates the error in the final solution 368 static PetscErrorCode ComputeErrorMax(OperatorApplyContext op_apply_ctx, 369 CeedOperator op_error, Vec X, 370 CeedVector target, PetscReal *max_error) { 371 PetscErrorCode ierr; 372 PetscScalar *x; 373 PetscMemType mem_type; 374 CeedVector collocated_error; 375 CeedSize length; 376 377 PetscFunctionBeginUser; 378 379 CeedVectorGetLength(target, &length); 380 CeedVectorCreate(op_apply_ctx->ceed, length, &collocated_error); 381 382 // Global-to-local 383 ierr = VecScatterBegin(op_apply_ctx->l_to_g, X, op_apply_ctx->X_loc, 384 INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 385 ierr = VecScatterEnd(op_apply_ctx->l_to_g, X, op_apply_ctx->X_loc, 386 INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 387 388 // Setup libCEED vector 389 ierr = VecGetArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x, 390 &mem_type); CHKERRQ(ierr); 391 CeedVectorSetArray(op_apply_ctx->x_ceed, MemTypeP2C(mem_type), 392 CEED_USE_POINTER, x); 393 394 // Apply libCEED operator 395 CeedOperatorApply(op_error, op_apply_ctx->x_ceed, collocated_error, 396 CEED_REQUEST_IMMEDIATE); 397 398 // Restore PETSc vector 399 CeedVectorTakeArray(op_apply_ctx->x_ceed, MemTypeP2C(mem_type), NULL); 400 ierr = VecRestoreArrayReadAndMemType(op_apply_ctx->X_loc, 401 (const PetscScalar **)&x); CHKERRQ(ierr); 402 403 // Reduce max error 404 *max_error = 0; 405 const CeedScalar *e; 406 CeedVectorGetArrayRead(collocated_error, CEED_MEM_HOST, &e); 407 for (CeedInt i=0; i<length; i++) { 408 *max_error = PetscMax(*max_error, PetscAbsScalar(e[i])); 409 } 410 CeedVectorRestoreArrayRead(collocated_error, &e); 411 ierr = MPI_Allreduce(MPI_IN_PLACE, max_error, 1, MPIU_REAL, MPIU_MAX, 412 op_apply_ctx->comm); CHKERRQ(ierr); 413 414 // Cleanup 415 CeedVectorDestroy(&collocated_error); 416 417 PetscFunctionReturn(0); 418 } 419 420 int main(int argc, char **argv) { 421 PetscInt ierr; 422 MPI_Comm comm; 423 char ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 424 double my_rt_start, my_rt, rt_min, rt_max; 425 PetscInt degree, q_extra, local_nodes, local_elem, mesh_elem[3], m_nodes[3], 426 p[3], 427 i_rank[3], l_nodes[3], l_size, num_comp_u = 1, ksp_max_it_clip[2]; 428 PetscScalar *r; 429 PetscBool test_mode, benchmark_mode, write_solution; 430 PetscMPIInt size, rank; 431 PetscLogStage solve_stage; 432 Vec X, X_loc, rhs, rhs_loc; 433 Mat mat; 434 KSP ksp; 435 VecScatter l_to_g, l_to_g_0, g_to_g_D; 436 PetscMemType mem_type; 437 OperatorApplyContext op_apply_ctx; 438 Ceed ceed; 439 CeedBasis basis_x, basis_u; 440 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i; 441 CeedQFunction qf_setup_geo, qf_setup_rhs, qf_apply, qf_error; 442 CeedOperator op_setup_geo, op_setup_rhs, op_apply, op_error; 443 CeedVector x_coord, q_data, rhs_ceed, target; 444 CeedInt P, Q; 445 const CeedInt dim = 3, num_comp_x = 3; 446 BPType bp_choice; 447 448 ierr = PetscInitialize(&argc, &argv, NULL, help); 449 if (ierr) return ierr; 450 comm = PETSC_COMM_WORLD; 451 452 // Read command line options 453 PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); 454 bp_choice = CEED_BP1; 455 ierr = PetscOptionsEnum("-problem", 456 "CEED benchmark problem to solve", NULL, 457 bp_types, (PetscEnum)bp_choice, (PetscEnum *)&bp_choice, 458 NULL); CHKERRQ(ierr); 459 num_comp_u = bp_options[bp_choice].num_comp_u; 460 test_mode = PETSC_FALSE; 461 ierr = PetscOptionsBool("-test", 462 "Testing mode (do not print unless error is large)", 463 NULL, test_mode, &test_mode, NULL); CHKERRQ(ierr); 464 benchmark_mode = PETSC_FALSE; 465 ierr = PetscOptionsBool("-benchmark", 466 "Benchmarking mode (prints benchmark statistics)", 467 NULL, benchmark_mode, &benchmark_mode, NULL); 468 CHKERRQ(ierr); 469 write_solution = PETSC_FALSE; 470 ierr = PetscOptionsBool("-write_solution", 471 "Write solution for visualization", 472 NULL, write_solution, &write_solution, NULL); 473 CHKERRQ(ierr); 474 degree = test_mode ? 3 : 1; 475 ierr = PetscOptionsInt("-degree", "Polynomial degree of tensor product basis", 476 NULL, degree, °ree, NULL); CHKERRQ(ierr); 477 q_extra = bp_options[bp_choice].q_extra; 478 ierr = PetscOptionsInt("-q_extra", "Number of extra quadrature points", 479 NULL, q_extra, &q_extra, NULL); CHKERRQ(ierr); 480 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 481 NULL, ceed_resource, ceed_resource, 482 sizeof(ceed_resource), NULL); CHKERRQ(ierr); 483 local_nodes = 1000; 484 ierr = PetscOptionsInt("-local", 485 "Target number of locally owned nodes per process", 486 NULL, local_nodes, &local_nodes, NULL); CHKERRQ(ierr); 487 PetscInt two = 2; 488 ksp_max_it_clip[0] = 5; 489 ksp_max_it_clip[1] = 20; 490 ierr = PetscOptionsIntArray("-ksp_max_it_clip", 491 "Min and max number of iterations to use during benchmarking", 492 NULL, ksp_max_it_clip, &two, NULL); 493 CHKERRQ(ierr); 494 PetscOptionsEnd(); 495 P = degree + 1; 496 Q = P + q_extra; 497 498 // Set up libCEED 499 CeedInit(ceed_resource, &ceed); 500 CeedMemType mem_type_backend; 501 CeedGetPreferredMemType(ceed, &mem_type_backend); 502 503 VecType default_vec_type = NULL, vec_type; 504 switch (mem_type_backend) { 505 case CEED_MEM_HOST: default_vec_type = VECSTANDARD; break; 506 case CEED_MEM_DEVICE: { 507 const char *resolved; 508 CeedGetResource(ceed, &resolved); 509 if (strstr(resolved, "/gpu/cuda")) default_vec_type = VECCUDA; 510 else if (strstr(resolved, "/gpu/hip/occa")) 511 default_vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 512 else if (strstr(resolved, "/gpu/hip")) default_vec_type = VECHIP; 513 else default_vec_type = VECSTANDARD; 514 } 515 } 516 517 // Determine size of process grid 518 ierr = MPI_Comm_size(comm, &size); CHKERRQ(ierr); 519 Split3(size, p, false); 520 521 // Find a nicely composite number of elements no less than local_nodes 522 for (local_elem = PetscMax(1, local_nodes / (degree*degree*degree)); ; 523 local_elem++) { 524 Split3(local_elem, mesh_elem, true); 525 if (Max3(mesh_elem) / Min3(mesh_elem) <= 2) break; 526 } 527 528 // Find my location in the process grid 529 ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr); 530 for (int d=0, rank_left=rank; d<dim; d++) { 531 const int pstride[3] = {p[1] *p[2], p[2], 1}; 532 i_rank[d] = rank_left / pstride[d]; 533 rank_left -= i_rank[d] * pstride[d]; 534 } 535 536 GlobalNodes(p, i_rank, degree, mesh_elem, m_nodes); 537 538 // Setup global vector 539 ierr = VecCreate(comm, &X); CHKERRQ(ierr); 540 ierr = VecSetType(X, default_vec_type); CHKERRQ(ierr); 541 ierr = VecSetSizes(X, m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 542 PETSC_DECIDE); 543 CHKERRQ(ierr); 544 ierr = VecSetFromOptions(X); CHKERRQ(ierr); 545 ierr = VecSetUp(X); CHKERRQ(ierr); 546 547 // Set up libCEED 548 CeedInit(ceed_resource, &ceed); 549 550 // Print summary 551 CeedInt gsize; 552 ierr = VecGetSize(X, &gsize); CHKERRQ(ierr); 553 if (!test_mode) { 554 const char *used_resource; 555 CeedGetResource(ceed, &used_resource); 556 557 ierr = VecGetType(X, &vec_type); CHKERRQ(ierr); 558 559 ierr = PetscPrintf(comm, 560 "\n-- CEED Benchmark Problem %" CeedInt_FMT " -- libCEED + PETSc --\n" 561 " PETSc:\n" 562 " PETSc Vec Type : %s\n" 563 " libCEED:\n" 564 " libCEED Backend : %s\n" 565 " libCEED Backend MemType : %s\n" 566 " Mesh:\n" 567 " Solution Order (P) : %" CeedInt_FMT "\n" 568 " Quadrature Order (Q) : %" CeedInt_FMT "\n" 569 " Global nodes : %" PetscInt_FMT "\n" 570 " Process Decomposition : %" PetscInt_FMT 571 " %" PetscInt_FMT " %" PetscInt_FMT "\n" 572 " Local Elements : %" PetscInt_FMT 573 " = %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n" 574 " Owned nodes : %" PetscInt_FMT 575 " = %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n" 576 " DoF per node : %" PetscInt_FMT "\n", 577 bp_choice+1, vec_type, used_resource, 578 CeedMemTypes[mem_type_backend], 579 P, Q, gsize/num_comp_u, p[0], p[1], p[2], local_elem, 580 mesh_elem[0], mesh_elem[1], mesh_elem[2], 581 m_nodes[0]*m_nodes[1]*m_nodes[2], m_nodes[0], m_nodes[1], 582 m_nodes[2], num_comp_u); CHKERRQ(ierr); 583 } 584 585 { 586 l_size = 1; 587 for (int d=0; d<dim; d++) { 588 l_nodes[d] = mesh_elem[d]*degree + 1; 589 l_size *= l_nodes[d]; 590 } 591 ierr = VecCreate(PETSC_COMM_SELF, &X_loc); CHKERRQ(ierr); 592 ierr = VecSetType(X_loc, default_vec_type); CHKERRQ(ierr); 593 ierr = VecSetSizes(X_loc, l_size*num_comp_u, PETSC_DECIDE); CHKERRQ(ierr); 594 ierr = VecSetFromOptions(X_loc); CHKERRQ(ierr); 595 ierr = VecSetUp(X_loc); CHKERRQ(ierr); 596 597 // Create local-to-global scatter 598 PetscInt *l_to_g_ind, *l_to_g_ind_0, *loc_ind, l_0_count; 599 IS l_to_g_is, l_to_g_is_0, loc_is; 600 PetscInt g_start[2][2][2], g_m_nodes[2][2][2][dim]; 601 602 for (int i=0; i<2; i++) { 603 for (int j=0; j<2; j++) { 604 for (int k=0; k<2; k++) { 605 PetscInt ijk_rank[3] = {i_rank[0]+i, i_rank[1]+j, i_rank[2]+k}; 606 g_start[i][j][k] = GlobalStart(p, ijk_rank, degree, mesh_elem); 607 GlobalNodes(p, ijk_rank, degree, mesh_elem, g_m_nodes[i][j][k]); 608 } 609 } 610 } 611 612 ierr = PetscMalloc1(l_size, &l_to_g_ind); CHKERRQ(ierr); 613 ierr = PetscMalloc1(l_size, &l_to_g_ind_0); CHKERRQ(ierr); 614 ierr = PetscMalloc1(l_size, &loc_ind); CHKERRQ(ierr); 615 l_0_count = 0; 616 for (PetscInt i=0,ir,ii; ir=i>=m_nodes[0], ii=i-ir*m_nodes[0], i<l_nodes[0]; 617 i++) 618 for (PetscInt j=0,jr,jj; jr=j>=m_nodes[1], jj=j-jr*m_nodes[1], j<l_nodes[1]; 619 j++) 620 for (PetscInt k=0,kr,kk; kr=k>=m_nodes[2], kk=k-kr*m_nodes[2], k<l_nodes[2]; 621 k++) { 622 PetscInt here = (i*l_nodes[1]+j)*l_nodes[2]+k; 623 l_to_g_ind[here] = 624 g_start[ir][jr][kr] + (ii*g_m_nodes[ir][jr][kr][1]+jj)*g_m_nodes[ir][jr][kr][2] 625 +kk; 626 if ((i_rank[0] == 0 && i == 0) 627 || (i_rank[1] == 0 && j == 0) 628 || (i_rank[2] == 0 && k == 0) 629 || (i_rank[0]+1 == p[0] && i+1 == l_nodes[0]) 630 || (i_rank[1]+1 == p[1] && j+1 == l_nodes[1]) 631 || (i_rank[2]+1 == p[2] && k+1 == l_nodes[2])) 632 continue; 633 l_to_g_ind_0[l_0_count] = l_to_g_ind[here]; 634 loc_ind[l_0_count++] = here; 635 } 636 ierr = ISCreateBlock(comm, num_comp_u, l_size, l_to_g_ind, PETSC_OWN_POINTER, 637 &l_to_g_is); CHKERRQ(ierr); 638 ierr = VecScatterCreate(X_loc, NULL, X, l_to_g_is, &l_to_g); CHKERRQ(ierr); 639 CHKERRQ(ierr); 640 ierr = ISCreateBlock(comm, num_comp_u, l_0_count, l_to_g_ind_0, 641 PETSC_OWN_POINTER, 642 &l_to_g_is_0); CHKERRQ(ierr); 643 ierr = ISCreateBlock(comm, num_comp_u, l_0_count, loc_ind, PETSC_OWN_POINTER, 644 &loc_is); CHKERRQ(ierr); 645 ierr = VecScatterCreate(X_loc, loc_is, X, l_to_g_is_0, &l_to_g_0); 646 CHKERRQ(ierr); 647 { 648 // Create global-to-global scatter for Dirichlet values (everything not in 649 // l_to_g_is_0, which is the range of l_to_g_0) 650 PetscInt x_start, x_end, *ind_D, count_D = 0; 651 IS is_D; 652 const PetscScalar *x; 653 ierr = VecZeroEntries(X_loc); CHKERRQ(ierr); 654 ierr = VecSet(X, 1.0); CHKERRQ(ierr); 655 ierr = VecScatterBegin(l_to_g_0, X_loc, X, INSERT_VALUES, SCATTER_FORWARD); 656 CHKERRQ(ierr); 657 ierr = VecScatterEnd(l_to_g_0, X_loc, X, INSERT_VALUES, SCATTER_FORWARD); 658 CHKERRQ(ierr); 659 ierr = VecGetOwnershipRange(X, &x_start, &x_end); CHKERRQ(ierr); 660 ierr = PetscMalloc1(x_end-x_start, &ind_D); CHKERRQ(ierr); 661 ierr = VecGetArrayRead(X, &x); CHKERRQ(ierr); 662 for (PetscInt i=0; i<x_end-x_start; i++) { 663 if (x[i] == 1.) ind_D[count_D++] = x_start + i; 664 } 665 ierr = VecRestoreArrayRead(X, &x); CHKERRQ(ierr); 666 ierr = ISCreateGeneral(comm, count_D, ind_D, PETSC_COPY_VALUES, &is_D); 667 CHKERRQ(ierr); 668 ierr = PetscFree(ind_D); CHKERRQ(ierr); 669 ierr = VecScatterCreate(X, is_D, X, is_D, &g_to_g_D); CHKERRQ(ierr); 670 ierr = ISDestroy(&is_D); CHKERRQ(ierr); 671 } 672 ierr = ISDestroy(&l_to_g_is); CHKERRQ(ierr); 673 ierr = ISDestroy(&l_to_g_is_0); CHKERRQ(ierr); 674 ierr = ISDestroy(&loc_is); CHKERRQ(ierr); 675 } 676 677 // CEED bases 678 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, 679 bp_options[bp_choice].q_mode, &basis_u); 680 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, 681 bp_options[bp_choice].q_mode, &basis_x); 682 683 // CEED restrictions 684 CreateRestriction(ceed, mesh_elem, P, num_comp_u, &elem_restr_u); 685 CreateRestriction(ceed, mesh_elem, 2, dim, &elem_restr_x); 686 CeedInt num_elem = mesh_elem[0]*mesh_elem[1]*mesh_elem[2]; 687 CeedElemRestrictionCreateStrided(ceed, num_elem, Q*Q*Q, num_comp_u, 688 num_comp_u*num_elem*Q*Q*Q, 689 CEED_STRIDES_BACKEND, &elem_restr_u_i); 690 CeedElemRestrictionCreateStrided(ceed, num_elem, Q*Q*Q, 691 bp_options[bp_choice].q_data_size, 692 bp_options[bp_choice].q_data_size*num_elem*Q*Q*Q, 693 CEED_STRIDES_BACKEND, &elem_restr_qd_i); 694 { 695 CeedScalar *x_loc; 696 CeedInt shape[3] = {mesh_elem[0]+1, mesh_elem[1]+1, mesh_elem[2]+1}, len = 697 shape[0]*shape[1]*shape[2]; 698 x_loc = malloc(len*num_comp_x*sizeof x_loc[0]); 699 for (CeedInt i=0; i<shape[0]; i++) { 700 for (CeedInt j=0; j<shape[1]; j++) { 701 for (CeedInt k=0; k<shape[2]; k++) { 702 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 0] = 1.*(i_rank[0]*mesh_elem[0]+i) / 703 (p[0]*mesh_elem[0]); 704 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 1] = 1.*(i_rank[1]*mesh_elem[1]+j) / 705 (p[1]*mesh_elem[1]); 706 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 2] = 1.*(i_rank[2]*mesh_elem[2]+k) / 707 (p[2]*mesh_elem[2]); 708 } 709 } 710 } 711 CeedVectorCreate(ceed, len*num_comp_x, &x_coord); 712 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_OWN_POINTER, x_loc); 713 } 714 715 // Create the QFunction that builds the operator quadrature data 716 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].setup_geo, 717 bp_options[bp_choice].setup_geo_loc, &qf_setup_geo); 718 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 719 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x*dim, CEED_EVAL_GRAD); 720 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 721 CeedQFunctionAddOutput(qf_setup_geo, "q_data", 722 bp_options[bp_choice].q_data_size, 723 CEED_EVAL_NONE); 724 725 // Create the QFunction that sets up the RHS and true solution 726 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].setup_rhs, 727 bp_options[bp_choice].setup_rhs_loc, &qf_setup_rhs); 728 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 729 CeedQFunctionAddInput(qf_setup_rhs, "q_data", bp_options[bp_choice].q_data_size, 730 CEED_EVAL_NONE); 731 CeedQFunctionAddOutput(qf_setup_rhs, "true_soln", num_comp_u, CEED_EVAL_NONE); 732 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 733 734 // Set up PDE operator 735 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].apply, 736 bp_options[bp_choice].apply_loc, &qf_apply); 737 // Add inputs and outputs 738 CeedInt in_scale = bp_options[bp_choice].in_mode==CEED_EVAL_GRAD ? 3 : 1; 739 CeedInt out_scale = bp_options[bp_choice].out_mode==CEED_EVAL_GRAD ? 3 : 1; 740 CeedQFunctionAddInput(qf_apply, "u", num_comp_u*in_scale, 741 bp_options[bp_choice].in_mode); 742 CeedQFunctionAddInput(qf_apply, "q_data", bp_options[bp_choice].q_data_size, 743 CEED_EVAL_NONE); 744 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u*out_scale, 745 bp_options[bp_choice].out_mode); 746 747 // Create the error qfunction 748 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].error, 749 bp_options[bp_choice].error_loc, &qf_error); 750 CeedQFunctionAddInput(qf_error, "u", num_comp_u, CEED_EVAL_INTERP); 751 CeedQFunctionAddInput(qf_error, "true_soln", num_comp_u, CEED_EVAL_NONE); 752 CeedQFunctionAddInput(qf_error, "qdata", bp_options[bp_choice].q_data_size, 753 CEED_EVAL_NONE); 754 CeedQFunctionAddOutput(qf_error, "error", num_comp_u, CEED_EVAL_NONE); 755 756 // Create the persistent vectors that will be needed in setup 757 CeedInt num_qpts; 758 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 759 CeedVectorCreate(ceed, bp_options[bp_choice].q_data_size*num_elem*num_qpts, 760 &q_data); 761 CeedVectorCreate(ceed, num_elem*num_qpts*num_comp_u, &target); 762 CeedVectorCreate(ceed, l_size*num_comp_u, &rhs_ceed); 763 764 // Create the operator that builds the quadrature data for the ceed operator 765 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, 766 CEED_QFUNCTION_NONE, &op_setup_geo); 767 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, 768 CEED_VECTOR_ACTIVE); 769 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, 770 CEED_VECTOR_ACTIVE); 771 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, 772 CEED_VECTOR_NONE); 773 CeedOperatorSetField(op_setup_geo, "q_data", elem_restr_qd_i, 774 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 775 776 // Create the operator that builds the RHS and true solution 777 CeedOperatorCreate(ceed, qf_setup_rhs, CEED_QFUNCTION_NONE, 778 CEED_QFUNCTION_NONE, &op_setup_rhs); 779 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, 780 CEED_VECTOR_ACTIVE); 781 CeedOperatorSetField(op_setup_rhs, "q_data", elem_restr_qd_i, 782 CEED_BASIS_COLLOCATED, 783 q_data); 784 CeedOperatorSetField(op_setup_rhs, "true_soln", elem_restr_u_i, 785 CEED_BASIS_COLLOCATED, target); 786 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, 787 CEED_VECTOR_ACTIVE); 788 789 // Create the mass or diff operator 790 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 791 &op_apply); 792 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 793 CeedOperatorSetField(op_apply, "q_data", elem_restr_qd_i, CEED_BASIS_COLLOCATED, 794 q_data); 795 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 796 797 // Create the error operator 798 CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 799 &op_error); 800 CeedOperatorSetField(op_error, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 801 CeedOperatorSetField(op_error, "true_soln", elem_restr_u_i, 802 CEED_BASIS_COLLOCATED, target); 803 CeedOperatorSetField(op_error, "qdata", elem_restr_qd_i, 804 CEED_BASIS_COLLOCATED, q_data); 805 CeedOperatorSetField(op_error, "error", elem_restr_u_i, CEED_BASIS_COLLOCATED, 806 CEED_VECTOR_ACTIVE); 807 808 // Set up Mat 809 ierr = PetscMalloc1(1, &op_apply_ctx); CHKERRQ(ierr); 810 op_apply_ctx->comm = comm; 811 op_apply_ctx->l_to_g = l_to_g; 812 if (bp_choice != CEED_BP1 && bp_choice != CEED_BP2) { 813 op_apply_ctx->l_to_g_0 = l_to_g_0; 814 op_apply_ctx->g_to_g_D = g_to_g_D; 815 } 816 op_apply_ctx->X_loc = X_loc; 817 ierr = VecDuplicate(X_loc, &op_apply_ctx->Y_loc); CHKERRQ(ierr); 818 CeedVectorCreate(ceed, l_size*num_comp_u, &op_apply_ctx->x_ceed); 819 CeedVectorCreate(ceed, l_size*num_comp_u, &op_apply_ctx->y_ceed); 820 op_apply_ctx->op = op_apply; 821 op_apply_ctx->q_data = q_data; 822 op_apply_ctx->ceed = ceed; 823 824 ierr = MatCreateShell(comm, m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 825 m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 826 PETSC_DECIDE, PETSC_DECIDE, op_apply_ctx, &mat); CHKERRQ(ierr); 827 if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) { 828 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Mass); 829 CHKERRQ(ierr); 830 } else { 831 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Diff); 832 CHKERRQ(ierr); 833 } 834 ierr = VecGetType(X, &vec_type); CHKERRQ(ierr); 835 ierr = MatShellSetVecType(mat, vec_type); CHKERRQ(ierr); 836 837 // Get RHS vector 838 ierr = VecDuplicate(X, &rhs); CHKERRQ(ierr); 839 ierr = VecDuplicate(X_loc, &rhs_loc); CHKERRQ(ierr); 840 ierr = VecZeroEntries(rhs_loc); CHKERRQ(ierr); 841 ierr = VecGetArrayAndMemType(rhs_loc, &r, &mem_type); CHKERRQ(ierr); 842 CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r); 843 844 // Setup q_data, rhs, and target 845 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 846 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 847 CeedVectorDestroy(&x_coord); 848 849 // Gather RHS 850 ierr = CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL); CHKERRQ(ierr); 851 ierr = VecRestoreArrayAndMemType(rhs_loc, &r); CHKERRQ(ierr); 852 ierr = VecZeroEntries(rhs); CHKERRQ(ierr); 853 ierr = VecScatterBegin(l_to_g, rhs_loc, rhs, ADD_VALUES, SCATTER_FORWARD); 854 CHKERRQ(ierr); 855 ierr = VecScatterEnd(l_to_g, rhs_loc, rhs, ADD_VALUES, SCATTER_FORWARD); 856 CHKERRQ(ierr); 857 CeedVectorDestroy(&rhs_ceed); 858 859 ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr); 860 { 861 PC pc; 862 ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); 863 if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) { 864 ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr); 865 ierr = PCJacobiSetType(pc, PC_JACOBI_ROWSUM); CHKERRQ(ierr); 866 } else { 867 ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr); 868 } 869 ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr); 870 ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr); 871 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 872 PETSC_DEFAULT); CHKERRQ(ierr); 873 } 874 ierr = KSPSetOperators(ksp, mat, mat); CHKERRQ(ierr); 875 // First run's performance log is not considered for benchmarking purposes 876 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1); 877 CHKERRQ(ierr); 878 my_rt_start = MPI_Wtime(); 879 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 880 my_rt = MPI_Wtime() - my_rt_start; 881 ierr = MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, comm); 882 CHKERRQ(ierr); 883 // Set maxits based on first iteration timing 884 if (my_rt > 0.02) { 885 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 886 ksp_max_it_clip[0]); CHKERRQ(ierr); 887 } else { 888 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 889 ksp_max_it_clip[1]); CHKERRQ(ierr); 890 } 891 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); 892 893 // Timed solve 894 ierr = VecZeroEntries(X); CHKERRQ(ierr); 895 ierr = PetscBarrier((PetscObject)ksp); CHKERRQ(ierr); 896 897 // -- Performance logging 898 ierr = PetscLogStageRegister("Solve Stage", &solve_stage); CHKERRQ(ierr); 899 ierr = PetscLogStagePush(solve_stage); CHKERRQ(ierr); 900 901 // -- Solve 902 my_rt_start = MPI_Wtime(); 903 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 904 my_rt = MPI_Wtime() - my_rt_start; 905 906 // -- Performance logging 907 ierr = PetscLogStagePop(); 908 909 // Output results 910 { 911 KSPType ksp_type; 912 KSPConvergedReason reason; 913 PetscReal rnorm; 914 PetscInt its; 915 ierr = KSPGetType(ksp, &ksp_type); CHKERRQ(ierr); 916 ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr); 917 ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr); 918 ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr); 919 if (!test_mode || reason < 0 || rnorm > 1e-8) { 920 ierr = PetscPrintf(comm, 921 " KSP:\n" 922 " KSP Type : %s\n" 923 " KSP Convergence : %s\n" 924 " Total KSP Iterations : %" PetscInt_FMT "\n" 925 " Final rnorm : %e\n", 926 ksp_type, KSPConvergedReasons[reason], its, 927 (double)rnorm); CHKERRQ(ierr); 928 } 929 if (!test_mode) { 930 ierr = PetscPrintf(comm," Performance:\n"); CHKERRQ(ierr); 931 } 932 { 933 PetscReal max_error; 934 ierr = ComputeErrorMax(op_apply_ctx, op_error, X, target, &max_error); 935 CHKERRQ(ierr); 936 PetscReal tol = 5e-2; 937 if (!test_mode || max_error > tol) { 938 ierr = MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, comm); 939 CHKERRQ(ierr); 940 ierr = MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, comm); 941 CHKERRQ(ierr); 942 ierr = PetscPrintf(comm, 943 " Pointwise Error (max) : %e\n" 944 " CG Solve Time : %g (%g) sec\n", 945 (double)max_error, rt_max, rt_min); CHKERRQ(ierr); 946 } 947 } 948 if (!test_mode) { 949 ierr = PetscPrintf(comm, 950 " DoFs/Sec in CG : %g (%g) million\n", 951 1e-6*gsize*its/rt_max, 952 1e-6*gsize*its/rt_min); CHKERRQ(ierr); 953 } 954 } 955 956 if (write_solution) { 957 PetscViewer vtk_viewer_soln; 958 959 ierr = PetscViewerCreate(comm, &vtk_viewer_soln); CHKERRQ(ierr); 960 ierr = PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK); CHKERRQ(ierr); 961 ierr = PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu"); CHKERRQ(ierr); 962 ierr = VecView(X, vtk_viewer_soln); CHKERRQ(ierr); 963 ierr = PetscViewerDestroy(&vtk_viewer_soln); CHKERRQ(ierr); 964 } 965 966 ierr = VecDestroy(&rhs); CHKERRQ(ierr); 967 ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr); 968 ierr = VecDestroy(&X); CHKERRQ(ierr); 969 ierr = VecDestroy(&op_apply_ctx->X_loc); CHKERRQ(ierr); 970 ierr = VecDestroy(&op_apply_ctx->Y_loc); CHKERRQ(ierr); 971 ierr = VecScatterDestroy(&l_to_g); CHKERRQ(ierr); 972 ierr = VecScatterDestroy(&l_to_g_0); CHKERRQ(ierr); 973 ierr = VecScatterDestroy(&g_to_g_D); CHKERRQ(ierr); 974 ierr = MatDestroy(&mat); CHKERRQ(ierr); 975 ierr = KSPDestroy(&ksp); CHKERRQ(ierr); 976 977 CeedVectorDestroy(&op_apply_ctx->x_ceed); 978 CeedVectorDestroy(&op_apply_ctx->y_ceed); 979 CeedVectorDestroy(&op_apply_ctx->q_data); 980 CeedVectorDestroy(&target); 981 CeedOperatorDestroy(&op_setup_geo); 982 CeedOperatorDestroy(&op_setup_rhs); 983 CeedOperatorDestroy(&op_apply); 984 CeedOperatorDestroy(&op_error); 985 CeedElemRestrictionDestroy(&elem_restr_u); 986 CeedElemRestrictionDestroy(&elem_restr_x); 987 CeedElemRestrictionDestroy(&elem_restr_u_i); 988 CeedElemRestrictionDestroy(&elem_restr_qd_i); 989 CeedQFunctionDestroy(&qf_setup_geo); 990 CeedQFunctionDestroy(&qf_setup_rhs); 991 CeedQFunctionDestroy(&qf_apply); 992 CeedQFunctionDestroy(&qf_error); 993 CeedBasisDestroy(&basis_u); 994 CeedBasisDestroy(&basis_x); 995 CeedDestroy(&ceed); 996 ierr = PetscFree(op_apply_ctx); CHKERRQ(ierr); 997 return PetscFinalize(); 998 } 999