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