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;
103: PetscFunctionBeginUser;
104: PetscCall(KSPGetPC(ksp, &pc));
105: PetscCall(PCFieldSplitSetIS(pc, "0", s->isg[0]));
106: PetscCall(PCFieldSplitSetIS(pc, "1", s->isg[1]));
107: if (s->userPC) PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS));
108: if (s->userKSP) {
109: PetscCall(PCSetUp(pc));
110: PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp));
111: PetscCall(KSPSetOperators(subksp[1], s->myS, s->myS));
112: PetscCall(PetscFree(subksp));
113: }
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: PetscErrorCode StokesWriteSolution(Stokes *s)
118: {
119: PetscMPIInt size;
120: PetscInt n, i, j;
121: const PetscScalar *array;
123: PetscFunctionBeginUser;
124: /* write data (*warning* only works sequential) */
125: PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
126: if (size == 1) {
127: PetscViewer viewer;
128: PetscCall(VecGetArrayRead(s->x, &array));
129: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer));
130: PetscCall(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: PetscCall(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: PetscCall(VecRestoreArrayRead(s->x, &array));
138: PetscCall(PetscViewerDestroy(&viewer));
139: }
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: PetscErrorCode StokesSetupIndexSets(Stokes *s)
144: {
145: PetscFunctionBeginUser;
146: /* the two index sets */
147: PetscCall(MatNestGetISs(s->A, s->isg, NULL));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: PetscErrorCode StokesSetupVectors(Stokes *s)
152: {
153: PetscFunctionBeginUser;
154: /* solution vector x */
155: PetscCall(VecCreate(PETSC_COMM_WORLD, &s->x));
156: PetscCall(VecSetSizes(s->x, PETSC_DECIDE, 3 * s->nx * s->ny));
157: PetscCall(VecSetType(s->x, VECMPI));
159: /* exact solution y */
160: PetscCall(VecDuplicate(s->x, &s->y));
161: PetscCall(StokesExactSolution(s));
163: /* rhs vector b */
164: PetscCall(VecDuplicate(s->x, &s->b));
165: PetscCall(StokesRhs(s));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
170: {
171: PetscInt n;
173: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: PetscErrorCode StokesExactSolution(Stokes *s)
182: {
183: PetscInt row, start, end, i, j;
184: PetscScalar val;
185: Vec y0, y1;
187: PetscFunctionBeginUser;
188: /* velocity part */
189: PetscCall(VecGetSubVector(s->y, s->isg[0], &y0));
190: PetscCall(VecGetOwnershipRange(y0, &start, &end));
191: for (row = start; row < end; row++) {
192: PetscCall(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: PetscCall(VecSetValue(y0, row, val, INSERT_VALUES));
199: }
200: PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));
202: /* pressure part */
203: PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
204: PetscCall(VecGetOwnershipRange(y1, &start, &end));
205: for (row = start; row < end; row++) {
206: PetscCall(StokesGetPosition(s, row, &i, &j));
207: val = StokesExactPressure(i * s->hx + s->hx / 2);
208: PetscCall(VecSetValue(y1, row, val, INSERT_VALUES));
209: }
210: PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: PetscErrorCode StokesRhs(Stokes *s)
215: {
216: PetscInt row, start, end, i, j;
217: PetscScalar val;
218: Vec b0, b1;
220: PetscFunctionBeginUser;
221: /* velocity part */
222: PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
223: PetscCall(VecGetOwnershipRange(b0, &start, &end));
224: for (row = start; row < end; row++) {
225: PetscCall(StokesGetPosition(s, row, &i, &j));
226: if (row < s->nx * s->ny) {
227: PetscCall(StokesRhsMomX(s, i, j, &val));
228: } else {
229: PetscCall(StokesRhsMomY(s, i, j, &val));
230: }
231: PetscCall(VecSetValue(b0, row, val, INSERT_VALUES));
232: }
233: PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));
235: /* pressure part */
236: PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
237: PetscCall(VecGetOwnershipRange(b1, &start, &end));
238: for (row = start; row < end; row++) {
239: PetscCall(StokesGetPosition(s, row, &i, &j));
240: PetscCall(StokesRhsMass(s, i, j, &val));
241: if (s->matsymmetric) val = -1.0 * val;
242: PetscCall(VecSetValue(b1, row, val, INSERT_VALUES));
243: }
244: PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
249: {
250: PetscInt row, start, end, sz, i, j;
251: PetscInt cols[5];
252: PetscScalar vals[5];
254: PetscFunctionBeginUser;
255: /* A[0] is 2N-by-2N */
256: PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[0]));
257: PetscCall(MatSetOptionsPrefix(s->subA[0], "a00_"));
258: PetscCall(MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny));
259: PetscCall(MatSetType(s->subA[0], MATMPIAIJ));
260: PetscCall(MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL));
261: PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end));
263: for (row = start; row < end; row++) {
264: PetscCall(StokesGetPosition(s, row, &i, &j));
265: /* first part: rows 0 to (nx*ny-1) */
266: PetscCall(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: PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES));
273: }
274: PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY));
275: PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
280: {
281: PetscInt row, start, end, sz, i, j;
282: PetscInt cols[5];
283: PetscScalar vals[5];
285: PetscFunctionBeginUser;
286: /* A[1] is 2N-by-N */
287: PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1]));
288: PetscCall(MatSetOptionsPrefix(s->subA[1], "a01_"));
289: PetscCall(MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny));
290: PetscCall(MatSetType(s->subA[1], MATMPIAIJ));
291: PetscCall(MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL));
292: PetscCall(MatGetOwnershipRange(s->subA[1], &start, &end));
294: PetscCall(MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
296: for (row = start; row < end; row++) {
297: PetscCall(StokesGetPosition(s, row, &i, &j));
298: /* first part: rows 0 to (nx*ny-1) */
299: if (row < s->nx * s->ny) {
300: PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals));
301: } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */
302: PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals));
303: }
304: PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES));
305: }
306: PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY));
307: PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY));
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
312: {
313: PetscFunctionBeginUser;
314: /* A[2] is minus transpose of A[1] */
315: PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]));
316: if (!s->matsymmetric) PetscCall(MatScale(s->subA[2], -1.0));
317: PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_"));
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
322: {
323: PetscFunctionBeginUser;
324: /* A[3] is N-by-N null matrix */
325: PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[3]));
326: PetscCall(MatSetOptionsPrefix(s->subA[3], "a11_"));
327: PetscCall(MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny));
328: PetscCall(MatSetType(s->subA[3], MATMPIAIJ));
329: PetscCall(MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL));
330: PetscCall(MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY));
331: PetscCall(MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
336: {
337: Vec diag;
339: PetscFunctionBeginUser;
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: PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
347: PetscCall(VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny));
348: PetscCall(VecSetType(diag, VECMPI));
349: PetscCall(MatGetDiagonal(s->subA[0], diag));
350: PetscCall(VecReciprocal(diag));
352: /* compute: - A10 inv(DIAGFORM(A00)) A01 */
353: PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); /* (*warning* overwrites subA[1]) */
354: PetscCall(MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DEFAULT, &s->myS));
355: PetscCall(MatScale(s->myS, -1.0));
357: /* restore A10 */
358: PetscCall(MatGetDiagonal(s->subA[0], diag));
359: PetscCall(MatDiagonalScale(s->subA[1], diag, NULL));
360: PetscCall(VecDestroy(&diag));
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: PetscErrorCode StokesSetupMatrix(Stokes *s)
365: {
366: PetscFunctionBeginUser;
367: PetscCall(StokesSetupMatBlock00(s));
368: PetscCall(StokesSetupMatBlock01(s));
369: PetscCall(StokesSetupMatBlock10(s));
370: PetscCall(StokesSetupMatBlock11(s));
371: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A));
372: PetscCall(StokesSetupApproxSchur(s));
373: PetscFunctionReturn(PETSC_SUCCESS);
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);
384: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
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;
479: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
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;
550: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
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);
620: PetscFunctionBeginUser;
621: if (i == 0) { /* west boundary */
622: *val = awb * StokesExactVelocityX(y);
623: } else {
624: *val = 0.0;
625: }
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
630: {
631: PetscFunctionBeginUser;
632: *val = 0.0;
633: PetscFunctionReturn(PETSC_SUCCESS);
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;
641: PetscFunctionBeginUser;
642: if (i == 0) { /* west boundary */
643: *val = aeb * StokesExactVelocityX(y);
644: } else {
645: *val = 0.0;
646: }
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: PetscErrorCode StokesCalcResidual(Stokes *s)
651: {
652: PetscReal val;
653: Vec b0, b1;
655: PetscFunctionBeginUser;
656: /* residual Ax-b (*warning* overwrites b) */
657: PetscCall(VecScale(s->b, -1.0));
658: PetscCall(MatMultAdd(s->A, s->x, s->b, s->b));
660: /* residual velocity */
661: PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
662: PetscCall(VecNorm(b0, NORM_2, &val));
663: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual u = %g\n", (double)val));
664: PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));
666: /* residual pressure */
667: PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
668: PetscCall(VecNorm(b1, NORM_2, &val));
669: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual p = %g\n", (double)val));
670: PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));
672: /* total residual */
673: PetscCall(VecNorm(s->b, NORM_2, &val));
674: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual [u,p] = %g\n", (double)val));
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: PetscErrorCode StokesCalcError(Stokes *s)
679: {
680: PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny);
681: PetscReal val;
682: Vec y0, y1;
684: PetscFunctionBeginUser;
685: /* error y-x */
686: PetscCall(VecAXPY(s->y, -1.0, s->x));
688: /* error in velocity */
689: PetscCall(VecGetSubVector(s->y, s->isg[0], &y0));
690: PetscCall(VecNorm(y0, NORM_2, &val));
691: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error u = %g\n", (double)(PetscRealPart(val / scale))));
692: PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));
694: /* error in pressure */
695: PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
696: PetscCall(VecNorm(y1, NORM_2, &val));
697: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error p = %g\n", (double)(PetscRealPart(val / scale))));
698: PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));
700: /* total error */
701: PetscCall(VecNorm(s->y, NORM_2, &val));
702: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error [u,p] = %g\n", (double)PetscRealPart((val / scale))));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: int main(int argc, char **argv)
707: {
708: Stokes s;
709: KSP ksp;
711: PetscFunctionBeginUser;
712: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
713: s.nx = 4;
714: s.ny = 6;
715: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL));
716: PetscCall(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: PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL));
721: s.userPC = s.userKSP = PETSC_FALSE;
722: PetscCall(PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC));
723: PetscCall(PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP));
725: PetscCall(StokesSetupMatrix(&s));
726: PetscCall(StokesSetupIndexSets(&s));
727: PetscCall(StokesSetupVectors(&s));
729: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
730: PetscCall(KSPSetOperators(ksp, s.A, s.A));
731: PetscCall(KSPSetFromOptions(ksp));
732: PetscCall(StokesSetupPC(&s, ksp));
733: PetscCall(KSPSolve(ksp, s.b, s.x));
735: /* don't trust, verify! */
736: PetscCall(StokesCalcResidual(&s));
737: PetscCall(StokesCalcError(&s));
738: PetscCall(StokesWriteSolution(&s));
740: PetscCall(KSPDestroy(&ksp));
741: PetscCall(MatDestroy(&s.subA[0]));
742: PetscCall(MatDestroy(&s.subA[1]));
743: PetscCall(MatDestroy(&s.subA[2]));
744: PetscCall(MatDestroy(&s.subA[3]));
745: PetscCall(MatDestroy(&s.A));
746: PetscCall(VecDestroy(&s.x));
747: PetscCall(VecDestroy(&s.b));
748: PetscCall(VecDestroy(&s.y));
749: PetscCall(MatDestroy(&s.myS));
750: PetscCall(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*/