Actual source code: ex7.c

petsc-3.4.4 2014-03-13
  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 <petscdmda.h>
 42: #include <petscsnes.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

152:   SNESSetComputeInitialGuess(snes, FormInitialGuess,NULL);

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

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

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

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

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

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

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

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

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

236:   SNESGetDM(snes,&da);
237:   DMGetApplicationContext(da,&user);
238:   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);

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

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

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

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

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

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

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

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

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

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

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

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

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

383:   /* Naive Jacobian calculation:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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: }