Actual source code: ex7.c

petsc-dev 2014-07-22
Report Typos and Errors
  2: static char help[] = "Solves the Stokes equation in a 2D rectangular\n\
  3: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";

  5: /*T
  6:    Concepts: SNES^parallel Stokes example
  7:    Concepts: DMDA^using distributed arrays;
  8:    Processors: n
  9: T*/

 11: /* ------------------------------------------------------------------------

 13:     The Stokes equation is given by the partial differential equation

 15:         -alpha*Laplacian u - \nabla p = f,  0 < x,y < 1,

 17:           \nabla \cdot u              = 0

 19:     with boundary conditions

 21:              u = 0  for  x = 0, x = 1, y = 0, y = 1.

 23:     A P2/P1 finite element approximation is used to discretize the boundary
 24:     value problem on the two triangles which make up each rectangle in the DMDA
 25:     to obtain a nonlinear system of equations.

 27:   ------------------------------------------------------------------------- */

 29: /*
 30:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 31:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 32:    file automatically includes:
 33:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 34:      petscmat.h - matrices
 35:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 36:      petscviewer.h - viewers               petscpc.h  - preconditioners
 37:      petscksp.h   - linear solvers
 38: */
 39: #include <petscsys.h>
 40: #include <petscbag.h>
 41: #include <petscdm.h>
 42: #include <petscdmda.h>
 43: #include <petscsnes.h>

 45: /*
 46:    User-defined application context - contains data needed by the
 47:    application-provided call-back routines, FormJacobianLocal() and
 48:    FormFunctionLocal().
 49: */
 50: typedef struct {
 51:   PetscReal alpha;          /* parameter controlling linearity */
 52:   PetscReal lambda;         /* parameter controlling nonlinearity */
 53: } AppCtx;

 55: typedef struct {
 56:   PetscScalar u;
 57:   PetscScalar v;
 58:   PetscScalar p;
 59: } Field;

 61: static PetscScalar Kref[36] = { 0.5,  0.5, -0.5,  0,  0, -0.5,
 62:                                 0.5,  0.5, -0.5,  0,  0, -0.5,
 63:                                -0.5, -0.5,  0.5,  0,  0,  0.5,
 64:                                   0,    0,    0,  0,  0,    0,
 65:                                   0,    0,    0,  0,  0,    0,
 66:                                -0.5, -0.5,  0.5,  0,  0,  0.5};

 68: static PetscScalar Gradient[18] = {-0.1666667, -0.1666667, -0.1666667,
 69:                                    -0.1666667, -0.1666667, -0.1666667,
 70:                                     0.1666667,  0.1666667,  0.1666667,
 71:                                     0.0,        0.0,        0.0,
 72:                                     0.0,        0.0,        0.0,
 73:                                     0.1666667,  0.1666667,  0.1666667};

 75: static PetscScalar Divergence[18] = {-0.1666667, 0.1666667, 0.0,
 76:                                      -0.1666667, 0.0,       0.1666667,
 77:                                      -0.1666667, 0.1666667, 0.0,
 78:                                      -0.1666667, 0.0,       0.1666667,
 79:                                      -0.1666667, 0.1666667, 0.0,
 80:                                      -0.1666667, 0.0,       0.1666667};

 82: /* These are */
 83: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
 84:                                     0.07503111, 0.64494897,
 85:                                     0.66639025, 0.15505103,
 86:                                     0.28001992, 0.64494897};
 87: static PetscScalar quadWeights[4] = {0.15902069,  0.09097931,  0.15902069,  0.09097931};

 89: /*
 90:    User-defined routines
 91: */
 92: extern PetscErrorCode CreateNullSpace(DM, Vec*);
 93: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
 94: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,AppCtx*);
 95: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,Field**,Mat,Mat,AppCtx*);
 96: extern PetscErrorCode L_2Error(DM, Vec, PetscReal*, AppCtx*);

100: int main(int argc,char **argv)
101: {
102:   DM                  da;
103:   SNES                snes;                    /* nonlinear solver */
104:   AppCtx              *user;                   /* user-defined work context */
105:   PetscBag            bag;
106:   PetscInt            its;                     /* iterations for convergence */
107:   PetscMPIInt         size;
108:   SNESConvergedReason reason;
109:   PetscErrorCode      ierr;
110:   PetscReal           lambda_max = 6.81, lambda_min = 0.0, error;
111:   Vec                 x;

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Initialize program
115:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116:   PetscInitialize(&argc,&argv,(char*)0,help);
117:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
118:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example only works for one process.");

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Initialize problem parameters
122:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123:   PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
124:   PetscBagGetData(bag, (void**) &user);
125:   PetscBagSetName(bag, "params", "Parameters for SNES example 4");
126:   PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
127:   PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
128:   PetscBagSetFromOptions(bag);
129:   PetscOptionsGetReal(NULL,"-alpha",&user->alpha,NULL);
130:   PetscOptionsGetReal(NULL,"-lambda",&user->lambda,NULL);
131:   if (user->lambda > lambda_max || user->lambda < lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda %g is out of range [%g, %g]", (double)user->lambda, (double)lambda_min, (double)lambda_max);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:      Create multilevel DM data structure (SNES) to manage hierarchical solvers
135:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136:   SNESCreate(PETSC_COMM_WORLD,&snes);

138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139:      Create distributed array (DMDA) to manage parallel grid and vectors
140:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
142:   DMDASetFieldName(da, 0, "ooblek");
143:   DMSetApplicationContext(da,user);
144:   SNESSetDM(snes, (DM) da);

146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147:      Set the discretization functions
148:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,user);
150:   DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,void*))FormJacobianLocal,user);
151:   SNESSetFromOptions(snes);

153:   SNESSetComputeInitialGuess(snes, FormInitialGuess,NULL);

155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156:      Solve nonlinear system
157:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158:   SNESSolve(snes,NULL,NULL);
159:   SNESGetIterationNumber(snes,&its);
160:   SNESGetConvergedReason(snes, &reason);

162:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
164:   DMDestroy(&da);
165:   SNESGetDM(snes,&da);
166:   SNESGetSolution(snes,&x);
167:   L_2Error(da, x, &error, user);
168:   PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %g\n", (double)error);

170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171:      Free work space.  All PETSc objects should be destroyed when they
172:      are no longer needed.
173:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174:   SNESDestroy(&snes);
175:   PetscBagDestroy(&bag);
176:   PetscFinalize();
177:   return(0);
178: }

182: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, Field *u)
183: {
185:   (*u).u = x;
186:   (*u).v = -y;
187:   (*u).p = 0.0;
188:   return(0);
189: }

193: PetscErrorCode CreateNullSpace(DM da, Vec *N)
194: {
195:   Field          **x;
196:   PetscInt       xs,ys,xm,ym,i,j;

200:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
201:   DMCreateGlobalVector(da,N);
202:   DMDAVecGetArray(da, *N, &x);
203:   for (j = ys; j < ys+ym; j++) {
204:     for (i = xs; i < xs+xm; i++) {
205:       x[j][i].u = 0.0;
206:       x[j][i].v = 0.0;
207:       x[j][i].p = 1.0;
208:     }
209:   }
210:   DMDAVecRestoreArray(da, *N, &x);
211:   return(0);
212: }

216: /*
217:    FormInitialGuess - Forms initial approximation.

219:    Input Parameters:
220:    dm - The DM context
221:    X - vector

223:    Output Parameter:
224:    X - vector
225: */
226: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
227: {
228:   AppCtx         *user;
229:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
231:   PetscReal      lambda,hx,hy;
232:   PETSC_UNUSED PetscReal temp1;
233:   Field          **x;
234:   DM             da;

237:   SNESGetDM(snes,&da);
238:   DMGetApplicationContext(da,&user);
239:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

241:   lambda = user->lambda;
242:   hx     = 1.0/(PetscReal)(Mx-1);
243:   hy     = 1.0/(PetscReal)(My-1);
244:   if (lambda == 0.0) temp1 = 0.0;
245:   else temp1 = lambda/(lambda + 1.0);

247:   /*
248:      Get a pointer to vector data.
249:        - For default PETSc vectors, VecGetArray() returns a pointer to
250:          the data array.  Otherwise, the routine is implementation dependent.
251:        - You MUST call VecRestoreArray() when you no longer need access to
252:          the array.
253:   */
254:   DMDAVecGetArray(da,X,&x);

256:   /*
257:      Get local grid boundaries (for 2-dimensional DMDA):
258:        xs, ys   - starting grid indices (no ghost points)
259:        xm, ym   - widths of local grid (no ghost points)

261:   */
262:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

264:   /*
265:      Compute initial guess over the locally owned part of the grid
266:   */
267:   for (j=ys; j<ys+ym; j++) {
268:     for (i=xs; i<xs+xm; i++) {
269: #define CHECK_SOLUTION
270: #if defined(CHECK_SOLUTION)
271:       ExactSolution(i*hx, j*hy, &x[j][i]);
272: #else
273:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
274:         /* Boundary conditions are usually zero Dirichlet */
275:         ExactSolution(i*hx, j*hy, &x[j][i]);
276:       } else {
277:         PetscReal temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
278:         x[j][i].u = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
279:         x[j][i].v = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
280:         x[j][i].p = 1.0;
281:       }
282: #endif
283:     }
284:   }

286:   /*
287:      Restore vector
288:   */
289:   DMDAVecRestoreArray(da,X,&x);
290:   return(0);
291: }

295: PetscErrorCode constantResidual(PetscReal lambda, PetscBool isLower, int i, int j, PetscReal hx, PetscReal hy, Field r[])
296: {
297:   Field       rLocal[3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
298:   PetscScalar phi[3]    = {0.0, 0.0, 0.0};
299:   PetscReal   xI = i*hx, yI = j*hy, hxhy = hx*hy;
300:   Field       res;
301:   PetscInt    q, k;

304:   for (q = 0; q < 4; q++) {
305:     PETSC_UNUSED PetscReal x, y;
306:     phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
307:     phi[1] = quadPoints[q*2];
308:     phi[2] = quadPoints[q*2+1];
309:     if (isLower) {
310:       x = xI + PetscAbsScalar(quadPoints[q*2])*hx;
311:       y = yI + PetscAbsScalar(quadPoints[q*2+1])*hy;
312:     } else {
313:       x = xI + 1.0 - PetscAbsScalar(quadPoints[q*2])*hx;
314:       y = yI + 1.0 - PetscAbsScalar(quadPoints[q*2+1])*hy;
315:     }
316:     res.u = quadWeights[q]*(0.0);
317:     res.v = quadWeights[q]*(0.0);
318:     res.p = quadWeights[q]*(0.0);
319:     for (k = 0; k < 3; k++) {
320:       rLocal[k].u += phi[k]*res.u;
321:       rLocal[k].v += phi[k]*res.v;
322:       rLocal[k].p += phi[k]*res.p;
323:     }
324:   }
325:   for (k = 0; k < 3; k++) {
326:     /*printf("  constLocal[%d] = (%g, %g, %g)\n", k, lambda*hxhy*rLocal[k].u, lambda*hxhy*rLocal[k].v, hxhy*rLocal[k].p); */
327:     r[k].u += lambda*hxhy*rLocal[k].u;
328:     r[k].v += lambda*hxhy*rLocal[k].v;
329:     r[k].p += hxhy*rLocal[k].p;
330:   }
331:   return(0);
332: }

336: PetscErrorCode nonlinearResidual(PetscReal lambda, Field u[], Field r[])
337: {
338:   Field       rLocal[3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
339:   PetscScalar phi[3]    = {0.0, 0.0, 0.0};
340:   Field       res;
341:   PetscInt    q;

344:   for (q = 0; q < 4; q++) {
345:     phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
346:     phi[1] = quadPoints[q*2];
347:     phi[2] = quadPoints[q*2+1];
348:     res.u  = quadWeights[q]*PetscExpScalar(u[0].u*phi[0] + u[1].u*phi[1] + u[2].u*phi[2]);
349:     res.v  = quadWeights[q]*PetscExpScalar(u[0].v*phi[0] + u[1].v*phi[1] + u[2].v*phi[2]);

351:     rLocal[0].u += phi[0]*res.u;
352:     rLocal[0].v += phi[0]*res.v;
353:     rLocal[1].u += phi[1]*res.u;
354:     rLocal[1].v += phi[1]*res.v;
355:     rLocal[2].u += phi[2]*res.u;
356:     rLocal[2].v += phi[2]*res.v;
357:   }
358:   r[0].u += lambda*rLocal[0].u;
359:   r[0].v += lambda*rLocal[0].v;
360:   r[1].u += lambda*rLocal[1].u;
361:   r[1].v += lambda*rLocal[1].v;
362:   r[2].u += lambda*rLocal[2].u;
363:   r[2].v += lambda*rLocal[2].v;
364:   return(0);
365: }

369: /*
370:    FormFunctionLocal - Evaluates nonlinear function, F(x).

372:  */
373: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, AppCtx *user)
374: {
375:   Field          uLocal[3];
376:   Field          rLocal[3];
377:   PetscScalar    G[4];
378:   Field          uExact;
379:   PetscReal      alpha,lambda,hx,hy,hxhy,sc,detJInv;
380:   PetscInt       i,j,k,l;

384:   /* Naive Jacobian calculation:

386:      J = / 1/hx  0   \ J^{-1} = / hx   0 \  1/|J| = hx*hy = |J^{-1}|
387:          \  0   1/hy /          \  0  hy /
388:    */
389:   alpha   = user->alpha;
390:   lambda  = user->lambda;
391:   hx      = 1.0/(PetscReal)(info->mx-1);
392:   hy      = 1.0/(PetscReal)(info->my-1);
393:   sc      = hx*hy*lambda;
394:   hxhy    = hx*hy;
395:   detJInv = hxhy;

397:   G[0] = (1.0/(hx*hx)) * detJInv;
398:   G[1] = 0.0;
399:   G[2] = G[1];
400:   G[3] = (1.0/(hy*hy)) * detJInv;
401:   for (k = 0; k < 4; k++) {
402:     /* printf("G[%d] = %g\n", k, G[k]);*/
403:   }

405:   /* Zero the vector */
406:   PetscMemzero((void*) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(Field));
407:   /* Compute function over the locally owned part of the grid. For each
408:      vertex (i,j), we consider the element below:

410:        2 (1)    (0)
411:      i,j+1 --- i+1,j+1
412:        |  \      |
413:        |   \     |
414:        |    \    |
415:        |     \   |
416:        |      \  |
417:       i,j  --- i+1,j
418:        0         1 (2)

420:      and therefore we do not loop over the last vertex in each dimension.
421:   */
422:   for (j = info->ys; j < info->ys+info->ym-1; j++) {
423:     for (i = info->xs; i < info->xs+info->xm-1; i++) {
424:       /* Lower element */
425:       uLocal[0] = x[j][i];
426:       uLocal[1] = x[j][i+1];
427:       uLocal[2] = x[j+1][i];
428:       /* printf("Lower Solution ElementVector for (%d, %d)\n", i, j);*/
429:       for (k = 0; k < 3; k++) {
430:         /* printf("  uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
431:       }
432:       for (k = 0; k < 3; k++) {
433:         rLocal[k].u = 0.0;
434:         rLocal[k].v = 0.0;
435:         rLocal[k].p = 0.0;
436:         for (l = 0; l < 3; l++) {
437:           rLocal[k].u += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].u;
438:           rLocal[k].v += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].v;
439:           /* rLocal[k].u += hxhy*Identity[k*3+l]*uLocal[l].u; */
440:           /* rLocal[k].p += hxhy*Identity[k*3+l]*uLocal[l].p; */
441:           /* Gradient */
442:           rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
443:           rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
444:           /* Divergence */
445:           rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
446:         }
447:       }
448:       /* printf("Lower Laplacian ElementVector for (%d, %d)\n", i, j);*/
449:       for (k = 0; k < 3; k++) {
450:         /* printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
451:       }
452:       constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
453:       /* printf("Lower Laplacian+Constant ElementVector for (%d, %d)\n", i, j);*/
454:       for (k = 0; k < 3; k++) {
455:         /* printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
456:       }
457:       nonlinearResidual(0.0*sc, uLocal, rLocal);
458:       /* printf("Lower Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
459:       for (k = 0; k < 3; k++) {
460:         /* printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
461:       }
462:       f[j][i].u   += rLocal[0].u;
463:       f[j][i].v   += rLocal[0].v;
464:       f[j][i].p   += rLocal[0].p;
465:       f[j][i+1].u += rLocal[1].u;
466:       f[j][i+1].v += rLocal[1].v;
467:       f[j][i+1].p += rLocal[1].p;
468:       f[j+1][i].u += rLocal[2].u;
469:       f[j+1][i].v += rLocal[2].v;
470:       f[j+1][i].p += rLocal[2].p;
471:       /* Upper element */
472:       uLocal[0] = x[j+1][i+1];
473:       uLocal[1] = x[j+1][i];
474:       uLocal[2] = x[j][i+1];
475:       /* printf("Upper Solution ElementVector for (%d, %d)\n", i, j);*/
476:       for (k = 0; k < 3; k++) {
477:         /* printf("  uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);*/
478:       }
479:       for (k = 0; k < 3; k++) {
480:         rLocal[k].u = 0.0;
481:         rLocal[k].v = 0.0;
482:         rLocal[k].p = 0.0;
483:         for (l = 0; l < 3; l++) {
484:           rLocal[k].u += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].u;
485:           rLocal[k].v += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l].v;
486:           /* rLocal[k].p += Identity[k*3+l]*uLocal[l].p; */
487:           /* Gradient */
488:           rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
489:           rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
490:           /* Divergence */
491:           rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
492:         }
493:       }
494:       /* printf("Upper Laplacian ElementVector for (%d, %d)\n", i, j);*/
495:       for (k = 0; k < 3; k++) {
496:         /* printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
497:       }
498:       constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
499:       /* printf("Upper Laplacian+Constant ElementVector for (%d, %d)\n", i, j);*/
500:       for (k = 0; k < 3; k++) {
501:         /* printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
502:       }
503:       nonlinearResidual(0.0*sc, uLocal, rLocal);
504:       /* printf("Upper Full nonlinear ElementVector for (%d, %d)\n", i, j);*/
505:       for (k = 0; k < 3; k++) {
506:         /* printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);*/
507:       }
508:       f[j+1][i+1].u += rLocal[0].u;
509:       f[j+1][i+1].v += rLocal[0].v;
510:       f[j+1][i+1].p += rLocal[0].p;
511:       f[j+1][i].u   += rLocal[1].u;
512:       f[j+1][i].v   += rLocal[1].v;
513:       f[j+1][i].p   += rLocal[1].p;
514:       f[j][i+1].u   += rLocal[2].u;
515:       f[j][i+1].v   += rLocal[2].v;
516:       f[j][i+1].p   += rLocal[2].p;
517:       /* Boundary conditions */
518:       if (i == 0 || j == 0) {
519:         ExactSolution(i*hx, j*hy, &uExact);

521:         f[j][i].u = x[j][i].u - uExact.u;
522:         f[j][i].v = x[j][i].v - uExact.v;
523:         f[j][i].p = x[j][i].p - uExact.p;
524:       }
525:       if ((i == info->mx-2) || (j == 0)) {
526:         ExactSolution((i+1)*hx, j*hy, &uExact);

528:         f[j][i+1].u = x[j][i+1].u - uExact.u;
529:         f[j][i+1].v = x[j][i+1].v - uExact.v;
530:         f[j][i+1].p = x[j][i+1].p - uExact.p;
531:       }
532:       if ((i == info->mx-2) || (j == info->my-2)) {
533:         ExactSolution((i+1)*hx, (j+1)*hy, &uExact);

535:         f[j+1][i+1].u = x[j+1][i+1].u - uExact.u;
536:         f[j+1][i+1].v = x[j+1][i+1].v - uExact.v;
537:         f[j+1][i+1].p = x[j+1][i+1].p - uExact.p;
538:       }
539:       if ((i == 0) || (j == info->my-2)) {
540:         ExactSolution(i*hx, (j+1)*hy, &uExact);

542:         f[j+1][i].u = x[j+1][i].u - uExact.u;
543:         f[j+1][i].v = x[j+1][i].v - uExact.v;
544:         f[j+1][i].p = x[j+1][i].p - uExact.p;
545:       }
546:     }
547:   }

549:   for (j = info->ys+info->ym-1; j >= info->ys; j--) {
550:     for (i = info->xs; i < info->xs+info->xm; i++) {
551:       /* printf("f[%d][%d] = (%g, %g, %g) ", j, i, f[j][i].u, f[j][i].v, f[j][i].p);*/
552:     }
553:     /*printf("\n");*/
554:   }
555:   PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
556:   return(0);
557: }

561: PetscErrorCode nonlinearJacobian(PetscReal lambda, Field u[], PetscScalar J[])
562: {
564:   return(0);
565: }

569: /*
570:    FormJacobianLocal - Evaluates Jacobian matrix.
571: */
572: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field **x, Mat A,Mat jac, AppCtx *user)
573: {
574:   Field          uLocal[4];
575:   PetscScalar    JLocal[144];
576:   MatStencil     rows[4*3], cols[4*3], ident;
577:   PetscInt       lowerRow[3] = {0, 1, 3};
578:   PetscInt       upperRow[3] = {2, 3, 1};
579:   PetscInt       hasLower[3], hasUpper[3], velocityRows[4], pressureRows[4];
580:   PetscScalar    alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
581:   PetscInt       i,j,k,l,numRows,dof = info->dof;
583:   MatNullSpace   nullspace;
584:   Vec            N;

587:   alpha   = user->alpha;
588:   lambda  = user->lambda;
589:   hx      = 1.0/(PetscReal)(info->mx-1);
590:   hy      = 1.0/(PetscReal)(info->my-1);
591:   sc      = hx*hy*lambda;
592:   hxhy    = hx*hy;
593:   detJInv = hxhy;

595:   G[0] = (1.0/(hx*hx)) * detJInv;
596:   G[1] = 0.0;
597:   G[2] = G[1];
598:   G[3] = (1.0/(hy*hy)) * detJInv;
599:   for (k = 0; k < 4; k++) {
600:     /* printf("G[%d] = %g\n", k, G[k]);*/
601:   }

603:   MatZeroEntries(jac);
604:   /*
605:      Compute entries for the locally owned part of the Jacobian.
606:       - Currently, all PETSc parallel matrix formats are partitioned by
607:         contiguous chunks of rows across the processors.
608:       - Each processor needs to insert only elements that it owns
609:         locally (but any non-local elements will be sent to the
610:         appropriate processor during matrix assembly).
611:       - Here, we set all entries for a particular row at once.
612:       - We can set matrix entries either using either
613:         MatSetValuesLocal() or MatSetValues(), as discussed above.
614:   */
615: #define NOT_PRES_BC 1
616:   for (j=info->ys; j<info->ys+info->ym-1; j++) {
617:     for (i=info->xs; i<info->xs+info->xm-1; i++) {
618:       PetscMemzero(JLocal, 144 * sizeof(PetscScalar));
619:       numRows = 0;
620:       /* Lower element */
621:       uLocal[0] = x[j][i];
622:       uLocal[1] = x[j][i+1];
623:       uLocal[2] = x[j+1][i+1];
624:       uLocal[3] = x[j+1][i];
625:       /* i,j */
626:       if (i == 0 || j == 0) {
627:         hasLower[0] = 0;

629:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
630:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);

632:         ident.i = i; ident.j = j; ident.c = 0;

634:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);

636:         ident.i = i; ident.j = j; ident.c = 1;

638:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
639: #if defined(PRES_BC)
640:         ident.i = i; ident.j = j; ident.c = 2;

642:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
643: #endif
644:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
645:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
646:       } else {
647:         hasLower[0]     = 1;
648:         velocityRows[0] = numRows;
649:         rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 0;
650:         numRows++;
651:         rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 1;
652:         numRows++;
653: #if defined(PRES_BC)
654:         pressureRows[0] = numRows;
655:         rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
656:         numRows++;
657: #endif
658:       }
659: #if !defined(PRES_BC)
660:       pressureRows[0] = numRows;
661:       rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
662:       numRows++;
663: #endif
664:       cols[0*dof+0].i = i; cols[0*dof+0].j = j; cols[0*dof+0].c = 0;
665:       cols[0*dof+1].i = i; cols[0*dof+1].j = j; cols[0*dof+1].c = 1;
666:       cols[0*dof+2].i = i; cols[0*dof+2].j = j; cols[0*dof+2].c = 2;
667:       /* i+1,j */
668:       if ((i == info->mx-2) || (j == 0)) {
669:         hasLower[1] = 0;
670:         hasUpper[2] = 0;

672:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
673:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);

675:         ident.i = i+1; ident.j = j; ident.c = 0;

677:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);

679:         ident.i = i+1; ident.j = j; ident.c = 1;

681:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
682: #if defined(PRES_BC)
683:         ident.i = i+1; ident.j = j; ident.c = 2;

685:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
686: #endif
687:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
688:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
689:       } else {
690:         hasLower[1]     = 1;
691:         hasUpper[2]     = 1;
692:         velocityRows[1] = numRows;
693:         rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 0;
694:         numRows++;
695:         rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 1;
696:         numRows++;
697: #if defined(PRES_BC)
698:         pressureRows[1] = numRows;
699:         rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
700:         numRows++;
701: #endif
702:       }
703: #if !defined(PRES_BC)
704:       pressureRows[1] = numRows;
705:       rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
706:       numRows++;
707: #endif
708:       cols[1*dof+0].i = i+1; cols[1*dof+0].j = j; cols[1*dof+0].c = 0;
709:       cols[1*dof+1].i = i+1; cols[1*dof+1].j = j; cols[1*dof+1].c = 1;
710:       cols[1*dof+2].i = i+1; cols[1*dof+2].j = j; cols[1*dof+2].c = 2;
711:       /* i+1,j+1 */
712:       if ((i == info->mx-2) || (j == info->my-2)) {
713:         hasUpper[0] = 0;
714:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
715:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);

717:         ident.i = i+1; ident.j = j+1; ident.c = 0;

719:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);

721:         ident.i = i+1; ident.j = j+1; ident.c = 1;

723:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
724: #if defined(PRES_BC)
725:         ident.i = i+1; ident.j = j+1; ident.c = 2;

727:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
728: #endif
729:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
730:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
731:       } else {
732:         hasUpper[0]     = 1;
733:         velocityRows[2] = numRows;
734:         rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 0;
735:         numRows++;
736:         rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 1;
737:         numRows++;
738: #if defined(PRES_BC)
739:         pressureRows[2] = numRows;
740:         rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
741:         numRows++;
742: #endif
743:       }
744: #if !defined(PRES_BC)
745:       pressureRows[2] = numRows;
746:       rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
747:       numRows++;
748: #endif
749:       cols[2*dof+0].i = i+1; cols[2*dof+0].j = j+1; cols[2*dof+0].c = 0;
750:       cols[2*dof+1].i = i+1; cols[2*dof+1].j = j+1; cols[2*dof+1].c = 1;
751:       cols[2*dof+2].i = i+1; cols[2*dof+2].j = j+1; cols[2*dof+2].c = 2;
752:       /* i,j+1 */
753:       if ((i == 0) || (j == info->my-2)) {
754:         hasLower[2] = 0;
755:         hasUpper[1] = 0;
756:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
757:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);

759:         ident.i = i; ident.j = j+1; ident.c = 0;

761:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);

763:         ident.i = i; ident.j = j+1; ident.c = 1;

765:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
766: #if defined(PRES_BC)
767:         ident.i = i; ident.j = j+1; ident.c = 2;

769:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
770: #endif
771:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
772:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
773:       } else {
774:         hasLower[2]     = 1;
775:         hasUpper[1]     = 1;
776:         velocityRows[3] = numRows;
777:         rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 0;
778:         numRows++;
779:         rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 1;
780:         numRows++;
781: #if defined(PRES_BC)
782:         pressureRows[3] = numRows;
783:         rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
784:         numRows++;
785: #endif
786:       }
787: #if !defined(PRES_BC)
788:       pressureRows[3] = numRows;
789:       rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
790:       numRows++;
791: #endif
792:       cols[3*dof+0].i = i; cols[3*dof+0].j = j+1; cols[3*dof+0].c = 0;
793:       cols[3*dof+1].i = i; cols[3*dof+1].j = j+1; cols[3*dof+1].c = 1;
794:       cols[3*dof+2].i = i; cols[3*dof+2].j = j+1; cols[3*dof+2].c = 2;

796:       /* Lower Element */
797:       for (k = 0; k < 3; k++) {
798: #if defined(PRES_BC)
799:         if (!hasLower[k]) continue;
800: #endif
801:         for (l = 0; l < 3; l++) {
802:           /* Divergence */
803:           JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
804:           JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
805: /*        JLocal[pressureRows[lowerRow[k]]*dof*4 + lowerRow[l]*dof+2] += Identity[k*3 + l]; */
806:         }
807:         if (!hasLower[k]) continue;
808:         for (l = 0; l < 3; l++) {
809:           /* Laplacian */
810:           JLocal[velocityRows[lowerRow[k]]*dof*4     + lowerRow[l]*dof+0] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
811:           JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+1] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
812: /*        JLocal[velocityRows[lowerRow[k]]*dof*4     + lowerRow[l]*dof+0] += Identity[k*3 + l]; */
813: /*        JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+1] += Identity[k*3 + l]; */
814:           /* Gradient */
815:           JLocal[velocityRows[lowerRow[k]]*dof*4     + lowerRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
816:           JLocal[(velocityRows[lowerRow[k]]+1)*dof*4 + lowerRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
817:         }
818:       }
819:       /* Upper Element */
820:       for (k = 0; k < 3; k++) {
821: #if defined(PRES_BC)
822:         if (!hasUpper[k]) continue;
823: #endif
824:         for (l = 0; l < 3; l++) {
825:           /* Divergence */
826:           JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+0] += hx*Divergence[(k*2+0)*3 + l];
827:           JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+1] += hy*Divergence[(k*2+1)*3 + l];
828: /*        JLocal[pressureRows[upperRow[k]]*dof*4 + upperRow[l]*dof+2] += Identity[k*3 + l]; */
829:         }
830:         if (!hasUpper[k]) continue;
831:         for (l = 0; l < 3; l++) {
832:           /* Laplacian */
833:           JLocal[velocityRows[upperRow[k]]*dof*4     + upperRow[l]*dof+0] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
834:           JLocal[(velocityRows[upperRow[k]]+1)*dof*4 + upperRow[l]*dof+1] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
835:           /* Gradient */
836:           JLocal[velocityRows[upperRow[k]]*dof*4     + upperRow[l]*dof+2] += hx*Gradient[(k*2+0)*3 + l];
837:           JLocal[(velocityRows[upperRow[k]]+1)*dof*4 + upperRow[l]*dof+2] += hy*Gradient[(k*2+1)*3 + l];
838:         }
839:       }

841:       nonlinearJacobian(-1.0*PetscAbsScalar(sc), uLocal, JLocal);
842:       /* printf("Element matrix for (%d, %d)\n", i, j);*/
843:       /* printf("   col         ");*/
844:       for (l = 0; l < 4*3; l++) {
845:         /* printf("(%d, %d, %d) ", cols[l].i, cols[l].j, cols[l].c);*/
846:       }
847:       /* printf("\n");*/
848:       for (k = 0; k < numRows; k++) {
849:         /* printf("row (%d, %d, %d): ", rows[k].i, rows[k].j, rows[k].c);*/
850:         for (l = 0; l < 4; l++) {
851:           /* printf("%9.6g %9.6g %9.6g ", JLocal[k*dof*4 + l*dof+0], JLocal[k*dof*4 + l*dof+1], JLocal[k*dof*4 + l*dof+2]);*/
852:         }
853:         /* printf("\n");*/
854:       }
855:       MatSetValuesStencil(jac,numRows,rows,4*dof,cols, JLocal,ADD_VALUES);
856:     }
857:   }

859:   /*
860:      Assemble matrix, using the 2-step process:
861:        MatAssemblyBegin(), MatAssemblyEnd().
862:   */
863:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
864:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
865:   if (A != jac) {
866:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
867:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
868:   }

870:   /*
871:      Tell the matrix we will never add a new nonzero location to the
872:      matrix. If we do, it will generate an error.
873:   */
874:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);

876:   CreateNullSpace(info->da,&N);
877:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&N,&nullspace);
878:   VecDestroy(&N);
879:   MatSetNullSpace(jac,nullspace);
880:   MatNullSpaceDestroy(&nullspace);
881:   return(0);
882: }

886: /*
887:   L_2Error - Integrate the L_2 error of our solution over each face
888: */
889: PetscErrorCode L_2Error(DM da, Vec fVec, PetscReal *error, AppCtx *user)
890: {
891:   DMDALocalInfo  info;
892:   Vec            fLocalVec;
893:   Field          **f;
894:   Field          u, uExact, uLocal[4];
895:   PetscScalar    hx, hy, hxhy, x, y, phi[3];
896:   PetscInt       i, j, q;

900:   DMDAGetLocalInfo(da, &info);
901:   DMGetLocalVector(da, &fLocalVec);
902:   DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
903:   DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
904:   DMDAVecGetArray(da, fLocalVec, &f);

906:   *error = 0.0;
907:   hx     = 1.0/(PetscReal)(info.mx-1);
908:   hy     = 1.0/(PetscReal)(info.my-1);
909:   hxhy   = hx*hy;
910:   for (j = info.ys; j < info.ys+info.ym-1; j++) {
911:     for (i = info.xs; i < info.xs+info.xm-1; i++) {
912:       uLocal[0] = f[j][i];
913:       uLocal[1] = f[j][i+1];
914:       uLocal[2] = f[j+1][i+1];
915:       uLocal[3] = f[j+1][i];
916:       /* Lower element */
917:       for (q = 0; q < 4; q++) {
918:         phi[0]  = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
919:         phi[1]  = quadPoints[q*2];
920:         phi[2]  = quadPoints[q*2+1];
921:         u.u     = uLocal[0].u*phi[0]+ uLocal[1].u*phi[1] + uLocal[3].u*phi[2];
922:         u.v     = uLocal[0].v*phi[0]+ uLocal[1].v*phi[1] + uLocal[3].v*phi[2];
923:         u.p     = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
924:         x       = (quadPoints[q*2] + (PetscReal)i)*hx;
925:         y       = (quadPoints[q*2+1] + (PetscReal)j)*hy;
926:         ExactSolution(PetscAbsScalar(x), PetscAbsScalar(y), &uExact);
927:         *error += PetscAbsScalar(hxhy*quadWeights[q]*((u.u - uExact.u)*(u.u - uExact.u) + (u.v - uExact.v)*(u.v - uExact.v) + (u.p - uExact.p)*(u.p - uExact.p)));
928:       }
929:       /* Upper element */
930:       /*
931:         The affine map from the lower to the upper is

933:         / x_U \ = / -1  0 \ / x_L \ + / hx \
934:         \ y_U /   \  0 -1 / \ y_L /   \ hy /
935:        */
936:       for (q = 0; q < 4; q++) {
937:         phi[0]  = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
938:         phi[1]  = quadPoints[q*2];
939:         phi[2]  = quadPoints[q*2+1];
940:         u.u     = uLocal[2].u*phi[0]+ uLocal[3].u*phi[1] + uLocal[1].u*phi[2];
941:         u.v     = uLocal[2].v*phi[0]+ uLocal[3].v*phi[1] + uLocal[1].v*phi[2];
942:         u.p     = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
943:         x       = (1.0 - quadPoints[q*2] + (PetscReal)i)*hx;
944:         y       = (1.0 - quadPoints[q*2+1] + (PetscReal)j)*hy;
945:         ExactSolution(PetscAbsScalar(x), PetscAbsScalar(y), &uExact);
946:         *error += PetscAbsScalar(hxhy*quadWeights[q]*((u.u - uExact.u)*(u.u - uExact.u) + (u.v - uExact.v)*(u.v - uExact.v) + (u.p - uExact.p)*(u.p - uExact.p)));
947:       }
948:     }
949:   }

951:   DMDAVecRestoreArray(da, fLocalVec, &f);
952:   /* DMLocalToGlobalBegin(da,xLocalVec,ADD_VALUES,xVec); */
953:   /* DMLocalToGlobalEnd(da,xLocalVec,ADD_VALUES,xVec); */
954:   DMRestoreLocalVector(da, &fLocalVec);
955:   return(0);
956: }