Actual source code: ex70.c

petsc-3.14.6 2021-03-30
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. 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: /*---------------------------------------------------------------------------- */


 56: #include <petscksp.h>

 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;

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

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

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

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

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

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

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

 96: /* exact solution for the pressure */
 97: PetscScalar StokesExactPressure(const PetscScalar x)
 98: {
 99:   return 8.0*(2.0-x);
100: }

102: PetscErrorCode StokesSetupPC(Stokes *s, KSP ksp)
103: {
104:   KSP            *subksp;
105:   PC             pc;
106:   PetscInt       n = 1;

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

125: PetscErrorCode StokesWriteSolution(Stokes *s)
126: {
127:   PetscMPIInt       size;
128:   PetscInt          n,i,j;
129:   const PetscScalar *array;
130:   PetscErrorCode    ierr;

133:   /* write data (*warning* only works sequential) */
134:   MPI_Comm_size(MPI_COMM_WORLD,&size);
135:   /*PetscPrintf(PETSC_COMM_WORLD," number of processors = %D\n",size);*/
136:   if (size == 1) {
137:     PetscViewer viewer;
138:     VecGetArrayRead(s->x, &array);
139:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer);
140:     PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n");
141:     for (j = 0; j < s->ny; j++) {
142:       for (i = 0; i < s->nx; i++) {
143:         n    = j*s->nx+i;
144:         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]));
145:       }
146:     }
147:     VecRestoreArrayRead(s->x, &array);
148:     PetscViewerDestroy(&viewer);
149:   }
150:   return(0);
151: }

153: PetscErrorCode StokesSetupIndexSets(Stokes *s)
154: {

158:   /* the two index sets */
159:   MatNestGetISs(s->A, s->isg, NULL);
160:   /*  ISView(isg[0],PETSC_VIEWER_STDOUT_WORLD); */
161:   /*  ISView(isg[1],PETSC_VIEWER_STDOUT_WORLD); */
162:   return(0);
163: }

165: PetscErrorCode StokesSetupVectors(Stokes *s)
166: {

170:   /* solution vector x */
171:   VecCreate(PETSC_COMM_WORLD, &s->x);
172:   VecSetSizes(s->x, PETSC_DECIDE, 3*s->nx*s->ny);
173:   VecSetType(s->x, VECMPI);
174:   /*  VecSetRandom(s->x, NULL); */
175:   /*  VecView(s->x, (PetscViewer) PETSC_VIEWER_DEFAULT); */

177:   /* exact solution y */
178:   VecDuplicate(s->x, &s->y);
179:   StokesExactSolution(s);
180:   /*  VecView(s->y, (PetscViewer) PETSC_VIEWER_DEFAULT); */

182:   /* rhs vector b */
183:   VecDuplicate(s->x, &s->b);
184:   StokesRhs(s);
185:   /*VecView(s->b, (PetscViewer) PETSC_VIEWER_DEFAULT);*/
186:   return(0);
187: }

189: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
190: {
191:   PetscInt n;

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:   return(0);
199: }

201: PetscErrorCode StokesExactSolution(Stokes *s)
202: {
203:   PetscInt       row, start, end, i, j;
204:   PetscScalar    val;
205:   Vec            y0,y1;

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

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

235: PetscErrorCode StokesRhs(Stokes *s)
236: {
237:   PetscInt       row, start, end, i, j;
238:   PetscScalar    val;
239:   Vec            b0,b1;

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

257:   /* pressure part */
258:   VecGetSubVector(s->b, s->isg[1], &b1);
259:   VecGetOwnershipRange(b1, &start, &end);
260:   for (row = start; row < end; row++) {
261:     StokesGetPosition(s, row, &i, &j);
262:     StokesRhsMass(s, i, j, &val);
263:     if (s->matsymmetric) {
264:       val = -1.0*val;
265:     }
266:     VecSetValue(b1, row, val, INSERT_VALUES);
267:   }
268:   VecRestoreSubVector(s->b, s->isg[1], &b1);
269:   return(0);
270: }

272: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
273: {
274:   PetscInt       row, start, end, sz, i, j;
275:   PetscInt       cols[5];
276:   PetscScalar    vals[5];

280:   /* A[0] is 2N-by-2N */
281:   MatCreate(PETSC_COMM_WORLD,&s->subA[0]);
282:   MatSetOptionsPrefix(s->subA[0],"a00_");
283:   MatSetSizes(s->subA[0],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,2*s->nx*s->ny);
284:   MatSetType(s->subA[0],MATMPIAIJ);
285:   MatMPIAIJSetPreallocation(s->subA[0],5,NULL,5,NULL);
286:   MatGetOwnershipRange(s->subA[0], &start, &end);

288:   for (row = start; row < end; row++) {
289:     StokesGetPosition(s, row, &i, &j);
290:     /* first part: rows 0 to (nx*ny-1) */
291:     StokesStencilLaplacian(s, i, j, &sz, cols, vals);
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:     MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES);
298:   }
299:   MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY);
300:   MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY);
301:   return(0);
302: }

304: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
305: {
306:   PetscInt       row, start, end, sz, i, j;
307:   PetscInt       cols[5];
308:   PetscScalar    vals[5];

312:   /* A[1] is 2N-by-N */
313:   MatCreate(PETSC_COMM_WORLD, &s->subA[1]);
314:   MatSetOptionsPrefix(s->subA[1],"a01_");
315:   MatSetSizes(s->subA[1],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,s->nx*s->ny);
316:   MatSetType(s->subA[1],MATMPIAIJ);
317:   MatMPIAIJSetPreallocation(s->subA[1],5,NULL,5,NULL);
318:   MatGetOwnershipRange(s->subA[1],&start,&end);

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

322:   for (row = start; row < end; row++) {
323:     StokesGetPosition(s, row, &i, &j);
324:     /* first part: rows 0 to (nx*ny-1) */
325:     if (row < s->nx*s->ny) {
326:       StokesStencilGradientX(s, i, j, &sz, cols, vals);
327:     } else {    /* second part: rows (nx*ny) to (2*nx*ny-1) */
328:       StokesStencilGradientY(s, i, j, &sz, cols, vals);
329:     }
330:     MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES);
331:   }
332:   MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY);
333:   MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY);
334:   return(0);
335: }

337: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
338: {

342:   /* A[2] is minus transpose of A[1] */
343:   MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]);
344:   if (!s->matsymmetric) {
345:     MatScale(s->subA[2], -1.0);
346:   }
347:   MatSetOptionsPrefix(s->subA[2], "a10_");
348:   return(0);
349: }

351: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
352: {

356:   /* A[3] is N-by-N null matrix */
357:   MatCreate(PETSC_COMM_WORLD, &s->subA[3]);
358:   MatSetOptionsPrefix(s->subA[3], "a11_");
359:   MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx*s->ny, s->nx*s->ny);
360:   MatSetType(s->subA[3], MATMPIAIJ);
361:   MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL);
362:   MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY);
363:   MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY);
364:   return(0);
365: }

367: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
368: {
369:   Vec            diag;

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

378:   /* inverse of diagonal of A00 */
379:   VecCreate(PETSC_COMM_WORLD,&diag);
380:   VecSetSizes(diag,PETSC_DECIDE,2*s->nx*s->ny);
381:   VecSetType(diag,VECMPI);
382:   MatGetDiagonal(s->subA[0],diag);
383:   VecReciprocal(diag);

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

390:   /* restore A10 */
391:   MatGetDiagonal(s->subA[0],diag);
392:   MatDiagonalScale(s->subA[1],diag,NULL);
393:   VecDestroy(&diag);
394:   return(0);
395: }

397: PetscErrorCode StokesSetupMatrix(Stokes *s)
398: {

402:   StokesSetupMatBlock00(s);
403:   StokesSetupMatBlock01(s);
404:   StokesSetupMatBlock10(s);
405:   StokesSetupMatBlock11(s);
406:   MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A);
407:   StokesSetupApproxSchur(s);
408:   return(0);
409: }

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);

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:   return(0);
473: }

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;

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:   return(0);
523: }

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;

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:   return(0);
573: }

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);

581:   if (i == 0) { /* west boundary */
582:     *val = awb*StokesExactVelocityX(y);
583:   } else {
584:     *val = 0.0;
585:   }
586:   return(0);
587: }

589: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
590: {
592:   *val = 0.0;
593:   return(0);
594: }

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;

602:   if (i == 0) { /* west boundary */
603:     *val = aeb*StokesExactVelocityX(y);
604:   } else {
605:     *val = 0.0;
606:   }
607:   return(0);
608: }

610: PetscErrorCode StokesCalcResidual(Stokes *s)
611: {
612:   PetscReal      val;
613:   Vec            b0, b1;

617:   /* residual Ax-b (*warning* overwrites b) */
618:   VecScale(s->b, -1.0);
619:   MatMultAdd(s->A, s->x, s->b, s->b);
620:   /*  VecView(s->b, (PetscViewer)PETSC_VIEWER_DEFAULT); */

622:   /* residual velocity */
623:   VecGetSubVector(s->b, s->isg[0], &b0);
624:   VecNorm(b0, NORM_2, &val);
625:   PetscPrintf(PETSC_COMM_WORLD," residual u = %g\n",(double)val);
626:   VecRestoreSubVector(s->b, s->isg[0], &b0);

628:   /* residual pressure */
629:   VecGetSubVector(s->b, s->isg[1], &b1);
630:   VecNorm(b1, NORM_2, &val);
631:   PetscPrintf(PETSC_COMM_WORLD," residual p = %g\n",(double)val);
632:   VecRestoreSubVector(s->b, s->isg[1], &b1);

634:   /* total residual */
635:   VecNorm(s->b, NORM_2, &val);
636:   PetscPrintf(PETSC_COMM_WORLD," residual [u,p] = %g\n", (double)val);
637:   return(0);
638: }

640: PetscErrorCode StokesCalcError(Stokes *s)
641: {
642:   PetscScalar    scale = PetscSqrtReal((double)s->nx*s->ny);
643:   PetscReal      val;
644:   Vec            y0, y1;

648:   /* error y-x */
649:   VecAXPY(s->y, -1.0, s->x);
650:   /* VecView(s->y, (PetscViewer)PETSC_VIEWER_DEFAULT); */

652:   /* error in velocity */
653:   VecGetSubVector(s->y, s->isg[0], &y0);
654:   VecNorm(y0, NORM_2, &val);
655:   PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(PetscRealPart(val/scale)));
656:   VecRestoreSubVector(s->y, s->isg[0], &y0);

658:   /* error in pressure */
659:   VecGetSubVector(s->y, s->isg[1], &y1);
660:   VecNorm(y1, NORM_2, &val);
661:   PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(PetscRealPart(val/scale)));
662:   VecRestoreSubVector(s->y, s->isg[1], &y1);

664:   /* total error */
665:   VecNorm(s->y, NORM_2, &val);
666:   PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)PetscRealPart((val/scale)));
667:   return(0);
668: }

670: int main(int argc, char **argv)
671: {
672:   Stokes         s;
673:   KSP            ksp;

676:   PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
677:   s.nx           = 4;
678:   s.ny           = 6;
679:   PetscOptionsGetInt(NULL,NULL, "-nx", &s.nx, NULL);
680:   PetscOptionsGetInt(NULL,NULL, "-ny", &s.ny, NULL);
681:   s.hx           = 2.0/s.nx;
682:   s.hy           = 1.0/s.ny;
683:   s.matsymmetric = PETSC_FALSE;
684:   PetscOptionsGetBool(NULL,NULL, "-mat_set_symmetric", &s.matsymmetric,NULL);
685:   s.userPC       = s.userKSP = PETSC_FALSE;
686:   PetscOptionsHasName(NULL,NULL, "-user_pc", &s.userPC);
687:   PetscOptionsHasName(NULL,NULL, "-user_ksp", &s.userKSP);

689:   StokesSetupMatrix(&s);
690:   StokesSetupIndexSets(&s);
691:   StokesSetupVectors(&s);

693:   KSPCreate(PETSC_COMM_WORLD, &ksp);
694:   KSPSetOperators(ksp, s.A, s.A);
695:   KSPSetFromOptions(ksp);
696:   StokesSetupPC(&s, ksp);
697:   KSPSolve(ksp, s.b, s.x);

699:   /* don't trust, verify! */
700:   StokesCalcResidual(&s);
701:   StokesCalcError(&s);
702:   StokesWriteSolution(&s);

704:   KSPDestroy(&ksp);
705:   MatDestroy(&s.subA[0]);
706:   MatDestroy(&s.subA[1]);
707:   MatDestroy(&s.subA[2]);
708:   MatDestroy(&s.subA[3]);
709:   MatDestroy(&s.A);
710:   VecDestroy(&s.x);
711:   VecDestroy(&s.b);
712:   VecDestroy(&s.y);
713:   MatDestroy(&s.myS);
714:   PetscFinalize();
715:   return ierr;
716: }


719: /*TEST

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

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

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

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

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

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

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

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


761: TEST*/