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) { 242 val = -1.0*val; 243 } 244 PetscCall(VecSetValue(b1, row, val, INSERT_VALUES)); 245 } 246 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); 247 PetscFunctionReturn(0); 248 } 249 250 PetscErrorCode StokesSetupMatBlock00(Stokes *s) 251 { 252 PetscInt row, start, end, sz, i, j; 253 PetscInt cols[5]; 254 PetscScalar vals[5]; 255 256 PetscFunctionBeginUser; 257 /* A[0] is 2N-by-2N */ 258 PetscCall(MatCreate(PETSC_COMM_WORLD,&s->subA[0])); 259 PetscCall(MatSetOptionsPrefix(s->subA[0],"a00_")); 260 PetscCall(MatSetSizes(s->subA[0],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,2*s->nx*s->ny)); 261 PetscCall(MatSetType(s->subA[0],MATMPIAIJ)); 262 PetscCall(MatMPIAIJSetPreallocation(s->subA[0],5,NULL,5,NULL)); 263 PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end)); 264 265 for (row = start; row < end; row++) { 266 PetscCall(StokesGetPosition(s, row, &i, &j)); 267 /* first part: rows 0 to (nx*ny-1) */ 268 PetscCall(StokesStencilLaplacian(s, i, j, &sz, cols, vals)); 269 /* second part: rows (nx*ny) to (2*nx*ny-1) */ 270 if (row >= s->nx*s->ny) { 271 for (i = 0; i < sz; i++) cols[i] += s->nx*s->ny; 272 } 273 for (i = 0; i < sz; i++) vals[i] = -1.0*vals[i]; /* dynamic viscosity coef mu=-1 */ 274 PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES)); 275 } 276 PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY)); 277 PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY)); 278 PetscFunctionReturn(0); 279 } 280 281 PetscErrorCode StokesSetupMatBlock01(Stokes *s) 282 { 283 PetscInt row, start, end, sz, i, j; 284 PetscInt cols[5]; 285 PetscScalar vals[5]; 286 287 PetscFunctionBeginUser; 288 /* A[1] is 2N-by-N */ 289 PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1])); 290 PetscCall(MatSetOptionsPrefix(s->subA[1],"a01_")); 291 PetscCall(MatSetSizes(s->subA[1],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,s->nx*s->ny)); 292 PetscCall(MatSetType(s->subA[1],MATMPIAIJ)); 293 PetscCall(MatMPIAIJSetPreallocation(s->subA[1],5,NULL,5,NULL)); 294 PetscCall(MatGetOwnershipRange(s->subA[1],&start,&end)); 295 296 PetscCall(MatSetOption(s->subA[1],MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 297 298 for (row = start; row < end; row++) { 299 PetscCall(StokesGetPosition(s, row, &i, &j)); 300 /* first part: rows 0 to (nx*ny-1) */ 301 if (row < s->nx*s->ny) { 302 PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals)); 303 } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */ 304 PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals)); 305 } 306 PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES)); 307 } 308 PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY)); 309 PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY)); 310 PetscFunctionReturn(0); 311 } 312 313 PetscErrorCode StokesSetupMatBlock10(Stokes *s) 314 { 315 PetscFunctionBeginUser; 316 /* A[2] is minus transpose of A[1] */ 317 PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2])); 318 if (!s->matsymmetric) { 319 PetscCall(MatScale(s->subA[2], -1.0)); 320 } 321 PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_")); 322 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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; vals[0]=-(ae+awb+asb+an); 392 cols[1]=e; vals[1]=ae; 393 cols[2]=n; vals[2]=an; 394 } else if (i==0 && j==s->ny-1) { /* north-west corner */ 395 *sz =3; 396 cols[0]=s2; vals[0]=as; 397 cols[1]=p; vals[1]=-(ae+awb+as+anb); 398 cols[2]=e; vals[2]=ae; 399 } else if (i==s->nx-1 && j==0) { /* south-east corner */ 400 *sz =3; 401 cols[0]=w; vals[0]=aw; 402 cols[1]=p; vals[1]=-(aeb+aw+asb+an); 403 cols[2]=n; vals[2]=an; 404 } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */ 405 *sz =3; 406 cols[0]=s2; vals[0]=as; 407 cols[1]=w; vals[1]=aw; 408 cols[2]=p; vals[2]=-(aeb+aw+as+anb); 409 } else if (i==0) { /* west boundary */ 410 *sz =4; 411 cols[0]=s2; vals[0]=as; 412 cols[1]=p; vals[1]=-(ae+awb+as+an); 413 cols[2]=e; vals[2]=ae; 414 cols[3]=n; vals[3]=an; 415 } else if (i==s->nx-1) { /* east boundary */ 416 *sz =4; 417 cols[0]=s2; vals[0]=as; 418 cols[1]=w; vals[1]=aw; 419 cols[2]=p; vals[2]=-(aeb+aw+as+an); 420 cols[3]=n; vals[3]=an; 421 } else if (j==0) { /* south boundary */ 422 *sz =4; 423 cols[0]=w; vals[0]=aw; 424 cols[1]=p; vals[1]=-(ae+aw+asb+an); 425 cols[2]=e; vals[2]=ae; 426 cols[3]=n; vals[3]=an; 427 } else if (j==s->ny-1) { /* north boundary */ 428 *sz =4; 429 cols[0]=s2; vals[0]=as; 430 cols[1]=w; vals[1]=aw; 431 cols[2]=p; vals[2]=-(ae+aw+as+anb); 432 cols[3]=e; vals[3]=ae; 433 } else { /* interior */ 434 *sz =5; 435 cols[0]=s2; vals[0]=as; 436 cols[1]=w; vals[1]=aw; 437 cols[2]=p; vals[2]=-(ae+aw+as+an); 438 cols[3]=e; vals[3]=ae; 439 cols[4]=n; vals[4]=an; 440 } 441 PetscFunctionReturn(0); 442 } 443 444 PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 445 { 446 PetscInt p =j*s->nx+i, w=p-1, e=p+1; 447 PetscScalar ae= s->hy/2, aeb=s->hy; 448 PetscScalar aw=-s->hy/2, awb=0; 449 450 PetscFunctionBeginUser; 451 if (i==0 && j==0) { /* south-west corner */ 452 *sz =2; 453 cols[0]=p; vals[0]=-(ae+awb); 454 cols[1]=e; vals[1]=ae; 455 } else if (i==0 && j==s->ny-1) { /* north-west corner */ 456 *sz =2; 457 cols[0]=p; vals[0]=-(ae+awb); 458 cols[1]=e; vals[1]=ae; 459 } else if (i==s->nx-1 && j==0) { /* south-east corner */ 460 *sz =2; 461 cols[0]=w; vals[0]=aw; 462 cols[1]=p; vals[1]=-(aeb+aw); 463 } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */ 464 *sz =2; 465 cols[0]=w; vals[0]=aw; 466 cols[1]=p; vals[1]=-(aeb+aw); 467 } else if (i==0) { /* west boundary */ 468 *sz =2; 469 cols[0]=p; vals[0]=-(ae+awb); 470 cols[1]=e; vals[1]=ae; 471 } else if (i==s->nx-1) { /* east boundary */ 472 *sz =2; 473 cols[0]=w; vals[0]=aw; 474 cols[1]=p; vals[1]=-(aeb+aw); 475 } else if (j==0) { /* south boundary */ 476 *sz =3; 477 cols[0]=w; vals[0]=aw; 478 cols[1]=p; vals[1]=-(ae+aw); 479 cols[2]=e; vals[2]=ae; 480 } else if (j==s->ny-1) { /* north boundary */ 481 *sz =3; 482 cols[0]=w; vals[0]=aw; 483 cols[1]=p; vals[1]=-(ae+aw); 484 cols[2]=e; vals[2]=ae; 485 } else { /* interior */ 486 *sz =3; 487 cols[0]=w; vals[0]=aw; 488 cols[1]=p; vals[1]=-(ae+aw); 489 cols[2]=e; vals[2]=ae; 490 } 491 PetscFunctionReturn(0); 492 } 493 494 PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 495 { 496 PetscInt p =j*s->nx+i, s2=p-s->nx, n=p+s->nx; 497 PetscScalar as=-s->hx/2, asb=0; 498 PetscScalar an= s->hx/2, anb=0; 499 500 PetscFunctionBeginUser; 501 if (i==0 && j==0) { /* south-west corner */ 502 *sz =2; 503 cols[0]=p; vals[0]=-(asb+an); 504 cols[1]=n; vals[1]=an; 505 } else if (i==0 && j==s->ny-1) { /* north-west corner */ 506 *sz =2; 507 cols[0]=s2; vals[0]=as; 508 cols[1]=p; vals[1]=-(as+anb); 509 } else if (i==s->nx-1 && j==0) { /* south-east corner */ 510 *sz =2; 511 cols[0]=p; vals[0]=-(asb+an); 512 cols[1]=n; vals[1]=an; 513 } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */ 514 *sz =2; 515 cols[0]=s2; vals[0]=as; 516 cols[1]=p; vals[1]=-(as+anb); 517 } else if (i==0) { /* west boundary */ 518 *sz =3; 519 cols[0]=s2; vals[0]=as; 520 cols[1]=p; vals[1]=-(as+an); 521 cols[2]=n; vals[2]=an; 522 } else if (i==s->nx-1) { /* east boundary */ 523 *sz =3; 524 cols[0]=s2; vals[0]=as; 525 cols[1]=p; vals[1]=-(as+an); 526 cols[2]=n; vals[2]=an; 527 } else if (j==0) { /* south boundary */ 528 *sz =2; 529 cols[0]=p; vals[0]=-(asb+an); 530 cols[1]=n; vals[1]=an; 531 } else if (j==s->ny-1) { /* north boundary */ 532 *sz =2; 533 cols[0]=s2; vals[0]=as; 534 cols[1]=p; vals[1]=-(as+anb); 535 } else { /* interior */ 536 *sz =3; 537 cols[0]=s2; vals[0]=as; 538 cols[1]=p; vals[1]=-(as+an); 539 cols[2]=n; vals[2]=an; 540 } 541 PetscFunctionReturn(0); 542 } 543 544 PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 545 { 546 PetscScalar y = j*s->hy+s->hy/2; 547 PetscScalar awb = s->hy/(s->hx/2); 548 549 PetscFunctionBeginUser; 550 if (i == 0) { /* west boundary */ 551 *val = awb*StokesExactVelocityX(y); 552 } else { 553 *val = 0.0; 554 } 555 PetscFunctionReturn(0); 556 } 557 558 PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 559 { 560 PetscFunctionBeginUser; 561 *val = 0.0; 562 PetscFunctionReturn(0); 563 } 564 565 PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 566 { 567 PetscScalar y = j*s->hy+s->hy/2; 568 PetscScalar aeb = s->hy; 569 570 PetscFunctionBeginUser; 571 if (i == 0) { /* west boundary */ 572 *val = aeb*StokesExactVelocityX(y); 573 } else { 574 *val = 0.0; 575 } 576 PetscFunctionReturn(0); 577 } 578 579 PetscErrorCode StokesCalcResidual(Stokes *s) 580 { 581 PetscReal val; 582 Vec b0, b1; 583 584 PetscFunctionBeginUser; 585 /* residual Ax-b (*warning* overwrites b) */ 586 PetscCall(VecScale(s->b, -1.0)); 587 PetscCall(MatMultAdd(s->A, s->x, s->b, s->b)); 588 589 /* residual velocity */ 590 PetscCall(VecGetSubVector(s->b, s->isg[0], &b0)); 591 PetscCall(VecNorm(b0, NORM_2, &val)); 592 PetscCall(PetscPrintf(PETSC_COMM_WORLD," residual u = %g\n",(double)val)); 593 PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0)); 594 595 /* residual pressure */ 596 PetscCall(VecGetSubVector(s->b, s->isg[1], &b1)); 597 PetscCall(VecNorm(b1, NORM_2, &val)); 598 PetscCall(PetscPrintf(PETSC_COMM_WORLD," residual p = %g\n",(double)val)); 599 PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1)); 600 601 /* total residual */ 602 PetscCall(VecNorm(s->b, NORM_2, &val)); 603 PetscCall(PetscPrintf(PETSC_COMM_WORLD," residual [u,p] = %g\n", (double)val)); 604 PetscFunctionReturn(0); 605 } 606 607 PetscErrorCode StokesCalcError(Stokes *s) 608 { 609 PetscScalar scale = PetscSqrtReal((double)s->nx*s->ny); 610 PetscReal val; 611 Vec y0, y1; 612 613 PetscFunctionBeginUser; 614 /* error y-x */ 615 PetscCall(VecAXPY(s->y, -1.0, s->x)); 616 617 /* error in velocity */ 618 PetscCall(VecGetSubVector(s->y, s->isg[0], &y0)); 619 PetscCall(VecNorm(y0, NORM_2, &val)); 620 PetscCall(PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(PetscRealPart(val/scale)))); 621 PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0)); 622 623 /* error in pressure */ 624 PetscCall(VecGetSubVector(s->y, s->isg[1], &y1)); 625 PetscCall(VecNorm(y1, NORM_2, &val)); 626 PetscCall(PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(PetscRealPart(val/scale)))); 627 PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1)); 628 629 /* total error */ 630 PetscCall(VecNorm(s->y, NORM_2, &val)); 631 PetscCall(PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)PetscRealPart((val/scale)))); 632 PetscFunctionReturn(0); 633 } 634 635 int main(int argc, char **argv) 636 { 637 Stokes s; 638 KSP ksp; 639 640 PetscCall(PetscInitialize(&argc, &argv, NULL,help)); 641 s.nx = 4; 642 s.ny = 6; 643 PetscCall(PetscOptionsGetInt(NULL,NULL, "-nx", &s.nx, NULL)); 644 PetscCall(PetscOptionsGetInt(NULL,NULL, "-ny", &s.ny, NULL)); 645 s.hx = 2.0/s.nx; 646 s.hy = 1.0/s.ny; 647 s.matsymmetric = PETSC_FALSE; 648 PetscCall(PetscOptionsGetBool(NULL,NULL, "-mat_set_symmetric", &s.matsymmetric,NULL)); 649 s.userPC = s.userKSP = PETSC_FALSE; 650 PetscCall(PetscOptionsHasName(NULL,NULL, "-user_pc", &s.userPC)); 651 PetscCall(PetscOptionsHasName(NULL,NULL, "-user_ksp", &s.userKSP)); 652 653 PetscCall(StokesSetupMatrix(&s)); 654 PetscCall(StokesSetupIndexSets(&s)); 655 PetscCall(StokesSetupVectors(&s)); 656 657 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 658 PetscCall(KSPSetOperators(ksp, s.A, s.A)); 659 PetscCall(KSPSetFromOptions(ksp)); 660 PetscCall(StokesSetupPC(&s, ksp)); 661 PetscCall(KSPSolve(ksp, s.b, s.x)); 662 663 /* don't trust, verify! */ 664 PetscCall(StokesCalcResidual(&s)); 665 PetscCall(StokesCalcError(&s)); 666 PetscCall(StokesWriteSolution(&s)); 667 668 PetscCall(KSPDestroy(&ksp)); 669 PetscCall(MatDestroy(&s.subA[0])); 670 PetscCall(MatDestroy(&s.subA[1])); 671 PetscCall(MatDestroy(&s.subA[2])); 672 PetscCall(MatDestroy(&s.subA[3])); 673 PetscCall(MatDestroy(&s.A)); 674 PetscCall(VecDestroy(&s.x)); 675 PetscCall(VecDestroy(&s.b)); 676 PetscCall(VecDestroy(&s.y)); 677 PetscCall(MatDestroy(&s.myS)); 678 PetscCall(PetscFinalize()); 679 return 0; 680 } 681 682 /*TEST 683 684 test: 685 nsize: 2 686 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 687 688 test: 689 suffix: 2 690 nsize: 2 691 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 692 693 test: 694 suffix: 3 695 nsize: 2 696 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 697 698 test: 699 suffix: 4 700 nsize: 2 701 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 702 703 test: 704 suffix: 4_pcksp 705 nsize: 2 706 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 707 708 test: 709 suffix: 5 710 nsize: 2 711 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 712 713 test: 714 suffix: 6 715 nsize: 2 716 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 717 718 test: 719 suffix: 7 720 nsize: 2 721 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 722 723 TEST*/ 724