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;
105: KSPGetPC(ksp, &pc);
106: PCFieldSplitSetIS(pc, "0", s->isg[0]);
107: PCFieldSplitSetIS(pc, "1", s->isg[1]);
108: if (s->userPC) {
109: PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS);
110: }
111: if (s->userKSP) {
112: PCSetUp(pc);
113: PCFieldSplitGetSubKSP(pc, &n, &subksp);
114: KSPSetOperators(subksp[1], s->myS, s->myS);
115: PetscFree(subksp);
116: }
117: return(0);
118: }
120: PetscErrorCode StokesWriteSolution(Stokes *s)
121: {
122: PetscMPIInt size;
123: PetscInt n,i,j;
124: const PetscScalar *array;
125: PetscErrorCode ierr;
128: /* write data (*warning* only works sequential) */
129: MPI_Comm_size(MPI_COMM_WORLD,&size);
130: if (size == 1) {
131: PetscViewer viewer;
132: VecGetArrayRead(s->x, &array);
133: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer);
134: PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n");
135: for (j = 0; j < s->ny; j++) {
136: for (i = 0; i < s->nx; i++) {
137: n = j*s->nx+i;
138: 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]));
139: }
140: }
141: VecRestoreArrayRead(s->x, &array);
142: PetscViewerDestroy(&viewer);
143: }
144: return(0);
145: }
147: PetscErrorCode StokesSetupIndexSets(Stokes *s)
148: {
152: /* the two index sets */
153: MatNestGetISs(s->A, s->isg, NULL);
154: /* ISView(isg[0],PETSC_VIEWER_STDOUT_WORLD); */
155: /* ISView(isg[1],PETSC_VIEWER_STDOUT_WORLD); */
156: return(0);
157: }
159: PetscErrorCode StokesSetupVectors(Stokes *s)
160: {
164: /* solution vector x */
165: VecCreate(PETSC_COMM_WORLD, &s->x);
166: VecSetSizes(s->x, PETSC_DECIDE, 3*s->nx*s->ny);
167: VecSetType(s->x, VECMPI);
169: /* exact solution y */
170: VecDuplicate(s->x, &s->y);
171: StokesExactSolution(s);
173: /* rhs vector b */
174: VecDuplicate(s->x, &s->b);
175: StokesRhs(s);
176: return(0);
177: }
179: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
180: {
181: PetscInt n;
184: /* cell number n=j*nx+i has position (i,j) in grid */
185: n = row%(s->nx*s->ny);
186: *i = n%s->nx;
187: *j = (n-(*i))/s->nx;
188: return(0);
189: }
191: PetscErrorCode StokesExactSolution(Stokes *s)
192: {
193: PetscInt row, start, end, i, j;
194: PetscScalar val;
195: Vec y0,y1;
199: /* velocity part */
200: VecGetSubVector(s->y, s->isg[0], &y0);
201: VecGetOwnershipRange(y0, &start, &end);
202: for (row = start; row < end; row++) {
203: StokesGetPosition(s, row,&i,&j);
204: if (row < s->nx*s->ny) {
205: val = StokesExactVelocityX(j*s->hy+s->hy/2);
206: } else {
207: val = 0;
208: }
209: VecSetValue(y0, row, val, INSERT_VALUES);
210: }
211: VecRestoreSubVector(s->y, s->isg[0], &y0);
213: /* pressure part */
214: VecGetSubVector(s->y, s->isg[1], &y1);
215: VecGetOwnershipRange(y1, &start, &end);
216: for (row = start; row < end; row++) {
217: StokesGetPosition(s, row, &i, &j);
218: val = StokesExactPressure(i*s->hx+s->hx/2);
219: VecSetValue(y1, row, val, INSERT_VALUES);
220: }
221: VecRestoreSubVector(s->y, s->isg[1], &y1);
222: return(0);
223: }
225: PetscErrorCode StokesRhs(Stokes *s)
226: {
227: PetscInt row, start, end, i, j;
228: PetscScalar val;
229: Vec b0,b1;
233: /* velocity part */
234: VecGetSubVector(s->b, s->isg[0], &b0);
235: VecGetOwnershipRange(b0, &start, &end);
236: for (row = start; row < end; row++) {
237: StokesGetPosition(s, row, &i, &j);
238: if (row < s->nx*s->ny) {
239: StokesRhsMomX(s, i, j, &val);
240: } else {
241: StokesRhsMomY(s, i, j, &val);
242: }
243: VecSetValue(b0, row, val, INSERT_VALUES);
244: }
245: VecRestoreSubVector(s->b, s->isg[0], &b0);
247: /* pressure part */
248: VecGetSubVector(s->b, s->isg[1], &b1);
249: VecGetOwnershipRange(b1, &start, &end);
250: for (row = start; row < end; row++) {
251: StokesGetPosition(s, row, &i, &j);
252: StokesRhsMass(s, i, j, &val);
253: if (s->matsymmetric) {
254: val = -1.0*val;
255: }
256: VecSetValue(b1, row, val, INSERT_VALUES);
257: }
258: VecRestoreSubVector(s->b, s->isg[1], &b1);
259: return(0);
260: }
262: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
263: {
264: PetscInt row, start, end, sz, i, j;
265: PetscInt cols[5];
266: PetscScalar vals[5];
270: /* A[0] is 2N-by-2N */
271: MatCreate(PETSC_COMM_WORLD,&s->subA[0]);
272: MatSetOptionsPrefix(s->subA[0],"a00_");
273: MatSetSizes(s->subA[0],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,2*s->nx*s->ny);
274: MatSetType(s->subA[0],MATMPIAIJ);
275: MatMPIAIJSetPreallocation(s->subA[0],5,NULL,5,NULL);
276: MatGetOwnershipRange(s->subA[0], &start, &end);
278: for (row = start; row < end; row++) {
279: StokesGetPosition(s, row, &i, &j);
280: /* first part: rows 0 to (nx*ny-1) */
281: StokesStencilLaplacian(s, i, j, &sz, cols, vals);
282: /* second part: rows (nx*ny) to (2*nx*ny-1) */
283: if (row >= s->nx*s->ny) {
284: for (i = 0; i < sz; i++) cols[i] += s->nx*s->ny;
285: }
286: for (i = 0; i < sz; i++) vals[i] = -1.0*vals[i]; /* dynamic viscosity coef mu=-1 */
287: MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES);
288: }
289: MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY);
290: MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY);
291: return(0);
292: }
294: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
295: {
296: PetscInt row, start, end, sz, i, j;
297: PetscInt cols[5];
298: PetscScalar vals[5];
302: /* A[1] is 2N-by-N */
303: MatCreate(PETSC_COMM_WORLD, &s->subA[1]);
304: MatSetOptionsPrefix(s->subA[1],"a01_");
305: MatSetSizes(s->subA[1],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,s->nx*s->ny);
306: MatSetType(s->subA[1],MATMPIAIJ);
307: MatMPIAIJSetPreallocation(s->subA[1],5,NULL,5,NULL);
308: MatGetOwnershipRange(s->subA[1],&start,&end);
310: MatSetOption(s->subA[1],MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
312: for (row = start; row < end; row++) {
313: StokesGetPosition(s, row, &i, &j);
314: /* first part: rows 0 to (nx*ny-1) */
315: if (row < s->nx*s->ny) {
316: StokesStencilGradientX(s, i, j, &sz, cols, vals);
317: } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */
318: StokesStencilGradientY(s, i, j, &sz, cols, vals);
319: }
320: MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES);
321: }
322: MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY);
323: MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY);
324: return(0);
325: }
327: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
328: {
332: /* A[2] is minus transpose of A[1] */
333: MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]);
334: if (!s->matsymmetric) {
335: MatScale(s->subA[2], -1.0);
336: }
337: MatSetOptionsPrefix(s->subA[2], "a10_");
338: return(0);
339: }
341: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
342: {
346: /* A[3] is N-by-N null matrix */
347: MatCreate(PETSC_COMM_WORLD, &s->subA[3]);
348: MatSetOptionsPrefix(s->subA[3], "a11_");
349: MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx*s->ny, s->nx*s->ny);
350: MatSetType(s->subA[3], MATMPIAIJ);
351: MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL);
352: MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY);
353: MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY);
354: return(0);
355: }
357: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
358: {
359: Vec diag;
363: /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */
364: /* note: A11 is zero */
365: /* note: in real life this matrix would be build directly, */
366: /* i.e. without MatMatMult */
368: /* inverse of diagonal of A00 */
369: VecCreate(PETSC_COMM_WORLD,&diag);
370: VecSetSizes(diag,PETSC_DECIDE,2*s->nx*s->ny);
371: VecSetType(diag,VECMPI);
372: MatGetDiagonal(s->subA[0],diag);
373: VecReciprocal(diag);
375: /* compute: - A10 inv(DIAGFORM(A00)) A01 */
376: MatDiagonalScale(s->subA[1],diag,NULL); /* (*warning* overwrites subA[1]) */
377: MatMatMult(s->subA[2],s->subA[1],MAT_INITIAL_MATRIX,PETSC_DEFAULT,&s->myS);
378: MatScale(s->myS,-1.0);
380: /* restore A10 */
381: MatGetDiagonal(s->subA[0],diag);
382: MatDiagonalScale(s->subA[1],diag,NULL);
383: VecDestroy(&diag);
384: return(0);
385: }
387: PetscErrorCode StokesSetupMatrix(Stokes *s)
388: {
392: StokesSetupMatBlock00(s);
393: StokesSetupMatBlock01(s);
394: StokesSetupMatBlock10(s);
395: StokesSetupMatBlock11(s);
396: MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A);
397: StokesSetupApproxSchur(s);
398: return(0);
399: }
401: PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
402: {
403: PetscInt p =j*s->nx+i, w=p-1, e=p+1, s2=p-s->nx, n=p+s->nx;
404: PetscScalar ae=s->hy/s->hx, aeb=0;
405: PetscScalar aw=s->hy/s->hx, awb=s->hy/(s->hx/2);
406: PetscScalar as=s->hx/s->hy, asb=s->hx/(s->hy/2);
407: PetscScalar an=s->hx/s->hy, anb=s->hx/(s->hy/2);
410: if (i==0 && j==0) { /* south-west corner */
411: *sz =3;
412: cols[0]=p; vals[0]=-(ae+awb+asb+an);
413: cols[1]=e; vals[1]=ae;
414: cols[2]=n; vals[2]=an;
415: } else if (i==0 && j==s->ny-1) { /* north-west corner */
416: *sz =3;
417: cols[0]=s2; vals[0]=as;
418: cols[1]=p; vals[1]=-(ae+awb+as+anb);
419: cols[2]=e; vals[2]=ae;
420: } else if (i==s->nx-1 && j==0) { /* south-east corner */
421: *sz =3;
422: cols[0]=w; vals[0]=aw;
423: cols[1]=p; vals[1]=-(aeb+aw+asb+an);
424: cols[2]=n; vals[2]=an;
425: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
426: *sz =3;
427: cols[0]=s2; vals[0]=as;
428: cols[1]=w; vals[1]=aw;
429: cols[2]=p; vals[2]=-(aeb+aw+as+anb);
430: } else if (i==0) { /* west boundary */
431: *sz =4;
432: cols[0]=s2; vals[0]=as;
433: cols[1]=p; vals[1]=-(ae+awb+as+an);
434: cols[2]=e; vals[2]=ae;
435: cols[3]=n; vals[3]=an;
436: } else if (i==s->nx-1) { /* east boundary */
437: *sz =4;
438: cols[0]=s2; vals[0]=as;
439: cols[1]=w; vals[1]=aw;
440: cols[2]=p; vals[2]=-(aeb+aw+as+an);
441: cols[3]=n; vals[3]=an;
442: } else if (j==0) { /* south boundary */
443: *sz =4;
444: cols[0]=w; vals[0]=aw;
445: cols[1]=p; vals[1]=-(ae+aw+asb+an);
446: cols[2]=e; vals[2]=ae;
447: cols[3]=n; vals[3]=an;
448: } else if (j==s->ny-1) { /* north boundary */
449: *sz =4;
450: cols[0]=s2; vals[0]=as;
451: cols[1]=w; vals[1]=aw;
452: cols[2]=p; vals[2]=-(ae+aw+as+anb);
453: cols[3]=e; vals[3]=ae;
454: } else { /* interior */
455: *sz =5;
456: cols[0]=s2; vals[0]=as;
457: cols[1]=w; vals[1]=aw;
458: cols[2]=p; vals[2]=-(ae+aw+as+an);
459: cols[3]=e; vals[3]=ae;
460: cols[4]=n; vals[4]=an;
461: }
462: return(0);
463: }
465: PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
466: {
467: PetscInt p =j*s->nx+i, w=p-1, e=p+1;
468: PetscScalar ae= s->hy/2, aeb=s->hy;
469: PetscScalar aw=-s->hy/2, awb=0;
472: if (i==0 && j==0) { /* south-west corner */
473: *sz =2;
474: cols[0]=p; vals[0]=-(ae+awb);
475: cols[1]=e; vals[1]=ae;
476: } else if (i==0 && j==s->ny-1) { /* north-west corner */
477: *sz =2;
478: cols[0]=p; vals[0]=-(ae+awb);
479: cols[1]=e; vals[1]=ae;
480: } else if (i==s->nx-1 && j==0) { /* south-east corner */
481: *sz =2;
482: cols[0]=w; vals[0]=aw;
483: cols[1]=p; vals[1]=-(aeb+aw);
484: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
485: *sz =2;
486: cols[0]=w; vals[0]=aw;
487: cols[1]=p; vals[1]=-(aeb+aw);
488: } else if (i==0) { /* west boundary */
489: *sz =2;
490: cols[0]=p; vals[0]=-(ae+awb);
491: cols[1]=e; vals[1]=ae;
492: } else if (i==s->nx-1) { /* east boundary */
493: *sz =2;
494: cols[0]=w; vals[0]=aw;
495: cols[1]=p; vals[1]=-(aeb+aw);
496: } else if (j==0) { /* south boundary */
497: *sz =3;
498: cols[0]=w; vals[0]=aw;
499: cols[1]=p; vals[1]=-(ae+aw);
500: cols[2]=e; vals[2]=ae;
501: } else if (j==s->ny-1) { /* north boundary */
502: *sz =3;
503: cols[0]=w; vals[0]=aw;
504: cols[1]=p; vals[1]=-(ae+aw);
505: cols[2]=e; vals[2]=ae;
506: } else { /* interior */
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: }
512: return(0);
513: }
515: PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
516: {
517: PetscInt p =j*s->nx+i, s2=p-s->nx, n=p+s->nx;
518: PetscScalar as=-s->hx/2, asb=0;
519: PetscScalar an= s->hx/2, anb=0;
522: if (i==0 && j==0) { /* south-west corner */
523: *sz =2;
524: cols[0]=p; vals[0]=-(asb+an);
525: cols[1]=n; vals[1]=an;
526: } else if (i==0 && j==s->ny-1) { /* north-west corner */
527: *sz =2;
528: cols[0]=s2; vals[0]=as;
529: cols[1]=p; vals[1]=-(as+anb);
530: } else if (i==s->nx-1 && j==0) { /* south-east corner */
531: *sz =2;
532: cols[0]=p; vals[0]=-(asb+an);
533: cols[1]=n; vals[1]=an;
534: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
535: *sz =2;
536: cols[0]=s2; vals[0]=as;
537: cols[1]=p; vals[1]=-(as+anb);
538: } else if (i==0) { /* west boundary */
539: *sz =3;
540: cols[0]=s2; vals[0]=as;
541: cols[1]=p; vals[1]=-(as+an);
542: cols[2]=n; vals[2]=an;
543: } else if (i==s->nx-1) { /* east boundary */
544: *sz =3;
545: cols[0]=s2; vals[0]=as;
546: cols[1]=p; vals[1]=-(as+an);
547: cols[2]=n; vals[2]=an;
548: } else if (j==0) { /* south boundary */
549: *sz =2;
550: cols[0]=p; vals[0]=-(asb+an);
551: cols[1]=n; vals[1]=an;
552: } else if (j==s->ny-1) { /* north boundary */
553: *sz =2;
554: cols[0]=s2; vals[0]=as;
555: cols[1]=p; vals[1]=-(as+anb);
556: } else { /* interior */
557: *sz =3;
558: cols[0]=s2; vals[0]=as;
559: cols[1]=p; vals[1]=-(as+an);
560: cols[2]=n; vals[2]=an;
561: }
562: return(0);
563: }
565: PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
566: {
567: PetscScalar y = j*s->hy+s->hy/2;
568: PetscScalar awb = s->hy/(s->hx/2);
571: if (i == 0) { /* west boundary */
572: *val = awb*StokesExactVelocityX(y);
573: } else {
574: *val = 0.0;
575: }
576: return(0);
577: }
579: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
580: {
582: *val = 0.0;
583: return(0);
584: }
586: PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
587: {
588: PetscScalar y = j*s->hy+s->hy/2;
589: PetscScalar aeb = s->hy;
592: if (i == 0) { /* west boundary */
593: *val = aeb*StokesExactVelocityX(y);
594: } else {
595: *val = 0.0;
596: }
597: return(0);
598: }
600: PetscErrorCode StokesCalcResidual(Stokes *s)
601: {
602: PetscReal val;
603: Vec b0, b1;
607: /* residual Ax-b (*warning* overwrites b) */
608: VecScale(s->b, -1.0);
609: MatMultAdd(s->A, s->x, s->b, s->b);
611: /* residual velocity */
612: VecGetSubVector(s->b, s->isg[0], &b0);
613: VecNorm(b0, NORM_2, &val);
614: PetscPrintf(PETSC_COMM_WORLD," residual u = %g\n",(double)val);
615: VecRestoreSubVector(s->b, s->isg[0], &b0);
617: /* residual pressure */
618: VecGetSubVector(s->b, s->isg[1], &b1);
619: VecNorm(b1, NORM_2, &val);
620: PetscPrintf(PETSC_COMM_WORLD," residual p = %g\n",(double)val);
621: VecRestoreSubVector(s->b, s->isg[1], &b1);
623: /* total residual */
624: VecNorm(s->b, NORM_2, &val);
625: PetscPrintf(PETSC_COMM_WORLD," residual [u,p] = %g\n", (double)val);
626: return(0);
627: }
629: PetscErrorCode StokesCalcError(Stokes *s)
630: {
631: PetscScalar scale = PetscSqrtReal((double)s->nx*s->ny);
632: PetscReal val;
633: Vec y0, y1;
637: /* error y-x */
638: VecAXPY(s->y, -1.0, s->x);
640: /* error in velocity */
641: VecGetSubVector(s->y, s->isg[0], &y0);
642: VecNorm(y0, NORM_2, &val);
643: PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(PetscRealPart(val/scale)));
644: VecRestoreSubVector(s->y, s->isg[0], &y0);
646: /* error in pressure */
647: VecGetSubVector(s->y, s->isg[1], &y1);
648: VecNorm(y1, NORM_2, &val);
649: PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(PetscRealPart(val/scale)));
650: VecRestoreSubVector(s->y, s->isg[1], &y1);
652: /* total error */
653: VecNorm(s->y, NORM_2, &val);
654: PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)PetscRealPart((val/scale)));
655: return(0);
656: }
658: int main(int argc, char **argv)
659: {
660: Stokes s;
661: KSP ksp;
664: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
665: s.nx = 4;
666: s.ny = 6;
667: PetscOptionsGetInt(NULL,NULL, "-nx", &s.nx, NULL);
668: PetscOptionsGetInt(NULL,NULL, "-ny", &s.ny, NULL);
669: s.hx = 2.0/s.nx;
670: s.hy = 1.0/s.ny;
671: s.matsymmetric = PETSC_FALSE;
672: PetscOptionsGetBool(NULL,NULL, "-mat_set_symmetric", &s.matsymmetric,NULL);
673: s.userPC = s.userKSP = PETSC_FALSE;
674: PetscOptionsHasName(NULL,NULL, "-user_pc", &s.userPC);
675: PetscOptionsHasName(NULL,NULL, "-user_ksp", &s.userKSP);
677: StokesSetupMatrix(&s);
678: StokesSetupIndexSets(&s);
679: StokesSetupVectors(&s);
681: KSPCreate(PETSC_COMM_WORLD, &ksp);
682: KSPSetOperators(ksp, s.A, s.A);
683: KSPSetFromOptions(ksp);
684: StokesSetupPC(&s, ksp);
685: KSPSolve(ksp, s.b, s.x);
687: /* don't trust, verify! */
688: StokesCalcResidual(&s);
689: StokesCalcError(&s);
690: StokesWriteSolution(&s);
692: KSPDestroy(&ksp);
693: MatDestroy(&s.subA[0]);
694: MatDestroy(&s.subA[1]);
695: MatDestroy(&s.subA[2]);
696: MatDestroy(&s.subA[3]);
697: MatDestroy(&s.A);
698: VecDestroy(&s.x);
699: VecDestroy(&s.b);
700: VecDestroy(&s.y);
701: MatDestroy(&s.myS);
702: PetscFinalize();
703: return ierr;
704: }
706: /*TEST
708: test:
709: nsize: 2
710: 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
712: test:
713: suffix: 2
714: nsize: 2
715: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
717: test:
718: suffix: 3
719: nsize: 2
720: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
722: test:
723: suffix: 4
724: nsize: 2
725: 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
727: test:
728: suffix: 4_pcksp
729: nsize: 2
730: 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
732: test:
733: suffix: 5
734: nsize: 2
735: 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
737: test:
738: suffix: 6
739: nsize: 2
740: 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
742: test:
743: suffix: 7
744: nsize: 2
745: 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
747: TEST*/