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: {
83: DM dmSol;
84: Vec sol,solRef,rhs;
85: Mat A;
86: KSP ksp;
87: PC pc;
88: PetscBool pinPressure;
90: /* Initialize PETSc and process command line arguments */
91: PetscInitialize(&argc,&argv,(char*)0,help);
92: pinPressure = PETSC_TRUE;
93: PetscOptionsGetBool(NULL,NULL,"-pinpressure",&pinPressure,NULL);
95: /* Create 3D DMStag for the solution, and set up. */
96: {
97: const PetscInt dof0 = 0, dof1 = 0,dof2 = 1, dof3 = 1; /* 1 dof on each face and element center */
98: const PetscInt stencilWidth = 1;
99: 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);
100: DMSetFromOptions(dmSol);
101: DMSetUp(dmSol);
102: DMStagSetUniformCoordinatesExplicit(dmSol,0.0,1.0,0.0,1.0,0.0,1.0);
103: /* Note: also see ex2.c, where another, more efficient option is demonstrated,
104: using DMStagSetUniformCoordinatesProduct() */
105: }
107: /* Compute (manufactured) reference solution */
108: CreateReferenceSolution(dmSol,&solRef);
110: /* Assemble system */
111: CreateSystem(dmSol,&A,&rhs,pinPressure);
113: /* Attach a constant-pressure nullspace to the operator */
114: if (!pinPressure) {
115: AttachNullspace(dmSol,A);
116: }
118: /* Solve */
119: DMCreateGlobalVector(dmSol,&sol);
120: KSPCreate(PETSC_COMM_WORLD,&ksp);
121: KSPSetType(ksp,KSPFGMRES);
122: KSPSetOperators(ksp,A,A);
123: KSPGetPC(ksp,&pc);
124: PCSetType(pc,PCFIELDSPLIT);
125: PCFieldSplitSetDetectSaddlePoint(pc,PETSC_TRUE);
126: KSPSetFromOptions(ksp);
127: KSPSolve(ksp,rhs,sol);
129: /* Check Solution */
130: CheckSolution(sol,solRef);
132: /* Clean up and finalize PETSc */
133: KSPDestroy(&ksp);
134: VecDestroy(&sol);
135: VecDestroy(&solRef);
136: VecDestroy(&rhs);
137: MatDestroy(&A);
138: DMDestroy(&dmSol);
139: PetscFinalize();
140: return 0;
141: }
143: static PetscErrorCode CreateSystem(DM dmSol,Mat *pA,Vec *pRhs, PetscBool pinPressure)
144: {
145: Vec rhs,coordLocal;
146: Mat A;
147: PetscInt startx,starty,startz,N[3],nx,ny,nz,ex,ey,ez,d;
148: PetscInt icp[3],icux[3],icuy[3],icuz[3],icux_right[3],icuy_up[3],icuz_front[3];
149: PetscReal hx,hy,hz;
150: DM dmCoord;
151: PetscScalar ****arrCoord;
154: DMCreateMatrix(dmSol,pA);
155: A = *pA;
156: DMCreateGlobalVector(dmSol,pRhs);
157: rhs = *pRhs;
159: DMStagGetCorners(dmSol,&startx,&starty,&startz,&nx,&ny,&nz,NULL,NULL,NULL);
160: DMStagGetGlobalSizes(dmSol,&N[0],&N[1],&N[2]);
162: hx = 1.0/N[0]; hy = 1.0/N[1]; hz = 1.0/N[2];
163: DMGetCoordinateDM(dmSol,&dmCoord);
164: DMGetCoordinatesLocal(dmSol,&coordLocal);
165: DMStagVecGetArrayRead(dmCoord,coordLocal,&arrCoord);
166: for (d=0; d<3; ++d) {
167: DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d]);
168: DMStagGetLocationSlot(dmCoord,LEFT, d,&icux[d]);
169: DMStagGetLocationSlot(dmCoord,DOWN, d,&icuy[d]);
170: DMStagGetLocationSlot(dmCoord,BACK, d,&icuz[d]);
171: DMStagGetLocationSlot(dmCoord,RIGHT, d,&icux_right[d]);
172: DMStagGetLocationSlot(dmCoord,UP, d,&icuy_up[d]);
173: DMStagGetLocationSlot(dmCoord,FRONT, d,&icuz_front[d]);
174: }
176: /* Loop over all local elements. Note that it may be more efficient in real
177: applications to loop over each boundary separately */
178: for (ez = startz; ez<startz+nz; ++ez) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
179: for (ey = starty; ey<starty+ny; ++ey) {
180: 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; row.j = ey; row.k = ez; row.loc = RIGHT; row.c = 0;
188: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
189: valRhs = uxRef(arrCoord[ez][ey][ex][icux_right[0]], arrCoord[ez][ey][ex][icux_right[1]], arrCoord[ez][ey][ex][icux_right[2]]);
190: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
191: }
192: if (ey == N[1]-1) {
193: /* Top boundary velocity Dirichlet */
194: DMStagStencil row;
195: PetscScalar valRhs;
196: const PetscScalar valA = 1.0;
197: row.i = ex; row.j = ey; row.k = ez; row.loc = UP; row.c = 0;
198: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
199: valRhs = uyRef(arrCoord[ez][ey][ex][icuy_up[0]],arrCoord[ez][ey][ex][icuy_up[1]],arrCoord[ez][ey][ex][icuy_up[2]]);
200: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
201: }
202: if (ez == N[2]-1) {
203: /* Front boundary velocity Dirichlet */
204: DMStagStencil row;
205: PetscScalar valRhs;
206: const PetscScalar valA = 1.0;
207: row.i = ex; row.j = ey; row.k = ez; row.loc = FRONT; row.c = 0;
208: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
209: valRhs = uzRef(arrCoord[ez][ey][ex][icuz_front[0]],arrCoord[ez][ey][ex][icuz_front[1]],arrCoord[ez][ey][ex][icuz_front[2]]);
210: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
211: }
213: /* Equation on left face of this element */
214: if (ex == 0) {
215: /* Left velocity Dirichlet */
216: DMStagStencil row;
217: PetscScalar valRhs;
218: const PetscScalar valA = 1.0;
219: row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
220: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
221: valRhs = uxRef(arrCoord[ez][ey][ex][icux[0]],arrCoord[ez][ey][ex][icux[1]],arrCoord[ez][ey][ex][icux[2]]);
222: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
223: } else {
224: /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
225: DMStagStencil row,col[9];
226: PetscScalar valA[9],valRhs;
227: PetscInt nEntries;
229: row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
230: if (ey == 0) {
231: if (ez == 0) {
232: nEntries = 7;
233: 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);
234: /* Missing down term */
235: 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);
236: 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);
237: 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);
238: /* Missing back term */
239: 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);
240: 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;
241: 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;
242: } else if (ez == N[2]-1) {
243: nEntries = 7;
244: 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);
245: /* Missing down term */
246: 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);
247: 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);
248: 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);
249: 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);
250: /* Missing front term */
251: 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;
252: 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;
253: } else {
254: nEntries = 8;
255: 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);
256: /* Missing down term */
257: 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);
258: 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);
259: 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);
260: 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);
261: 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);
262: 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;
263: 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;
264: }
265: } else if (ey == N[1]-1) {
266: if (ez == 0) {
267: nEntries = 7;
268: 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);
269: 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);
270: /* Missing up term */
271: 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);
272: 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);
273: /* Missing back entry */
274: 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);
275: 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;
276: 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;
277: } else if (ez == N[2]-1) {
278: nEntries = 7;
279: 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);
280: 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);
281: /* Missing up term */
282: 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);
283: 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);
284: 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);
285: /* Missing front term */
286: 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;
287: 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;
288: } else {
289: nEntries = 8;
290: 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);
291: 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);
292: /* Missing up term */
293: 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);
294: 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);
295: 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);
296: 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);
297: 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;
298: 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;
299: }
300: } else if (ez == 0) {
301: nEntries = 8;
302: 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);
303: 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);
304: 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);
305: 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);
306: 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);
307: /* Missing back term */
308: 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);
309: 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;
310: 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;
311: } else if (ez == N[2]-1) {
312: nEntries = 8;
313: 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);
314: 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);
315: 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);
316: 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);
317: 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);
318: 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);
319: /* Missing front term */
320: 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;
321: 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;
322: } else {
323: nEntries = 9;
324: 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);
325: 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);
326: 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);
327: 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);
328: 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);
329: 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);
330: 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);
331: 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;
332: 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;
333: }
334: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
335: valRhs = fx(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
336: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
337: }
339: /* Equation on bottom face of this element */
340: if (ey == 0) {
341: /* Bottom boundary velocity Dirichlet */
342: DMStagStencil row;
343: PetscScalar valRhs;
344: const PetscScalar valA = 1.0;
345: row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
346: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
347: valRhs = uyRef(arrCoord[ez][ey][ex][icuy[0]],arrCoord[ez][ey][ex][icuy[1]],arrCoord[ez][ey][ex][icuy[2]]);
348: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
349: } else {
350: /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
351: DMStagStencil row,col[9];
352: PetscScalar valA[9],valRhs;
353: PetscInt nEntries;
355: row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
356: if (ex ==0) {
357: if (ez == 0) {
358: nEntries = 7;
359: 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);
360: 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);
361: 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);
362: /* Left term missing */
363: 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);
364: /* Back term missing */
365: 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);
366: 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;
367: 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;
368: } else if (ez == N[2]-1) {
369: nEntries = 7;
370: 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);
371: 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);
372: 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);
373: /* Left term missing */
374: 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);
375: 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);
376: /* Front term missing */
377: 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;
378: 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;
379: } else {
380: nEntries = 8;
381: 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);
382: 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);
383: 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);
384: /* Left term missing */
385: 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);
386: 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);
387: 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);
388: 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;
389: 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;
390: }
391: } else if (ex == N[0]-1) {
392: if (ez == 0) {
393: nEntries = 7;
394: 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);
395: 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);
396: 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);
397: 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);
398: /* Right term missing */
399: /* Back term missing */
400: 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);
401: 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;
402: 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;
403: } else if (ez == N[2]-1) {
404: nEntries = 7;
405: 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);
406: 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);
407: 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);
408: 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);
409: /* Right term missing */
410: 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);
411: /* Front term missing */
412: 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;
413: 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;
414: } else {
415: nEntries = 8;
416: 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);
417: 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);
418: 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);
419: 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);
420: /* Right term missing */
421: 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);
422: 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);
423: 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;
424: 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;
425: }
426: } else if (ez == 0) {
427: nEntries = 8;
428: 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);
429: 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);
430: 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);
431: 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);
432: 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);
433: /* Back term missing */
434: 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);
435: 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;
436: 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;
437: } else if (ez == N[2]-1) {
438: nEntries = 8;
439: 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);
440: 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);
441: 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);
442: 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);
443: 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);
444: 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);
445: /* Front term missing */
446: 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;
447: 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;
448: } else {
449: nEntries = 9;
450: 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);
451: 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);
452: 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);
453: 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);
454: 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);
455: 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);
456: 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);
457: 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;
458: 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;
459: }
460: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
461: valRhs = fy(arrCoord[ez][ey][ex][icuy[0]],arrCoord[ez][ey][ex][icuy[1]],arrCoord[ez][ey][ex][icuy[2]]);
462: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
463: }
465: /* Equation on back face of this element */
466: if (ez == 0) {
467: /* Back boundary velocity Dirichlet */
468: DMStagStencil row;
469: PetscScalar valRhs;
470: const PetscScalar valA = 1.0;
471: row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
472: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
473: valRhs = uzRef(arrCoord[ez][ey][ex][icuz[0]],arrCoord[ez][ey][ex][icuz[1]],arrCoord[ez][ey][ex][icuz[2]]);
474: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
475: } else {
476: /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
477: DMStagStencil row,col[9];
478: PetscScalar valA[9],valRhs;
479: PetscInt nEntries;
481: row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
482: if (ex == 0) {
483: if (ey == 0) {
484: nEntries = 7;
485: 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);
486: /* Down term missing */
487: 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);
488: /* Left term missing */
489: 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);
490: 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);
491: 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);
492: 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;
493: 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;
494: } else if (ey == N[1]-1) {
495: nEntries = 7;
496: 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);
497: 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);
498: /* Up term missing */
499: /* Left term missing */
500: 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);
501: 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);
502: 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);
503: 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;
504: 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;
505: } else {
506: nEntries = 8;
507: 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);
508: 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);
509: 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);
510: /* Left term missing */
511: 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);
512: 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);
513: 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);
514: 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;
515: 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;
516: }
517: } else if (ex == N[0]-1) {
518: if (ey == 0) {
519: nEntries = 7;
520: 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);
521: /* Down term missing */
522: 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);
523: 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);
524: /* Right term missing */
525: 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);
526: 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);
527: 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;
528: 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;
529: } else if (ey == N[1]-1) {
530: nEntries = 7;
531: 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);
532: 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);
533: /* Up term missing */
534: 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);
535: /* Right term missing */
536: 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);
537: 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);
538: 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;
539: 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;
540: } else {
541: nEntries = 8;
542: 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);
543: 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);
544: 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);
545: 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);
546: /* Right term missing */
547: 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);
548: 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);
549: 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;
550: 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;
551: }
552: } else if (ey == 0) {
553: nEntries = 8;
554: 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);
555: /* Down term missing */
556: 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);
557: 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);
558: 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);
559: 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);
560: 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);
561: 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;
562: 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;
563: } else if (ey == N[1]-1) {
564: nEntries = 8;
565: 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);
566: 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);
567: /* Up term missing */
568: 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);
569: 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);
570: 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);
571: 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);
572: 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;
573: 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;
574: } else {
575: nEntries = 9;
576: 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);
577: 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);
578: 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);
579: 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);
580: 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);
581: 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);
582: 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);
583: 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;
584: 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;
585: }
586: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
587: valRhs = fz(arrCoord[ez][ey][ex][icuz[0]],arrCoord[ez][ey][ex][icuz[1]],arrCoord[ez][ey][ex][icuz[2]]);
588: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
589: }
591: /* P equation : u_x + v_y + w_z = g
592: Note that this includes an explicit zero on the diagonal. This is only needed for
593: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
594: if (pinPressure && ex == 0 && ey == 0 && ez == 0) { /* Pin the first pressure node, if requested */
595: DMStagStencil row;
596: PetscScalar valA,valRhs;
597: row.i = ex; row.j = ey; row.k = ez; row.loc = ELEMENT; row.c = 0;
598: valA = 1.0;
599: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
600: valRhs = pRef(arrCoord[ez][ey][ex][icp[0]],arrCoord[ez][ey][ex][icp[1]],arrCoord[ez][ey][ex][icp[2]]);
601: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
602: } else {
603: DMStagStencil row,col[7];
604: PetscScalar valA[7],valRhs;
606: row.i = ex; row.j = ey; row.k = ez; row.loc = ELEMENT; row.c = 0;
607: 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;
608: 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;
609: 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;
610: 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;
611: 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;
612: 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;
613: col[6] = row; valA[6] = 0.0;
614: DMStagMatSetValuesStencil(dmSol,A,1,&row,7,col,valA,INSERT_VALUES);
615: valRhs = g(arrCoord[ez][ey][ex][icp[0]],arrCoord[ez][ey][ex][icp[1]],arrCoord[ez][ey][ex][icp[2]]);
616: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
617: }
618: }
619: }
620: }
621: DMStagVecRestoreArrayRead(dmCoord,coordLocal,&arrCoord);
622: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
623: VecAssemblyBegin(rhs);
624: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
625: VecAssemblyEnd(rhs);
627: return 0;
628: }
630: /* Create a pressure-only DMStag and use it to generate a nullspace vector
631: - Create a compatible DMStag with one dof per element (and nothing else).
632: - Create a constant vector and normalize it
633: - Migrate it to a vector on the original dmSol (making use of the fact
634: that this will fill in zeros for "extra" dof)
635: - Set the nullspace for the operator
636: - Destroy everything (the operator keeps the references it needs) */
637: static PetscErrorCode AttachNullspace(DM dmSol,Mat A)
638: {
639: DM dmPressure;
640: Vec constantPressure,basis;
641: PetscReal nrm;
642: MatNullSpace matNullSpace;
645: DMStagCreateCompatibleDMStag(dmSol,0,0,0,1,&dmPressure);
646: DMGetGlobalVector(dmPressure,&constantPressure);
647: VecSet(constantPressure,1.0);
648: VecNorm(constantPressure,NORM_2,&nrm);
649: VecScale(constantPressure,1.0/nrm);
650: DMCreateGlobalVector(dmSol,&basis);
651: DMStagMigrateVec(dmPressure,constantPressure,dmSol,basis);
652: MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol),PETSC_FALSE,1,&basis,&matNullSpace);
653: VecDestroy(&basis);
654: VecDestroy(&constantPressure);
655: MatSetNullSpace(A,matNullSpace);
656: MatNullSpaceDestroy(&matNullSpace);
657: return 0;
658: }
660: static PetscErrorCode CreateReferenceSolution(DM dmSol,Vec *pSolRef)
661: {
662: PetscInt start[3],n[3],nExtra[3],ex,ey,ez,d;
663: PetscInt ip,iux,iuy,iuz,icp[3],icux[3],icuy[3],icuz[3];
664: Vec solRef,solRefLocal,coord,coordLocal;
665: DM dmCoord;
666: PetscScalar ****arrSol,****arrCoord;
669: DMCreateGlobalVector(dmSol,pSolRef);
670: solRef = *pSolRef;
671: DMStagGetCorners(dmSol,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
672: DMGetCoordinateDM(dmSol,&dmCoord);
673: DMGetCoordinates(dmSol,&coord);
674: DMGetLocalVector(dmCoord,&coordLocal);
675: DMGlobalToLocal(dmCoord,coord,INSERT_VALUES,coordLocal);
676: DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip);
677: DMStagGetLocationSlot(dmSol,LEFT, 0,&iux);
678: DMStagGetLocationSlot(dmSol,DOWN, 0,&iuy);
679: DMStagGetLocationSlot(dmSol,BACK, 0,&iuz);
680: for (d=0; d<3; ++d) {
681: DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d]);
682: DMStagGetLocationSlot(dmCoord,LEFT, d,&icux[d]);
683: DMStagGetLocationSlot(dmCoord,DOWN, d,&icuy[d]);
684: DMStagGetLocationSlot(dmCoord,BACK, d,&icuz[d]);
685: }
686: DMStagVecGetArrayRead(dmCoord,coordLocal,&arrCoord);
687: DMGetLocalVector(dmSol,&solRefLocal);
688: DMStagVecGetArray(dmSol,solRefLocal,&arrSol);
689: for (ez=start[2]; ez<start[2] + n[2] + nExtra[2]; ++ez) {
690: for (ey=start[1]; ey<start[1] + n[1] + nExtra[1]; ++ey) {
691: for (ex=start[0]; ex<start[0] + n[0] + nExtra[0]; ++ex) {
692: if (ex < start[0] + n[0] && ey < start[1] + n[1]) {
693: arrSol[ez][ey][ex][iuz] = uzRef(
694: arrCoord[ez][ey][ex][icuz[0]],
695: arrCoord[ez][ey][ex][icuz[1]],
696: arrCoord[ez][ey][ex][icuz[2]]);
697: }
698: if (ex < start[0] + n[0] && ey < start[2] + n[2]) {
699: arrSol[ez][ey][ex][iuy] = uyRef(
700: arrCoord[ez][ey][ex][icuy[0]],
701: arrCoord[ez][ey][ex][icuy[1]],
702: arrCoord[ez][ey][ex][icuy[2]]);
703: }
704: if (ex < start[1] + n[1] && ey < start[2] + n[2]) {
705: arrSol[ez][ey][ex][iux] = uxRef(
706: arrCoord[ez][ey][ex][icux[0]],
707: arrCoord[ez][ey][ex][icux[1]],
708: arrCoord[ez][ey][ex][icux[2]]);
709: }
710: if (ex < start[0] + n[0] && ey < start[1] + n[1] && ez < start[2] + n[2]) {
711: arrSol[ez][ey][ex][ip] = pRef(
712: arrCoord[ez][ey][ex][icp[0]],
713: arrCoord[ez][ey][ex][icp[1]],
714: arrCoord[ez][ey][ex][icp[2]]);
715: }
716: }
717: }
718: }
719: DMStagVecRestoreArrayRead(dmCoord,coordLocal,&arrCoord);
720: DMStagVecRestoreArray(dmSol,solRefLocal,&arrSol);
721: DMLocalToGlobal(dmSol,solRefLocal,INSERT_VALUES,solRef);
722: DMRestoreLocalVector(dmCoord,&coordLocal);
723: DMRestoreLocalVector(dmSol,&solRefLocal);
724: return 0;
725: }
727: static PetscErrorCode CheckSolution(Vec sol,Vec solRef)
728: {
729: Vec diff;
730: PetscReal normsolRef,errAbs,errRel;
733: VecDuplicate(sol,&diff);
734: VecCopy(sol,diff);
735: VecAXPY(diff,-1.0,solRef);
736: VecNorm(diff,NORM_2,&errAbs);
737: VecNorm(solRef,NORM_2,&normsolRef);
738: errRel = errAbs/normsolRef;
739: PetscPrintf(PETSC_COMM_WORLD,"Error (abs): %g\nError (rel): %g\n",(double)errAbs,(double)errRel);
740: VecDestroy(&diff);
741: return 0;
742: }
744: /*TEST
746: test:
747: suffix: 1
748: requires: mumps
749: nsize: 27
750: 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
752: test:
753: suffix: 2
754: requires: !single
755: nsize: 4
756: 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_ksp_type gmres -fieldsplit_0_pc_gamg_esteig_ksp_max_it 5
758: test:
759: suffix: direct_umfpack
760: requires: suitesparse
761: nsize: 1
762: 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
764: TEST*/