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: author : Christiaan M. Klaij
8: Poiseuille flow problem.
10: Viscous, laminar flow in a 2D channel with parabolic velocity
11: profile and linear pressure drop, exact solution of the 2D Stokes
12: equations.
14: Discretized with the cell-centered finite-volume method on a
15: Cartesian grid with co-located variables. Variables ordered as
16: [u1...uN v1...vN p1...pN]^T. Matrix [A00 A01; A10, A11] solved with
17: PCFIELDSPLIT. Lower factorization is used to mimic the Semi-Implicit
18: Method for Pressure Linked Equations (SIMPLE) used as preconditioner
19: instead of solver.
21: Disclaimer: does not contain the pressure-weighed interpolation
22: method needed to suppress spurious pressure modes in real-life
23: problems.
25: Usage:
26: 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
28: Runs with PCFIELDSPLIT on 32x48 grid, no PC for the Schur
29: complement because A11 is zero. FGMRES is needed because
30: PCFIELDSPLIT is a variable preconditioner.
32: 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
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).
38: 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
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.
44: 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
46: Out-of-the-box SIMPLE-type preconditioning. The major advantage
47: is that the user neither needs to provide the approximation of
48: the Schur complement, nor the corresponding preconditioner.
49: */
51: #include <petscksp.h>
53: typedef struct {
54: PetscBool userPC, userKSP, matsymmetric; /* 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;
104: KSPGetPC(ksp, &pc);
105: PCFieldSplitSetIS(pc, "0", s->isg[0]);
106: PCFieldSplitSetIS(pc, "1", s->isg[1]);
107: if (s->userPC) {
108: PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS);
109: }
110: if (s->userKSP) {
111: PCSetUp(pc);
112: PCFieldSplitGetSubKSP(pc, &n, &subksp);
113: KSPSetOperators(subksp[1], s->myS, s->myS);
114: PetscFree(subksp);
115: }
116: return 0;
117: }
119: PetscErrorCode StokesWriteSolution(Stokes *s)
120: {
121: PetscMPIInt size;
122: PetscInt n,i,j;
123: const PetscScalar *array;
126: /* write data (*warning* only works sequential) */
127: MPI_Comm_size(MPI_COMM_WORLD,&size);
128: if (size == 1) {
129: PetscViewer viewer;
130: VecGetArrayRead(s->x, &array);
131: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer);
132: PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n");
133: for (j = 0; j < s->ny; j++) {
134: for (i = 0; i < s->nx; i++) {
135: n = j*s->nx+i;
136: 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]));
137: }
138: }
139: VecRestoreArrayRead(s->x, &array);
140: PetscViewerDestroy(&viewer);
141: }
142: return 0;
143: }
145: PetscErrorCode StokesSetupIndexSets(Stokes *s)
146: {
148: /* the two index sets */
149: MatNestGetISs(s->A, s->isg, NULL);
150: return 0;
151: }
153: PetscErrorCode StokesSetupVectors(Stokes *s)
154: {
156: /* solution vector x */
157: VecCreate(PETSC_COMM_WORLD, &s->x);
158: VecSetSizes(s->x, PETSC_DECIDE, 3*s->nx*s->ny);
159: VecSetType(s->x, VECMPI);
161: /* exact solution y */
162: VecDuplicate(s->x, &s->y);
163: StokesExactSolution(s);
165: /* rhs vector b */
166: VecDuplicate(s->x, &s->b);
167: StokesRhs(s);
168: return 0;
169: }
171: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
172: {
173: PetscInt n;
176: /* cell number n=j*nx+i has position (i,j) in grid */
177: n = row%(s->nx*s->ny);
178: *i = n%s->nx;
179: *j = (n-(*i))/s->nx;
180: return 0;
181: }
183: PetscErrorCode StokesExactSolution(Stokes *s)
184: {
185: PetscInt row, start, end, i, j;
186: PetscScalar val;
187: Vec y0,y1;
190: /* velocity part */
191: VecGetSubVector(s->y, s->isg[0], &y0);
192: VecGetOwnershipRange(y0, &start, &end);
193: for (row = start; row < end; row++) {
194: StokesGetPosition(s, row,&i,&j);
195: if (row < s->nx*s->ny) {
196: val = StokesExactVelocityX(j*s->hy+s->hy/2);
197: } else {
198: val = 0;
199: }
200: VecSetValue(y0, row, val, INSERT_VALUES);
201: }
202: VecRestoreSubVector(s->y, s->isg[0], &y0);
204: /* pressure part */
205: VecGetSubVector(s->y, s->isg[1], &y1);
206: VecGetOwnershipRange(y1, &start, &end);
207: for (row = start; row < end; row++) {
208: StokesGetPosition(s, row, &i, &j);
209: val = StokesExactPressure(i*s->hx+s->hx/2);
210: VecSetValue(y1, row, val, INSERT_VALUES);
211: }
212: VecRestoreSubVector(s->y, s->isg[1], &y1);
213: return 0;
214: }
216: PetscErrorCode StokesRhs(Stokes *s)
217: {
218: PetscInt row, start, end, i, j;
219: PetscScalar val;
220: Vec b0,b1;
223: /* velocity part */
224: VecGetSubVector(s->b, s->isg[0], &b0);
225: VecGetOwnershipRange(b0, &start, &end);
226: for (row = start; row < end; row++) {
227: StokesGetPosition(s, row, &i, &j);
228: if (row < s->nx*s->ny) {
229: StokesRhsMomX(s, i, j, &val);
230: } else {
231: StokesRhsMomY(s, i, j, &val);
232: }
233: VecSetValue(b0, row, val, INSERT_VALUES);
234: }
235: VecRestoreSubVector(s->b, s->isg[0], &b0);
237: /* pressure part */
238: VecGetSubVector(s->b, s->isg[1], &b1);
239: VecGetOwnershipRange(b1, &start, &end);
240: for (row = start; row < end; row++) {
241: StokesGetPosition(s, row, &i, &j);
242: StokesRhsMass(s, i, j, &val);
243: if (s->matsymmetric) {
244: val = -1.0*val;
245: }
246: VecSetValue(b1, row, val, INSERT_VALUES);
247: }
248: VecRestoreSubVector(s->b, s->isg[1], &b1);
249: return 0;
250: }
252: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
253: {
254: PetscInt row, start, end, sz, i, j;
255: PetscInt cols[5];
256: PetscScalar vals[5];
259: /* A[0] is 2N-by-2N */
260: MatCreate(PETSC_COMM_WORLD,&s->subA[0]);
261: MatSetOptionsPrefix(s->subA[0],"a00_");
262: MatSetSizes(s->subA[0],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,2*s->nx*s->ny);
263: MatSetType(s->subA[0],MATMPIAIJ);
264: MatMPIAIJSetPreallocation(s->subA[0],5,NULL,5,NULL);
265: MatGetOwnershipRange(s->subA[0], &start, &end);
267: for (row = start; row < end; row++) {
268: StokesGetPosition(s, row, &i, &j);
269: /* first part: rows 0 to (nx*ny-1) */
270: StokesStencilLaplacian(s, i, j, &sz, cols, vals);
271: /* second part: rows (nx*ny) to (2*nx*ny-1) */
272: if (row >= s->nx*s->ny) {
273: for (i = 0; i < sz; i++) cols[i] += s->nx*s->ny;
274: }
275: for (i = 0; i < sz; i++) vals[i] = -1.0*vals[i]; /* dynamic viscosity coef mu=-1 */
276: MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES);
277: }
278: MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY);
279: MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY);
280: return 0;
281: }
283: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
284: {
285: PetscInt row, start, end, sz, i, j;
286: PetscInt cols[5];
287: PetscScalar vals[5];
290: /* A[1] is 2N-by-N */
291: MatCreate(PETSC_COMM_WORLD, &s->subA[1]);
292: MatSetOptionsPrefix(s->subA[1],"a01_");
293: MatSetSizes(s->subA[1],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,s->nx*s->ny);
294: MatSetType(s->subA[1],MATMPIAIJ);
295: MatMPIAIJSetPreallocation(s->subA[1],5,NULL,5,NULL);
296: MatGetOwnershipRange(s->subA[1],&start,&end);
298: MatSetOption(s->subA[1],MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
300: for (row = start; row < end; row++) {
301: StokesGetPosition(s, row, &i, &j);
302: /* first part: rows 0 to (nx*ny-1) */
303: if (row < s->nx*s->ny) {
304: StokesStencilGradientX(s, i, j, &sz, cols, vals);
305: } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */
306: StokesStencilGradientY(s, i, j, &sz, cols, vals);
307: }
308: MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES);
309: }
310: MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY);
311: MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY);
312: return 0;
313: }
315: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
316: {
318: /* A[2] is minus transpose of A[1] */
319: MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]);
320: if (!s->matsymmetric) {
321: MatScale(s->subA[2], -1.0);
322: }
323: MatSetOptionsPrefix(s->subA[2], "a10_");
324: return 0;
325: }
327: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
328: {
330: /* A[3] is N-by-N null matrix */
331: MatCreate(PETSC_COMM_WORLD, &s->subA[3]);
332: MatSetOptionsPrefix(s->subA[3], "a11_");
333: MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx*s->ny, s->nx*s->ny);
334: MatSetType(s->subA[3], MATMPIAIJ);
335: MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL);
336: MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY);
337: MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY);
338: return 0;
339: }
341: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
342: {
343: Vec diag;
346: /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */
347: /* note: A11 is zero */
348: /* note: in real life this matrix would be build directly, */
349: /* i.e. without MatMatMult */
351: /* inverse of diagonal of A00 */
352: VecCreate(PETSC_COMM_WORLD,&diag);
353: VecSetSizes(diag,PETSC_DECIDE,2*s->nx*s->ny);
354: VecSetType(diag,VECMPI);
355: MatGetDiagonal(s->subA[0],diag);
356: VecReciprocal(diag);
358: /* compute: - A10 inv(DIAGFORM(A00)) A01 */
359: MatDiagonalScale(s->subA[1],diag,NULL); /* (*warning* overwrites subA[1]) */
360: MatMatMult(s->subA[2],s->subA[1],MAT_INITIAL_MATRIX,PETSC_DEFAULT,&s->myS);
361: MatScale(s->myS,-1.0);
363: /* restore A10 */
364: MatGetDiagonal(s->subA[0],diag);
365: MatDiagonalScale(s->subA[1],diag,NULL);
366: VecDestroy(&diag);
367: return 0;
368: }
370: PetscErrorCode StokesSetupMatrix(Stokes *s)
371: {
373: StokesSetupMatBlock00(s);
374: StokesSetupMatBlock01(s);
375: StokesSetupMatBlock10(s);
376: StokesSetupMatBlock11(s);
377: MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A);
378: StokesSetupApproxSchur(s);
379: return 0;
380: }
382: PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
383: {
384: PetscInt p =j*s->nx+i, w=p-1, e=p+1, s2=p-s->nx, n=p+s->nx;
385: PetscScalar ae=s->hy/s->hx, aeb=0;
386: PetscScalar aw=s->hy/s->hx, awb=s->hy/(s->hx/2);
387: PetscScalar as=s->hx/s->hy, asb=s->hx/(s->hy/2);
388: PetscScalar an=s->hx/s->hy, anb=s->hx/(s->hy/2);
391: if (i==0 && j==0) { /* south-west corner */
392: *sz =3;
393: cols[0]=p; vals[0]=-(ae+awb+asb+an);
394: cols[1]=e; vals[1]=ae;
395: cols[2]=n; vals[2]=an;
396: } else if (i==0 && j==s->ny-1) { /* north-west corner */
397: *sz =3;
398: cols[0]=s2; vals[0]=as;
399: cols[1]=p; vals[1]=-(ae+awb+as+anb);
400: cols[2]=e; vals[2]=ae;
401: } else if (i==s->nx-1 && j==0) { /* south-east corner */
402: *sz =3;
403: cols[0]=w; vals[0]=aw;
404: cols[1]=p; vals[1]=-(aeb+aw+asb+an);
405: cols[2]=n; vals[2]=an;
406: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
407: *sz =3;
408: cols[0]=s2; vals[0]=as;
409: cols[1]=w; vals[1]=aw;
410: cols[2]=p; vals[2]=-(aeb+aw+as+anb);
411: } else if (i==0) { /* west boundary */
412: *sz =4;
413: cols[0]=s2; vals[0]=as;
414: cols[1]=p; vals[1]=-(ae+awb+as+an);
415: cols[2]=e; vals[2]=ae;
416: cols[3]=n; vals[3]=an;
417: } else if (i==s->nx-1) { /* east boundary */
418: *sz =4;
419: cols[0]=s2; vals[0]=as;
420: cols[1]=w; vals[1]=aw;
421: cols[2]=p; vals[2]=-(aeb+aw+as+an);
422: cols[3]=n; vals[3]=an;
423: } else if (j==0) { /* south boundary */
424: *sz =4;
425: cols[0]=w; vals[0]=aw;
426: cols[1]=p; vals[1]=-(ae+aw+asb+an);
427: cols[2]=e; vals[2]=ae;
428: cols[3]=n; vals[3]=an;
429: } else if (j==s->ny-1) { /* north boundary */
430: *sz =4;
431: cols[0]=s2; vals[0]=as;
432: cols[1]=w; vals[1]=aw;
433: cols[2]=p; vals[2]=-(ae+aw+as+anb);
434: cols[3]=e; vals[3]=ae;
435: } else { /* interior */
436: *sz =5;
437: cols[0]=s2; vals[0]=as;
438: cols[1]=w; vals[1]=aw;
439: cols[2]=p; vals[2]=-(ae+aw+as+an);
440: cols[3]=e; vals[3]=ae;
441: cols[4]=n; vals[4]=an;
442: }
443: return 0;
444: }
446: PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
447: {
448: PetscInt p =j*s->nx+i, w=p-1, e=p+1;
449: PetscScalar ae= s->hy/2, aeb=s->hy;
450: PetscScalar aw=-s->hy/2, awb=0;
453: if (i==0 && j==0) { /* south-west corner */
454: *sz =2;
455: cols[0]=p; vals[0]=-(ae+awb);
456: cols[1]=e; vals[1]=ae;
457: } else if (i==0 && j==s->ny-1) { /* north-west corner */
458: *sz =2;
459: cols[0]=p; vals[0]=-(ae+awb);
460: cols[1]=e; vals[1]=ae;
461: } else if (i==s->nx-1 && j==0) { /* south-east corner */
462: *sz =2;
463: cols[0]=w; vals[0]=aw;
464: cols[1]=p; vals[1]=-(aeb+aw);
465: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
466: *sz =2;
467: cols[0]=w; vals[0]=aw;
468: cols[1]=p; vals[1]=-(aeb+aw);
469: } else if (i==0) { /* west boundary */
470: *sz =2;
471: cols[0]=p; vals[0]=-(ae+awb);
472: cols[1]=e; vals[1]=ae;
473: } else if (i==s->nx-1) { /* east boundary */
474: *sz =2;
475: cols[0]=w; vals[0]=aw;
476: cols[1]=p; vals[1]=-(aeb+aw);
477: } else if (j==0) { /* south boundary */
478: *sz =3;
479: cols[0]=w; vals[0]=aw;
480: cols[1]=p; vals[1]=-(ae+aw);
481: cols[2]=e; vals[2]=ae;
482: } else if (j==s->ny-1) { /* north boundary */
483: *sz =3;
484: cols[0]=w; vals[0]=aw;
485: cols[1]=p; vals[1]=-(ae+aw);
486: cols[2]=e; vals[2]=ae;
487: } else { /* interior */
488: *sz =3;
489: cols[0]=w; vals[0]=aw;
490: cols[1]=p; vals[1]=-(ae+aw);
491: cols[2]=e; vals[2]=ae;
492: }
493: return 0;
494: }
496: PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
497: {
498: PetscInt p =j*s->nx+i, s2=p-s->nx, n=p+s->nx;
499: PetscScalar as=-s->hx/2, asb=0;
500: PetscScalar an= s->hx/2, anb=0;
503: if (i==0 && j==0) { /* south-west corner */
504: *sz =2;
505: cols[0]=p; vals[0]=-(asb+an);
506: cols[1]=n; vals[1]=an;
507: } else if (i==0 && j==s->ny-1) { /* north-west corner */
508: *sz =2;
509: cols[0]=s2; vals[0]=as;
510: cols[1]=p; vals[1]=-(as+anb);
511: } else if (i==s->nx-1 && j==0) { /* south-east corner */
512: *sz =2;
513: cols[0]=p; vals[0]=-(asb+an);
514: cols[1]=n; vals[1]=an;
515: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
516: *sz =2;
517: cols[0]=s2; vals[0]=as;
518: cols[1]=p; vals[1]=-(as+anb);
519: } else if (i==0) { /* west boundary */
520: *sz =3;
521: cols[0]=s2; vals[0]=as;
522: cols[1]=p; vals[1]=-(as+an);
523: cols[2]=n; vals[2]=an;
524: } else if (i==s->nx-1) { /* east boundary */
525: *sz =3;
526: cols[0]=s2; vals[0]=as;
527: cols[1]=p; vals[1]=-(as+an);
528: cols[2]=n; vals[2]=an;
529: } else if (j==0) { /* south boundary */
530: *sz =2;
531: cols[0]=p; vals[0]=-(asb+an);
532: cols[1]=n; vals[1]=an;
533: } else if (j==s->ny-1) { /* north boundary */
534: *sz =2;
535: cols[0]=s2; vals[0]=as;
536: cols[1]=p; vals[1]=-(as+anb);
537: } else { /* interior */
538: *sz =3;
539: cols[0]=s2; vals[0]=as;
540: cols[1]=p; vals[1]=-(as+an);
541: cols[2]=n; vals[2]=an;
542: }
543: return 0;
544: }
546: PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
547: {
548: PetscScalar y = j*s->hy+s->hy/2;
549: PetscScalar awb = s->hy/(s->hx/2);
552: if (i == 0) { /* west boundary */
553: *val = awb*StokesExactVelocityX(y);
554: } else {
555: *val = 0.0;
556: }
557: return 0;
558: }
560: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
561: {
563: *val = 0.0;
564: return 0;
565: }
567: PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
568: {
569: PetscScalar y = j*s->hy+s->hy/2;
570: PetscScalar aeb = s->hy;
573: if (i == 0) { /* west boundary */
574: *val = aeb*StokesExactVelocityX(y);
575: } else {
576: *val = 0.0;
577: }
578: return 0;
579: }
581: PetscErrorCode StokesCalcResidual(Stokes *s)
582: {
583: PetscReal val;
584: Vec b0, b1;
587: /* residual Ax-b (*warning* overwrites b) */
588: VecScale(s->b, -1.0);
589: MatMultAdd(s->A, s->x, s->b, s->b);
591: /* residual velocity */
592: VecGetSubVector(s->b, s->isg[0], &b0);
593: VecNorm(b0, NORM_2, &val);
594: PetscPrintf(PETSC_COMM_WORLD," residual u = %g\n",(double)val);
595: VecRestoreSubVector(s->b, s->isg[0], &b0);
597: /* residual pressure */
598: VecGetSubVector(s->b, s->isg[1], &b1);
599: VecNorm(b1, NORM_2, &val);
600: PetscPrintf(PETSC_COMM_WORLD," residual p = %g\n",(double)val);
601: VecRestoreSubVector(s->b, s->isg[1], &b1);
603: /* total residual */
604: VecNorm(s->b, NORM_2, &val);
605: PetscPrintf(PETSC_COMM_WORLD," residual [u,p] = %g\n", (double)val);
606: return 0;
607: }
609: PetscErrorCode StokesCalcError(Stokes *s)
610: {
611: PetscScalar scale = PetscSqrtReal((double)s->nx*s->ny);
612: PetscReal val;
613: Vec y0, y1;
616: /* error y-x */
617: VecAXPY(s->y, -1.0, s->x);
619: /* error in velocity */
620: VecGetSubVector(s->y, s->isg[0], &y0);
621: VecNorm(y0, NORM_2, &val);
622: PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(PetscRealPart(val/scale)));
623: VecRestoreSubVector(s->y, s->isg[0], &y0);
625: /* error in pressure */
626: VecGetSubVector(s->y, s->isg[1], &y1);
627: VecNorm(y1, NORM_2, &val);
628: PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(PetscRealPart(val/scale)));
629: VecRestoreSubVector(s->y, s->isg[1], &y1);
631: /* total error */
632: VecNorm(s->y, NORM_2, &val);
633: PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)PetscRealPart((val/scale)));
634: return 0;
635: }
637: int main(int argc, char **argv)
638: {
639: Stokes s;
640: KSP ksp;
642: PetscInitialize(&argc, &argv, NULL,help);
643: s.nx = 4;
644: s.ny = 6;
645: PetscOptionsGetInt(NULL,NULL, "-nx", &s.nx, NULL);
646: PetscOptionsGetInt(NULL,NULL, "-ny", &s.ny, NULL);
647: s.hx = 2.0/s.nx;
648: s.hy = 1.0/s.ny;
649: s.matsymmetric = PETSC_FALSE;
650: PetscOptionsGetBool(NULL,NULL, "-mat_set_symmetric", &s.matsymmetric,NULL);
651: s.userPC = s.userKSP = PETSC_FALSE;
652: PetscOptionsHasName(NULL,NULL, "-user_pc", &s.userPC);
653: PetscOptionsHasName(NULL,NULL, "-user_ksp", &s.userKSP);
655: StokesSetupMatrix(&s);
656: StokesSetupIndexSets(&s);
657: StokesSetupVectors(&s);
659: KSPCreate(PETSC_COMM_WORLD, &ksp);
660: KSPSetOperators(ksp, s.A, s.A);
661: KSPSetFromOptions(ksp);
662: StokesSetupPC(&s, ksp);
663: KSPSolve(ksp, s.b, s.x);
665: /* don't trust, verify! */
666: StokesCalcResidual(&s);
667: StokesCalcError(&s);
668: StokesWriteSolution(&s);
670: KSPDestroy(&ksp);
671: MatDestroy(&s.subA[0]);
672: MatDestroy(&s.subA[1]);
673: MatDestroy(&s.subA[2]);
674: MatDestroy(&s.subA[3]);
675: MatDestroy(&s.A);
676: VecDestroy(&s.x);
677: VecDestroy(&s.b);
678: VecDestroy(&s.y);
679: MatDestroy(&s.myS);
680: PetscFinalize();
681: return 0;
682: }
684: /*TEST
686: test:
687: nsize: 2
688: 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
690: test:
691: suffix: 2
692: nsize: 2
693: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
695: test:
696: suffix: 3
697: nsize: 2
698: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
700: test:
701: suffix: 4
702: nsize: 2
703: 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
705: test:
706: suffix: 4_pcksp
707: nsize: 2
708: 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
710: test:
711: suffix: 5
712: nsize: 2
713: 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
715: test:
716: suffix: 6
717: nsize: 2
718: 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
720: test:
721: suffix: 7
722: nsize: 2
723: 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
725: TEST*/