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(VecAssemblyBegin(y0));
201: PetscCall(VecAssemblyEnd(y0));
202: PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));
204: /* pressure part */
205: PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
206: PetscCall(VecGetOwnershipRange(y1, &start, &end));
207: for (row = start; row < end; row++) {
208: PetscCall(StokesGetPosition(s, row, &i, &j));
209: val = StokesExactPressure(i * s->hx + s->hx / 2);
210: PetscCall(VecSetValue(y1, row, val, INSERT_VALUES));
211: }
212: PetscCall(VecAssemblyBegin(y1));
213: PetscCall(VecAssemblyEnd(y1));
214: PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: PetscErrorCode StokesRhs(Stokes *s)
219: {
220: PetscInt row, start, end, i, j;
221: PetscScalar val;
222: Vec b0, b1;
224: PetscFunctionBeginUser;
225: /* velocity part */
226: PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
227: PetscCall(VecGetOwnershipRange(b0, &start, &end));
228: for (row = start; row < end; row++) {
229: PetscCall(StokesGetPosition(s, row, &i, &j));
230: if (row < s->nx * s->ny) {
231: PetscCall(StokesRhsMomX(s, i, j, &val));
232: } else {
233: PetscCall(StokesRhsMomY(s, i, j, &val));
234: }
235: PetscCall(VecSetValue(b0, row, val, INSERT_VALUES));
236: }
237: PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));
239: /* pressure part */
240: PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
241: PetscCall(VecGetOwnershipRange(b1, &start, &end));
242: for (row = start; row < end; row++) {
243: PetscCall(StokesGetPosition(s, row, &i, &j));
244: PetscCall(StokesRhsMass(s, i, j, &val));
245: if (s->matsymmetric) val = -1.0 * val;
246: PetscCall(VecSetValue(b1, row, val, INSERT_VALUES));
247: }
248: PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
253: {
254: PetscInt row, start, end, sz, i, j;
255: PetscInt cols[5];
256: PetscScalar vals[5];
258: PetscFunctionBeginUser;
259: /* A[0] is 2N-by-2N */
260: PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[0]));
261: PetscCall(MatSetOptionsPrefix(s->subA[0], "a00_"));
262: PetscCall(MatSetSizes(s->subA[0], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, 2 * s->nx * s->ny));
263: PetscCall(MatSetType(s->subA[0], MATMPIAIJ));
264: PetscCall(MatMPIAIJSetPreallocation(s->subA[0], 5, NULL, 5, NULL));
265: PetscCall(MatGetOwnershipRange(s->subA[0], &start, &end));
267: for (row = start; row < end; row++) {
268: PetscCall(StokesGetPosition(s, row, &i, &j));
269: /* first part: rows 0 to (nx*ny-1) */
270: PetscCall(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: PetscCall(MatSetValues(s->subA[0], 1, &row, sz, cols, vals, INSERT_VALUES));
277: }
278: PetscCall(MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY));
279: PetscCall(MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
284: {
285: PetscInt row, start, end, sz, i, j;
286: PetscInt cols[5];
287: PetscScalar vals[5];
289: PetscFunctionBeginUser;
290: /* A[1] is 2N-by-N */
291: PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[1]));
292: PetscCall(MatSetOptionsPrefix(s->subA[1], "a01_"));
293: PetscCall(MatSetSizes(s->subA[1], PETSC_DECIDE, PETSC_DECIDE, 2 * s->nx * s->ny, s->nx * s->ny));
294: PetscCall(MatSetType(s->subA[1], MATMPIAIJ));
295: PetscCall(MatMPIAIJSetPreallocation(s->subA[1], 5, NULL, 5, NULL));
296: PetscCall(MatGetOwnershipRange(s->subA[1], &start, &end));
298: PetscCall(MatSetOption(s->subA[1], MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE));
300: for (row = start; row < end; row++) {
301: PetscCall(StokesGetPosition(s, row, &i, &j));
302: /* first part: rows 0 to (nx*ny-1) */
303: if (row < s->nx * s->ny) {
304: PetscCall(StokesStencilGradientX(s, i, j, &sz, cols, vals));
305: } else { /* second part: rows (nx*ny) to (2*nx*ny-1) */
306: PetscCall(StokesStencilGradientY(s, i, j, &sz, cols, vals));
307: }
308: PetscCall(MatSetValues(s->subA[1], 1, &row, sz, cols, vals, INSERT_VALUES));
309: }
310: PetscCall(MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY));
311: PetscCall(MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY));
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
316: {
317: PetscFunctionBeginUser;
318: /* A[2] is minus transpose of A[1] */
319: PetscCall(MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]));
320: if (!s->matsymmetric) PetscCall(MatScale(s->subA[2], -1.0));
321: PetscCall(MatSetOptionsPrefix(s->subA[2], "a10_"));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
326: {
327: PetscFunctionBeginUser;
328: /* A[3] is N-by-N null matrix */
329: PetscCall(MatCreate(PETSC_COMM_WORLD, &s->subA[3]));
330: PetscCall(MatSetOptionsPrefix(s->subA[3], "a11_"));
331: PetscCall(MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx * s->ny, s->nx * s->ny));
332: PetscCall(MatSetType(s->subA[3], MATMPIAIJ));
333: PetscCall(MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL));
334: PetscCall(MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY));
335: PetscCall(MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
340: {
341: Vec diag;
343: PetscFunctionBeginUser;
344: /* Schur complement approximation: myS = A11 - A10 inv(DIAGFORM(A00)) A01 */
345: /* note: A11 is zero */
346: /* note: in real life this matrix would be build directly, */
347: /* i.e. without MatMatMult */
349: /* inverse of diagonal of A00 */
350: PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
351: PetscCall(VecSetSizes(diag, PETSC_DECIDE, 2 * s->nx * s->ny));
352: PetscCall(VecSetType(diag, VECMPI));
353: PetscCall(MatGetDiagonal(s->subA[0], diag));
354: PetscCall(VecReciprocal(diag));
356: /* compute: - A10 inv(DIAGFORM(A00)) A01 */
357: PetscCall(MatDiagonalScale(s->subA[1], diag, NULL)); /* (*warning* overwrites subA[1]) */
358: PetscCall(MatMatMult(s->subA[2], s->subA[1], MAT_INITIAL_MATRIX, PETSC_DETERMINE, &s->myS));
359: PetscCall(MatScale(s->myS, -1.0));
361: /* restore A10 */
362: PetscCall(MatGetDiagonal(s->subA[0], diag));
363: PetscCall(MatDiagonalScale(s->subA[1], diag, NULL));
364: PetscCall(VecDestroy(&diag));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: PetscErrorCode StokesSetupMatrix(Stokes *s)
369: {
370: PetscFunctionBeginUser;
371: PetscCall(StokesSetupMatBlock00(s));
372: PetscCall(StokesSetupMatBlock01(s));
373: PetscCall(StokesSetupMatBlock10(s));
374: PetscCall(StokesSetupMatBlock11(s));
375: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A));
376: PetscCall(StokesSetupApproxSchur(s));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
381: {
382: PetscInt p = j * s->nx + i, w = p - 1, e = p + 1, s2 = p - s->nx, n = p + s->nx;
383: PetscScalar ae = s->hy / s->hx, aeb = 0;
384: PetscScalar aw = s->hy / s->hx, awb = s->hy / (s->hx / 2);
385: PetscScalar as = s->hx / s->hy, asb = s->hx / (s->hy / 2);
386: PetscScalar an = s->hx / s->hy, anb = s->hx / (s->hy / 2);
388: PetscFunctionBeginUser;
389: if (i == 0 && j == 0) { /* south-west corner */
390: *sz = 3;
391: cols[0] = p;
392: vals[0] = -(ae + awb + asb + an);
393: cols[1] = e;
394: vals[1] = ae;
395: cols[2] = n;
396: vals[2] = an;
397: } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
398: *sz = 3;
399: cols[0] = s2;
400: vals[0] = as;
401: cols[1] = p;
402: vals[1] = -(ae + awb + as + anb);
403: cols[2] = e;
404: vals[2] = ae;
405: } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
406: *sz = 3;
407: cols[0] = w;
408: vals[0] = aw;
409: cols[1] = p;
410: vals[1] = -(aeb + aw + asb + an);
411: cols[2] = n;
412: vals[2] = an;
413: } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
414: *sz = 3;
415: cols[0] = s2;
416: vals[0] = as;
417: cols[1] = w;
418: vals[1] = aw;
419: cols[2] = p;
420: vals[2] = -(aeb + aw + as + anb);
421: } else if (i == 0) { /* west boundary */
422: *sz = 4;
423: cols[0] = s2;
424: vals[0] = as;
425: cols[1] = p;
426: vals[1] = -(ae + awb + as + an);
427: cols[2] = e;
428: vals[2] = ae;
429: cols[3] = n;
430: vals[3] = an;
431: } else if (i == s->nx - 1) { /* east boundary */
432: *sz = 4;
433: cols[0] = s2;
434: vals[0] = as;
435: cols[1] = w;
436: vals[1] = aw;
437: cols[2] = p;
438: vals[2] = -(aeb + aw + as + an);
439: cols[3] = n;
440: vals[3] = an;
441: } else if (j == 0) { /* south boundary */
442: *sz = 4;
443: cols[0] = w;
444: vals[0] = aw;
445: cols[1] = p;
446: vals[1] = -(ae + aw + asb + an);
447: cols[2] = e;
448: vals[2] = ae;
449: cols[3] = n;
450: vals[3] = an;
451: } else if (j == s->ny - 1) { /* north boundary */
452: *sz = 4;
453: cols[0] = s2;
454: vals[0] = as;
455: cols[1] = w;
456: vals[1] = aw;
457: cols[2] = p;
458: vals[2] = -(ae + aw + as + anb);
459: cols[3] = e;
460: vals[3] = ae;
461: } else { /* interior */
462: *sz = 5;
463: cols[0] = s2;
464: vals[0] = as;
465: cols[1] = w;
466: vals[1] = aw;
467: cols[2] = p;
468: vals[2] = -(ae + aw + as + an);
469: cols[3] = e;
470: vals[3] = ae;
471: cols[4] = n;
472: vals[4] = an;
473: }
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
478: {
479: PetscInt p = j * s->nx + i, w = p - 1, e = p + 1;
480: PetscScalar ae = s->hy / 2, aeb = s->hy;
481: PetscScalar aw = -s->hy / 2, awb = 0;
483: PetscFunctionBeginUser;
484: if (i == 0 && j == 0) { /* south-west corner */
485: *sz = 2;
486: cols[0] = p;
487: vals[0] = -(ae + awb);
488: cols[1] = e;
489: vals[1] = ae;
490: } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
491: *sz = 2;
492: cols[0] = p;
493: vals[0] = -(ae + awb);
494: cols[1] = e;
495: vals[1] = ae;
496: } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
497: *sz = 2;
498: cols[0] = w;
499: vals[0] = aw;
500: cols[1] = p;
501: vals[1] = -(aeb + aw);
502: } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
503: *sz = 2;
504: cols[0] = w;
505: vals[0] = aw;
506: cols[1] = p;
507: vals[1] = -(aeb + aw);
508: } else if (i == 0) { /* west boundary */
509: *sz = 2;
510: cols[0] = p;
511: vals[0] = -(ae + awb);
512: cols[1] = e;
513: vals[1] = ae;
514: } else if (i == s->nx - 1) { /* east boundary */
515: *sz = 2;
516: cols[0] = w;
517: vals[0] = aw;
518: cols[1] = p;
519: vals[1] = -(aeb + aw);
520: } else if (j == 0) { /* south boundary */
521: *sz = 3;
522: cols[0] = w;
523: vals[0] = aw;
524: cols[1] = p;
525: vals[1] = -(ae + aw);
526: cols[2] = e;
527: vals[2] = ae;
528: } else if (j == s->ny - 1) { /* north boundary */
529: *sz = 3;
530: cols[0] = w;
531: vals[0] = aw;
532: cols[1] = p;
533: vals[1] = -(ae + aw);
534: cols[2] = e;
535: vals[2] = ae;
536: } else { /* interior */
537: *sz = 3;
538: cols[0] = w;
539: vals[0] = aw;
540: cols[1] = p;
541: vals[1] = -(ae + aw);
542: cols[2] = e;
543: vals[2] = ae;
544: }
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *sz, PetscInt *cols, PetscScalar *vals)
549: {
550: PetscInt p = j * s->nx + i, s2 = p - s->nx, n = p + s->nx;
551: PetscScalar as = -s->hx / 2, asb = 0;
552: PetscScalar an = s->hx / 2, anb = 0;
554: PetscFunctionBeginUser;
555: if (i == 0 && j == 0) { /* south-west corner */
556: *sz = 2;
557: cols[0] = p;
558: vals[0] = -(asb + an);
559: cols[1] = n;
560: vals[1] = an;
561: } else if (i == 0 && j == s->ny - 1) { /* north-west corner */
562: *sz = 2;
563: cols[0] = s2;
564: vals[0] = as;
565: cols[1] = p;
566: vals[1] = -(as + anb);
567: } else if (i == s->nx - 1 && j == 0) { /* south-east corner */
568: *sz = 2;
569: cols[0] = p;
570: vals[0] = -(asb + an);
571: cols[1] = n;
572: vals[1] = an;
573: } else if (i == s->nx - 1 && j == s->ny - 1) { /* north-east corner */
574: *sz = 2;
575: cols[0] = s2;
576: vals[0] = as;
577: cols[1] = p;
578: vals[1] = -(as + anb);
579: } else if (i == 0) { /* west boundary */
580: *sz = 3;
581: cols[0] = s2;
582: vals[0] = as;
583: cols[1] = p;
584: vals[1] = -(as + an);
585: cols[2] = n;
586: vals[2] = an;
587: } else if (i == s->nx - 1) { /* east boundary */
588: *sz = 3;
589: cols[0] = s2;
590: vals[0] = as;
591: cols[1] = p;
592: vals[1] = -(as + an);
593: cols[2] = n;
594: vals[2] = an;
595: } else if (j == 0) { /* south boundary */
596: *sz = 2;
597: cols[0] = p;
598: vals[0] = -(asb + an);
599: cols[1] = n;
600: vals[1] = an;
601: } else if (j == s->ny - 1) { /* north boundary */
602: *sz = 2;
603: cols[0] = s2;
604: vals[0] = as;
605: cols[1] = p;
606: vals[1] = -(as + anb);
607: } else { /* interior */
608: *sz = 3;
609: cols[0] = s2;
610: vals[0] = as;
611: cols[1] = p;
612: vals[1] = -(as + an);
613: cols[2] = n;
614: vals[2] = an;
615: }
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
620: {
621: PetscScalar y = j * s->hy + s->hy / 2;
622: PetscScalar awb = s->hy / (s->hx / 2);
624: PetscFunctionBeginUser;
625: if (i == 0) { /* west boundary */
626: *val = awb * StokesExactVelocityX(y);
627: } else {
628: *val = 0.0;
629: }
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
634: {
635: PetscFunctionBeginUser;
636: *val = 0.0;
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
641: {
642: PetscScalar y = j * s->hy + s->hy / 2;
643: PetscScalar aeb = s->hy;
645: PetscFunctionBeginUser;
646: if (i == 0) { /* west boundary */
647: *val = aeb * StokesExactVelocityX(y);
648: } else {
649: *val = 0.0;
650: }
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: PetscErrorCode StokesCalcResidual(Stokes *s)
655: {
656: PetscReal val;
657: Vec b0, b1;
659: PetscFunctionBeginUser;
660: /* residual Ax-b (*warning* overwrites b) */
661: PetscCall(VecScale(s->b, -1.0));
662: PetscCall(MatMultAdd(s->A, s->x, s->b, s->b));
664: /* residual velocity */
665: PetscCall(VecGetSubVector(s->b, s->isg[0], &b0));
666: PetscCall(VecNorm(b0, NORM_2, &val));
667: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual u = %g\n", (double)val));
668: PetscCall(VecRestoreSubVector(s->b, s->isg[0], &b0));
670: /* residual pressure */
671: PetscCall(VecGetSubVector(s->b, s->isg[1], &b1));
672: PetscCall(VecNorm(b1, NORM_2, &val));
673: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual p = %g\n", (double)val));
674: PetscCall(VecRestoreSubVector(s->b, s->isg[1], &b1));
676: /* total residual */
677: PetscCall(VecNorm(s->b, NORM_2, &val));
678: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual [u,p] = %g\n", (double)val));
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: PetscErrorCode StokesCalcError(Stokes *s)
683: {
684: PetscScalar scale = PetscSqrtReal((double)s->nx * s->ny);
685: PetscReal val;
686: Vec y0, y1;
688: PetscFunctionBeginUser;
689: /* error y-x */
690: PetscCall(VecAXPY(s->y, -1.0, s->x));
692: /* error in velocity */
693: PetscCall(VecGetSubVector(s->y, s->isg[0], &y0));
694: PetscCall(VecNorm(y0, NORM_2, &val));
695: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error u = %g\n", (double)(PetscRealPart(val / scale))));
696: PetscCall(VecRestoreSubVector(s->y, s->isg[0], &y0));
698: /* error in pressure */
699: PetscCall(VecGetSubVector(s->y, s->isg[1], &y1));
700: PetscCall(VecNorm(y1, NORM_2, &val));
701: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error p = %g\n", (double)(PetscRealPart(val / scale))));
702: PetscCall(VecRestoreSubVector(s->y, s->isg[1], &y1));
704: /* total error */
705: PetscCall(VecNorm(s->y, NORM_2, &val));
706: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " discretization error [u,p] = %g\n", (double)PetscRealPart(val / scale)));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: int main(int argc, char **argv)
711: {
712: Stokes s;
713: KSP ksp;
715: PetscFunctionBeginUser;
716: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
717: s.nx = 4;
718: s.ny = 6;
719: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &s.nx, NULL));
720: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ny", &s.ny, NULL));
721: s.hx = 2.0 / s.nx;
722: s.hy = 1.0 / s.ny;
723: s.matsymmetric = PETSC_FALSE;
724: PetscCall(PetscOptionsGetBool(NULL, NULL, "-mat_set_symmetric", &s.matsymmetric, NULL));
725: s.userPC = s.userKSP = PETSC_FALSE;
726: PetscCall(PetscOptionsHasName(NULL, NULL, "-user_pc", &s.userPC));
727: PetscCall(PetscOptionsHasName(NULL, NULL, "-user_ksp", &s.userKSP));
729: PetscCall(StokesSetupMatrix(&s));
730: PetscCall(StokesSetupIndexSets(&s));
731: PetscCall(StokesSetupVectors(&s));
733: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
734: PetscCall(KSPSetOperators(ksp, s.A, s.A));
735: PetscCall(KSPSetFromOptions(ksp));
736: PetscCall(StokesSetupPC(&s, ksp));
737: PetscCall(KSPSolve(ksp, s.b, s.x));
739: /* don't trust, verify! */
740: PetscCall(StokesCalcResidual(&s));
741: PetscCall(StokesCalcError(&s));
742: PetscCall(StokesWriteSolution(&s));
744: PetscCall(KSPDestroy(&ksp));
745: PetscCall(MatDestroy(&s.subA[0]));
746: PetscCall(MatDestroy(&s.subA[1]));
747: PetscCall(MatDestroy(&s.subA[2]));
748: PetscCall(MatDestroy(&s.subA[3]));
749: PetscCall(MatDestroy(&s.A));
750: PetscCall(VecDestroy(&s.x));
751: PetscCall(VecDestroy(&s.b));
752: PetscCall(VecDestroy(&s.y));
753: PetscCall(MatDestroy(&s.myS));
754: PetscCall(PetscFinalize());
755: return 0;
756: }
758: /*TEST
760: test:
761: nsize: 2
762: 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
764: test:
765: suffix: 2
766: nsize: 2
767: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
769: test:
770: suffix: 3
771: nsize: 2
772: args: -nx 16 -ny 24 -ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -user_pc
774: test:
775: suffix: 4
776: nsize: 2
777: 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
779: test:
780: suffix: 4_pcksp
781: nsize: 2
782: 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
784: test:
785: suffix: 5
786: nsize: 2
787: 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
789: test:
790: suffix: 6
791: nsize: 2
792: 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
794: test:
795: suffix: 7
796: nsize: 2
797: 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
799: TEST*/