Actual source code: ex3.c
1: static char help[] = "Solve a toy 3D 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 3D domain.
7: u_{xx} + u_{yy} + u_{zz} - p_x = f^x
8: v_{xx} + v_{yy} + u_{zz} - p_y = f^y
9: w_{xx} + w_{yy} + w_{zz} - p_y = f^z
10: u_x + v_y + w_z = g
12: g = 0 for the physical case.
14: Boundary conditions give prescribed flow perpendicular to the boundaries,
15: and zero derivative perpendicular to them (free slip). This involves
16: using a modifed stencil at the boundaries. Another option would be to
17: use DM_BOUNDARY_GHOSTED in DMStagCreate3d() and a matrix-free operator (MATSHELL)
18: making use of the uniformly-available ghost layer.
20: Use the -pinpressure option to fix a pressure node, instead of providing
21: a constant-pressure nullspace. This allows use of direct solvers, e.g. to
22: use UMFPACK,
24: ./ex3 -pinpressure 1 -pc_type lu -pc_factor_mat_solver_type umfpack
26: */
27: #include <petscdm.h>
28: #include <petscksp.h>
29: #include <petscdmstag.h>
31: /* Shorter, more convenient names for DMStagStencilLocation entries */
32: #define BACK_DOWN_LEFT DMSTAG_BACK_DOWN_LEFT
33: #define BACK_DOWN DMSTAG_BACK_DOWN
34: #define BACK_DOWN_RIGHT DMSTAG_BACK_DOWN_RIGHT
35: #define BACK_LEFT DMSTAG_BACK_LEFT
36: #define BACK DMSTAG_BACK
37: #define BACK_RIGHT DMSTAG_BACK_RIGHT
38: #define BACK_UP_LEFT DMSTAG_BACK_UP_LEFT
39: #define BACK_UP DMSTAG_BACK_UP
40: #define BACK_UP_RIGHT DMSTAG_BACK_UP_RIGHT
41: #define DOWN_LEFT DMSTAG_DOWN_LEFT
42: #define DOWN DMSTAG_DOWN
43: #define DOWN_RIGHT DMSTAG_DOWN_RIGHT
44: #define LEFT DMSTAG_LEFT
45: #define ELEMENT DMSTAG_ELEMENT
46: #define RIGHT DMSTAG_RIGHT
47: #define UP_LEFT DMSTAG_UP_LEFT
48: #define UP DMSTAG_UP
49: #define UP_RIGHT DMSTAG_UP_RIGHT
50: #define FRONT_DOWN_LEFT DMSTAG_FRONT_DOWN_LEFT
51: #define FRONT_DOWN DMSTAG_FRONT_DOWN
52: #define FRONT_DOWN_RIGHT DMSTAG_FRONT_DOWN_RIGHT
53: #define FRONT_LEFT DMSTAG_FRONT_LEFT
54: #define FRONT DMSTAG_FRONT
55: #define FRONT_RIGHT DMSTAG_FRONT_RIGHT
56: #define FRONT_UP_LEFT DMSTAG_FRONT_UP_LEFT
57: #define FRONT_UP DMSTAG_FRONT_UP
58: #define FRONT_UP_RIGHT DMSTAG_FRONT_UP_RIGHT
60: static PetscErrorCode CreateReferenceSolution(DM,Vec*);
61: static PetscErrorCode CreateSystem(DM,Mat*,Vec*,PetscBool);
62: static PetscErrorCode AttachNullspace(DM,Mat);
63: static PetscErrorCode CheckSolution(Vec,Vec);
65: /* Manufactured solution. Chosen to be higher order than can be solved exactly,
66: and to have a zero derivative for flow parallel to the boundaries. That is,
67: d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
68: and left boundaries.
69: These expressions could be made more interesting with more z terms,
70: and convergence could be confirmed. */
72: static PetscScalar uxRef(PetscScalar x,PetscScalar y, PetscScalar z) {return 0.0*x + y*y - 2.0*y*y*y + y*y*y*y + 0.0*z;}
73: static PetscScalar uyRef(PetscScalar x,PetscScalar y, PetscScalar z) {return x*x - 2.0*x*x*x + x*x*x*x +0.0*y + 0.0*z;}
74: static PetscScalar uzRef(PetscScalar x,PetscScalar y, PetscScalar z) {return 0.0*x + 0.0*y + 0.0*z + 1.0;}
75: static PetscScalar pRef (PetscScalar x,PetscScalar y, PetscScalar z) {return -1.0*(x-0.5) + -3.0/2.0*y*y + 0.5 -2.0*(z-1.0);} /* zero integral */
76: static PetscScalar fx (PetscScalar x,PetscScalar y, PetscScalar z) {return 0.0*x + 2.0 -12.0*y + 12.0*y*y + 0.0*z + 1.0;}
77: static PetscScalar fy (PetscScalar x,PetscScalar y, PetscScalar z) {return 2.0 -12.0*x + 12.0*x*x + 3.0*y + 0.0*z;}
78: static PetscScalar fz (PetscScalar x,PetscScalar y, PetscScalar z) {return 0.0*x + 0.0*y + 0.0*z + 2.0;}
79: static PetscScalar g (PetscScalar x,PetscScalar y, PetscScalar z) {return 0.0*x + 0.0*y + 0.0*z + 0.0;}
81: int main(int argc,char **argv)
82: {
84: DM dmSol;
85: Vec sol,solRef,rhs;
86: Mat A;
87: KSP ksp;
88: PC pc;
89: PetscBool pinPressure;
91: /* Initialize PETSc and process command line arguments */
92: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
93: pinPressure = PETSC_TRUE;
94: PetscOptionsGetBool(NULL,NULL,"-pinpressure",&pinPressure,NULL);
96: /* Create 3D DMStag for the solution, and set up. */
97: {
98: const PetscInt dof0 = 0, dof1 = 0,dof2 = 1, dof3 = 1; /* 1 dof on each face and element center */
99: const PetscInt stencilWidth = 1;
100: DMStagCreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,4,5,6,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof0,dof1,dof2,dof3,DMSTAG_STENCIL_BOX,stencilWidth,NULL,NULL,NULL,&dmSol);
101: DMSetFromOptions(dmSol);
102: DMSetUp(dmSol);
103: DMStagSetUniformCoordinatesExplicit(dmSol,0.0,1.0,0.0,1.0,0.0,1.0);
104: /* Note: also see ex2.c, where another, more efficient option is demonstrated,
105: using DMStagSetUniformCoordinatesProduct() */
106: }
108: /* Compute (manufactured) reference solution */
109: CreateReferenceSolution(dmSol,&solRef);
111: /* Assemble system */
112: CreateSystem(dmSol,&A,&rhs,pinPressure);
114: /* Attach a constant-pressure nullspace to the operator */
115: if (!pinPressure) {
116: AttachNullspace(dmSol,A);
117: }
119: /* Solve */
120: DMCreateGlobalVector(dmSol,&sol);
121: KSPCreate(PETSC_COMM_WORLD,&ksp);
122: KSPSetType(ksp,KSPFGMRES);
123: KSPSetOperators(ksp,A,A);
124: KSPGetPC(ksp,&pc);
125: PCSetType(pc,PCFIELDSPLIT);
126: PCFieldSplitSetDetectSaddlePoint(pc,PETSC_TRUE);
127: KSPSetFromOptions(ksp);
128: KSPSolve(ksp,rhs,sol);
130: /* Check Solution */
131: CheckSolution(sol,solRef);
133: /* Clean up and finalize PETSc */
134: KSPDestroy(&ksp);
135: VecDestroy(&sol);
136: VecDestroy(&solRef);
137: VecDestroy(&rhs);
138: MatDestroy(&A);
139: DMDestroy(&dmSol);
140: PetscFinalize();
141: return ierr;
142: }
144: static PetscErrorCode CreateSystem(DM dmSol,Mat *pA,Vec *pRhs, PetscBool pinPressure)
145: {
146: PetscErrorCode ierr;
147: Vec rhs,coordLocal;
148: Mat A;
149: PetscInt startx,starty,startz,N[3],nx,ny,nz,ex,ey,ez,d;
150: PetscInt icp[3],icux[3],icuy[3],icuz[3],icux_right[3],icuy_up[3],icuz_front[3];
151: PetscReal hx,hy,hz;
152: DM dmCoord;
153: PetscScalar ****arrCoord;
156: DMCreateMatrix(dmSol,pA);
157: A = *pA;
158: DMCreateGlobalVector(dmSol,pRhs);
159: rhs = *pRhs;
161: DMStagGetCorners(dmSol,&startx,&starty,&startz,&nx,&ny,&nz,NULL,NULL,NULL);
162: DMStagGetGlobalSizes(dmSol,&N[0],&N[1],&N[2]);
163: if (N[0] < 2 || N[1] < 2 || N[2] < 2) SETERRQ(PetscObjectComm((PetscObject)dmSol),PETSC_ERR_ARG_SIZ,"This example requires at least two elements in each dimensions");
164: hx = 1.0/N[0]; hy = 1.0/N[1]; hz = 1.0/N[2];
165: DMGetCoordinateDM(dmSol,&dmCoord);
166: DMGetCoordinatesLocal(dmSol,&coordLocal);
167: DMStagVecGetArrayRead(dmCoord,coordLocal,&arrCoord);
168: for (d=0; d<3; ++d) {
169: DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d]);
170: DMStagGetLocationSlot(dmCoord,LEFT, d,&icux[d]);
171: DMStagGetLocationSlot(dmCoord,DOWN, d,&icuy[d]);
172: DMStagGetLocationSlot(dmCoord,BACK, d,&icuz[d]);
173: DMStagGetLocationSlot(dmCoord,RIGHT, d,&icux_right[d]);
174: DMStagGetLocationSlot(dmCoord,UP, d,&icuy_up[d]);
175: DMStagGetLocationSlot(dmCoord,FRONT, d,&icuz_front[d]);
176: }
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 (ez = startz; ez<startz+nz; ++ez) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
181: for (ey = starty; ey<starty+ny; ++ey) {
182: for (ex = startx; ex<startx+nx; ++ex) {
184: if (ex == N[0]-1) {
185: /* Right Boundary velocity Dirichlet */
186: DMStagStencil row;
187: PetscScalar valRhs;
188: const PetscScalar valA = 1.0;
189: row.i = ex; row.j = ey; row.k = ez; row.loc = RIGHT; row.c = 0;
190: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
191: valRhs = uxRef(arrCoord[ez][ey][ex][icux_right[0]], arrCoord[ez][ey][ex][icux_right[1]], arrCoord[ez][ey][ex][icux_right[2]]);
192: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
193: }
194: if (ey == N[1]-1) {
195: /* Top boundary velocity Dirichlet */
196: DMStagStencil row;
197: PetscScalar valRhs;
198: const PetscScalar valA = 1.0;
199: row.i = ex; row.j = ey; row.k = ez; row.loc = UP; row.c = 0;
200: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
201: valRhs = uyRef(arrCoord[ez][ey][ex][icuy_up[0]],arrCoord[ez][ey][ex][icuy_up[1]],arrCoord[ez][ey][ex][icuy_up[2]]);
202: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
203: }
204: if (ez == N[2]-1) {
205: /* Front boundary velocity Dirichlet */
206: DMStagStencil row;
207: PetscScalar valRhs;
208: const PetscScalar valA = 1.0;
209: row.i = ex; row.j = ey; row.k = ez; row.loc = FRONT; row.c = 0;
210: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
211: valRhs = uzRef(arrCoord[ez][ey][ex][icuz_front[0]],arrCoord[ez][ey][ex][icuz_front[1]],arrCoord[ez][ey][ex][icuz_front[2]]);
212: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
213: }
215: /* Equation on left face of this element */
216: if (ex == 0) {
217: /* Left velocity Dirichlet */
218: DMStagStencil row;
219: PetscScalar valRhs;
220: const PetscScalar valA = 1.0;
221: row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
222: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
223: valRhs = uxRef(arrCoord[ez][ey][ex][icux[0]],arrCoord[ez][ey][ex][icux[1]],arrCoord[ez][ey][ex][icux[2]]);
224: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
225: } else {
226: /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
227: DMStagStencil row,col[9];
228: PetscScalar valA[9],valRhs;
229: PetscInt nEntries;
231: row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
232: if (ey == 0) {
233: if (ez == 0) {
234: nEntries = 7;
235: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -1.0 / (hy*hy) -1.0 / (hz*hz);
236: /* Missing down term */
237: col[1].i = ex ; col[1].j = ey+1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
238: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
239: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
240: /* Missing back term */
241: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
242: col[5].i = ex-1; col[5].j = ey ; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hx;
243: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hx;
244: } else if (ez == N[2]-1) {
245: nEntries = 7;
246: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -1.0 / (hy*hy) -1.0 / (hz*hz);
247: /* Missing down term */
248: col[1].i = ex ; col[1].j = ey+1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
249: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
250: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
251: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
252: /* Missing front term */
253: col[5].i = ex-1; col[5].j = ey ; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hx;
254: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hx;
255: } else {
256: nEntries = 8;
257: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
258: /* Missing down term */
259: col[1].i = ex ; col[1].j = ey+1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
260: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
261: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
262: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
263: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = LEFT; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
264: col[6].i = ex-1; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hx;
265: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hx;
266: }
267: } else if (ey == N[1]-1) {
268: if (ez == 0) {
269: nEntries = 7;
270: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
271: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
272: /* Missing up term */
273: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; 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].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
275: /* Missing back entry */
276: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
277: col[5].i = ex-1; col[5].j = ey ; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hx;
278: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hx;
279: } else if (ez == N[2]-1) {
280: nEntries = 7;
281: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
282: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
283: /* Missing up term */
284: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
285: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
286: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
287: /* Missing front term */
288: col[5].i = ex-1; col[5].j = ey ; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hx;
289: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hx;
290: } else {
291: nEntries = 8;
292: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
293: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
294: /* Missing up term */
295: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
296: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
297: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
298: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = LEFT; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
299: col[6].i = ex-1; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hx;
300: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hx;
301: }
302: } else if (ez == 0) {
303: nEntries = 8;
304: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
305: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
306: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
307: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
308: col[4].i = ex+1; col[4].j = ey ; col[4].k = ez ; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
309: /* Missing back term */
310: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = LEFT; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
311: col[6].i = ex-1; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hx;
312: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hx;
313: } else if (ez == N[2]-1) {
314: nEntries = 8;
315: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
316: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
317: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
318: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
319: col[4].i = ex+1; col[4].j = ey ; col[4].k = ez ; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
320: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = LEFT; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
321: /* Missing front term */
322: col[6].i = ex-1; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hx;
323: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hx;
324: } else {
325: nEntries = 9;
326: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = LEFT; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
327: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = LEFT; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
328: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = LEFT; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
329: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = LEFT; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
330: col[4].i = ex+1; col[4].j = ey ; col[4].k = ez ; col[4].loc = LEFT; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
331: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = LEFT; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
332: col[6].i = ex ; col[6].j = ey ; col[6].k = ez+1; col[6].loc = LEFT; col[6].c = 0; valA[6] = 1.0 / (hz*hz);
333: col[7].i = ex-1; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = 1.0 / hx;
334: col[8].i = ex ; col[8].j = ey ; col[8].k = ez ; col[8].loc = ELEMENT; col[8].c = 0; valA[8] = -1.0 / hx;
335: }
336: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
337: valRhs = fx(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
338: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
339: }
341: /* Equation on bottom face of this element */
342: if (ey == 0) {
343: /* Bottom boundary velocity Dirichlet */
344: DMStagStencil row;
345: PetscScalar valRhs;
346: const PetscScalar valA = 1.0;
347: row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
348: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
349: valRhs = uyRef(arrCoord[ez][ey][ex][icuy[0]],arrCoord[ez][ey][ex][icuy[1]],arrCoord[ez][ey][ex][icuy[2]]);
350: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
351: } else {
352: /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
353: DMStagStencil row,col[9];
354: PetscScalar valA[9],valRhs;
355: PetscInt nEntries;
357: row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
358: if (ex ==0) {
359: if (ez == 0) {
360: nEntries = 7;
361: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
362: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
363: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
364: /* Left term missing */
365: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
366: /* Back term missing */
367: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
368: col[5].i = ex ; col[5].j = ey-1; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hy;
369: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hy;
370: } else if (ez == N[2]-1) {
371: nEntries = 7;
372: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
373: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
374: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
375: /* Left term missing */
376: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
377: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
378: /* Front term missing */
379: col[5].i = ex ; col[5].j = ey-1; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hy;
380: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hy;
381: } else {
382: nEntries = 8;
383: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
384: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
385: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
386: /* Left term missing */
387: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
388: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
389: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = DOWN; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
390: col[6].i = ex ; col[6].j = ey-1; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hy;
391: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hy;
392: }
393: } else if (ex == N[0]-1) {
394: if (ez == 0) {
395: nEntries = 7;
396: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
397: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
398: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
399: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
400: /* Right term missing */
401: /* Back term missing */
402: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
403: col[5].i = ex ; col[5].j = ey-1; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hy;
404: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hy;
405: } else if (ez == N[2]-1) {
406: nEntries = 7;
407: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
408: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
409: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
410: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
411: /* Right term missing */
412: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
413: /* Front term missing */
414: col[5].i = ex ; col[5].j = ey-1; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hy;
415: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hy;
416: } else {
417: nEntries = 8;
418: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
419: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
420: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
421: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
422: /* Right term missing */
423: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
424: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = DOWN; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
425: col[6].i = ex ; col[6].j = ey-1; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hy;
426: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hy;
427: }
428: } else if (ez == 0) {
429: nEntries = 8;
430: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
431: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
432: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
433: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
434: col[4].i = ex+1; col[4].j = ey ; col[4].k = ez ; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
435: /* Back term missing */
436: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = DOWN; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
437: col[6].i = ex ; col[6].j = ey-1; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hy;
438: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hy;
439: } else if (ez == N[2]-1) {
440: nEntries = 8;
441: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -1.0 / (hz*hz);
442: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
443: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
444: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
445: col[4].i = ex+1; col[4].j = ey ; col[4].k = ez ; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
446: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = DOWN; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
447: /* Front term missing */
448: col[6].i = ex ; col[6].j = ey-1; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hy;
449: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hy;
450: } else {
451: nEntries = 9;
452: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = DOWN; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
453: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = DOWN; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
454: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = DOWN; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
455: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = DOWN; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
456: col[4].i = ex+1; col[4].j = ey ; col[4].k = ez ; col[4].loc = DOWN; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
457: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = DOWN; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
458: col[6].i = ex ; col[6].j = ey ; col[6].k = ez+1; col[6].loc = DOWN; col[6].c = 0; valA[6] = 1.0 / (hz*hz);
459: col[7].i = ex ; col[7].j = ey-1; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = 1.0 / hy;
460: col[8].i = ex ; col[8].j = ey ; col[8].k = ez ; col[8].loc = ELEMENT; col[8].c = 0; valA[8] = -1.0 / hy;
461: }
462: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
463: valRhs = fy(arrCoord[ez][ey][ex][icuy[0]],arrCoord[ez][ey][ex][icuy[1]],arrCoord[ez][ey][ex][icuy[2]]);
464: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
465: }
467: /* Equation on back face of this element */
468: if (ez == 0) {
469: /* Back boundary velocity Dirichlet */
470: DMStagStencil row;
471: PetscScalar valRhs;
472: const PetscScalar valA = 1.0;
473: row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
474: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
475: valRhs = uzRef(arrCoord[ez][ey][ex][icuz[0]],arrCoord[ez][ey][ex][icuz[1]],arrCoord[ez][ey][ex][icuz[2]]);
476: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
477: } else {
478: /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
479: DMStagStencil row,col[9];
480: PetscScalar valA[9],valRhs;
481: PetscInt nEntries;
483: row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
484: if (ex == 0) {
485: if (ey == 0) {
486: nEntries = 7;
487: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
488: /* Down term missing */
489: col[1].i = ex ; col[1].j = ey+1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
490: /* Left term missing */
491: col[2].i = ex+1; col[2].j = ey ; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
492: col[3].i = ex ; col[3].j = ey ; col[3].k = ez-1; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hz*hz);
493: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
494: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hz;
495: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hz;
496: } else if (ey == N[1]-1) {
497: nEntries = 7;
498: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
499: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
500: /* Up term missing */
501: /* Left term missing */
502: col[2].i = ex+1; col[2].j = ey ; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
503: col[3].i = ex ; col[3].j = ey ; col[3].k = ez-1; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hz*hz);
504: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
505: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hz;
506: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hz;
507: } else {
508: nEntries = 8;
509: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
510: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
511: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
512: /* Left term missing */
513: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
514: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
515: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = BACK; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
516: col[6].i = ex ; col[6].j = ey ; col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hz;
517: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hz;
518: }
519: } else if (ex == N[0]-1) {
520: if (ey == 0) {
521: nEntries = 7;
522: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
523: /* Down term missing */
524: col[1].i = ex ; col[1].j = ey+1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
525: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
526: /* Right term missing */
527: col[3].i = ex ; col[3].j = ey ; col[3].k = ez-1; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hz*hz);
528: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
529: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hz;
530: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hz;
531: } else if (ey == N[1]-1) {
532: nEntries = 7;
533: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
534: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
535: /* Up term missing */
536: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
537: /* Right term missing */
538: col[3].i = ex ; col[3].j = ey ; col[3].k = ez-1; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hz*hz);
539: col[4].i = ex ; col[4].j = ey ; col[4].k = ez+1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
540: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hz;
541: col[6].i = ex ; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hz;
542: } else {
543: nEntries = 8;
544: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
545: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
546: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
547: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
548: /* Right term missing */
549: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
550: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = BACK; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
551: col[6].i = ex ; col[6].j = ey ; col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hz;
552: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hz;
553: }
554: } else if (ey == 0) {
555: nEntries = 8;
556: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
557: /* Down term missing */
558: col[1].i = ex ; col[1].j = ey+1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
559: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
560: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
561: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
562: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = BACK; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
563: col[6].i = ex ; col[6].j = ey ; col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hz;
564: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hz;
565: } else if (ey == N[1]-1) {
566: nEntries = 8;
567: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
568: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
569: /* Up term missing */
570: col[2].i = ex-1; col[2].j = ey ; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hx*hx);
571: col[3].i = ex+1; col[3].j = ey ; col[3].k = ez ; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
572: col[4].i = ex ; col[4].j = ey ; col[4].k = ez-1; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hz*hz);
573: col[5].i = ex ; col[5].j = ey ; col[5].k = ez+1; col[5].loc = BACK; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
574: col[6].i = ex ; col[6].j = ey ; col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hz;
575: col[7].i = ex ; col[7].j = ey ; col[7].k = ez ; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = -1.0 / hz;
576: } else {
577: nEntries = 9;
578: col[0].i = ex ; col[0].j = ey ; col[0].k = ez ; col[0].loc = BACK; col[0].c = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
579: col[1].i = ex ; col[1].j = ey-1; col[1].k = ez ; col[1].loc = BACK; col[1].c = 0; valA[1] = 1.0 / (hy*hy);
580: col[2].i = ex ; col[2].j = ey+1; col[2].k = ez ; col[2].loc = BACK; col[2].c = 0; valA[2] = 1.0 / (hy*hy);
581: col[3].i = ex-1; col[3].j = ey ; col[3].k = ez ; col[3].loc = BACK; col[3].c = 0; valA[3] = 1.0 / (hx*hx);
582: col[4].i = ex+1; col[4].j = ey ; col[4].k = ez ; col[4].loc = BACK; col[4].c = 0; valA[4] = 1.0 / (hx*hx);
583: col[5].i = ex ; col[5].j = ey ; col[5].k = ez-1; col[5].loc = BACK; col[5].c = 0; valA[5] = 1.0 / (hz*hz);
584: col[6].i = ex ; col[6].j = ey ; col[6].k = ez+1; col[6].loc = BACK; col[6].c = 0; valA[6] = 1.0 / (hz*hz);
585: col[7].i = ex ; col[7].j = ey ; col[7].k = ez-1; col[7].loc = ELEMENT; col[7].c = 0; valA[7] = 1.0 / hz;
586: col[8].i = ex ; col[8].j = ey ; col[8].k = ez ; col[8].loc = ELEMENT; col[8].c = 0; valA[8] = -1.0 / hz;
587: }
588: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
589: valRhs = fz(arrCoord[ez][ey][ex][icuz[0]],arrCoord[ez][ey][ex][icuz[1]],arrCoord[ez][ey][ex][icuz[2]]);
590: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
591: }
593: /* P equation : u_x + v_y + w_z = g
594: Note that this includes an explicit zero on the diagonal. This is only needed for
595: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
596: if (pinPressure && ex == 0 && ey == 0 && ez == 0) { /* Pin the first pressure node, if requested */
597: DMStagStencil row;
598: PetscScalar valA,valRhs;
599: row.i = ex; row.j = ey; row.k = ez; row.loc = ELEMENT; row.c = 0;
600: valA = 1.0;
601: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
602: valRhs = pRef(arrCoord[ez][ey][ex][icp[0]],arrCoord[ez][ey][ex][icp[1]],arrCoord[ez][ey][ex][icp[2]]);
603: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
604: } else {
605: DMStagStencil row,col[7];
606: PetscScalar valA[7],valRhs;
608: row.i = ex; row.j = ey; row.k = ez; row.loc = ELEMENT; row.c = 0;
609: col[0].i = ex; col[0].j = ey; col[0].k = ez; col[0].loc = LEFT; col[0].c = 0; valA[0] = -1.0 / hx;
610: col[1].i = ex; col[1].j = ey; col[1].k = ez; col[1].loc = RIGHT; col[1].c = 0; valA[1] = 1.0 / hx;
611: col[2].i = ex; col[2].j = ey; col[2].k = ez; col[2].loc = DOWN; col[2].c = 0; valA[2] = -1.0 / hy;
612: col[3].i = ex; col[3].j = ey; col[3].k = ez; col[3].loc = UP; col[3].c = 0; valA[3] = 1.0 / hy;
613: col[4].i = ex; col[4].j = ey; col[4].k = ez; col[4].loc = BACK; col[4].c = 0; valA[4] = -1.0 / hz;
614: col[5].i = ex; col[5].j = ey; col[5].k = ez; col[5].loc = FRONT; col[5].c = 0; valA[5] = 1.0 / hz;
615: col[6] = row; valA[6] = 0.0;
616: DMStagMatSetValuesStencil(dmSol,A,1,&row,7,col,valA,INSERT_VALUES);
617: valRhs = g(arrCoord[ez][ey][ex][icp[0]],arrCoord[ez][ey][ex][icp[1]],arrCoord[ez][ey][ex][icp[2]]);
618: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
619: }
620: }
621: }
622: }
623: DMStagVecRestoreArrayRead(dmCoord,coordLocal,&arrCoord);
624: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
625: VecAssemblyBegin(rhs);
626: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
627: VecAssemblyEnd(rhs);
629: return(0);
630: }
632: /* Create a pressure-only DMStag and use it to generate a nullspace vector
633: - Create a compatible DMStag with one dof per element (and nothing else).
634: - Create a constant vector and normalize it
635: - Migrate it to a vector on the original dmSol (making use of the fact
636: that this will fill in zeros for "extra" dof)
637: - Set the nullspace for the operator
638: - Destroy everything (the operator keeps the references it needs) */
639: static PetscErrorCode AttachNullspace(DM dmSol,Mat A)
640: {
642: DM dmPressure;
643: Vec constantPressure,basis;
644: PetscReal nrm;
645: MatNullSpace matNullSpace;
648: DMStagCreateCompatibleDMStag(dmSol,0,0,0,1,&dmPressure);
649: DMGetGlobalVector(dmPressure,&constantPressure);
650: VecSet(constantPressure,1.0);
651: VecNorm(constantPressure,NORM_2,&nrm);
652: VecScale(constantPressure,1.0/nrm);
653: DMCreateGlobalVector(dmSol,&basis);
654: DMStagMigrateVec(dmPressure,constantPressure,dmSol,basis);
655: MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol),PETSC_FALSE,1,&basis,&matNullSpace);
656: VecDestroy(&basis);
657: VecDestroy(&constantPressure);
658: MatSetNullSpace(A,matNullSpace);
659: MatNullSpaceDestroy(&matNullSpace);
660: return(0);
661: }
663: static PetscErrorCode CreateReferenceSolution(DM dmSol,Vec *pSolRef)
664: {
665: PetscErrorCode ierr;
666: PetscInt start[3],n[3],nExtra[3],ex,ey,ez,d;
667: PetscInt ip,iux,iuy,iuz,icp[3],icux[3],icuy[3],icuz[3];
668: Vec solRef,solRefLocal,coord,coordLocal;
669: DM dmCoord;
670: PetscScalar ****arrSol,****arrCoord;
673: DMCreateGlobalVector(dmSol,pSolRef);
674: solRef = *pSolRef;
675: DMStagGetCorners(dmSol,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
676: DMGetCoordinateDM(dmSol,&dmCoord);
677: DMGetCoordinates(dmSol,&coord);
678: DMGetLocalVector(dmCoord,&coordLocal);
679: DMGlobalToLocal(dmCoord,coord,INSERT_VALUES,coordLocal);
680: DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip);
681: DMStagGetLocationSlot(dmSol,LEFT, 0,&iux);
682: DMStagGetLocationSlot(dmSol,DOWN, 0,&iuy);
683: DMStagGetLocationSlot(dmSol,BACK, 0,&iuz);
684: for (d=0; d<3; ++d) {
685: DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d]);
686: DMStagGetLocationSlot(dmCoord,LEFT, d,&icux[d]);
687: DMStagGetLocationSlot(dmCoord,DOWN, d,&icuy[d]);
688: DMStagGetLocationSlot(dmCoord,BACK, d,&icuz[d]);
689: }
690: DMStagVecGetArrayRead(dmCoord,coordLocal,&arrCoord);
691: DMGetLocalVector(dmSol,&solRefLocal);
692: DMStagVecGetArray(dmSol,solRefLocal,&arrSol);
693: for (ez=start[2]; ez<start[2] + n[2] + nExtra[2]; ++ez) {
694: for (ey=start[1]; ey<start[1] + n[1] + nExtra[1]; ++ey) {
695: for (ex=start[0]; ex<start[0] + n[0] + nExtra[0]; ++ex) {
696: if (ex < start[0] + n[0] && ey < start[1] + n[1]) {
697: arrSol[ez][ey][ex][iuz] = uzRef(
698: arrCoord[ez][ey][ex][icuz[0]],
699: arrCoord[ez][ey][ex][icuz[1]],
700: arrCoord[ez][ey][ex][icuz[2]]);
701: }
702: if (ex < start[0] + n[0] && ey < start[2] + n[2]) {
703: arrSol[ez][ey][ex][iuy] = uyRef(
704: arrCoord[ez][ey][ex][icuy[0]],
705: arrCoord[ez][ey][ex][icuy[1]],
706: arrCoord[ez][ey][ex][icuy[2]]);
707: }
708: if (ex < start[1] + n[1] && ey < start[2] + n[2]) {
709: arrSol[ez][ey][ex][iux] = uxRef(
710: arrCoord[ez][ey][ex][icux[0]],
711: arrCoord[ez][ey][ex][icux[1]],
712: arrCoord[ez][ey][ex][icux[2]]);
713: }
714: if (ex < start[0] + n[0] && ey < start[1] + n[1] && ez < start[2] + n[2]) {
715: arrSol[ez][ey][ex][ip] = pRef(
716: arrCoord[ez][ey][ex][icp[0]],
717: arrCoord[ez][ey][ex][icp[1]],
718: arrCoord[ez][ey][ex][icp[2]]);
719: }
720: }
721: }
722: }
723: DMStagVecRestoreArrayRead(dmCoord,coordLocal,&arrCoord);
724: DMStagVecRestoreArray(dmSol,solRefLocal,&arrSol);
725: DMLocalToGlobal(dmSol,solRefLocal,INSERT_VALUES,solRef);
726: DMRestoreLocalVector(dmCoord,&coordLocal);
727: DMRestoreLocalVector(dmSol,&solRefLocal);
728: return(0);
729: }
731: static PetscErrorCode CheckSolution(Vec sol,Vec solRef)
732: {
734: Vec diff;
735: PetscReal normsolRef,errAbs,errRel;
738: VecDuplicate(sol,&diff);
739: VecCopy(sol,diff);
740: VecAXPY(diff,-1.0,solRef);
741: VecNorm(diff,NORM_2,&errAbs);
742: VecNorm(solRef,NORM_2,&normsolRef);
743: errRel = errAbs/normsolRef;
744: PetscPrintf(PETSC_COMM_WORLD,"Error (abs): %g\nError (rel): %g\n",(double)errAbs,(double)errRel);
745: VecDestroy(&diff);
746: return(0);
747: }
749: /*TEST
751: test:
752: suffix: 1
753: requires: mumps
754: nsize: 27
755: args: -ksp_monitor_short -ksp_converged_reason -stag_ranks_x 3 -stag_ranks_y 3 -stag_ranks_z 3 -pc_fieldsplit_schur_fact_type diag -fieldsplit_0_ksp_type preonly -fieldsplit_1_pc_type none -fieldsplit_0_pc_type lu -fieldsplit_1_ksp_type gmres -fieldsplit_1_ksp_max_it 20
757: test:
758: suffix: 2
759: requires: !single
760: nsize: 4
761: args: -ksp_monitor_short -ksp_converged_reason -pc_fieldsplit_schur_fact_type diag -fieldsplit_0_ksp_type preonly -fieldsplit_1_pc_type none -fieldsplit_0_pc_type gamg -fieldsplit_0_mg_levels_ksp_max_it 3 -fieldsplit_1_ksp_type gmres -fieldsplit_1_ksp_max_it 20 -fieldsplit_0_pc_gamg_esteig_ksp_type gmres
763: test:
764: suffix: direct_umfpack
765: requires: suitesparse
766: nsize: 1
767: args: -pinpressure 1 -stag_grid_x 5 -stag_grid_y 3 -stag_grid_z 4 -ksp_monitor_short -pc_type lu -pc_factor_mat_solver_type umfpack
769: TEST*/