Actual source code: ex18.c
1: static char help[] = "Test: Solve a toy 2D problem on a staggered grid using 2-level multigrid. \n"
2: "The solve is only applied to the u-u block.\n\n";
3: /*
5: Solves only the velocity block of a isoviscous incompressible Stokes problem on a rectangular 2D
6: domain.
8: u_xx + u_yy - p_x = f^x
9: v_xx + v_yy - p_y = f^y
10: u_x + v_y = g
12: g is 0 in the physical case.
14: Boundary conditions give prescribed flow perpendicular to the boundaries,
15: and zero derivative perpendicular to them (free slip).
17: Supply the -analyze flag to activate a custom KSP monitor. Note that
18: this does an additional direct solve of the velocity block of the system to have an "exact"
19: solution to the discrete system available (see KSPSetOptionsPrefix
20: call below for the prefix).
22: This is for testing purposes, and uses some routines to make
23: sure that transfer operators are consistent with extracting submatrices.
25: -extractTransferOperators (true by default) defines transfer operators for
26: the velocity-velocity system by extracting submatrices from the operators for
27: the full system.
29: */
30: #include <petscdm.h>
31: #include <petscksp.h>
32: #include <petscdmstag.h>
34: /* Shorter, more convenient names for DMStagStencilLocation entries */
35: #define DOWN_LEFT DMSTAG_DOWN_LEFT
36: #define DOWN DMSTAG_DOWN
37: #define DOWN_RIGHT DMSTAG_DOWN_RIGHT
38: #define LEFT DMSTAG_LEFT
39: #define ELEMENT DMSTAG_ELEMENT
40: #define RIGHT DMSTAG_RIGHT
41: #define UP_LEFT DMSTAG_UP_LEFT
42: #define UP DMSTAG_UP
43: #define UP_RIGHT DMSTAG_UP_RIGHT
45: static PetscErrorCode CreateSystem(DM, Mat *, Vec *);
46: static PetscErrorCode CreateNumericalReferenceSolution(Mat, Vec, Vec *);
47: static PetscErrorCode DMStagAnalysisKSPMonitor(KSP, PetscInt, PetscReal, void *);
48: typedef struct {
49: DM dm;
50: Vec solRef, solRefNum, solPrev;
51: } DMStagAnalysisKSPMonitorContext;
53: /* Manufactured solution. Chosen to be higher order than can be solved exactly,
54: and to have a zero derivative for flow parallel to the boundaries. That is,
55: d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
56: and left boundaries. */
57: static PetscScalar uxRef(PetscScalar x, PetscScalar y)
58: {
59: return 0.0 * x + y * y - 2.0 * y * y * y + y * y * y * y;
60: } /* no x-dependence */
61: static PetscScalar uyRef(PetscScalar x, PetscScalar y)
62: {
63: return x * x - 2.0 * x * x * x + x * x * x * x + 0.0 * y;
64: } /* no y-dependence */
65: static PetscScalar fx(PetscScalar x, PetscScalar y)
66: {
67: return 0.0 * x + 2.0 - 12.0 * y + 12.0 * y * y + 1.0;
68: } /* no x-dependence */
69: static PetscScalar fy(PetscScalar x, PetscScalar y)
70: {
71: return 2.0 - 12.0 * x + 12.0 * x * x + 3.0 * y;
72: }
73: static PetscScalar g(PetscScalar x, PetscScalar y)
74: {
75: return 0.0 * x * y;
76: } /* identically zero */
78: int main(int argc, char **argv)
79: {
80: DM dmSol, dmSolc, dmuu, dmuuc;
81: KSP ksp, kspc;
82: PC pc;
83: Mat IIu, Ru, Auu, Auuc;
84: Vec su, xu, fu;
85: IS isuf, ispf, isuc, ispc;
86: PetscInt cnt = 0;
87: DMStagStencil stencil_set[1 + 1 + 1];
88: PetscBool extractTransferOperators, extractSystem;
90: DMStagAnalysisKSPMonitorContext mctx;
91: PetscBool analyze;
93: /* Initialize PETSc and process command line arguments */
95: PetscInitialize(&argc, &argv, (char *)0, help);
96: analyze = PETSC_FALSE;
97: PetscOptionsGetBool(NULL, NULL, "-analyze", &analyze, NULL);
98: extractTransferOperators = PETSC_TRUE;
99: PetscOptionsGetBool(NULL, NULL, "-extractTransferOperators", &extractTransferOperators, NULL);
100: extractSystem = PETSC_TRUE;
101: PetscOptionsGetBool(NULL, NULL, "-extractSystem", &extractSystem, NULL);
103: /* Create 2D DMStag for the solution, and set up. */
104: {
105: const PetscInt dof0 = 0, dof1 = 1, dof2 = 1; /* 1 dof on each edge and element center */
106: const PetscInt stencilWidth = 1;
107: DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 8, 8, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, &dmSol);
108: DMSetFromOptions(dmSol);
109: DMSetUp(dmSol);
110: }
112: /* Create coarse DMStag (which we may or may not directly use) */
113: DMCoarsen(dmSol, MPI_COMM_NULL, &dmSolc);
114: DMSetUp(dmSolc);
116: /* Create compatible DMStags with only velocity dof (which we may or may not
117: directly use) */
118: DMStagCreateCompatibleDMStag(dmSol, 0, 1, 0, 0, &dmuu); /* vel-only */
119: DMSetUp(dmuu);
120: DMCoarsen(dmuu, MPI_COMM_NULL, &dmuuc);
121: DMSetUp(dmuuc);
123: /* Define uniform coordinates as a product of 1D arrays */
124: DMStagSetUniformCoordinatesProduct(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
125: DMStagSetUniformCoordinatesProduct(dmSolc, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
126: DMStagSetUniformCoordinatesProduct(dmuu, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
127: DMStagSetUniformCoordinatesProduct(dmuuc, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
129: /* Create ISs for the velocity and pressure blocks */
130: /* i,j will be ignored */
131: stencil_set[cnt].loc = DMSTAG_DOWN;
132: stencil_set[cnt].c = 0;
133: cnt++; /* u */
134: stencil_set[cnt].loc = DMSTAG_LEFT;
135: stencil_set[cnt].c = 0;
136: cnt++; /* v */
137: stencil_set[cnt].loc = DMSTAG_ELEMENT;
138: stencil_set[cnt].c = 0;
139: cnt++; /* p */
141: DMStagCreateISFromStencils(dmSol, 2, stencil_set, &isuf);
142: DMStagCreateISFromStencils(dmSol, 1, &stencil_set[2], &ispf);
143: DMStagCreateISFromStencils(dmSolc, 2, stencil_set, &isuc);
144: DMStagCreateISFromStencils(dmSolc, 1, &stencil_set[2], &ispc);
146: /* Assemble velocity-velocity system */
147: if (extractSystem) {
148: Mat A, Ac;
149: Vec tmp, rhs;
151: CreateSystem(dmSol, &A, &rhs);
152: CreateSystem(dmSolc, &Ac, NULL);
153: MatCreateSubMatrix(A, isuf, isuf, MAT_INITIAL_MATRIX, &Auu);
154: MatCreateSubMatrix(Ac, isuc, isuc, MAT_INITIAL_MATRIX, &Auuc);
155: MatCreateVecs(Auu, &xu, &fu);
156: VecGetSubVector(rhs, isuf, &tmp);
157: VecCopy(tmp, fu);
158: VecRestoreSubVector(rhs, isuf, &tmp);
159: MatDestroy(&Ac);
160: VecDestroy(&rhs);
161: MatDestroy(&A);
162: } else {
163: CreateSystem(dmuu, &Auu, &fu);
164: CreateSystem(dmuuc, &Auuc, NULL);
165: MatCreateVecs(Auu, &xu, NULL);
166: }
167: PetscObjectSetName((PetscObject)Auu, "Auu");
168: PetscObjectSetName((PetscObject)Auuc, "Auuc");
170: /* Create Transfer Operators and scaling for the velocity-velocity block */
171: if (extractTransferOperators) {
172: Mat II, R;
173: Vec s, tmp;
175: DMCreateInterpolation(dmSolc, dmSol, &II, NULL);
176: DMCreateRestriction(dmSolc, dmSol, &R);
177: MatCreateSubMatrix(II, isuf, isuc, MAT_INITIAL_MATRIX, &IIu);
178: MatCreateSubMatrix(R, isuc, isuf, MAT_INITIAL_MATRIX, &Ru);
179: DMCreateInterpolationScale(dmSolc, dmSol, II, &s);
180: MatCreateVecs(IIu, &su, NULL);
181: VecGetSubVector(s, isuc, &tmp);
182: VecCopy(tmp, su);
183: VecRestoreSubVector(s, isuc, &tmp);
184: MatDestroy(&R);
185: MatDestroy(&II);
186: VecDestroy(&s);
187: } else {
188: DMCreateInterpolation(dmuuc, dmuu, &IIu, NULL);
189: DMCreateInterpolationScale(dmuuc, dmuu, IIu, &su);
190: DMCreateRestriction(dmuuc, dmuu, &Ru);
191: }
193: /* Create and configure solver */
194: KSPCreate(PETSC_COMM_WORLD, &ksp);
195: KSPSetOperators(ksp, Auu, Auu);
196: KSPGetPC(ksp, &pc);
197: PCSetType(pc, PCMG);
198: PCMGSetLevels(pc, 2, NULL);
199: PCMGSetInterpolation(pc, 1, IIu);
200: /* PCMGSetRestriction(pc,1,Ru); */
201: PCMGSetRScale(pc, 1, su);
203: PCMGGetCoarseSolve(pc, &kspc);
204: KSPSetOperators(kspc, Auuc, Auuc);
205: KSPSetFromOptions(ksp);
207: if (analyze) {
208: mctx.dm = dmuu;
209: mctx.solRef = NULL; /* Reference solution not computed for u-u only */
210: mctx.solPrev = NULL; /* Populated automatically */
211: CreateNumericalReferenceSolution(Auu, fu, &mctx.solRefNum);
212: KSPMonitorSet(ksp, DMStagAnalysisKSPMonitor, &mctx, NULL);
213: }
215: /* Solve */
216: KSPSolve(ksp, fu, xu);
218: /* Clean up and finalize PETSc */
219: if (analyze) {
220: VecDestroy(&mctx.solPrev);
221: VecDestroy(&mctx.solRef);
222: VecDestroy(&mctx.solRefNum);
223: }
224: DMDestroy(&dmuu);
225: DMDestroy(&dmuuc);
226: VecDestroy(&su);
227: ISDestroy(&ispc);
228: ISDestroy(&ispf);
229: ISDestroy(&isuc);
230: ISDestroy(&isuf);
231: MatDestroy(&Ru);
232: MatDestroy(&IIu);
233: MatDestroy(&Auuc);
234: MatDestroy(&Auu);
235: VecDestroy(&xu);
236: VecDestroy(&fu);
237: KSPDestroy(&ksp);
238: DMDestroy(&dmSolc);
239: DMDestroy(&dmSol);
240: PetscFinalize();
241: return 0;
242: }
244: /*
245: Note: in this system all stencil coefficients which are not related to the Dirichlet boundary are scaled by dv = dx*dy.
246: This scaling is necessary for multigrid to converge.
247: */
248: static PetscErrorCode CreateSystem(DM dm, Mat *pA, Vec *pRhs)
249: {
250: PetscInt N[2], dof[3];
251: PetscBool isLastRankx, isLastRanky, isFirstRankx, isFirstRanky;
252: PetscInt ex, ey, startx, starty, nx, ny;
253: PetscInt iprev, icenter, inext;
254: Mat A;
255: Vec rhs;
256: PetscReal hx, hy, dv, bogusScale;
257: PetscScalar **cArrX, **cArrY;
258: PetscBool hasPressure;
262: /* Determine whether or not to create system including pressure dof (on elements) */
263: DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL);
265: hasPressure = (PetscBool)(dof[2] == 1);
267: bogusScale = 1.0;
268: PetscOptionsGetReal(NULL, NULL, "-bogus", &bogusScale, NULL); /* Use this to break MG (try small values) */
270: DMCreateMatrix(dm, pA);
271: A = *pA;
272: if (pRhs) {
273: DMCreateGlobalVector(dm, pRhs);
274: rhs = *pRhs;
275: } else {
276: rhs = NULL;
277: }
278: DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL);
279: DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL);
280: DMStagGetIsLastRank(dm, &isLastRankx, &isLastRanky, NULL);
281: DMStagGetIsFirstRank(dm, &isFirstRankx, &isFirstRanky, NULL);
282: hx = 1.0 / N[0];
283: hy = 1.0 / N[1];
284: dv = hx * hy;
285: DMStagGetProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL);
286: DMStagGetProductCoordinateLocationSlot(dm, LEFT, &iprev);
287: DMStagGetProductCoordinateLocationSlot(dm, RIGHT, &inext);
288: DMStagGetProductCoordinateLocationSlot(dm, ELEMENT, &icenter);
290: /* Loop over all local elements. Note that it may be more efficient in real
291: applications to loop over each boundary separately */
292: for (ey = starty; ey < starty + ny; ++ey) {
293: for (ex = startx; ex < startx + nx; ++ex) {
294: if (ex == N[0] - 1) {
295: /* Right Boundary velocity Dirichlet */
296: DMStagStencil row;
297: PetscScalar valRhs;
299: const PetscScalar valA = bogusScale * 1.0;
300: row.i = ex;
301: row.j = ey;
302: row.loc = RIGHT;
303: row.c = 0;
304: DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
305: if (rhs) {
306: valRhs = bogusScale * uxRef(cArrX[ex][inext], cArrY[ey][icenter]);
307: DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
308: }
309: }
310: if (ey == N[1] - 1) {
311: /* Top boundary velocity Dirichlet */
312: DMStagStencil row;
313: PetscScalar valRhs;
314: const PetscScalar valA = 1.0;
315: row.i = ex;
316: row.j = ey;
317: row.loc = UP;
318: row.c = 0;
319: DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
320: if (rhs) {
321: valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][inext]);
322: DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
323: }
324: }
326: if (ey == 0) {
327: /* Bottom boundary velocity Dirichlet */
328: DMStagStencil row;
329: PetscScalar valRhs;
330: const PetscScalar valA = 1.0;
331: row.i = ex;
332: row.j = ey;
333: row.loc = DOWN;
334: row.c = 0;
335: DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
336: if (rhs) {
337: valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
338: DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
339: }
340: } else {
341: /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
342: DMStagStencil row, col[7];
343: PetscScalar valA[7], valRhs;
344: PetscInt nEntries;
346: row.i = ex;
347: row.j = ey;
348: row.loc = DOWN;
349: row.c = 0;
350: if (ex == 0) {
351: nEntries = 4;
352: col[0].i = ex;
353: col[0].j = ey;
354: col[0].loc = DOWN;
355: col[0].c = 0;
356: valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
357: col[1].i = ex;
358: col[1].j = ey - 1;
359: col[1].loc = DOWN;
360: col[1].c = 0;
361: valA[1] = dv * 1.0 / (hy * hy);
362: col[2].i = ex;
363: col[2].j = ey + 1;
364: col[2].loc = DOWN;
365: col[2].c = 0;
366: valA[2] = dv * 1.0 / (hy * hy);
367: /* Missing left element */
368: col[3].i = ex + 1;
369: col[3].j = ey;
370: col[3].loc = DOWN;
371: col[3].c = 0;
372: valA[3] = dv * 1.0 / (hx * hx);
373: if (hasPressure) {
374: nEntries += 2;
375: col[4].i = ex;
376: col[4].j = ey - 1;
377: col[4].loc = ELEMENT;
378: col[4].c = 0;
379: valA[4] = dv * 1.0 / hy;
380: col[5].i = ex;
381: col[5].j = ey;
382: col[5].loc = ELEMENT;
383: col[5].c = 0;
384: valA[5] = -dv * 1.0 / hy;
385: }
386: } else if (ex == N[0] - 1) {
387: /* Right boundary y velocity stencil */
388: nEntries = 4;
389: col[0].i = ex;
390: col[0].j = ey;
391: col[0].loc = DOWN;
392: col[0].c = 0;
393: valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
394: col[1].i = ex;
395: col[1].j = ey - 1;
396: col[1].loc = DOWN;
397: col[1].c = 0;
398: valA[1] = dv * 1.0 / (hy * hy);
399: col[2].i = ex;
400: col[2].j = ey + 1;
401: col[2].loc = DOWN;
402: col[2].c = 0;
403: valA[2] = dv * 1.0 / (hy * hy);
404: col[3].i = ex - 1;
405: col[3].j = ey;
406: col[3].loc = DOWN;
407: col[3].c = 0;
408: valA[3] = dv * 1.0 / (hx * hx);
409: /* Missing right element */
410: if (hasPressure) {
411: nEntries += 2;
412: col[4].i = ex;
413: col[4].j = ey - 1;
414: col[4].loc = ELEMENT;
415: col[4].c = 0;
416: valA[4] = dv * 1.0 / hy;
417: col[5].i = ex;
418: col[5].j = ey;
419: col[5].loc = ELEMENT;
420: col[5].c = 0;
421: valA[5] = -dv * 1.0 / hy;
422: }
423: } else {
424: nEntries = 5;
425: col[0].i = ex;
426: col[0].j = ey;
427: col[0].loc = DOWN;
428: col[0].c = 0;
429: valA[0] = -dv * 2.0 / (hx * hx) - dv * 2.0 / (hy * hy);
430: col[1].i = ex;
431: col[1].j = ey - 1;
432: col[1].loc = DOWN;
433: col[1].c = 0;
434: valA[1] = dv * 1.0 / (hy * hy);
435: col[2].i = ex;
436: col[2].j = ey + 1;
437: col[2].loc = DOWN;
438: col[2].c = 0;
439: valA[2] = dv * 1.0 / (hy * hy);
440: col[3].i = ex - 1;
441: col[3].j = ey;
442: col[3].loc = DOWN;
443: col[3].c = 0;
444: valA[3] = dv * 1.0 / (hx * hx);
445: col[4].i = ex + 1;
446: col[4].j = ey;
447: col[4].loc = DOWN;
448: col[4].c = 0;
449: valA[4] = dv * 1.0 / (hx * hx);
450: if (hasPressure) {
451: nEntries += 2;
452: col[5].i = ex;
453: col[5].j = ey - 1;
454: col[5].loc = ELEMENT;
455: col[5].c = 0;
456: valA[5] = dv * 1.0 / hy;
457: col[6].i = ex;
458: col[6].j = ey;
459: col[6].loc = ELEMENT;
460: col[6].c = 0;
461: valA[6] = -dv * 1.0 / hy;
462: }
463: }
464: DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
465: if (rhs) {
466: valRhs = dv * fy(cArrX[ex][icenter], cArrY[ey][iprev]);
467: DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
468: }
469: }
471: if (ex == 0) {
472: /* Left velocity Dirichlet */
473: DMStagStencil row;
474: PetscScalar valRhs;
475: const PetscScalar valA = 1.0;
476: row.i = ex;
477: row.j = ey;
478: row.loc = LEFT;
479: row.c = 0;
480: DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
481: if (rhs) {
482: valRhs = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
483: DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
484: }
485: } else {
486: /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
487: DMStagStencil row, col[7];
488: PetscScalar valA[7], valRhs;
489: PetscInt nEntries;
490: row.i = ex;
491: row.j = ey;
492: row.loc = LEFT;
493: row.c = 0;
495: if (ey == 0) {
496: nEntries = 4;
497: col[0].i = ex;
498: col[0].j = ey;
499: col[0].loc = LEFT;
500: col[0].c = 0;
501: valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
502: /* missing term from element below */
503: col[1].i = ex;
504: col[1].j = ey + 1;
505: col[1].loc = LEFT;
506: col[1].c = 0;
507: valA[1] = dv * 1.0 / (hy * hy);
508: col[2].i = ex - 1;
509: col[2].j = ey;
510: col[2].loc = LEFT;
511: col[2].c = 0;
512: valA[2] = dv * 1.0 / (hx * hx);
513: col[3].i = ex + 1;
514: col[3].j = ey;
515: col[3].loc = LEFT;
516: col[3].c = 0;
517: valA[3] = dv * 1.0 / (hx * hx);
518: if (hasPressure) {
519: nEntries += 2;
520: col[4].i = ex - 1;
521: col[4].j = ey;
522: col[4].loc = ELEMENT;
523: col[4].c = 0;
524: valA[4] = dv * 1.0 / hx;
525: col[5].i = ex;
526: col[5].j = ey;
527: col[5].loc = ELEMENT;
528: col[5].c = 0;
529: valA[5] = -dv * 1.0 / hx;
530: }
531: } else if (ey == N[1] - 1) {
532: /* Top boundary x velocity stencil */
533: nEntries = 4;
534: row.i = ex;
535: row.j = ey;
536: row.loc = LEFT;
537: row.c = 0;
538: col[0].i = ex;
539: col[0].j = ey;
540: col[0].loc = LEFT;
541: col[0].c = 0;
542: valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
543: col[1].i = ex;
544: col[1].j = ey - 1;
545: col[1].loc = LEFT;
546: col[1].c = 0;
547: valA[1] = dv * 1.0 / (hy * hy);
548: /* Missing element above term */
549: col[2].i = ex - 1;
550: col[2].j = ey;
551: col[2].loc = LEFT;
552: col[2].c = 0;
553: valA[2] = dv * 1.0 / (hx * hx);
554: col[3].i = ex + 1;
555: col[3].j = ey;
556: col[3].loc = LEFT;
557: col[3].c = 0;
558: valA[3] = dv * 1.0 / (hx * hx);
559: if (hasPressure) {
560: nEntries += 2;
561: col[4].i = ex - 1;
562: col[4].j = ey;
563: col[4].loc = ELEMENT;
564: col[4].c = 0;
565: valA[4] = dv * 1.0 / hx;
566: col[5].i = ex;
567: col[5].j = ey;
568: col[5].loc = ELEMENT;
569: col[5].c = 0;
570: valA[5] = -dv * 1.0 / hx;
571: }
572: } else {
573: /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
574: nEntries = 5;
575: col[0].i = ex;
576: col[0].j = ey;
577: col[0].loc = LEFT;
578: col[0].c = 0;
579: valA[0] = -dv * 2.0 / (hx * hx) + -dv * 2.0 / (hy * hy);
580: col[1].i = ex;
581: col[1].j = ey - 1;
582: col[1].loc = LEFT;
583: col[1].c = 0;
584: valA[1] = dv * 1.0 / (hy * hy);
585: col[2].i = ex;
586: col[2].j = ey + 1;
587: col[2].loc = LEFT;
588: col[2].c = 0;
589: valA[2] = dv * 1.0 / (hy * hy);
590: col[3].i = ex - 1;
591: col[3].j = ey;
592: col[3].loc = LEFT;
593: col[3].c = 0;
594: valA[3] = dv * 1.0 / (hx * hx);
595: col[4].i = ex + 1;
596: col[4].j = ey;
597: col[4].loc = LEFT;
598: col[4].c = 0;
599: valA[4] = dv * 1.0 / (hx * hx);
600: if (hasPressure) {
601: nEntries += 2;
602: col[5].i = ex - 1;
603: col[5].j = ey;
604: col[5].loc = ELEMENT;
605: col[5].c = 0;
606: valA[5] = dv * 1.0 / hx;
607: col[6].i = ex;
608: col[6].j = ey;
609: col[6].loc = ELEMENT;
610: col[6].c = 0;
611: valA[6] = -dv * 1.0 / hx;
612: }
613: }
614: DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
615: if (rhs) {
616: valRhs = dv * fx(cArrX[ex][iprev], cArrY[ey][icenter]);
617: DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
618: }
619: }
621: /* P equation : u_x + v_y = g
622: Note that this includes an explicit zero on the diagonal. This is only needed for
623: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
624: if (hasPressure) {
625: DMStagStencil row, col[5];
626: PetscScalar valA[5], valRhs;
628: row.i = ex;
629: row.j = ey;
630: row.loc = ELEMENT;
631: row.c = 0;
632: /* Note: the scaling by dv here may not be optimal (but this test isn't concerned with these equations) */
633: col[0].i = ex;
634: col[0].j = ey;
635: col[0].loc = LEFT;
636: col[0].c = 0;
637: valA[0] = -dv * 1.0 / hx;
638: col[1].i = ex;
639: col[1].j = ey;
640: col[1].loc = RIGHT;
641: col[1].c = 0;
642: valA[1] = dv * 1.0 / hx;
643: col[2].i = ex;
644: col[2].j = ey;
645: col[2].loc = DOWN;
646: col[2].c = 0;
647: valA[2] = -dv * 1.0 / hy;
648: col[3].i = ex;
649: col[3].j = ey;
650: col[3].loc = UP;
651: col[3].c = 0;
652: valA[3] = dv * 1.0 / hy;
653: col[4] = row;
654: valA[4] = dv * 0.0;
655: DMStagMatSetValuesStencil(dm, A, 1, &row, 5, col, valA, INSERT_VALUES);
656: if (rhs) {
657: valRhs = dv * g(cArrX[ex][icenter], cArrY[ey][icenter]);
658: DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES);
659: }
660: }
661: }
662: }
663: DMStagRestoreProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL);
664: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
665: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
666: if (rhs) {
667: VecAssemblyBegin(rhs);
668: VecAssemblyEnd(rhs);
669: }
671: return 0;
672: }
674: /* A custom monitor function for analysis purposes. Computes and dumps
675: residuals and errors for each KSP iteration */
676: PetscErrorCode DMStagAnalysisKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
677: {
678: DM dm;
679: Vec r, sol;
680: DMStagAnalysisKSPMonitorContext *ctx = (DMStagAnalysisKSPMonitorContext *)mctx;
683: KSPBuildSolution(ksp, NULL, &sol); /* don't destroy sol */
684: KSPBuildResidual(ksp, NULL, NULL, &r);
685: dm = ctx->dm; /* Would typically get this with VecGetDM(), KSPGetDM() */
687: PetscPrintf(PetscObjectComm((PetscObject)dm), " *** DMStag Analysis KSP Monitor (it. %" PetscInt_FMT ") ***\n", it);
688: PetscPrintf(PetscObjectComm((PetscObject)dm), " Residual Norm: %g\n", (double)rnorm);
689: PetscPrintf(PetscObjectComm((PetscObject)dm), " Dumping files...\n");
691: /* Note: these blocks are almost entirely duplicated */
692: {
693: const DMStagStencilLocation loc = DMSTAG_LEFT;
694: const PetscInt c = 0;
695: Vec vec;
696: DM da;
697: PetscViewer viewer;
698: char filename[PETSC_MAX_PATH_LEN];
700: DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec);
701: PetscSNPrintf(filename, sizeof(filename), "res_vx_%" PetscInt_FMT ".vtr", it);
702: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
703: VecView(vec, viewer);
704: PetscViewerDestroy(&viewer);
705: VecDestroy(&vec);
706: DMDestroy(&da);
707: }
709: {
710: const DMStagStencilLocation loc = DMSTAG_DOWN;
711: const PetscInt c = 0;
712: Vec vec;
713: DM da;
714: PetscViewer viewer;
715: char filename[PETSC_MAX_PATH_LEN];
717: DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec);
718: PetscSNPrintf(filename, sizeof(filename), "res_vy_%" PetscInt_FMT ".vtr", it);
719: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
720: VecView(vec, viewer);
721: PetscViewerDestroy(&viewer);
722: VecDestroy(&vec);
723: DMDestroy(&da);
724: }
726: if (ctx->solRef) {
727: Vec e;
728: VecDuplicate(ctx->solRef, &e);
729: VecCopy(ctx->solRef, e);
730: VecAXPY(e, -1.0, sol);
732: {
733: const DMStagStencilLocation loc = DMSTAG_LEFT;
734: const PetscInt c = 0;
735: Vec vec;
736: DM da;
737: PetscViewer viewer;
738: char filename[PETSC_MAX_PATH_LEN];
740: DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
741: PetscSNPrintf(filename, sizeof(filename), "err_vx_%" PetscInt_FMT ".vtr", it);
742: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
743: VecView(vec, viewer);
744: PetscViewerDestroy(&viewer);
745: VecDestroy(&vec);
746: DMDestroy(&da);
747: }
749: {
750: const DMStagStencilLocation loc = DMSTAG_DOWN;
751: const PetscInt c = 0;
752: Vec vec;
753: DM da;
754: PetscViewer viewer;
755: char filename[PETSC_MAX_PATH_LEN];
757: DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
758: PetscSNPrintf(filename, sizeof(filename), "err_vy_%" PetscInt_FMT ".vtr", it);
759: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
760: VecView(vec, viewer);
761: PetscViewerDestroy(&viewer);
762: VecDestroy(&vec);
763: DMDestroy(&da);
764: }
766: VecDestroy(&e);
767: }
769: /* Repeat error computations wrt an "exact" solution to the discrete equations */
770: if (ctx->solRefNum) {
771: Vec e;
772: VecDuplicate(ctx->solRefNum, &e);
773: VecCopy(ctx->solRefNum, e);
774: VecAXPY(e, -1.0, sol);
776: {
777: const DMStagStencilLocation loc = DMSTAG_LEFT;
778: const PetscInt c = 0;
779: Vec vec;
780: DM da;
781: PetscViewer viewer;
782: char filename[PETSC_MAX_PATH_LEN];
784: DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
785: PetscSNPrintf(filename, sizeof(filename), "err_num_vx_%" PetscInt_FMT ".vtr", it);
786: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
787: VecView(vec, viewer);
788: PetscViewerDestroy(&viewer);
789: VecDestroy(&vec);
790: DMDestroy(&da);
791: }
793: {
794: const DMStagStencilLocation loc = DMSTAG_DOWN;
795: const PetscInt c = 0;
796: Vec vec;
797: DM da;
798: PetscViewer viewer;
799: char filename[PETSC_MAX_PATH_LEN];
801: DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec);
802: PetscSNPrintf(filename, sizeof(filename), "err_num_vy_%" PetscInt_FMT ".vtr", it);
803: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
804: VecView(vec, viewer);
805: PetscViewerDestroy(&viewer);
806: VecDestroy(&vec);
807: DMDestroy(&da);
808: }
810: VecDestroy(&e);
811: }
813: /* Step */
814: if (!ctx->solPrev) {
815: VecDuplicate(sol, &ctx->solPrev);
816: } else {
817: Vec diff;
818: VecDuplicate(sol, &diff);
819: VecCopy(sol, diff);
820: VecAXPY(diff, -1.0, ctx->solPrev);
821: {
822: const DMStagStencilLocation loc = DMSTAG_LEFT;
823: const PetscInt c = 0;
824: Vec vec;
825: DM da;
826: PetscViewer viewer;
827: char filename[PETSC_MAX_PATH_LEN];
829: DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec);
830: PetscSNPrintf(filename, sizeof(filename), "diff_vx_%" PetscInt_FMT ".vtr", it);
831: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
832: VecView(vec, viewer);
833: PetscViewerDestroy(&viewer);
834: VecDestroy(&vec);
835: DMDestroy(&da);
836: }
838: {
839: const DMStagStencilLocation loc = DMSTAG_DOWN;
840: const PetscInt c = 0;
841: Vec vec;
842: DM da;
843: PetscViewer viewer;
844: char filename[PETSC_MAX_PATH_LEN];
846: DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec);
847: PetscSNPrintf(filename, sizeof(filename), "diff_vy_%" PetscInt_FMT ".vtr", it);
848: PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer);
849: VecView(vec, viewer);
850: PetscViewerDestroy(&viewer);
851: VecDestroy(&vec);
852: DMDestroy(&da);
853: }
854: VecDestroy(&diff);
855: }
856: VecCopy(sol, ctx->solPrev);
858: VecDestroy(&r);
859: return 0;
860: }
862: /* Use a direct solver to create an "exact" solution to the discrete system
863: useful for testing solvers (in that it doesn't include discretization error) */
864: static PetscErrorCode CreateNumericalReferenceSolution(Mat A, Vec rhs, Vec *px)
865: {
866: KSP ksp;
867: PC pc;
868: Vec x;
871: KSPCreate(PetscObjectComm((PetscObject)A), &ksp);
872: KSPSetType(ksp, KSPPREONLY);
873: KSPGetPC(ksp, &pc);
874: PCSetType(pc, PCLU);
875: PCFactorSetMatSolverType(pc, MATSOLVERUMFPACK);
876: KSPSetOptionsPrefix(ksp, "numref_");
877: KSPSetFromOptions(ksp);
878: KSPSetOperators(ksp, A, A);
879: VecDuplicate(rhs, px);
880: x = *px;
881: KSPSolve(ksp, rhs, x);
882: KSPDestroy(&ksp);
883: return 0;
884: }
886: /*TEST
888: test:
889: suffix: gmg_1
890: nsize: 1
891: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason
893: test:
894: suffix: gmg_1_b
895: nsize: 1
896: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false
898: test:
899: suffix: gmg_1_c
900: nsize: 1
901: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false
903: test:
904: suffix: gmg_1_d
905: nsize: 1
906: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false
908: test:
909: suffix: gmg_1_bigger
910: requires: !complex
911: nsize: 1
912: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32
914: test:
915: suffix: gmg_8
916: requires: !single !complex
917: nsize: 8
918: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason
920: test:
921: suffix: gmg_8_b
922: nsize: 8
923: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false
925: test:
926: suffix: gmg_8_c
927: nsize: 8
928: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false
930: test:
931: suffix: gmg_8_d
932: nsize: 8
933: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false
935: test:
936: suffix: gmg_8_galerkin
937: nsize: 8
938: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false -pc_mg_galerkin
940: test:
941: suffix: gmg_8_bigger
942: requires: !complex
943: nsize: 8
944: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32
946: TEST*/