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