Actual source code: ex70.c
petsc-3.5.4 2015-05-23
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: /*---------------------------------------------------------------------------- */
7: /* author : Christiaan M. Klaij */
8: /*---------------------------------------------------------------------------- */
9: /* */
10: /* Poiseuille flow problem. */
11: /* */
12: /* Viscous, laminar flow in a 2D channel with parabolic velocity */
13: /* profile and linear pressure drop, exact solution of the 2D Stokes */
14: /* equations. */
15: /* */
16: /* Discretized with the cell-centered finite-volume method on a */
17: /* Cartesian grid with co-located variables. Variables ordered as */
18: /* [u1...uN v1...vN p1...pN]^T. Matrix [A00 A01; A10, A11] solved with */
19: /* PCFIELDSPLIT. */
20: /* */
21: /* Disclaimer: does not contain the pressure-weighed interpolation */
22: /* method needed to suppress spurious pressure modes in real-life */
23: /* problems. */
24: /* */
25: /* usage: */
26: /* */
27: /* mpiexec -n 2 ./stokes -nx 32 -ny 48 */
28: /* */
29: /* Runs with PETSc defaults on 32x48 grid, no PC for the Schur */
30: /* complement because A11 is zero. */
31: /* */
32: /* mpiexec -n 2 ./stokes -nx 32 -ny 48 -fieldsplit_1_user_pc */
33: /* */
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). */
37: /* */
38: /* mpiexec -n 2 ./stokes -nx 32 -ny 48 -fieldsplit_1_user_ksp */
39: /* */
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. */
43: /* */
44: /* mpiexec -n 2 ./stokes -nx 32 -ny 48 -fieldsplit_1_user_ksp -fieldsplit_1_ksp_rtol 0.01 -fieldsplit_0_ksp_rtol 0.01 */
45: /* */
46: /* SIMPLE-type approximations are crude, there's no benefit in */
47: /* solving the subsystems in the preconditioner very accurately. */
48: /* */
49: /*---------------------------------------------------------------------------- */
51: #include <petscksp.h>
53: typedef struct {
54: PetscBool userPC, userKSP; /* 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: PetscScalar *array;
128: /* write data (*warning* only works sequential) */
129: MPI_Comm_size(MPI_COMM_WORLD,&size);
130: /*PetscPrintf(PETSC_COMM_WORLD," number of processors = %D\n",size);*/
131: if (size == 1) {
132: PetscViewer viewer;
133: VecGetArray(s->x, &array);
134: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "solution.dat", &viewer);
135: PetscViewerASCIIPrintf(viewer, "# x, y, u, v, p\n");
136: for (j = 0; j < s->ny; j++) {
137: for (i = 0; i < s->nx; i++) {
138: n = j*s->nx+i;
139: PetscViewerASCIIPrintf(viewer, "%.12g %.12g %.12g %.12g %.12g\n", i*s->hx+s->hx/2, j*s->hy+s->hy/2, array[n], array[n+s->nx*s->ny], array[n+2*s->nx*s->ny]);
140: }
141: }
142: VecRestoreArray(s->x, &array);
143: PetscViewerDestroy(&viewer);
144: }
145: return(0);
146: }
148: PetscErrorCode StokesSetupIndexSets(Stokes *s)
149: {
153: /* the two index sets */
154: MatNestGetISs(s->A, s->isg, NULL);
155: /* ISView(isg[0],PETSC_VIEWER_STDOUT_WORLD); */
156: /* ISView(isg[1],PETSC_VIEWER_STDOUT_WORLD); */
157: return(0);
158: }
160: PetscErrorCode StokesSetupVectors(Stokes *s)
161: {
165: /* solution vector x */
166: VecCreate(PETSC_COMM_WORLD, &s->x);
167: VecSetSizes(s->x, PETSC_DECIDE, 3*s->nx*s->ny);
168: VecSetType(s->x, VECMPI);
169: /* VecSetRandom(s->x, NULL); */
170: /* VecView(s->x, (PetscViewer) PETSC_VIEWER_DEFAULT); */
172: /* exact solution y */
173: VecDuplicate(s->x, &s->y);
174: StokesExactSolution(s);
175: /* VecView(s->y, (PetscViewer) PETSC_VIEWER_DEFAULT); */
177: /* rhs vector b */
178: VecDuplicate(s->x, &s->b);
179: StokesRhs(s);
180: /*VecView(s->b, (PetscViewer) PETSC_VIEWER_DEFAULT);*/
181: return(0);
182: }
184: PetscErrorCode StokesGetPosition(Stokes *s, PetscInt row, PetscInt *i, PetscInt *j)
185: {
186: PetscInt n;
189: /* cell number n=j*nx+i has position (i,j) in grid */
190: n = row%(s->nx*s->ny);
191: *i = n%s->nx;
192: *j = (n-(*i))/s->nx;
193: return(0);
194: }
196: PetscErrorCode StokesExactSolution(Stokes *s)
197: {
198: PetscInt row, start, end, i, j;
199: PetscScalar val;
200: Vec y0,y1;
204: /* velocity part */
205: VecGetSubVector(s->y, s->isg[0], &y0);
206: VecGetOwnershipRange(y0, &start, &end);
207: for (row = start; row < end; row++) {
208: StokesGetPosition(s, row,&i,&j);
209: if (row < s->nx*s->ny) {
210: val = StokesExactVelocityX(j*s->hy+s->hy/2);
211: } else {
212: val = 0;
213: }
214: VecSetValue(y0, row, val, INSERT_VALUES);
215: }
216: VecRestoreSubVector(s->y, s->isg[0], &y0);
218: /* pressure part */
219: VecGetSubVector(s->y, s->isg[1], &y1);
220: VecGetOwnershipRange(y1, &start, &end);
221: for (row = start; row < end; row++) {
222: StokesGetPosition(s, row, &i, &j);
223: val = StokesExactPressure(i*s->hx+s->hx/2);
224: VecSetValue(y1, row, val, INSERT_VALUES);
225: }
226: VecRestoreSubVector(s->y, s->isg[1], &y1);
227: return(0);
228: }
230: PetscErrorCode StokesRhs(Stokes *s)
231: {
232: PetscInt row, start, end, i, j;
233: PetscScalar val;
234: Vec b0,b1;
238: /* velocity part */
239: VecGetSubVector(s->b, s->isg[0], &b0);
240: VecGetOwnershipRange(b0, &start, &end);
241: for (row = start; row < end; row++) {
242: StokesGetPosition(s, row, &i, &j);
243: if (row < s->nx*s->ny) {
244: StokesRhsMomX(s, i, j, &val);
245: } else if (row < 2*s->nx*s->ny) {
246: StokesRhsMomY(s, i, j, &val);
247: }
248: VecSetValue(b0, row, val, INSERT_VALUES);
249: }
250: VecRestoreSubVector(s->b, s->isg[0], &b0);
252: /* pressure part */
253: VecGetSubVector(s->b, s->isg[1], &b1);
254: VecGetOwnershipRange(b1, &start, &end);
255: for (row = start; row < end; row++) {
256: StokesGetPosition(s, row, &i, &j);
257: StokesRhsMass(s, i, j, &val);
258: VecSetValue(b1, row, val, INSERT_VALUES);
259: }
260: VecRestoreSubVector(s->b, s->isg[1], &b1);
261: return(0);
262: }
264: PetscErrorCode StokesSetupMatBlock00(Stokes *s)
265: {
266: PetscInt row, start, end, size, i, j;
267: PetscInt cols[5];
268: PetscScalar vals[5];
272: /* A[0] is 2N-by-2N */
273: MatCreate(PETSC_COMM_WORLD,&s->subA[0]);
274: MatSetOptionsPrefix(s->subA[0],"a00_");
275: MatSetSizes(s->subA[0],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,2*s->nx*s->ny);
276: MatSetType(s->subA[0],MATMPIAIJ);
277: MatMPIAIJSetPreallocation(s->subA[0],5,NULL,5,NULL);
278: MatGetOwnershipRange(s->subA[0], &start, &end);
280: for (row = start; row < end; row++) {
281: StokesGetPosition(s, row, &i, &j);
282: /* first part: rows 0 to (nx*ny-1) */
283: StokesStencilLaplacian(s, i, j, &size, cols, vals);
284: /* second part: rows (nx*ny) to (2*nx*ny-1) */
285: if (row >= s->nx*s->ny) {
286: for (i = 0; i < 5; i++) cols[i] = cols[i] + s->nx*s->ny;
287: }
288: for (i = 0; i < 5; i++) vals[i] = -1.0*vals[i]; /* dynamic viscosity coef mu=-1 */
289: MatSetValues(s->subA[0], 1, &row, size, cols, vals, INSERT_VALUES);
290: }
291: MatAssemblyBegin(s->subA[0], MAT_FINAL_ASSEMBLY);
292: MatAssemblyEnd(s->subA[0], MAT_FINAL_ASSEMBLY);
293: return(0);
294: }
296: PetscErrorCode StokesSetupMatBlock01(Stokes *s)
297: {
298: PetscInt row, start, end, size, i, j;
299: PetscInt cols[5];
300: PetscScalar vals[5];
304: /* A[1] is 2N-by-N */
305: MatCreate(PETSC_COMM_WORLD, &s->subA[1]);
306: MatSetOptionsPrefix(s->subA[1],"a01_");
307: MatSetSizes(s->subA[1],PETSC_DECIDE,PETSC_DECIDE,2*s->nx*s->ny,s->nx*s->ny);
308: MatSetType(s->subA[1],MATMPIAIJ);
309: MatMPIAIJSetPreallocation(s->subA[1],5,NULL,5,NULL);
310: MatGetOwnershipRange(s->subA[1],&start,&end);
312: MatSetOption(s->subA[1],MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
314: for (row = start; row < end; row++) {
315: StokesGetPosition(s, row, &i, &j);
316: /* first part: rows 0 to (nx*ny-1) */
317: if (row < s->nx*s->ny) {
318: StokesStencilGradientX(s, i, j, &size, cols, vals);
319: }
320: /* second part: rows (nx*ny) to (2*nx*ny-1) */
321: else {
322: StokesStencilGradientY(s, i, j, &size, cols, vals);
323: }
324: MatSetValues(s->subA[1], 1, &row, size, cols, vals, INSERT_VALUES);
325: }
326: MatAssemblyBegin(s->subA[1], MAT_FINAL_ASSEMBLY);
327: MatAssemblyEnd(s->subA[1], MAT_FINAL_ASSEMBLY);
328: return(0);
329: }
331: PetscErrorCode StokesSetupMatBlock10(Stokes *s)
332: {
336: /* A[2] is minus transpose of A[1] */
337: MatTranspose(s->subA[1], MAT_INITIAL_MATRIX, &s->subA[2]);
338: MatScale(s->subA[2], -1.0);
339: MatSetOptionsPrefix(s->subA[2], "a10_");
340: return(0);
341: }
343: PetscErrorCode StokesSetupMatBlock11(Stokes *s)
344: {
348: /* A[3] is N-by-N null matrix */
349: MatCreate(PETSC_COMM_WORLD, &s->subA[3]);
350: MatSetOptionsPrefix(s->subA[3], "a11_");
351: MatSetSizes(s->subA[3], PETSC_DECIDE, PETSC_DECIDE, s->nx*s->ny, s->nx*s->ny);
352: MatSetType(s->subA[3], MATMPIAIJ);
353: MatMPIAIJSetPreallocation(s->subA[3], 0, NULL, 0, NULL);
354: MatAssemblyBegin(s->subA[3], MAT_FINAL_ASSEMBLY);
355: MatAssemblyEnd(s->subA[3], MAT_FINAL_ASSEMBLY);
356: return(0);
357: }
359: PetscErrorCode StokesSetupApproxSchur(Stokes *s)
360: {
361: Vec diag;
365: /* Schur complement approximation: myS = A11 - A10 diag(A00)^(-1) A01 */
366: /* note: A11 is zero */
367: /* note: in real life this matrix would be build directly, */
368: /* i.e. without MatMatMult */
370: /* inverse of diagonal of A00 */
371: VecCreate(PETSC_COMM_WORLD,&diag);
372: VecSetSizes(diag,PETSC_DECIDE,2*s->nx*s->ny);
373: VecSetType(diag,VECMPI);
374: MatGetDiagonal(s->subA[0],diag);
375: VecReciprocal(diag);
377: /* compute: - A10 diag(A00)^(-1) A01 */
378: MatDiagonalScale(s->subA[1],diag,NULL); /* (*warning* overwrites subA[1]) */
379: MatMatMult(s->subA[2],s->subA[1],MAT_INITIAL_MATRIX,PETSC_DEFAULT,&s->myS);
380: MatScale(s->myS,-1.0);
382: /* restore A10 */
383: MatGetDiagonal(s->subA[0],diag);
384: MatDiagonalScale(s->subA[1],diag,NULL);
385: VecDestroy(&diag);
386: return(0);
387: }
389: PetscErrorCode StokesSetupMatrix(Stokes *s)
390: {
394: StokesSetupMatBlock00(s);
395: StokesSetupMatBlock01(s);
396: StokesSetupMatBlock10(s);
397: StokesSetupMatBlock11(s);
398: MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, s->subA, &s->A);
399: StokesSetupApproxSchur(s);
400: return(0);
401: }
403: PetscErrorCode StokesStencilLaplacian(Stokes *s, PetscInt i, PetscInt j, PetscInt *size, PetscInt *cols, PetscScalar *vals)
404: {
405: PetscInt p =j*s->nx+i, w=p-1, e=p+1, s2=p-s->nx, n=p+s->nx;
406: PetscScalar ae=s->hy/s->hx, aeb=0;
407: PetscScalar aw=s->hy/s->hx, awb=s->hy/(s->hx/2);
408: PetscScalar as=s->hx/s->hy, asb=s->hx/(s->hy/2);
409: PetscScalar an=s->hx/s->hy, anb=s->hx/(s->hy/2);
412: if (i==0 && j==0) { /* south-west corner */
413: *size =3;
414: cols[0]=p; vals[0]=-(ae+awb+asb+an);
415: cols[1]=e; vals[1]=ae;
416: cols[2]=n; vals[2]=an;
417: } else if (i==0 && j==s->ny-1) { /* north-west corner */
418: *size =3;
419: cols[0]=s2; vals[0]=as;
420: cols[1]=p; vals[1]=-(ae+awb+as+anb);
421: cols[2]=e; vals[2]=ae;
422: } else if (i==s->nx-1 && j==0) { /* south-east corner */
423: *size =3;
424: cols[0]=w; vals[0]=aw;
425: cols[1]=p; vals[1]=-(aeb+aw+asb+an);
426: cols[2]=n; vals[2]=an;
427: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
428: *size =3;
429: cols[0]=s2; vals[0]=as;
430: cols[1]=w; vals[1]=aw;
431: cols[2]=p; vals[2]=-(aeb+aw+as+anb);
432: } else if (i==0) { /* west boundary */
433: *size =4;
434: cols[0]=s2; vals[0]=as;
435: cols[1]=p; vals[1]=-(ae+awb+as+an);
436: cols[2]=e; vals[2]=ae;
437: cols[3]=n; vals[3]=an;
438: } else if (i==s->nx-1) { /* east boundary */
439: *size =4;
440: cols[0]=s2; vals[0]=as;
441: cols[1]=w; vals[1]=aw;
442: cols[2]=p; vals[2]=-(aeb+aw+as+an);
443: cols[3]=n; vals[3]=an;
444: } else if (j==0) { /* south boundary */
445: *size =4;
446: cols[0]=w; vals[0]=aw;
447: cols[1]=p; vals[1]=-(ae+aw+asb+an);
448: cols[2]=e; vals[2]=ae;
449: cols[3]=n; vals[3]=an;
450: } else if (j==s->ny-1) { /* north boundary */
451: *size =4;
452: cols[0]=s2; vals[0]=as;
453: cols[1]=w; vals[1]=aw;
454: cols[2]=p; vals[2]=-(ae+aw+as+anb);
455: cols[3]=e; vals[3]=ae;
456: } else { /* interior */
457: *size =5;
458: cols[0]=s2; vals[0]=as;
459: cols[1]=w; vals[1]=aw;
460: cols[2]=p; vals[2]=-(ae+aw+as+an);
461: cols[3]=e; vals[3]=ae;
462: cols[4]=n; vals[4]=an;
463: }
464: return(0);
465: }
467: PetscErrorCode StokesStencilGradientX(Stokes *s, PetscInt i, PetscInt j, PetscInt *size, PetscInt *cols, PetscScalar *vals)
468: {
469: PetscInt p =j*s->nx+i, w=p-1, e=p+1;
470: PetscScalar ae= s->hy/2, aeb=s->hy;
471: PetscScalar aw=-s->hy/2, awb=0;
474: if (i==0 && j==0) { /* south-west corner */
475: *size =2;
476: cols[0]=p; vals[0]=-(ae+awb);
477: cols[1]=e; vals[1]=ae;
478: } else if (i==0 && j==s->ny-1) { /* north-west corner */
479: *size =2;
480: cols[0]=p; vals[0]=-(ae+awb);
481: cols[1]=e; vals[1]=ae;
482: } else if (i==s->nx-1 && j==0) { /* south-east corner */
483: *size =2;
484: cols[0]=w; vals[0]=aw;
485: cols[1]=p; vals[1]=-(aeb+aw);
486: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
487: *size =2;
488: cols[0]=w; vals[0]=aw;
489: cols[1]=p; vals[1]=-(aeb+aw);
490: } else if (i==0) { /* west boundary */
491: *size =2;
492: cols[0]=p; vals[0]=-(ae+awb);
493: cols[1]=e; vals[1]=ae;
494: } else if (i==s->nx-1) { /* east boundary */
495: *size =2;
496: cols[0]=w; vals[0]=aw;
497: cols[1]=p; vals[1]=-(aeb+aw);
498: } else if (j==0) { /* south boundary */
499: *size =3;
500: cols[0]=w; vals[0]=aw;
501: cols[1]=p; vals[1]=-(ae+aw);
502: cols[2]=e; vals[2]=ae;
503: } else if (j==s->ny-1) { /* north boundary */
504: *size =3;
505: cols[0]=w; vals[0]=aw;
506: cols[1]=p; vals[1]=-(ae+aw);
507: cols[2]=e; vals[2]=ae;
508: } else { /* interior */
509: *size =3;
510: cols[0]=w; vals[0]=aw;
511: cols[1]=p; vals[1]=-(ae+aw);
512: cols[2]=e; vals[2]=ae;
513: }
514: return(0);
515: }
517: PetscErrorCode StokesStencilGradientY(Stokes *s, PetscInt i, PetscInt j, PetscInt *size, PetscInt *cols, PetscScalar *vals)
518: {
519: PetscInt p =j*s->nx+i, s2=p-s->nx, n=p+s->nx;
520: PetscScalar as=-s->hx/2, asb=0;
521: PetscScalar an= s->hx/2, anb=0;
524: if (i==0 && j==0) { /* south-west corner */
525: *size =2;
526: cols[0]=p; vals[0]=-(asb+an);
527: cols[1]=n; vals[1]=an;
528: } else if (i==0 && j==s->ny-1) { /* north-west corner */
529: *size =2;
530: cols[0]=s2; vals[0]=as;
531: cols[1]=p; vals[1]=-(as+anb);
532: } else if (i==s->nx-1 && j==0) { /* south-east corner */
533: *size =2;
534: cols[0]=p; vals[0]=-(asb+an);
535: cols[1]=n; vals[1]=an;
536: } else if (i==s->nx-1 && j==s->ny-1) { /* north-east corner */
537: *size =2;
538: cols[0]=s2; vals[0]=as;
539: cols[1]=p; vals[1]=-(as+anb);
540: } else if (i==0) { /* west boundary */
541: *size =3;
542: cols[0]=s2; vals[0]=as;
543: cols[1]=p; vals[1]=-(as+an);
544: cols[2]=n; vals[2]=an;
545: } else if (i==s->nx-1) { /* east boundary */
546: *size =3;
547: cols[0]=s2; vals[0]=as;
548: cols[1]=p; vals[1]=-(as+an);
549: cols[2]=n; vals[2]=an;
550: } else if (j==0) { /* south boundary */
551: *size =2;
552: cols[0]=p; vals[0]=-(asb+an);
553: cols[1]=n; vals[1]=an;
554: } else if (j==s->ny-1) { /* north boundary */
555: *size =2;
556: cols[0]=s2; vals[0]=as;
557: cols[1]=p; vals[1]=-(as+anb);
558: } else { /* interior */
559: *size =3;
560: cols[0]=s2; vals[0]=as;
561: cols[1]=p; vals[1]=-(as+an);
562: cols[2]=n; vals[2]=an;
563: }
564: return(0);
565: }
567: PetscErrorCode StokesRhsMomX(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
568: {
569: PetscScalar y = j*s->hy+s->hy/2;
570: PetscScalar awb = s->hy/(s->hx/2);
573: if (i == 0) { /* west boundary */
574: *val = awb*StokesExactVelocityX(y);
575: } else {
576: *val = 0.0;
577: }
578: return(0);
579: }
581: PetscErrorCode StokesRhsMomY(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
582: {
584: *val = 0.0;
585: return(0);
586: }
588: PetscErrorCode StokesRhsMass(Stokes *s, PetscInt i, PetscInt j, PetscScalar *val)
589: {
590: PetscScalar y = j*s->hy+s->hy/2;
591: PetscScalar aeb = s->hy;
594: if (i == 0) { /* west boundary */
595: *val = aeb*StokesExactVelocityX(y);
596: } else {
597: *val = 0.0;
598: }
599: return(0);
600: }
602: PetscErrorCode StokesCalcResidual(Stokes *s)
603: {
604: PetscReal val;
605: Vec b0, b1;
609: /* residual Ax-b (*warning* overwrites b) */
610: VecScale(s->b, -1.0);
611: MatMultAdd(s->A, s->x, s->b, s->b);
612: /* VecView(s->b, (PetscViewer)PETSC_VIEWER_DEFAULT); */
614: /* residual velocity */
615: VecGetSubVector(s->b, s->isg[0], &b0);
616: VecNorm(b0, NORM_2, &val);
617: PetscPrintf(PETSC_COMM_WORLD," residual u = %g\n",(double)val);
618: VecRestoreSubVector(s->b, s->isg[0], &b0);
620: /* residual pressure */
621: VecGetSubVector(s->b, s->isg[1], &b1);
622: VecNorm(b1, NORM_2, &val);
623: PetscPrintf(PETSC_COMM_WORLD," residual p = %g\n",(double)val);
624: VecRestoreSubVector(s->b, s->isg[1], &b1);
626: /* total residual */
627: VecNorm(s->b, NORM_2, &val);
628: PetscPrintf(PETSC_COMM_WORLD," residual [u,p] = %g\n", (double)val);
629: return(0);
630: }
632: PetscErrorCode StokesCalcError(Stokes *s)
633: {
634: PetscScalar scale = PetscSqrtScalar(s->nx*s->ny);
635: PetscReal val;
636: Vec y0, y1;
640: /* error y-x */
641: VecAXPY(s->y, -1.0, s->x);
642: /* VecView(s->y, (PetscViewer)PETSC_VIEWER_DEFAULT); */
644: /* error in velocity */
645: VecGetSubVector(s->y, s->isg[0], &y0);
646: VecNorm(y0, NORM_2, &val);
647: PetscPrintf(PETSC_COMM_WORLD," discretization error u = %g\n",(double)(val/scale));
648: VecRestoreSubVector(s->y, s->isg[0], &y0);
650: /* error in pressure */
651: VecGetSubVector(s->y, s->isg[1], &y1);
652: VecNorm(y1, NORM_2, &val);
653: PetscPrintf(PETSC_COMM_WORLD," discretization error p = %g\n",(double)(val/scale));
654: VecRestoreSubVector(s->y, s->isg[1], &y1);
656: /* total error */
657: VecNorm(s->y, NORM_2, &val);
658: PetscPrintf(PETSC_COMM_WORLD," discretization error [u,p] = %g\n", (double)(val/scale));
659: return(0);
660: }
662: int main(int argc, char **argv)
663: {
664: Stokes s;
665: KSP ksp;
668: PetscInitialize(&argc, &argv, NULL, help);
669: s.nx = 4;
670: s.ny = 6;
671: PetscOptionsGetInt(NULL, "-nx", &s.nx, NULL);
672: PetscOptionsGetInt(NULL, "-ny", &s.ny, NULL);
673: s.hx = 2.0/s.nx;
674: s.hy = 1.0/s.ny;
675: s.userPC = s.userKSP = PETSC_FALSE;
676: PetscOptionsHasName(NULL, "-user_pc", &s.userPC);
677: PetscOptionsHasName(NULL, "-user_ksp", &s.userKSP);
679: StokesSetupMatrix(&s);
680: StokesSetupIndexSets(&s);
681: StokesSetupVectors(&s);
683: KSPCreate(PETSC_COMM_WORLD, &ksp);
684: KSPSetOperators(ksp, s.A, s.A);
685: KSPSetFromOptions(ksp);
686: StokesSetupPC(&s, ksp);
687: KSPSolve(ksp, s.b, s.x);
689: /* don't trust, verify! */
690: StokesCalcResidual(&s);
691: StokesCalcError(&s);
692: StokesWriteSolution(&s);
694: KSPDestroy(&ksp);
695: MatDestroy(&s.subA[0]);
696: MatDestroy(&s.subA[1]);
697: MatDestroy(&s.subA[2]);
698: MatDestroy(&s.subA[3]);
699: MatDestroy(&s.A);
700: VecDestroy(&s.x);
701: VecDestroy(&s.b);
702: VecDestroy(&s.y);
703: MatDestroy(&s.myS);
704: PetscFinalize();
705: return 0;
706: }