Actual source code: ex3.c
petsc-3.11.4 2019-09-28
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: PetscBool isLastRankx,isLastRanky,isLastRankz,isFirstRankx,isFirstRanky,isFirstRankz;
152: PetscReal hx,hy,hz;
153: DM dmCoord;
154: PetscScalar ****arrCoord;
157: DMCreateMatrix(dmSol,pA);
158: A = *pA;
159: DMCreateGlobalVector(dmSol,pRhs);
160: rhs = *pRhs;
162: DMStagGetCorners(dmSol,&startx,&starty,&startz,&nx,&ny,&nz,NULL,NULL,NULL);
163: DMStagGetGlobalSizes(dmSol,&N[0],&N[1],&N[2]);
164: 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");
165: DMStagGetIsLastRank(dmSol,&isLastRankx,&isLastRanky,&isLastRankz);
166: DMStagGetIsFirstRank(dmSol,&isFirstRankx,&isFirstRanky,&isFirstRankz);
167: hx = 1.0/N[0]; hy = 1.0/N[1]; hz = 1.0/N[2];
168: DMGetCoordinateDM(dmSol,&dmCoord);
169: DMGetCoordinatesLocal(dmSol,&coordLocal);
170: DMStagVecGetArrayDOFRead(dmCoord,coordLocal,&arrCoord);
171: for (d=0; d<3; ++d) {
172: DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d] );
173: DMStagGetLocationSlot(dmCoord,LEFT, d,&icux[d] );
174: DMStagGetLocationSlot(dmCoord,DOWN, d,&icuy[d] );
175: DMStagGetLocationSlot(dmCoord,BACK, d,&icuz[d] );
176: DMStagGetLocationSlot(dmCoord,RIGHT, d,&icux_right[d]);
177: DMStagGetLocationSlot(dmCoord,UP, d,&icuy_up[d] );
178: DMStagGetLocationSlot(dmCoord,FRONT, d,&icuz_front[d]);
179: }
181: /* Loop over all local elements. Note that it may be more efficient in real
182: applications to loop over each boundary separately */
183: for (ez = startz; ez<startz+nz; ++ez) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
184: for (ey = starty; ey<starty+ny; ++ey) {
185: for (ex = startx; ex<startx+nx; ++ex) {
187: if (ex == N[0]-1) {
188: /* Right Boundary velocity Dirichlet */
189: DMStagStencil row;
190: PetscScalar valRhs;
191: const PetscScalar valA = 1.0;
192: row.i = ex; row.j = ey; row.k = ez; row.loc = RIGHT; row.c = 0;
193: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
194: valRhs = uxRef(
195: arrCoord[ez][ey][ex][icux_right[0]],
196: arrCoord[ez][ey][ex][icux_right[1]],
197: arrCoord[ez][ey][ex][icux_right[2]]
198: );
199: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
200: }
201: if (ey == N[1]-1) {
202: /* Top boundary velocity Dirichlet */
203: DMStagStencil row;
204: PetscScalar valRhs;
205: const PetscScalar valA = 1.0;
206: row.i = ex; row.j = ey; row.k = ez; row.loc = UP; row.c = 0;
207: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
208: valRhs = uyRef(
209: arrCoord[ez][ey][ex][icuy_up[0]],
210: arrCoord[ez][ey][ex][icuy_up[1]],
211: arrCoord[ez][ey][ex][icuy_up[2]]
212: );
213: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
214: }
215: if (ez == N[2]-1) {
216: /* Top boundary velocity Dirichlet */
217: DMStagStencil row;
218: PetscScalar valRhs;
219: const PetscScalar valA = 1.0;
220: row.i = ex; row.j = ey; row.k = ez; row.loc = FRONT; row.c = 0;
221: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
222: valRhs = uzRef(
223: arrCoord[ez][ey][ex][icuz_front[0]],
224: arrCoord[ez][ey][ex][icuz_front[1]],
225: arrCoord[ez][ey][ex][icuz_front[2]]
226: );
227: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
228: }
230: /* Equation on left face of this element */
231: if (ex == 0) {
232: /* Left velocity Dirichlet */
233: DMStagStencil row;
234: PetscScalar valRhs;
235: const PetscScalar valA = 1.0;
236: row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
237: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
238: valRhs = uxRef(
239: arrCoord[ez][ey][ex][icux[0]],
240: arrCoord[ez][ey][ex][icux[1]],
241: arrCoord[ez][ey][ex][icux[2]]
242: );
243: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
244: } else {
245: /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
246: DMStagStencil row,col[9];
247: PetscScalar valA[9],valRhs;
248: PetscInt nEntries;
250: row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
251: if (ey == 0) {
252: if (ez == 0) {
253: nEntries = 7;
254: 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);
255: /* Missing down term */
256: 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);
257: 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);
258: 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);
259: /* Missing back term */
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-1; col[5].j = ey ; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hx;
262: 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;
263: } else if (ez == N[2]-1) {
264: nEntries = 7;
265: 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);
266: /* Missing down term */
267: 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);
268: 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);
269: 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);
270: 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);
271: /* Missing front term */
272: 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;
273: 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;
274: } else {
275: nEntries = 8;
276: 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);
277: /* Missing down term */
278: 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);
279: 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);
280: 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);
281: 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);
282: 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);
283: 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;
284: 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;
285: }
286: } else if (ey == N[1]-1) {
287: if (ez == 0) {
288: nEntries = 7;
289: 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);
290: 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);
291: /* Missing up term */
292: 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);
293: 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);
294: /* Missing back entry */
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-1; col[5].j = ey ; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hx;
297: 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;
298: } else if (ez == N[2]-1) {
299: nEntries = 7;
300: 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);
301: 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);
302: /* Missing up term */
303: 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);
304: 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);
305: 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);
306: /* Missing front term */
307: 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;
308: 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;
309: } else {
310: nEntries = 8;
311: 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);
312: 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);
313: /* Missing up term */
314: 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);
315: 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);
316: 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);
317: 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);
318: 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;
319: 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;
320: }
321: } else if (ez == 0) {
322: nEntries = 8;
323: 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);
324: 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);
325: 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);
326: 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);
327: 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);
328: /* Missing back term */
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-1; col[6].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = 1.0 / hx;
331: 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;
332: } else if (ez == N[2]-1) {
333: nEntries = 8;
334: 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);
335: 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);
336: 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);
337: 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);
338: 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);
339: 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);
340: /* Missing front term */
341: 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;
342: 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;
343: } else {
344: nEntries = 9;
345: 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);
346: 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);
347: 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);
348: 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);
349: 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);
350: 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);
351: 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);
352: 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;
353: 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;
354: }
355: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
356: valRhs = fx(
357: arrCoord[ez][ey][ex][icux[0]],
358: arrCoord[ez][ey][ex][icux[1]],
359: arrCoord[ez][ey][ex][icux[2]]
360: );
361: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
362: }
364: /* Equation on bottom face of this element */
365: if (ey == 0) {
366: /* Bottom boundary velocity Dirichlet */
367: DMStagStencil row;
368: PetscScalar valRhs;
369: const PetscScalar valA = 1.0;
370: row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
371: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
372: valRhs = uyRef(
373: arrCoord[ez][ey][ex][icuy[0]],
374: arrCoord[ez][ey][ex][icuy[1]],
375: arrCoord[ez][ey][ex][icuy[2]]
376: );
377: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
378: } else {
379: /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
380: DMStagStencil row,col[9];
381: PetscScalar valA[9],valRhs;
382: PetscInt nEntries;
384: row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
385: if (ex ==0) {
386: if (ez == 0) {
387: nEntries = 7;
388: 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);
389: 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);
390: 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);
391: /* Left term missing */
392: 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);
393: /* Back term missing */
394: 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);
395: 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;
396: 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;
397: } else if (ez == N[2]-1) {
398: nEntries = 7;
399: 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);
400: 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);
401: 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);
402: /* Left term missing */
403: 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);
404: 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);
405: /* Front term missing */
406: col[5].i = ex ; col[4].j = ey-1; col[5].k = ez ; col[5].loc = ELEMENT; col[5].c = 0; valA[5] = 1.0 / hy;
407: col[6].i = ex ; col[5].j = ey ; col[6].k = ez ; col[6].loc = ELEMENT; col[6].c = 0; valA[6] = -1.0 / hy;
408: } else {
409: nEntries = 8;
410: 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);
411: 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);
412: 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);
413: /* Left term missing */
414: 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);
415: 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);
416: 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);
417: 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;
418: 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;
419: }
420: } else if (ex == N[0]-1) {
421: if (ez == 0) {
422: nEntries = 7;
423: 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);
424: 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);
425: 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);
426: 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);
427: /* Right term missing */
428: /* Back term missing */
429: 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);
430: 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;
431: 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;
432: } else if (ez == N[2]-1) {
433: nEntries = 7;
434: 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);
435: 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);
436: 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);
437: 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);
438: /* Right term missing */
439: 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);
440: /* Front term missing */
441: 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;
442: 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;
443: } else {
444: nEntries = 8;
445: 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);
446: 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);
447: 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);
448: 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);
449: /* Right term missing */
450: 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);
451: 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);
452: 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;
453: 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;
454: }
455: } else if (ez == 0) {
456: nEntries = 8;
457: 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);
458: 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);
459: 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);
460: 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);
461: 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);
462: /* Back term missing */
463: 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);
464: 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;
465: 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;
466: } else if (ez == N[2]-1) {
467: nEntries = 8;
468: 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);
469: 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);
470: 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);
471: 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);
472: 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);
473: 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);
474: /* Front term missing */
475: 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;
476: 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;
477: } else {
478: nEntries = 9;
479: 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);
480: 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);
481: 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);
482: 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);
483: 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);
484: 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);
485: 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);
486: 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;
487: 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;
488: }
489: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
490: valRhs = fy(
491: arrCoord[ez][ey][ex][icuy[0]],
492: arrCoord[ez][ey][ex][icuy[1]],
493: arrCoord[ez][ey][ex][icuy[2]]
494: );
495: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
496: }
498: /* Equation on back face of this element */
499: if (ez == 0) {
500: /* Back boundary velocity Dirichlet */
501: DMStagStencil row;
502: PetscScalar valRhs;
503: const PetscScalar valA = 1.0;
504: row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
505: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
506: valRhs = uzRef(
507: arrCoord[ez][ey][ex][icuz[0]],
508: arrCoord[ez][ey][ex][icuz[1]],
509: arrCoord[ez][ey][ex][icuz[2]]
510: );
511: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
512: } else {
513: /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
514: DMStagStencil row,col[9];
515: PetscScalar valA[9],valRhs;
516: PetscInt nEntries;
518: row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
519: if (ex == 0) {
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: /* Left term missing */
526: 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);
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: /* Left term missing */
537: 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);
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: /* Left term missing */
548: 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);
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 (ex == N[0]-1) {
555: if (ey == 0) {
556: nEntries = 7;
557: 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);
558: /* Down term missing */
559: 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);
560: 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);
561: /* Right term missing */
562: 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);
563: 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);
564: 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;
565: 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;
566: } else if (ey == N[1]-1) {
567: nEntries = 7;
568: 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);
569: 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);
570: /* Up term missing */
571: 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);
572: /* Right term missing */
573: 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);
574: 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);
575: 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;
576: 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;
577: } else {
578: nEntries = 8;
579: 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);
580: 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);
581: 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);
582: 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);
583: /* Right term missing */
584: 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);
585: 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);
586: 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;
587: 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;
588: }
589: } else if (ey == 0) {
590: nEntries = 8;
591: 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);
592: /* Down term missing */
593: 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);
594: 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);
595: 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);
596: 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);
597: 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);
598: 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;
599: 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;
600: } else if (ey == N[1]-1) {
601: nEntries = 8;
602: 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);
603: 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);
604: /* Up term missing */
605: 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);
606: 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);
607: 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);
608: 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);
609: 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;
610: 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;
611: } else {
612: nEntries = 9;
613: 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);
614: 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);
615: 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);
616: 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);
617: 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);
618: 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);
619: 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);
620: 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;
621: 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;
622: }
623: DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
624: valRhs = fz(
625: arrCoord[ez][ey][ex][icuz[0]],
626: arrCoord[ez][ey][ex][icuz[1]],
627: arrCoord[ez][ey][ex][icuz[2]]
628: );
629: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
630: }
632: /* P equation : u_x + v_y + w_z = g
633: Note that this includes an explicit zero on the diagonal. This is only needed for
634: direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
635: if (pinPressure && ex == 0 && ey == 0 && ez == 0) { /* Pin the first pressure node, if requested */
636: DMStagStencil row;
637: PetscScalar valA,valRhs;
638: row.i = ex; row.j = ey; row.k = ez; row.loc = ELEMENT; row.c = 0;
639: valA = 1.0;
640: DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
641: valRhs = pRef(
642: arrCoord[ez][ey][ex][icp[0]],
643: arrCoord[ez][ey][ex][icp[1]],
644: arrCoord[ez][ey][ex][icp[2]]
645: );
646: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
647: } else {
648: DMStagStencil row,col[7];
649: PetscScalar valA[7],valRhs;
651: row.i = ex; row.j = ey; row.k = ez; row.loc = ELEMENT; row.c = 0;
652: 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;
653: 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;
654: 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;
655: 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;
656: 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;
657: 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;
658: col[6] = row; valA[6] = 0.0;
659: DMStagMatSetValuesStencil(dmSol,A,1,&row,7,col,valA,INSERT_VALUES);
660: valRhs = g(
661: arrCoord[ez][ey][ex][icp[0]],
662: arrCoord[ez][ey][ex][icp[1]],
663: arrCoord[ez][ey][ex][icp[2]]
664: );
665: DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
666: }
667: }
668: }
669: }
670: DMStagVecRestoreArrayDOFRead(dmCoord,coordLocal,&arrCoord);
671: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
672: VecAssemblyBegin(rhs);
673: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
674: VecAssemblyEnd(rhs);
676: return(0);
677: }
679: /* Create a pressure-only DMStag and use it to generate a nullspace vector
680: - Create a compatible DMStag with one dof per element (and nothing else).
681: - Create a constant vector and normalize it
682: - Migrate it to a vector on the original dmSol (making use of the fact
683: that this will fill in zeros for "extra" dof)
684: - Set the nullspace for the operator
685: - Destroy everything (the operator keeps the references it needs) */
686: static PetscErrorCode AttachNullspace(DM dmSol,Mat A)
687: {
689: DM dmPressure;
690: Vec constantPressure,basis;
691: PetscReal nrm;
692: MatNullSpace matNullSpace;
695: DMStagCreateCompatibleDMStag(dmSol,0,0,0,1,&dmPressure);
696: DMGetGlobalVector(dmPressure,&constantPressure);
697: VecSet(constantPressure,1.0);
698: VecNorm(constantPressure,NORM_2,&nrm);
699: VecScale(constantPressure,1.0/nrm);
700: DMCreateGlobalVector(dmSol,&basis);
701: DMStagMigrateVec(dmPressure,constantPressure,dmSol,basis);
702: MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol),PETSC_FALSE,1,&basis,&matNullSpace);
703: VecDestroy(&basis);
704: VecDestroy(&constantPressure);
705: MatSetNullSpace(A,matNullSpace);
706: MatNullSpaceDestroy(&matNullSpace);
707: return(0);
708: }
710: static PetscErrorCode CreateReferenceSolution(DM dmSol,Vec *pSolRef)
711: {
712: PetscErrorCode ierr;
713: PetscInt start[3],n[3],nExtra[3],ex,ey,ez,d;
714: PetscInt ip,iux,iuy,iuz,icp[3],icux[3],icuy[3],icuz[3];
715: Vec solRef,solRefLocal,coord,coordLocal;
716: DM dmCoord;
717: PetscScalar ****arrSol,****arrCoord;
720: DMCreateGlobalVector(dmSol,pSolRef);
721: solRef = *pSolRef;
722: DMStagGetCorners(dmSol,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
723: DMGetCoordinateDM(dmSol,&dmCoord);
724: DMGetCoordinates(dmSol,&coord);
725: DMGetLocalVector(dmCoord,&coordLocal);
726: DMGlobalToLocal(dmCoord,coord,INSERT_VALUES,coordLocal);
727: DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip );
728: DMStagGetLocationSlot(dmSol,LEFT, 0,&iux);
729: DMStagGetLocationSlot(dmSol,DOWN, 0,&iuy);
730: DMStagGetLocationSlot(dmSol,BACK, 0,&iuz);
731: for (d=0; d<3; ++d) {
732: DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d] );
733: DMStagGetLocationSlot(dmCoord,LEFT, d,&icux[d]);
734: DMStagGetLocationSlot(dmCoord,DOWN, d,&icuy[d]);
735: DMStagGetLocationSlot(dmCoord,BACK, d,&icuz[d]);
736: }
737: DMStagVecGetArrayDOFRead(dmCoord,coordLocal,&arrCoord);
738: DMGetLocalVector(dmSol,&solRefLocal);
739: DMStagVecGetArrayDOF(dmSol,solRefLocal,&arrSol);
740: for (ez=start[2]; ez<start[2] + n[2] + nExtra[2]; ++ez) {
741: for (ey=start[1]; ey<start[1] + n[1] + nExtra[1]; ++ey) {
742: for (ex=start[0]; ex<start[0] + n[0] + nExtra[0]; ++ex) {
743: if (ex < start[0] + n[0] && ey < start[1] + n[1]) {
744: arrSol[ez][ey][ex][iuz] = uzRef(
745: arrCoord[ez][ey][ex][icuz[0]],
746: arrCoord[ez][ey][ex][icuz[1]],
747: arrCoord[ez][ey][ex][icuz[2]]);
748: }
749: if (ex < start[0] + n[0] && ey < start[2] + n[2]) {
750: arrSol[ez][ey][ex][iuy] = uyRef(
751: arrCoord[ez][ey][ex][icuy[0]],
752: arrCoord[ez][ey][ex][icuy[1]],
753: arrCoord[ez][ey][ex][icuy[2]]);
754: }
755: if (ex < start[1] + n[1] && ey < start[2] + n[2]) {
756: arrSol[ez][ey][ex][iux] = uxRef(
757: arrCoord[ez][ey][ex][icux[0]],
758: arrCoord[ez][ey][ex][icux[1]],
759: arrCoord[ez][ey][ex][icux[2]]);
760: }
761: if (ex < start[0] + n[0] && ey < start[1] + n[1] && ez < start[2] + n[2]) {
762: arrSol[ez][ey][ex][ip] = pRef(
763: arrCoord[ez][ey][ex][icp[0]],
764: arrCoord[ez][ey][ex][icp[1]],
765: arrCoord[ez][ey][ex][icp[2]]);
766: }
767: }
768: }
769: }
770: DMStagVecRestoreArrayDOFRead(dmCoord,coordLocal,&arrCoord);
771: DMStagVecRestoreArrayDOF(dmSol,solRefLocal,&arrSol);
772: DMLocalToGlobal(dmSol,solRefLocal,INSERT_VALUES,solRef);
773: DMRestoreLocalVector(dmCoord,&coordLocal);
774: DMRestoreLocalVector(dmSol,&solRefLocal);
775: return(0);
776: }
778: static PetscErrorCode CheckSolution(Vec sol,Vec solRef)
779: {
781: Vec diff;
782: PetscReal normsolRef,errAbs,errRel;
785: VecDuplicate(sol,&diff);
786: VecCopy(sol,diff);
787: VecAXPY(diff,-1.0,solRef);
788: VecNorm(diff,NORM_2,&errAbs);
789: VecNorm(solRef,NORM_2,&normsolRef);
790: errRel = errAbs/normsolRef;
791: PetscPrintf(PETSC_COMM_WORLD,"Error (abs): %g\nError (rel): %g\n",(double)errAbs,(double)errRel);
792: VecDestroy(&diff);
793: return(0);
794: }
796: /*TEST
798: test:
799: suffix: 1
800: requires: mumps
801: nsize: 27
802: 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
804: test:
805: suffix: 2
806: requires: !single
807: nsize: 4
808: 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
810: test:
811: suffix: direct_umfpack
812: requires: suitesparse
813: nsize: 1
814: 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
816: TEST*/