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