Actual source code: ex2.c
1: static char help[] = "Solve a toy 2D problem on a staggered grid\n\n";
2: /*
4: To demonstrate the basic functionality of DMStag, solves an isoviscous
5: incompressible Stokes problem on a rectangular 2D domain, using a manufactured
6: solution.
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: Use the -pinpressure option to fix a pressure node, instead of providing
18: a constant-pressure nullspace. This allows use of direct solvers, e.g. to
19: use UMFPACK,
21: ./ex2 -pinpressure 1 -pc_type lu -pc_factor_mat_solver_type umfpack
23: This example demonstrates the use of DMProduct to efficiently store coordinates
24: on an orthogonal grid.
26: */
27: #include <petscdm.h>
28: #include <petscksp.h>
29: #include <petscdmstag.h>
31: /* Shorter, more convenient names for DMStagStencilLocation entries */
32: #define DOWN_LEFT DMSTAG_DOWN_LEFT
33: #define DOWN DMSTAG_DOWN
34: #define DOWN_RIGHT DMSTAG_DOWN_RIGHT
35: #define LEFT DMSTAG_LEFT
36: #define ELEMENT DMSTAG_ELEMENT
37: #define RIGHT DMSTAG_RIGHT
38: #define UP_LEFT DMSTAG_UP_LEFT
39: #define UP DMSTAG_UP
40: #define UP_RIGHT DMSTAG_UP_RIGHT
42: static PetscErrorCode CreateSystem(DM, Mat *, Vec *, PetscBool);
43: static PetscErrorCode CreateReferenceSolution(DM, Vec *);
44: static PetscErrorCode AttachNullspace(DM, Mat);
45: static PetscErrorCode CheckSolution(Vec, Vec);
47: /* Manufactured solution. Chosen to be higher order than can be solved exactly,
48: and to have a zero derivative for flow parallel to the boundaries. That is,
49: d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
50: and left boundaries. */
51: static PetscScalar uxRef(PetscScalar x, PetscScalar y)
52: {
53: return 0.0 * x + y * y - 2.0 * y * y * y + y * y * y * y;
54: } /* no x-dependence */
55: static PetscScalar uyRef(PetscScalar x, PetscScalar y)
56: {
57: return x * x - 2.0 * x * x * x + x * x * x * x + 0.0 * y;
58: } /* no y-dependence */
59: static PetscScalar pRef(PetscScalar x, PetscScalar y)
60: {
61: return -1.0 * (x - 0.5) + -3.0 / 2.0 * y * y + 0.5;
62: } /* zero integral */
63: static PetscScalar fx(PetscScalar x, PetscScalar y)
64: {
65: return 0.0 * x + 2.0 - 12.0 * y + 12.0 * y * y + 1.0;
66: } /* no x-dependence */
67: static PetscScalar fy(PetscScalar x, PetscScalar y)
68: {
69: return 2.0 - 12.0 * x + 12.0 * x * x + 3.0 * y;
70: }
71: static PetscScalar g(PetscScalar x, PetscScalar y)
72: {
73: return 0.0 * x * y;
74: } /* identically zero */
76: int main(int argc, char **argv)
77: {
78: DM dmSol;
79: Vec sol, solRef, rhs;
80: Mat A;
81: KSP ksp;
82: PC pc;
83: PetscBool pinPressure;
85: /* Initialize PETSc and process command line arguments */
86: PetscFunctionBeginUser;
87: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
88: pinPressure = PETSC_TRUE;
89: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pinpressure", &pinPressure, NULL));
91: /* Create 2D DMStag for the solution, and set up. */
92: {
93: const PetscInt dof0 = 0, dof1 = 1, dof2 = 1; /* 1 dof on each edge and element center */
94: const PetscInt stencilWidth = 1;
95: PetscCall(DMStagCreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 7, 9, PETSC_DECIDE, PETSC_DECIDE, dof0, dof1, dof2, DMSTAG_STENCIL_BOX, stencilWidth, NULL, NULL, &dmSol));
96: PetscCall(DMSetFromOptions(dmSol));
97: PetscCall(DMSetUp(dmSol));
98: }
100: /* Define uniform coordinates as a product of 1D arrays */
101: PetscCall(DMStagSetUniformCoordinatesProduct(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
103: /* Compute (manufactured) reference solution */
104: PetscCall(CreateReferenceSolution(dmSol, &solRef));
106: /* Assemble system */
107: PetscCall(CreateSystem(dmSol, &A, &rhs, pinPressure));
109: /* Attach a constant-pressure nullspace to the operator
110: (logically, this should be in CreateSystem, but we separate it here to highlight
111: the features of DMStag exposed, in particular DMStagMigrateVec()) */
112: if (!pinPressure) PetscCall(AttachNullspace(dmSol, A));
114: /* Solve, using the default FieldSplit (Approximate Block Factorization) Preconditioner
115: This is not intended to be an example of a good solver! */
116: PetscCall(DMCreateGlobalVector(dmSol, &sol));
117: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
118: PetscCall(KSPSetType(ksp, KSPFGMRES));
119: PetscCall(KSPSetOperators(ksp, A, A));
120: PetscCall(KSPGetPC(ksp, &pc));
121: PetscCall(PCSetType(pc, PCFIELDSPLIT));
122: PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, PETSC_TRUE));
123: PetscCall(KSPSetFromOptions(ksp));
124: PetscCall(KSPSolve(ksp, rhs, sol));
126: /* Check Solution */
127: PetscCall(CheckSolution(sol, solRef));
129: /* Clean up and finalize PETSc */
130: PetscCall(KSPDestroy(&ksp));
131: PetscCall(VecDestroy(&sol));
132: PetscCall(VecDestroy(&solRef));
133: PetscCall(VecDestroy(&rhs));
134: PetscCall(MatDestroy(&A));
135: PetscCall(DMDestroy(&dmSol));
136: PetscCall(PetscFinalize());
137: return 0;
138: }
140: /*
141: Note: this system is not well-scaled! Generally one would adjust the equations
142: to try to get matrix entries to be of comparable order, regardless of grid spacing
143: or choice of coefficients.
144: */
145: static PetscErrorCode CreateSystem(DM dmSol, Mat *pA, Vec *pRhs, PetscBool pinPressure)
146: {
147: PetscInt N[2];
148: PetscInt ex, ey, startx, starty, nx, ny;
149: PetscInt iprev, icenter, inext;
150: Mat A;
151: Vec rhs;
152: PetscReal hx, hy;
153: PetscScalar **cArrX, **cArrY;
155: /* Here, we showcase two different methods for manipulating local vector entries.
156: One can use DMStagStencil objects with DMStagVecSetValuesStencil(),
157: making sure to call VecAssemble[Begin/End]() after all values are set.
158: Alternately, one can use DMStagVecGetArray[Read]() and DMStagVecRestoreArray[Read]().
159: The first approach is used to build the rhs, and the second is used to
160: obtain coordinate values. Working with the array is almost certainly more efficient,
161: but only allows setting local entries, requires understanding which "slot" to use,
162: and doesn't correspond as precisely to the matrix assembly process using DMStagStencil objects */
164: PetscFunctionBeginUser;
165: PetscCall(DMCreateMatrix(dmSol, pA));
166: A = *pA;
167: PetscCall(DMCreateGlobalVector(dmSol, pRhs));
168: rhs = *pRhs;
169: PetscCall(DMStagGetCorners(dmSol, &startx, &starty, NULL, &nx, &ny, NULL, NULL, NULL, NULL));
170: PetscCall(DMStagGetGlobalSizes(dmSol, &N[0], &N[1], NULL));
171: hx = 1.0 / N[0];
172: hy = 1.0 / N[1];
173: PetscCall(DMStagGetProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
174: PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, ELEMENT, &icenter));
175: PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, LEFT, &iprev));
176: PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, RIGHT, &inext));
178: /* Loop over all local elements. Note that it may be more efficient in real
179: applications to loop over each boundary separately */
180: for (ey = starty; ey < starty + ny; ++ey) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
181: for (ex = startx; ex < startx + nx; ++ex) {
182: if (ex == N[0] - 1) {
183: /* Right Boundary velocity Dirichlet */
184: DMStagStencil row;
185: PetscScalar valRhs;
186: const PetscScalar valA = 1.0;
187: row.i = ex;
188: row.j = ey;
189: row.loc = RIGHT;
190: row.c = 0;
191: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
192: valRhs = uxRef(cArrX[ex][inext], cArrY[ey][icenter]);
193: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
194: }
195: if (ey == N[1] - 1) {
196: /* Top boundary velocity Dirichlet */
197: DMStagStencil row;
198: PetscScalar valRhs;
199: const PetscScalar valA = 1.0;
200: row.i = ex;
201: row.j = ey;
202: row.loc = UP;
203: row.c = 0;
204: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
205: valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][inext]);
206: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
207: }
209: if (ey == 0) {
210: /* Bottom boundary velocity Dirichlet */
211: DMStagStencil row;
212: PetscScalar valRhs;
213: const PetscScalar valA = 1.0;
214: row.i = ex;
215: row.j = ey;
216: row.loc = DOWN;
217: row.c = 0;
218: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
219: valRhs = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
220: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
221: } else {
222: /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
223: DMStagStencil row, col[7];
224: PetscScalar valA[7], valRhs;
225: PetscInt nEntries;
227: row.i = ex;
228: row.j = ey;
229: row.loc = DOWN;
230: row.c = 0;
231: if (ex == 0) {
232: nEntries = 6;
233: col[0].i = ex;
234: col[0].j = ey;
235: col[0].loc = DOWN;
236: col[0].c = 0;
237: valA[0] = -1.0 / (hx * hx) - 2.0 / (hy * hy);
238: col[1].i = ex;
239: col[1].j = ey - 1;
240: col[1].loc = DOWN;
241: col[1].c = 0;
242: valA[1] = 1.0 / (hy * hy);
243: col[2].i = ex;
244: col[2].j = ey + 1;
245: col[2].loc = DOWN;
246: col[2].c = 0;
247: valA[2] = 1.0 / (hy * hy);
248: /* Missing left element */
249: col[3].i = ex + 1;
250: col[3].j = ey;
251: col[3].loc = DOWN;
252: col[3].c = 0;
253: valA[3] = 1.0 / (hx * hx);
254: col[4].i = ex;
255: col[4].j = ey - 1;
256: col[4].loc = ELEMENT;
257: col[4].c = 0;
258: valA[4] = 1.0 / hy;
259: col[5].i = ex;
260: col[5].j = ey;
261: col[5].loc = ELEMENT;
262: col[5].c = 0;
263: valA[5] = -1.0 / hy;
264: } else if (ex == N[0] - 1) {
265: /* Right boundary y velocity stencil */
266: nEntries = 6;
267: col[0].i = ex;
268: col[0].j = ey;
269: col[0].loc = DOWN;
270: col[0].c = 0;
271: valA[0] = -1.0 / (hx * hx) - 2.0 / (hy * hy);
272: col[1].i = ex;
273: col[1].j = ey - 1;
274: col[1].loc = DOWN;
275: col[1].c = 0;
276: valA[1] = 1.0 / (hy * hy);
277: col[2].i = ex;
278: col[2].j = ey + 1;
279: col[2].loc = DOWN;
280: col[2].c = 0;
281: valA[2] = 1.0 / (hy * hy);
282: col[3].i = ex - 1;
283: col[3].j = ey;
284: col[3].loc = DOWN;
285: col[3].c = 0;
286: valA[3] = 1.0 / (hx * hx);
287: /* Missing right element */
288: col[4].i = ex;
289: col[4].j = ey - 1;
290: col[4].loc = ELEMENT;
291: col[4].c = 0;
292: valA[4] = 1.0 / hy;
293: col[5].i = ex;
294: col[5].j = ey;
295: col[5].loc = ELEMENT;
296: col[5].c = 0;
297: valA[5] = -1.0 / hy;
298: } else {
299: nEntries = 7;
300: col[0].i = ex;
301: col[0].j = ey;
302: col[0].loc = DOWN;
303: col[0].c = 0;
304: valA[0] = -2.0 / (hx * hx) - 2.0 / (hy * hy);
305: col[1].i = ex;
306: col[1].j = ey - 1;
307: col[1].loc = DOWN;
308: col[1].c = 0;
309: valA[1] = 1.0 / (hy * hy);
310: col[2].i = ex;
311: col[2].j = ey + 1;
312: col[2].loc = DOWN;
313: col[2].c = 0;
314: valA[2] = 1.0 / (hy * hy);
315: col[3].i = ex - 1;
316: col[3].j = ey;
317: col[3].loc = DOWN;
318: col[3].c = 0;
319: valA[3] = 1.0 / (hx * hx);
320: col[4].i = ex + 1;
321: col[4].j = ey;
322: col[4].loc = DOWN;
323: col[4].c = 0;
324: valA[4] = 1.0 / (hx * hx);
325: col[5].i = ex;
326: col[5].j = ey - 1;
327: col[5].loc = ELEMENT;
328: col[5].c = 0;
329: valA[5] = 1.0 / hy;
330: col[6].i = ex;
331: col[6].j = ey;
332: col[6].loc = ELEMENT;
333: col[6].c = 0;
334: valA[6] = -1.0 / hy;
335: }
336: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
337: valRhs = fy(cArrX[ex][icenter], cArrY[ey][iprev]);
338: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
339: }
341: if (ex == 0) {
342: /* Left velocity Dirichlet */
343: DMStagStencil row;
344: PetscScalar valRhs;
345: const PetscScalar valA = 1.0;
346: row.i = ex;
347: row.j = ey;
348: row.loc = LEFT;
349: row.c = 0;
350: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
351: valRhs = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
352: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
353: } else {
354: /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
355: DMStagStencil row, col[7];
356: PetscScalar valA[7], valRhs;
357: PetscInt nEntries;
358: row.i = ex;
359: row.j = ey;
360: row.loc = LEFT;
361: row.c = 0;
363: if (ey == 0) {
364: nEntries = 6;
365: col[0].i = ex;
366: col[0].j = ey;
367: col[0].loc = LEFT;
368: col[0].c = 0;
369: valA[0] = -2.0 / (hx * hx) - 1.0 / (hy * hy);
370: /* missing term from element below */
371: col[1].i = ex;
372: col[1].j = ey + 1;
373: col[1].loc = LEFT;
374: col[1].c = 0;
375: valA[1] = 1.0 / (hy * hy);
376: col[2].i = ex - 1;
377: col[2].j = ey;
378: col[2].loc = LEFT;
379: col[2].c = 0;
380: valA[2] = 1.0 / (hx * hx);
381: col[3].i = ex + 1;
382: col[3].j = ey;
383: col[3].loc = LEFT;
384: col[3].c = 0;
385: valA[3] = 1.0 / (hx * hx);
386: col[4].i = ex - 1;
387: col[4].j = ey;
388: col[4].loc = ELEMENT;
389: col[4].c = 0;
390: valA[4] = 1.0 / hx;
391: col[5].i = ex;
392: col[5].j = ey;
393: col[5].loc = ELEMENT;
394: col[5].c = 0;
395: valA[5] = -1.0 / hx;
396: } else if (ey == N[1] - 1) {
397: /* Top boundary x velocity stencil */
398: nEntries = 6;
399: row.i = ex;
400: row.j = ey;
401: row.loc = LEFT;
402: row.c = 0;
403: col[0].i = ex;
404: col[0].j = ey;
405: col[0].loc = LEFT;
406: col[0].c = 0;
407: valA[0] = -2.0 / (hx * hx) - 1.0 / (hy * hy);
408: col[1].i = ex;
409: col[1].j = ey - 1;
410: col[1].loc = LEFT;
411: col[1].c = 0;
412: valA[1] = 1.0 / (hy * hy);
413: /* Missing element above term */
414: col[2].i = ex - 1;
415: col[2].j = ey;
416: col[2].loc = LEFT;
417: col[2].c = 0;
418: valA[2] = 1.0 / (hx * hx);
419: col[3].i = ex + 1;
420: col[3].j = ey;
421: col[3].loc = LEFT;
422: col[3].c = 0;
423: valA[3] = 1.0 / (hx * hx);
424: col[4].i = ex - 1;
425: col[4].j = ey;
426: col[4].loc = ELEMENT;
427: col[4].c = 0;
428: valA[4] = 1.0 / hx;
429: col[5].i = ex;
430: col[5].j = ey;
431: col[5].loc = ELEMENT;
432: col[5].c = 0;
433: valA[5] = -1.0 / hx;
434: } else {
435: /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
436: nEntries = 7;
437: col[0].i = ex;
438: col[0].j = ey;
439: col[0].loc = LEFT;
440: col[0].c = 0;
441: valA[0] = -2.0 / (hx * hx) + -2.0 / (hy * hy);
442: col[1].i = ex;
443: col[1].j = ey - 1;
444: col[1].loc = LEFT;
445: col[1].c = 0;
446: valA[1] = 1.0 / (hy * hy);
447: col[2].i = ex;
448: col[2].j = ey + 1;
449: col[2].loc = LEFT;
450: col[2].c = 0;
451: valA[2] = 1.0 / (hy * hy);
452: col[3].i = ex - 1;
453: col[3].j = ey;
454: col[3].loc = LEFT;
455: col[3].c = 0;
456: valA[3] = 1.0 / (hx * hx);
457: col[4].i = ex + 1;
458: col[4].j = ey;
459: col[4].loc = LEFT;
460: col[4].c = 0;
461: valA[4] = 1.0 / (hx * hx);
462: col[5].i = ex - 1;
463: col[5].j = ey;
464: col[5].loc = ELEMENT;
465: col[5].c = 0;
466: valA[5] = 1.0 / hx;
467: col[6].i = ex;
468: col[6].j = ey;
469: col[6].loc = ELEMENT;
470: col[6].c = 0;
471: valA[6] = -1.0 / hx;
472: }
473: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES));
474: valRhs = fx(cArrX[ex][iprev], cArrY[ey][icenter]);
475: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
476: }
478: /* P equation : u_x + v_y = g
479: Note that this includes an explicit zero on the diagonal. This is only needed for
480: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
481: if (pinPressure && ex == 0 && ey == 0) { /* Pin the first pressure node, if requested */
482: DMStagStencil row;
483: PetscScalar valA, valRhs;
484: row.i = ex;
485: row.j = ey;
486: row.loc = ELEMENT;
487: row.c = 0;
488: valA = 1.0;
489: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES));
490: valRhs = pRef(cArrX[ex][icenter], cArrY[ey][icenter]);
491: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
492: } else {
493: DMStagStencil row, col[5];
494: PetscScalar valA[5], valRhs;
496: row.i = ex;
497: row.j = ey;
498: row.loc = ELEMENT;
499: row.c = 0;
500: col[0].i = ex;
501: col[0].j = ey;
502: col[0].loc = LEFT;
503: col[0].c = 0;
504: valA[0] = -1.0 / hx;
505: col[1].i = ex;
506: col[1].j = ey;
507: col[1].loc = RIGHT;
508: col[1].c = 0;
509: valA[1] = 1.0 / hx;
510: col[2].i = ex;
511: col[2].j = ey;
512: col[2].loc = DOWN;
513: col[2].c = 0;
514: valA[2] = -1.0 / hy;
515: col[3].i = ex;
516: col[3].j = ey;
517: col[3].loc = UP;
518: col[3].c = 0;
519: valA[3] = 1.0 / hy;
520: col[4] = row;
521: valA[4] = 0.0;
522: PetscCall(DMStagMatSetValuesStencil(dmSol, A, 1, &row, 5, col, valA, INSERT_VALUES));
523: valRhs = g(cArrX[ex][icenter], cArrY[ey][icenter]);
524: PetscCall(DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES));
525: }
526: }
527: }
528: PetscCall(DMStagRestoreProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
529: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
530: PetscCall(VecAssemblyBegin(rhs));
531: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
532: PetscCall(VecAssemblyEnd(rhs));
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: /* Create a pressure-only DMStag and use it to generate a nullspace vector
537: - Create a compatible DMStag with one dof per element (and nothing else).
538: - Create a constant vector and normalize it
539: - Migrate it to a vector on the original dmSol (making use of the fact
540: that this will fill in zeros for "extra" dof)
541: - Set the nullspace for the operator
542: - Destroy everything (the operator keeps the references it needs) */
543: static PetscErrorCode AttachNullspace(DM dmSol, Mat A)
544: {
545: DM dmPressure;
546: Vec constantPressure, basis;
547: PetscReal nrm;
548: MatNullSpace matNullSpace;
550: PetscFunctionBeginUser;
551: PetscCall(DMStagCreateCompatibleDMStag(dmSol, 0, 0, 1, 0, &dmPressure));
552: PetscCall(DMGetGlobalVector(dmPressure, &constantPressure));
553: PetscCall(VecSet(constantPressure, 1.0));
554: PetscCall(VecNorm(constantPressure, NORM_2, &nrm));
555: PetscCall(VecScale(constantPressure, 1.0 / nrm));
556: PetscCall(DMCreateGlobalVector(dmSol, &basis));
557: PetscCall(DMStagMigrateVec(dmPressure, constantPressure, dmSol, basis));
558: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol), PETSC_FALSE, 1, &basis, &matNullSpace));
559: PetscCall(VecDestroy(&basis));
560: PetscCall(VecDestroy(&constantPressure));
561: PetscCall(MatSetNullSpace(A, matNullSpace));
562: PetscCall(MatNullSpaceDestroy(&matNullSpace));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /* Create a reference solution.
567: Here, we use the more direct method of iterating over arrays. */
568: static PetscErrorCode CreateReferenceSolution(DM dmSol, Vec *pSolRef)
569: {
570: PetscInt startx, starty, nx, ny, nExtra[2], ex, ey;
571: PetscInt iuy, iux, ip, iprev, icenter;
572: PetscScalar ***arrSol, **cArrX, **cArrY;
573: Vec solRefLocal;
575: PetscFunctionBeginUser;
576: PetscCall(DMCreateGlobalVector(dmSol, pSolRef));
577: PetscCall(DMGetLocalVector(dmSol, &solRefLocal));
579: /* Obtain indices to use in the raw arrays */
580: PetscCall(DMStagGetLocationSlot(dmSol, DOWN, 0, &iuy));
581: PetscCall(DMStagGetLocationSlot(dmSol, LEFT, 0, &iux));
582: PetscCall(DMStagGetLocationSlot(dmSol, ELEMENT, 0, &ip));
584: /* Use high-level convenience functions to get raw arrays and indices for 1d coordinates */
585: PetscCall(DMStagGetProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
586: PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, ELEMENT, &icenter));
587: PetscCall(DMStagGetProductCoordinateLocationSlot(dmSol, LEFT, &iprev));
589: PetscCall(DMStagVecGetArray(dmSol, solRefLocal, &arrSol));
590: PetscCall(DMStagGetCorners(dmSol, &startx, &starty, NULL, &nx, &ny, NULL, &nExtra[0], &nExtra[1], NULL));
591: for (ey = starty; ey < starty + ny + nExtra[1]; ++ey) {
592: for (ex = startx; ex < startx + nx + nExtra[0]; ++ex) {
593: arrSol[ey][ex][iuy] = uyRef(cArrX[ex][icenter], cArrY[ey][iprev]);
594: arrSol[ey][ex][iux] = uxRef(cArrX[ex][iprev], cArrY[ey][icenter]);
595: if (ey < starty + ny && ex < startx + nx) { /* Don't fill on the dummy elements (though you could, and these values would just be ignored) */
596: arrSol[ey][ex][ip] = pRef(cArrX[ex][icenter], cArrY[ey][icenter]);
597: }
598: }
599: }
600: PetscCall(DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol));
601: PetscCall(DMStagRestoreProductCoordinateArraysRead(dmSol, &cArrX, &cArrY, NULL));
602: PetscCall(DMLocalToGlobal(dmSol, solRefLocal, INSERT_VALUES, *pSolRef));
603: PetscCall(DMRestoreLocalVector(dmSol, &solRefLocal));
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: static PetscErrorCode CheckSolution(Vec sol, Vec solRef)
608: {
609: Vec diff;
610: PetscReal normsolRef, errAbs, errRel;
612: PetscFunctionBeginUser;
613: PetscCall(VecDuplicate(sol, &diff));
614: PetscCall(VecCopy(sol, diff));
615: PetscCall(VecAXPY(diff, -1.0, solRef));
616: PetscCall(VecNorm(diff, NORM_2, &errAbs));
617: PetscCall(VecNorm(solRef, NORM_2, &normsolRef));
618: errRel = errAbs / normsolRef;
619: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel));
620: PetscCall(VecDestroy(&diff));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: /*TEST
626: test:
627: suffix: 1
628: nsize: 4
629: args: -ksp_monitor_short -ksp_converged_reason
631: test:
632: suffix: direct_umfpack
633: requires: suitesparse
634: nsize: 1
635: args: -pinpressure 1 -stag_grid_x 8 -stag_grid_y 6 -ksp_monitor_short -pc_type lu -pc_factor_mat_solver_type umfpack
637: TEST*/