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