Actual source code: ex2.c

      1: static char help[] = "Solve a toy 2D 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 2D domain, using a manufactured
      6:   solution.

      8:   u_xx + u_yy - p_x = f^x
      9:   v_xx + v_yy - p_y = f^y
     10:   u_x + v_y         = g

     12:   g is 0 in the physical case.

     14:   Boundary conditions give prescribed flow perpendicular to the boundaries,
     15:   and zero derivative perpendicular to them (free slip).

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

     21:      ./ex2 -pinpressure 1 -pc_type lu -pc_factor_mat_solver_type umfpack

     23:   This example demonstrates the use of DMProduct to efficiently store coordinates
     24:   on an orthogonal grid.

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

     31: /* Shorter, more convenient names for DMStagStencilLocation entries */
     32: #define DOWN_LEFT  DMSTAG_DOWN_LEFT
     33: #define DOWN       DMSTAG_DOWN
     34: #define DOWN_RIGHT DMSTAG_DOWN_RIGHT
     35: #define LEFT       DMSTAG_LEFT
     36: #define ELEMENT    DMSTAG_ELEMENT
     37: #define RIGHT      DMSTAG_RIGHT
     38: #define UP_LEFT    DMSTAG_UP_LEFT
     39: #define UP         DMSTAG_UP
     40: #define UP_RIGHT   DMSTAG_UP_RIGHT

     42: static PetscErrorCode CreateSystem(DM,Mat*,Vec*,PetscBool);
     43: static PetscErrorCode CreateReferenceSolution(DM,Vec*);
     44: static PetscErrorCode AttachNullspace(DM,Mat);
     45: static PetscErrorCode CheckSolution(Vec,Vec);

     47: /* Manufactured solution. Chosen to be higher order than can be solved exactly,
     48: and to have a zero derivative for flow parallel to the boundaries. That is,
     49: d(ux)/dy = 0 at the top and bottom boundaries, and d(uy)/dx = 0 at the right
     50: and left boundaries. */
     51: static PetscScalar uxRef(PetscScalar x,PetscScalar y){return 0.0*x + y*y - 2.0*y*y*y + y*y*y*y;}      /* no x-dependence  */
     52: static PetscScalar uyRef(PetscScalar x,PetscScalar y) {return x*x - 2.0*x*x*x + x*x*x*x + 0.0*y;}      /* no y-dependence  */
     53: static PetscScalar pRef (PetscScalar x,PetscScalar y) {return -1.0*(x-0.5) + -3.0/2.0*y*y + 0.5;}    /* zero integral    */
     54: static PetscScalar fx   (PetscScalar x,PetscScalar y) {return 0.0*x + 2.0 -12.0*y + 12.0*y*y + 1.0;} /* no x-dependence  */
     55: static PetscScalar fy   (PetscScalar x,PetscScalar y) {return 2.0 -12.0*x + 12.0*x*x + 3.0*y;}
     56: static PetscScalar g    (PetscScalar x,PetscScalar y) {return 0.0*x*y;}                              /* identically zero */

     58: int main(int argc,char **argv)
     59: {
     61:   DM             dmSol;
     62:   Vec            sol,solRef,rhs;
     63:   Mat            A;
     64:   KSP            ksp;
     65:   PC             pc;
     66:   PetscBool      pinPressure;

     68:   /* Initialize PETSc and process command line arguments */
     69:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     70:   pinPressure = PETSC_TRUE;
     71:   PetscOptionsGetBool(NULL,NULL,"-pinpressure",&pinPressure,NULL);

     73:   /* Create 2D DMStag for the solution, and set up. */
     74:   {
     75:     const PetscInt dof0 = 0, dof1 = 1,dof2 = 1; /* 1 dof on each edge and element center */
     76:     const PetscInt stencilWidth = 1;
     77:     DMStagCreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,7,9,PETSC_DECIDE,PETSC_DECIDE,dof0,dof1,dof2,DMSTAG_STENCIL_BOX,stencilWidth,NULL,NULL,&dmSol);
     78:     DMSetFromOptions(dmSol);
     79:     DMSetUp(dmSol);
     80:   }

     82:   /* Define uniform coordinates as a product of 1D arrays */
     83:   DMStagSetUniformCoordinatesProduct(dmSol,0.0,1.0,0.0,1.0,0.0,0.0);

     85:   /* Compute (manufactured) reference solution */
     86:   CreateReferenceSolution(dmSol,&solRef);

     88:   /* Assemble system */
     89:   CreateSystem(dmSol,&A,&rhs,pinPressure);

     91:   /* Attach a constant-pressure nullspace to the operator
     92:   (logically, this should be in CreateSystem, but we separate it here to highlight
     93:    the features of DMStag exposed, in particular DMStagMigrateVec()) */
     94:   if (!pinPressure) {
     95:     AttachNullspace(dmSol,A);
     96:   }

     98:   /* Solve, using the default FieldSplit (Approximate Block Factorization) Preconditioner
     99:      This is not intended to be an example of a good solver!  */
    100:   DMCreateGlobalVector(dmSol,&sol);
    101:   KSPCreate(PETSC_COMM_WORLD,&ksp);
    102:   KSPSetType(ksp,KSPFGMRES);
    103:   KSPSetOperators(ksp,A,A);
    104:   KSPGetPC(ksp,&pc);
    105:   PCSetType(pc,PCFIELDSPLIT);
    106:   PCFieldSplitSetDetectSaddlePoint(pc,PETSC_TRUE);
    107:   KSPSetFromOptions(ksp);
    108:   KSPSolve(ksp,rhs,sol);

    110:   /* Check Solution */
    111:   CheckSolution(sol,solRef);

    113:   /* Clean up and finalize PETSc */
    114:   KSPDestroy(&ksp);
    115:   VecDestroy(&sol);
    116:   VecDestroy(&solRef);
    117:   VecDestroy(&rhs);
    118:   MatDestroy(&A);
    119:   DMDestroy(&dmSol);
    120:   PetscFinalize();
    121:   return ierr;
    122: }

    124: /*
    125: Note: this system is not well-scaled! Generally one would adjust the equations
    126:  to try to get matrix entries to be of comparable order, regardless of grid spacing
    127:  or choice of coefficients.
    128: */
    129: static PetscErrorCode CreateSystem(DM dmSol,Mat *pA,Vec *pRhs, PetscBool pinPressure)
    130: {
    132:   PetscInt       N[2];
    133:   PetscInt       ex,ey,startx,starty,nx,ny;
    134:   PetscInt       iprev,icenter,inext;
    135:   Mat            A;
    136:   Vec            rhs;
    137:   PetscReal      hx,hy;
    138:   PetscScalar    **cArrX,**cArrY;

    140:   /* Here, we showcase two different methods for manipulating local vector entries.
    141:      One can use DMStagStencil objects with DMStagVecSetValuesStencil(),
    142:      making sure to call VecAssemble[Begin/End]() after all values are set.
    143:      Alternately, one can use DMStagVecGetArray[Read]() and DMStagVecRestoreArray[Read]().
    144:      The first approach is used to build the rhs, and the second is used to
    145:      obtain coordinate values. Working with the array is almost certainly more efficient,
    146:      but only allows setting local entries, requires understanding which "slot" to use,
    147:      and doesn't correspond as precisely to the matrix assembly process using DMStagStencil objects */

    150:   DMCreateMatrix(dmSol,pA);
    151:   A = *pA;
    152:   DMCreateGlobalVector(dmSol,pRhs);
    153:   rhs = *pRhs;
    154:   DMStagGetCorners(dmSol,&startx,&starty,NULL,&nx,&ny,NULL,NULL,NULL,NULL);
    155:   DMStagGetGlobalSizes(dmSol,&N[0],&N[1],NULL);
    156:   hx = 1.0/N[0]; hy = 1.0/N[1];
    157:   DMStagGetProductCoordinateArraysRead(dmSol,&cArrX,&cArrY,NULL);
    158:   DMStagGetProductCoordinateLocationSlot(dmSol,ELEMENT,&icenter);
    159:   DMStagGetProductCoordinateLocationSlot(dmSol,LEFT,&iprev);
    160:   DMStagGetProductCoordinateLocationSlot(dmSol,RIGHT,&inext);

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

    167:       if (ex == N[0]-1) {
    168:         /* Right Boundary velocity Dirichlet */
    169:         DMStagStencil row;
    170:         PetscScalar   valRhs;
    171:         const PetscScalar valA = 1.0;
    172:         row.i = ex; row.j = ey; row.loc = RIGHT; row.c = 0;
    173:         DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
    174:         valRhs = uxRef(cArrX[ex][inext],cArrY[ey][icenter]);
    175:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    176:       }
    177:       if (ey == N[1]-1) {
    178:         /* Top boundary velocity Dirichlet */
    179:         DMStagStencil row;
    180:         PetscScalar   valRhs;
    181:         const PetscScalar valA = 1.0;
    182:         row.i = ex; row.j = ey; row.loc = UP; row.c = 0;
    183:         DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
    184:         valRhs = uyRef(cArrX[ex][icenter],cArrY[ey][inext]);
    185:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    186:       }

    188:       if (ey == 0) {
    189:         /* Bottom boundary velocity Dirichlet */
    190:         DMStagStencil row;
    191:         PetscScalar   valRhs;
    192:         const PetscScalar valA = 1.0;
    193:         row.i = ex; row.j = ey; row.loc = DOWN; row.c = 0;
    194:         DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
    195:         valRhs = uyRef(cArrX[ex][icenter],cArrY[ey][iprev]);
    196:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    197:       } else {
    198:         /* Y-momentum equation : (u_xx + u_yy) - p_y = f^y */
    199:         DMStagStencil row,col[7];
    200:         PetscScalar   valA[7],valRhs;
    201:         PetscInt      nEntries;

    203:         row.i    = ex  ; row.j    = ey  ; row.loc    = DOWN;    row.c     = 0;
    204:         if (ex == 0) {
    205:           nEntries = 6;
    206:           col[0].i = ex  ; col[0].j = ey  ; col[0].loc = DOWN;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) -2.0 / (hy*hy);
    207:           col[1].i = ex  ; col[1].j = ey-1; col[1].loc = DOWN;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
    208:           col[2].i = ex  ; col[2].j = ey+1; col[2].loc = DOWN;    col[2].c  = 0; valA[2] =  1.0 / (hy*hy);
    209:           /* Missing left element */
    210:           col[3].i = ex+1; col[3].j = ey  ; col[3].loc = DOWN;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
    211:           col[4].i = ex  ; col[4].j = ey-1; col[4].loc = ELEMENT; col[4].c  = 0; valA[4] =  1.0 / hy;
    212:           col[5].i = ex  ; col[5].j = ey  ; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] = -1.0 / hy;
    213:         } else if (ex == N[0]-1) {
    214:           /* Right boundary y velocity stencil */
    215:           nEntries = 6;
    216:           col[0].i = ex  ; col[0].j = ey  ; col[0].loc = DOWN;    col[0].c  = 0; valA[0] = -1.0 / (hx*hx) -2.0 / (hy*hy);
    217:           col[1].i = ex  ; col[1].j = ey-1; col[1].loc = DOWN;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
    218:           col[2].i = ex  ; col[2].j = ey+1; col[2].loc = DOWN;    col[2].c  = 0; valA[2] =  1.0 / (hy*hy);
    219:           col[3].i = ex-1; col[3].j = ey  ; col[3].loc = DOWN;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
    220:           /* Missing right element */
    221:           col[4].i = ex  ; col[4].j = ey-1; col[4].loc = ELEMENT; col[4].c  = 0; valA[4] =  1.0 / hy;
    222:           col[5].i = ex  ; col[5].j = ey  ; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] = -1.0 / hy;
    223:         } else {
    224:           nEntries = 7;
    225:           col[0].i = ex  ; col[0].j = ey  ; col[0].loc = DOWN;    col[0].c  = 0; valA[0] = -2.0 / (hx*hx) -2.0 / (hy*hy);
    226:           col[1].i = ex  ; col[1].j = ey-1; col[1].loc = DOWN;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
    227:           col[2].i = ex  ; col[2].j = ey+1; col[2].loc = DOWN;    col[2].c  = 0; valA[2] =  1.0 / (hy*hy);
    228:           col[3].i = ex-1; col[3].j = ey  ; col[3].loc = DOWN;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
    229:           col[4].i = ex+1; col[4].j = ey  ; col[4].loc = DOWN;    col[4].c  = 0; valA[4] =  1.0 / (hx*hx);
    230:           col[5].i = ex  ; col[5].j = ey-1; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] =  1.0 / hy;
    231:           col[6].i = ex  ; col[6].j = ey  ; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] = -1.0 / hy;
    232:         }
    233:         DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
    234:         valRhs = fy(cArrX[ex][icenter],cArrY[ey][iprev]);
    235:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    236:       }

    238:       if (ex == 0) {
    239:         /* Left velocity Dirichlet */
    240:         DMStagStencil row;
    241:         PetscScalar   valRhs;
    242:         const PetscScalar valA = 1.0;
    243:         row.i = ex; row.j = ey; row.loc = LEFT; row.c = 0;
    244:         DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
    245:         valRhs = uxRef(cArrX[ex][iprev],cArrY[ey][icenter]);
    246:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    247:       } else {
    248:         /* X-momentum equation : (u_xx + u_yy) - p_x = f^x */
    249:         DMStagStencil row,col[7];
    250:         PetscScalar   valA[7],valRhs;
    251:         PetscInt      nEntries;
    252:         row.i    = ex  ; row.j    = ey  ; row.loc    = LEFT;    row.c     = 0;

    254:         if (ey == 0) {
    255:           nEntries = 6;
    256:           col[0].i = ex  ; col[0].j = ey  ; col[0].loc = LEFT;    col[0].c  = 0; valA[0] = -2.0 /(hx*hx) -1.0 /(hy*hy);
    257:           /* missing term from element below */
    258:           col[1].i = ex  ; col[1].j = ey+1; col[1].loc = LEFT;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
    259:           col[2].i = ex-1; col[2].j = ey  ; col[2].loc = LEFT;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
    260:           col[3].i = ex+1; col[3].j = ey  ; col[3].loc = LEFT;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
    261:           col[4].i = ex-1; col[4].j = ey  ; col[4].loc = ELEMENT; col[4].c  = 0; valA[4] =  1.0 / hx;
    262:           col[5].i = ex  ; col[5].j = ey  ; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] = -1.0 / hx;
    263:         } else if (ey == N[1]-1) {
    264:           /* Top boundary x velocity stencil */
    265:           nEntries = 6;
    266:           row.i    = ex  ; row.j    = ey  ; row.loc    = LEFT;    row.c     = 0;
    267:           col[0].i = ex  ; col[0].j = ey  ; col[0].loc = LEFT;    col[0].c  = 0; valA[0] = -2.0 / (hx*hx) -1.0 / (hy*hy);
    268:           col[1].i = ex  ; col[1].j = ey-1; col[1].loc = LEFT;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
    269:           /* Missing element above term */
    270:           col[2].i = ex-1; col[2].j = ey  ; col[2].loc = LEFT;    col[2].c  = 0; valA[2] =  1.0 / (hx*hx);
    271:           col[3].i = ex+1; col[3].j = ey  ; col[3].loc = LEFT;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
    272:           col[4].i = ex-1; col[4].j = ey  ; col[4].loc = ELEMENT; col[4].c  = 0; valA[4] =  1.0 / hx;
    273:           col[5].i = ex  ; col[5].j = ey  ; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] = -1.0 / hx;
    274:         } else {
    275:           /* Note how this is identical to the stencil for U_y, with "DOWN" replaced by "LEFT" and the pressure derivative in the other direction */
    276:           nEntries = 7;
    277:           col[0].i = ex  ; col[0].j = ey  ; col[0].loc = LEFT;    col[0].c  = 0; valA[0] = -2.0 / (hx*hx) + -2.0 / (hy*hy);
    278:           col[1].i = ex  ; col[1].j = ey-1; col[1].loc = LEFT;    col[1].c  = 0; valA[1] =  1.0 / (hy*hy);
    279:           col[2].i = ex  ; col[2].j = ey+1; col[2].loc = LEFT;    col[2].c  = 0; valA[2] =  1.0 / (hy*hy);
    280:           col[3].i = ex-1; col[3].j = ey  ; col[3].loc = LEFT;    col[3].c  = 0; valA[3] =  1.0 / (hx*hx);
    281:           col[4].i = ex+1; col[4].j = ey  ; col[4].loc = LEFT;    col[4].c  = 0; valA[4] =  1.0 / (hx*hx);
    282:           col[5].i = ex-1; col[5].j = ey  ; col[5].loc = ELEMENT; col[5].c  = 0; valA[5] =  1.0 / hx;
    283:           col[6].i = ex  ; col[6].j = ey  ; col[6].loc = ELEMENT; col[6].c  = 0; valA[6] = -1.0 / hx;

    285:         }
    286:         DMStagMatSetValuesStencil(dmSol,A,1,&row,nEntries,col,valA,INSERT_VALUES);
    287:         valRhs = fx(cArrX[ex][iprev],cArrY[ey][icenter]);
    288:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    289:       }

    291:       /* P equation : u_x + v_y = g
    292:          Note that this includes an explicit zero on the diagonal. This is only needed for
    293:          direct solvers (not required if using an iterative solver and setting the constant-pressure nullspace) */
    294:       if (pinPressure && ex == 0 && ey == 0) { /* Pin the first pressure node, if requested */
    295:         DMStagStencil row;
    296:         PetscScalar valA,valRhs;
    297:         row.i = ex; row.j = ey; row.loc  = ELEMENT; row.c = 0;
    298:         valA = 1.0;
    299:         DMStagMatSetValuesStencil(dmSol,A,1,&row,1,&row,&valA,INSERT_VALUES);
    300:         valRhs = pRef(cArrX[ex][icenter],cArrY[ey][icenter]);
    301:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    302:       } else {
    303:         DMStagStencil row,col[5];
    304:         PetscScalar   valA[5],valRhs;

    306:         row.i    = ex; row.j    = ey; row.loc    = ELEMENT; row.c    = 0;
    307:         col[0].i = ex; col[0].j = ey; col[0].loc = LEFT;    col[0].c = 0; valA[0] = -1.0 / hx;
    308:         col[1].i = ex; col[1].j = ey; col[1].loc = RIGHT;   col[1].c = 0; valA[1] =  1.0 / hx;
    309:         col[2].i = ex; col[2].j = ey; col[2].loc = DOWN;    col[2].c = 0; valA[2] = -1.0 / hy;
    310:         col[3].i = ex; col[3].j = ey; col[3].loc = UP;      col[3].c = 0; valA[3] =  1.0 / hy;
    311:         col[4] = row;                                                     valA[4] = 0.0;
    312:         DMStagMatSetValuesStencil(dmSol,A,1,&row,5,col,valA,INSERT_VALUES);
    313:         valRhs = g(cArrX[ex][icenter],cArrY[ey][icenter]);
    314:         DMStagVecSetValuesStencil(dmSol,rhs,1,&row,&valRhs,INSERT_VALUES);
    315:       }
    316:     }
    317:   }
    318:   DMStagRestoreProductCoordinateArraysRead(dmSol,&cArrX,&cArrY,NULL);
    319:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    320:   VecAssemblyBegin(rhs);
    321:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    322:   VecAssemblyEnd(rhs);
    323:   return(0);
    324: }

    326: /* Create a pressure-only DMStag and use it to generate a nullspace vector
    327:    - Create a compatible DMStag with one dof per element (and nothing else).
    328:    - Create a constant vector and normalize it
    329:    - Migrate it to a vector on the original dmSol (making use of the fact
    330:    that this will fill in zeros for "extra" dof)
    331:    - Set the nullspace for the operator
    332:    - Destroy everything (the operator keeps the references it needs) */
    333: static PetscErrorCode AttachNullspace(DM dmSol,Mat A)
    334: {
    336:   DM             dmPressure;
    337:   Vec            constantPressure,basis;
    338:   PetscReal      nrm;
    339:   MatNullSpace   matNullSpace;

    342:   DMStagCreateCompatibleDMStag(dmSol,0,0,1,0,&dmPressure);
    343:   DMGetGlobalVector(dmPressure,&constantPressure);
    344:   VecSet(constantPressure,1.0);
    345:   VecNorm(constantPressure,NORM_2,&nrm);
    346:   VecScale(constantPressure,1.0/nrm);
    347:   DMCreateGlobalVector(dmSol,&basis);
    348:   DMStagMigrateVec(dmPressure,constantPressure,dmSol,basis);
    349:   MatNullSpaceCreate(PetscObjectComm((PetscObject)dmSol),PETSC_FALSE,1,&basis,&matNullSpace);
    350:   VecDestroy(&basis);
    351:   VecDestroy(&constantPressure);
    352:   MatSetNullSpace(A,matNullSpace);
    353:   MatNullSpaceDestroy(&matNullSpace);
    354:   return(0);
    355: }

    357: /* Create a reference solution.
    358:    Here, we use the more direct method of iterating over arrays.  */
    359: static PetscErrorCode CreateReferenceSolution(DM dmSol,Vec *pSolRef)
    360: {
    362:   PetscInt       startx,starty,nx,ny,nExtra[2],ex,ey;
    363:   PetscInt       iuy,iux,ip,iprev,icenter;
    364:   PetscScalar    ***arrSol,**cArrX,**cArrY;
    365:   Vec            solRefLocal;

    368:   DMCreateGlobalVector(dmSol,pSolRef);
    369:   DMGetLocalVector(dmSol,&solRefLocal);

    371:   /* Obtain indices to use in the raw arrays */
    372:   DMStagGetLocationSlot(dmSol,DOWN,0,&iuy);
    373:   DMStagGetLocationSlot(dmSol,LEFT,0,&iux);
    374:   DMStagGetLocationSlot(dmSol,ELEMENT,0,&ip);

    376:   /* Use high-level convenience functions to get raw arrays and indices for 1d coordinates */
    377:   DMStagGetProductCoordinateArraysRead(dmSol,&cArrX,&cArrY,NULL);
    378:   DMStagGetProductCoordinateLocationSlot(dmSol,ELEMENT,&icenter);
    379:   DMStagGetProductCoordinateLocationSlot(dmSol,LEFT,&iprev);

    381:   DMStagVecGetArray(dmSol,solRefLocal,&arrSol);
    382:   DMStagGetCorners(dmSol,&startx,&starty,NULL,&nx,&ny,NULL,&nExtra[0],&nExtra[1],NULL);
    383:   for (ey=starty; ey<starty + ny + nExtra[1]; ++ey) {
    384:     for (ex=startx; ex<startx + nx + nExtra[0]; ++ex) {
    385:       arrSol[ey][ex][iuy] = uyRef(cArrX[ex][icenter],cArrY[ey][iprev]);
    386:       arrSol[ey][ex][iux] = uxRef(cArrX[ex][iprev],cArrY[ey][icenter]);
    387:       if (ey < starty+ny && ex < startx+nx) { /* Don't fill on the dummy elements (though you could, and these values would just be ignored) */
    388:         arrSol[ey][ex][ip]  = pRef(cArrX[ex][icenter],cArrY[ey][icenter]);
    389:       }
    390:     }
    391:   }
    392:   DMStagVecRestoreArray(dmSol,solRefLocal,&arrSol);
    393:   DMStagRestoreProductCoordinateArraysRead(dmSol,&cArrX,&cArrY,NULL);
    394:   DMLocalToGlobal(dmSol,solRefLocal,INSERT_VALUES,*pSolRef);
    395:   DMRestoreLocalVector(dmSol,&solRefLocal);
    396:   return(0);
    397: }

    399: static PetscErrorCode CheckSolution(Vec sol,Vec solRef)
    400: {
    402:   Vec            diff;
    403:   PetscReal      normsolRef,errAbs,errRel;

    406:   VecDuplicate(sol,&diff);
    407:   VecCopy(sol,diff);
    408:   VecAXPY(diff,-1.0,solRef);
    409:   VecNorm(diff,NORM_2,&errAbs);
    410:   VecNorm(solRef,NORM_2,&normsolRef);
    411:   errRel = errAbs/normsolRef;
    412:   PetscPrintf(PETSC_COMM_WORLD,"Error (abs): %g\nError (rel): %g\n",(double)errAbs,(double)errRel);
    413:   VecDestroy(&diff);
    414:   return(0);
    415: }

    417: /*TEST

    419:    test:
    420:       suffix: 1
    421:       nsize: 4
    422:       args: -ksp_monitor_short -ksp_converged_reason

    424:    test:
    425:       suffix: direct_umfpack
    426:       requires: suitesparse
    427:       nsize: 1
    428:       args: -pinpressure 1 -stag_grid_x 8 -stag_grid_y 6 -ksp_monitor_short -pc_type lu -pc_factor_mat_solver_type umfpack

    430: TEST*/