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