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 : %D\n" 559 " Process Decomposition : %D %D %D\n" 560 " Local Elements : %D = %D %D %D\n" 561 " Owned nodes : %D = %D %D %D\n" 562 " DoF per node : %D\n", 563 bp_choice+1, vec_type, used_resource, 564 CeedMemTypes[mem_type_backend], 565 P, Q, gsize/num_comp_u, p[0], p[1], p[2], local_elem, 566 mesh_elem[0], mesh_elem[1], mesh_elem[2], 567 m_nodes[0]*m_nodes[1]*m_nodes[2], m_nodes[0], m_nodes[1], 568 m_nodes[2], num_comp_u); CHKERRQ(ierr); 569 } 570 571 { 572 l_size = 1; 573 for (int d=0; d<dim; d++) { 574 l_nodes[d] = mesh_elem[d]*degree + 1; 575 l_size *= l_nodes[d]; 576 } 577 ierr = VecCreate(PETSC_COMM_SELF, &X_loc); CHKERRQ(ierr); 578 ierr = VecSetType(X_loc, default_vec_type); CHKERRQ(ierr); 579 ierr = VecSetSizes(X_loc, l_size*num_comp_u, PETSC_DECIDE); CHKERRQ(ierr); 580 ierr = VecSetFromOptions(X_loc); CHKERRQ(ierr); 581 ierr = VecSetUp(X_loc); CHKERRQ(ierr); 582 583 // Create local-to-global scatter 584 PetscInt *l_to_g_ind, *l_to_g_ind_0, *loc_ind, l_0_count; 585 IS l_to_g_is, l_to_g_is_0, loc_is; 586 PetscInt g_start[2][2][2], g_m_nodes[2][2][2][dim]; 587 588 for (int i=0; i<2; i++) { 589 for (int j=0; j<2; j++) { 590 for (int k=0; k<2; k++) { 591 PetscInt ijk_rank[3] = {i_rank[0]+i, i_rank[1]+j, i_rank[2]+k}; 592 g_start[i][j][k] = GlobalStart(p, ijk_rank, degree, mesh_elem); 593 GlobalNodes(p, ijk_rank, degree, mesh_elem, g_m_nodes[i][j][k]); 594 } 595 } 596 } 597 598 ierr = PetscMalloc1(l_size, &l_to_g_ind); CHKERRQ(ierr); 599 ierr = PetscMalloc1(l_size, &l_to_g_ind_0); CHKERRQ(ierr); 600 ierr = PetscMalloc1(l_size, &loc_ind); CHKERRQ(ierr); 601 l_0_count = 0; 602 for (PetscInt i=0,ir,ii; ir=i>=m_nodes[0], ii=i-ir*m_nodes[0], i<l_nodes[0]; 603 i++) 604 for (PetscInt j=0,jr,jj; jr=j>=m_nodes[1], jj=j-jr*m_nodes[1], j<l_nodes[1]; 605 j++) 606 for (PetscInt k=0,kr,kk; kr=k>=m_nodes[2], kk=k-kr*m_nodes[2], k<l_nodes[2]; 607 k++) { 608 PetscInt here = (i*l_nodes[1]+j)*l_nodes[2]+k; 609 l_to_g_ind[here] = 610 g_start[ir][jr][kr] + (ii*g_m_nodes[ir][jr][kr][1]+jj)*g_m_nodes[ir][jr][kr][2] 611 +kk; 612 if ((i_rank[0] == 0 && i == 0) 613 || (i_rank[1] == 0 && j == 0) 614 || (i_rank[2] == 0 && k == 0) 615 || (i_rank[0]+1 == p[0] && i+1 == l_nodes[0]) 616 || (i_rank[1]+1 == p[1] && j+1 == l_nodes[1]) 617 || (i_rank[2]+1 == p[2] && k+1 == l_nodes[2])) 618 continue; 619 l_to_g_ind_0[l_0_count] = l_to_g_ind[here]; 620 loc_ind[l_0_count++] = here; 621 } 622 ierr = ISCreateBlock(comm, num_comp_u, l_size, l_to_g_ind, PETSC_OWN_POINTER, 623 &l_to_g_is); CHKERRQ(ierr); 624 ierr = VecScatterCreate(X_loc, NULL, X, l_to_g_is, &l_to_g); CHKERRQ(ierr); 625 CHKERRQ(ierr); 626 ierr = ISCreateBlock(comm, num_comp_u, l_0_count, l_to_g_ind_0, 627 PETSC_OWN_POINTER, 628 &l_to_g_is_0); CHKERRQ(ierr); 629 ierr = ISCreateBlock(comm, num_comp_u, l_0_count, loc_ind, PETSC_OWN_POINTER, 630 &loc_is); CHKERRQ(ierr); 631 ierr = VecScatterCreate(X_loc, loc_is, X, l_to_g_is_0, &l_to_g_0); 632 CHKERRQ(ierr); 633 { 634 // Create global-to-global scatter for Dirichlet values (everything not in 635 // l_to_g_is_0, which is the range of l_to_g_0) 636 PetscInt x_start, x_end, *ind_D, count_D = 0; 637 IS is_D; 638 const PetscScalar *x; 639 ierr = VecZeroEntries(X_loc); CHKERRQ(ierr); 640 ierr = VecSet(X, 1.0); CHKERRQ(ierr); 641 ierr = VecScatterBegin(l_to_g_0, X_loc, X, INSERT_VALUES, SCATTER_FORWARD); 642 CHKERRQ(ierr); 643 ierr = VecScatterEnd(l_to_g_0, X_loc, X, INSERT_VALUES, SCATTER_FORWARD); 644 CHKERRQ(ierr); 645 ierr = VecGetOwnershipRange(X, &x_start, &x_end); CHKERRQ(ierr); 646 ierr = PetscMalloc1(x_end-x_start, &ind_D); CHKERRQ(ierr); 647 ierr = VecGetArrayRead(X, &x); CHKERRQ(ierr); 648 for (PetscInt i=0; i<x_end-x_start; i++) { 649 if (x[i] == 1.) ind_D[count_D++] = x_start + i; 650 } 651 ierr = VecRestoreArrayRead(X, &x); CHKERRQ(ierr); 652 ierr = ISCreateGeneral(comm, count_D, ind_D, PETSC_COPY_VALUES, &is_D); 653 CHKERRQ(ierr); 654 ierr = PetscFree(ind_D); CHKERRQ(ierr); 655 ierr = VecScatterCreate(X, is_D, X, is_D, &g_to_g_D); CHKERRQ(ierr); 656 ierr = ISDestroy(&is_D); CHKERRQ(ierr); 657 } 658 ierr = ISDestroy(&l_to_g_is); CHKERRQ(ierr); 659 ierr = ISDestroy(&l_to_g_is_0); CHKERRQ(ierr); 660 ierr = ISDestroy(&loc_is); CHKERRQ(ierr); 661 } 662 663 // CEED bases 664 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, 665 bp_options[bp_choice].q_mode, &basis_u); 666 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, 667 bp_options[bp_choice].q_mode, &basis_x); 668 669 // CEED restrictions 670 CreateRestriction(ceed, mesh_elem, P, num_comp_u, &elem_restr_u); 671 CreateRestriction(ceed, mesh_elem, 2, dim, &elem_restr_x); 672 CeedInt num_elem = mesh_elem[0]*mesh_elem[1]*mesh_elem[2]; 673 CeedElemRestrictionCreateStrided(ceed, num_elem, Q*Q*Q, num_comp_u, 674 num_comp_u*num_elem*Q*Q*Q, 675 CEED_STRIDES_BACKEND, &elem_restr_u_i); 676 CeedElemRestrictionCreateStrided(ceed, num_elem, Q*Q*Q, 677 bp_options[bp_choice].q_data_size, 678 bp_options[bp_choice].q_data_size*num_elem*Q*Q*Q, 679 CEED_STRIDES_BACKEND, &elem_restr_qd_i); 680 { 681 CeedScalar *x_loc; 682 CeedInt shape[3] = {mesh_elem[0]+1, mesh_elem[1]+1, mesh_elem[2]+1}, len = 683 shape[0]*shape[1]*shape[2]; 684 x_loc = malloc(len*num_comp_x*sizeof x_loc[0]); 685 for (CeedInt i=0; i<shape[0]; i++) { 686 for (CeedInt j=0; j<shape[1]; j++) { 687 for (CeedInt k=0; k<shape[2]; k++) { 688 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 0] = 1.*(i_rank[0]*mesh_elem[0]+i) / 689 (p[0]*mesh_elem[0]); 690 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 1] = 1.*(i_rank[1]*mesh_elem[1]+j) / 691 (p[1]*mesh_elem[1]); 692 x_loc[dim*((i*shape[1]+j)*shape[2]+k) + 2] = 1.*(i_rank[2]*mesh_elem[2]+k) / 693 (p[2]*mesh_elem[2]); 694 } 695 } 696 } 697 CeedVectorCreate(ceed, len*num_comp_x, &x_coord); 698 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_OWN_POINTER, x_loc); 699 } 700 701 // Create the QFunction that builds the operator quadrature data 702 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].setup_geo, 703 bp_options[bp_choice].setup_geo_loc, &qf_setup_geo); 704 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 705 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x*dim, CEED_EVAL_GRAD); 706 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 707 CeedQFunctionAddOutput(qf_setup_geo, "q_data", 708 bp_options[bp_choice].q_data_size, 709 CEED_EVAL_NONE); 710 711 // Create the QFunction that sets up the RHS and true solution 712 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].setup_rhs, 713 bp_options[bp_choice].setup_rhs_loc, &qf_setup_rhs); 714 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 715 CeedQFunctionAddInput(qf_setup_rhs, "q_data", bp_options[bp_choice].q_data_size, 716 CEED_EVAL_NONE); 717 CeedQFunctionAddOutput(qf_setup_rhs, "true_soln", num_comp_u, CEED_EVAL_NONE); 718 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 719 720 // Set up PDE operator 721 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].apply, 722 bp_options[bp_choice].apply_loc, &qf_apply); 723 // Add inputs and outputs 724 CeedInt in_scale = bp_options[bp_choice].in_mode==CEED_EVAL_GRAD ? 3 : 1; 725 CeedInt out_scale = bp_options[bp_choice].out_mode==CEED_EVAL_GRAD ? 3 : 1; 726 CeedQFunctionAddInput(qf_apply, "u", num_comp_u*in_scale, 727 bp_options[bp_choice].in_mode); 728 CeedQFunctionAddInput(qf_apply, "q_data", bp_options[bp_choice].q_data_size, 729 CEED_EVAL_NONE); 730 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u*out_scale, 731 bp_options[bp_choice].out_mode); 732 733 // Create the error qfunction 734 CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].error, 735 bp_options[bp_choice].error_loc, &qf_error); 736 CeedQFunctionAddInput(qf_error, "u", num_comp_u, CEED_EVAL_INTERP); 737 CeedQFunctionAddInput(qf_error, "true_soln", num_comp_u, CEED_EVAL_NONE); 738 CeedQFunctionAddOutput(qf_error, "error", num_comp_u, CEED_EVAL_NONE); 739 740 // Create the persistent vectors that will be needed in setup 741 CeedInt num_qpts; 742 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 743 CeedVectorCreate(ceed, bp_options[bp_choice].q_data_size*num_elem*num_qpts, 744 &q_data); 745 CeedVectorCreate(ceed, num_elem*num_qpts*num_comp_u, &target); 746 CeedVectorCreate(ceed, l_size*num_comp_u, &rhs_ceed); 747 748 // Create the operator that builds the quadrature data for the ceed operator 749 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, 750 CEED_QFUNCTION_NONE, &op_setup_geo); 751 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, 752 CEED_VECTOR_ACTIVE); 753 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, 754 CEED_VECTOR_ACTIVE); 755 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, 756 CEED_VECTOR_NONE); 757 CeedOperatorSetField(op_setup_geo, "q_data", elem_restr_qd_i, 758 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 759 760 // Create the operator that builds the RHS and true solution 761 CeedOperatorCreate(ceed, qf_setup_rhs, CEED_QFUNCTION_NONE, 762 CEED_QFUNCTION_NONE, &op_setup_rhs); 763 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, 764 CEED_VECTOR_ACTIVE); 765 CeedOperatorSetField(op_setup_rhs, "q_data", elem_restr_qd_i, 766 CEED_BASIS_COLLOCATED, 767 q_data); 768 CeedOperatorSetField(op_setup_rhs, "true_soln", elem_restr_u_i, 769 CEED_BASIS_COLLOCATED, target); 770 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, 771 CEED_VECTOR_ACTIVE); 772 773 // Create the mass or diff operator 774 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 775 &op_apply); 776 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 777 CeedOperatorSetField(op_apply, "q_data", elem_restr_qd_i, CEED_BASIS_COLLOCATED, 778 q_data); 779 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 780 781 // Create the error operator 782 CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 783 &op_error); 784 CeedOperatorSetField(op_error, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 785 CeedOperatorSetField(op_error, "true_soln", elem_restr_u_i, 786 CEED_BASIS_COLLOCATED, target); 787 CeedOperatorSetField(op_error, "error", elem_restr_u_i, CEED_BASIS_COLLOCATED, 788 CEED_VECTOR_ACTIVE); 789 790 // Set up Mat 791 ierr = PetscMalloc1(1, &user); CHKERRQ(ierr); 792 user->comm = comm; 793 user->l_to_g = l_to_g; 794 if (bp_choice != CEED_BP1 && bp_choice != CEED_BP2) { 795 user->l_to_g_0 = l_to_g_0; 796 user->g_to_g_D = g_to_g_D; 797 } 798 user->X_loc = X_loc; 799 ierr = VecDuplicate(X_loc, &user->Y_loc); CHKERRQ(ierr); 800 CeedVectorCreate(ceed, l_size*num_comp_u, &user->x_ceed); 801 CeedVectorCreate(ceed, l_size*num_comp_u, &user->y_ceed); 802 user->op = op_apply; 803 user->q_data = q_data; 804 user->ceed = ceed; 805 806 ierr = MatCreateShell(comm, m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 807 m_nodes[0]*m_nodes[1]*m_nodes[2]*num_comp_u, 808 PETSC_DECIDE, PETSC_DECIDE, user, &mat); CHKERRQ(ierr); 809 if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) { 810 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Mass); 811 CHKERRQ(ierr); 812 } else { 813 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Diff); 814 CHKERRQ(ierr); 815 } 816 ierr = VecGetType(X, &vec_type); CHKERRQ(ierr); 817 ierr = MatShellSetVecType(mat, vec_type); CHKERRQ(ierr); 818 819 // Get RHS vector 820 ierr = VecDuplicate(X, &rhs); CHKERRQ(ierr); 821 ierr = VecDuplicate(X_loc, &rhs_loc); CHKERRQ(ierr); 822 ierr = VecZeroEntries(rhs_loc); CHKERRQ(ierr); 823 ierr = VecGetArrayAndMemType(rhs_loc, &r, &mem_type); CHKERRQ(ierr); 824 CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r); 825 826 // Setup q_data, rhs, and target 827 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 828 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 829 CeedVectorDestroy(&x_coord); 830 831 // Gather RHS 832 ierr = CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL); CHKERRQ(ierr); 833 ierr = VecRestoreArrayAndMemType(rhs_loc, &r); CHKERRQ(ierr); 834 ierr = VecZeroEntries(rhs); CHKERRQ(ierr); 835 ierr = VecScatterBegin(l_to_g, rhs_loc, rhs, ADD_VALUES, SCATTER_FORWARD); 836 CHKERRQ(ierr); 837 ierr = VecScatterEnd(l_to_g, rhs_loc, rhs, ADD_VALUES, SCATTER_FORWARD); 838 CHKERRQ(ierr); 839 CeedVectorDestroy(&rhs_ceed); 840 841 ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr); 842 { 843 PC pc; 844 ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); 845 if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) { 846 ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr); 847 ierr = PCJacobiSetType(pc, PC_JACOBI_ROWSUM); CHKERRQ(ierr); 848 } else { 849 ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr); 850 } 851 ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr); 852 ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr); 853 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 854 PETSC_DEFAULT); CHKERRQ(ierr); 855 } 856 ierr = KSPSetOperators(ksp, mat, mat); CHKERRQ(ierr); 857 // First run's performance log is not considered for benchmarking purposes 858 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1); 859 CHKERRQ(ierr); 860 my_rt_start = MPI_Wtime(); 861 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 862 my_rt = MPI_Wtime() - my_rt_start; 863 ierr = MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, comm); 864 CHKERRQ(ierr); 865 // Set maxits based on first iteration timing 866 if (my_rt > 0.02) { 867 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 868 ksp_max_it_clip[0]); CHKERRQ(ierr); 869 } else { 870 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 871 ksp_max_it_clip[1]); CHKERRQ(ierr); 872 } 873 ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); 874 875 // Timed solve 876 ierr = VecZeroEntries(X); CHKERRQ(ierr); 877 ierr = PetscBarrier((PetscObject)ksp); CHKERRQ(ierr); 878 879 // -- Performance logging 880 ierr = PetscLogStageRegister("Solve Stage", &solve_stage); CHKERRQ(ierr); 881 ierr = PetscLogStagePush(solve_stage); CHKERRQ(ierr); 882 883 // -- Solve 884 my_rt_start = MPI_Wtime(); 885 ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr); 886 my_rt = MPI_Wtime() - my_rt_start; 887 888 // -- Performance logging 889 ierr = PetscLogStagePop(); 890 891 // Output results 892 { 893 KSPType ksp_type; 894 KSPConvergedReason reason; 895 PetscReal rnorm; 896 PetscInt its; 897 ierr = KSPGetType(ksp, &ksp_type); CHKERRQ(ierr); 898 ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr); 899 ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr); 900 ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr); 901 if (!test_mode || reason < 0 || rnorm > 1e-8) { 902 ierr = PetscPrintf(comm, 903 " KSP:\n" 904 " KSP Type : %s\n" 905 " KSP Convergence : %s\n" 906 " Total KSP Iterations : %D\n" 907 " Final rnorm : %e\n", 908 ksp_type, KSPConvergedReasons[reason], its, 909 (double)rnorm); CHKERRQ(ierr); 910 } 911 if (!test_mode) { 912 ierr = PetscPrintf(comm," Performance:\n"); CHKERRQ(ierr); 913 } 914 { 915 PetscReal max_error; 916 ierr = ComputeErrorMax(user, op_error, X, target, &max_error); 917 CHKERRQ(ierr); 918 PetscReal tol = 5e-2; 919 if (!test_mode || max_error > tol) { 920 ierr = MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, comm); 921 CHKERRQ(ierr); 922 ierr = MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, comm); 923 CHKERRQ(ierr); 924 ierr = PetscPrintf(comm, 925 " Pointwise Error (max) : %e\n" 926 " CG Solve Time : %g (%g) sec\n", 927 (double)max_error, rt_max, rt_min); CHKERRQ(ierr); 928 } 929 } 930 if (!test_mode) { 931 ierr = PetscPrintf(comm, 932 " DoFs/Sec in CG : %g (%g) million\n", 933 1e-6*gsize*its/rt_max, 934 1e-6*gsize*its/rt_min); CHKERRQ(ierr); 935 } 936 } 937 938 if (write_solution) { 939 PetscViewer vtk_viewer_soln; 940 941 ierr = PetscViewerCreate(comm, &vtk_viewer_soln); CHKERRQ(ierr); 942 ierr = PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK); CHKERRQ(ierr); 943 ierr = PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu"); CHKERRQ(ierr); 944 ierr = VecView(X, vtk_viewer_soln); CHKERRQ(ierr); 945 ierr = PetscViewerDestroy(&vtk_viewer_soln); CHKERRQ(ierr); 946 } 947 948 ierr = VecDestroy(&rhs); CHKERRQ(ierr); 949 ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr); 950 ierr = VecDestroy(&X); CHKERRQ(ierr); 951 ierr = VecDestroy(&user->X_loc); CHKERRQ(ierr); 952 ierr = VecDestroy(&user->Y_loc); CHKERRQ(ierr); 953 ierr = VecScatterDestroy(&l_to_g); CHKERRQ(ierr); 954 ierr = VecScatterDestroy(&l_to_g_0); CHKERRQ(ierr); 955 ierr = VecScatterDestroy(&g_to_g_D); CHKERRQ(ierr); 956 ierr = MatDestroy(&mat); CHKERRQ(ierr); 957 ierr = KSPDestroy(&ksp); CHKERRQ(ierr); 958 959 CeedVectorDestroy(&user->x_ceed); 960 CeedVectorDestroy(&user->y_ceed); 961 CeedVectorDestroy(&user->q_data); 962 CeedVectorDestroy(&target); 963 CeedOperatorDestroy(&op_setup_geo); 964 CeedOperatorDestroy(&op_setup_rhs); 965 CeedOperatorDestroy(&op_apply); 966 CeedOperatorDestroy(&op_error); 967 CeedElemRestrictionDestroy(&elem_restr_u); 968 CeedElemRestrictionDestroy(&elem_restr_x); 969 CeedElemRestrictionDestroy(&elem_restr_u_i); 970 CeedElemRestrictionDestroy(&elem_restr_qd_i); 971 CeedQFunctionDestroy(&qf_setup_geo); 972 CeedQFunctionDestroy(&qf_setup_rhs); 973 CeedQFunctionDestroy(&qf_apply); 974 CeedQFunctionDestroy(&qf_error); 975 CeedBasisDestroy(&basis_u); 976 CeedBasisDestroy(&basis_x); 977 CeedDestroy(&ceed); 978 ierr = PetscFree(user); CHKERRQ(ierr); 979 return PetscFinalize(); 980 } 981