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) PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS);
108: if (s->userKSP) {
109: PCSetUp(pc);
110: PCFieldSplitGetSubKSP(pc, &n, &subksp);
111: KSPSetOperators(subksp[1], s->myS, s->myS);
112: PetscFree(subksp);
113: }
114: return 0;
115: }
117: PetscErrorCode StokesWriteSolution(Stokes *s)
118: {
119: PetscMPIInt size;
120: PetscInt n, i, j;
121: const PetscScalar *array;
124: /* write data (*warning* only works sequential) */
125: MPI_Comm_size(MPI_COMM_WORLD, &size);
126: if (size == 1) {
127: PetscViewer viewer;
128: VecGetArrayRead(s->x, &array);
129: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer);
130: PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n");
131: for (j = 0; j < s->ny; j++) {
132: for (i = 0; i < s->nx; i++) {
133: n = j * s->nx + i;
134: 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]));
135: }
136: }
137: VecRestoreArrayRead(s->x, &array);
138: PetscViewerDestroy(&viewer);
139: }
140: return 0;
141: }
143: PetscErrorCode StokesSetupIndexSets(Stokes *s)
144: {
146: /* the two index sets */
147: MatNestGetISs(s->A, s->isg, NULL);
148: return 0;
149: }
151: PetscErrorCode StokesSetupVectors(Stokes *s)
152: {
154: /* solution vector x */
155: VecCreate(PETSC_COMM_WORLD, &s->x);
156: VecSetSizes(s->x, PETSC_DECIDE, 3 * s->nx * s->ny);
157: VecSetType(s->x, VECMPI);
159: /* exact solution y */
160: VecDuplicate(s->x, &s->y);
161: StokesExactSolution(s);
163: /* rhs vector b */
164: VecDuplicate(s->x, &s->b);
165: StokesRhs(s);
166: return 0;
167: }
169: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
170: {
171: PetscInt n;
174: /* cell number n=j*nx+i has position (i,j) in grid */
175: n = row % (s->nx * s->ny);
176: *i = n % s->nx;
177: *j = (n - (*i)) / s->nx;
178: return 0;
179: }
181: PetscErrorCode StokesExactSolution(Stokes *s)
182: {
183: PetscInt row, start, end, i, j;
184: PetscScalar val;
185: Vec y0, y1;
188: /* velocity part */
189: VecGetSubVector(s->y, s->isg[0], &y0);
190: VecGetOwnershipRange(y0, &start, &end);
191: for (row = start; row < end; row++) {
192: StokesGetPosition(s, row, &i, &j);
193: if (row < s->nx * s->ny) {
194: val = StokesExactVelocityX(j * s->hy + s->hy / 2);
195: } else {
196: val = 0;
197: }
198: VecSetValue(y0, row, val, INSERT_VALUES);
199: }
200: VecRestoreSubVector(s->y, s->isg[0], &y0);
202: /* pressure part */
203: VecGetSubVector(s->y, s->isg[1], &y1);
204: VecGetOwnershipRange(y1, &start, &end);
205: for (row = start; row < end; row++) {
206: StokesGetPosition(s, row, &i, &j);
207: val = StokesExactPressure(i * s->hx + s->hx / 2);
208: VecSetValue(y1, row, val, INSERT_VALUES);
209: }
210: VecRestoreSubVector(s->y, s->isg[1], &y1);
211: return 0;
212: }
214: PetscErrorCode StokesRhs(Stokes *s)
215: {
216: PetscInt row, start, end, i, j;
217: PetscScalar val;
218: Vec b0, b1;
221: /* velocity part */
222: VecGetSubVector(s->b, s->isg[0], &b0);
223: VecGetOwnershipRange(b0, &start, &end);
224: for (row = start; row < end; row++) {
225: StokesGetPosition(s, row, &i, &j);
226: if (row < s->nx * s->ny) {
227: StokesRhsMomX(s, i, j, &val);
228: } else {
229: StokesRhsMomY(s, i, j, &val);
230: }
231: VecSetValue(b0, row, val, INSERT_VALUES);
232: }
233: VecRestoreSubVector(s->b, s->isg[0], &b0);
235: /* pressure part */
236: VecGetSubVector(s->b, s->isg[1], &b1);
237: VecGetOwnershipRange(b1, &start, &end);
238: for (row = start; row < end; row++) {
239: StokesGetPosition(s, row, &i, &j);
240: StokesRhsMass(s, i, j, &val);
241: if (s->matsymmetric) val = -1.0 * val;
242: VecSetValue(b1, row, val, INSERT_VALUES);
243: }
244: VecRestoreSubVector(s->b, s->isg[1], &b1);
245: return 0;
246: }
248: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
249: {
250: PetscInt row, start, end, sz, i, j;
251: PetscInt cols[5];
252: PetscScalar vals[5];
255: /* A[0] is 2N-by-2N */
256: MatCreate(PETSC_COMM_WORLD, &s->subA[0]);
257: MatSetOptionsPrefix(s->subA[0], "a00_");
258: MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny);
259: MatSetType(s->subA[0], MATMPIAIJ);
260: MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL);
261: MatGetOwnershipRange(s->subA[0], &start, &end);
263: for (row = start; row < end; row++) {
264: StokesGetPosition(s, row, &i, &j);
265: /* first part: rows 0 to (nx*ny-1) */
266: StokesStencilLaplacian(s, i, j, &sz, cols, vals);
267: /* second part: rows (nx*ny) to (2*nx*ny-1) */
268: if (row >= s->nx * s->ny) {
269: for (i = 0; i < sz; i++) cols[i] += s->nx * s->ny;
270: }
271: for (i = 0; i < sz; i++) vals[i] = -1.0 * vals[i]; /* dynamic viscosity coef mu=-1 */
272: MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES);
273: }
274: MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY);
275: MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY);
276: return 0;
277: }
279: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
280: {
281: PetscInt row, start, end, sz, i, j;
282: PetscInt cols[5];
283: PetscScalar vals[5];
286: /* A[1] is 2N-by-N */
287: MatCreate(PETSC_COMM_WORLD, &s->subA[1]);
288: MatSetOptionsPrefix(s->subA[1], "a01_");
289: MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny);
290: MatSetType(s->subA[1], MATMPIAIJ);
291: MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL);
292: MatGetOwnershipRange(s->subA[1], &start, &end);
294: MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
296: for (row = start; row < end; row++) {
297: StokesGetPosition(s, row, &i, &j);
298: /* first part: rows 0 to (nx*ny-1) */
299: if (row < s->nx * s->ny) {
300: StokesStencilGradientX(s, i, j, &sz, cols, vals);
301: } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */
302: StokesStencilGradientY(s, i, j, &sz, cols, vals);
303: }
304: MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES);
305: }
306: MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY);
307: MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY);
308: return 0;
309: }
311: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
312: {
314: /* A[2] is minus transpose of A[1] */
315: MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]);
316: if (!s->matsymmetric) MatScale(s->subA[2], -1.0);
317: MatSetOptionsPrefix(s->subA[2], "a10_");
318: return 0;
319: }
321: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
322: {
324: /* A[3] is N-by-N null matrix */
325: MatCreate(PETSC_COMM_WORLD, &s->subA[3]);
326: MatSetOptionsPrefix(s->subA[3], "a11_");
327: MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny);
328: MatSetType(s->subA[3], MATMPIAIJ);
329: MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL);
330: MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY);
331: MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY);
332: return 0;
333: }
335: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
336: {
337: Vec diag;
340: /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */
341: /* note: A11 is zero */
342: /* note: in real life this matrix would be build directly, */
343: /* i.e. without MatMatMult */
345: /* inverse of diagonal of A00 */
346: VecCreate(PETSC_COMM_WORLD, &diag);
347: VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny);
348: VecSetType(diag, VECMPI);
349: MatGetDiagonal(s->subA[0], diag);
350: VecReciprocal(diag);
352: /* compute: - A10 inv(DIAGFORM(A00)) A01 */
353: MatDiagonalScale(s->subA[1], diag, NULL); /* (*warning* overwrites subA[1]) */
354: MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DEFAULT, &s->myS);
355: MatScale(s->myS, -1.0);
357: /* restore A10 */
358: MatGetDiagonal(s->subA[0], diag);
359: MatDiagonalScale(s->subA[1], diag, NULL);
360: VecDestroy(&diag);
361: return 0;
362: }
364: PetscErrorCode StokesSetupMatrix(Stokes *s)
365: {
367: StokesSetupMatBlock00(s);
368: StokesSetupMatBlock01(s);
369: StokesSetupMatBlock10(s);
370: StokesSetupMatBlock11(s);
371: MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A);
372: StokesSetupApproxSchur(s);
373: return 0;
374: }
376: PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
377: {
378: PetscInt p = j * s->nx + i, w = p - 1, e = p + 1, s2 = p - s->nx, n = p + s->nx;
379: PetscScalar ae = s->hy / s->hx, aeb = 0;
380: PetscScalar aw = s->hy / s->hx, awb = s->hy / (s->hx / 2);
381: PetscScalar as = s->hx / s->hy, asb = s->hx / (s->hy / 2);
382: PetscScalar an = s->hx / s->hy, anb = s->hx / (s->hy / 2);
385: if (i == 0 && j == 0) { /* south-west corner */
386: *sz = 3;
387: cols[0] = p;
388: vals[0] = -(ae + awb + asb + an);
389: cols[1] = e;
390: vals[1] = ae;
391: cols[2] = n;
392: vals[2] = an;
393: } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
394: *sz = 3;
395: cols[0] = s2;
396: vals[0] = as;
397: cols[1] = p;
398: vals[1] = -(ae + awb + as + anb);
399: cols[2] = e;
400: vals[2] = ae;
401: } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
402: *sz = 3;
403: cols[0] = w;
404: vals[0] = aw;
405: cols[1] = p;
406: vals[1] = -(aeb + aw + asb + an);
407: cols[2] = n;
408: vals[2] = an;
409: } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
410: *sz = 3;
411: cols[0] = s2;
412: vals[0] = as;
413: cols[1] = w;
414: vals[1] = aw;
415: cols[2] = p;
416: vals[2] = -(aeb + aw + as + anb);
417: } else if (i == 0) { /* west boundary */
418: *sz = 4;
419: cols[0] = s2;
420: vals[0] = as;
421: cols[1] = p;
422: vals[1] = -(ae + awb + as + an);
423: cols[2] = e;
424: vals[2] = ae;
425: cols[3] = n;
426: vals[3] = an;
427: } else if (i == s->nx - 1) { /* east boundary */
428: *sz = 4;
429: cols[0] = s2;
430: vals[0] = as;
431: cols[1] = w;
432: vals[1] = aw;
433: cols[2] = p;
434: vals[2] = -(aeb + aw + as + an);
435: cols[3] = n;
436: vals[3] = an;
437: } else if (j == 0) { /* south boundary */
438: *sz = 4;
439: cols[0] = w;
440: vals[0] = aw;
441: cols[1] = p;
442: vals[1] = -(ae + aw + asb + an);
443: cols[2] = e;
444: vals[2] = ae;
445: cols[3] = n;
446: vals[3] = an;
447: } else if (j == s->ny - 1) { /* north boundary */
448: *sz = 4;
449: cols[0] = s2;
450: vals[0] = as;
451: cols[1] = w;
452: vals[1] = aw;
453: cols[2] = p;
454: vals[2] = -(ae + aw + as + anb);
455: cols[3] = e;
456: vals[3] = ae;
457: } else { /* interior */
458: *sz = 5;
459: cols[0] = s2;
460: vals[0] = as;
461: cols[1] = w;
462: vals[1] = aw;
463: cols[2] = p;
464: vals[2] = -(ae + aw + as + an);
465: cols[3] = e;
466: vals[3] = ae;
467: cols[4] = n;
468: vals[4] = an;
469: }
470: return 0;
471: }
473: PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
474: {
475: PetscInt p = j * s->nx + i, w = p - 1, e = p + 1;
476: PetscScalar ae = s->hy / 2, aeb = s->hy;
477: PetscScalar aw = -s->hy / 2, awb = 0;
480: if (i == 0 && j == 0) { /* south-west corner */
481: *sz = 2;
482: cols[0] = p;
483: vals[0] = -(ae + awb);
484: cols[1] = e;
485: vals[1] = ae;
486: } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
487: *sz = 2;
488: cols[0] = p;
489: vals[0] = -(ae + awb);
490: cols[1] = e;
491: vals[1] = ae;
492: } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
493: *sz = 2;
494: cols[0] = w;
495: vals[0] = aw;
496: cols[1] = p;
497: vals[1] = -(aeb + aw);
498: } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
499: *sz = 2;
500: cols[0] = w;
501: vals[0] = aw;
502: cols[1] = p;
503: vals[1] = -(aeb + aw);
504: } else if (i == 0) { /* west boundary */
505: *sz = 2;
506: cols[0] = p;
507: vals[0] = -(ae + awb);
508: cols[1] = e;
509: vals[1] = ae;
510: } else if (i == s->nx - 1) { /* east boundary */
511: *sz = 2;
512: cols[0] = w;
513: vals[0] = aw;
514: cols[1] = p;
515: vals[1] = -(aeb + aw);
516: } else if (j == 0) { /* south boundary */
517: *sz = 3;
518: cols[0] = w;
519: vals[0] = aw;
520: cols[1] = p;
521: vals[1] = -(ae + aw);
522: cols[2] = e;
523: vals[2] = ae;
524: } else if (j == s->ny - 1) { /* north boundary */
525: *sz = 3;
526: cols[0] = w;
527: vals[0] = aw;
528: cols[1] = p;
529: vals[1] = -(ae + aw);
530: cols[2] = e;
531: vals[2] = ae;
532: } else { /* interior */
533: *sz = 3;
534: cols[0] = w;
535: vals[0] = aw;
536: cols[1] = p;
537: vals[1] = -(ae + aw);
538: cols[2] = e;
539: vals[2] = ae;
540: }
541: return 0;
542: }
544: PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
545: {
546: PetscInt p = j * s->nx + i, s2 = p - s->nx, n = p + s->nx;
547: PetscScalar as = -s->hx / 2, asb = 0;
548: PetscScalar an = s->hx / 2, anb = 0;
551: if (i == 0 && j == 0) { /* south-west corner */
552: *sz = 2;
553: cols[0] = p;
554: vals[0] = -(asb + an);
555: cols[1] = n;
556: vals[1] = an;
557: } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
558: *sz = 2;
559: cols[0] = s2;
560: vals[0] = as;
561: cols[1] = p;
562: vals[1] = -(as + anb);
563: } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
564: *sz = 2;
565: cols[0] = p;
566: vals[0] = -(asb + an);
567: cols[1] = n;
568: vals[1] = an;
569: } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
570: *sz = 2;
571: cols[0] = s2;
572: vals[0] = as;
573: cols[1] = p;
574: vals[1] = -(as + anb);
575: } else if (i == 0) { /* west boundary */
576: *sz = 3;
577: cols[0] = s2;
578: vals[0] = as;
579: cols[1] = p;
580: vals[1] = -(as + an);
581: cols[2] = n;
582: vals[2] = an;
583: } else if (i == s->nx - 1) { /* east boundary */
584: *sz = 3;
585: cols[0] = s2;
586: vals[0] = as;
587: cols[1] = p;
588: vals[1] = -(as + an);
589: cols[2] = n;
590: vals[2] = an;
591: } else if (j == 0) { /* south boundary */
592: *sz = 2;
593: cols[0] = p;
594: vals[0] = -(asb + an);
595: cols[1] = n;
596: vals[1] = an;
597: } else if (j == s->ny - 1) { /* north boundary */
598: *sz = 2;
599: cols[0] = s2;
600: vals[0] = as;
601: cols[1] = p;
602: vals[1] = -(as + anb);
603: } else { /* interior */
604: *sz = 3;
605: cols[0] = s2;
606: vals[0] = as;
607: cols[1] = p;
608: vals[1] = -(as + an);
609: cols[2] = n;
610: vals[2] = an;
611: }
612: return 0;
613: }
615: PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
616: {
617: PetscScalar y = j * s->hy + s->hy / 2;
618: PetscScalar awb = s->hy / (s->hx / 2);
621: if (i == 0) { /* west boundary */
622: *val = awb * StokesExactVelocityX(y);
623: } else {
624: *val = 0.0;
625: }
626: return 0;
627: }
629: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
630: {
632: *val = 0.0;
633: return 0;
634: }
636: PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
637: {
638: PetscScalar y = j * s->hy + s->hy / 2;
639: PetscScalar aeb = s->hy;
642: if (i == 0) { /* west boundary */
643: *val = aeb * StokesExactVelocityX(y);
644: } else {
645: *val = 0.0;
646: }
647: return 0;
648: }
650: PetscErrorCode StokesCalcResidual(Stokes *s)
651: {
652: PetscReal val;
653: Vec b0, b1;
656: /* residual Ax-b (*warning* overwrites b) */
657: VecScale(s->b, -1.0);
658: MatMultAdd(s->A, s->x, s->b, s->b);
660: /* residual velocity */
661: VecGetSubVector(s->b, s->isg[0], &b0);
662: VecNorm(b0, NORM_2, &val);
663: PetscPrintf(PETSC_COMM_WORLD, " residual u = %g\n", (double)val);
664: VecRestoreSubVector(s->b, s->isg[0], &b0);
666: /* residual pressure */
667: VecGetSubVector(s->b, s->isg[1], &b1);
668: VecNorm(b1, NORM_2, &val);
669: PetscPrintf(PETSC_COMM_WORLD, " residual p = %g\n", (double)val);
670: VecRestoreSubVector(s->b, s->isg[1], &b1);
672: /* total residual */
673: VecNorm(s->b, NORM_2, &val);
674: PetscPrintf(PETSC_COMM_WORLD, " residual [u,p] = %g\n", (double)val);
675: return 0;
676: }
678: PetscErrorCode StokesCalcError(Stokes *s)
679: {
680: PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny);
681: PetscReal val;
682: Vec y0, y1;
685: /* error y-x */
686: VecAXPY(s->y, -1.0, s->x);
688: /* error in velocity */
689: VecGetSubVector(s->y, s->isg[0], &y0);
690: VecNorm(y0, NORM_2, &val);
691: PetscPrintf(PETSC_COMM_WORLD, " discretization error u = %g\n", (double)(PetscRealPart(val / scale)));
692: VecRestoreSubVector(s->y, s->isg[0], &y0);
694: /* error in pressure */
695: VecGetSubVector(s->y, s->isg[1], &y1);
696: VecNorm(y1, NORM_2, &val);
697: PetscPrintf(PETSC_COMM_WORLD, " discretization error p = %g\n", (double)(PetscRealPart(val / scale)));
698: VecRestoreSubVector(s->y, s->isg[1], &y1);
700: /* total error */
701: VecNorm(s->y, NORM_2, &val);
702: PetscPrintf(PETSC_COMM_WORLD, " discretization error [u,p] = %g\n", (double)PetscRealPart((val / scale)));
703: return 0;
704: }
706: int main(int argc, char **argv)
707: {
708: Stokes s;
709: KSP ksp;
712: PetscInitialize(&argc, &argv, NULL, help);
713: s.nx = 4;
714: s.ny = 6;
715: PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL);
716: PetscOptionsGetInt(NULL, NULL, "-ny", &s.ny, NULL);
717: s.hx = 2.0 / s.nx;
718: s.hy = 1.0 / s.ny;
719: s.matsymmetric = PETSC_FALSE;
720: PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL);
721: s.userPC = s.userKSP = PETSC_FALSE;
722: PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC);
723: PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP);
725: StokesSetupMatrix(&s);
726: StokesSetupIndexSets(&s);
727: StokesSetupVectors(&s);
729: KSPCreate(PETSC_COMM_WORLD, &ksp);
730: KSPSetOperators(ksp, s.A, s.A);
731: KSPSetFromOptions(ksp);
732: StokesSetupPC(&s, ksp);
733: KSPSolve(ksp, s.b, s.x);
735: /* don't trust, verify! */
736: StokesCalcResidual(&s);
737: StokesCalcError(&s);
738: StokesWriteSolution(&s);
740: KSPDestroy(&ksp);
741: MatDestroy(&s.subA[0]);
742: MatDestroy(&s.subA[1]);
743: MatDestroy(&s.subA[2]);
744: MatDestroy(&s.subA[3]);
745: MatDestroy(&s.A);
746: VecDestroy(&s.x);
747: VecDestroy(&s.b);
748: VecDestroy(&s.y);
749: MatDestroy(&s.myS);
750: PetscFinalize();
751: return 0;
752: }
754: /*TEST
756: test:
757: nsize: 2
758: 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
760: test:
761: suffix: 2
762: nsize: 2
763: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
765: test:
766: suffix: 3
767: nsize: 2
768: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
770: test:
771: suffix: 4
772: nsize: 2
773: 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
775: test:
776: suffix: 4_pcksp
777: nsize: 2
778: 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
780: test:
781: suffix: 5
782: nsize: 2
783: 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
785: test:
786: suffix: 6
787: nsize: 2
788: 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
790: test:
791: suffix: 7
792: nsize: 2
793: 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
795: TEST*/