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 /*---------------------------------------------------------------------------- */ 7 /* author : Christiaan M. Klaij */ 8 /*---------------------------------------------------------------------------- */ 9 /* */ 10 /* Poiseuille flow problem. */ 11 /* */ 12 /* Viscous, laminar flow in a 2D channel with parabolic velocity */ 13 /* profile and linear pressure drop, exact solution of the 2D Stokes */ 14 /* equations. */ 15 /* */ 16 /* Discretized with the cell-centered finite-volume method on a */ 17 /* Cartesian grid with co-located variables. Variables ordered as */ 18 /* [u1...uN v1...vN p1...pN]^T. Matrix [A00 A01; A10, A11] solved with */ 19 /* PCFIELDSPLIT. Lower factorization is used to mimic the Semi-Implicit */ 20 /* Method for Pressure Linked Equations (SIMPLE) used as preconditioner */ 21 /* instead of solver. */ 22 /* */ 23 /* Disclaimer: does not contain the pressure-weighed interpolation */ 24 /* method needed to suppress spurious pressure modes in real-life */ 25 /* problems. */ 26 /* */ 27 /* Usage: */ 28 /* */ 29 /* 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 */ 30 /* */ 31 /* Runs with PCFIELDSPLIT on 32x48 grid, no PC for the Schur */ 32 /* complement because A11 is zero. FGMRES is needed because */ 33 /* PCFIELDSPLIT is a variable preconditioner. */ 34 /* */ 35 /* 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 */ 36 /* */ 37 /* Same as above but with user defined PC for the true Schur */ 38 /* complement. PC based on the SIMPLE-type approximation (inverse of */ 39 /* A00 approximated by inverse of its diagonal). */ 40 /* */ 41 /* 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 */ 42 /* */ 43 /* Replace the true Schur complement with a user defined Schur */ 44 /* complement based on the SIMPLE-type approximation. Same matrix is */ 45 /* used as PC. */ 46 /* */ 47 /* 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 */ 48 /* */ 49 /* Out-of-the-box SIMPLE-type preconditioning. The major advantage */ 50 /* is that the user neither needs to provide the approximation of */ 51 /* the Schur complement, nor the corresponding preconditioner. */ 52 /* */ 53 /*---------------------------------------------------------------------------- */ 54 55 #include <petscksp.h> 56 57 typedef struct { 58 PetscBool userPC, userKSP, matsymmetric; /* user defined preconditioner and matrix for the Schur complement */ 59 PetscInt nx, ny; /* nb of cells in x- and y-direction */ 60 PetscReal hx, hy; /* mesh size in x- and y-direction */ 61 Mat A; /* block matrix */ 62 Mat subA[4]; /* the four blocks */ 63 Mat myS; /* the approximation of the Schur complement */ 64 Vec x, b, y; /* solution, rhs and temporary vector */ 65 IS isg[2]; /* index sets of split "0" and "1" */ 66 } Stokes; 67 68 PetscErrorCode StokesSetupMatBlock00(Stokes*); /* setup the block Q */ 69 PetscErrorCode StokesSetupMatBlock01(Stokes*); /* setup the block G */ 70 PetscErrorCode StokesSetupMatBlock10(Stokes*); /* setup the block D (equal to the transpose of G) */ 71 PetscErrorCode StokesSetupMatBlock11(Stokes*); /* setup the block C (equal to zero) */ 72 73 PetscErrorCode StokesGetPosition(Stokes*, PetscInt, PetscInt*, PetscInt*); /* row number j*nx+i corresponds to position (i,j) in grid */ 74 75 PetscErrorCode StokesStencilLaplacian(Stokes*, PetscInt, PetscInt, PetscInt*, PetscInt*, PetscScalar*); /* stencil of the Laplacian operator */ 76 PetscErrorCode StokesStencilGradientX(Stokes*, PetscInt, PetscInt, PetscInt*, PetscInt*, PetscScalar*); /* stencil of the Gradient operator (x-component) */ 77 PetscErrorCode StokesStencilGradientY(Stokes*, PetscInt, PetscInt, PetscInt*, PetscInt*, PetscScalar*); /* stencil of the Gradient operator (y-component) */ 78 79 PetscErrorCode StokesRhs(Stokes*); /* rhs vector */ 80 PetscErrorCode StokesRhsMomX(Stokes*, PetscInt, PetscInt, PetscScalar*); /* right hand side of velocity (x-component) */ 81 PetscErrorCode StokesRhsMomY(Stokes*, PetscInt, PetscInt, PetscScalar*); /* right hand side of velocity (y-component) */ 82 PetscErrorCode StokesRhsMass(Stokes*, PetscInt, PetscInt, PetscScalar*); /* right hand side of pressure */ 83 84 PetscErrorCode StokesSetupApproxSchur(Stokes*); /* approximation of the Schur complement */ 85 86 PetscErrorCode StokesExactSolution(Stokes*); /* exact solution vector */ 87 PetscErrorCode StokesWriteSolution(Stokes*); /* write solution to file */ 88 89 /* exact solution for the velocity (x-component, y-component is zero) */ 90 PetscScalar StokesExactVelocityX(const PetscScalar y) 91 { 92 return 4.0*y*(1.0-y); 93 } 94 95 /* exact solution for the pressure */ 96 PetscScalar StokesExactPressure(const PetscScalar x) 97 { 98 return 8.0*(2.0-x); 99 } 100 101 PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp) 102 { 103 KSP *subksp; 104 PC pc; 105 PetscInt n = 1; 106 PetscErrorCode ierr; 107 108 PetscFunctionBeginUser; 109 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 110 ierr = PCFieldSplitSetIS(pc, "0", s->isg[0]);CHKERRQ(ierr); 111 ierr = PCFieldSplitSetIS(pc, "1", s->isg[1]);CHKERRQ(ierr); 112 if (s->userPC) { 113 ierr = PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS);CHKERRQ(ierr); 114 } 115 if (s->userKSP) { 116 ierr = PCSetUp(pc);CHKERRQ(ierr); 117 ierr = PCFieldSplitGetSubKSP(pc, &n, &subksp);CHKERRQ(ierr); 118 ierr = KSPSetOperators(subksp[1], s->myS, s->myS);CHKERRQ(ierr); 119 ierr = PetscFree(subksp);CHKERRQ(ierr); 120 } 121 PetscFunctionReturn(0); 122 } 123 124 PetscErrorCode StokesWriteSolution(Stokes *s) 125 { 126 PetscMPIInt size; 127 PetscInt n,i,j; 128 const PetscScalar *array; 129 PetscErrorCode ierr; 130 131 PetscFunctionBeginUser; 132 /* write data (*warning* only works sequential) */ 133 ierr = MPI_Comm_size(MPI_COMM_WORLD,&size);CHKERRMPI(ierr); 134 /*ierr = PetscPrintf(PETSC_COMM_WORLD," number of processors = %D\n",size);CHKERRQ(ierr);*/ 135 if (size == 1) { 136 PetscViewer viewer; 137 ierr = VecGetArrayRead(s->x, &array);CHKERRQ(ierr); 138 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer);CHKERRQ(ierr); 139 ierr = PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n");CHKERRQ(ierr); 140 for (j = 0; j < s->ny; j++) { 141 for (i = 0; i < s->nx; i++) { 142 n = j*s->nx+i; 143 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); 144 } 145 } 146 ierr = VecRestoreArrayRead(s->x, &array);CHKERRQ(ierr); 147 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 148 } 149 PetscFunctionReturn(0); 150 } 151 152 PetscErrorCode StokesSetupIndexSets(Stokes *s) 153 { 154 PetscErrorCode ierr; 155 156 PetscFunctionBeginUser; 157 /* the two index sets */ 158 ierr = MatNestGetISs(s->A, s->isg, NULL);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 PetscErrorCode StokesSetupVectors(Stokes *s) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBeginUser; 167 /* solution vector x */ 168 ierr = VecCreate(PETSC_COMM_WORLD, &s->x);CHKERRQ(ierr); 169 ierr = VecSetSizes(s->x, PETSC_DECIDE, 3*s->nx*s->ny);CHKERRQ(ierr); 170 ierr = VecSetType(s->x, VECMPI);CHKERRQ(ierr); 171 /* ierr = VecSetRandom(s->x, NULL);CHKERRQ(ierr); */ 172 /* ierr = VecView(s->x, (PetscViewer) PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); */ 173 174 /* exact solution y */ 175 ierr = VecDuplicate(s->x, &s->y);CHKERRQ(ierr); 176 ierr = StokesExactSolution(s);CHKERRQ(ierr); 177 /* ierr = VecView(s->y, (PetscViewer) PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); */ 178 179 /* rhs vector b */ 180 ierr = VecDuplicate(s->x, &s->b);CHKERRQ(ierr); 181 ierr = StokesRhs(s);CHKERRQ(ierr); 182 /*ierr = VecView(s->b, (PetscViewer) PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);*/ 183 PetscFunctionReturn(0); 184 } 185 186 PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j) 187 { 188 PetscInt n; 189 190 PetscFunctionBeginUser; 191 /* cell number n=j*nx+i has position (i,j) in grid */ 192 n = row%(s->nx*s->ny); 193 *i = n%s->nx; 194 *j = (n-(*i))/s->nx; 195 PetscFunctionReturn(0); 196 } 197 198 PetscErrorCode StokesExactSolution(Stokes *s) 199 { 200 PetscInt row, start, end, i, j; 201 PetscScalar val; 202 Vec y0,y1; 203 PetscErrorCode ierr; 204 205 PetscFunctionBeginUser; 206 /* velocity part */ 207 ierr = VecGetSubVector(s->y, s->isg[0], &y0);CHKERRQ(ierr); 208 ierr = VecGetOwnershipRange(y0, &start, &end);CHKERRQ(ierr); 209 for (row = start; row < end; row++) { 210 ierr = StokesGetPosition(s, row,&i,&j);CHKERRQ(ierr); 211 if (row < s->nx*s->ny) { 212 val = StokesExactVelocityX(j*s->hy+s->hy/2); 213 } else { 214 val = 0; 215 } 216 ierr = VecSetValue(y0, row, val, INSERT_VALUES);CHKERRQ(ierr); 217 } 218 ierr = VecRestoreSubVector(s->y, s->isg[0], &y0);CHKERRQ(ierr); 219 220 /* pressure part */ 221 ierr = VecGetSubVector(s->y, s->isg[1], &y1);CHKERRQ(ierr); 222 ierr = VecGetOwnershipRange(y1, &start, &end);CHKERRQ(ierr); 223 for (row = start; row < end; row++) { 224 ierr = StokesGetPosition(s, row, &i, &j);CHKERRQ(ierr); 225 val = StokesExactPressure(i*s->hx+s->hx/2); 226 ierr = VecSetValue(y1, row, val, INSERT_VALUES);CHKERRQ(ierr); 227 } 228 ierr = VecRestoreSubVector(s->y, s->isg[1], &y1);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 PetscErrorCode StokesRhs(Stokes *s) 233 { 234 PetscInt row, start, end, i, j; 235 PetscScalar val; 236 Vec b0,b1; 237 PetscErrorCode ierr; 238 239 PetscFunctionBeginUser; 240 /* velocity part */ 241 ierr = VecGetSubVector(s->b, s->isg[0], &b0);CHKERRQ(ierr); 242 ierr = VecGetOwnershipRange(b0, &start, &end);CHKERRQ(ierr); 243 for (row = start; row < end; row++) { 244 ierr = StokesGetPosition(s, row, &i, &j);CHKERRQ(ierr); 245 if (row < s->nx*s->ny) { 246 ierr = StokesRhsMomX(s, i, j, &val);CHKERRQ(ierr); 247 } else { 248 ierr = StokesRhsMomY(s, i, j, &val);CHKERRQ(ierr); 249 } 250 ierr = VecSetValue(b0, row, val, INSERT_VALUES);CHKERRQ(ierr); 251 } 252 ierr = VecRestoreSubVector(s->b, s->isg[0], &b0);CHKERRQ(ierr); 253 254 /* pressure part */ 255 ierr = VecGetSubVector(s->b, s->isg[1], &b1);CHKERRQ(ierr); 256 ierr = VecGetOwnershipRange(b1, &start, &end);CHKERRQ(ierr); 257 for (row = start; row < end; row++) { 258 ierr = StokesGetPosition(s, row, &i, &j);CHKERRQ(ierr); 259 ierr = StokesRhsMass(s, i, j, &val);CHKERRQ(ierr); 260 if (s->matsymmetric) { 261 val = -1.0*val; 262 } 263 ierr = VecSetValue(b1, row, val, INSERT_VALUES);CHKERRQ(ierr); 264 } 265 ierr = VecRestoreSubVector(s->b, s->isg[1], &b1);CHKERRQ(ierr); 266 PetscFunctionReturn(0); 267 } 268 269 PetscErrorCode StokesSetupMatBlock00(Stokes *s) 270 { 271 PetscInt row, start, end, sz, i, j; 272 PetscInt cols[5]; 273 PetscScalar vals[5]; 274 PetscErrorCode ierr; 275 276 PetscFunctionBeginUser; 277 /* A[0] is 2N-by-2N */ 278 ierr = MatCreate(PETSC_COMM_WORLD,&s->subA[0]);CHKERRQ(ierr); 279 ierr = MatSetOptionsPrefix(s->subA[0],"a00_");CHKERRQ(ierr); 280 ierr = MatSetSizes(s->subA[0],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,2*s->nx*s->ny);CHKERRQ(ierr); 281 ierr = MatSetType(s->subA[0],MATMPIAIJ);CHKERRQ(ierr); 282 ierr = MatMPIAIJSetPreallocation(s->subA[0],5,NULL,5,NULL);CHKERRQ(ierr); 283 ierr = MatGetOwnershipRange(s->subA[0], &start, &end);CHKERRQ(ierr); 284 285 for (row = start; row < end; row++) { 286 ierr = StokesGetPosition(s, row, &i, &j);CHKERRQ(ierr); 287 /* first part: rows 0 to (nx*ny-1) */ 288 ierr = StokesStencilLaplacian(s, i, j, &sz, cols, vals);CHKERRQ(ierr); 289 /* second part: rows (nx*ny) to (2*nx*ny-1) */ 290 if (row >= s->nx*s->ny) { 291 for (i = 0; i < sz; i++) cols[i] += s->nx*s->ny; 292 } 293 for (i = 0; i < sz; i++) vals[i] = -1.0*vals[i]; /* dynamic viscosity coef mu=-1 */ 294 ierr = MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES);CHKERRQ(ierr); 295 } 296 ierr = MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 297 ierr = MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 298 PetscFunctionReturn(0); 299 } 300 301 PetscErrorCode StokesSetupMatBlock01(Stokes *s) 302 { 303 PetscInt row, start, end, sz, i, j; 304 PetscInt cols[5]; 305 PetscScalar vals[5]; 306 PetscErrorCode ierr; 307 308 PetscFunctionBeginUser; 309 /* A[1] is 2N-by-N */ 310 ierr = MatCreate(PETSC_COMM_WORLD, &s->subA[1]);CHKERRQ(ierr); 311 ierr = MatSetOptionsPrefix(s->subA[1],"a01_");CHKERRQ(ierr); 312 ierr = MatSetSizes(s->subA[1],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,s->nx*s->ny);CHKERRQ(ierr); 313 ierr = MatSetType(s->subA[1],MATMPIAIJ);CHKERRQ(ierr); 314 ierr = MatMPIAIJSetPreallocation(s->subA[1],5,NULL,5,NULL);CHKERRQ(ierr); 315 ierr = MatGetOwnershipRange(s->subA[1],&start,&end);CHKERRQ(ierr); 316 317 ierr = MatSetOption(s->subA[1],MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 318 319 for (row = start; row < end; row++) { 320 ierr = StokesGetPosition(s, row, &i, &j);CHKERRQ(ierr); 321 /* first part: rows 0 to (nx*ny-1) */ 322 if (row < s->nx*s->ny) { 323 ierr = StokesStencilGradientX(s, i, j, &sz, cols, vals);CHKERRQ(ierr); 324 } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */ 325 ierr = StokesStencilGradientY(s, i, j, &sz, cols, vals);CHKERRQ(ierr); 326 } 327 ierr = MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES);CHKERRQ(ierr); 328 } 329 ierr = MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 330 ierr = MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 331 PetscFunctionReturn(0); 332 } 333 334 PetscErrorCode StokesSetupMatBlock10(Stokes *s) 335 { 336 PetscErrorCode ierr; 337 338 PetscFunctionBeginUser; 339 /* A[2] is minus transpose of A[1] */ 340 ierr = MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]);CHKERRQ(ierr); 341 if (!s->matsymmetric) { 342 ierr = MatScale(s->subA[2], -1.0);CHKERRQ(ierr); 343 } 344 ierr = MatSetOptionsPrefix(s->subA[2], "a10_");CHKERRQ(ierr); 345 PetscFunctionReturn(0); 346 } 347 348 PetscErrorCode StokesSetupMatBlock11(Stokes *s) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBeginUser; 353 /* A[3] is N-by-N null matrix */ 354 ierr = MatCreate(PETSC_COMM_WORLD, &s->subA[3]);CHKERRQ(ierr); 355 ierr = MatSetOptionsPrefix(s->subA[3], "a11_");CHKERRQ(ierr); 356 ierr = MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx*s->ny, s->nx*s->ny);CHKERRQ(ierr); 357 ierr = MatSetType(s->subA[3], MATMPIAIJ);CHKERRQ(ierr); 358 ierr = MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL);CHKERRQ(ierr); 359 ierr = MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 360 ierr = MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 PetscErrorCode StokesSetupApproxSchur(Stokes *s) 365 { 366 Vec diag; 367 PetscErrorCode ierr; 368 369 PetscFunctionBeginUser; 370 /* Schur complement approximation: myS = A11 - A10 diag(A00)^(-1) A01 */ 371 /* note: A11 is zero */ 372 /* note: in real life this matrix would be build directly, */ 373 /* i.e. without MatMatMult */ 374 375 /* inverse of diagonal of A00 */ 376 ierr = VecCreate(PETSC_COMM_WORLD,&diag);CHKERRQ(ierr); 377 ierr = VecSetSizes(diag,PETSC_DECIDE,2*s->nx*s->ny);CHKERRQ(ierr); 378 ierr = VecSetType(diag,VECMPI);CHKERRQ(ierr); 379 ierr = MatGetDiagonal(s->subA[0],diag);CHKERRQ(ierr); 380 ierr = VecReciprocal(diag);CHKERRQ(ierr); 381 382 /* compute: - A10 diag(A00)^(-1) A01 */ 383 ierr = MatDiagonalScale(s->subA[1],diag,NULL);CHKERRQ(ierr); /* (*warning* overwrites subA[1]) */ 384 ierr = MatMatMult(s->subA[2],s->subA[1],MAT_INITIAL_MATRIX,PETSC_DEFAULT,&s->myS);CHKERRQ(ierr); 385 ierr = MatScale(s->myS,-1.0);CHKERRQ(ierr); 386 387 /* restore A10 */ 388 ierr = MatGetDiagonal(s->subA[0],diag);CHKERRQ(ierr); 389 ierr = MatDiagonalScale(s->subA[1],diag,NULL);CHKERRQ(ierr); 390 ierr = VecDestroy(&diag);CHKERRQ(ierr); 391 PetscFunctionReturn(0); 392 } 393 394 PetscErrorCode StokesSetupMatrix(Stokes *s) 395 { 396 PetscErrorCode ierr; 397 398 PetscFunctionBeginUser; 399 ierr = StokesSetupMatBlock00(s);CHKERRQ(ierr); 400 ierr = StokesSetupMatBlock01(s);CHKERRQ(ierr); 401 ierr = StokesSetupMatBlock10(s);CHKERRQ(ierr); 402 ierr = StokesSetupMatBlock11(s);CHKERRQ(ierr); 403 ierr = MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A);CHKERRQ(ierr); 404 ierr = StokesSetupApproxSchur(s);CHKERRQ(ierr); 405 PetscFunctionReturn(0); 406 } 407 408 PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 409 { 410 PetscInt p =j*s->nx+i, w=p-1, e=p+1, s2=p-s->nx, n=p+s->nx; 411 PetscScalar ae=s->hy/s->hx, aeb=0; 412 PetscScalar aw=s->hy/s->hx, awb=s->hy/(s->hx/2); 413 PetscScalar as=s->hx/s->hy, asb=s->hx/(s->hy/2); 414 PetscScalar an=s->hx/s->hy, anb=s->hx/(s->hy/2); 415 416 PetscFunctionBeginUser; 417 if (i==0 && j==0) { /* south-west corner */ 418 *sz =3; 419 cols[0]=p; vals[0]=-(ae+awb+asb+an); 420 cols[1]=e; vals[1]=ae; 421 cols[2]=n; vals[2]=an; 422 } else if (i==0 && j==s->ny-1) { /* north-west corner */ 423 *sz =3; 424 cols[0]=s2; vals[0]=as; 425 cols[1]=p; vals[1]=-(ae+awb+as+anb); 426 cols[2]=e; vals[2]=ae; 427 } else if (i==s->nx-1 && j==0) { /* south-east corner */ 428 *sz =3; 429 cols[0]=w; vals[0]=aw; 430 cols[1]=p; vals[1]=-(aeb+aw+asb+an); 431 cols[2]=n; vals[2]=an; 432 } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */ 433 *sz =3; 434 cols[0]=s2; vals[0]=as; 435 cols[1]=w; vals[1]=aw; 436 cols[2]=p; vals[2]=-(aeb+aw+as+anb); 437 } else if (i==0) { /* west boundary */ 438 *sz =4; 439 cols[0]=s2; vals[0]=as; 440 cols[1]=p; vals[1]=-(ae+awb+as+an); 441 cols[2]=e; vals[2]=ae; 442 cols[3]=n; vals[3]=an; 443 } else if (i==s->nx-1) { /* east boundary */ 444 *sz =4; 445 cols[0]=s2; vals[0]=as; 446 cols[1]=w; vals[1]=aw; 447 cols[2]=p; vals[2]=-(aeb+aw+as+an); 448 cols[3]=n; vals[3]=an; 449 } else if (j==0) { /* south boundary */ 450 *sz =4; 451 cols[0]=w; vals[0]=aw; 452 cols[1]=p; vals[1]=-(ae+aw+asb+an); 453 cols[2]=e; vals[2]=ae; 454 cols[3]=n; vals[3]=an; 455 } else if (j==s->ny-1) { /* north boundary */ 456 *sz =4; 457 cols[0]=s2; vals[0]=as; 458 cols[1]=w; vals[1]=aw; 459 cols[2]=p; vals[2]=-(ae+aw+as+anb); 460 cols[3]=e; vals[3]=ae; 461 } else { /* interior */ 462 *sz =5; 463 cols[0]=s2; vals[0]=as; 464 cols[1]=w; vals[1]=aw; 465 cols[2]=p; vals[2]=-(ae+aw+as+an); 466 cols[3]=e; vals[3]=ae; 467 cols[4]=n; vals[4]=an; 468 } 469 PetscFunctionReturn(0); 470 } 471 472 PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 473 { 474 PetscInt p =j*s->nx+i, w=p-1, e=p+1; 475 PetscScalar ae= s->hy/2, aeb=s->hy; 476 PetscScalar aw=-s->hy/2, awb=0; 477 478 PetscFunctionBeginUser; 479 if (i==0 && j==0) { /* south-west corner */ 480 *sz =2; 481 cols[0]=p; vals[0]=-(ae+awb); 482 cols[1]=e; vals[1]=ae; 483 } else if (i==0 && j==s->ny-1) { /* north-west corner */ 484 *sz =2; 485 cols[0]=p; vals[0]=-(ae+awb); 486 cols[1]=e; vals[1]=ae; 487 } else if (i==s->nx-1 && j==0) { /* south-east corner */ 488 *sz =2; 489 cols[0]=w; vals[0]=aw; 490 cols[1]=p; vals[1]=-(aeb+aw); 491 } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */ 492 *sz =2; 493 cols[0]=w; vals[0]=aw; 494 cols[1]=p; vals[1]=-(aeb+aw); 495 } else if (i==0) { /* west boundary */ 496 *sz =2; 497 cols[0]=p; vals[0]=-(ae+awb); 498 cols[1]=e; vals[1]=ae; 499 } else if (i==s->nx-1) { /* east boundary */ 500 *sz =2; 501 cols[0]=w; vals[0]=aw; 502 cols[1]=p; vals[1]=-(aeb+aw); 503 } else if (j==0) { /* south boundary */ 504 *sz =3; 505 cols[0]=w; vals[0]=aw; 506 cols[1]=p; vals[1]=-(ae+aw); 507 cols[2]=e; vals[2]=ae; 508 } else if (j==s->ny-1) { /* north boundary */ 509 *sz =3; 510 cols[0]=w; vals[0]=aw; 511 cols[1]=p; vals[1]=-(ae+aw); 512 cols[2]=e; vals[2]=ae; 513 } else { /* interior */ 514 *sz =3; 515 cols[0]=w; vals[0]=aw; 516 cols[1]=p; vals[1]=-(ae+aw); 517 cols[2]=e; vals[2]=ae; 518 } 519 PetscFunctionReturn(0); 520 } 521 522 PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals) 523 { 524 PetscInt p =j*s->nx+i, s2=p-s->nx, n=p+s->nx; 525 PetscScalar as=-s->hx/2, asb=0; 526 PetscScalar an= s->hx/2, anb=0; 527 528 PetscFunctionBeginUser; 529 if (i==0 && j==0) { /* south-west corner */ 530 *sz =2; 531 cols[0]=p; vals[0]=-(asb+an); 532 cols[1]=n; vals[1]=an; 533 } else if (i==0 && j==s->ny-1) { /* north-west corner */ 534 *sz =2; 535 cols[0]=s2; vals[0]=as; 536 cols[1]=p; vals[1]=-(as+anb); 537 } else if (i==s->nx-1 && j==0) { /* south-east corner */ 538 *sz =2; 539 cols[0]=p; vals[0]=-(asb+an); 540 cols[1]=n; vals[1]=an; 541 } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */ 542 *sz =2; 543 cols[0]=s2; vals[0]=as; 544 cols[1]=p; vals[1]=-(as+anb); 545 } else if (i==0) { /* west boundary */ 546 *sz =3; 547 cols[0]=s2; vals[0]=as; 548 cols[1]=p; vals[1]=-(as+an); 549 cols[2]=n; vals[2]=an; 550 } else if (i==s->nx-1) { /* east boundary */ 551 *sz =3; 552 cols[0]=s2; vals[0]=as; 553 cols[1]=p; vals[1]=-(as+an); 554 cols[2]=n; vals[2]=an; 555 } else if (j==0) { /* south boundary */ 556 *sz =2; 557 cols[0]=p; vals[0]=-(asb+an); 558 cols[1]=n; vals[1]=an; 559 } else if (j==s->ny-1) { /* north boundary */ 560 *sz =2; 561 cols[0]=s2; vals[0]=as; 562 cols[1]=p; vals[1]=-(as+anb); 563 } else { /* interior */ 564 *sz =3; 565 cols[0]=s2; vals[0]=as; 566 cols[1]=p; vals[1]=-(as+an); 567 cols[2]=n; vals[2]=an; 568 } 569 PetscFunctionReturn(0); 570 } 571 572 PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 573 { 574 PetscScalar y = j*s->hy+s->hy/2; 575 PetscScalar awb = s->hy/(s->hx/2); 576 577 PetscFunctionBeginUser; 578 if (i == 0) { /* west boundary */ 579 *val = awb*StokesExactVelocityX(y); 580 } else { 581 *val = 0.0; 582 } 583 PetscFunctionReturn(0); 584 } 585 586 PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 587 { 588 PetscFunctionBeginUser; 589 *val = 0.0; 590 PetscFunctionReturn(0); 591 } 592 593 PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val) 594 { 595 PetscScalar y = j*s->hy+s->hy/2; 596 PetscScalar aeb = s->hy; 597 598 PetscFunctionBeginUser; 599 if (i == 0) { /* west boundary */ 600 *val = aeb*StokesExactVelocityX(y); 601 } else { 602 *val = 0.0; 603 } 604 PetscFunctionReturn(0); 605 } 606 607 PetscErrorCode StokesCalcResidual(Stokes *s) 608 { 609 PetscReal val; 610 Vec b0, b1; 611 PetscErrorCode ierr; 612 613 PetscFunctionBeginUser; 614 /* residual Ax-b (*warning* overwrites b) */ 615 ierr = VecScale(s->b, -1.0);CHKERRQ(ierr); 616 ierr = MatMultAdd(s->A, s->x, s->b, s->b);CHKERRQ(ierr); 617 /* ierr = VecView(s->b, (PetscViewer)PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); */ 618 619 /* residual velocity */ 620 ierr = VecGetSubVector(s->b, s->isg[0], &b0);CHKERRQ(ierr); 621 ierr = VecNorm(b0, NORM_2, &val);CHKERRQ(ierr); 622 ierr = PetscPrintf(PETSC_COMM_WORLD," residual u = %g\n",(double)val);CHKERRQ(ierr); 623 ierr = VecRestoreSubVector(s->b, s->isg[0], &b0);CHKERRQ(ierr); 624 625 /* residual pressure */ 626 ierr = VecGetSubVector(s->b, s->isg[1], &b1);CHKERRQ(ierr); 627 ierr = VecNorm(b1, NORM_2, &val);CHKERRQ(ierr); 628 ierr = PetscPrintf(PETSC_COMM_WORLD," residual p = %g\n",(double)val);CHKERRQ(ierr); 629 ierr = VecRestoreSubVector(s->b, s->isg[1], &b1);CHKERRQ(ierr); 630 631 /* total residual */ 632 ierr = VecNorm(s->b, NORM_2, &val);CHKERRQ(ierr); 633 ierr = PetscPrintf(PETSC_COMM_WORLD," residual [u,p] = %g\n", (double)val);CHKERRQ(ierr); 634 PetscFunctionReturn(0); 635 } 636 637 PetscErrorCode StokesCalcError(Stokes *s) 638 { 639 PetscScalar scale = PetscSqrtReal((double)s->nx*s->ny); 640 PetscReal val; 641 Vec y0, y1; 642 PetscErrorCode ierr; 643 644 PetscFunctionBeginUser; 645 /* error y-x */ 646 ierr = VecAXPY(s->y, -1.0, s->x);CHKERRQ(ierr); 647 /* ierr = VecView(s->y, (PetscViewer)PETSC_VIEWER_DEFAULT);CHKERRQ(ierr); */ 648 649 /* error in velocity */ 650 ierr = VecGetSubVector(s->y, s->isg[0], &y0);CHKERRQ(ierr); 651 ierr = VecNorm(y0, NORM_2, &val);CHKERRQ(ierr); 652 ierr = PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(PetscRealPart(val/scale)));CHKERRQ(ierr); 653 ierr = VecRestoreSubVector(s->y, s->isg[0], &y0);CHKERRQ(ierr); 654 655 /* error in pressure */ 656 ierr = VecGetSubVector(s->y, s->isg[1], &y1);CHKERRQ(ierr); 657 ierr = VecNorm(y1, NORM_2, &val);CHKERRQ(ierr); 658 ierr = PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(PetscRealPart(val/scale)));CHKERRQ(ierr); 659 ierr = VecRestoreSubVector(s->y, s->isg[1], &y1);CHKERRQ(ierr); 660 661 /* total error */ 662 ierr = VecNorm(s->y, NORM_2, &val);CHKERRQ(ierr); 663 ierr = PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)PetscRealPart((val/scale)));CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 int main(int argc, char **argv) 668 { 669 Stokes s; 670 KSP ksp; 671 PetscErrorCode ierr; 672 673 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 674 s.nx = 4; 675 s.ny = 6; 676 ierr = PetscOptionsGetInt(NULL,NULL, "-nx", &s.nx, NULL);CHKERRQ(ierr); 677 ierr = PetscOptionsGetInt(NULL,NULL, "-ny", &s.ny, NULL);CHKERRQ(ierr); 678 s.hx = 2.0/s.nx; 679 s.hy = 1.0/s.ny; 680 s.matsymmetric = PETSC_FALSE; 681 ierr = PetscOptionsGetBool(NULL,NULL, "-mat_set_symmetric", &s.matsymmetric,NULL);CHKERRQ(ierr); 682 s.userPC = s.userKSP = PETSC_FALSE; 683 ierr = PetscOptionsHasName(NULL,NULL, "-user_pc", &s.userPC);CHKERRQ(ierr); 684 ierr = PetscOptionsHasName(NULL,NULL, "-user_ksp", &s.userKSP);CHKERRQ(ierr); 685 686 ierr = StokesSetupMatrix(&s);CHKERRQ(ierr); 687 ierr = StokesSetupIndexSets(&s);CHKERRQ(ierr); 688 ierr = StokesSetupVectors(&s);CHKERRQ(ierr); 689 690 ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);CHKERRQ(ierr); 691 ierr = KSPSetOperators(ksp, s.A, s.A);CHKERRQ(ierr); 692 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 693 ierr = StokesSetupPC(&s, ksp);CHKERRQ(ierr); 694 ierr = KSPSolve(ksp, s.b, s.x);CHKERRQ(ierr); 695 696 /* don't trust, verify! */ 697 ierr = StokesCalcResidual(&s);CHKERRQ(ierr); 698 ierr = StokesCalcError(&s);CHKERRQ(ierr); 699 ierr = StokesWriteSolution(&s);CHKERRQ(ierr); 700 701 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 702 ierr = MatDestroy(&s.subA[0]);CHKERRQ(ierr); 703 ierr = MatDestroy(&s.subA[1]);CHKERRQ(ierr); 704 ierr = MatDestroy(&s.subA[2]);CHKERRQ(ierr); 705 ierr = MatDestroy(&s.subA[3]);CHKERRQ(ierr); 706 ierr = MatDestroy(&s.A);CHKERRQ(ierr); 707 ierr = VecDestroy(&s.x);CHKERRQ(ierr); 708 ierr = VecDestroy(&s.b);CHKERRQ(ierr); 709 ierr = VecDestroy(&s.y);CHKERRQ(ierr); 710 ierr = MatDestroy(&s.myS);CHKERRQ(ierr); 711 ierr = PetscFinalize(); 712 return ierr; 713 } 714 715 /*TEST 716 717 test: 718 nsize: 2 719 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 720 721 test: 722 suffix: 2 723 nsize: 2 724 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 725 726 test: 727 suffix: 3 728 nsize: 2 729 args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc 730 731 test: 732 suffix: 4 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_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 735 736 test: 737 suffix: 4_pcksp 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 -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 740 741 test: 742 suffix: 5 743 nsize: 2 744 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 745 746 test: 747 suffix: 6 748 nsize: 2 749 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 750 751 test: 752 suffix: 7 753 nsize: 2 754 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 755 756 TEST*/ 757