Actual source code: ex70.c

petsc-dev 2014-04-19
Report Typos and Errors
  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";

  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.                                                               */
 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: /*                                                                             */
 27: /* mpiexec -n 2 ./stokes -nx 32 -ny 48                                         */
 28: /*                                                                             */
 29: /*   Runs with PETSc defaults on 32x48 grid, no PC for the Schur               */
 30: /*   complement because A11 is zero.                                           */
 31: /*                                                                             */
 32: /* mpiexec -n 2 ./stokes -nx 32 -ny 48 -fieldsplit_1_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 ./stokes -nx 32 -ny 48 -fieldsplit_1_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 ./stokes -nx 32 -ny 48 -fieldsplit_1_user_ksp -fieldsplit_1_ksp_rtol 0.01 -fieldsplit_0_ksp_rtol 0.01 */
 45: /*                                                                             */
 46: /*   SIMPLE-type approximations are crude, there's no benefit in               */
 47: /*   solving the subsystems in the preconditioner very accurately.             */
 48: /*                                                                             */
 49: /*---------------------------------------------------------------------------- */

 51: #include <petscksp.h>

 53: typedef struct {
 54:   PetscBool userPC, userKSP; /* 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;

 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) */

 69: PetscErrorCode StokesGetPosition(Stokes*, PetscInt, PetscInt*, PetscInt*); /* row number j*nx+i corresponds to position (i,j) in grid */

 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) */

 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 */

 80: PetscErrorCode StokesSetupApproxSchur(Stokes*);  /* approximation of the Schur complement */

 82: PetscErrorCode StokesExactSolution(Stokes*); /* exact solution vector */
 83: PetscErrorCode StokesWriteSolution(Stokes*); /* write solution to file */

 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: }

 91: /* exact solution for the pressure */
 92: PetscScalar StokesExactPressure(const PetscScalar x)
 93: {
 94:   return 8.0*(2.0-x);
 95: }

 97: PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp)
 98: {
 99:   KSP            *subksp;
100:   PC             pc;
101:   PetscInt       n = 1;

105:   KSPGetPC(ksp, &pc);
106:   PCFieldSplitSetIS(pc, "0", s->isg[0]);
107:   PCFieldSplitSetIS(pc, "1", s->isg[1]);
108:   if (s->userPC) {
109:     PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS);
110:   }
111:   if (s->userKSP) {
112:     PCSetUp(pc);
113:     PCFieldSplitGetSubKSP(pc, &n, &subksp);
114:     KSPSetOperators(subksp[1], s->myS, s->myS);
115:     PetscFree(subksp);
116:   }
117:   return(0);
118: }

120: PetscErrorCode StokesWriteSolution(Stokes *s)
121: {
122:   PetscMPIInt    size;
123:   PetscInt       n,i,j;
124:   PetscScalar    *array;

128:   /* write data (*warning* only works sequential) */
129:   MPI_Comm_size(MPI_COMM_WORLD,&size);
130:   /*PetscPrintf(PETSC_COMM_WORLD," number of processors = %D\n",size);*/
131:   if (size == 1) {
132:     PetscViewer viewer;
133:     VecGetArray(s->x, &array);
134:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer);
135:     PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n");
136:     for (j = 0; j < s->ny; j++) {
137:       for (i = 0; i < s->nx; i++) {
138:         n    = j*s->nx+i;
139:         PetscViewerASCIIPrintf(viewer, "%.12g %.12g %.12g %.12g %.12g\n", i*s->hx+s->hx/2, j*s->hy+s->hy/2, array[n], array[n+s->nx*s->ny], array[n+2*s->nx*s->ny]);
140:       }
141:     }
142:     VecRestoreArray(s->x, &array);
143:     PetscViewerDestroy(&viewer);
144:   }
145:   return(0);
146: }

148: PetscErrorCode StokesSetupIndexSets(Stokes *s)
149: {

153:   /* the two index sets */
154:   MatNestGetISs(s->A, s->isg, NULL);
155:   /*  ISView(isg[0],PETSC_VIEWER_STDOUT_WORLD); */
156:   /*  ISView(isg[1],PETSC_VIEWER_STDOUT_WORLD); */
157:   return(0);
158: }

160: PetscErrorCode StokesSetupVectors(Stokes *s)
161: {

165:   /* solution vector x */
166:   VecCreate(PETSC_COMM_WORLD, &s->x);
167:   VecSetSizes(s->x, PETSC_DECIDE, 3*s->nx*s->ny);
168:   VecSetType(s->x, VECMPI);
169:   /*  VecSetRandom(s->x, NULL); */
170:   /*  VecView(s->x, (PetscViewer) PETSC_VIEWER_DEFAULT); */

172:   /* exact solution y */
173:   VecDuplicate(s->x, &s->y);
174:   StokesExactSolution(s);
175:   /*  VecView(s->y, (PetscViewer) PETSC_VIEWER_DEFAULT); */

177:   /* rhs vector b */
178:   VecDuplicate(s->x, &s->b);
179:   StokesRhs(s);
180:   /*VecView(s->b, (PetscViewer) PETSC_VIEWER_DEFAULT);*/
181:   return(0);
182: }

184: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
185: {
186:   PetscInt n;

189:   /* cell number n=j*nx+i has position (i,j) in grid */
190:   n  = row%(s->nx*s->ny);
191:   *i = n%s->nx;
192:   *j = (n-(*i))/s->nx;
193:   return(0);
194: }

196: PetscErrorCode StokesExactSolution(Stokes *s)
197: {
198:   PetscInt       row, start, end, i, j;
199:   PetscScalar    val;
200:   Vec            y0,y1;

204:   /* velocity part */
205:   VecGetSubVector(s->y, s->isg[0], &y0);
206:   VecGetOwnershipRange(y0, &start, &end);
207:   for (row = start; row < end; row++) {
208:     StokesGetPosition(s, row,&i,&j);
209:     if (row < s->nx*s->ny) {
210:       val = StokesExactVelocityX(j*s->hy+s->hy/2);
211:     } else {
212:       val = 0;
213:     }
214:     VecSetValue(y0, row, val, INSERT_VALUES);
215:   }
216:   VecRestoreSubVector(s->y, s->isg[0], &y0);

218:   /* pressure part */
219:   VecGetSubVector(s->y, s->isg[1], &y1);
220:   VecGetOwnershipRange(y1, &start, &end);
221:   for (row = start; row < end; row++) {
222:     StokesGetPosition(s, row, &i, &j);
223:     val  = StokesExactPressure(i*s->hx+s->hx/2);
224:     VecSetValue(y1, row, val, INSERT_VALUES);
225:   }
226:   VecRestoreSubVector(s->y, s->isg[1], &y1);
227:   return(0);
228: }

230: PetscErrorCode StokesRhs(Stokes *s)
231: {
232:   PetscInt       row, start, end, i, j;
233:   PetscScalar    val;
234:   Vec            b0,b1;

238:   /* velocity part */
239:   VecGetSubVector(s->b, s->isg[0], &b0);
240:   VecGetOwnershipRange(b0, &start, &end);
241:   for (row = start; row < end; row++) {
242:     StokesGetPosition(s, row, &i, &j);
243:     if (row < s->nx*s->ny) {
244:       StokesRhsMomX(s, i, j, &val);
245:     } else if (row < 2*s->nx*s->ny) {
246:       StokesRhsMomY(s, i, j, &val);
247:     }
248:     VecSetValue(b0, row, val, INSERT_VALUES);
249:   }
250:   VecRestoreSubVector(s->b, s->isg[0], &b0);

252:   /* pressure part */
253:   VecGetSubVector(s->b, s->isg[1], &b1);
254:   VecGetOwnershipRange(b1, &start, &end);
255:   for (row = start; row < end; row++) {
256:     StokesGetPosition(s, row, &i, &j);
257:     StokesRhsMass(s, i, j, &val);
258:     VecSetValue(b1, row, val, INSERT_VALUES);
259:   }
260:   VecRestoreSubVector(s->b, s->isg[1], &b1);
261:   return(0);
262: }

264: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
265: {
266:   PetscInt       row, start, end, size, i, j;
267:   PetscInt       cols[5];
268:   PetscScalar    vals[5];

272:   /* A[0] is 2N-by-2N */
273:   MatCreate(PETSC_COMM_WORLD,&s->subA[0]);
274:   MatSetOptionsPrefix(s->subA[0],"a00_");
275:   MatSetSizes(s->subA[0],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,2*s->nx*s->ny);
276:   MatSetType(s->subA[0],MATMPIAIJ);
277:   MatMPIAIJSetPreallocation(s->subA[0],5,NULL,5,NULL);
278:   MatGetOwnershipRange(s->subA[0], &start, &end);

280:   for (row = start; row < end; row++) {
281:     StokesGetPosition(s, row, &i, &j);
282:     /* first part: rows 0 to (nx*ny-1) */
283:     StokesStencilLaplacian(s, i, j, &size, cols, vals);
284:     /* second part: rows (nx*ny) to (2*nx*ny-1) */
285:     if (row >= s->nx*s->ny) {
286:       for (i = 0; i < 5; i++) cols[i] = cols[i] + s->nx*s->ny;
287:     }
288:     for (i = 0; i < 5; i++) vals[i] = -1.0*vals[i]; /* dynamic viscosity coef mu=-1 */
289:     MatSetValues(s->subA[0], 1, &row, size, cols, vals, INSERT_VALUES);
290:   }
291:   MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY);
292:   MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY);
293:   return(0);
294: }

296: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
297: {
298:   PetscInt       row, start, end, size, i, j;
299:   PetscInt       cols[5];
300:   PetscScalar    vals[5];

304:   /* A[1] is 2N-by-N */
305:   MatCreate(PETSC_COMM_WORLD, &s->subA[1]);
306:   MatSetOptionsPrefix(s->subA[1],"a01_");
307:   MatSetSizes(s->subA[1],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,s->nx*s->ny);
308:   MatSetType(s->subA[1],MATMPIAIJ);
309:   MatMPIAIJSetPreallocation(s->subA[1],5,NULL,5,NULL);
310:   MatGetOwnershipRange(s->subA[1],&start,&end);

312:   MatSetOption(s->subA[1],MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);

314:   for (row = start; row < end; row++) {
315:     StokesGetPosition(s, row, &i, &j);
316:     /* first part: rows 0 to (nx*ny-1) */
317:     if (row < s->nx*s->ny) {
318:       StokesStencilGradientX(s, i, j, &size, cols, vals);
319:     }
320:     /* second part: rows (nx*ny) to (2*nx*ny-1) */
321:     else {
322:       StokesStencilGradientY(s, i, j, &size, cols, vals);
323:     }
324:     MatSetValues(s->subA[1], 1, &row, size, cols, vals, INSERT_VALUES);
325:   }
326:   MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY);
327:   MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY);
328:   return(0);
329: }

331: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
332: {

336:   /* A[2] is minus transpose of A[1] */
337:   MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]);
338:   MatScale(s->subA[2], -1.0);
339:   MatSetOptionsPrefix(s->subA[2], "a10_");
340:   return(0);
341: }

343: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
344: {

348:   /* A[3] is N-by-N null matrix */
349:   MatCreate(PETSC_COMM_WORLD, &s->subA[3]);
350:   MatSetOptionsPrefix(s->subA[3], "a11_");
351:   MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx*s->ny, s->nx*s->ny);
352:   MatSetType(s->subA[3], MATMPIAIJ);
353:   MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL);
354:   MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY);
355:   MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY);
356:   return(0);
357: }

359: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
360: {
361:   Vec            diag;

365:   /* Schur complement approximation: myS = A11 - A10 diag(A00)^(-1) A01 */
366:   /* note: A11 is zero */
367:   /* note: in real life this matrix would be build directly, */
368:   /* i.e. without MatMatMult */

370:   /* inverse of diagonal of A00 */
371:   VecCreate(PETSC_COMM_WORLD,&diag);
372:   VecSetSizes(diag,PETSC_DECIDE,2*s->nx*s->ny);
373:   VecSetType(diag,VECMPI);
374:   MatGetDiagonal(s->subA[0],diag);
375:   VecReciprocal(diag);

377:   /* compute: - A10 diag(A00)^(-1) A01 */
378:   MatDiagonalScale(s->subA[1],diag,NULL); /* (*warning* overwrites subA[1]) */
379:   MatMatMult(s->subA[2],s->subA[1],MAT_INITIAL_MATRIX,PETSC_DEFAULT,&s->myS);
380:   MatScale(s->myS,-1.0);

382:   /* restore A10 */
383:   MatGetDiagonal(s->subA[0],diag);
384:   MatDiagonalScale(s->subA[1],diag,NULL);
385:   VecDestroy(&diag);
386:   return(0);
387: }

389: PetscErrorCode StokesSetupMatrix(Stokes *s)
390: {

394:   StokesSetupMatBlock00(s);
395:   StokesSetupMatBlock01(s);
396:   StokesSetupMatBlock10(s);
397:   StokesSetupMatBlock11(s);
398:   MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A);
399:   StokesSetupApproxSchur(s);
400:   return(0);
401: }

403: PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *size, PetscInt *cols, PetscScalar *vals)
404: {
405:   PetscInt    p =j*s->nx+i, w=p-1, e=p+1, s2=p-s->nx, n=p+s->nx;
406:   PetscScalar ae=s->hy/s->hx, aeb=0;
407:   PetscScalar aw=s->hy/s->hx, awb=s->hy/(s->hx/2);
408:   PetscScalar as=s->hx/s->hy, asb=s->hx/(s->hy/2);
409:   PetscScalar an=s->hx/s->hy, anb=s->hx/(s->hy/2);

412:   if (i==0 && j==0) { /* south-west corner */
413:     *size  =3;
414:     cols[0]=p; vals[0]=-(ae+awb+asb+an);
415:     cols[1]=e; vals[1]=ae;
416:     cols[2]=n; vals[2]=an;
417:   } else if (i==0 && j==s->ny-1) { /* north-west corner */
418:     *size  =3;
419:     cols[0]=s2; vals[0]=as;
420:     cols[1]=p; vals[1]=-(ae+awb+as+anb);
421:     cols[2]=e; vals[2]=ae;
422:   } else if (i==s->nx-1 && j==0) { /* south-east corner */
423:     *size  =3;
424:     cols[0]=w; vals[0]=aw;
425:     cols[1]=p; vals[1]=-(aeb+aw+asb+an);
426:     cols[2]=n; vals[2]=an;
427:   } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
428:     *size  =3;
429:     cols[0]=s2; vals[0]=as;
430:     cols[1]=w; vals[1]=aw;
431:     cols[2]=p; vals[2]=-(aeb+aw+as+anb);
432:   } else if (i==0) { /* west boundary */
433:     *size  =4;
434:     cols[0]=s2; vals[0]=as;
435:     cols[1]=p; vals[1]=-(ae+awb+as+an);
436:     cols[2]=e; vals[2]=ae;
437:     cols[3]=n; vals[3]=an;
438:   } else if (i==s->nx-1) { /* east boundary */
439:     *size  =4;
440:     cols[0]=s2; vals[0]=as;
441:     cols[1]=w; vals[1]=aw;
442:     cols[2]=p; vals[2]=-(aeb+aw+as+an);
443:     cols[3]=n; vals[3]=an;
444:   } else if (j==0) { /* south boundary */
445:     *size  =4;
446:     cols[0]=w; vals[0]=aw;
447:     cols[1]=p; vals[1]=-(ae+aw+asb+an);
448:     cols[2]=e; vals[2]=ae;
449:     cols[3]=n; vals[3]=an;
450:   } else if (j==s->ny-1) { /* north boundary */
451:     *size  =4;
452:     cols[0]=s2; vals[0]=as;
453:     cols[1]=w; vals[1]=aw;
454:     cols[2]=p; vals[2]=-(ae+aw+as+anb);
455:     cols[3]=e; vals[3]=ae;
456:   } else { /* interior */
457:     *size  =5;
458:     cols[0]=s2; vals[0]=as;
459:     cols[1]=w; vals[1]=aw;
460:     cols[2]=p; vals[2]=-(ae+aw+as+an);
461:     cols[3]=e; vals[3]=ae;
462:     cols[4]=n; vals[4]=an;
463:   }
464:   return(0);
465: }

467: PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *size, PetscInt *cols, PetscScalar *vals)
468: {
469:   PetscInt    p =j*s->nx+i, w=p-1, e=p+1;
470:   PetscScalar ae= s->hy/2, aeb=s->hy;
471:   PetscScalar aw=-s->hy/2, awb=0;

474:   if (i==0 && j==0) { /* south-west corner */
475:     *size  =2;
476:     cols[0]=p; vals[0]=-(ae+awb);
477:     cols[1]=e; vals[1]=ae;
478:   } else if (i==0 && j==s->ny-1) { /* north-west corner */
479:     *size  =2;
480:     cols[0]=p; vals[0]=-(ae+awb);
481:     cols[1]=e; vals[1]=ae;
482:   } else if (i==s->nx-1 && j==0) { /* south-east corner */
483:     *size  =2;
484:     cols[0]=w; vals[0]=aw;
485:     cols[1]=p; vals[1]=-(aeb+aw);
486:   } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
487:     *size  =2;
488:     cols[0]=w; vals[0]=aw;
489:     cols[1]=p; vals[1]=-(aeb+aw);
490:   } else if (i==0) { /* west boundary */
491:     *size  =2;
492:     cols[0]=p; vals[0]=-(ae+awb);
493:     cols[1]=e; vals[1]=ae;
494:   } else if (i==s->nx-1) { /* east boundary */
495:     *size  =2;
496:     cols[0]=w; vals[0]=aw;
497:     cols[1]=p; vals[1]=-(aeb+aw);
498:   } else if (j==0) { /* south boundary */
499:     *size  =3;
500:     cols[0]=w; vals[0]=aw;
501:     cols[1]=p; vals[1]=-(ae+aw);
502:     cols[2]=e; vals[2]=ae;
503:   } else if (j==s->ny-1) { /* north boundary */
504:     *size  =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 { /* interior */
509:     *size  =3;
510:     cols[0]=w; vals[0]=aw;
511:     cols[1]=p; vals[1]=-(ae+aw);
512:     cols[2]=e; vals[2]=ae;
513:   }
514:   return(0);
515: }

517: PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *size, PetscInt *cols, PetscScalar *vals)
518: {
519:   PetscInt    p =j*s->nx+i, s2=p-s->nx, n=p+s->nx;
520:   PetscScalar as=-s->hx/2, asb=0;
521:   PetscScalar an= s->hx/2, anb=0;

524:   if (i==0 && j==0) { /* south-west corner */
525:     *size  =2;
526:     cols[0]=p; vals[0]=-(asb+an);
527:     cols[1]=n; vals[1]=an;
528:   } else if (i==0 && j==s->ny-1) { /* north-west corner */
529:     *size  =2;
530:     cols[0]=s2; vals[0]=as;
531:     cols[1]=p; vals[1]=-(as+anb);
532:   } else if (i==s->nx-1 && j==0) { /* south-east corner */
533:     *size  =2;
534:     cols[0]=p; vals[0]=-(asb+an);
535:     cols[1]=n; vals[1]=an;
536:   } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
537:     *size  =2;
538:     cols[0]=s2; vals[0]=as;
539:     cols[1]=p; vals[1]=-(as+anb);
540:   } else if (i==0) { /* west boundary */
541:     *size  =3;
542:     cols[0]=s2; vals[0]=as;
543:     cols[1]=p; vals[1]=-(as+an);
544:     cols[2]=n; vals[2]=an;
545:   } else if (i==s->nx-1) { /* east boundary */
546:     *size  =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 (j==0) { /* south boundary */
551:     *size  =2;
552:     cols[0]=p; vals[0]=-(asb+an);
553:     cols[1]=n; vals[1]=an;
554:   } else if (j==s->ny-1) { /* north boundary */
555:     *size  =2;
556:     cols[0]=s2; vals[0]=as;
557:     cols[1]=p; vals[1]=-(as+anb);
558:   } else { /* interior */
559:     *size  =3;
560:     cols[0]=s2; vals[0]=as;
561:     cols[1]=p; vals[1]=-(as+an);
562:     cols[2]=n; vals[2]=an;
563:   }
564:   return(0);
565: }

567: PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
568: {
569:   PetscScalar y   = j*s->hy+s->hy/2;
570:   PetscScalar awb = s->hy/(s->hx/2);

573:   if (i == 0) { /* west boundary */
574:     *val = awb*StokesExactVelocityX(y);
575:   } else {
576:     *val = 0.0;
577:   }
578:   return(0);
579: }

581: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
582: {
584:   *val = 0.0;
585:   return(0);
586: }

588: PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
589: {
590:   PetscScalar y   = j*s->hy+s->hy/2;
591:   PetscScalar aeb = s->hy;

594:   if (i == 0) { /* west boundary */
595:     *val = aeb*StokesExactVelocityX(y);
596:   } else {
597:     *val = 0.0;
598:   }
599:   return(0);
600: }

602: PetscErrorCode StokesCalcResidual(Stokes *s)
603: {
604:   PetscReal      val;
605:   Vec            b0, b1;

609:   /* residual Ax-b (*warning* overwrites b) */
610:   VecScale(s->b, -1.0);
611:   MatMultAdd(s->A, s->x, s->b, s->b);
612:   /*  VecView(s->b, (PetscViewer)PETSC_VIEWER_DEFAULT); */

614:   /* residual velocity */
615:   VecGetSubVector(s->b, s->isg[0], &b0);
616:   VecNorm(b0, NORM_2, &val);
617:   PetscPrintf(PETSC_COMM_WORLD," residual u = %g\n",(double)val);
618:   VecRestoreSubVector(s->b, s->isg[0], &b0);

620:   /* residual pressure */
621:   VecGetSubVector(s->b, s->isg[1], &b1);
622:   VecNorm(b1, NORM_2, &val);
623:   PetscPrintf(PETSC_COMM_WORLD," residual p = %g\n",(double)val);
624:   VecRestoreSubVector(s->b, s->isg[1], &b1);

626:   /* total residual */
627:   VecNorm(s->b, NORM_2, &val);
628:   PetscPrintf(PETSC_COMM_WORLD," residual [u,p] = %g\n", (double)val);
629:   return(0);
630: }

632: PetscErrorCode StokesCalcError(Stokes *s)
633: {
634:   PetscScalar    scale = PetscSqrtScalar(s->nx*s->ny);
635:   PetscReal      val;
636:   Vec            y0, y1;

640:   /* error y-x */
641:   VecAXPY(s->y, -1.0, s->x);
642:   /* VecView(s->y, (PetscViewer)PETSC_VIEWER_DEFAULT); */

644:   /* error in velocity */
645:   VecGetSubVector(s->y, s->isg[0], &y0);
646:   VecNorm(y0, NORM_2, &val);
647:   PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(val/scale));
648:   VecRestoreSubVector(s->y, s->isg[0], &y0);

650:   /* error in pressure */
651:   VecGetSubVector(s->y, s->isg[1], &y1);
652:   VecNorm(y1, NORM_2, &val);
653:   PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(val/scale));
654:   VecRestoreSubVector(s->y, s->isg[1], &y1);

656:   /* total error */
657:   VecNorm(s->y, NORM_2, &val);
658:   PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)(val/scale));
659:   return(0);
660: }

662: int main(int argc, char **argv)
663: {
664:   Stokes         s;
665:   KSP            ksp;

668:   PetscInitialize(&argc, &argv, NULL, help);
669:   s.nx     = 4;
670:   s.ny     = 6;
671:   PetscOptionsGetInt(NULL, "-nx", &s.nx, NULL);
672:   PetscOptionsGetInt(NULL, "-ny", &s.ny, NULL);
673:   s.hx     = 2.0/s.nx;
674:   s.hy     = 1.0/s.ny;
675:   s.userPC = s.userKSP = PETSC_FALSE;
676:   PetscOptionsHasName(NULL, "-user_pc", &s.userPC);
677:   PetscOptionsHasName(NULL, "-user_ksp", &s.userKSP);

679:   StokesSetupMatrix(&s);
680:   StokesSetupIndexSets(&s);
681:   StokesSetupVectors(&s);

683:   KSPCreate(PETSC_COMM_WORLD, &ksp);
684:   KSPSetOperators(ksp, s.A, s.A);
685:   KSPSetFromOptions(ksp);
686:   StokesSetupPC(&s, ksp);
687:   KSPSolve(ksp, s.b, s.x);

689:   /* don't trust, verify! */
690:   StokesCalcResidual(&s);
691:   StokesCalcError(&s);
692:   StokesWriteSolution(&s);

694:   KSPDestroy(&ksp);
695:   MatDestroy(&s.subA[0]);
696:   MatDestroy(&s.subA[1]);
697:   MatDestroy(&s.subA[2]);
698:   MatDestroy(&s.subA[3]);
699:   MatDestroy(&s.A);
700:   VecDestroy(&s.x);
701:   VecDestroy(&s.b);
702:   VecDestroy(&s.y);
703:   MatDestroy(&s.myS);
704:   PetscFinalize();
705:   return 0;
706: }