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 */
94: PetscFunctionBeginUser;
95: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
96: analyze = PETSC_FALSE;
97: PetscCall(PetscOptionsGetBool(NULL, NULL, "-analyze", &analyze, NULL));
98: extractTransferOperators = PETSC_TRUE;
99: PetscCall(PetscOptionsGetBool(NULL, NULL, "-extractTransferOperators", &extractTransferOperators, NULL));
100: extractSystem = PETSC_TRUE;
101: PetscCall(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: PetscCall(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: PetscCall(DMSetFromOptions(dmSol));
109: PetscCall(DMSetUp(dmSol));
110: }
112: /* Create coarse DMStag (which we may or may not directly use) */
113: PetscCall(DMCoarsen(dmSol, MPI_COMM_NULL, &dmSolc));
114: PetscCall(DMSetUp(dmSolc));
116: /* Create compatible DMStags with only velocity dof (which we may or may not
117: directly use) */
118: PetscCall(DMStagCreateCompatibleDMStag(dmSol, 0, 1, 0, 0, &dmuu)); /* vel-only */
119: PetscCall(DMSetUp(dmuu));
120: PetscCall(DMCoarsen(dmuu, MPI_COMM_NULL, &dmuuc));
121: PetscCall(DMSetUp(dmuuc));
123: /* Define uniform coordinates as a product of 1D arrays */
124: PetscCall(DMStagSetUniformCoordinatesProduct(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
125: PetscCall(DMStagSetUniformCoordinatesProduct(dmSolc, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
126: PetscCall(DMStagSetUniformCoordinatesProduct(dmuu, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
127: PetscCall(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: PetscCall(DMStagCreateISFromStencils(dmSol, 2, stencil_set, &isuf));
142: PetscCall(DMStagCreateISFromStencils(dmSol, 1, &stencil_set[2], &ispf));
143: PetscCall(DMStagCreateISFromStencils(dmSolc, 2, stencil_set, &isuc));
144: PetscCall(DMStagCreateISFromStencils(dmSolc, 1, &stencil_set[2], &ispc));
146: /* Assemble velocity-velocity system */
147: if (extractSystem) {
148: Mat A, Ac;
149: Vec tmp, rhs;
151: PetscCall(CreateSystem(dmSol, &A, &rhs));
152: PetscCall(CreateSystem(dmSolc, &Ac, NULL));
153: PetscCall(MatCreateSubMatrix(A, isuf, isuf, MAT_INITIAL_MATRIX, &Auu));
154: PetscCall(MatCreateSubMatrix(Ac, isuc, isuc, MAT_INITIAL_MATRIX, &Auuc));
155: PetscCall(MatCreateVecs(Auu, &xu, &fu));
156: PetscCall(VecGetSubVector(rhs, isuf, &tmp));
157: PetscCall(VecCopy(tmp, fu));
158: PetscCall(VecRestoreSubVector(rhs, isuf, &tmp));
159: PetscCall(MatDestroy(&Ac));
160: PetscCall(VecDestroy(&rhs));
161: PetscCall(MatDestroy(&A));
162: } else {
163: PetscCall(CreateSystem(dmuu, &Auu, &fu));
164: PetscCall(CreateSystem(dmuuc, &Auuc, NULL));
165: PetscCall(MatCreateVecs(Auu, &xu, NULL));
166: }
167: PetscCall(PetscObjectSetName((PetscObject)Auu, "Auu"));
168: PetscCall(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: PetscCall(DMCreateInterpolation(dmSolc, dmSol, &II, NULL));
176: PetscCall(DMCreateRestriction(dmSolc, dmSol, &R));
177: PetscCall(MatCreateSubMatrix(II, isuf, isuc, MAT_INITIAL_MATRIX, &IIu));
178: PetscCall(MatCreateSubMatrix(R, isuc, isuf, MAT_INITIAL_MATRIX, &Ru));
179: PetscCall(DMCreateInterpolationScale(dmSolc, dmSol, II, &s));
180: PetscCall(MatCreateVecs(IIu, &su, NULL));
181: PetscCall(VecGetSubVector(s, isuc, &tmp));
182: PetscCall(VecCopy(tmp, su));
183: PetscCall(VecRestoreSubVector(s, isuc, &tmp));
184: PetscCall(MatDestroy(&R));
185: PetscCall(MatDestroy(&II));
186: PetscCall(VecDestroy(&s));
187: } else {
188: PetscCall(DMCreateInterpolation(dmuuc, dmuu, &IIu, NULL));
189: PetscCall(DMCreateInterpolationScale(dmuuc, dmuu, IIu, &su));
190: PetscCall(DMCreateRestriction(dmuuc, dmuu, &Ru));
191: }
193: /* Create and configure solver */
194: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
195: PetscCall(KSPSetOperators(ksp, Auu, Auu));
196: PetscCall(KSPGetPC(ksp, &pc));
197: PetscCall(PCSetType(pc, PCMG));
198: PetscCall(PCMGSetLevels(pc, 2, NULL));
199: PetscCall(PCMGSetInterpolation(pc, 1, IIu));
200: /* PetscCall(PCMGSetRestriction(pc,1,Ru)); */
201: PetscCall(PCMGSetRScale(pc, 1, su));
203: PetscCall(PCMGGetCoarseSolve(pc, &kspc));
204: PetscCall(KSPSetOperators(kspc, Auuc, Auuc));
205: PetscCall(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: PetscCall(CreateNumericalReferenceSolution(Auu, fu, &mctx.solRefNum));
212: PetscCall(KSPMonitorSet(ksp, DMStagAnalysisKSPMonitor, &mctx, NULL));
213: }
215: /* Solve */
216: PetscCall(KSPSolve(ksp, fu, xu));
218: /* Clean up and finalize PETSc */
219: if (analyze) {
220: PetscCall(VecDestroy(&mctx.solPrev));
221: PetscCall(VecDestroy(&mctx.solRef));
222: PetscCall(VecDestroy(&mctx.solRefNum));
223: }
224: PetscCall(DMDestroy(&dmuu));
225: PetscCall(DMDestroy(&dmuuc));
226: PetscCall(VecDestroy(&su));
227: PetscCall(ISDestroy(&ispc));
228: PetscCall(ISDestroy(&ispf));
229: PetscCall(ISDestroy(&isuc));
230: PetscCall(ISDestroy(&isuf));
231: PetscCall(MatDestroy(&Ru));
232: PetscCall(MatDestroy(&IIu));
233: PetscCall(MatDestroy(&Auuc));
234: PetscCall(MatDestroy(&Auu));
235: PetscCall(VecDestroy(&xu));
236: PetscCall(VecDestroy(&fu));
237: PetscCall(KSPDestroy(&ksp));
238: PetscCall(DMDestroy(&dmSolc));
239: PetscCall(DMDestroy(&dmSol));
240: PetscCall(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;
260: PetscFunctionBeginUser;
261: /* Determine whether or not to create system including pressure dof (on elements) */
262: PetscCall(DMStagGetDOF(dm, &dof[0], &dof[1], &dof[2], NULL));
263: PetscCheck(dof[0] == 0 && dof[1] == 1 && (dof[2] == 1 || dof[2] == 0), PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "CreateSystem only implemented for velocity-only or velocity+pressure grids");
264: hasPressure = (PetscBool)(dof[2] == 1);
266: bogusScale = 1.0;
267: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bogus", &bogusScale, NULL)); /* Use this to break MG (try small values) */
269: PetscCall(DMCreateMatrix(dm, pA));
270: A = *pA;
271: if (pRhs) {
272: PetscCall(DMCreateGlobalVector(dm, pRhs));
273: rhs = *pRhs;
274: } else {
275: rhs = NULL;
276: }
277: PetscCall(DMStagGetCorners(dm, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
278: PetscCall(DMStagGetGlobalSizes(dm, &N[0], &N[1], NULL));
279: PetscCall(DMStagGetIsLastRank(dm, &isLastRankx, &isLastRanky, NULL));
280: PetscCall(DMStagGetIsFirstRank(dm, &isFirstRankx, &isFirstRanky, NULL));
281: hx = 1.0 / N[0];
282: hy = 1.0 / N[1];
283: dv = hx * hy;
284: PetscCall(DMStagGetProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL));
285: PetscCall(DMStagGetProductCoordinateLocationSlot(dm, LEFT, &iprev));
286: PetscCall(DMStagGetProductCoordinateLocationSlot(dm, RIGHT, &inext));
287: PetscCall(DMStagGetProductCoordinateLocationSlot(dm, ELEMENT, &icenter));
289: /* Loop over all local elements. Note that it may be more efficient in real
290: applications to loop over each boundary separately */
291: for (ey = starty; ey < starty + ny; ++ey) {
292: for (ex = startx; ex < startx + nx; ++ex) {
293: if (ex == N[0] - 1) {
294: /* Right Boundary velocity Dirichlet */
295: DMStagStencil row;
296: PetscScalar valRhs;
298: const PetscScalar valA = bogusScale * 1.0;
299: row.i = ex;
300: row.j = ey;
301: row.loc = RIGHT;
302: row.c = 0;
303: PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
304: if (rhs) {
305: valRhs = bogusScale * uxRef(cArrX[ex][inext], cArrY[ey][icenter]);
306: PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
307: }
308: }
309: if (ey == N[1] - 1) {
310: /* Top boundary velocity Dirichlet */
311: DMStagStencil row;
312: PetscScalar valRhs;
313: const PetscScalar valA = 1.0;
314: row.i = ex;
315: row.j = ey;
316: row.loc = UP;
317: row.c = 0;
318: PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
319: if (rhs) {
320: valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][inext]);
321: PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
322: }
323: }
325: if (ey == 0) {
326: /* Bottom boundary velocity Dirichlet */
327: DMStagStencil row;
328: PetscScalar valRhs;
329: const PetscScalar valA = 1.0;
330: row.i = ex;
331: row.j = ey;
332: row.loc = DOWN;
333: row.c = 0;
334: PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
335: if (rhs) {
336: valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
337: PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
338: }
339: } else {
340: /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
341: DMStagStencil row, col[7];
342: PetscScalar valA[7], valRhs;
343: PetscInt nEntries;
345: row.i = ex;
346: row.j = ey;
347: row.loc = DOWN;
348: row.c = 0;
349: if (ex == 0) {
350: nEntries = 4;
351: col[0].i = ex;
352: col[0].j = ey;
353: col[0].loc = DOWN;
354: col[0].c = 0;
355: valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
356: col[1].i = ex;
357: col[1].j = ey - 1;
358: col[1].loc = DOWN;
359: col[1].c = 0;
360: valA[1] = dv * 1.0 / (hy * hy);
361: col[2].i = ex;
362: col[2].j = ey + 1;
363: col[2].loc = DOWN;
364: col[2].c = 0;
365: valA[2] = dv * 1.0 / (hy * hy);
366: /* Missing left element */
367: col[3].i = ex + 1;
368: col[3].j = ey;
369: col[3].loc = DOWN;
370: col[3].c = 0;
371: valA[3] = dv * 1.0 / (hx * hx);
372: if (hasPressure) {
373: nEntries += 2;
374: col[4].i = ex;
375: col[4].j = ey - 1;
376: col[4].loc = ELEMENT;
377: col[4].c = 0;
378: valA[4] = dv * 1.0 / hy;
379: col[5].i = ex;
380: col[5].j = ey;
381: col[5].loc = ELEMENT;
382: col[5].c = 0;
383: valA[5] = -dv * 1.0 / hy;
384: }
385: } else if (ex == N[0] - 1) {
386: /* Right boundary y velocity stencil */
387: nEntries = 4;
388: col[0].i = ex;
389: col[0].j = ey;
390: col[0].loc = DOWN;
391: col[0].c = 0;
392: valA[0] = -dv * 1.0 / (hx * hx) - dv * 2.0 / (hy * hy);
393: col[1].i = ex;
394: col[1].j = ey - 1;
395: col[1].loc = DOWN;
396: col[1].c = 0;
397: valA[1] = dv * 1.0 / (hy * hy);
398: col[2].i = ex;
399: col[2].j = ey + 1;
400: col[2].loc = DOWN;
401: col[2].c = 0;
402: valA[2] = dv * 1.0 / (hy * hy);
403: col[3].i = ex - 1;
404: col[3].j = ey;
405: col[3].loc = DOWN;
406: col[3].c = 0;
407: valA[3] = dv * 1.0 / (hx * hx);
408: /* Missing right element */
409: if (hasPressure) {
410: nEntries += 2;
411: col[4].i = ex;
412: col[4].j = ey - 1;
413: col[4].loc = ELEMENT;
414: col[4].c = 0;
415: valA[4] = dv * 1.0 / hy;
416: col[5].i = ex;
417: col[5].j = ey;
418: col[5].loc = ELEMENT;
419: col[5].c = 0;
420: valA[5] = -dv * 1.0 / hy;
421: }
422: } else {
423: nEntries = 5;
424: col[0].i = ex;
425: col[0].j = ey;
426: col[0].loc = DOWN;
427: col[0].c = 0;
428: valA[0] = -dv * 2.0 / (hx * hx) - dv * 2.0 / (hy * hy);
429: col[1].i = ex;
430: col[1].j = ey - 1;
431: col[1].loc = DOWN;
432: col[1].c = 0;
433: valA[1] = dv * 1.0 / (hy * hy);
434: col[2].i = ex;
435: col[2].j = ey + 1;
436: col[2].loc = DOWN;
437: col[2].c = 0;
438: valA[2] = dv * 1.0 / (hy * hy);
439: col[3].i = ex - 1;
440: col[3].j = ey;
441: col[3].loc = DOWN;
442: col[3].c = 0;
443: valA[3] = dv * 1.0 / (hx * hx);
444: col[4].i = ex + 1;
445: col[4].j = ey;
446: col[4].loc = DOWN;
447: col[4].c = 0;
448: valA[4] = dv * 1.0 / (hx * hx);
449: if (hasPressure) {
450: nEntries += 2;
451: col[5].i = ex;
452: col[5].j = ey - 1;
453: col[5].loc = ELEMENT;
454: col[5].c = 0;
455: valA[5] = dv * 1.0 / hy;
456: col[6].i = ex;
457: col[6].j = ey;
458: col[6].loc = ELEMENT;
459: col[6].c = 0;
460: valA[6] = -dv * 1.0 / hy;
461: }
462: }
463: PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
464: if (rhs) {
465: valRhs = dv * fy(cArrX[ex][icenter], cArrY[ey][iprev]);
466: PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
467: }
468: }
470: if (ex == 0) {
471: /* Left velocity Dirichlet */
472: DMStagStencil row;
473: PetscScalar valRhs;
474: const PetscScalar valA = 1.0;
475: row.i = ex;
476: row.j = ey;
477: row.loc = LEFT;
478: row.c = 0;
479: PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
480: if (rhs) {
481: valRhs = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
482: PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
483: }
484: } else {
485: /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
486: DMStagStencil row, col[7];
487: PetscScalar valA[7], valRhs;
488: PetscInt nEntries;
489: row.i = ex;
490: row.j = ey;
491: row.loc = LEFT;
492: row.c = 0;
494: if (ey == 0) {
495: nEntries = 4;
496: col[0].i = ex;
497: col[0].j = ey;
498: col[0].loc = LEFT;
499: col[0].c = 0;
500: valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
501: /* missing term from element below */
502: col[1].i = ex;
503: col[1].j = ey + 1;
504: col[1].loc = LEFT;
505: col[1].c = 0;
506: valA[1] = dv * 1.0 / (hy * hy);
507: col[2].i = ex - 1;
508: col[2].j = ey;
509: col[2].loc = LEFT;
510: col[2].c = 0;
511: valA[2] = dv * 1.0 / (hx * hx);
512: col[3].i = ex + 1;
513: col[3].j = ey;
514: col[3].loc = LEFT;
515: col[3].c = 0;
516: valA[3] = dv * 1.0 / (hx * hx);
517: if (hasPressure) {
518: nEntries += 2;
519: col[4].i = ex - 1;
520: col[4].j = ey;
521: col[4].loc = ELEMENT;
522: col[4].c = 0;
523: valA[4] = dv * 1.0 / hx;
524: col[5].i = ex;
525: col[5].j = ey;
526: col[5].loc = ELEMENT;
527: col[5].c = 0;
528: valA[5] = -dv * 1.0 / hx;
529: }
530: } else if (ey == N[1] - 1) {
531: /* Top boundary x velocity stencil */
532: nEntries = 4;
533: row.i = ex;
534: row.j = ey;
535: row.loc = LEFT;
536: row.c = 0;
537: col[0].i = ex;
538: col[0].j = ey;
539: col[0].loc = LEFT;
540: col[0].c = 0;
541: valA[0] = -dv * 2.0 / (hx * hx) - dv * 1.0 / (hy * hy);
542: col[1].i = ex;
543: col[1].j = ey - 1;
544: col[1].loc = LEFT;
545: col[1].c = 0;
546: valA[1] = dv * 1.0 / (hy * hy);
547: /* Missing element above term */
548: col[2].i = ex - 1;
549: col[2].j = ey;
550: col[2].loc = LEFT;
551: col[2].c = 0;
552: valA[2] = dv * 1.0 / (hx * hx);
553: col[3].i = ex + 1;
554: col[3].j = ey;
555: col[3].loc = LEFT;
556: col[3].c = 0;
557: valA[3] = dv * 1.0 / (hx * hx);
558: if (hasPressure) {
559: nEntries += 2;
560: col[4].i = ex - 1;
561: col[4].j = ey;
562: col[4].loc = ELEMENT;
563: col[4].c = 0;
564: valA[4] = dv * 1.0 / hx;
565: col[5].i = ex;
566: col[5].j = ey;
567: col[5].loc = ELEMENT;
568: col[5].c = 0;
569: valA[5] = -dv * 1.0 / hx;
570: }
571: } else {
572: /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
573: nEntries = 5;
574: col[0].i = ex;
575: col[0].j = ey;
576: col[0].loc = LEFT;
577: col[0].c = 0;
578: valA[0] = -dv * 2.0 / (hx * hx) + -dv * 2.0 / (hy * hy);
579: col[1].i = ex;
580: col[1].j = ey - 1;
581: col[1].loc = LEFT;
582: col[1].c = 0;
583: valA[1] = dv * 1.0 / (hy * hy);
584: col[2].i = ex;
585: col[2].j = ey + 1;
586: col[2].loc = LEFT;
587: col[2].c = 0;
588: valA[2] = dv * 1.0 / (hy * hy);
589: col[3].i = ex - 1;
590: col[3].j = ey;
591: col[3].loc = LEFT;
592: col[3].c = 0;
593: valA[3] = dv * 1.0 / (hx * hx);
594: col[4].i = ex + 1;
595: col[4].j = ey;
596: col[4].loc = LEFT;
597: col[4].c = 0;
598: valA[4] = dv * 1.0 / (hx * hx);
599: if (hasPressure) {
600: nEntries += 2;
601: col[5].i = ex - 1;
602: col[5].j = ey;
603: col[5].loc = ELEMENT;
604: col[5].c = 0;
605: valA[5] = dv * 1.0 / hx;
606: col[6].i = ex;
607: col[6].j = ey;
608: col[6].loc = ELEMENT;
609: col[6].c = 0;
610: valA[6] = -dv * 1.0 / hx;
611: }
612: }
613: PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
614: if (rhs) {
615: valRhs = dv * fx(cArrX[ex][iprev], cArrY[ey][icenter]);
616: PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
617: }
618: }
620: /* P equation : u_x + v_y = g
621: Note that this includes an explicit zero on the diagonal. This is only needed for
622: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
623: if (hasPressure) {
624: DMStagStencil row, col[5];
625: PetscScalar valA[5], valRhs;
627: row.i = ex;
628: row.j = ey;
629: row.loc = ELEMENT;
630: row.c = 0;
631: /* Note: the scaling by dv here may not be optimal (but this test isn't concerned with these equations) */
632: col[0].i = ex;
633: col[0].j = ey;
634: col[0].loc = LEFT;
635: col[0].c = 0;
636: valA[0] = -dv * 1.0 / hx;
637: col[1].i = ex;
638: col[1].j = ey;
639: col[1].loc = RIGHT;
640: col[1].c = 0;
641: valA[1] = dv * 1.0 / hx;
642: col[2].i = ex;
643: col[2].j = ey;
644: col[2].loc = DOWN;
645: col[2].c = 0;
646: valA[2] = -dv * 1.0 / hy;
647: col[3].i = ex;
648: col[3].j = ey;
649: col[3].loc = UP;
650: col[3].c = 0;
651: valA[3] = dv * 1.0 / hy;
652: col[4] = row;
653: valA[4] = dv * 0.0;
654: PetscCall(DMStagMatSetValuesStencil(dm, A, 1, &row, 5, col, valA, INSERT_VALUES));
655: if (rhs) {
656: valRhs = dv * g(cArrX[ex][icenter], cArrY[ey][icenter]);
657: PetscCall(DMStagVecSetValuesStencil(dm, rhs, 1, &row, &valRhs, INSERT_VALUES));
658: }
659: }
660: }
661: }
662: PetscCall(DMStagRestoreProductCoordinateArraysRead(dm, &cArrX, &cArrY, NULL));
663: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
664: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
665: if (rhs) {
666: PetscCall(VecAssemblyBegin(rhs));
667: PetscCall(VecAssemblyEnd(rhs));
668: }
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /* A custom monitor function for analysis purposes. Computes and dumps
673: residuals and errors for each KSP iteration */
674: PetscErrorCode DMStagAnalysisKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
675: {
676: DM dm;
677: Vec r, sol;
678: DMStagAnalysisKSPMonitorContext *ctx = (DMStagAnalysisKSPMonitorContext *)mctx;
680: PetscFunctionBeginUser;
681: PetscCall(KSPBuildSolution(ksp, NULL, &sol)); /* don't destroy sol */
682: PetscCall(KSPBuildResidual(ksp, NULL, NULL, &r));
683: dm = ctx->dm; /* Would typically get this with VecGetDM(), KSPGetDM() */
684: PetscCheck(dm, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Analaysis monitor requires a DM which is properly associated with the solution Vec");
685: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " *** DMStag Analysis KSP Monitor (it. %" PetscInt_FMT ") ***\n", it));
686: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " Residual Norm: %g\n", (double)rnorm));
687: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " Dumping files...\n"));
689: /* Note: these blocks are almost entirely duplicated */
690: {
691: const DMStagStencilLocation loc = DMSTAG_LEFT;
692: const PetscInt c = 0;
693: Vec vec;
694: DM da;
695: PetscViewer viewer;
696: char filename[PETSC_MAX_PATH_LEN];
698: PetscCall(DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec));
699: PetscCall(PetscSNPrintf(filename, sizeof(filename), "res_vx_%" PetscInt_FMT ".vtr", it));
700: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
701: PetscCall(VecView(vec, viewer));
702: PetscCall(PetscViewerDestroy(&viewer));
703: PetscCall(VecDestroy(&vec));
704: PetscCall(DMDestroy(&da));
705: }
707: {
708: const DMStagStencilLocation loc = DMSTAG_DOWN;
709: const PetscInt c = 0;
710: Vec vec;
711: DM da;
712: PetscViewer viewer;
713: char filename[PETSC_MAX_PATH_LEN];
715: PetscCall(DMStagVecSplitToDMDA(dm, r, loc, c, &da, &vec));
716: PetscCall(PetscSNPrintf(filename, sizeof(filename), "res_vy_%" PetscInt_FMT ".vtr", it));
717: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
718: PetscCall(VecView(vec, viewer));
719: PetscCall(PetscViewerDestroy(&viewer));
720: PetscCall(VecDestroy(&vec));
721: PetscCall(DMDestroy(&da));
722: }
724: if (ctx->solRef) {
725: Vec e;
726: PetscCall(VecDuplicate(ctx->solRef, &e));
727: PetscCall(VecCopy(ctx->solRef, e));
728: PetscCall(VecAXPY(e, -1.0, sol));
730: {
731: const DMStagStencilLocation loc = DMSTAG_LEFT;
732: const PetscInt c = 0;
733: Vec vec;
734: DM da;
735: PetscViewer viewer;
736: char filename[PETSC_MAX_PATH_LEN];
738: PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
739: PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_vx_%" PetscInt_FMT ".vtr", it));
740: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
741: PetscCall(VecView(vec, viewer));
742: PetscCall(PetscViewerDestroy(&viewer));
743: PetscCall(VecDestroy(&vec));
744: PetscCall(DMDestroy(&da));
745: }
747: {
748: const DMStagStencilLocation loc = DMSTAG_DOWN;
749: const PetscInt c = 0;
750: Vec vec;
751: DM da;
752: PetscViewer viewer;
753: char filename[PETSC_MAX_PATH_LEN];
755: PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
756: PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_vy_%" PetscInt_FMT ".vtr", it));
757: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
758: PetscCall(VecView(vec, viewer));
759: PetscCall(PetscViewerDestroy(&viewer));
760: PetscCall(VecDestroy(&vec));
761: PetscCall(DMDestroy(&da));
762: }
764: PetscCall(VecDestroy(&e));
765: }
767: /* Repeat error computations wrt an "exact" solution to the discrete equations */
768: if (ctx->solRefNum) {
769: Vec e;
770: PetscCall(VecDuplicate(ctx->solRefNum, &e));
771: PetscCall(VecCopy(ctx->solRefNum, e));
772: PetscCall(VecAXPY(e, -1.0, sol));
774: {
775: const DMStagStencilLocation loc = DMSTAG_LEFT;
776: const PetscInt c = 0;
777: Vec vec;
778: DM da;
779: PetscViewer viewer;
780: char filename[PETSC_MAX_PATH_LEN];
782: PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
783: PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_num_vx_%" PetscInt_FMT ".vtr", it));
784: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
785: PetscCall(VecView(vec, viewer));
786: PetscCall(PetscViewerDestroy(&viewer));
787: PetscCall(VecDestroy(&vec));
788: PetscCall(DMDestroy(&da));
789: }
791: {
792: const DMStagStencilLocation loc = DMSTAG_DOWN;
793: const PetscInt c = 0;
794: Vec vec;
795: DM da;
796: PetscViewer viewer;
797: char filename[PETSC_MAX_PATH_LEN];
799: PetscCall(DMStagVecSplitToDMDA(dm, e, loc, c, &da, &vec));
800: PetscCall(PetscSNPrintf(filename, sizeof(filename), "err_num_vy_%" PetscInt_FMT ".vtr", it));
801: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
802: PetscCall(VecView(vec, viewer));
803: PetscCall(PetscViewerDestroy(&viewer));
804: PetscCall(VecDestroy(&vec));
805: PetscCall(DMDestroy(&da));
806: }
808: PetscCall(VecDestroy(&e));
809: }
811: /* Step */
812: if (!ctx->solPrev) {
813: PetscCall(VecDuplicate(sol, &ctx->solPrev));
814: } else {
815: Vec diff;
816: PetscCall(VecDuplicate(sol, &diff));
817: PetscCall(VecCopy(sol, diff));
818: PetscCall(VecAXPY(diff, -1.0, ctx->solPrev));
819: {
820: const DMStagStencilLocation loc = DMSTAG_LEFT;
821: const PetscInt c = 0;
822: Vec vec;
823: DM da;
824: PetscViewer viewer;
825: char filename[PETSC_MAX_PATH_LEN];
827: PetscCall(DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec));
828: PetscCall(PetscSNPrintf(filename, sizeof(filename), "diff_vx_%" PetscInt_FMT ".vtr", it));
829: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
830: PetscCall(VecView(vec, viewer));
831: PetscCall(PetscViewerDestroy(&viewer));
832: PetscCall(VecDestroy(&vec));
833: PetscCall(DMDestroy(&da));
834: }
836: {
837: const DMStagStencilLocation loc = DMSTAG_DOWN;
838: const PetscInt c = 0;
839: Vec vec;
840: DM da;
841: PetscViewer viewer;
842: char filename[PETSC_MAX_PATH_LEN];
844: PetscCall(DMStagVecSplitToDMDA(dm, diff, loc, c, &da, &vec));
845: PetscCall(PetscSNPrintf(filename, sizeof(filename), "diff_vy_%" PetscInt_FMT ".vtr", it));
846: PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)r), filename, FILE_MODE_WRITE, &viewer));
847: PetscCall(VecView(vec, viewer));
848: PetscCall(PetscViewerDestroy(&viewer));
849: PetscCall(VecDestroy(&vec));
850: PetscCall(DMDestroy(&da));
851: }
852: PetscCall(VecDestroy(&diff));
853: }
854: PetscCall(VecCopy(sol, ctx->solPrev));
856: PetscCall(VecDestroy(&r));
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: /* Use a direct solver to create an "exact" solution to the discrete system
861: useful for testing solvers (in that it doesn't include discretization error) */
862: static PetscErrorCode CreateNumericalReferenceSolution(Mat A, Vec rhs, Vec *px)
863: {
864: KSP ksp;
865: PC pc;
866: Vec x;
868: PetscFunctionBeginUser;
869: PetscCall(KSPCreate(PetscObjectComm((PetscObject)A), &ksp));
870: PetscCall(KSPSetType(ksp, KSPPREONLY));
871: PetscCall(KSPGetPC(ksp, &pc));
872: PetscCall(PCSetType(pc, PCLU));
873: PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERUMFPACK));
874: PetscCall(KSPSetOptionsPrefix(ksp, "numref_"));
875: PetscCall(KSPSetFromOptions(ksp));
876: PetscCall(KSPSetOperators(ksp, A, A));
877: PetscCall(VecDuplicate(rhs, px));
878: x = *px;
879: PetscCall(KSPSolve(ksp, rhs, x));
880: PetscCall(KSPDestroy(&ksp));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: /*TEST
886: test:
887: suffix: gmg_1
888: nsize: 1
889: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason
891: test:
892: suffix: gmg_1_b
893: nsize: 1
894: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false
896: test:
897: suffix: gmg_1_c
898: nsize: 1
899: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false
901: test:
902: suffix: gmg_1_d
903: nsize: 1
904: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false
906: test:
907: suffix: gmg_1_bigger
908: requires: !complex
909: nsize: 1
910: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32
911: filter: sed -e "s/ iterations 8/ iterations 7/g"
913: test:
914: suffix: gmg_8
915: requires: !single !complex
916: nsize: 8
917: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason
919: test:
920: suffix: gmg_8_b
921: nsize: 8
922: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractTransferOperators false
924: test:
925: suffix: gmg_8_c
926: nsize: 8
927: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false
929: test:
930: suffix: gmg_8_d
931: nsize: 8
932: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false
934: test:
935: suffix: gmg_8_galerkin
936: nsize: 8
937: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -extractSystem false -extractTransferOperators false -pc_mg_galerkin
939: test:
940: suffix: gmg_8_bigger
941: requires: !complex
942: nsize: 8
943: args: -ksp_type fgmres -mg_levels_pc_type jacobi -ksp_converged_reason -stag_grid_x 32 -stag_grid_y 32
945: TEST*/