Actual source code: ex2.c
petsc-3.11.4 2019-09-28
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 ){return 0.0*x + y*y - 2.0*y*y*y + y*y*y*y;} /* no x-dependence */
52: static PetscScalar uyRef(PetscScalar x,PetscScalar y) {return x*x - 2.0*x*x*x + x*x*x*x + 0.0*y;} /* no y-dependence */
53: static PetscScalar pRef (PetscScalar x,PetscScalar y) {return -1.0*(x-0.5) + -3.0/2.0*y*y + 0.5;} /* zero integral */
54: static PetscScalar fx (PetscScalar x,PetscScalar y) {return 0.0*x + 2.0 -12.0*y + 12.0*y*y + 1.0;} /* no x-dependence */
55: static PetscScalar fy (PetscScalar x,PetscScalar y) {return 2.0 -12.0*x + 12.0*x*x + 3.0*y;}
56: static PetscScalar g (PetscScalar x,PetscScalar y) {return 0.0*x*y;} /* identically zero */
58: int main(int argc,char **argv)
59: {
61: DM dmSol;
62: Vec sol,solRef,rhs;
63: Mat A;
64: KSP ksp;
65: PC pc;
66: PetscBool pinPressure;
68: /* Initialize PETSc and process command line arguments */
69: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
70: pinPressure = PETSC_TRUE;
71: PetscOptionsGetBool(NULL,NULL,"-pinpressure",&pinPressure,NULL);
73: /* Create 2D DMStag for the solution, and set up. */
74: {
75: const PetscInt dof0 = 0, dof1 = 1,dof2 = 1; /* 1 dof on each edge and element center */
76: const PetscInt stencilWidth = 1;
77: 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);
78: DMSetFromOptions(dmSol);
79: DMSetUp(dmSol);
80: }
82: /* Define uniform coordinates as a product of 1D arrays */
83: DMStagSetUniformCoordinatesProduct(dmSol,0.0,1.0,0.0,1.0,0.0,0.0);
85: /* Compute (manufactured) reference solution */
86: CreateReferenceSolution(dmSol,&solRef);
88: /* Assemble system */
89: CreateSystem(dmSol,&A,&rhs,pinPressure);
91: /* Attach a constant-pressure nullspace to the operator
92: (logically, this should be in CreateSystem, but we separate it here to highlight
93: the features of DMStag exposed, in particular DMStagMigrateVec() ) */
94: if (!pinPressure) {
95: AttachNullspace(dmSol,A);
96: }
98: /* Solve, using the default FieldSplit (Approximate Block Factorization) Preconditioner
99: This is not intended to be an example of a good solver! */
100: DMCreateGlobalVector(dmSol,&sol);
101: KSPCreate(PETSC_COMM_WORLD,&ksp);
102: KSPSetType(ksp,KSPFGMRES);
103: KSPSetOperators(ksp,A,A);
104: KSPGetPC(ksp,&pc);
105: PCSetType(pc,PCFIELDSPLIT);
106: PCFieldSplitSetDetectSaddlePoint(pc,PETSC_TRUE);
107: KSPSetFromOptions(ksp);
108: KSPSolve(ksp,rhs,sol);
110: /* Check Solution */
111: CheckSolution(sol,solRef);
113: /* Clean up and finalize PETSc */
114: KSPDestroy(&ksp);
115: VecDestroy(&sol);
116: VecDestroy(&solRef);
117: VecDestroy(&rhs);
118: MatDestroy(&A);
119: DMDestroy(&dmSol);
120: PetscFinalize();
121: return ierr;
122: }
124: /*
125: Note: this system is not well-scaled! Generally one would adjust the equations
126: to try to get matrix entries to be of comparable order, regardless of grid spacing
127: or choice of coefficients.
128: */
129: static PetscErrorCode CreateSystem(DM dmSol,Mat *pA,Vec *pRhs, PetscBool pinPressure)
130: {
132: PetscInt N[2];
133: PetscBool isLastRankx,isLastRanky,isFirstRankx,isFirstRanky;
134: PetscInt ex,ey,startx,starty,nx,ny;
135: PetscInt iprev,icenter,inext;
136: Mat A;
137: Vec rhs;
138: PetscReal hx,hy;
139: PetscScalar **cArrX,**cArrY;
141: /* Here, we showcase two different methods for manipulating local vector entries.
142: One can use DMStagStencil objects with DMStagVecSetValuesStencil(),
143: making sure to call VecAssemble[Begin/End]() after all values are set.
144: Alternately, one can use DMStagVecGetArrayDOF[Read]() and DMStagVecRestoreArrayDOF[Read]().
145: The first approach is used to build the rhs, and the second is used to
146: obtain coordinate values. Working with the array is almost certainly more efficient,
147: but only allows setting local entries, requires understanding which "slot" to use,
148: and doesn't correspond as precisely to the matrix assembly process using DMStagStencil objects */
151: DMCreateMatrix(dmSol,pA);
152: A = *pA;
153: DMCreateGlobalVector(dmSol,pRhs);
154: rhs = *pRhs;
155: DMStagGetCorners(dmSol,&startx,&starty,NULL,&nx,&ny,NULL,NULL,NULL,NULL);
156: DMStagGetGlobalSizes(dmSol,&N[0],&N[1],NULL);
157: DMStagGetIsLastRank(dmSol,&isLastRankx,&isLastRanky,NULL);
158: DMStagGetIsFirstRank(dmSol,&isFirstRankx,&isFirstRanky,NULL);
159: hx = 1.0/N[0]; hy = 1.0/N[1];
160: DMStagGet1dCoordinateArraysDOFRead(dmSol,&cArrX,&cArrY,NULL);
161: DMStagGet1dCoordinateLocationSlot(dmSol,ELEMENT,&icenter);
162: DMStagGet1dCoordinateLocationSlot(dmSol,LEFT,&iprev);
163: DMStagGet1dCoordinateLocationSlot(dmSol,RIGHT,&inext);
165: /* Loop over all local elements. Note that it may be more efficient in real
166: applications to loop over each boundary separately */
167: for (ey = starty; ey<starty+ny; ++ey) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
168: for (ex = startx; ex<startx+nx; ++ex) {
170: if (ex == N[0]-1) {
171: /* Right Boundary velocity Dirichlet */
172: DMStagStencil row;
173: PetscScalar valRhs;
174: const PetscScalar valA = 1.0;
175: row.i = ex; row.j = ey; row.loc = RIGHT; row.c = 0;
176: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
177: valRhs = uxRef(cArrX[ex][inext],cArrY[ey][icenter]);
178: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
179: }
180: if (ey == N[1]-1) {
181: /* Top boundary velocity Dirichlet */
182: DMStagStencil row;
183: PetscScalar valRhs;
184: const PetscScalar valA = 1.0;
185: row.i = ex; row.j = ey; row.loc = UP; row.c = 0;
186: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
187: valRhs = uyRef(cArrX[ex][icenter],cArrY[ey][inext]);
188: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
189: }
191: if (ey == 0) {
192: /* Bottom boundary velocity Dirichlet */
193: DMStagStencil row;
194: PetscScalar valRhs;
195: const PetscScalar valA = 1.0;
196: row.i = ex; row.j = ey; row.loc = DOWN; row.c = 0;
197: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
198: valRhs = uyRef(cArrX[ex][icenter],cArrY[ey][iprev]);
199: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
200: } else {
201: /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
202: DMStagStencil row,col[7];
203: PetscScalar valA[7],valRhs;
204: PetscInt nEntries;
206: row.i = ex ; row.j = ey ; row.loc = DOWN; row.c = 0;
207: if (ex == 0) {
208: nEntries = 6;
209: col[0].i = ex ; col[0].j = ey ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) -2.0 / (hy*hy);
210: col[1].i = ex ; col[1].j = ey-1; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
211: col[2].i = ex ; col[2].j = ey+1; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
212: /* Missing left element */
213: col[3].i = ex+1; col[3].j = ey ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
214: col[4].i = ex ; col[4].j = ey-1; col[4].loc = ELEMENT; col[4].c = 0; valA[4] = 1.0 / hy;
215: col[5].i = ex ; col[5].j = ey ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = -1.0 / hy;
216: } else if (ex == N[0]-1) {
217: /* Right boundary y velocity stencil */
218: nEntries = 6;
219: col[0].i = ex ; col[0].j = ey ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) -2.0 / (hy*hy);
220: col[1].i = ex ; col[1].j = ey-1; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
221: col[2].i = ex ; col[2].j = ey+1; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
222: col[3].i = ex-1; col[3].j = ey ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
223: /* Missing right element */
224: col[4].i = ex ; col[4].j = ey-1; col[4].loc = ELEMENT; col[4].c = 0; valA[4] = 1.0 / hy;
225: col[5].i = ex ; col[5].j = ey ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = -1.0 / hy;
226: } else {
227: nEntries = 7;
228: col[0].i = ex ; col[0].j = ey ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -2.0 / (hx*hx) -2.0 / (hy*hy);
229: col[1].i = ex ; col[1].j = ey-1; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
230: col[2].i = ex ; col[2].j = ey+1; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
231: col[3].i = ex-1; col[3].j = ey ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
232: col[4].i = ex+1; col[4].j = ey ; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
233: col[5].i = ex ; col[5].j = ey-1; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hy;
234: col[6].i = ex ; col[6].j = ey ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hy;
235: }
236: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
237: valRhs = fy(cArrX[ex][icenter],cArrY[ey][iprev]);
238: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
239: }
241: if (ex == 0) {
242: /* Left velocity Dirichlet */
243: DMStagStencil row;
244: PetscScalar valRhs;
245: const PetscScalar valA = 1.0;
246: row.i = ex; row.j = ey; row.loc = LEFT; row.c = 0;
247: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
248: valRhs = uxRef(cArrX[ex][iprev],cArrY[ey][icenter]);
249: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
250: } else {
251: /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
252: DMStagStencil row,col[7];
253: PetscScalar valA[7],valRhs;
254: PetscInt nEntries;
255: row.i = ex ; row.j = ey ; row.loc = LEFT; row.c = 0;
257: if (ey == 0) {
258: nEntries = 6;
259: col[0].i = ex ; col[0].j = ey ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 /(hx*hx) -1.0 /(hy*hy);
260: /* missing term from element below */
261: col[1].i = ex ; col[1].j = ey+1; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
262: col[2].i = ex-1; col[2].j = ey ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
263: col[3].i = ex+1; col[3].j = ey ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
264: col[4].i = ex-1; col[4].j = ey ; col[4].loc = ELEMENT; col[4].c = 0; valA[4] = 1.0 / hx;
265: col[5].i = ex ; col[5].j = ey ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = -1.0 / hx;
266: } else if (ey == N[1]-1) {
267: /* Top boundary x velocity stencil */
268: nEntries = 6;
269: row.i = ex ; row.j = ey ; row.loc = LEFT; row.c = 0;
270: col[0].i = ex ; col[0].j = ey ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) -1.0 / (hy*hy);
271: col[1].i = ex ; col[1].j = ey-1; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
272: /* Missing element above term */
273: col[2].i = ex-1; col[2].j = ey ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
274: col[3].i = ex+1; col[3].j = ey ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
275: col[4].i = ex-1; col[4].j = ey ; col[4].loc = ELEMENT; col[4].c = 0; valA[4] = 1.0 / hx;
276: col[5].i = ex ; col[5].j = ey ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = -1.0 / hx;
277: } else {
278: /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
279: nEntries = 7;
280: col[0].i = ex ; col[0].j = ey ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy);
281: col[1].i = ex ; col[1].j = ey-1; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
282: col[2].i = ex ; col[2].j = ey+1; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
283: col[3].i = ex-1; col[3].j = ey ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
284: col[4].i = ex+1; col[4].j = ey ; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
285: col[5].i = ex-1; col[5].j = ey ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hx;
286: col[6].i = ex ; col[6].j = ey ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hx;
288: }
289: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
290: valRhs = fx(cArrX[ex][iprev],cArrY[ey][icenter]);
291: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
292: }
294: /* P equation : u_x + v_y = g
295: Note that this includes an explicit zero on the diagonal. This is only needed for
296: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
297: if (pinPressure && ex == 0 && ey == 0) { /* Pin the first pressure node, if requested */
298: DMStagStencil row;
299: PetscScalar valA,valRhs;
300: row.i = ex; row.j = ey; row.loc = ELEMENT; row.c = 0;
301: valA = 1.0;
302: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
303: valRhs = pRef(cArrX[ex][icenter],cArrY[ey][icenter]);
304: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
305: } else {
306: DMStagStencil row,col[5];
307: PetscScalar valA[5],valRhs;
309: row.i = ex; row.j = ey; row.loc = ELEMENT; row.c = 0;
310: col[0].i = ex; col[0].j = ey; col[0].loc = LEFT; col[0].c = 0; valA[0] = -1.0 / hx;
311: col[1].i = ex; col[1].j = ey; col[1].loc = RIGHT; col[1].c = 0; valA[1] = 1.0 / hx;
312: col[2].i = ex; col[2].j = ey; col[2].loc = DOWN; col[2].c = 0; valA[2] = -1.0 / hy;
313: col[3].i = ex; col[3].j = ey; col[3].loc = UP; col[3].c = 0; valA[3] = 1.0 / hy;
314: col[4] = row; valA[4] = 0.0;
315: DMStagMatSetValuesStencil(dmSol,A,1,&row,5,col,valA,INSERT_VALUES);
316: valRhs = g(cArrX[ex][icenter],cArrY[ey][icenter]);
317: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
318: }
319: }
320: }
321: DMStagRestore1dCoordinateArraysDOFRead(dmSol,&cArrX,&cArrY,NULL);
322: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
323: VecAssemblyBegin(rhs);
324: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
325: VecAssemblyEnd(rhs);
326: return(0);
327: }
329: /* Create a pressure-only DMStag and use it to generate a nullspace vector
330: - Create a compatible DMStag with one dof per element (and nothing else).
331: - Create a constant vector and normalize it
332: - Migrate it to a vector on the original dmSol (making use of the fact
333: that this will fill in zeros for "extra" dof)
334: - Set the nullspace for the operator
335: - Destroy everything (the operator keeps the references it needs) */
336: static PetscErrorCode AttachNullspace(DM dmSol,Mat A)
337: {
339: DM dmPressure;
340: Vec constantPressure,basis;
341: PetscReal nrm;
342: MatNullSpace matNullSpace;
345: DMStagCreateCompatibleDMStag(dmSol,0,0,1,0,&dmPressure);
346: DMGetGlobalVector(dmPressure,&constantPressure);
347: VecSet(constantPressure,1.0);
348: VecNorm(constantPressure,NORM_2,&nrm);
349: VecScale(constantPressure,1.0/nrm);
350: DMCreateGlobalVector(dmSol,&basis);
351: DMStagMigrateVec(dmPressure,constantPressure,dmSol,basis);
352: MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol),PETSC_FALSE,1,&basis,&matNullSpace);
353: VecDestroy(&basis);
354: VecDestroy(&constantPressure);
355: MatSetNullSpace(A,matNullSpace);
356: MatNullSpaceDestroy(&matNullSpace);
357: return(0);
358: }
360: /* Create a reference solution.
361: Here, we use the more direct method of iterating over arrays. */
362: static PetscErrorCode CreateReferenceSolution(DM dmSol,Vec *pSolRef)
363: {
365: PetscInt startx,starty,nx,ny,nExtra[2],ex,ey;
366: PetscInt iuy,iux,ip,iprev,icenter;
367: PetscScalar ***arrSol,**cArrX,**cArrY;
368: Vec solRefLocal;
371: DMCreateGlobalVector(dmSol,pSolRef);
372: DMGetLocalVector(dmSol,&solRefLocal);
374: /* Obtain indices to use in the raw arrays */
375: DMStagGetLocationSlot(dmSol,DOWN,0,&iuy);
376: DMStagGetLocationSlot(dmSol,LEFT,0,&iux);
377: DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip);
379: /* Use high-level convenience functions to get raw arrays and indices for 1d coordinates */
380: DMStagGet1dCoordinateArraysDOFRead(dmSol,&cArrX,&cArrY,NULL);
381: DMStagGet1dCoordinateLocationSlot(dmSol,ELEMENT,&icenter);
382: DMStagGet1dCoordinateLocationSlot(dmSol,LEFT,&iprev);
384: DMStagVecGetArrayDOF(dmSol,solRefLocal,&arrSol);
385: DMStagGetCorners(dmSol,&startx,&starty,NULL,&nx,&ny,NULL,&nExtra[0],&nExtra[1],NULL);
386: for (ey=starty; ey<starty + ny + nExtra[1]; ++ey) {
387: for (ex=startx; ex<startx + nx + nExtra[0]; ++ex) {
388: arrSol[ey][ex][iuy] = uyRef(cArrX[ex][icenter],cArrY[ey][iprev]);
389: arrSol[ey][ex][iux] = uxRef(cArrX[ex][iprev],cArrY[ey][icenter]);
390: if (ey < starty+ny && ex < startx+nx) { /* Don't fill on the dummy elements (though you could, and these values would just be ignored) */
391: arrSol[ey][ex][ip] = pRef(cArrX[ex][icenter],cArrY[ey][icenter]);
392: }
393: }
394: }
395: DMStagVecRestoreArrayDOF(dmSol,solRefLocal,&arrSol);
396: DMStagRestore1dCoordinateArraysDOFRead(dmSol,&cArrX,&cArrY,NULL);
397: DMLocalToGlobal(dmSol,solRefLocal,INSERT_VALUES,*pSolRef);
398: DMRestoreLocalVector(dmSol,&solRefLocal);
399: return(0);
400: }
402: static PetscErrorCode CheckSolution(Vec sol,Vec solRef)
403: {
405: Vec diff;
406: PetscReal normsolRef,errAbs,errRel;
409: VecDuplicate(sol,&diff);
410: VecCopy(sol,diff);
411: VecAXPY(diff,-1.0,solRef);
412: VecNorm(diff,NORM_2,&errAbs);
413: VecNorm(solRef,NORM_2,&normsolRef);
414: errRel = errAbs/normsolRef;
415: PetscPrintf(PETSC_COMM_WORLD,"Error (abs): %g\nError (rel): %g\n",(double)errAbs,(double)errRel);
416: VecDestroy(&diff);
417: return(0);
418: }
420: /*TEST
422: test:
423: suffix: 1
424: nsize: 4
425: args: -ksp_monitor_short -ksp_converged_reason
427: test:
428: suffix: direct_umfpack
429: requires: suitesparse
430: nsize: 1
431: args: -pinpressure 1 -stag_grid_x 8 -stag_grid_y 6 -ksp_monitor_short -pc_type lu -pc_factor_mat_solver_type umfpack
433: TEST*/