Actual source code: ex3.c

petsc-master 2020-07-10
Report Typos and Errors
  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(
192:               arrCoord[ez][ey][ex][icux_right[0]],
193:               arrCoord[ez][ey][ex][icux_right[1]],
194:               arrCoord[ez][ey][ex][icux_right[2]]
195:               );
196:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
197:         }
198:         if (ey == N[1]-1) {
199:           /* Top boundary velocity Dirichlet */
200:           DMStagStencil row;
201:           PetscScalar   valRhs;
202:           const PetscScalar valA = 1.0;
203:           row.i = ex; row.j = ey; row.k = ez; row.loc = UP; row.c = 0;
204:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
205:           valRhs = uyRef(
206:               arrCoord[ez][ey][ex][icuy_up[0]],
207:               arrCoord[ez][ey][ex][icuy_up[1]],
208:               arrCoord[ez][ey][ex][icuy_up[2]]
209:               );
210:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
211:         }
212:         if (ez == N[2]-1) {
213:           /* Front boundary velocity Dirichlet */
214:           DMStagStencil row;
215:           PetscScalar   valRhs;
216:           const PetscScalar valA = 1.0;
217:           row.i = ex; row.j = ey; row.k = ez; row.loc = FRONT; row.c = 0;
218:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
219:           valRhs = uzRef(
220:               arrCoord[ez][ey][ex][icuz_front[0]],
221:               arrCoord[ez][ey][ex][icuz_front[1]],
222:               arrCoord[ez][ey][ex][icuz_front[2]]
223:               );
224:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
225:         }

227:         /* Equation on left face of this element */
228:         if (ex == 0) {
229:           /* Left velocity Dirichlet */
230:           DMStagStencil row;
231:           PetscScalar   valRhs;
232:           const PetscScalar valA = 1.0;
233:           row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
234:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
235:           valRhs = uxRef(
236:               arrCoord[ez][ey][ex][icux[0]],
237:               arrCoord[ez][ey][ex][icux[1]],
238:               arrCoord[ez][ey][ex][icux[2]]
239:               );
240:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
241:         } else {
242:           /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
243:           DMStagStencil row,col[9];
244:           PetscScalar   valA[9],valRhs;
245:           PetscInt      nEntries;

247:           row.i = ex; row.j = ey; row.k = ez; row.loc = LEFT; row.c = 0;
248:           if (ey == 0) {
249:             if (ez == 0) {
250:               nEntries = 7;
251:               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);
252:               /* Missing down term */
253:               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);
254:               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);
255:               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);
256:               /* Missing back term */
257:               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);
258:               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;
259:               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;
260:             } else if (ez == N[2]-1) {
261:               nEntries = 7;
262:               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);
263:               /* Missing down term */
264:               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);
265:               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);
266:               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);
267:               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);
268:               /* Missing front term */
269:               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;
270:               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;
271:             } else {
272:               nEntries = 8;
273:               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);
274:               /* Missing down term */
275:               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);
276:               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);
277:               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);
278:               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);
279:               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);
280:               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;
281:               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;
282:             }
283:           } else if (ey == N[1]-1) {
284:             if (ez == 0) {
285:               nEntries = 7;
286:               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);
287:               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);
288:               /* Missing up term */
289:               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);
290:               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);
291:               /* Missing back entry */
292:               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);
293:               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;
294:               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;
295:             } else if (ez == N[2]-1) {
296:               nEntries = 7;
297:               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);
298:               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);
299:               /* Missing up term */
300:               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);
301:               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);
302:               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);
303:               /* Missing front term */
304:               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;
305:               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;
306:             } else {
307:               nEntries = 8;
308:               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);
309:               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);
310:               /* Missing up term */
311:               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);
312:               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);
313:               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);
314:               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);
315:               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;
316:               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;
317:             }
318:           } else if (ez == 0) {
319:             nEntries = 8;
320:             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);
321:             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);
322:             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);
323:             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);
324:             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);
325:             /* Missing back term */
326:             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);
327:             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;
328:             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;
329:           } else if (ez == N[2]-1) {
330:             nEntries = 8;
331:             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);
332:             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);
333:             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);
334:             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);
335:             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);
336:             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);
337:             /* Missing front term */
338:             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;
339:             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;
340:           } else {
341:             nEntries = 9;
342:             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);
343:             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);
344:             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);
345:             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);
346:             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);
347:             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);
348:             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);
349:             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;
350:             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;
351:           }
352:           DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
353:           valRhs = fx(
354:               arrCoord[ez][ey][ex][icux[0]],
355:               arrCoord[ez][ey][ex][icux[1]],
356:               arrCoord[ez][ey][ex][icux[2]]
357:               );
358:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
359:         }

361:         /* Equation on bottom face of this element */
362:         if (ey == 0) {
363:           /* Bottom boundary velocity Dirichlet */
364:           DMStagStencil row;
365:           PetscScalar   valRhs;
366:           const PetscScalar valA = 1.0;
367:           row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
368:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
369:           valRhs = uyRef(
370:               arrCoord[ez][ey][ex][icuy[0]],
371:               arrCoord[ez][ey][ex][icuy[1]],
372:               arrCoord[ez][ey][ex][icuy[2]]
373:               );
374:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
375:         } else {
376:           /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
377:           DMStagStencil row,col[9];
378:           PetscScalar   valA[9],valRhs;
379:           PetscInt      nEntries;

381:           row.i = ex; row.j = ey; row.k = ez; row.loc = DOWN; row.c = 0;
382:           if (ex ==0) {
383:             if (ez == 0) {
384:               nEntries = 7;
385:               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);
386:               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);
387:               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);
388:               /* Left term missing */
389:               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);
390:               /* Back term missing */
391:               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);
392:               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;
393:               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;
394:             } else if (ez == N[2]-1) {
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:               /* Left term missing */
400:               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);
401:               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);
402:               /* Front term missing */
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 {
406:               nEntries = 8;
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) -2.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:               /* Left term missing */
411:               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);
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:               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);
414:               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;
415:               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;
416:             }
417:           } else if (ex == N[0]-1) {
418:             if (ez == 0) {
419:               nEntries = 7;
420:               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);
421:               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);
422:               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);
423:               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);
424:               /* Right term missing */
425:               /* Back term missing */
426:               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);
427:               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;
428:               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;
429:             } else if (ez == N[2]-1) {
430:               nEntries = 7;
431:               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);
432:               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);
433:               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);
434:               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);
435:               /* Right term missing */
436:               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);
437:               /* Front term missing */
438:               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;
439:               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;
440:             } else {
441:               nEntries = 8;
442:               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);
443:               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);
444:               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);
445:               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);
446:               /* Right term missing */
447:               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);
448:               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);
449:               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;
450:               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;
451:             }
452:           } else if (ez == 0) {
453:             nEntries = 8;
454:             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);
455:             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);
456:             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);
457:             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);
458:             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);
459:             /* Back term missing */
460:             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);
461:             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;
462:             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;
463:           } else if (ez == N[2]-1) {
464:             nEntries = 8;
465:             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);
466:             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);
467:             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);
468:             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);
469:             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);
470:             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);
471:             /* Front term missing */
472:             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;
473:             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;
474:           } else {
475:             nEntries = 9;
476:             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);
477:             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);
478:             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);
479:             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);
480:             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);
481:             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);
482:             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);
483:             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;
484:             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;
485:           }
486:           DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
487:           valRhs = fy(
488:               arrCoord[ez][ey][ex][icuy[0]],
489:               arrCoord[ez][ey][ex][icuy[1]],
490:               arrCoord[ez][ey][ex][icuy[2]]
491:               );
492:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
493:         }

495:         /* Equation on back face of this element */
496:         if (ez == 0) {
497:           /* Back boundary velocity Dirichlet */
498:           DMStagStencil row;
499:           PetscScalar   valRhs;
500:           const PetscScalar valA = 1.0;
501:           row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
502:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
503:           valRhs = uzRef(
504:               arrCoord[ez][ey][ex][icuz[0]],
505:               arrCoord[ez][ey][ex][icuz[1]],
506:               arrCoord[ez][ey][ex][icuz[2]]
507:               );
508:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
509:         } else {
510:           /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
511:           DMStagStencil row,col[9];
512:           PetscScalar   valA[9],valRhs;
513:           PetscInt      nEntries;

515:           row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
516:           if (ex == 0) {
517:             if (ey == 0) {
518:               nEntries = 7;
519:               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);
520:               /* Down term missing */
521:               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);
522:               /* Left term missing */
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:               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);
525:               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);
526:               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;
527:               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;
528:             } else if (ey == N[1]-1) {
529:               nEntries = 7;
530:               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);
531:               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);
532:               /* Up term missing */
533:               /* Left 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:               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);
536:               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);
537:               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;
538:               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;
539:             } else {
540:               nEntries = 8;
541:               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);
542:               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);
543:               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);
544:               /* Left term missing */
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:               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);
547:               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);
548:               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;
549:               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;
550:             }
551:           } else if (ex == N[0]-1) {
552:             if (ey == 0) {
553:               nEntries = 7;
554:               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);
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:               /* Right term missing */
559:               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);
560:               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);
561:               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;
562:               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;
563:             } else if (ey == N[1]-1) {
564:               nEntries = 7;
565:               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);
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:               /* Right term missing */
570:               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);
571:               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);
572:               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;
573:               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;
574:             } else {
575:               nEntries = 8;
576:               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);
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:               /* Right term missing */
581:               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);
582:               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);
583:               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;
584:               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;
585:             }
586:           } else if (ey == 0) {
587:             nEntries = 8;
588:             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);
589:             /* Down term missing */
590:             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);
591:             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);
592:             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);
593:             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);
594:             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);
595:             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;
596:             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;
597:           } else if (ey == N[1]-1) {
598:             nEntries = 8;
599:             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);
600:             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);
601:             /* Up term missing */
602:             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);
603:             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);
604:             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);
605:             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);
606:             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;
607:             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;
608:           } else {
609:             nEntries = 9;
610:             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);
611:             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);
612:             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);
613:             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);
614:             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);
615:             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);
616:             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);
617:             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;
618:             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;
619:           }
620:           DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
621:           valRhs = fz(
622:               arrCoord[ez][ey][ex][icuz[0]],
623:               arrCoord[ez][ey][ex][icuz[1]],
624:               arrCoord[ez][ey][ex][icuz[2]]
625:               );
626:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
627:         }

629:         /* P equation : u_x + v_y + w_z = g
630:            Note that this includes an explicit zero on the diagonal. This is only needed for
631:            direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
632:         if (pinPressure && ex == 0 && ey == 0 && ez == 0) { /* Pin the first pressure node, if requested */
633:           DMStagStencil row;
634:           PetscScalar valA,valRhs;
635:           row.i = ex; row.j = ey; row.k = ez; row.loc  = ELEMENT; row.c = 0;
636:           valA = 1.0;
637:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
638:           valRhs = pRef(
639:               arrCoord[ez][ey][ex][icp[0]],
640:               arrCoord[ez][ey][ex][icp[1]],
641:               arrCoord[ez][ey][ex][icp[2]]
642:               );
643:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
644:         } else {
645:           DMStagStencil row,col[7];
646:           PetscScalar   valA[7],valRhs;

648:           row.i    = ex; row.j    = ey; row.k    = ez; row.loc    = ELEMENT; row.c    = 0;
649:           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;
650:           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;
651:           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;
652:           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;
653:           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;
654:           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;
655:           col[6]   = row;                                                                  valA[6] =  0.0;
656:           DMStagMatSetValuesStencil(dmSol,A,1,&row,7,col,valA,INSERT_VALUES);
657:           valRhs = g(
658:               arrCoord[ez][ey][ex][icp[0]],
659:               arrCoord[ez][ey][ex][icp[1]],
660:               arrCoord[ez][ey][ex][icp[2]]
661:               );
662:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
663:         }
664:       }
665:     }
666:   }
667:   DMStagVecRestoreArrayRead(dmCoord,coordLocal,&arrCoord);
668:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
669:   VecAssemblyBegin(rhs);
670:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
671:   VecAssemblyEnd(rhs);

673:   return(0);
674: }

676: /* Create a pressure-only DMStag and use it to generate a nullspace vector
677:    - Create a compatible DMStag with one dof per element (and nothing else).
678:    - Create a constant vector and normalize it
679:    - Migrate it to a vector on the original dmSol (making use of the fact
680:    that this will fill in zeros for "extra" dof)
681:    - Set the nullspace for the operator
682:    - Destroy everything (the operator keeps the references it needs) */
683: static PetscErrorCode AttachNullspace(DM dmSol,Mat A)
684: {
686:   DM             dmPressure;
687:   Vec            constantPressure,basis;
688:   PetscReal      nrm;
689:   MatNullSpace   matNullSpace;

692:   DMStagCreateCompatibleDMStag(dmSol,0,0,0,1,&dmPressure);
693:   DMGetGlobalVector(dmPressure,&constantPressure);
694:   VecSet(constantPressure,1.0);
695:   VecNorm(constantPressure,NORM_2,&nrm);
696:   VecScale(constantPressure,1.0/nrm);
697:   DMCreateGlobalVector(dmSol,&basis);
698:   DMStagMigrateVec(dmPressure,constantPressure,dmSol,basis);
699:   MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol),PETSC_FALSE,1,&basis,&matNullSpace);
700:   VecDestroy(&basis);
701:   VecDestroy(&constantPressure);
702:   MatSetNullSpace(A,matNullSpace);
703:   MatNullSpaceDestroy(&matNullSpace);
704:   return(0);
705: }

707: static PetscErrorCode CreateReferenceSolution(DM dmSol,Vec *pSolRef)
708: {
709:   PetscErrorCode    ierr;
710:   PetscInt          start[3],n[3],nExtra[3],ex,ey,ez,d;
711:   PetscInt          ip,iux,iuy,iuz,icp[3],icux[3],icuy[3],icuz[3];
712:   Vec               solRef,solRefLocal,coord,coordLocal;
713:   DM                dmCoord;
714:   PetscScalar       ****arrSol,****arrCoord;

717:   DMCreateGlobalVector(dmSol,pSolRef);
718:   solRef = *pSolRef;
719:   DMStagGetCorners(dmSol,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);
720:   DMGetCoordinateDM(dmSol,&dmCoord);
721:   DMGetCoordinates(dmSol,&coord);
722:   DMGetLocalVector(dmCoord,&coordLocal);
723:   DMGlobalToLocal(dmCoord,coord,INSERT_VALUES,coordLocal);
724:   DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip );
725:   DMStagGetLocationSlot(dmSol,LEFT,   0,&iux);
726:   DMStagGetLocationSlot(dmSol,DOWN,   0,&iuy);
727:   DMStagGetLocationSlot(dmSol,BACK,   0,&iuz);
728:   for (d=0; d<3; ++d) {
729:     DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d] );
730:     DMStagGetLocationSlot(dmCoord,LEFT,   d,&icux[d]);
731:     DMStagGetLocationSlot(dmCoord,DOWN,   d,&icuy[d]);
732:     DMStagGetLocationSlot(dmCoord,BACK,   d,&icuz[d]);
733:   }
734:   DMStagVecGetArrayRead(dmCoord,coordLocal,&arrCoord);
735:   DMGetLocalVector(dmSol,&solRefLocal);
736:   DMStagVecGetArray(dmSol,solRefLocal,&arrSol);
737:   for (ez=start[2]; ez<start[2] + n[2] + nExtra[2]; ++ez) {
738:     for (ey=start[1]; ey<start[1] + n[1] + nExtra[1]; ++ey) {
739:       for (ex=start[0]; ex<start[0] + n[0] + nExtra[0]; ++ex) {
740:         if (ex < start[0] + n[0] && ey < start[1] + n[1]) {
741:           arrSol[ez][ey][ex][iuz] = uzRef(
742:               arrCoord[ez][ey][ex][icuz[0]],
743:               arrCoord[ez][ey][ex][icuz[1]],
744:               arrCoord[ez][ey][ex][icuz[2]]);
745:         }
746:         if (ex < start[0] + n[0] && ey < start[2] + n[2]) {
747:           arrSol[ez][ey][ex][iuy] = uyRef(
748:               arrCoord[ez][ey][ex][icuy[0]],
749:               arrCoord[ez][ey][ex][icuy[1]],
750:               arrCoord[ez][ey][ex][icuy[2]]);
751:         }
752:         if (ex < start[1] + n[1] && ey < start[2] + n[2]) {
753:           arrSol[ez][ey][ex][iux] = uxRef(
754:               arrCoord[ez][ey][ex][icux[0]],
755:               arrCoord[ez][ey][ex][icux[1]],
756:               arrCoord[ez][ey][ex][icux[2]]);
757:         }
758:         if (ex < start[0] + n[0] && ey < start[1] + n[1] && ez < start[2] + n[2]) {
759:           arrSol[ez][ey][ex][ip]  = pRef(
760:               arrCoord[ez][ey][ex][icp[0]],
761:               arrCoord[ez][ey][ex][icp[1]],
762:               arrCoord[ez][ey][ex][icp[2]]);
763:         }
764:       }
765:     }
766:   }
767:   DMStagVecRestoreArrayRead(dmCoord,coordLocal,&arrCoord);
768:   DMStagVecRestoreArray(dmSol,solRefLocal,&arrSol);
769:   DMLocalToGlobal(dmSol,solRefLocal,INSERT_VALUES,solRef);
770:   DMRestoreLocalVector(dmCoord,&coordLocal);
771:   DMRestoreLocalVector(dmSol,&solRefLocal);
772:   return(0);
773: }

775: static PetscErrorCode CheckSolution(Vec sol,Vec solRef)
776: {
778:   Vec            diff;
779:   PetscReal      normsolRef,errAbs,errRel;

782:   VecDuplicate(sol,&diff);
783:   VecCopy(sol,diff);
784:   VecAXPY(diff,-1.0,solRef);
785:   VecNorm(diff,NORM_2,&errAbs);
786:   VecNorm(solRef,NORM_2,&normsolRef);
787:   errRel = errAbs/normsolRef;
788:   PetscPrintf(PETSC_COMM_WORLD,"Error (abs): %g\nError (rel): %g\n",(double)errAbs,(double)errRel);
789:   VecDestroy(&diff);
790:   return(0);
791: }

793: /*TEST

795:    test:
796:       suffix: 1
797:       requires: mumps
798:       nsize: 27
799:       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

801:    test:
802:       suffix: 2
803:       requires: !single
804:       nsize: 4
805:       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

807:    test:
808:       suffix: direct_umfpack
809:       requires: suitesparse
810:       nsize: 1
811:       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

813: TEST*/