Actual source code: ex3.c

  1: static char help[] = "Solve a toy 3D problem on a staggered grid\n\n";
  2: /*

  4:   To demonstrate the basic functionality of DMStag, solves an isoviscous
  5:   incompressible Stokes problem on a rectangular 3D domain.

  7:   u_{xx} + u_{yy} + u_{zz} - p_x = f^x
  8:   v_{xx} + v_{yy} + u_{zz} - p_y = f^y
  9:   w_{xx} + w_{yy} + w_{zz} - p_y = f^z
 10:   u_x    + v_y    + w_z          = g

 12:   g = 0 for the physical case.

 14:   Boundary conditions give prescribed flow perpendicular to the boundaries,
 15:   and zero derivative perpendicular to them (free slip). This involves
 16:   using a modifed stencil at the boundaries. Another option would be to
 17:   use DM_BOUNDARY_GHOSTED in DMStagCreate3d() and a matrix-free operator (MATSHELL)
 18:   making use of the uniformly-available ghost layer.

 20:   Use the -pinpressure option to fix a pressure node, instead of providing
 21:   a constant-pressure nullspace. This allows use of direct solvers, e.g. to
 22:   use UMFPACK,

 24:      ./ex3 -pinpressure 1 -pc_type lu -pc_factor_mat_solver_type umfpack

 26: */
 27: #include <petscdm.h>
 28: #include <petscksp.h>
 29: #include <petscdmstag.h>

 31: /* Shorter, more convenient names for DMStagStencilLocation entries */
 32: #define BACK_DOWN_LEFT   DMSTAG_BACK_DOWN_LEFT
 33: #define BACK_DOWN        DMSTAG_BACK_DOWN
 34: #define BACK_DOWN_RIGHT  DMSTAG_BACK_DOWN_RIGHT
 35: #define BACK_LEFT        DMSTAG_BACK_LEFT
 36: #define BACK             DMSTAG_BACK
 37: #define BACK_RIGHT       DMSTAG_BACK_RIGHT
 38: #define BACK_UP_LEFT     DMSTAG_BACK_UP_LEFT
 39: #define BACK_UP          DMSTAG_BACK_UP
 40: #define BACK_UP_RIGHT    DMSTAG_BACK_UP_RIGHT
 41: #define DOWN_LEFT        DMSTAG_DOWN_LEFT
 42: #define DOWN             DMSTAG_DOWN
 43: #define DOWN_RIGHT       DMSTAG_DOWN_RIGHT
 44: #define LEFT             DMSTAG_LEFT
 45: #define ELEMENT          DMSTAG_ELEMENT
 46: #define RIGHT            DMSTAG_RIGHT
 47: #define UP_LEFT          DMSTAG_UP_LEFT
 48: #define UP               DMSTAG_UP
 49: #define UP_RIGHT         DMSTAG_UP_RIGHT
 50: #define FRONT_DOWN_LEFT  DMSTAG_FRONT_DOWN_LEFT
 51: #define FRONT_DOWN       DMSTAG_FRONT_DOWN
 52: #define FRONT_DOWN_RIGHT DMSTAG_FRONT_DOWN_RIGHT
 53: #define FRONT_LEFT       DMSTAG_FRONT_LEFT
 54: #define FRONT            DMSTAG_FRONT
 55: #define FRONT_RIGHT      DMSTAG_FRONT_RIGHT
 56: #define FRONT_UP_LEFT    DMSTAG_FRONT_UP_LEFT
 57: #define FRONT_UP         DMSTAG_FRONT_UP
 58: #define FRONT_UP_RIGHT   DMSTAG_FRONT_UP_RIGHT

 60: static PetscErrorCode CreateReferenceSolution(DM, Vec *);
 61: static PetscErrorCode CreateSystem(DM, Mat *, Vec *, PetscBool);
 62: static PetscErrorCode AttachNullspace(DM, Mat);
 63: static PetscErrorCode CheckSolution(Vec, Vec);

 65: /* Manufactured solution. Chosen to be higher order than can be solved exactly,
 66: and to have a zero derivative for flow parallel to the boundaries. That is,
 67: d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
 68: and left boundaries.
 69: These expressions could be made more interesting with more z terms,
 70: and convergence could be confirmed.  */

 72: static PetscScalar uxRef(PetscScalar x, PetscScalar y, PetscScalar z)
 73: {
 74:   return 0.0 * x + y * y - 2.0 * y * y * y + y * y * y * y + 0.0 * z;
 75: }
 76: static PetscScalar uyRef(PetscScalar x, PetscScalar y, PetscScalar z)
 77: {
 78:   return x * x - 2.0 * x * x * x + x * x * x * x + 0.0 * y + 0.0 * z;
 79: }
 80: static PetscScalar uzRef(PetscScalar x, PetscScalar y, PetscScalar z)
 81: {
 82:   return 0.0 * x + 0.0 * y + 0.0 * z + 1.0;
 83: }
 84: static PetscScalar pRef(PetscScalar x, PetscScalar y, PetscScalar z)
 85: {
 86:   return -1.0 * (x - 0.5) + -3.0 / 2.0 * y * y + 0.5 - 2.0 * (z - 1.0);
 87: } /* zero integral */
 88: static PetscScalar fx(PetscScalar x, PetscScalar y, PetscScalar z)
 89: {
 90:   return 0.0 * x + 2.0 - 12.0 * y + 12.0 * y * y + 0.0 * z + 1.0;
 91: }
 92: static PetscScalar fy(PetscScalar x, PetscScalar y, PetscScalar z)
 93: {
 94:   return 2.0 - 12.0 * x + 12.0 * x * x + 3.0 * y + 0.0 * z;
 95: }
 96: static PetscScalar fz(PetscScalar x, PetscScalar y, PetscScalar z)
 97: {
 98:   return 0.0 * x + 0.0 * y + 0.0 * z + 2.0;
 99: }
100: static PetscScalar g(PetscScalar x, PetscScalar y, PetscScalar z)
101: {
102:   return 0.0 * x + 0.0 * y + 0.0 * z + 0.0;
103: }

105: int main(int argc, char **argv)
106: {
107:   DM        dmSol;
108:   Vec       sol, solRef, rhs;
109:   Mat       A;
110:   KSP       ksp;
111:   PC        pc;
112:   PetscBool pinPressure;

114:   /* Initialize PETSc and process command line arguments */
116:   PetscInitialize(&argc, &argv, (char *)0, help);
117:   pinPressure = PETSC_TRUE;
118:   PetscOptionsGetBool(NULL, NULL, "-pinpressure", &pinPressure, NULL);

120:   /* Create 3D DMStag for the solution, and set up. */
121:   {
122:     const PetscInt dof0 = 0, dof1 = 0, dof2 = 1, dof3 = 1; /* 1 dof on each face and element center */
123:     const PetscInt stencilWidth = 1;
124:     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);
125:     DMSetFromOptions(dmSol);
126:     DMSetUp(dmSol);
127:     DMStagSetUniformCoordinatesExplicit(dmSol, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
128:     /* Note: also see ex2.c, where another, more efficient option is demonstrated,
129:        using DMStagSetUniformCoordinatesProduct() */
130:   }

132:   /* Compute (manufactured) reference solution */
133:   CreateReferenceSolution(dmSol, &solRef);

135:   /* Assemble system */
136:   CreateSystem(dmSol, &A, &rhs, pinPressure);

138:   /* Attach a constant-pressure nullspace to the operator */
139:   if (!pinPressure) AttachNullspace(dmSol, A);

141:   /* Solve */
142:   DMCreateGlobalVector(dmSol, &sol);
143:   KSPCreate(PETSC_COMM_WORLD, &ksp);
144:   KSPSetType(ksp, KSPFGMRES);
145:   KSPSetOperators(ksp, A, A);
146:   KSPGetPC(ksp, &pc);
147:   PCSetType(pc, PCFIELDSPLIT);
148:   PCFieldSplitSetDetectSaddlePoint(pc, PETSC_TRUE);
149:   KSPSetFromOptions(ksp);
150:   KSPSolve(ksp, rhs, sol);

152:   /* Check Solution */
153:   CheckSolution(sol, solRef);

155:   /* Clean up and finalize PETSc */
156:   KSPDestroy(&ksp);
157:   VecDestroy(&sol);
158:   VecDestroy(&solRef);
159:   VecDestroy(&rhs);
160:   MatDestroy(&A);
161:   DMDestroy(&dmSol);
162:   PetscFinalize();
163:   return 0;
164: }

166: static PetscErrorCode CreateSystem(DM dmSol, Mat *pA, Vec *pRhs, PetscBool pinPressure)
167: {
168:   Vec             rhs, coordLocal;
169:   Mat             A;
170:   PetscInt        startx, starty, startz, N[3], nx, ny, nz, ex, ey, ez, d;
171:   PetscInt        icp[3], icux[3], icuy[3], icuz[3], icux_right[3], icuy_up[3], icuz_front[3];
172:   PetscReal       hx, hy, hz;
173:   DM              dmCoord;
174:   PetscScalar ****arrCoord;

177:   DMCreateMatrix(dmSol, pA);
178:   A = *pA;
179:   DMCreateGlobalVector(dmSol, pRhs);
180:   rhs = *pRhs;

182:   DMStagGetCorners(dmSol, &startx, &starty, &startz, &nx, &ny, &nz, NULL, NULL, NULL);
183:   DMStagGetGlobalSizes(dmSol, &N[0], &N[1], &N[2]);
185:   hx = 1.0 / N[0];
186:   hy = 1.0 / N[1];
187:   hz = 1.0 / N[2];
188:   DMGetCoordinateDM(dmSol, &dmCoord);
189:   DMGetCoordinatesLocal(dmSol, &coordLocal);
190:   DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord);
191:   for (d = 0; d < 3; ++d) {
192:     DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]);
193:     DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]);
194:     DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]);
195:     DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]);
196:     DMStagGetLocationSlot(dmCoord, RIGHT, d, &icux_right[d]);
197:     DMStagGetLocationSlot(dmCoord, UP, d, &icuy_up[d]);
198:     DMStagGetLocationSlot(dmCoord, FRONT, d, &icuz_front[d]);
199:   }

201:   /* Loop over all local elements. Note that it may be more efficient in real
202:      applications to loop over each boundary separately */
203:   for (ez = startz; ez < startz + nz; ++ez) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
204:     for (ey = starty; ey < starty + ny; ++ey) {
205:       for (ex = startx; ex < startx + nx; ++ex) {
206:         if (ex == N[0] - 1) {
207:           /* Right Boundary velocity Dirichlet */
208:           DMStagStencil     row;
209:           PetscScalar       valRhs;
210:           const PetscScalar valA = 1.0;
211:           row.i                  = ex;
212:           row.j                  = ey;
213:           row.k                  = ez;
214:           row.loc                = RIGHT;
215:           row.c                  = 0;
216:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
217:           valRhs = uxRef(arrCoord[ez][ey][ex][icux_right[0]], arrCoord[ez][ey][ex][icux_right[1]], arrCoord[ez][ey][ex][icux_right[2]]);
218:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
219:         }
220:         if (ey == N[1] - 1) {
221:           /* Top boundary velocity Dirichlet */
222:           DMStagStencil     row;
223:           PetscScalar       valRhs;
224:           const PetscScalar valA = 1.0;
225:           row.i                  = ex;
226:           row.j                  = ey;
227:           row.k                  = ez;
228:           row.loc                = UP;
229:           row.c                  = 0;
230:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
231:           valRhs = uyRef(arrCoord[ez][ey][ex][icuy_up[0]], arrCoord[ez][ey][ex][icuy_up[1]], arrCoord[ez][ey][ex][icuy_up[2]]);
232:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
233:         }
234:         if (ez == N[2] - 1) {
235:           /* Front boundary velocity Dirichlet */
236:           DMStagStencil     row;
237:           PetscScalar       valRhs;
238:           const PetscScalar valA = 1.0;
239:           row.i                  = ex;
240:           row.j                  = ey;
241:           row.k                  = ez;
242:           row.loc                = FRONT;
243:           row.c                  = 0;
244:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
245:           valRhs = uzRef(arrCoord[ez][ey][ex][icuz_front[0]], arrCoord[ez][ey][ex][icuz_front[1]], arrCoord[ez][ey][ex][icuz_front[2]]);
246:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
247:         }

249:         /* Equation on left face of this element */
250:         if (ex == 0) {
251:           /* Left velocity Dirichlet */
252:           DMStagStencil     row;
253:           PetscScalar       valRhs;
254:           const PetscScalar valA = 1.0;
255:           row.i                  = ex;
256:           row.j                  = ey;
257:           row.k                  = ez;
258:           row.loc                = LEFT;
259:           row.c                  = 0;
260:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
261:           valRhs = uxRef(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
262:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
263:         } else {
264:           /* X-momentum interior equation : (u_xx + u_yy + u_zz) - p_x = f^x */
265:           DMStagStencil row, col[9];
266:           PetscScalar   valA[9], valRhs;
267:           PetscInt      nEntries;

269:           row.i   = ex;
270:           row.j   = ey;
271:           row.k   = ez;
272:           row.loc = LEFT;
273:           row.c   = 0;
274:           if (ey == 0) {
275:             if (ez == 0) {
276:               nEntries   = 7;
277:               col[0].i   = ex;
278:               col[0].j   = ey;
279:               col[0].k   = ez;
280:               col[0].loc = LEFT;
281:               col[0].c   = 0;
282:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
283:               /* Missing down term */
284:               col[1].i   = ex;
285:               col[1].j   = ey + 1;
286:               col[1].k   = ez;
287:               col[1].loc = LEFT;
288:               col[1].c   = 0;
289:               valA[1]    = 1.0 / (hy * hy);
290:               col[2].i   = ex - 1;
291:               col[2].j   = ey;
292:               col[2].k   = ez;
293:               col[2].loc = LEFT;
294:               col[2].c   = 0;
295:               valA[2]    = 1.0 / (hx * hx);
296:               col[3].i   = ex + 1;
297:               col[3].j   = ey;
298:               col[3].k   = ez;
299:               col[3].loc = LEFT;
300:               col[3].c   = 0;
301:               valA[3]    = 1.0 / (hx * hx);
302:               /* Missing back term */
303:               col[4].i   = ex;
304:               col[4].j   = ey;
305:               col[4].k   = ez + 1;
306:               col[4].loc = LEFT;
307:               col[4].c   = 0;
308:               valA[4]    = 1.0 / (hz * hz);
309:               col[5].i   = ex - 1;
310:               col[5].j   = ey;
311:               col[5].k   = ez;
312:               col[5].loc = ELEMENT;
313:               col[5].c   = 0;
314:               valA[5]    = 1.0 / hx;
315:               col[6].i   = ex;
316:               col[6].j   = ey;
317:               col[6].k   = ez;
318:               col[6].loc = ELEMENT;
319:               col[6].c   = 0;
320:               valA[6]    = -1.0 / hx;
321:             } else if (ez == N[2] - 1) {
322:               nEntries   = 7;
323:               col[0].i   = ex;
324:               col[0].j   = ey;
325:               col[0].k   = ez;
326:               col[0].loc = LEFT;
327:               col[0].c   = 0;
328:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 1.0 / (hz * hz);
329:               /* Missing down term */
330:               col[1].i   = ex;
331:               col[1].j   = ey + 1;
332:               col[1].k   = ez;
333:               col[1].loc = LEFT;
334:               col[1].c   = 0;
335:               valA[1]    = 1.0 / (hy * hy);
336:               col[2].i   = ex - 1;
337:               col[2].j   = ey;
338:               col[2].k   = ez;
339:               col[2].loc = LEFT;
340:               col[2].c   = 0;
341:               valA[2]    = 1.0 / (hx * hx);
342:               col[3].i   = ex + 1;
343:               col[3].j   = ey;
344:               col[3].k   = ez;
345:               col[3].loc = LEFT;
346:               col[3].c   = 0;
347:               valA[3]    = 1.0 / (hx * hx);
348:               col[4].i   = ex;
349:               col[4].j   = ey;
350:               col[4].k   = ez - 1;
351:               col[4].loc = LEFT;
352:               col[4].c   = 0;
353:               valA[4]    = 1.0 / (hz * hz);
354:               /* Missing front term */
355:               col[5].i   = ex - 1;
356:               col[5].j   = ey;
357:               col[5].k   = ez;
358:               col[5].loc = ELEMENT;
359:               col[5].c   = 0;
360:               valA[5]    = 1.0 / hx;
361:               col[6].i   = ex;
362:               col[6].j   = ey;
363:               col[6].k   = ez;
364:               col[6].loc = ELEMENT;
365:               col[6].c   = 0;
366:               valA[6]    = -1.0 / hx;
367:             } else {
368:               nEntries   = 8;
369:               col[0].i   = ex;
370:               col[0].j   = ey;
371:               col[0].k   = ez;
372:               col[0].loc = LEFT;
373:               col[0].c   = 0;
374:               valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
375:               /* Missing down term */
376:               col[1].i   = ex;
377:               col[1].j   = ey + 1;
378:               col[1].k   = ez;
379:               col[1].loc = LEFT;
380:               col[1].c   = 0;
381:               valA[1]    = 1.0 / (hy * hy);
382:               col[2].i   = ex - 1;
383:               col[2].j   = ey;
384:               col[2].k   = ez;
385:               col[2].loc = LEFT;
386:               col[2].c   = 0;
387:               valA[2]    = 1.0 / (hx * hx);
388:               col[3].i   = ex + 1;
389:               col[3].j   = ey;
390:               col[3].k   = ez;
391:               col[3].loc = LEFT;
392:               col[3].c   = 0;
393:               valA[3]    = 1.0 / (hx * hx);
394:               col[4].i   = ex;
395:               col[4].j   = ey;
396:               col[4].k   = ez - 1;
397:               col[4].loc = LEFT;
398:               col[4].c   = 0;
399:               valA[4]    = 1.0 / (hz * hz);
400:               col[5].i   = ex;
401:               col[5].j   = ey;
402:               col[5].k   = ez + 1;
403:               col[5].loc = LEFT;
404:               col[5].c   = 0;
405:               valA[5]    = 1.0 / (hz * hz);
406:               col[6].i   = ex - 1;
407:               col[6].j   = ey;
408:               col[6].k   = ez;
409:               col[6].loc = ELEMENT;
410:               col[6].c   = 0;
411:               valA[6]    = 1.0 / hx;
412:               col[7].i   = ex;
413:               col[7].j   = ey;
414:               col[7].k   = ez;
415:               col[7].loc = ELEMENT;
416:               col[7].c   = 0;
417:               valA[7]    = -1.0 / hx;
418:             }
419:           } else if (ey == N[1] - 1) {
420:             if (ez == 0) {
421:               nEntries   = 7;
422:               col[0].i   = ex;
423:               col[0].j   = ey;
424:               col[0].k   = ez;
425:               col[0].loc = LEFT;
426:               col[0].c   = 0;
427:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
428:               col[1].i   = ex;
429:               col[1].j   = ey - 1;
430:               col[1].k   = ez;
431:               col[1].loc = LEFT;
432:               col[1].c   = 0;
433:               valA[1]    = 1.0 / (hy * hy);
434:               /* Missing up term */
435:               col[2].i   = ex - 1;
436:               col[2].j   = ey;
437:               col[2].k   = ez;
438:               col[2].loc = LEFT;
439:               col[2].c   = 0;
440:               valA[2]    = 1.0 / (hx * hx);
441:               col[3].i   = ex + 1;
442:               col[3].j   = ey;
443:               col[3].k   = ez;
444:               col[3].loc = LEFT;
445:               col[3].c   = 0;
446:               valA[3]    = 1.0 / (hx * hx);
447:               /* Missing back entry */
448:               col[4].i   = ex;
449:               col[4].j   = ey;
450:               col[4].k   = ez + 1;
451:               col[4].loc = LEFT;
452:               col[4].c   = 0;
453:               valA[4]    = 1.0 / (hz * hz);
454:               col[5].i   = ex - 1;
455:               col[5].j   = ey;
456:               col[5].k   = ez;
457:               col[5].loc = ELEMENT;
458:               col[5].c   = 0;
459:               valA[5]    = 1.0 / hx;
460:               col[6].i   = ex;
461:               col[6].j   = ey;
462:               col[6].k   = ez;
463:               col[6].loc = ELEMENT;
464:               col[6].c   = 0;
465:               valA[6]    = -1.0 / hx;
466:             } else if (ez == N[2] - 1) {
467:               nEntries   = 7;
468:               col[0].i   = ex;
469:               col[0].j   = ey;
470:               col[0].k   = ez;
471:               col[0].loc = LEFT;
472:               col[0].c   = 0;
473:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
474:               col[1].i   = ex;
475:               col[1].j   = ey - 1;
476:               col[1].k   = ez;
477:               col[1].loc = LEFT;
478:               col[1].c   = 0;
479:               valA[1]    = 1.0 / (hy * hy);
480:               /* Missing up term */
481:               col[2].i   = ex - 1;
482:               col[2].j   = ey;
483:               col[2].k   = ez;
484:               col[2].loc = LEFT;
485:               col[2].c   = 0;
486:               valA[2]    = 1.0 / (hx * hx);
487:               col[3].i   = ex + 1;
488:               col[3].j   = ey;
489:               col[3].k   = ez;
490:               col[3].loc = LEFT;
491:               col[3].c   = 0;
492:               valA[3]    = 1.0 / (hx * hx);
493:               col[4].i   = ex;
494:               col[4].j   = ey;
495:               col[4].k   = ez - 1;
496:               col[4].loc = LEFT;
497:               col[4].c   = 0;
498:               valA[4]    = 1.0 / (hz * hz);
499:               /* Missing front term */
500:               col[5].i   = ex - 1;
501:               col[5].j   = ey;
502:               col[5].k   = ez;
503:               col[5].loc = ELEMENT;
504:               col[5].c   = 0;
505:               valA[5]    = 1.0 / hx;
506:               col[6].i   = ex;
507:               col[6].j   = ey;
508:               col[6].k   = ez;
509:               col[6].loc = ELEMENT;
510:               col[6].c   = 0;
511:               valA[6]    = -1.0 / hx;
512:             } else {
513:               nEntries   = 8;
514:               col[0].i   = ex;
515:               col[0].j   = ey;
516:               col[0].k   = ez;
517:               col[0].loc = LEFT;
518:               col[0].c   = 0;
519:               valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
520:               col[1].i   = ex;
521:               col[1].j   = ey - 1;
522:               col[1].k   = ez;
523:               col[1].loc = LEFT;
524:               col[1].c   = 0;
525:               valA[1]    = 1.0 / (hy * hy);
526:               /* Missing up term */
527:               col[2].i   = ex - 1;
528:               col[2].j   = ey;
529:               col[2].k   = ez;
530:               col[2].loc = LEFT;
531:               col[2].c   = 0;
532:               valA[2]    = 1.0 / (hx * hx);
533:               col[3].i   = ex + 1;
534:               col[3].j   = ey;
535:               col[3].k   = ez;
536:               col[3].loc = LEFT;
537:               col[3].c   = 0;
538:               valA[3]    = 1.0 / (hx * hx);
539:               col[4].i   = ex;
540:               col[4].j   = ey;
541:               col[4].k   = ez - 1;
542:               col[4].loc = LEFT;
543:               col[4].c   = 0;
544:               valA[4]    = 1.0 / (hz * hz);
545:               col[5].i   = ex;
546:               col[5].j   = ey;
547:               col[5].k   = ez + 1;
548:               col[5].loc = LEFT;
549:               col[5].c   = 0;
550:               valA[5]    = 1.0 / (hz * hz);
551:               col[6].i   = ex - 1;
552:               col[6].j   = ey;
553:               col[6].k   = ez;
554:               col[6].loc = ELEMENT;
555:               col[6].c   = 0;
556:               valA[6]    = 1.0 / hx;
557:               col[7].i   = ex;
558:               col[7].j   = ey;
559:               col[7].k   = ez;
560:               col[7].loc = ELEMENT;
561:               col[7].c   = 0;
562:               valA[7]    = -1.0 / hx;
563:             }
564:           } else if (ez == 0) {
565:             nEntries   = 8;
566:             col[0].i   = ex;
567:             col[0].j   = ey;
568:             col[0].k   = ez;
569:             col[0].loc = LEFT;
570:             col[0].c   = 0;
571:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
572:             col[1].i   = ex;
573:             col[1].j   = ey - 1;
574:             col[1].k   = ez;
575:             col[1].loc = LEFT;
576:             col[1].c   = 0;
577:             valA[1]    = 1.0 / (hy * hy);
578:             col[2].i   = ex;
579:             col[2].j   = ey + 1;
580:             col[2].k   = ez;
581:             col[2].loc = LEFT;
582:             col[2].c   = 0;
583:             valA[2]    = 1.0 / (hy * hy);
584:             col[3].i   = ex - 1;
585:             col[3].j   = ey;
586:             col[3].k   = ez;
587:             col[3].loc = LEFT;
588:             col[3].c   = 0;
589:             valA[3]    = 1.0 / (hx * hx);
590:             col[4].i   = ex + 1;
591:             col[4].j   = ey;
592:             col[4].k   = ez;
593:             col[4].loc = LEFT;
594:             col[4].c   = 0;
595:             valA[4]    = 1.0 / (hx * hx);
596:             /* Missing back term */
597:             col[5].i   = ex;
598:             col[5].j   = ey;
599:             col[5].k   = ez + 1;
600:             col[5].loc = LEFT;
601:             col[5].c   = 0;
602:             valA[5]    = 1.0 / (hz * hz);
603:             col[6].i   = ex - 1;
604:             col[6].j   = ey;
605:             col[6].k   = ez;
606:             col[6].loc = ELEMENT;
607:             col[6].c   = 0;
608:             valA[6]    = 1.0 / hx;
609:             col[7].i   = ex;
610:             col[7].j   = ey;
611:             col[7].k   = ez;
612:             col[7].loc = ELEMENT;
613:             col[7].c   = 0;
614:             valA[7]    = -1.0 / hx;
615:           } else if (ez == N[2] - 1) {
616:             nEntries   = 8;
617:             col[0].i   = ex;
618:             col[0].j   = ey;
619:             col[0].k   = ez;
620:             col[0].loc = LEFT;
621:             col[0].c   = 0;
622:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
623:             col[1].i   = ex;
624:             col[1].j   = ey - 1;
625:             col[1].k   = ez;
626:             col[1].loc = LEFT;
627:             col[1].c   = 0;
628:             valA[1]    = 1.0 / (hy * hy);
629:             col[2].i   = ex;
630:             col[2].j   = ey + 1;
631:             col[2].k   = ez;
632:             col[2].loc = LEFT;
633:             col[2].c   = 0;
634:             valA[2]    = 1.0 / (hy * hy);
635:             col[3].i   = ex - 1;
636:             col[3].j   = ey;
637:             col[3].k   = ez;
638:             col[3].loc = LEFT;
639:             col[3].c   = 0;
640:             valA[3]    = 1.0 / (hx * hx);
641:             col[4].i   = ex + 1;
642:             col[4].j   = ey;
643:             col[4].k   = ez;
644:             col[4].loc = LEFT;
645:             col[4].c   = 0;
646:             valA[4]    = 1.0 / (hx * hx);
647:             col[5].i   = ex;
648:             col[5].j   = ey;
649:             col[5].k   = ez - 1;
650:             col[5].loc = LEFT;
651:             col[5].c   = 0;
652:             valA[5]    = 1.0 / (hz * hz);
653:             /* Missing front term */
654:             col[6].i   = ex - 1;
655:             col[6].j   = ey;
656:             col[6].k   = ez;
657:             col[6].loc = ELEMENT;
658:             col[6].c   = 0;
659:             valA[6]    = 1.0 / hx;
660:             col[7].i   = ex;
661:             col[7].j   = ey;
662:             col[7].k   = ez;
663:             col[7].loc = ELEMENT;
664:             col[7].c   = 0;
665:             valA[7]    = -1.0 / hx;
666:           } else {
667:             nEntries   = 9;
668:             col[0].i   = ex;
669:             col[0].j   = ey;
670:             col[0].k   = ez;
671:             col[0].loc = LEFT;
672:             col[0].c   = 0;
673:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
674:             col[1].i   = ex;
675:             col[1].j   = ey - 1;
676:             col[1].k   = ez;
677:             col[1].loc = LEFT;
678:             col[1].c   = 0;
679:             valA[1]    = 1.0 / (hy * hy);
680:             col[2].i   = ex;
681:             col[2].j   = ey + 1;
682:             col[2].k   = ez;
683:             col[2].loc = LEFT;
684:             col[2].c   = 0;
685:             valA[2]    = 1.0 / (hy * hy);
686:             col[3].i   = ex - 1;
687:             col[3].j   = ey;
688:             col[3].k   = ez;
689:             col[3].loc = LEFT;
690:             col[3].c   = 0;
691:             valA[3]    = 1.0 / (hx * hx);
692:             col[4].i   = ex + 1;
693:             col[4].j   = ey;
694:             col[4].k   = ez;
695:             col[4].loc = LEFT;
696:             col[4].c   = 0;
697:             valA[4]    = 1.0 / (hx * hx);
698:             col[5].i   = ex;
699:             col[5].j   = ey;
700:             col[5].k   = ez - 1;
701:             col[5].loc = LEFT;
702:             col[5].c   = 0;
703:             valA[5]    = 1.0 / (hz * hz);
704:             col[6].i   = ex;
705:             col[6].j   = ey;
706:             col[6].k   = ez + 1;
707:             col[6].loc = LEFT;
708:             col[6].c   = 0;
709:             valA[6]    = 1.0 / (hz * hz);
710:             col[7].i   = ex - 1;
711:             col[7].j   = ey;
712:             col[7].k   = ez;
713:             col[7].loc = ELEMENT;
714:             col[7].c   = 0;
715:             valA[7]    = 1.0 / hx;
716:             col[8].i   = ex;
717:             col[8].j   = ey;
718:             col[8].k   = ez;
719:             col[8].loc = ELEMENT;
720:             col[8].c   = 0;
721:             valA[8]    = -1.0 / hx;
722:           }
723:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
724:           valRhs = fx(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
725:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
726:         }

728:         /* Equation on bottom face of this element */
729:         if (ey == 0) {
730:           /* Bottom boundary velocity Dirichlet */
731:           DMStagStencil     row;
732:           PetscScalar       valRhs;
733:           const PetscScalar valA = 1.0;
734:           row.i                  = ex;
735:           row.j                  = ey;
736:           row.k                  = ez;
737:           row.loc                = DOWN;
738:           row.c                  = 0;
739:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
740:           valRhs = uyRef(arrCoord[ez][ey][ex][icuy[0]], arrCoord[ez][ey][ex][icuy[1]], arrCoord[ez][ey][ex][icuy[2]]);
741:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
742:         } else {
743:           /* Y-momentum equation, (v_xx + v_yy + v_zz) - p_y = f^y */
744:           DMStagStencil row, col[9];
745:           PetscScalar   valA[9], valRhs;
746:           PetscInt      nEntries;

748:           row.i   = ex;
749:           row.j   = ey;
750:           row.k   = ez;
751:           row.loc = DOWN;
752:           row.c   = 0;
753:           if (ex == 0) {
754:             if (ez == 0) {
755:               nEntries   = 7;
756:               col[0].i   = ex;
757:               col[0].j   = ey;
758:               col[0].k   = ez;
759:               col[0].loc = DOWN;
760:               col[0].c   = 0;
761:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
762:               col[1].i   = ex;
763:               col[1].j   = ey - 1;
764:               col[1].k   = ez;
765:               col[1].loc = DOWN;
766:               col[1].c   = 0;
767:               valA[1]    = 1.0 / (hy * hy);
768:               col[2].i   = ex;
769:               col[2].j   = ey + 1;
770:               col[2].k   = ez;
771:               col[2].loc = DOWN;
772:               col[2].c   = 0;
773:               valA[2]    = 1.0 / (hy * hy);
774:               /* Left term missing */
775:               col[3].i   = ex + 1;
776:               col[3].j   = ey;
777:               col[3].k   = ez;
778:               col[3].loc = DOWN;
779:               col[3].c   = 0;
780:               valA[3]    = 1.0 / (hx * hx);
781:               /* Back term missing */
782:               col[4].i   = ex;
783:               col[4].j   = ey;
784:               col[4].k   = ez + 1;
785:               col[4].loc = DOWN;
786:               col[4].c   = 0;
787:               valA[4]    = 1.0 / (hz * hz);
788:               col[5].i   = ex;
789:               col[5].j   = ey - 1;
790:               col[5].k   = ez;
791:               col[5].loc = ELEMENT;
792:               col[5].c   = 0;
793:               valA[5]    = 1.0 / hy;
794:               col[6].i   = ex;
795:               col[6].j   = ey;
796:               col[6].k   = ez;
797:               col[6].loc = ELEMENT;
798:               col[6].c   = 0;
799:               valA[6]    = -1.0 / hy;
800:             } else if (ez == N[2] - 1) {
801:               nEntries   = 7;
802:               col[0].i   = ex;
803:               col[0].j   = ey;
804:               col[0].k   = ez;
805:               col[0].loc = DOWN;
806:               col[0].c   = 0;
807:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
808:               col[1].i   = ex;
809:               col[1].j   = ey - 1;
810:               col[1].k   = ez;
811:               col[1].loc = DOWN;
812:               col[1].c   = 0;
813:               valA[1]    = 1.0 / (hy * hy);
814:               col[2].i   = ex;
815:               col[2].j   = ey + 1;
816:               col[2].k   = ez;
817:               col[2].loc = DOWN;
818:               col[2].c   = 0;
819:               valA[2]    = 1.0 / (hy * hy);
820:               /* Left term missing */
821:               col[3].i   = ex + 1;
822:               col[3].j   = ey;
823:               col[3].k   = ez;
824:               col[3].loc = DOWN;
825:               col[3].c   = 0;
826:               valA[3]    = 1.0 / (hx * hx);
827:               col[4].i   = ex;
828:               col[4].j   = ey;
829:               col[4].k   = ez - 1;
830:               col[4].loc = DOWN;
831:               col[4].c   = 0;
832:               valA[4]    = 1.0 / (hz * hz);
833:               /* Front term missing */
834:               col[5].i   = ex;
835:               col[5].j   = ey - 1;
836:               col[5].k   = ez;
837:               col[5].loc = ELEMENT;
838:               col[5].c   = 0;
839:               valA[5]    = 1.0 / hy;
840:               col[6].i   = ex;
841:               col[6].j   = ey;
842:               col[6].k   = ez;
843:               col[6].loc = ELEMENT;
844:               col[6].c   = 0;
845:               valA[6]    = -1.0 / hy;
846:             } else {
847:               nEntries   = 8;
848:               col[0].i   = ex;
849:               col[0].j   = ey;
850:               col[0].k   = ez;
851:               col[0].loc = DOWN;
852:               col[0].c   = 0;
853:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
854:               col[1].i   = ex;
855:               col[1].j   = ey - 1;
856:               col[1].k   = ez;
857:               col[1].loc = DOWN;
858:               col[1].c   = 0;
859:               valA[1]    = 1.0 / (hy * hy);
860:               col[2].i   = ex;
861:               col[2].j   = ey + 1;
862:               col[2].k   = ez;
863:               col[2].loc = DOWN;
864:               col[2].c   = 0;
865:               valA[2]    = 1.0 / (hy * hy);
866:               /* Left term missing */
867:               col[3].i   = ex + 1;
868:               col[3].j   = ey;
869:               col[3].k   = ez;
870:               col[3].loc = DOWN;
871:               col[3].c   = 0;
872:               valA[3]    = 1.0 / (hx * hx);
873:               col[4].i   = ex;
874:               col[4].j   = ey;
875:               col[4].k   = ez - 1;
876:               col[4].loc = DOWN;
877:               col[4].c   = 0;
878:               valA[4]    = 1.0 / (hz * hz);
879:               col[5].i   = ex;
880:               col[5].j   = ey;
881:               col[5].k   = ez + 1;
882:               col[5].loc = DOWN;
883:               col[5].c   = 0;
884:               valA[5]    = 1.0 / (hz * hz);
885:               col[6].i   = ex;
886:               col[6].j   = ey - 1;
887:               col[6].k   = ez;
888:               col[6].loc = ELEMENT;
889:               col[6].c   = 0;
890:               valA[6]    = 1.0 / hy;
891:               col[7].i   = ex;
892:               col[7].j   = ey;
893:               col[7].k   = ez;
894:               col[7].loc = ELEMENT;
895:               col[7].c   = 0;
896:               valA[7]    = -1.0 / hy;
897:             }
898:           } else if (ex == N[0] - 1) {
899:             if (ez == 0) {
900:               nEntries   = 7;
901:               col[0].i   = ex;
902:               col[0].j   = ey;
903:               col[0].k   = ez;
904:               col[0].loc = DOWN;
905:               col[0].c   = 0;
906:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
907:               col[1].i   = ex;
908:               col[1].j   = ey - 1;
909:               col[1].k   = ez;
910:               col[1].loc = DOWN;
911:               col[1].c   = 0;
912:               valA[1]    = 1.0 / (hy * hy);
913:               col[2].i   = ex;
914:               col[2].j   = ey + 1;
915:               col[2].k   = ez;
916:               col[2].loc = DOWN;
917:               col[2].c   = 0;
918:               valA[2]    = 1.0 / (hy * hy);
919:               col[3].i   = ex - 1;
920:               col[3].j   = ey;
921:               col[3].k   = ez;
922:               col[3].loc = DOWN;
923:               col[3].c   = 0;
924:               valA[3]    = 1.0 / (hx * hx);
925:               /* Right term missing */
926:               /* Back term missing */
927:               col[4].i   = ex;
928:               col[4].j   = ey;
929:               col[4].k   = ez + 1;
930:               col[4].loc = DOWN;
931:               col[4].c   = 0;
932:               valA[4]    = 1.0 / (hz * hz);
933:               col[5].i   = ex;
934:               col[5].j   = ey - 1;
935:               col[5].k   = ez;
936:               col[5].loc = ELEMENT;
937:               col[5].c   = 0;
938:               valA[5]    = 1.0 / hy;
939:               col[6].i   = ex;
940:               col[6].j   = ey;
941:               col[6].k   = ez;
942:               col[6].loc = ELEMENT;
943:               col[6].c   = 0;
944:               valA[6]    = -1.0 / hy;
945:             } else if (ez == N[2] - 1) {
946:               nEntries   = 7;
947:               col[0].i   = ex;
948:               col[0].j   = ey;
949:               col[0].k   = ez;
950:               col[0].loc = DOWN;
951:               col[0].c   = 0;
952:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
953:               col[1].i   = ex;
954:               col[1].j   = ey - 1;
955:               col[1].k   = ez;
956:               col[1].loc = DOWN;
957:               col[1].c   = 0;
958:               valA[1]    = 1.0 / (hy * hy);
959:               col[2].i   = ex;
960:               col[2].j   = ey + 1;
961:               col[2].k   = ez;
962:               col[2].loc = DOWN;
963:               col[2].c   = 0;
964:               valA[2]    = 1.0 / (hy * hy);
965:               col[3].i   = ex - 1;
966:               col[3].j   = ey;
967:               col[3].k   = ez;
968:               col[3].loc = DOWN;
969:               col[3].c   = 0;
970:               valA[3]    = 1.0 / (hx * hx);
971:               /* Right term missing */
972:               col[4].i   = ex;
973:               col[4].j   = ey;
974:               col[4].k   = ez - 1;
975:               col[4].loc = DOWN;
976:               col[4].c   = 0;
977:               valA[4]    = 1.0 / (hz * hz);
978:               /* Front term missing */
979:               col[5].i   = ex;
980:               col[5].j   = ey - 1;
981:               col[5].k   = ez;
982:               col[5].loc = ELEMENT;
983:               col[5].c   = 0;
984:               valA[5]    = 1.0 / hy;
985:               col[6].i   = ex;
986:               col[6].j   = ey;
987:               col[6].k   = ez;
988:               col[6].loc = ELEMENT;
989:               col[6].c   = 0;
990:               valA[6]    = -1.0 / hy;
991:             } else {
992:               nEntries   = 8;
993:               col[0].i   = ex;
994:               col[0].j   = ey;
995:               col[0].k   = ez;
996:               col[0].loc = DOWN;
997:               col[0].c   = 0;
998:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
999:               col[1].i   = ex;
1000:               col[1].j   = ey - 1;
1001:               col[1].k   = ez;
1002:               col[1].loc = DOWN;
1003:               col[1].c   = 0;
1004:               valA[1]    = 1.0 / (hy * hy);
1005:               col[2].i   = ex;
1006:               col[2].j   = ey + 1;
1007:               col[2].k   = ez;
1008:               col[2].loc = DOWN;
1009:               col[2].c   = 0;
1010:               valA[2]    = 1.0 / (hy * hy);
1011:               col[3].i   = ex - 1;
1012:               col[3].j   = ey;
1013:               col[3].k   = ez;
1014:               col[3].loc = DOWN;
1015:               col[3].c   = 0;
1016:               valA[3]    = 1.0 / (hx * hx);
1017:               /* Right term missing */
1018:               col[4].i   = ex;
1019:               col[4].j   = ey;
1020:               col[4].k   = ez - 1;
1021:               col[4].loc = DOWN;
1022:               col[4].c   = 0;
1023:               valA[4]    = 1.0 / (hz * hz);
1024:               col[5].i   = ex;
1025:               col[5].j   = ey;
1026:               col[5].k   = ez + 1;
1027:               col[5].loc = DOWN;
1028:               col[5].c   = 0;
1029:               valA[5]    = 1.0 / (hz * hz);
1030:               col[6].i   = ex;
1031:               col[6].j   = ey - 1;
1032:               col[6].k   = ez;
1033:               col[6].loc = ELEMENT;
1034:               col[6].c   = 0;
1035:               valA[6]    = 1.0 / hy;
1036:               col[7].i   = ex;
1037:               col[7].j   = ey;
1038:               col[7].k   = ez;
1039:               col[7].loc = ELEMENT;
1040:               col[7].c   = 0;
1041:               valA[7]    = -1.0 / hy;
1042:             }
1043:           } else if (ez == 0) {
1044:             nEntries   = 8;
1045:             col[0].i   = ex;
1046:             col[0].j   = ey;
1047:             col[0].k   = ez;
1048:             col[0].loc = DOWN;
1049:             col[0].c   = 0;
1050:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1051:             col[1].i   = ex;
1052:             col[1].j   = ey - 1;
1053:             col[1].k   = ez;
1054:             col[1].loc = DOWN;
1055:             col[1].c   = 0;
1056:             valA[1]    = 1.0 / (hy * hy);
1057:             col[2].i   = ex;
1058:             col[2].j   = ey + 1;
1059:             col[2].k   = ez;
1060:             col[2].loc = DOWN;
1061:             col[2].c   = 0;
1062:             valA[2]    = 1.0 / (hy * hy);
1063:             col[3].i   = ex - 1;
1064:             col[3].j   = ey;
1065:             col[3].k   = ez;
1066:             col[3].loc = DOWN;
1067:             col[3].c   = 0;
1068:             valA[3]    = 1.0 / (hx * hx);
1069:             col[4].i   = ex + 1;
1070:             col[4].j   = ey;
1071:             col[4].k   = ez;
1072:             col[4].loc = DOWN;
1073:             col[4].c   = 0;
1074:             valA[4]    = 1.0 / (hx * hx);
1075:             /* Back term missing */
1076:             col[5].i   = ex;
1077:             col[5].j   = ey;
1078:             col[5].k   = ez + 1;
1079:             col[5].loc = DOWN;
1080:             col[5].c   = 0;
1081:             valA[5]    = 1.0 / (hz * hz);
1082:             col[6].i   = ex;
1083:             col[6].j   = ey - 1;
1084:             col[6].k   = ez;
1085:             col[6].loc = ELEMENT;
1086:             col[6].c   = 0;
1087:             valA[6]    = 1.0 / hy;
1088:             col[7].i   = ex;
1089:             col[7].j   = ey;
1090:             col[7].k   = ez;
1091:             col[7].loc = ELEMENT;
1092:             col[7].c   = 0;
1093:             valA[7]    = -1.0 / hy;
1094:           } else if (ez == N[2] - 1) {
1095:             nEntries   = 8;
1096:             col[0].i   = ex;
1097:             col[0].j   = ey;
1098:             col[0].k   = ez;
1099:             col[0].loc = DOWN;
1100:             col[0].c   = 0;
1101:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 1.0 / (hz * hz);
1102:             col[1].i   = ex;
1103:             col[1].j   = ey - 1;
1104:             col[1].k   = ez;
1105:             col[1].loc = DOWN;
1106:             col[1].c   = 0;
1107:             valA[1]    = 1.0 / (hy * hy);
1108:             col[2].i   = ex;
1109:             col[2].j   = ey + 1;
1110:             col[2].k   = ez;
1111:             col[2].loc = DOWN;
1112:             col[2].c   = 0;
1113:             valA[2]    = 1.0 / (hy * hy);
1114:             col[3].i   = ex - 1;
1115:             col[3].j   = ey;
1116:             col[3].k   = ez;
1117:             col[3].loc = DOWN;
1118:             col[3].c   = 0;
1119:             valA[3]    = 1.0 / (hx * hx);
1120:             col[4].i   = ex + 1;
1121:             col[4].j   = ey;
1122:             col[4].k   = ez;
1123:             col[4].loc = DOWN;
1124:             col[4].c   = 0;
1125:             valA[4]    = 1.0 / (hx * hx);
1126:             col[5].i   = ex;
1127:             col[5].j   = ey;
1128:             col[5].k   = ez - 1;
1129:             col[5].loc = DOWN;
1130:             col[5].c   = 0;
1131:             valA[5]    = 1.0 / (hz * hz);
1132:             /* Front term missing */
1133:             col[6].i   = ex;
1134:             col[6].j   = ey - 1;
1135:             col[6].k   = ez;
1136:             col[6].loc = ELEMENT;
1137:             col[6].c   = 0;
1138:             valA[6]    = 1.0 / hy;
1139:             col[7].i   = ex;
1140:             col[7].j   = ey;
1141:             col[7].k   = ez;
1142:             col[7].loc = ELEMENT;
1143:             col[7].c   = 0;
1144:             valA[7]    = -1.0 / hy;
1145:           } else {
1146:             nEntries   = 9;
1147:             col[0].i   = ex;
1148:             col[0].j   = ey;
1149:             col[0].k   = ez;
1150:             col[0].loc = DOWN;
1151:             col[0].c   = 0;
1152:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1153:             col[1].i   = ex;
1154:             col[1].j   = ey - 1;
1155:             col[1].k   = ez;
1156:             col[1].loc = DOWN;
1157:             col[1].c   = 0;
1158:             valA[1]    = 1.0 / (hy * hy);
1159:             col[2].i   = ex;
1160:             col[2].j   = ey + 1;
1161:             col[2].k   = ez;
1162:             col[2].loc = DOWN;
1163:             col[2].c   = 0;
1164:             valA[2]    = 1.0 / (hy * hy);
1165:             col[3].i   = ex - 1;
1166:             col[3].j   = ey;
1167:             col[3].k   = ez;
1168:             col[3].loc = DOWN;
1169:             col[3].c   = 0;
1170:             valA[3]    = 1.0 / (hx * hx);
1171:             col[4].i   = ex + 1;
1172:             col[4].j   = ey;
1173:             col[4].k   = ez;
1174:             col[4].loc = DOWN;
1175:             col[4].c   = 0;
1176:             valA[4]    = 1.0 / (hx * hx);
1177:             col[5].i   = ex;
1178:             col[5].j   = ey;
1179:             col[5].k   = ez - 1;
1180:             col[5].loc = DOWN;
1181:             col[5].c   = 0;
1182:             valA[5]    = 1.0 / (hz * hz);
1183:             col[6].i   = ex;
1184:             col[6].j   = ey;
1185:             col[6].k   = ez + 1;
1186:             col[6].loc = DOWN;
1187:             col[6].c   = 0;
1188:             valA[6]    = 1.0 / (hz * hz);
1189:             col[7].i   = ex;
1190:             col[7].j   = ey - 1;
1191:             col[7].k   = ez;
1192:             col[7].loc = ELEMENT;
1193:             col[7].c   = 0;
1194:             valA[7]    = 1.0 / hy;
1195:             col[8].i   = ex;
1196:             col[8].j   = ey;
1197:             col[8].k   = ez;
1198:             col[8].loc = ELEMENT;
1199:             col[8].c   = 0;
1200:             valA[8]    = -1.0 / hy;
1201:           }
1202:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
1203:           valRhs = fy(arrCoord[ez][ey][ex][icuy[0]], arrCoord[ez][ey][ex][icuy[1]], arrCoord[ez][ey][ex][icuy[2]]);
1204:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
1205:         }

1207:         /* Equation on back face of this element */
1208:         if (ez == 0) {
1209:           /* Back boundary velocity Dirichlet */
1210:           DMStagStencil     row;
1211:           PetscScalar       valRhs;
1212:           const PetscScalar valA = 1.0;
1213:           row.i                  = ex;
1214:           row.j                  = ey;
1215:           row.k                  = ez;
1216:           row.loc                = BACK;
1217:           row.c                  = 0;
1218:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
1219:           valRhs = uzRef(arrCoord[ez][ey][ex][icuz[0]], arrCoord[ez][ey][ex][icuz[1]], arrCoord[ez][ey][ex][icuz[2]]);
1220:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
1221:         } else {
1222:           /* Z-momentum equation, (w_xx + w_yy + w_zz) - p_z = f^z */
1223:           DMStagStencil row, col[9];
1224:           PetscScalar   valA[9], valRhs;
1225:           PetscInt      nEntries;

1227:           row.i   = ex;
1228:           row.j   = ey;
1229:           row.k   = ez;
1230:           row.loc = BACK;
1231:           row.c   = 0;
1232:           if (ex == 0) {
1233:             if (ey == 0) {
1234:               nEntries   = 7;
1235:               col[0].i   = ex;
1236:               col[0].j   = ey;
1237:               col[0].k   = ez;
1238:               col[0].loc = BACK;
1239:               col[0].c   = 0;
1240:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1241:               /* Down term missing */
1242:               col[1].i   = ex;
1243:               col[1].j   = ey + 1;
1244:               col[1].k   = ez;
1245:               col[1].loc = BACK;
1246:               col[1].c   = 0;
1247:               valA[1]    = 1.0 / (hy * hy);
1248:               /* Left term missing */
1249:               col[2].i   = ex + 1;
1250:               col[2].j   = ey;
1251:               col[2].k   = ez;
1252:               col[2].loc = BACK;
1253:               col[2].c   = 0;
1254:               valA[2]    = 1.0 / (hx * hx);
1255:               col[3].i   = ex;
1256:               col[3].j   = ey;
1257:               col[3].k   = ez - 1;
1258:               col[3].loc = BACK;
1259:               col[3].c   = 0;
1260:               valA[3]    = 1.0 / (hz * hz);
1261:               col[4].i   = ex;
1262:               col[4].j   = ey;
1263:               col[4].k   = ez + 1;
1264:               col[4].loc = BACK;
1265:               col[4].c   = 0;
1266:               valA[4]    = 1.0 / (hz * hz);
1267:               col[5].i   = ex;
1268:               col[5].j   = ey;
1269:               col[5].k   = ez - 1;
1270:               col[5].loc = ELEMENT;
1271:               col[5].c   = 0;
1272:               valA[5]    = 1.0 / hz;
1273:               col[6].i   = ex;
1274:               col[6].j   = ey;
1275:               col[6].k   = ez;
1276:               col[6].loc = ELEMENT;
1277:               col[6].c   = 0;
1278:               valA[6]    = -1.0 / hz;
1279:             } else if (ey == N[1] - 1) {
1280:               nEntries   = 7;
1281:               col[0].i   = ex;
1282:               col[0].j   = ey;
1283:               col[0].k   = ez;
1284:               col[0].loc = BACK;
1285:               col[0].c   = 0;
1286:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1287:               col[1].i   = ex;
1288:               col[1].j   = ey - 1;
1289:               col[1].k   = ez;
1290:               col[1].loc = BACK;
1291:               col[1].c   = 0;
1292:               valA[1]    = 1.0 / (hy * hy);
1293:               /* Up term missing */
1294:               /* Left term missing */
1295:               col[2].i   = ex + 1;
1296:               col[2].j   = ey;
1297:               col[2].k   = ez;
1298:               col[2].loc = BACK;
1299:               col[2].c   = 0;
1300:               valA[2]    = 1.0 / (hx * hx);
1301:               col[3].i   = ex;
1302:               col[3].j   = ey;
1303:               col[3].k   = ez - 1;
1304:               col[3].loc = BACK;
1305:               col[3].c   = 0;
1306:               valA[3]    = 1.0 / (hz * hz);
1307:               col[4].i   = ex;
1308:               col[4].j   = ey;
1309:               col[4].k   = ez + 1;
1310:               col[4].loc = BACK;
1311:               col[4].c   = 0;
1312:               valA[4]    = 1.0 / (hz * hz);
1313:               col[5].i   = ex;
1314:               col[5].j   = ey;
1315:               col[5].k   = ez - 1;
1316:               col[5].loc = ELEMENT;
1317:               col[5].c   = 0;
1318:               valA[5]    = 1.0 / hz;
1319:               col[6].i   = ex;
1320:               col[6].j   = ey;
1321:               col[6].k   = ez;
1322:               col[6].loc = ELEMENT;
1323:               col[6].c   = 0;
1324:               valA[6]    = -1.0 / hz;
1325:             } else {
1326:               nEntries   = 8;
1327:               col[0].i   = ex;
1328:               col[0].j   = ey;
1329:               col[0].k   = ez;
1330:               col[0].loc = BACK;
1331:               col[0].c   = 0;
1332:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1333:               col[1].i   = ex;
1334:               col[1].j   = ey - 1;
1335:               col[1].k   = ez;
1336:               col[1].loc = BACK;
1337:               col[1].c   = 0;
1338:               valA[1]    = 1.0 / (hy * hy);
1339:               col[2].i   = ex;
1340:               col[2].j   = ey + 1;
1341:               col[2].k   = ez;
1342:               col[2].loc = BACK;
1343:               col[2].c   = 0;
1344:               valA[2]    = 1.0 / (hy * hy);
1345:               /* Left term missing */
1346:               col[3].i   = ex + 1;
1347:               col[3].j   = ey;
1348:               col[3].k   = ez;
1349:               col[3].loc = BACK;
1350:               col[3].c   = 0;
1351:               valA[3]    = 1.0 / (hx * hx);
1352:               col[4].i   = ex;
1353:               col[4].j   = ey;
1354:               col[4].k   = ez - 1;
1355:               col[4].loc = BACK;
1356:               col[4].c   = 0;
1357:               valA[4]    = 1.0 / (hz * hz);
1358:               col[5].i   = ex;
1359:               col[5].j   = ey;
1360:               col[5].k   = ez + 1;
1361:               col[5].loc = BACK;
1362:               col[5].c   = 0;
1363:               valA[5]    = 1.0 / (hz * hz);
1364:               col[6].i   = ex;
1365:               col[6].j   = ey;
1366:               col[6].k   = ez - 1;
1367:               col[6].loc = ELEMENT;
1368:               col[6].c   = 0;
1369:               valA[6]    = 1.0 / hz;
1370:               col[7].i   = ex;
1371:               col[7].j   = ey;
1372:               col[7].k   = ez;
1373:               col[7].loc = ELEMENT;
1374:               col[7].c   = 0;
1375:               valA[7]    = -1.0 / hz;
1376:             }
1377:           } else if (ex == N[0] - 1) {
1378:             if (ey == 0) {
1379:               nEntries   = 7;
1380:               col[0].i   = ex;
1381:               col[0].j   = ey;
1382:               col[0].k   = ez;
1383:               col[0].loc = BACK;
1384:               col[0].c   = 0;
1385:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1386:               /* Down term missing */
1387:               col[1].i   = ex;
1388:               col[1].j   = ey + 1;
1389:               col[1].k   = ez;
1390:               col[1].loc = BACK;
1391:               col[1].c   = 0;
1392:               valA[1]    = 1.0 / (hy * hy);
1393:               col[2].i   = ex - 1;
1394:               col[2].j   = ey;
1395:               col[2].k   = ez;
1396:               col[2].loc = BACK;
1397:               col[2].c   = 0;
1398:               valA[2]    = 1.0 / (hx * hx);
1399:               /* Right term missing */
1400:               col[3].i   = ex;
1401:               col[3].j   = ey;
1402:               col[3].k   = ez - 1;
1403:               col[3].loc = BACK;
1404:               col[3].c   = 0;
1405:               valA[3]    = 1.0 / (hz * hz);
1406:               col[4].i   = ex;
1407:               col[4].j   = ey;
1408:               col[4].k   = ez + 1;
1409:               col[4].loc = BACK;
1410:               col[4].c   = 0;
1411:               valA[4]    = 1.0 / (hz * hz);
1412:               col[5].i   = ex;
1413:               col[5].j   = ey;
1414:               col[5].k   = ez - 1;
1415:               col[5].loc = ELEMENT;
1416:               col[5].c   = 0;
1417:               valA[5]    = 1.0 / hz;
1418:               col[6].i   = ex;
1419:               col[6].j   = ey;
1420:               col[6].k   = ez;
1421:               col[6].loc = ELEMENT;
1422:               col[6].c   = 0;
1423:               valA[6]    = -1.0 / hz;
1424:             } else if (ey == N[1] - 1) {
1425:               nEntries   = 7;
1426:               col[0].i   = ex;
1427:               col[0].j   = ey;
1428:               col[0].k   = ez;
1429:               col[0].loc = BACK;
1430:               col[0].c   = 0;
1431:               valA[0]    = -1.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1432:               col[1].i   = ex;
1433:               col[1].j   = ey - 1;
1434:               col[1].k   = ez;
1435:               col[1].loc = BACK;
1436:               col[1].c   = 0;
1437:               valA[1]    = 1.0 / (hy * hy);
1438:               /* Up term missing */
1439:               col[2].i   = ex - 1;
1440:               col[2].j   = ey;
1441:               col[2].k   = ez;
1442:               col[2].loc = BACK;
1443:               col[2].c   = 0;
1444:               valA[2]    = 1.0 / (hx * hx);
1445:               /* Right term missing */
1446:               col[3].i   = ex;
1447:               col[3].j   = ey;
1448:               col[3].k   = ez - 1;
1449:               col[3].loc = BACK;
1450:               col[3].c   = 0;
1451:               valA[3]    = 1.0 / (hz * hz);
1452:               col[4].i   = ex;
1453:               col[4].j   = ey;
1454:               col[4].k   = ez + 1;
1455:               col[4].loc = BACK;
1456:               col[4].c   = 0;
1457:               valA[4]    = 1.0 / (hz * hz);
1458:               col[5].i   = ex;
1459:               col[5].j   = ey;
1460:               col[5].k   = ez - 1;
1461:               col[5].loc = ELEMENT;
1462:               col[5].c   = 0;
1463:               valA[5]    = 1.0 / hz;
1464:               col[6].i   = ex;
1465:               col[6].j   = ey;
1466:               col[6].k   = ez;
1467:               col[6].loc = ELEMENT;
1468:               col[6].c   = 0;
1469:               valA[6]    = -1.0 / hz;
1470:             } else {
1471:               nEntries   = 8;
1472:               col[0].i   = ex;
1473:               col[0].j   = ey;
1474:               col[0].k   = ez;
1475:               col[0].loc = BACK;
1476:               col[0].c   = 0;
1477:               valA[0]    = -1.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1478:               col[1].i   = ex;
1479:               col[1].j   = ey - 1;
1480:               col[1].k   = ez;
1481:               col[1].loc = BACK;
1482:               col[1].c   = 0;
1483:               valA[1]    = 1.0 / (hy * hy);
1484:               col[2].i   = ex;
1485:               col[2].j   = ey + 1;
1486:               col[2].k   = ez;
1487:               col[2].loc = BACK;
1488:               col[2].c   = 0;
1489:               valA[2]    = 1.0 / (hy * hy);
1490:               col[3].i   = ex - 1;
1491:               col[3].j   = ey;
1492:               col[3].k   = ez;
1493:               col[3].loc = BACK;
1494:               col[3].c   = 0;
1495:               valA[3]    = 1.0 / (hx * hx);
1496:               /* Right term missing */
1497:               col[4].i   = ex;
1498:               col[4].j   = ey;
1499:               col[4].k   = ez - 1;
1500:               col[4].loc = BACK;
1501:               col[4].c   = 0;
1502:               valA[4]    = 1.0 / (hz * hz);
1503:               col[5].i   = ex;
1504:               col[5].j   = ey;
1505:               col[5].k   = ez + 1;
1506:               col[5].loc = BACK;
1507:               col[5].c   = 0;
1508:               valA[5]    = 1.0 / (hz * hz);
1509:               col[6].i   = ex;
1510:               col[6].j   = ey;
1511:               col[6].k   = ez - 1;
1512:               col[6].loc = ELEMENT;
1513:               col[6].c   = 0;
1514:               valA[6]    = 1.0 / hz;
1515:               col[7].i   = ex;
1516:               col[7].j   = ey;
1517:               col[7].k   = ez;
1518:               col[7].loc = ELEMENT;
1519:               col[7].c   = 0;
1520:               valA[7]    = -1.0 / hz;
1521:             }
1522:           } else if (ey == 0) {
1523:             nEntries   = 8;
1524:             col[0].i   = ex;
1525:             col[0].j   = ey;
1526:             col[0].k   = ez;
1527:             col[0].loc = BACK;
1528:             col[0].c   = 0;
1529:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1530:             /* Down term missing */
1531:             col[1].i   = ex;
1532:             col[1].j   = ey + 1;
1533:             col[1].k   = ez;
1534:             col[1].loc = BACK;
1535:             col[1].c   = 0;
1536:             valA[1]    = 1.0 / (hy * hy);
1537:             col[2].i   = ex - 1;
1538:             col[2].j   = ey;
1539:             col[2].k   = ez;
1540:             col[2].loc = BACK;
1541:             col[2].c   = 0;
1542:             valA[2]    = 1.0 / (hx * hx);
1543:             col[3].i   = ex + 1;
1544:             col[3].j   = ey;
1545:             col[3].k   = ez;
1546:             col[3].loc = BACK;
1547:             col[3].c   = 0;
1548:             valA[3]    = 1.0 / (hx * hx);
1549:             col[4].i   = ex;
1550:             col[4].j   = ey;
1551:             col[4].k   = ez - 1;
1552:             col[4].loc = BACK;
1553:             col[4].c   = 0;
1554:             valA[4]    = 1.0 / (hz * hz);
1555:             col[5].i   = ex;
1556:             col[5].j   = ey;
1557:             col[5].k   = ez + 1;
1558:             col[5].loc = BACK;
1559:             col[5].c   = 0;
1560:             valA[5]    = 1.0 / (hz * hz);
1561:             col[6].i   = ex;
1562:             col[6].j   = ey;
1563:             col[6].k   = ez - 1;
1564:             col[6].loc = ELEMENT;
1565:             col[6].c   = 0;
1566:             valA[6]    = 1.0 / hz;
1567:             col[7].i   = ex;
1568:             col[7].j   = ey;
1569:             col[7].k   = ez;
1570:             col[7].loc = ELEMENT;
1571:             col[7].c   = 0;
1572:             valA[7]    = -1.0 / hz;
1573:           } else if (ey == N[1] - 1) {
1574:             nEntries   = 8;
1575:             col[0].i   = ex;
1576:             col[0].j   = ey;
1577:             col[0].k   = ez;
1578:             col[0].loc = BACK;
1579:             col[0].c   = 0;
1580:             valA[0]    = -2.0 / (hx * hx) + -1.0 / (hy * hy) - 2.0 / (hz * hz);
1581:             col[1].i   = ex;
1582:             col[1].j   = ey - 1;
1583:             col[1].k   = ez;
1584:             col[1].loc = BACK;
1585:             col[1].c   = 0;
1586:             valA[1]    = 1.0 / (hy * hy);
1587:             /* Up term missing */
1588:             col[2].i   = ex - 1;
1589:             col[2].j   = ey;
1590:             col[2].k   = ez;
1591:             col[2].loc = BACK;
1592:             col[2].c   = 0;
1593:             valA[2]    = 1.0 / (hx * hx);
1594:             col[3].i   = ex + 1;
1595:             col[3].j   = ey;
1596:             col[3].k   = ez;
1597:             col[3].loc = BACK;
1598:             col[3].c   = 0;
1599:             valA[3]    = 1.0 / (hx * hx);
1600:             col[4].i   = ex;
1601:             col[4].j   = ey;
1602:             col[4].k   = ez - 1;
1603:             col[4].loc = BACK;
1604:             col[4].c   = 0;
1605:             valA[4]    = 1.0 / (hz * hz);
1606:             col[5].i   = ex;
1607:             col[5].j   = ey;
1608:             col[5].k   = ez + 1;
1609:             col[5].loc = BACK;
1610:             col[5].c   = 0;
1611:             valA[5]    = 1.0 / (hz * hz);
1612:             col[6].i   = ex;
1613:             col[6].j   = ey;
1614:             col[6].k   = ez - 1;
1615:             col[6].loc = ELEMENT;
1616:             col[6].c   = 0;
1617:             valA[6]    = 1.0 / hz;
1618:             col[7].i   = ex;
1619:             col[7].j   = ey;
1620:             col[7].k   = ez;
1621:             col[7].loc = ELEMENT;
1622:             col[7].c   = 0;
1623:             valA[7]    = -1.0 / hz;
1624:           } else {
1625:             nEntries   = 9;
1626:             col[0].i   = ex;
1627:             col[0].j   = ey;
1628:             col[0].k   = ez;
1629:             col[0].loc = BACK;
1630:             col[0].c   = 0;
1631:             valA[0]    = -2.0 / (hx * hx) + -2.0 / (hy * hy) - 2.0 / (hz * hz);
1632:             col[1].i   = ex;
1633:             col[1].j   = ey - 1;
1634:             col[1].k   = ez;
1635:             col[1].loc = BACK;
1636:             col[1].c   = 0;
1637:             valA[1]    = 1.0 / (hy * hy);
1638:             col[2].i   = ex;
1639:             col[2].j   = ey + 1;
1640:             col[2].k   = ez;
1641:             col[2].loc = BACK;
1642:             col[2].c   = 0;
1643:             valA[2]    = 1.0 / (hy * hy);
1644:             col[3].i   = ex - 1;
1645:             col[3].j   = ey;
1646:             col[3].k   = ez;
1647:             col[3].loc = BACK;
1648:             col[3].c   = 0;
1649:             valA[3]    = 1.0 / (hx * hx);
1650:             col[4].i   = ex + 1;
1651:             col[4].j   = ey;
1652:             col[4].k   = ez;
1653:             col[4].loc = BACK;
1654:             col[4].c   = 0;
1655:             valA[4]    = 1.0 / (hx * hx);
1656:             col[5].i   = ex;
1657:             col[5].j   = ey;
1658:             col[5].k   = ez - 1;
1659:             col[5].loc = BACK;
1660:             col[5].c   = 0;
1661:             valA[5]    = 1.0 / (hz * hz);
1662:             col[6].i   = ex;
1663:             col[6].j   = ey;
1664:             col[6].k   = ez + 1;
1665:             col[6].loc = BACK;
1666:             col[6].c   = 0;
1667:             valA[6]    = 1.0 / (hz * hz);
1668:             col[7].i   = ex;
1669:             col[7].j   = ey;
1670:             col[7].k   = ez - 1;
1671:             col[7].loc = ELEMENT;
1672:             col[7].c   = 0;
1673:             valA[7]    = 1.0 / hz;
1674:             col[8].i   = ex;
1675:             col[8].j   = ey;
1676:             col[8].k   = ez;
1677:             col[8].loc = ELEMENT;
1678:             col[8].c   = 0;
1679:             valA[8]    = -1.0 / hz;
1680:           }
1681:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, nEntries, col, valA, INSERT_VALUES);
1682:           valRhs = fz(arrCoord[ez][ey][ex][icuz[0]], arrCoord[ez][ey][ex][icuz[1]], arrCoord[ez][ey][ex][icuz[2]]);
1683:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
1684:         }

1686:         /* P equation : u_x + v_y + w_z = g
1687:            Note that this includes an explicit zero on the diagonal. This is only needed for
1688:            direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
1689:         if (pinPressure && ex == 0 && ey == 0 && ez == 0) { /* Pin the first pressure node, if requested */
1690:           DMStagStencil row;
1691:           PetscScalar   valA, valRhs;
1692:           row.i   = ex;
1693:           row.j   = ey;
1694:           row.k   = ez;
1695:           row.loc = ELEMENT;
1696:           row.c   = 0;
1697:           valA    = 1.0;
1698:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 1, &row, &valA, INSERT_VALUES);
1699:           valRhs = pRef(arrCoord[ez][ey][ex][icp[0]], arrCoord[ez][ey][ex][icp[1]], arrCoord[ez][ey][ex][icp[2]]);
1700:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
1701:         } else {
1702:           DMStagStencil row, col[7];
1703:           PetscScalar   valA[7], valRhs;

1705:           row.i      = ex;
1706:           row.j      = ey;
1707:           row.k      = ez;
1708:           row.loc    = ELEMENT;
1709:           row.c      = 0;
1710:           col[0].i   = ex;
1711:           col[0].j   = ey;
1712:           col[0].k   = ez;
1713:           col[0].loc = LEFT;
1714:           col[0].c   = 0;
1715:           valA[0]    = -1.0 / hx;
1716:           col[1].i   = ex;
1717:           col[1].j   = ey;
1718:           col[1].k   = ez;
1719:           col[1].loc = RIGHT;
1720:           col[1].c   = 0;
1721:           valA[1]    = 1.0 / hx;
1722:           col[2].i   = ex;
1723:           col[2].j   = ey;
1724:           col[2].k   = ez;
1725:           col[2].loc = DOWN;
1726:           col[2].c   = 0;
1727:           valA[2]    = -1.0 / hy;
1728:           col[3].i   = ex;
1729:           col[3].j   = ey;
1730:           col[3].k   = ez;
1731:           col[3].loc = UP;
1732:           col[3].c   = 0;
1733:           valA[3]    = 1.0 / hy;
1734:           col[4].i   = ex;
1735:           col[4].j   = ey;
1736:           col[4].k   = ez;
1737:           col[4].loc = BACK;
1738:           col[4].c   = 0;
1739:           valA[4]    = -1.0 / hz;
1740:           col[5].i   = ex;
1741:           col[5].j   = ey;
1742:           col[5].k   = ez;
1743:           col[5].loc = FRONT;
1744:           col[5].c   = 0;
1745:           valA[5]    = 1.0 / hz;
1746:           col[6]     = row;
1747:           valA[6]    = 0.0;
1748:           DMStagMatSetValuesStencil(dmSol, A, 1, &row, 7, col, valA, INSERT_VALUES);
1749:           valRhs = g(arrCoord[ez][ey][ex][icp[0]], arrCoord[ez][ey][ex][icp[1]], arrCoord[ez][ey][ex][icp[2]]);
1750:           DMStagVecSetValuesStencil(dmSol, rhs, 1, &row, &valRhs, INSERT_VALUES);
1751:         }
1752:       }
1753:     }
1754:   }
1755:   DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord);
1756:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1757:   VecAssemblyBegin(rhs);
1758:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1759:   VecAssemblyEnd(rhs);

1761:   return 0;
1762: }

1764: /* Create a pressure-only DMStag and use it to generate a nullspace vector
1765:    - Create a compatible DMStag with one dof per element (and nothing else).
1766:    - Create a constant vector and normalize it
1767:    - Migrate it to a vector on the original dmSol (making use of the fact
1768:    that this will fill in zeros for "extra" dof)
1769:    - Set the nullspace for the operator
1770:    - Destroy everything (the operator keeps the references it needs) */
1771: static PetscErrorCode AttachNullspace(DM dmSol, Mat A)
1772: {
1773:   DM           dmPressure;
1774:   Vec          constantPressure, basis;
1775:   PetscReal    nrm;
1776:   MatNullSpace matNullSpace;

1779:   DMStagCreateCompatibleDMStag(dmSol, 0, 0, 0, 1, &dmPressure);
1780:   DMGetGlobalVector(dmPressure, &constantPressure);
1781:   VecSet(constantPressure, 1.0);
1782:   VecNorm(constantPressure, NORM_2, &nrm);
1783:   VecScale(constantPressure, 1.0 / nrm);
1784:   DMCreateGlobalVector(dmSol, &basis);
1785:   DMStagMigrateVec(dmPressure, constantPressure, dmSol, basis);
1786:   MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol), PETSC_FALSE, 1, &basis, &matNullSpace);
1787:   VecDestroy(&basis);
1788:   VecDestroy(&constantPressure);
1789:   MatSetNullSpace(A, matNullSpace);
1790:   MatNullSpaceDestroy(&matNullSpace);
1791:   return 0;
1792: }

1794: static PetscErrorCode CreateReferenceSolution(DM dmSol, Vec *pSolRef)
1795: {
1796:   PetscInt        start[3], n[3], nExtra[3], ex, ey, ez, d;
1797:   PetscInt        ip, iux, iuy, iuz, icp[3], icux[3], icuy[3], icuz[3];
1798:   Vec             solRef, solRefLocal, coord, coordLocal;
1799:   DM              dmCoord;
1800:   PetscScalar ****arrSol, ****arrCoord;

1803:   DMCreateGlobalVector(dmSol, pSolRef);
1804:   solRef = *pSolRef;
1805:   DMStagGetCorners(dmSol, &start[0], &start[1], &start[2], &n[0], &n[1], &n[2], &nExtra[0], &nExtra[1], &nExtra[2]);
1806:   DMGetCoordinateDM(dmSol, &dmCoord);
1807:   DMGetCoordinates(dmSol, &coord);
1808:   DMGetLocalVector(dmCoord, &coordLocal);
1809:   DMGlobalToLocal(dmCoord, coord, INSERT_VALUES, coordLocal);
1810:   DMStagGetLocationSlot(dmSol, ELEMENT, 0, &ip);
1811:   DMStagGetLocationSlot(dmSol, LEFT, 0, &iux);
1812:   DMStagGetLocationSlot(dmSol, DOWN, 0, &iuy);
1813:   DMStagGetLocationSlot(dmSol, BACK, 0, &iuz);
1814:   for (d = 0; d < 3; ++d) {
1815:     DMStagGetLocationSlot(dmCoord, ELEMENT, d, &icp[d]);
1816:     DMStagGetLocationSlot(dmCoord, LEFT, d, &icux[d]);
1817:     DMStagGetLocationSlot(dmCoord, DOWN, d, &icuy[d]);
1818:     DMStagGetLocationSlot(dmCoord, BACK, d, &icuz[d]);
1819:   }
1820:   DMStagVecGetArrayRead(dmCoord, coordLocal, &arrCoord);
1821:   DMGetLocalVector(dmSol, &solRefLocal);
1822:   DMStagVecGetArray(dmSol, solRefLocal, &arrSol);
1823:   for (ez = start[2]; ez < start[2] + n[2] + nExtra[2]; ++ez) {
1824:     for (ey = start[1]; ey < start[1] + n[1] + nExtra[1]; ++ey) {
1825:       for (ex = start[0]; ex < start[0] + n[0] + nExtra[0]; ++ex) {
1826:         if (ex < start[0] + n[0] && ey < start[1] + n[1]) arrSol[ez][ey][ex][iuz] = uzRef(arrCoord[ez][ey][ex][icuz[0]], arrCoord[ez][ey][ex][icuz[1]], arrCoord[ez][ey][ex][icuz[2]]);
1827:         if (ex < start[0] + n[0] && ey < start[2] + n[2]) arrSol[ez][ey][ex][iuy] = uyRef(arrCoord[ez][ey][ex][icuy[0]], arrCoord[ez][ey][ex][icuy[1]], arrCoord[ez][ey][ex][icuy[2]]);
1828:         if (ex < start[1] + n[1] && ey < start[2] + n[2]) arrSol[ez][ey][ex][iux] = uxRef(arrCoord[ez][ey][ex][icux[0]], arrCoord[ez][ey][ex][icux[1]], arrCoord[ez][ey][ex][icux[2]]);
1829:         if (ex < start[0] + n[0] && ey < start[1] + n[1] && ez < start[2] + n[2]) arrSol[ez][ey][ex][ip] = pRef(arrCoord[ez][ey][ex][icp[0]], arrCoord[ez][ey][ex][icp[1]], arrCoord[ez][ey][ex][icp[2]]);
1830:       }
1831:     }
1832:   }
1833:   DMStagVecRestoreArrayRead(dmCoord, coordLocal, &arrCoord);
1834:   DMStagVecRestoreArray(dmSol, solRefLocal, &arrSol);
1835:   DMLocalToGlobal(dmSol, solRefLocal, INSERT_VALUES, solRef);
1836:   DMRestoreLocalVector(dmCoord, &coordLocal);
1837:   DMRestoreLocalVector(dmSol, &solRefLocal);
1838:   return 0;
1839: }

1841: static PetscErrorCode CheckSolution(Vec sol, Vec solRef)
1842: {
1843:   Vec       diff;
1844:   PetscReal normsolRef, errAbs, errRel;

1847:   VecDuplicate(sol, &diff);
1848:   VecCopy(sol, diff);
1849:   VecAXPY(diff, -1.0, solRef);
1850:   VecNorm(diff, NORM_2, &errAbs);
1851:   VecNorm(solRef, NORM_2, &normsolRef);
1852:   errRel = errAbs / normsolRef;
1853:   PetscPrintf(PETSC_COMM_WORLD, "Error (abs): %g\nError (rel): %g\n", (double)errAbs, (double)errRel);
1854:   VecDestroy(&diff);
1855:   return 0;
1856: }

1858: /*TEST

1860:    test:
1861:       suffix: 1
1862:       requires: mumps
1863:       nsize: 27
1864:       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

1866:    test:
1867:       suffix: 2
1868:       requires: !single
1869:       nsize: 4
1870:       args: -ksp_monitor_short -ksp_converged_reason -pc_fieldsplit_schur_fact_type diag -fieldsplit_0_ksp_type preonly -fieldsplit_1_pc_type none -fieldsplit_0_pc_type gamg -fieldsplit_0_ksp_type gmres -fieldsplit_0_pc_gamg_esteig_ksp_max_it 5

1872:    test:
1873:       suffix: direct_umfpack
1874:       requires: suitesparse
1875:       nsize: 1
1876:       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

1878: TEST*/