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 User_ *User; 133 struct User_ { 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 User user; 266 PetscScalar *x, *y; 267 PetscMemType x_mem_type, y_mem_type; 268 269 PetscFunctionBeginUser; 270 271 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 272 273 // Global-to-local 274 ierr = VecScatterBegin(user->l_to_g, X, user->X_loc, INSERT_VALUES, 275 SCATTER_REVERSE); CHKERRQ(ierr); 276 ierr = VecScatterEnd(user->l_to_g, X, user->X_loc, INSERT_VALUES, 277 SCATTER_REVERSE); CHKERRQ(ierr); 278 279 // Setup libCEED vectors 280 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 281 &x_mem_type); CHKERRQ(ierr); 282 ierr = VecGetArrayAndMemType(user->Y_loc, &y, &y_mem_type); CHKERRQ(ierr); 283 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 284 CeedVectorSetArray(user->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 285 286 // Apply libCEED operator 287 CeedOperatorApply(user->op, user->x_ceed, user->y_ceed, 288 CEED_REQUEST_IMMEDIATE); 289 290 // Restore PETSc vectors 291 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 292 CeedVectorTakeArray(user->y_ceed, MemTypeP2C(y_mem_type), NULL); 293 ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); 294 CHKERRQ(ierr); 295 ierr = VecRestoreArrayAndMemType(user->Y_loc, &y); CHKERRQ(ierr); 296 297 // Local-to-global 298 if (Y) { 299 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 300 ierr = VecScatterBegin(user->l_to_g, user->Y_loc, Y, ADD_VALUES, 301 SCATTER_FORWARD); CHKERRQ(ierr); 302 ierr = VecScatterEnd(user->l_to_g, user->Y_loc, Y, ADD_VALUES, 303 SCATTER_FORWARD); CHKERRQ(ierr); 304 } 305 PetscFunctionReturn(0); 306 } 307 308 // This function uses libCEED to compute the action of the Laplacian with 309 // Dirichlet boundary conditions 310 static PetscErrorCode MatMult_Diff(Mat A, Vec X, Vec Y) { 311 PetscErrorCode ierr; 312 User user; 313 PetscScalar *x, *y; 314 PetscMemType x_mem_type, y_mem_type; 315 316 PetscFunctionBeginUser; 317 318 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 319 320 // Global-to-local 321 ierr = VecScatterBegin(user->l_to_g_0, X, user->X_loc, INSERT_VALUES, 322 SCATTER_REVERSE); CHKERRQ(ierr); 323 ierr = VecScatterEnd(user->l_to_g_0, X, user->X_loc, INSERT_VALUES, 324 SCATTER_REVERSE); 325 CHKERRQ(ierr); 326 327 // Setup libCEED vectors 328 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 329 &x_mem_type); CHKERRQ(ierr); 330 ierr = VecGetArrayAndMemType(user->Y_loc, &y, &y_mem_type); CHKERRQ(ierr); 331 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 332 CeedVectorSetArray(user->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 333 334 // Apply libCEED operator 335 CeedOperatorApply(user->op, user->x_ceed, user->y_ceed, 336 CEED_REQUEST_IMMEDIATE); 337 338 // Restore PETSc vectors 339 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 340 CeedVectorTakeArray(user->y_ceed, MemTypeP2C(y_mem_type), NULL); 341 ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); 342 CHKERRQ(ierr); 343 ierr = VecRestoreArrayAndMemType(user->Y_loc, &y); CHKERRQ(ierr); 344 345 // Local-to-global 346 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 347 ierr = VecScatterBegin(user->g_to_g_D, X, Y, INSERT_VALUES, SCATTER_FORWARD); 348 CHKERRQ(ierr); 349 ierr = VecScatterEnd(user->g_to_g_D, X, Y, INSERT_VALUES, SCATTER_FORWARD); 350 CHKERRQ(ierr); 351 ierr = VecScatterBegin(user->l_to_g_0, user->Y_loc, Y, ADD_VALUES, 352 SCATTER_FORWARD); CHKERRQ(ierr); 353 ierr = VecScatterEnd(user->l_to_g_0, user->Y_loc, Y, ADD_VALUES, 354 SCATTER_FORWARD); 355 CHKERRQ(ierr); 356 357 PetscFunctionReturn(0); 358 } 359 360 // This function calculates the error in the final solution 361 static PetscErrorCode ComputeErrorMax(User user, CeedOperator op_error, Vec X, 362 CeedVector target, PetscReal *max_error) { 363 PetscErrorCode ierr; 364 PetscScalar *x; 365 PetscMemType mem_type; 366 CeedVector collocated_error; 367 CeedSize length; 368 369 PetscFunctionBeginUser; 370 371 CeedVectorGetLength(target, &length); 372 CeedVectorCreate(user->ceed, length, &collocated_error); 373 374 // Global-to-local 375 ierr = VecScatterBegin(user->l_to_g, X, user->X_loc, INSERT_VALUES, 376 SCATTER_REVERSE); CHKERRQ(ierr); 377 ierr = VecScatterEnd(user->l_to_g, X, user->X_loc, INSERT_VALUES, 378 SCATTER_REVERSE); CHKERRQ(ierr); 379 380 // Setup libCEED vector 381 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 382 &mem_type); CHKERRQ(ierr); 383 CeedVectorSetArray(user->x_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, x); 384 385 // Apply libCEED operator 386 CeedOperatorApply(op_error, user->x_ceed, collocated_error, 387 CEED_REQUEST_IMMEDIATE); 388 389 // Restore PETSc vector 390 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(mem_type), NULL); 391 ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); 392 CHKERRQ(ierr); 393 394 // Reduce max error 395 *max_error = 0; 396 const CeedScalar *e; 397 CeedVectorGetArrayRead(collocated_error, CEED_MEM_HOST, &e); 398 for (CeedInt i=0; i<length; i++) { 399 *max_error = PetscMax(*max_error, PetscAbsScalar(e[i])); 400 } 401 CeedVectorRestoreArrayRead(collocated_error, &e); 402 ierr = MPI_Allreduce(MPI_IN_PLACE, max_error, 1, MPIU_REAL, MPIU_MAX, 403 user->comm); CHKERRQ(ierr); 404 405 // Cleanup 406 CeedVectorDestroy(&collocated_error); 407 408 PetscFunctionReturn(0); 409 } 410 411 int main(int argc, char **argv) { 412 PetscInt ierr; 413 MPI_Comm comm; 414 char ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 415 double my_rt_start, my_rt, rt_min, rt_max; 416 PetscInt degree, q_extra, local_nodes, local_elem, mesh_elem[3], m_nodes[3], 417 p[3], 418 i_rank[3], l_nodes[3], l_size, num_comp_u = 1, ksp_max_it_clip[2]; 419 PetscScalar *r; 420 PetscBool test_mode, benchmark_mode, write_solution; 421 PetscMPIInt size, rank; 422 PetscLogStage solve_stage; 423 Vec X, X_loc, rhs, rhs_loc; 424 Mat mat; 425 KSP ksp; 426 VecScatter l_to_g, l_to_g_0, g_to_g_D; 427 PetscMemType mem_type; 428 User user; 429 Ceed ceed; 430 CeedBasis basis_x, basis_u; 431 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i; 432 CeedQFunction qf_setup_geo, qf_setup_rhs, qf_apply, qf_error; 433 CeedOperator op_setup_geo, op_setup_rhs, op_apply, op_error; 434 CeedVector x_coord, q_data, rhs_ceed, target; 435 CeedInt P, Q; 436 const CeedInt dim = 3, num_comp_x = 3; 437 BPType bp_choice; 438 439 ierr = PetscInitialize(&argc, &argv, NULL, help); 440 if (ierr) return ierr; 441 comm = PETSC_COMM_WORLD; 442 443 // Read command line options 444 PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); 445 bp_choice = CEED_BP1; 446 ierr = PetscOptionsEnum("-problem", 447 "CEED benchmark problem to solve", NULL, 448 bp_types, (PetscEnum)bp_choice, (PetscEnum *)&bp_choice, 449 NULL); CHKERRQ(ierr); 450 num_comp_u = bp_options[bp_choice].num_comp_u; 451 test_mode = PETSC_FALSE; 452 ierr = PetscOptionsBool("-test", 453 "Testing mode (do not print unless error is large)", 454 NULL, test_mode, &test_mode, NULL); CHKERRQ(ierr); 455 benchmark_mode = PETSC_FALSE; 456 ierr = PetscOptionsBool("-benchmark", 457 "Benchmarking mode (prints benchmark statistics)", 458 NULL, benchmark_mode, &benchmark_mode, NULL); 459 CHKERRQ(ierr); 460 write_solution = PETSC_FALSE; 461 ierr = PetscOptionsBool("-write_solution", 462 "Write solution for visualization", 463 NULL, write_solution, &write_solution, NULL); 464 CHKERRQ(ierr); 465 degree = test_mode ? 3 : 1; 466 ierr = PetscOptionsInt("-degree", "Polynomial degree of tensor product basis", 467 NULL, degree, °ree, NULL); CHKERRQ(ierr); 468 q_extra = bp_options[bp_choice].q_extra; 469 ierr = PetscOptionsInt("-q_extra", "Number of extra quadrature points", 470 NULL, q_extra, &q_extra, NULL); CHKERRQ(ierr); 471 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 472 NULL, ceed_resource, ceed_resource, 473 sizeof(ceed_resource), NULL); CHKERRQ(ierr); 474 local_nodes = 1000; 475 ierr = PetscOptionsInt("-local", 476 "Target number of locally owned nodes per process", 477 NULL, local_nodes, &local_nodes, NULL); CHKERRQ(ierr); 478 PetscInt two = 2; 479 ksp_max_it_clip[0] = 5; 480 ksp_max_it_clip[1] = 20; 481 ierr = PetscOptionsIntArray("-ksp_max_it_clip", 482 "Min and max number of iterations to use during benchmarking", 483 NULL, ksp_max_it_clip, &two, NULL); 484 CHKERRQ(ierr); 485 PetscOptionsEnd(); 486 P = degree + 1; 487 Q = P + q_extra; 488 489 // Set up libCEED 490 CeedInit(ceed_resource, &ceed); 491 CeedMemType mem_type_backend; 492 CeedGetPreferredMemType(ceed, &mem_type_backend); 493 494 VecType default_vec_type = NULL, vec_type; 495 switch (mem_type_backend) { 496 case CEED_MEM_HOST: default_vec_type = VECSTANDARD; break; 497 case CEED_MEM_DEVICE: { 498 const char *resolved; 499 CeedGetResource(ceed, &resolved); 500 if (strstr(resolved, "/gpu/cuda")) default_vec_type = VECCUDA; 501 else if (strstr(resolved, "/gpu/hip/occa")) 502 default_vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 503 else if (strstr(resolved, "/gpu/hip")) default_vec_type = VECHIP; 504 else default_vec_type = VECSTANDARD; 505 } 506 } 507 508 // Determine size of process grid 509 ierr = MPI_Comm_size(comm, &size); CHKERRQ(ierr); 510 Split3(size, p, false); 511 512 // Find a nicely composite number of elements no less than local_nodes 513 for (local_elem = PetscMax(1, local_nodes / (degree*degree*degree)); ; 514 local_elem++) { 515 Split3(local_elem, mesh_elem, true); 516 if (Max3(mesh_elem) / Min3(mesh_elem) <= 2) break; 517 } 518 519 // Find my location in the process grid 520 ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr); 521 for (int d=0, rank_left=rank; d<dim; d++) { 522 const int pstride[3] = {p[1] *p[2], p[2], 1}; 523 i_rank[d] = rank_left / pstride[d]; 524 rank_left -= i_rank[d] * pstride[d]; 525 } 526 527 GlobalNodes(p, i_rank, degree, mesh_elem, m_nodes); 528 529 // Setup global vector 530 ierr = VecCreate(comm, &X); CHKERRQ(ierr); 531 ierr = VecSetType(X, default_vec_type); CHKERRQ(ierr); 532 ierr = VecSetSizes(X, m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 533 PETSC_DECIDE); 534 CHKERRQ(ierr); 535 ierr = VecSetFromOptions(X); CHKERRQ(ierr); 536 ierr = VecSetUp(X); CHKERRQ(ierr); 537 538 // Set up libCEED 539 CeedInit(ceed_resource, &ceed); 540 541 // Print summary 542 CeedInt gsize; 543 ierr = VecGetSize(X, &gsize); CHKERRQ(ierr); 544 if (!test_mode) { 545 const char *used_resource; 546 CeedGetResource(ceed, &used_resource); 547 548 ierr = VecGetType(X, &vec_type); CHKERRQ(ierr); 549 550 ierr = PetscPrintf(comm, 551 "\n-- CEED Benchmark Problem %" CeedInt_FMT " -- libCEED + PETSc --\n" 552 " PETSc:\n" 553 " PETSc Vec Type : %s\n" 554 " libCEED:\n" 555 " libCEED Backend : %s\n" 556 " libCEED Backend MemType : %s\n" 557 " Mesh:\n" 558 " Number of 1D Basis Nodes (P) : %" CeedInt_FMT "\n" 559 " Number of 1D Quadrature Points (Q) : %" CeedInt_FMT "\n" 560 " Global nodes : %" PetscInt_FMT "\n" 561 " Process Decomposition : %" PetscInt_FMT 562 " %" PetscInt_FMT " %" PetscInt_FMT "\n" 563 " Local Elements : %" PetscInt_FMT 564 " = %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n" 565 " Owned nodes : %" PetscInt_FMT 566 " = %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n" 567 " DoF per node : %" PetscInt_FMT "\n", 568 bp_choice+1, vec_type, used_resource, 569 CeedMemTypes[mem_type_backend], 570 P, Q, gsize/num_comp_u, p[0], p[1], p[2], local_elem, 571 mesh_elem[0], mesh_elem[1], mesh_elem[2], 572 m_nodes[0]*m_nodes[1]*m_nodes[2], m_nodes[0], m_nodes[1], 573 m_nodes[2], num_comp_u); CHKERRQ(ierr); 574 } 575 576 { 577 l_size = 1; 578 for (int d=0; d<dim; d++) { 579 l_nodes[d] = mesh_elem[d]*degree + 1; 580 l_size *= l_nodes[d]; 581 } 582 ierr = VecCreate(PETSC_COMM_SELF, &X_loc); CHKERRQ(ierr); 583 ierr = VecSetType(X_loc, default_vec_type); CHKERRQ(ierr); 584 ierr = VecSetSizes(X_loc, l_size*num_comp_u, PETSC_DECIDE); CHKERRQ(ierr); 585 ierr = VecSetFromOptions(X_loc); CHKERRQ(ierr); 586 ierr = VecSetUp(X_loc); CHKERRQ(ierr); 587 588 // Create local-to-global scatter 589 PetscInt *l_to_g_ind, *l_to_g_ind_0, *loc_ind, l_0_count; 590 IS l_to_g_is, l_to_g_is_0, loc_is; 591 PetscInt g_start[2][2][2], g_m_nodes[2][2][2][dim]; 592 593 for (int i=0; i<2; i++) { 594 for (int j=0; j<2; j++) { 595 for (int k=0; k<2; k++) { 596 PetscInt ijk_rank[3] = {i_rank[0]+i, i_rank[1]+j, i_rank[2]+k}; 597 g_start[i][j][k] = GlobalStart(p, ijk_rank, degree, mesh_elem); 598 GlobalNodes(p, ijk_rank, degree, mesh_elem, g_m_nodes[i][j][k]); 599 } 600 } 601 } 602 603 ierr = PetscMalloc1(l_size, &l_to_g_ind); CHKERRQ(ierr); 604 ierr = PetscMalloc1(l_size, &l_to_g_ind_0); CHKERRQ(ierr); 605 ierr = PetscMalloc1(l_size, &loc_ind); CHKERRQ(ierr); 606 l_0_count = 0; 607 for (PetscInt i=0,ir,ii; ir=i>=m_nodes[0], ii=i-ir*m_nodes[0], i<l_nodes[0]; 608 i++) 609 for (PetscInt j=0,jr,jj; jr=j>=m_nodes[1], jj=j-jr*m_nodes[1], j<l_nodes[1]; 610 j++) 611 for (PetscInt k=0,kr,kk; kr=k>=m_nodes[2], kk=k-kr*m_nodes[2], k<l_nodes[2]; 612 k++) { 613 PetscInt here = (i*l_nodes[1]+j)*l_nodes[2]+k; 614 l_to_g_ind[here] = 615 g_start[ir][jr][kr] + (ii*g_m_nodes[ir][jr][kr][1]+jj)*g_m_nodes[ir][jr][kr][2] 616 +kk; 617 if ((i_rank[0] == 0 && i == 0) 618 || (i_rank[1] == 0 && j == 0) 619 || (i_rank[2] == 0 && k == 0) 620 || (i_rank[0]+1 == p[0] && i+1 == l_nodes[0]) 621 || (i_rank[1]+1 == p[1] && j+1 == l_nodes[1]) 622 || (i_rank[2]+1 == p[2] && k+1 == l_nodes[2])) 623 continue; 624 l_to_g_ind_0[l_0_count] = l_to_g_ind[here]; 625 loc_ind[l_0_count++] = here; 626 } 627 ierr = ISCreateBlock(comm, num_comp_u, l_size, l_to_g_ind, PETSC_OWN_POINTER, 628 &l_to_g_is); CHKERRQ(ierr); 629 ierr = VecScatterCreate(X_loc, NULL, X, l_to_g_is, &l_to_g); CHKERRQ(ierr); 630 CHKERRQ(ierr); 631 ierr = ISCreateBlock(comm, num_comp_u, l_0_count, l_to_g_ind_0, 632 PETSC_OWN_POINTER, 633 &l_to_g_is_0); CHKERRQ(ierr); 634 ierr = ISCreateBlock(comm, num_comp_u, l_0_count, loc_ind, PETSC_OWN_POINTER, 635 &loc_is); CHKERRQ(ierr); 636 ierr = VecScatterCreate(X_loc, loc_is, X, l_to_g_is_0, &l_to_g_0); 637 CHKERRQ(ierr); 638 { 639 // Create global-to-global scatter for Dirichlet values (everything not in 640 // l_to_g_is_0, which is the range of l_to_g_0) 641 PetscInt x_start, x_end, *ind_D, count_D = 0; 642 IS is_D; 643 const PetscScalar *x; 644 ierr = VecZeroEntries(X_loc); CHKERRQ(ierr); 645 ierr = VecSet(X, 1.0); CHKERRQ(ierr); 646 ierr = VecScatterBegin(l_to_g_0, X_loc, X, INSERT_VALUES, SCATTER_FORWARD); 647 CHKERRQ(ierr); 648 ierr = VecScatterEnd(l_to_g_0, X_loc, X, INSERT_VALUES, SCATTER_FORWARD); 649 CHKERRQ(ierr); 650 ierr = VecGetOwnershipRange(X, &x_start, &x_end); CHKERRQ(ierr); 651 ierr = PetscMalloc1(x_end-x_start, &ind_D); CHKERRQ(ierr); 652 ierr = VecGetArrayRead(X, &x); CHKERRQ(ierr); 653 for (PetscInt i=0; i<x_end-x_start; i++) { 654 if (x[i] == 1.) ind_D[count_D++] = x_start + i; 655 } 656 ierr = VecRestoreArrayRead(X, &x); CHKERRQ(ierr); 657 ierr = ISCreateGeneral(comm, count_D, ind_D, PETSC_COPY_VALUES, &is_D); 658 CHKERRQ(ierr); 659 ierr = PetscFree(ind_D); CHKERRQ(ierr); 660 ierr = VecScatterCreate(X, is_D, X, is_D, &g_to_g_D); CHKERRQ(ierr); 661 ierr = ISDestroy(&is_D); CHKERRQ(ierr); 662 } 663 ierr = ISDestroy(&l_to_g_is); CHKERRQ(ierr); 664 ierr = ISDestroy(&l_to_g_is_0); CHKERRQ(ierr); 665 ierr = ISDestroy(&loc_is); CHKERRQ(ierr); 666 } 667 668 // CEED bases 669 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, 670 bp_options[bp_choice].q_mode, &basis_u); 671 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, 672 bp_options[bp_choice].q_mode, &basis_x); 673 674 // CEED restrictions 675 CreateRestriction(ceed, mesh_elem, P, num_comp_u, &elem_restr_u); 676 CreateRestriction(ceed, mesh_elem, 2, dim, &elem_restr_x); 677 CeedInt num_elem = mesh_elem[0]*mesh_elem[1]*mesh_elem[2]; 678 CeedElemRestrictionCreateStrided(ceed, num_elem, Q*Q*Q, num_comp_u, 679 num_comp_u*num_elem*Q*Q*Q, 680 CEED_STRIDES_BACKEND, &elem_restr_u_i); 681 CeedElemRestrictionCreateStrided(ceed, num_elem, Q*Q*Q, 682 bp_options[bp_choice].q_data_size, 683 bp_options[bp_choice].q_data_size*num_elem*Q*Q*Q, 684 CEED_STRIDES_BACKEND, &elem_restr_qd_i); 685 { 686 CeedScalar *x_loc; 687 CeedInt shape[3] = {mesh_elem[0]+1, mesh_elem[1]+1, mesh_elem[2]+1}, len = 688 shape[0]*shape[1]*shape[2]; 689 x_loc = malloc(len*num_comp_x*sizeof x_loc[0]); 690 for (CeedInt i=0; i<shape[0]; i++) { 691 for (CeedInt j=0; j<shape[1]; j++) { 692 for (CeedInt k=0; k<shape[2]; k++) { 693 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 0] = 1.*(i_rank[0]*mesh_elem[0]+i) / 694 (p[0]*mesh_elem[0]); 695 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 1] = 1.*(i_rank[1]*mesh_elem[1]+j) / 696 (p[1]*mesh_elem[1]); 697 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 2] = 1.*(i_rank[2]*mesh_elem[2]+k) / 698 (p[2]*mesh_elem[2]); 699 } 700 } 701 } 702 CeedVectorCreate(ceed, len*num_comp_x, &x_coord); 703 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_OWN_POINTER, x_loc); 704 } 705 706 // Create the QFunction that builds the operator quadrature data 707 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].setup_geo, 708 bp_options[bp_choice].setup_geo_loc, &qf_setup_geo); 709 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 710 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x*dim, CEED_EVAL_GRAD); 711 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 712 CeedQFunctionAddOutput(qf_setup_geo, "q_data", 713 bp_options[bp_choice].q_data_size, 714 CEED_EVAL_NONE); 715 716 // Create the QFunction that sets up the RHS and true solution 717 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].setup_rhs, 718 bp_options[bp_choice].setup_rhs_loc, &qf_setup_rhs); 719 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 720 CeedQFunctionAddInput(qf_setup_rhs, "q_data", bp_options[bp_choice].q_data_size, 721 CEED_EVAL_NONE); 722 CeedQFunctionAddOutput(qf_setup_rhs, "true_soln", num_comp_u, CEED_EVAL_NONE); 723 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 724 725 // Set up PDE operator 726 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].apply, 727 bp_options[bp_choice].apply_loc, &qf_apply); 728 // Add inputs and outputs 729 CeedInt in_scale = bp_options[bp_choice].in_mode==CEED_EVAL_GRAD ? 3 : 1; 730 CeedInt out_scale = bp_options[bp_choice].out_mode==CEED_EVAL_GRAD ? 3 : 1; 731 CeedQFunctionAddInput(qf_apply, "u", num_comp_u*in_scale, 732 bp_options[bp_choice].in_mode); 733 CeedQFunctionAddInput(qf_apply, "q_data", bp_options[bp_choice].q_data_size, 734 CEED_EVAL_NONE); 735 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u*out_scale, 736 bp_options[bp_choice].out_mode); 737 738 // Create the error qfunction 739 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].error, 740 bp_options[bp_choice].error_loc, &qf_error); 741 CeedQFunctionAddInput(qf_error, "u", num_comp_u, CEED_EVAL_INTERP); 742 CeedQFunctionAddInput(qf_error, "true_soln", num_comp_u, CEED_EVAL_NONE); 743 CeedQFunctionAddOutput(qf_error, "error", num_comp_u, CEED_EVAL_NONE); 744 745 // Create the persistent vectors that will be needed in setup 746 CeedInt num_qpts; 747 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 748 CeedVectorCreate(ceed, bp_options[bp_choice].q_data_size*num_elem*num_qpts, 749 &q_data); 750 CeedVectorCreate(ceed, num_elem*num_qpts*num_comp_u, &target); 751 CeedVectorCreate(ceed, l_size*num_comp_u, &rhs_ceed); 752 753 // Create the operator that builds the quadrature data for the ceed operator 754 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, 755 CEED_QFUNCTION_NONE, &op_setup_geo); 756 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, 757 CEED_VECTOR_ACTIVE); 758 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, 759 CEED_VECTOR_ACTIVE); 760 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, 761 CEED_VECTOR_NONE); 762 CeedOperatorSetField(op_setup_geo, "q_data", elem_restr_qd_i, 763 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 764 765 // Create the operator that builds the RHS and true solution 766 CeedOperatorCreate(ceed, qf_setup_rhs, CEED_QFUNCTION_NONE, 767 CEED_QFUNCTION_NONE, &op_setup_rhs); 768 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, 769 CEED_VECTOR_ACTIVE); 770 CeedOperatorSetField(op_setup_rhs, "q_data", elem_restr_qd_i, 771 CEED_BASIS_COLLOCATED, 772 q_data); 773 CeedOperatorSetField(op_setup_rhs, "true_soln", elem_restr_u_i, 774 CEED_BASIS_COLLOCATED, target); 775 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, 776 CEED_VECTOR_ACTIVE); 777 778 // Create the mass or diff operator 779 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 780 &op_apply); 781 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 782 CeedOperatorSetField(op_apply, "q_data", elem_restr_qd_i, CEED_BASIS_COLLOCATED, 783 q_data); 784 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 785 786 // Create the error operator 787 CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 788 &op_error); 789 CeedOperatorSetField(op_error, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 790 CeedOperatorSetField(op_error, "true_soln", elem_restr_u_i, 791 CEED_BASIS_COLLOCATED, target); 792 CeedOperatorSetField(op_error, "error", elem_restr_u_i, CEED_BASIS_COLLOCATED, 793 CEED_VECTOR_ACTIVE); 794 795 // Set up Mat 796 ierr = PetscMalloc1(1, &user); CHKERRQ(ierr); 797 user->comm = comm; 798 user->l_to_g = l_to_g; 799 if (bp_choice != CEED_BP1 && bp_choice != CEED_BP2) { 800 user->l_to_g_0 = l_to_g_0; 801 user->g_to_g_D = g_to_g_D; 802 } 803 user->X_loc = X_loc; 804 ierr = VecDuplicate(X_loc, &user->Y_loc); CHKERRQ(ierr); 805 CeedVectorCreate(ceed, l_size*num_comp_u, &user->x_ceed); 806 CeedVectorCreate(ceed, l_size*num_comp_u, &user->y_ceed); 807 user->op = op_apply; 808 user->q_data = q_data; 809 user->ceed = ceed; 810 811 ierr = MatCreateShell(comm, m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 812 m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 813 PETSC_DECIDE, PETSC_DECIDE, user, &mat); CHKERRQ(ierr); 814 if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) { 815 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Mass); 816 CHKERRQ(ierr); 817 } else { 818 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Diff); 819 CHKERRQ(ierr); 820 } 821 ierr = VecGetType(X, &vec_type); CHKERRQ(ierr); 822 ierr = MatShellSetVecType(mat, vec_type); CHKERRQ(ierr); 823 824 // Get RHS vector 825 ierr = VecDuplicate(X, &rhs); CHKERRQ(ierr); 826 ierr = VecDuplicate(X_loc, &rhs_loc); CHKERRQ(ierr); 827 ierr = VecZeroEntries(rhs_loc); CHKERRQ(ierr); 828 ierr = VecGetArrayAndMemType(rhs_loc, &r, &mem_type); CHKERRQ(ierr); 829 CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r); 830 831 // Setup q_data, rhs, and target 832 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 833 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 834 CeedVectorDestroy(&x_coord); 835 836 // Gather RHS 837 ierr = CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL); CHKERRQ(ierr); 838 ierr = VecRestoreArrayAndMemType(rhs_loc, &r); CHKERRQ(ierr); 839 ierr = VecZeroEntries(rhs); CHKERRQ(ierr); 840 ierr = VecScatterBegin(l_to_g, rhs_loc, rhs, ADD_VALUES, SCATTER_FORWARD); 841 CHKERRQ(ierr); 842 ierr = VecScatterEnd(l_to_g, rhs_loc, rhs, ADD_VALUES, SCATTER_FORWARD); 843 CHKERRQ(ierr); 844 CeedVectorDestroy(&rhs_ceed); 845 846 ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr); 847 { 848 PC pc; 849 ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); 850 if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) { 851 ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr); 852 ierr = PCJacobiSetType(pc, PC_JACOBI_ROWSUM); CHKERRQ(ierr); 853 } else { 854 ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr); 855 } 856 ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr); 857 ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr); 858 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 859 PETSC_DEFAULT); CHKERRQ(ierr); 860 } 861 ierr = KSPSetOperators(ksp, mat, mat); CHKERRQ(ierr); 862 // First run's performance log is not considered for benchmarking purposes 863 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1); 864 CHKERRQ(ierr); 865 my_rt_start = MPI_Wtime(); 866 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 867 my_rt = MPI_Wtime() - my_rt_start; 868 ierr = MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, comm); 869 CHKERRQ(ierr); 870 // Set maxits based on first iteration timing 871 if (my_rt > 0.02) { 872 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 873 ksp_max_it_clip[0]); CHKERRQ(ierr); 874 } else { 875 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 876 ksp_max_it_clip[1]); CHKERRQ(ierr); 877 } 878 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); 879 880 // Timed solve 881 ierr = VecZeroEntries(X); CHKERRQ(ierr); 882 ierr = PetscBarrier((PetscObject)ksp); CHKERRQ(ierr); 883 884 // -- Performance logging 885 ierr = PetscLogStageRegister("Solve Stage", &solve_stage); CHKERRQ(ierr); 886 ierr = PetscLogStagePush(solve_stage); CHKERRQ(ierr); 887 888 // -- Solve 889 my_rt_start = MPI_Wtime(); 890 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 891 my_rt = MPI_Wtime() - my_rt_start; 892 893 // -- Performance logging 894 ierr = PetscLogStagePop(); 895 896 // Output results 897 { 898 KSPType ksp_type; 899 KSPConvergedReason reason; 900 PetscReal rnorm; 901 PetscInt its; 902 ierr = KSPGetType(ksp, &ksp_type); CHKERRQ(ierr); 903 ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr); 904 ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr); 905 ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr); 906 if (!test_mode || reason < 0 || rnorm > 1e-8) { 907 ierr = PetscPrintf(comm, 908 " KSP:\n" 909 " KSP Type : %s\n" 910 " KSP Convergence : %s\n" 911 " Total KSP Iterations : %" PetscInt_FMT "\n" 912 " Final rnorm : %e\n", 913 ksp_type, KSPConvergedReasons[reason], its, 914 (double)rnorm); CHKERRQ(ierr); 915 } 916 if (!test_mode) { 917 ierr = PetscPrintf(comm," Performance:\n"); CHKERRQ(ierr); 918 } 919 { 920 PetscReal max_error; 921 ierr = ComputeErrorMax(user, op_error, X, target, &max_error); 922 CHKERRQ(ierr); 923 PetscReal tol = 5e-2; 924 if (!test_mode || max_error > tol) { 925 ierr = MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, comm); 926 CHKERRQ(ierr); 927 ierr = MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, comm); 928 CHKERRQ(ierr); 929 ierr = PetscPrintf(comm, 930 " Pointwise Error (max) : %e\n" 931 " CG Solve Time : %g (%g) sec\n", 932 (double)max_error, rt_max, rt_min); CHKERRQ(ierr); 933 } 934 } 935 if (!test_mode) { 936 ierr = PetscPrintf(comm, 937 " DoFs/Sec in CG : %g (%g) million\n", 938 1e-6*gsize*its/rt_max, 939 1e-6*gsize*its/rt_min); CHKERRQ(ierr); 940 } 941 } 942 943 if (write_solution) { 944 PetscViewer vtk_viewer_soln; 945 946 ierr = PetscViewerCreate(comm, &vtk_viewer_soln); CHKERRQ(ierr); 947 ierr = PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK); CHKERRQ(ierr); 948 ierr = PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu"); CHKERRQ(ierr); 949 ierr = VecView(X, vtk_viewer_soln); CHKERRQ(ierr); 950 ierr = PetscViewerDestroy(&vtk_viewer_soln); CHKERRQ(ierr); 951 } 952 953 ierr = VecDestroy(&rhs); CHKERRQ(ierr); 954 ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr); 955 ierr = VecDestroy(&X); CHKERRQ(ierr); 956 ierr = VecDestroy(&user->X_loc); CHKERRQ(ierr); 957 ierr = VecDestroy(&user->Y_loc); CHKERRQ(ierr); 958 ierr = VecScatterDestroy(&l_to_g); CHKERRQ(ierr); 959 ierr = VecScatterDestroy(&l_to_g_0); CHKERRQ(ierr); 960 ierr = VecScatterDestroy(&g_to_g_D); CHKERRQ(ierr); 961 ierr = MatDestroy(&mat); CHKERRQ(ierr); 962 ierr = KSPDestroy(&ksp); CHKERRQ(ierr); 963 964 CeedVectorDestroy(&user->x_ceed); 965 CeedVectorDestroy(&user->y_ceed); 966 CeedVectorDestroy(&user->q_data); 967 CeedVectorDestroy(&target); 968 CeedOperatorDestroy(&op_setup_geo); 969 CeedOperatorDestroy(&op_setup_rhs); 970 CeedOperatorDestroy(&op_apply); 971 CeedOperatorDestroy(&op_error); 972 CeedElemRestrictionDestroy(&elem_restr_u); 973 CeedElemRestrictionDestroy(&elem_restr_x); 974 CeedElemRestrictionDestroy(&elem_restr_u_i); 975 CeedElemRestrictionDestroy(&elem_restr_qd_i); 976 CeedQFunctionDestroy(&qf_setup_geo); 977 CeedQFunctionDestroy(&qf_setup_rhs); 978 CeedQFunctionDestroy(&qf_apply); 979 CeedQFunctionDestroy(&qf_error); 980 CeedBasisDestroy(&basis_u); 981 CeedBasisDestroy(&basis_x); 982 CeedDestroy(&ceed); 983 ierr = PetscFree(user); CHKERRQ(ierr); 984 return PetscFinalize(); 985 } 986