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(0); 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(0); 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(0); 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(0); 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(0); 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(VecRestoreSubVector(s->y, s->isg[0], &y0)); 201 202 /* pressure part */ 203 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1)); 204 PetscCall(VecGetOwnershipRange(y1, &start, &end)); 205 for (row = start; row < end; row++) { 206 PetscCall(StokesGetPosition(s, row, &i, &j)); 207 val = StokesExactPressure(i * s->hx + s->hx / 2); 208 PetscCall(VecSetValue(y1, row, val, INSERT_VALUES)); 209 } 210 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1)); 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode StokesRhs(Stokes *s) 215 { 216 PetscInt row, start, end, i, j; 217 PetscScalar val; 218 Vec b0, b1; 219 220 PetscFunctionBeginUser; 221 /* velocity part */ 222 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0)); 223 PetscCall(VecGetOwnershipRange(b0, &start, &end)); 224 for (row = start; row < end; row++) { 225 PetscCall(StokesGetPosition(s, row, &i, &j)); 226 if (row < s->nx * s->ny) { 227 PetscCall(StokesRhsMomX(s, i, j, &val)); 228 } else { 229 PetscCall(StokesRhsMomY(s, i, j, &val)); 230 } 231 PetscCall(VecSetValue(b0, row, val, INSERT_VALUES)); 232 } 233 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0)); 234 235 /* pressure part */ 236 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1)); 237 PetscCall(VecGetOwnershipRange(b1, &start, &end)); 238 for (row = start; row < end; row++) { 239 PetscCall(StokesGetPosition(s, row, &i, &j)); 240 PetscCall(StokesRhsMass(s, i, j, &val)); 241 if (s->matsymmetric) val = -1.0 * val; 242 PetscCall(VecSetValue(b1, row, val, INSERT_VALUES)); 243 } 244 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); 245 PetscFunctionReturn(0); 246 } 247 248 PetscErrorCode StokesSetupMatBlock00(Stokes *s) 249 { 250 PetscInt row, start, end, sz, i, j; 251 PetscInt cols[5]; 252 PetscScalar vals[5]; 253 254 PetscFunctionBeginUser; 255 /* A[0] is 2N-by-2N */ 256 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[0])); 257 PetscCall(MatSetOptionsPrefix(s->subA[0], "a00_")); 258 PetscCall(MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny)); 259 PetscCall(MatSetType(s->subA[0], MATMPIAIJ)); 260 PetscCall(MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL)); 261 PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end)); 262 263 for (row = start; row < end; row++) { 264 PetscCall(StokesGetPosition(s, row, &i, &j)); 265 /* first part: rows 0 to (nx*ny-1) */ 266 PetscCall(StokesStencilLaplacian(s, i, j, &sz, cols, vals)); 267 /* second part: rows (nx*ny) to (2*nx*ny-1) */ 268 if (row >= s->nx * s->ny) { 269 for (i = 0; i < sz; i++) cols[i] += s->nx * s->ny; 270 } 271 for (i = 0; i < sz; i++) vals[i] = -1.0 * vals[i]; /* dynamic viscosity coef mu=-1 */ 272 PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES)); 273 } 274 PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY)); 275 PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY)); 276 PetscFunctionReturn(0); 277 } 278 279 PetscErrorCode StokesSetupMatBlock01(Stokes *s) 280 { 281 PetscInt row, start, end, sz, i, j; 282 PetscInt cols[5]; 283 PetscScalar vals[5]; 284 285 PetscFunctionBeginUser; 286 /* A[1] is 2N-by-N */ 287 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1])); 288 PetscCall(MatSetOptionsPrefix(s->subA[1], "a01_")); 289 PetscCall(MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny)); 290 PetscCall(MatSetType(s->subA[1], MATMPIAIJ)); 291 PetscCall(MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL)); 292 PetscCall(MatGetOwnershipRange(s->subA[1], &start, &end)); 293 294 PetscCall(MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); 295 296 for (row = start; row < end; row++) { 297 PetscCall(StokesGetPosition(s, row, &i, &j)); 298 /* first part: rows 0 to (nx*ny-1) */ 299 if (row < s->nx * s->ny) { 300 PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals)); 301 } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */ 302 PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals)); 303 } 304 PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES)); 305 } 306 PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY)); 307 PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY)); 308 PetscFunctionReturn(0); 309 } 310 311 PetscErrorCode StokesSetupMatBlock10(Stokes *s) 312 { 313 PetscFunctionBeginUser; 314 /* A[2] is minus transpose of A[1] */ 315 PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2])); 316 if (!s->matsymmetric) PetscCall(MatScale(s->subA[2], -1.0)); 317 PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_")); 318 PetscFunctionReturn(0); 319 } 320 321 PetscErrorCode StokesSetupMatBlock11(Stokes *s) 322 { 323 PetscFunctionBeginUser; 324 /* A[3] is N-by-N null matrix */ 325 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[3])); 326 PetscCall(MatSetOptionsPrefix(s->subA[3], "a11_")); 327 PetscCall(MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny)); 328 PetscCall(MatSetType(s->subA[3], MATMPIAIJ)); 329 PetscCall(MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL)); 330 PetscCall(MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY)); 331 PetscCall(MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY)); 332 PetscFunctionReturn(0); 333 } 334 335 PetscErrorCode StokesSetupApproxSchur(Stokes *s) 336 { 337 Vec diag; 338 339 PetscFunctionBeginUser; 340 /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */ 341 /* note: A11 is zero */ 342 /* note: in real life this matrix would be build directly, */ 343 /* i.e. without MatMatMult */ 344 345 /* inverse of diagonal of A00 */ 346 PetscCall(VecCreate(PETSC_COMM_WORLD, &diag)); 347 PetscCall(VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny)); 348 PetscCall(VecSetType(diag, VECMPI)); 349 PetscCall(MatGetDiagonal(s->subA[0], diag)); 350 PetscCall(VecReciprocal(diag)); 351 352 /* compute: - A10 inv(DIAGFORM(A00)) A01 */ 353 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); /* (*warning* overwrites subA[1]) */ 354 PetscCall(MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DEFAULT, &s->myS)); 355 PetscCall(MatScale(s->myS, -1.0)); 356 357 /* restore A10 */ 358 PetscCall(MatGetDiagonal(s->subA[0], diag)); 359 PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); 360 PetscCall(VecDestroy(&diag)); 361 PetscFunctionReturn(0); 362 } 363 364 PetscErrorCode StokesSetupMatrix(Stokes *s) 365 { 366 PetscFunctionBeginUser; 367 PetscCall(StokesSetupMatBlock00(s)); 368 PetscCall(StokesSetupMatBlock01(s)); 369 PetscCall(StokesSetupMatBlock10(s)); 370 PetscCall(StokesSetupMatBlock11(s)); 371 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A)); 372 PetscCall(StokesSetupApproxSchur(s)); 373 PetscFunctionReturn(0); 374 } 375 376 PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 377 { 378 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1, s2 = p - s->nx, n = p + s->nx; 379 PetscScalar ae = s->hy / s->hx, aeb = 0; 380 PetscScalar aw = s->hy / s->hx, awb = s->hy / (s->hx / 2); 381 PetscScalar as = s->hx / s->hy, asb = s->hx / (s->hy / 2); 382 PetscScalar an = s->hx / s->hy, anb = s->hx / (s->hy / 2); 383 384 PetscFunctionBeginUser; 385 if (i == 0 && j == 0) { /* south-west corner */ 386 *sz = 3; 387 cols[0] = p; 388 vals[0] = -(ae + awb + asb + an); 389 cols[1] = e; 390 vals[1] = ae; 391 cols[2] = n; 392 vals[2] = an; 393 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ 394 *sz = 3; 395 cols[0] = s2; 396 vals[0] = as; 397 cols[1] = p; 398 vals[1] = -(ae + awb + as + anb); 399 cols[2] = e; 400 vals[2] = ae; 401 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ 402 *sz = 3; 403 cols[0] = w; 404 vals[0] = aw; 405 cols[1] = p; 406 vals[1] = -(aeb + aw + asb + an); 407 cols[2] = n; 408 vals[2] = an; 409 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ 410 *sz = 3; 411 cols[0] = s2; 412 vals[0] = as; 413 cols[1] = w; 414 vals[1] = aw; 415 cols[2] = p; 416 vals[2] = -(aeb + aw + as + anb); 417 } else if (i == 0) { /* west boundary */ 418 *sz = 4; 419 cols[0] = s2; 420 vals[0] = as; 421 cols[1] = p; 422 vals[1] = -(ae + awb + as + an); 423 cols[2] = e; 424 vals[2] = ae; 425 cols[3] = n; 426 vals[3] = an; 427 } else if (i == s->nx - 1) { /* east boundary */ 428 *sz = 4; 429 cols[0] = s2; 430 vals[0] = as; 431 cols[1] = w; 432 vals[1] = aw; 433 cols[2] = p; 434 vals[2] = -(aeb + aw + as + an); 435 cols[3] = n; 436 vals[3] = an; 437 } else if (j == 0) { /* south boundary */ 438 *sz = 4; 439 cols[0] = w; 440 vals[0] = aw; 441 cols[1] = p; 442 vals[1] = -(ae + aw + asb + an); 443 cols[2] = e; 444 vals[2] = ae; 445 cols[3] = n; 446 vals[3] = an; 447 } else if (j == s->ny - 1) { /* north boundary */ 448 *sz = 4; 449 cols[0] = s2; 450 vals[0] = as; 451 cols[1] = w; 452 vals[1] = aw; 453 cols[2] = p; 454 vals[2] = -(ae + aw + as + anb); 455 cols[3] = e; 456 vals[3] = ae; 457 } else { /* interior */ 458 *sz = 5; 459 cols[0] = s2; 460 vals[0] = as; 461 cols[1] = w; 462 vals[1] = aw; 463 cols[2] = p; 464 vals[2] = -(ae + aw + as + an); 465 cols[3] = e; 466 vals[3] = ae; 467 cols[4] = n; 468 vals[4] = an; 469 } 470 PetscFunctionReturn(0); 471 } 472 473 PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 474 { 475 PetscInt p = j * s->nx + i, w = p - 1, e = p + 1; 476 PetscScalar ae = s->hy / 2, aeb = s->hy; 477 PetscScalar aw = -s->hy / 2, awb = 0; 478 479 PetscFunctionBeginUser; 480 if (i == 0 && j == 0) { /* south-west corner */ 481 *sz = 2; 482 cols[0] = p; 483 vals[0] = -(ae + awb); 484 cols[1] = e; 485 vals[1] = ae; 486 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ 487 *sz = 2; 488 cols[0] = p; 489 vals[0] = -(ae + awb); 490 cols[1] = e; 491 vals[1] = ae; 492 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ 493 *sz = 2; 494 cols[0] = w; 495 vals[0] = aw; 496 cols[1] = p; 497 vals[1] = -(aeb + aw); 498 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ 499 *sz = 2; 500 cols[0] = w; 501 vals[0] = aw; 502 cols[1] = p; 503 vals[1] = -(aeb + aw); 504 } else if (i == 0) { /* west boundary */ 505 *sz = 2; 506 cols[0] = p; 507 vals[0] = -(ae + awb); 508 cols[1] = e; 509 vals[1] = ae; 510 } else if (i == s->nx - 1) { /* east boundary */ 511 *sz = 2; 512 cols[0] = w; 513 vals[0] = aw; 514 cols[1] = p; 515 vals[1] = -(aeb + aw); 516 } else if (j == 0) { /* south boundary */ 517 *sz = 3; 518 cols[0] = w; 519 vals[0] = aw; 520 cols[1] = p; 521 vals[1] = -(ae + aw); 522 cols[2] = e; 523 vals[2] = ae; 524 } else if (j == s->ny - 1) { /* north boundary */ 525 *sz = 3; 526 cols[0] = w; 527 vals[0] = aw; 528 cols[1] = p; 529 vals[1] = -(ae + aw); 530 cols[2] = e; 531 vals[2] = ae; 532 } else { /* interior */ 533 *sz = 3; 534 cols[0] = w; 535 vals[0] = aw; 536 cols[1] = p; 537 vals[1] = -(ae + aw); 538 cols[2] = e; 539 vals[2] = ae; 540 } 541 PetscFunctionReturn(0); 542 } 543 544 PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 545 { 546 PetscInt p = j * s->nx + i, s2 = p - s->nx, n = p + s->nx; 547 PetscScalar as = -s->hx / 2, asb = 0; 548 PetscScalar an = s->hx / 2, anb = 0; 549 550 PetscFunctionBeginUser; 551 if (i == 0 && j == 0) { /* south-west corner */ 552 *sz = 2; 553 cols[0] = p; 554 vals[0] = -(asb + an); 555 cols[1] = n; 556 vals[1] = an; 557 } else if (i == 0 && j == s->ny - 1) { /* north-west corner */ 558 *sz = 2; 559 cols[0] = s2; 560 vals[0] = as; 561 cols[1] = p; 562 vals[1] = -(as + anb); 563 } else if (i == s->nx - 1 && j == 0) { /* south-east corner */ 564 *sz = 2; 565 cols[0] = p; 566 vals[0] = -(asb + an); 567 cols[1] = n; 568 vals[1] = an; 569 } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */ 570 *sz = 2; 571 cols[0] = s2; 572 vals[0] = as; 573 cols[1] = p; 574 vals[1] = -(as + anb); 575 } else if (i == 0) { /* west boundary */ 576 *sz = 3; 577 cols[0] = s2; 578 vals[0] = as; 579 cols[1] = p; 580 vals[1] = -(as + an); 581 cols[2] = n; 582 vals[2] = an; 583 } else if (i == s->nx - 1) { /* east boundary */ 584 *sz = 3; 585 cols[0] = s2; 586 vals[0] = as; 587 cols[1] = p; 588 vals[1] = -(as + an); 589 cols[2] = n; 590 vals[2] = an; 591 } else if (j == 0) { /* south boundary */ 592 *sz = 2; 593 cols[0] = p; 594 vals[0] = -(asb + an); 595 cols[1] = n; 596 vals[1] = an; 597 } else if (j == s->ny - 1) { /* north boundary */ 598 *sz = 2; 599 cols[0] = s2; 600 vals[0] = as; 601 cols[1] = p; 602 vals[1] = -(as + anb); 603 } else { /* interior */ 604 *sz = 3; 605 cols[0] = s2; 606 vals[0] = as; 607 cols[1] = p; 608 vals[1] = -(as + an); 609 cols[2] = n; 610 vals[2] = an; 611 } 612 PetscFunctionReturn(0); 613 } 614 615 PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 616 { 617 PetscScalar y = j * s->hy + s->hy / 2; 618 PetscScalar awb = s->hy / (s->hx / 2); 619 620 PetscFunctionBeginUser; 621 if (i == 0) { /* west boundary */ 622 *val = awb * StokesExactVelocityX(y); 623 } else { 624 *val = 0.0; 625 } 626 PetscFunctionReturn(0); 627 } 628 629 PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 630 { 631 PetscFunctionBeginUser; 632 *val = 0.0; 633 PetscFunctionReturn(0); 634 } 635 636 PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 637 { 638 PetscScalar y = j * s->hy + s->hy / 2; 639 PetscScalar aeb = s->hy; 640 641 PetscFunctionBeginUser; 642 if (i == 0) { /* west boundary */ 643 *val = aeb * StokesExactVelocityX(y); 644 } else { 645 *val = 0.0; 646 } 647 PetscFunctionReturn(0); 648 } 649 650 PetscErrorCode StokesCalcResidual(Stokes *s) 651 { 652 PetscReal val; 653 Vec b0, b1; 654 655 PetscFunctionBeginUser; 656 /* residual Ax-b (*warning* overwrites b) */ 657 PetscCall(VecScale(s->b, -1.0)); 658 PetscCall(MatMultAdd(s->A, s->x, s->b, s->b)); 659 660 /* residual velocity */ 661 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0)); 662 PetscCall(VecNorm(b0, NORM_2, &val)); 663 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual u = %g\n", (double)val)); 664 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0)); 665 666 /* residual pressure */ 667 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1)); 668 PetscCall(VecNorm(b1, NORM_2, &val)); 669 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual p = %g\n", (double)val)); 670 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); 671 672 /* total residual */ 673 PetscCall(VecNorm(s->b, NORM_2, &val)); 674 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual [u,p] = %g\n", (double)val)); 675 PetscFunctionReturn(0); 676 } 677 678 PetscErrorCode StokesCalcError(Stokes *s) 679 { 680 PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny); 681 PetscReal val; 682 Vec y0, y1; 683 684 PetscFunctionBeginUser; 685 /* error y-x */ 686 PetscCall(VecAXPY(s->y, -1.0, s->x)); 687 688 /* error in velocity */ 689 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0)); 690 PetscCall(VecNorm(y0, NORM_2, &val)); 691 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error u = %g\n", (double)(PetscRealPart(val / scale)))); 692 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0)); 693 694 /* error in pressure */ 695 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1)); 696 PetscCall(VecNorm(y1, NORM_2, &val)); 697 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error p = %g\n", (double)(PetscRealPart(val / scale)))); 698 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1)); 699 700 /* total error */ 701 PetscCall(VecNorm(s->y, NORM_2, &val)); 702 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error [u,p] = %g\n", (double)PetscRealPart((val / scale)))); 703 PetscFunctionReturn(0); 704 } 705 706 int main(int argc, char **argv) 707 { 708 Stokes s; 709 KSP ksp; 710 711 PetscFunctionBeginUser; 712 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 713 s.nx = 4; 714 s.ny = 6; 715 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL)); 716 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ny", &s.ny, NULL)); 717 s.hx = 2.0 / s.nx; 718 s.hy = 1.0 / s.ny; 719 s.matsymmetric = PETSC_FALSE; 720 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL)); 721 s.userPC = s.userKSP = PETSC_FALSE; 722 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC)); 723 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP)); 724 725 PetscCall(StokesSetupMatrix(&s)); 726 PetscCall(StokesSetupIndexSets(&s)); 727 PetscCall(StokesSetupVectors(&s)); 728 729 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 730 PetscCall(KSPSetOperators(ksp, s.A, s.A)); 731 PetscCall(KSPSetFromOptions(ksp)); 732 PetscCall(StokesSetupPC(&s, ksp)); 733 PetscCall(KSPSolve(ksp, s.b, s.x)); 734 735 /* don't trust, verify! */ 736 PetscCall(StokesCalcResidual(&s)); 737 PetscCall(StokesCalcError(&s)); 738 PetscCall(StokesWriteSolution(&s)); 739 740 PetscCall(KSPDestroy(&ksp)); 741 PetscCall(MatDestroy(&s.subA[0])); 742 PetscCall(MatDestroy(&s.subA[1])); 743 PetscCall(MatDestroy(&s.subA[2])); 744 PetscCall(MatDestroy(&s.subA[3])); 745 PetscCall(MatDestroy(&s.A)); 746 PetscCall(VecDestroy(&s.x)); 747 PetscCall(VecDestroy(&s.b)); 748 PetscCall(VecDestroy(&s.y)); 749 PetscCall(MatDestroy(&s.myS)); 750 PetscCall(PetscFinalize()); 751 return 0; 752 } 753 754 /*TEST 755 756 test: 757 nsize: 2 758 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 759 760 test: 761 suffix: 2 762 nsize: 2 763 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 764 765 test: 766 suffix: 3 767 nsize: 2 768 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 769 770 test: 771 suffix: 4 772 nsize: 2 773 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 774 775 test: 776 suffix: 4_pcksp 777 nsize: 2 778 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 779 780 test: 781 suffix: 5 782 nsize: 2 783 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 784 785 test: 786 suffix: 6 787 nsize: 2 788 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 789 790 test: 791 suffix: 7 792 nsize: 2 793 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 794 795 TEST*/ 796