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