1 static char help[] = "Poiseuille flow problem. Viscous, laminar flow in a 2D channel with parabolic velocity\n\ 2 profile and linear pressure drop, exact solution of the 2D Stokes\n"; 3 4 /* 5 M A R I T I M E R E S E A R C H I N S T I T U T E N E T H E R L A N D S 6 author : Christiaan M. Klaij 7 8 Poiseuille flow problem. 9 10 Viscous, laminar flow in a 2D channel with parabolic velocity 11 profile and linear pressure drop, exact solution of the 2D Stokes 12 equations. 13 14 Discretized with the cell-centered finite-volume method on a 15 Cartesian grid with co-located variables. Variables ordered as 16 [u1...uN v1...vN p1...pN]^T. Matrix [A00 A01; A10, A11] solved with 17 PCFIELDSPLIT. Lower factorization is used to mimic the Semi-Implicit 18 Method for Pressure Linked Equations (SIMPLE) used as preconditioner 19 instead of solver. 20 21 Disclaimer: does not contain the pressure-weighed interpolation 22 method needed to suppress spurious pressure modes in real-life 23 problems. 24 25 Usage: 26 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_1_pc_type none 27 28 Runs with PCFIELDSPLIT on 32x48 grid, no PC for the Schur 29 complement because A11 is zero. FGMRES is needed because 30 PCFIELDSPLIT is a variable preconditioner. 31 32 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 33 34 Same as above but with user defined PC for the true Schur 35 complement. PC based on the SIMPLE-type approximation (inverse of 36 A00 approximated by inverse of its diagonal). 37 38 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_ksp 39 40 Replace the true Schur complement with a user defined Schur 41 complement based on the SIMPLE-type approximation. Same matrix is 42 used as PC. 43 44 mpiexec -n 2 ./ex70 -nx 32 -ny 48 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi 45 46 Out-of-the-box SIMPLE-type preconditioning. The major advantage 47 is that the user neither needs to provide the approximation of 48 the Schur complement, nor the corresponding preconditioner. 49 */ 50 51 #include <petscksp.h> 52 53 typedef struct { 54 PetscBool userPC, userKSP, matsymmetric; /* user defined preconditioner and matrix for the Schur complement */ 55 PetscInt nx, ny; /* nb of cells in x- and y-direction */ 56 PetscReal hx, hy; /* mesh size in x- and y-direction */ 57 Mat A; /* block matrix */ 58 Mat subA[4]; /* the four blocks */ 59 Mat myS; /* the approximation of the Schur complement */ 60 Vec x, b, y; /* solution, rhs and temporary vector */ 61 IS isg[2]; /* index sets of split "0" and "1" */ 62 } Stokes; 63 64 PetscErrorCode StokesSetupMatBlock00(Stokes *); /* setup the block Q */ 65 PetscErrorCode StokesSetupMatBlock01(Stokes *); /* setup the block G */ 66 PetscErrorCode StokesSetupMatBlock10(Stokes *); /* setup the block D (equal to the transpose of G) */ 67 PetscErrorCode StokesSetupMatBlock11(Stokes *); /* setup the block C (equal to zero) */ 68 69 PetscErrorCode StokesGetPosition(Stokes *, PetscInt, PetscInt *, PetscInt *); /* row number j*nx+i corresponds to position (i,j) in grid */ 70 71 PetscErrorCode StokesStencilLaplacian(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Laplacian operator */ 72 PetscErrorCode StokesStencilGradientX(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Gradient operator (x-component) */ 73 PetscErrorCode StokesStencilGradientY(Stokes *, PetscInt, PetscInt, PetscInt *, PetscInt *, PetscScalar *); /* stencil of the Gradient operator (y-component) */ 74 75 PetscErrorCode StokesRhs(Stokes *); /* rhs vector */ 76 PetscErrorCode StokesRhsMomX(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of velocity (x-component) */ 77 PetscErrorCode StokesRhsMomY(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of velocity (y-component) */ 78 PetscErrorCode StokesRhsMass(Stokes *, PetscInt, PetscInt, PetscScalar *); /* right-hand side of pressure */ 79 80 PetscErrorCode StokesSetupApproxSchur(Stokes *); /* approximation of the Schur complement */ 81 82 PetscErrorCode StokesExactSolution(Stokes *); /* exact solution vector */ 83 PetscErrorCode StokesWriteSolution(Stokes *); /* write solution to file */ 84 85 /* exact solution for the velocity (x-component, y-component is zero) */ 86 PetscScalar StokesExactVelocityX(const PetscScalar y) 87 { 88 return 4.0 * y * (1.0 - y); 89 } 90 91 /* exact solution for the pressure */ 92 PetscScalar StokesExactPressure(const PetscScalar x) 93 { 94 return 8.0 * (2.0 - x); 95 } 96 97 PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp) 98 { 99 KSP *subksp; 100 PC pc; 101 PetscInt n = 1; 102 103 PetscFunctionBeginUser; 104 PetscCall(KSPGetPC(ksp, &pc)); 105 PetscCall(PCFieldSplitSetIS(pc, "0", s->isg[0])); 106 PetscCall(PCFieldSplitSetIS(pc, "1", s->isg[1])); 107 if (s->userPC) PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS)); 108 if (s->userKSP) { 109 PetscCall(PCSetUp(pc)); 110 PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp)); 111 PetscCall(KSPSetOperators(subksp[1], s->myS, s->myS)); 112 PetscCall(PetscFree(subksp)); 113 } 114 PetscFunctionReturn(PETSC_SUCCESS); 115 } 116 117 PetscErrorCode StokesWriteSolution(Stokes *s) 118 { 119 PetscMPIInt size; 120 PetscInt n, i, j; 121 const PetscScalar *array; 122 123 PetscFunctionBeginUser; 124 /* write data (*warning* only works sequential) */ 125 PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size)); 126 if (size == 1) { 127 PetscViewer viewer; 128 PetscCall(VecGetArrayRead(s->x, &array)); 129 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer)); 130 PetscCall(PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n")); 131 for (j = 0; j < s->ny; j++) { 132 for (i = 0; i < s->nx; i++) { 133 n = j * s->nx + i; 134 PetscCall(PetscViewerASCIIPrintf(viewer, "%.12g %.12g %.12g %.12g %.12g\n", (double)(i * s->hx + s->hx / 2), (double)(j * s->hy + s->hy / 2), (double)PetscRealPart(array[n]), (double)PetscRealPart(array[n + s->nx * s->ny]), (double)PetscRealPart(array[n + 2 * s->nx * s->ny]))); 135 } 136 } 137 PetscCall(VecRestoreArrayRead(s->x, &array)); 138 PetscCall(PetscViewerDestroy(&viewer)); 139 } 140 PetscFunctionReturn(PETSC_SUCCESS); 141 } 142 143 PetscErrorCode StokesSetupIndexSets(Stokes *s) 144 { 145 PetscFunctionBeginUser; 146 /* the two index sets */ 147 PetscCall(MatNestGetISs(s->A, s->isg, NULL)); 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150 151 PetscErrorCode StokesSetupVectors(Stokes *s) 152 { 153 PetscFunctionBeginUser; 154 /* solution vector x */ 155 PetscCall(VecCreate(PETSC_COMM_WORLD, &s->x)); 156 PetscCall(VecSetSizes(s->x, PETSC_DECIDE, 3 * s->nx * s->ny)); 157 PetscCall(VecSetType(s->x, VECMPI)); 158 159 /* exact solution y */ 160 PetscCall(VecDuplicate(s->x, &s->y)); 161 PetscCall(StokesExactSolution(s)); 162 163 /* rhs vector b */ 164 PetscCall(VecDuplicate(s->x, &s->b)); 165 PetscCall(StokesRhs(s)); 166 PetscFunctionReturn(PETSC_SUCCESS); 167 } 168 169 PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j) 170 { 171 PetscInt n; 172 173 PetscFunctionBeginUser; 174 /* cell number n=j*nx+i has position (i,j) in grid */ 175 n = row % (s->nx * s->ny); 176 *i = n % s->nx; 177 *j = (n - (*i)) / s->nx; 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 PetscErrorCode StokesExactSolution(Stokes *s) 182 { 183 PetscInt row, start, end, i, j; 184 PetscScalar val; 185 Vec y0, y1; 186 187 PetscFunctionBeginUser; 188 /* velocity part */ 189 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0)); 190 PetscCall(VecGetOwnershipRange(y0, &start, &end)); 191 for (row = start; row < end; row++) { 192 PetscCall(StokesGetPosition(s, row, &i, &j)); 193 if (row < s->nx * s->ny) { 194 val = StokesExactVelocityX(j * s->hy + s->hy / 2); 195 } else { 196 val = 0; 197 } 198 PetscCall(VecSetValue(y0, row, val, INSERT_VALUES)); 199 } 200 PetscCall(VecAssemblyBegin(y0)); 201 PetscCall(VecAssemblyEnd(y0)); 202 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0)); 203 204 /* pressure part */ 205 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1)); 206 PetscCall(VecGetOwnershipRange(y1, &start, &end)); 207 for (row = start; row < end; row++) { 208 PetscCall(StokesGetPosition(s, row, &i, &j)); 209 val = StokesExactPressure(i * s->hx + s->hx / 2); 210 PetscCall(VecSetValue(y1, row, val, INSERT_VALUES)); 211 } 212 PetscCall(VecAssemblyBegin(y1)); 213 PetscCall(VecAssemblyEnd(y1)); 214 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1)); 215 PetscFunctionReturn(PETSC_SUCCESS); 216 } 217 218 PetscErrorCode StokesRhs(Stokes *s) 219 { 220 PetscInt row, start, end, i, j; 221 PetscScalar val; 222 Vec b0, b1; 223 224 PetscFunctionBeginUser; 225 /* velocity part */ 226 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0)); 227 PetscCall(VecGetOwnershipRange(b0, &start, &end)); 228 for (row = start; row < end; row++) { 229 PetscCall(StokesGetPosition(s, row, &i, &j)); 230 if (row < s->nx * s->ny) { 231 PetscCall(StokesRhsMomX(s, i, j, &val)); 232 } else { 233 PetscCall(StokesRhsMomY(s, i, j, &val)); 234 } 235 PetscCall(VecSetValue(b0, row, val, INSERT_VALUES)); 236 } 237 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0)); 238 239 /* pressure part */ 240 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1)); 241 PetscCall(VecGetOwnershipRange(b1, &start, &end)); 242 for (row = start; row < end; row++) { 243 PetscCall(StokesGetPosition(s, row, &i, &j)); 244 PetscCall(StokesRhsMass(s, i, j, &val)); 245 if (s->matsymmetric) val = -1.0 * val; 246 PetscCall(VecSetValue(b1, row, val, INSERT_VALUES)); 247 } 248 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); 249 PetscFunctionReturn(PETSC_SUCCESS); 250 } 251 252 PetscErrorCode StokesSetupMatBlock00(Stokes *s) 253 { 254 PetscInt row, start, end, sz, i, j; 255 PetscInt cols[5]; 256 PetscScalar vals[5]; 257 258 PetscFunctionBeginUser; 259 /* A[0] is 2N-by-2N */ 260 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[0])); 261 PetscCall(MatSetOptionsPrefix(s->subA[0], "a00_")); 262 PetscCall(MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny)); 263 PetscCall(MatSetType(s->subA[0], MATMPIAIJ)); 264 PetscCall(MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL)); 265 PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end)); 266 267 for (row = start; row < end; row++) { 268 PetscCall(StokesGetPosition(s, row, &i, &j)); 269 /* first part: rows 0 to (nx*ny-1) */ 270 PetscCall(StokesStencilLaplacian(s, i, j, &sz, cols, vals)); 271 /* second part: rows (nx*ny) to (2*nx*ny-1) */ 272 if (row >= s->nx * s->ny) { 273 for (i = 0; i < sz; i++) cols[i] += s->nx * s->ny; 274 } 275 for (i = 0; i < sz; i++) vals[i] = -1.0 * vals[i]; /* dynamic viscosity coef mu=-1 */ 276 PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES)); 277 } 278 PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY)); 279 PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY)); 280 PetscFunctionReturn(PETSC_SUCCESS); 281 } 282 283 PetscErrorCode StokesSetupMatBlock01(Stokes *s) 284 { 285 PetscInt row, start, end, sz, i, j; 286 PetscInt cols[5]; 287 PetscScalar vals[5]; 288 289 PetscFunctionBeginUser; 290 /* A[1] is 2N-by-N */ 291 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1])); 292 PetscCall(MatSetOptionsPrefix(s->subA[1], "a01_")); 293 PetscCall(MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny)); 294 PetscCall(MatSetType(s->subA[1], MATMPIAIJ)); 295 PetscCall(MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL)); 296 PetscCall(MatGetOwnershipRange(s->subA[1], &start, &end)); 297 298 PetscCall(MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); 299 300 for (row = start; row < end; row++) { 301 PetscCall(StokesGetPosition(s, row, &i, &j)); 302 /* first part: rows 0 to (nx*ny-1) */ 303 if (row < s->nx * s->ny) { 304 PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals)); 305 } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */ 306 PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals)); 307 } 308 PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES)); 309 } 310 PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY)); 311 PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY)); 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314 315 PetscErrorCode StokesSetupMatBlock10(Stokes *s) 316 { 317 PetscFunctionBeginUser; 318 /* A[2] is minus transpose of A[1] */ 319 PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2])); 320 if (!s->matsymmetric) PetscCall(MatScale(s->subA[2], -1.0)); 321 PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_")); 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 PetscErrorCode StokesSetupMatBlock11(Stokes *s) 326 { 327 PetscFunctionBeginUser; 328 /* A[3] is N-by-N null matrix */ 329 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[3])); 330 PetscCall(MatSetOptionsPrefix(s->subA[3], "a11_")); 331 PetscCall(MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny)); 332 PetscCall(MatSetType(s->subA[3], MATMPIAIJ)); 333 PetscCall(MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL)); 334 PetscCall(MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY)); 335 PetscCall(MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY)); 336 PetscFunctionReturn(PETSC_SUCCESS); 337 } 338 339 PetscErrorCode StokesSetupApproxSchur(Stokes *s) 340 { 341 Vec diag; 342 343 PetscFunctionBeginUser; 344 /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */ 345 /* note: A11 is zero */ 346 /* note: in real life this matrix would be build directly, */ 347 /* i.e. without MatMatMult */ 348 349 /* inverse of diagonal of A00 */ 350 PetscCall(VecCreate(PETSC_COMM_WORLD, &diag)); 351 PetscCall(VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny)); 352 PetscCall(VecSetType(diag, VECMPI)); 353 PetscCall(MatGetDiagonal(s->subA[0], diag)); 354 PetscCall(VecReciprocal(diag)); 355 356 /* compute: - A10 inv(DIAGFORM(A00)) A01 */ 357 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); /* (*warning* overwrites subA[1]) */ 358 PetscCall(MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DEFAULT, &s->myS)); 359 PetscCall(MatScale(s->myS, -1.0)); 360 361 /* restore A10 */ 362 PetscCall(MatGetDiagonal(s->subA[0], diag)); 363 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); 364 PetscCall(VecDestroy(&diag)); 365 PetscFunctionReturn(PETSC_SUCCESS); 366 } 367 368 PetscErrorCode StokesSetupMatrix(Stokes *s) 369 { 370 PetscFunctionBeginUser; 371 PetscCall(StokesSetupMatBlock00(s)); 372 PetscCall(StokesSetupMatBlock01(s)); 373 PetscCall(StokesSetupMatBlock10(s)); 374 PetscCall(StokesSetupMatBlock11(s)); 375 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A)); 376 PetscCall(StokesSetupApproxSchur(s)); 377 PetscFunctionReturn(PETSC_SUCCESS); 378 } 379 380 PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 381 { 382 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1, s2 = p - s->nx, n = p + s->nx; 383 PetscScalar ae = s->hy / s->hx, aeb = 0; 384 PetscScalar aw = s->hy / s->hx, awb = s->hy / (s->hx / 2); 385 PetscScalar as = s->hx / s->hy, asb = s->hx / (s->hy / 2); 386 PetscScalar an = s->hx / s->hy, anb = s->hx / (s->hy / 2); 387 388 PetscFunctionBeginUser; 389 if (i == 0 && j == 0) { /* south-west corner */ 390 *sz = 3; 391 cols[0] = p; 392 vals[0] = -(ae + awb + asb + an); 393 cols[1] = e; 394 vals[1] = ae; 395 cols[2] = n; 396 vals[2] = an; 397 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ 398 *sz = 3; 399 cols[0] = s2; 400 vals[0] = as; 401 cols[1] = p; 402 vals[1] = -(ae + awb + as + anb); 403 cols[2] = e; 404 vals[2] = ae; 405 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ 406 *sz = 3; 407 cols[0] = w; 408 vals[0] = aw; 409 cols[1] = p; 410 vals[1] = -(aeb + aw + asb + an); 411 cols[2] = n; 412 vals[2] = an; 413 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ 414 *sz = 3; 415 cols[0] = s2; 416 vals[0] = as; 417 cols[1] = w; 418 vals[1] = aw; 419 cols[2] = p; 420 vals[2] = -(aeb + aw + as + anb); 421 } else if (i == 0) { /* west boundary */ 422 *sz = 4; 423 cols[0] = s2; 424 vals[0] = as; 425 cols[1] = p; 426 vals[1] = -(ae + awb + as + an); 427 cols[2] = e; 428 vals[2] = ae; 429 cols[3] = n; 430 vals[3] = an; 431 } else if (i == s->nx - 1) { /* east boundary */ 432 *sz = 4; 433 cols[0] = s2; 434 vals[0] = as; 435 cols[1] = w; 436 vals[1] = aw; 437 cols[2] = p; 438 vals[2] = -(aeb + aw + as + an); 439 cols[3] = n; 440 vals[3] = an; 441 } else if (j == 0) { /* south boundary */ 442 *sz = 4; 443 cols[0] = w; 444 vals[0] = aw; 445 cols[1] = p; 446 vals[1] = -(ae + aw + asb + an); 447 cols[2] = e; 448 vals[2] = ae; 449 cols[3] = n; 450 vals[3] = an; 451 } else if (j == s->ny - 1) { /* north boundary */ 452 *sz = 4; 453 cols[0] = s2; 454 vals[0] = as; 455 cols[1] = w; 456 vals[1] = aw; 457 cols[2] = p; 458 vals[2] = -(ae + aw + as + anb); 459 cols[3] = e; 460 vals[3] = ae; 461 } else { /* interior */ 462 *sz = 5; 463 cols[0] = s2; 464 vals[0] = as; 465 cols[1] = w; 466 vals[1] = aw; 467 cols[2] = p; 468 vals[2] = -(ae + aw + as + an); 469 cols[3] = e; 470 vals[3] = ae; 471 cols[4] = n; 472 vals[4] = an; 473 } 474 PetscFunctionReturn(PETSC_SUCCESS); 475 } 476 477 PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 478 { 479 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1; 480 PetscScalar ae = s->hy / 2, aeb = s->hy; 481 PetscScalar aw = -s->hy / 2, awb = 0; 482 483 PetscFunctionBeginUser; 484 if (i == 0 && j == 0) { /* south-west corner */ 485 *sz = 2; 486 cols[0] = p; 487 vals[0] = -(ae + awb); 488 cols[1] = e; 489 vals[1] = ae; 490 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ 491 *sz = 2; 492 cols[0] = p; 493 vals[0] = -(ae + awb); 494 cols[1] = e; 495 vals[1] = ae; 496 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ 497 *sz = 2; 498 cols[0] = w; 499 vals[0] = aw; 500 cols[1] = p; 501 vals[1] = -(aeb + aw); 502 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ 503 *sz = 2; 504 cols[0] = w; 505 vals[0] = aw; 506 cols[1] = p; 507 vals[1] = -(aeb + aw); 508 } else if (i == 0) { /* west boundary */ 509 *sz = 2; 510 cols[0] = p; 511 vals[0] = -(ae + awb); 512 cols[1] = e; 513 vals[1] = ae; 514 } else if (i == s->nx - 1) { /* east boundary */ 515 *sz = 2; 516 cols[0] = w; 517 vals[0] = aw; 518 cols[1] = p; 519 vals[1] = -(aeb + aw); 520 } else if (j == 0) { /* south boundary */ 521 *sz = 3; 522 cols[0] = w; 523 vals[0] = aw; 524 cols[1] = p; 525 vals[1] = -(ae + aw); 526 cols[2] = e; 527 vals[2] = ae; 528 } else if (j == s->ny - 1) { /* north boundary */ 529 *sz = 3; 530 cols[0] = w; 531 vals[0] = aw; 532 cols[1] = p; 533 vals[1] = -(ae + aw); 534 cols[2] = e; 535 vals[2] = ae; 536 } else { /* interior */ 537 *sz = 3; 538 cols[0] = w; 539 vals[0] = aw; 540 cols[1] = p; 541 vals[1] = -(ae + aw); 542 cols[2] = e; 543 vals[2] = ae; 544 } 545 PetscFunctionReturn(PETSC_SUCCESS); 546 } 547 548 PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 549 { 550 PetscInt p = j * s->nx + i, s2 = p - s->nx, n = p + s->nx; 551 PetscScalar as = -s->hx / 2, asb = 0; 552 PetscScalar an = s->hx / 2, anb = 0; 553 554 PetscFunctionBeginUser; 555 if (i == 0 && j == 0) { /* south-west corner */ 556 *sz = 2; 557 cols[0] = p; 558 vals[0] = -(asb + an); 559 cols[1] = n; 560 vals[1] = an; 561 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ 562 *sz = 2; 563 cols[0] = s2; 564 vals[0] = as; 565 cols[1] = p; 566 vals[1] = -(as + anb); 567 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ 568 *sz = 2; 569 cols[0] = p; 570 vals[0] = -(asb + an); 571 cols[1] = n; 572 vals[1] = an; 573 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ 574 *sz = 2; 575 cols[0] = s2; 576 vals[0] = as; 577 cols[1] = p; 578 vals[1] = -(as + anb); 579 } else if (i == 0) { /* west boundary */ 580 *sz = 3; 581 cols[0] = s2; 582 vals[0] = as; 583 cols[1] = p; 584 vals[1] = -(as + an); 585 cols[2] = n; 586 vals[2] = an; 587 } else if (i == s->nx - 1) { /* east boundary */ 588 *sz = 3; 589 cols[0] = s2; 590 vals[0] = as; 591 cols[1] = p; 592 vals[1] = -(as + an); 593 cols[2] = n; 594 vals[2] = an; 595 } else if (j == 0) { /* south boundary */ 596 *sz = 2; 597 cols[0] = p; 598 vals[0] = -(asb + an); 599 cols[1] = n; 600 vals[1] = an; 601 } else if (j == s->ny - 1) { /* north boundary */ 602 *sz = 2; 603 cols[0] = s2; 604 vals[0] = as; 605 cols[1] = p; 606 vals[1] = -(as + anb); 607 } else { /* interior */ 608 *sz = 3; 609 cols[0] = s2; 610 vals[0] = as; 611 cols[1] = p; 612 vals[1] = -(as + an); 613 cols[2] = n; 614 vals[2] = an; 615 } 616 PetscFunctionReturn(PETSC_SUCCESS); 617 } 618 619 PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 620 { 621 PetscScalar y = j * s->hy + s->hy / 2; 622 PetscScalar awb = s->hy / (s->hx / 2); 623 624 PetscFunctionBeginUser; 625 if (i == 0) { /* west boundary */ 626 *val = awb * StokesExactVelocityX(y); 627 } else { 628 *val = 0.0; 629 } 630 PetscFunctionReturn(PETSC_SUCCESS); 631 } 632 633 PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 634 { 635 PetscFunctionBeginUser; 636 *val = 0.0; 637 PetscFunctionReturn(PETSC_SUCCESS); 638 } 639 640 PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 641 { 642 PetscScalar y = j * s->hy + s->hy / 2; 643 PetscScalar aeb = s->hy; 644 645 PetscFunctionBeginUser; 646 if (i == 0) { /* west boundary */ 647 *val = aeb * StokesExactVelocityX(y); 648 } else { 649 *val = 0.0; 650 } 651 PetscFunctionReturn(PETSC_SUCCESS); 652 } 653 654 PetscErrorCode StokesCalcResidual(Stokes *s) 655 { 656 PetscReal val; 657 Vec b0, b1; 658 659 PetscFunctionBeginUser; 660 /* residual Ax-b (*warning* overwrites b) */ 661 PetscCall(VecScale(s->b, -1.0)); 662 PetscCall(MatMultAdd(s->A, s->x, s->b, s->b)); 663 664 /* residual velocity */ 665 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0)); 666 PetscCall(VecNorm(b0, NORM_2, &val)); 667 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual u = %g\n", (double)val)); 668 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0)); 669 670 /* residual pressure */ 671 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1)); 672 PetscCall(VecNorm(b1, NORM_2, &val)); 673 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual p = %g\n", (double)val)); 674 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); 675 676 /* total residual */ 677 PetscCall(VecNorm(s->b, NORM_2, &val)); 678 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual [u,p] = %g\n", (double)val)); 679 PetscFunctionReturn(PETSC_SUCCESS); 680 } 681 682 PetscErrorCode StokesCalcError(Stokes *s) 683 { 684 PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny); 685 PetscReal val; 686 Vec y0, y1; 687 688 PetscFunctionBeginUser; 689 /* error y-x */ 690 PetscCall(VecAXPY(s->y, -1.0, s->x)); 691 692 /* error in velocity */ 693 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0)); 694 PetscCall(VecNorm(y0, NORM_2, &val)); 695 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error u = %g\n", (double)(PetscRealPart(val / scale)))); 696 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0)); 697 698 /* error in pressure */ 699 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1)); 700 PetscCall(VecNorm(y1, NORM_2, &val)); 701 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error p = %g\n", (double)(PetscRealPart(val / scale)))); 702 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1)); 703 704 /* total error */ 705 PetscCall(VecNorm(s->y, NORM_2, &val)); 706 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error [u,p] = %g\n", (double)PetscRealPart(val / scale))); 707 PetscFunctionReturn(PETSC_SUCCESS); 708 } 709 710 int main(int argc, char **argv) 711 { 712 Stokes s; 713 KSP ksp; 714 715 PetscFunctionBeginUser; 716 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 717 s.nx = 4; 718 s.ny = 6; 719 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL)); 720 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ny", &s.ny, NULL)); 721 s.hx = 2.0 / s.nx; 722 s.hy = 1.0 / s.ny; 723 s.matsymmetric = PETSC_FALSE; 724 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL)); 725 s.userPC = s.userKSP = PETSC_FALSE; 726 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC)); 727 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP)); 728 729 PetscCall(StokesSetupMatrix(&s)); 730 PetscCall(StokesSetupIndexSets(&s)); 731 PetscCall(StokesSetupVectors(&s)); 732 733 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 734 PetscCall(KSPSetOperators(ksp, s.A, s.A)); 735 PetscCall(KSPSetFromOptions(ksp)); 736 PetscCall(StokesSetupPC(&s, ksp)); 737 PetscCall(KSPSolve(ksp, s.b, s.x)); 738 739 /* don't trust, verify! */ 740 PetscCall(StokesCalcResidual(&s)); 741 PetscCall(StokesCalcError(&s)); 742 PetscCall(StokesWriteSolution(&s)); 743 744 PetscCall(KSPDestroy(&ksp)); 745 PetscCall(MatDestroy(&s.subA[0])); 746 PetscCall(MatDestroy(&s.subA[1])); 747 PetscCall(MatDestroy(&s.subA[2])); 748 PetscCall(MatDestroy(&s.subA[3])); 749 PetscCall(MatDestroy(&s.A)); 750 PetscCall(VecDestroy(&s.x)); 751 PetscCall(VecDestroy(&s.b)); 752 PetscCall(VecDestroy(&s.y)); 753 PetscCall(MatDestroy(&s.myS)); 754 PetscCall(PetscFinalize()); 755 return 0; 756 } 757 758 /*TEST 759 760 test: 761 nsize: 2 762 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_1_pc_type none 763 764 test: 765 suffix: 2 766 nsize: 2 767 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 768 769 test: 770 suffix: 3 771 nsize: 2 772 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 773 774 test: 775 suffix: 4 776 nsize: 2 777 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_inner_pc_type jacobi -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi 778 779 test: 780 suffix: 4_pcksp 781 nsize: 2 782 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_type bjacobi -fieldsplit_1_pc_type jacobi -fieldsplit_1_inner_ksp_type preonly -fieldsplit_1_upper_ksp_type preonly -fieldsplit_1_upper_pc_type jacobi 783 784 test: 785 suffix: 5 786 nsize: 2 787 args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type gkb -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-10 788 789 test: 790 suffix: 6 791 nsize: 2 792 args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type gkb -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-10 793 794 test: 795 suffix: 7 796 nsize: 2 797 args: -nx 4 -ny 8 -mat_set_symmetric -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type gkb -pc_fieldsplit_gkb_tol 1e-4 -pc_fieldsplit_gkb_nu 5 -fieldsplit_0_ksp_type cg -fieldsplit_0_pc_type jacobi -fieldsplit_0_ksp_rtol 1e-6 798 799 TEST*/ 800