Actual source code: ex70.c
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*/