Actual source code: ex3.c

petsc-master 2019-07-21
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:   PetscBool         isLastRankx,isLastRanky,isLastRankz,isFirstRankx,isFirstRanky,isFirstRankz;
152:   PetscReal         hx,hy,hz;
153:   DM                dmCoord;
154:   PetscScalar       ****arrCoord;

157:   DMCreateMatrix(dmSol,pA);
158:   A = *pA;
159:   DMCreateGlobalVector(dmSol,pRhs);
160:   rhs = *pRhs;

162:   DMStagGetCorners(dmSol,&startx,&starty,&startz,&nx,&ny,&nz,NULL,NULL,NULL);
163:   DMStagGetGlobalSizes(dmSol,&N[0],&N[1],&N[2]);
164:   if (N[0] < 2 || N[1] < 2 || N[2] < 2) SETERRQ(PetscObjectComm((PetscObject)dmSol),PETSC_ERR_ARG_SIZ,"This example requires at least two elements in each dimensions");
165:   DMStagGetIsLastRank(dmSol,&isLastRankx,&isLastRanky,&isLastRankz);
166:   DMStagGetIsFirstRank(dmSol,&isFirstRankx,&isFirstRanky,&isFirstRankz);
167:   hx = 1.0/N[0]; hy = 1.0/N[1]; hz = 1.0/N[2];
168:   DMGetCoordinateDM(dmSol,&dmCoord);
169:   DMGetCoordinatesLocal(dmSol,&coordLocal);
170:   DMStagVecGetArrayDOFRead(dmCoord,coordLocal,&arrCoord);
171:   for (d=0; d<3; ++d) {
172:     DMStagGetLocationSlot(dmCoord,ELEMENT,d,&icp[d]       );
173:     DMStagGetLocationSlot(dmCoord,LEFT,   d,&icux[d]      );
174:     DMStagGetLocationSlot(dmCoord,DOWN,   d,&icuy[d]      );
175:     DMStagGetLocationSlot(dmCoord,BACK,   d,&icuz[d]      );
176:     DMStagGetLocationSlot(dmCoord,RIGHT,  d,&icux_right[d]);
177:     DMStagGetLocationSlot(dmCoord,UP,     d,&icuy_up[d]   );
178:     DMStagGetLocationSlot(dmCoord,FRONT,  d,&icuz_front[d]);
179:   }

181:   /* Loop over all local elements. Note that it may be more efficient in real
182:      applications to loop over each boundary separately */
183:   for (ez = startz; ez<startz+nz; ++ez) { /* With DMStag, always iterate x fastest, y second fastest, z slowest */
184:     for (ey = starty; ey<starty+ny; ++ey) {
185:       for (ex = startx; ex<startx+nx; ++ex) {

187:         if (ex == N[0]-1) {
188:           /* Right Boundary velocity Dirichlet */
189:           DMStagStencil row;
190:           PetscScalar   valRhs;
191:           const PetscScalar valA = 1.0;
192:           row.i = ex; row.j = ey; row.k = ez; row.loc = RIGHT; row.c = 0;
193:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
194:           valRhs = uxRef(
195:               arrCoord[ez][ey][ex][icux_right[0]],
196:               arrCoord[ez][ey][ex][icux_right[1]],
197:               arrCoord[ez][ey][ex][icux_right[2]]
198:               );
199:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
200:         }
201:         if (ey == N[1]-1) {
202:           /* Top boundary velocity Dirichlet */
203:           DMStagStencil row;
204:           PetscScalar   valRhs;
205:           const PetscScalar valA = 1.0;
206:           row.i = ex; row.j = ey; row.k = ez; row.loc = UP; row.c = 0;
207:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
208:           valRhs = uyRef(
209:               arrCoord[ez][ey][ex][icuy_up[0]],
210:               arrCoord[ez][ey][ex][icuy_up[1]],
211:               arrCoord[ez][ey][ex][icuy_up[2]]
212:               );
213:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
214:         }
215:         if (ez == N[2]-1) {
216:           /* Top boundary velocity Dirichlet */
217:           DMStagStencil row;
218:           PetscScalar   valRhs;
219:           const PetscScalar valA = 1.0;
220:           row.i = ex; row.j = ey; row.k = ez; row.loc = FRONT; row.c = 0;
221:           DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
222:           valRhs = uzRef(
223:               arrCoord[ez][ey][ex][icuz_front[0]],
224:               arrCoord[ez][ey][ex][icuz_front[1]],
225:               arrCoord[ez][ey][ex][icuz_front[2]]
226:               );
227:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
228:         }

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

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

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

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

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

518:           row.i = ex; row.j = ey; row.k = ez; row.loc = BACK; row.c = 0;
519:           if (ex == 0) {
520:             if (ey == 0) {
521:               nEntries = 7;
522:               col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
523:               /* Down term missing */
524:               col[1].i = ex  ; col[1].j = ey+1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
525:               /* Left term missing */
526:               col[2].i = ex+1; col[2].j = ey  ;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
527:               col[3].i = ex  ; col[3].j = ey  ;  col[3].k = ez-1; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hz*hz);
528:               col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez+1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
529:               col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] =  1.0 / hz;
530:               col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez  ; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] = -1.0 / hz;
531:             } else if (ey == N[1]-1) {
532:               nEntries = 7;
533:               col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
534:               col[1].i = ex  ; col[1].j = ey-1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
535:               /* Up term missing */
536:               /* Left term missing */
537:               col[2].i = ex+1; col[2].j = ey  ;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
538:               col[3].i = ex  ; col[3].j = ey  ;  col[3].k = ez-1; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hz*hz);
539:               col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez+1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
540:               col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] =  1.0 / hz;
541:               col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez  ; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] = -1.0 / hz;
542:             } else {
543:               nEntries = 8;
544:               col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
545:               col[1].i = ex  ; col[1].j = ey-1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
546:               col[2].i = ex  ; col[2].j = ey+1;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hy*hy);
547:               /* Left term missing */
548:               col[3].i = ex+1; col[3].j = ey  ;  col[3].k = ez  ; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
549:               col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez-1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
550:               col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez+1; col[5].loc = BACK;    col[5].c  = 0; valA[5] =  1.0 / (hz*hz);
551:               col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] =  1.0 / hz;
552:               col[7].i = ex  ; col[7].j = ey  ;  col[7].k = ez  ; col[7].loc = ELEMENT; col[7].c  = 0; valA[7] = -1.0 / hz;
553:             }
554:           } else if (ex == N[0]-1) {
555:             if (ey == 0) {
556:               nEntries = 7;
557:               col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
558:               /* Down term missing */
559:               col[1].i = ex  ; col[1].j = ey+1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
560:               col[2].i = ex-1; col[2].j = ey  ;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
561:               /* Right term missing */
562:               col[3].i = ex  ; col[3].j = ey  ;  col[3].k = ez-1; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hz*hz);
563:               col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez+1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
564:               col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] =  1.0 / hz;
565:               col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez  ; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] = -1.0 / hz;
566:             } else if (ey == N[1]-1) {
567:               nEntries = 7;
568:               col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
569:               col[1].i = ex  ; col[1].j = ey-1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
570:               /* Up term missing */
571:               col[2].i = ex-1; col[2].j = ey  ;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
572:               /* Right term missing */
573:               col[3].i = ex  ; col[3].j = ey  ;  col[3].k = ez-1; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hz*hz);
574:               col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez+1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
575:               col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez-1; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] =  1.0 / hz;
576:               col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez  ; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] = -1.0 / hz;
577:             } else {
578:               nEntries = 8;
579:               col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
580:               col[1].i = ex  ; col[1].j = ey-1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
581:               col[2].i = ex  ; col[2].j = ey+1;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hy*hy);
582:               col[3].i = ex-1; col[3].j = ey  ;  col[3].k = ez  ; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
583:               /* Right term missing */
584:               col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez-1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
585:               col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez+1; col[5].loc = BACK;    col[5].c  = 0; valA[5] =  1.0 / (hz*hz);
586:               col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] =  1.0 / hz;
587:               col[7].i = ex  ; col[7].j = ey  ;  col[7].k = ez  ; col[7].loc = ELEMENT; col[7].c  = 0; valA[7] = -1.0 / hz;
588:             }
589:           } else if (ey == 0) {
590:             nEntries = 8;
591:             col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -2.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
592:             /* Down term missing */
593:             col[1].i = ex  ; col[1].j = ey+1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
594:             col[2].i = ex-1; col[2].j = ey  ;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
595:             col[3].i = ex+1; col[3].j = ey  ;  col[3].k = ez  ; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
596:             col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez-1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
597:             col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez+1; col[5].loc = BACK;    col[5].c  = 0; valA[5] =  1.0 / (hz*hz);
598:             col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] =  1.0 / hz;
599:             col[7].i = ex  ; col[7].j = ey  ;  col[7].k = ez  ; col[7].loc = ELEMENT; col[7].c  = 0; valA[7] = -1.0 / hz;
600:           } else if (ey == N[1]-1) {
601:             nEntries = 8;
602:             col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -2.0 / (hx*hx) + -1.0 / (hy*hy) -2.0 / (hz*hz);
603:             col[1].i = ex  ; col[1].j = ey-1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
604:             /* Up term missing */
605:             col[2].i = ex-1; col[2].j = ey  ;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
606:             col[3].i = ex+1; col[3].j = ey  ;  col[3].k = ez  ; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
607:             col[4].i = ex  ; col[4].j = ey  ;  col[4].k = ez-1; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hz*hz);
608:             col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez+1; col[5].loc = BACK;    col[5].c  = 0; valA[5] =  1.0 / (hz*hz);
609:             col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez-1; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] =  1.0 / hz;
610:             col[7].i = ex  ; col[7].j = ey  ;  col[7].k = ez  ; col[7].loc = ELEMENT; col[7].c  = 0; valA[7] = -1.0 / hz;
611:           } else {
612:             nEntries = 9;
613:             col[0].i = ex  ; col[0].j = ey  ;  col[0].k = ez  ; col[0].loc = BACK;    col[0].c  = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy) -2.0 / (hz*hz);
614:             col[1].i = ex  ; col[1].j = ey-1;  col[1].k = ez  ; col[1].loc = BACK;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
615:             col[2].i = ex  ; col[2].j = ey+1;  col[2].k = ez  ; col[2].loc = BACK;    col[2].c  = 0; valA[2] =  1.0 / (hy*hy);
616:             col[3].i = ex-1; col[3].j = ey  ;  col[3].k = ez  ; col[3].loc = BACK;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
617:             col[4].i = ex+1; col[4].j = ey  ;  col[4].k = ez  ; col[4].loc = BACK;    col[4].c  = 0; valA[4] =  1.0 / (hx*hx);
618:             col[5].i = ex  ; col[5].j = ey  ;  col[5].k = ez-1; col[5].loc = BACK;    col[5].c  = 0; valA[5] =  1.0 / (hz*hz);
619:             col[6].i = ex  ; col[6].j = ey  ;  col[6].k = ez+1; col[6].loc = BACK;    col[6].c  = 0; valA[6] =  1.0 / (hz*hz);
620:             col[7].i = ex  ; col[7].j = ey  ;  col[7].k = ez-1; col[7].loc = ELEMENT; col[7].c  = 0; valA[7] =  1.0 / hz;
621:             col[8].i = ex  ; col[8].j = ey  ;  col[8].k = ez  ; col[8].loc = ELEMENT; col[8].c  = 0; valA[8] = -1.0 / hz;
622:           }
623:           DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
624:           valRhs = fz(
625:               arrCoord[ez][ey][ex][icuz[0]],
626:               arrCoord[ez][ey][ex][icuz[1]],
627:               arrCoord[ez][ey][ex][icuz[2]]
628:               );
629:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
630:         }

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

651:           row.i    = ex; row.j    = ey; row.k    = ez; row.loc    = ELEMENT; row.c    = 0;
652:           col[0].i = ex; col[0].j = ey; col[0].k = ez; col[0].loc = LEFT;    col[0].c = 0; valA[0] = -1.0 / hx;
653:           col[1].i = ex; col[1].j = ey; col[1].k = ez; col[1].loc = RIGHT;   col[1].c = 0; valA[1] =  1.0 / hx;
654:           col[2].i = ex; col[2].j = ey; col[2].k = ez; col[2].loc = DOWN;    col[2].c = 0; valA[2] = -1.0 / hy;
655:           col[3].i = ex; col[3].j = ey; col[3].k = ez; col[3].loc = UP;      col[3].c = 0; valA[3] =  1.0 / hy;
656:           col[4].i = ex; col[4].j = ey; col[4].k = ez; col[4].loc = BACK;    col[4].c = 0; valA[4] = -1.0 / hz;
657:           col[5].i = ex; col[5].j = ey; col[5].k = ez; col[5].loc = FRONT;   col[5].c = 0; valA[5] =  1.0 / hz;
658:           col[6]   = row;                                                                  valA[6] =  0.0;
659:           DMStagMatSetValuesStencil(dmSol,A,1,&row,7,col,valA,INSERT_VALUES);
660:           valRhs = g(
661:               arrCoord[ez][ey][ex][icp[0]],
662:               arrCoord[ez][ey][ex][icp[1]],
663:               arrCoord[ez][ey][ex][icp[2]]
664:               );
665:           DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
666:         }
667:       }
668:     }
669:   }
670:   DMStagVecRestoreArrayDOFRead(dmCoord,coordLocal,&arrCoord);
671:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
672:   VecAssemblyBegin(rhs);
673:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
674:   VecAssemblyEnd(rhs);

676:   return(0);
677: }

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

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

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

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

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

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

796: /*TEST

798:    test:
799:       suffix: 1
800:       requires: mumps
801:       nsize: 27
802:       args: -ksp_monitor_short -ksp_converged_reason -stag_ranks_x 3 -stag_ranks_y 3 -stag_ranks_z 3 -pc_fieldsplit_schur_fact_type diag -fieldsplit_0_ksp_type preonly -fieldsplit_1_pc_type none -fieldsplit_0_pc_type lu -fieldsplit_1_ksp_type gmres -fieldsplit_1_ksp_max_it 20

804:    test:
805:       suffix: 2
806:       requires: !single
807:       nsize: 4
808:       args: -ksp_monitor_short -ksp_converged_reason -pc_fieldsplit_schur_fact_type diag -fieldsplit_0_ksp_type preonly -fieldsplit_1_pc_type none -fieldsplit_0_pc_type gamg -fieldsplit_0_mg_levels_ksp_max_it 3 -fieldsplit_1_ksp_type gmres -fieldsplit_1_ksp_max_it 20

810:    test:
811:       suffix: direct_umfpack
812:       requires: suitesparse
813:       nsize: 1
814:       args: -pinpressure 1 -stag_grid_x 5 -stag_grid_y 3 -stag_grid_z 4 -ksp_monitor_short -pc_type lu -pc_factor_mat_solver_type umfpack

816: TEST*/