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