Actual source code: ex7.c

petsc-3.7.3 2016-08-01
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,NULL,"-alpha",&user->alpha,NULL);
130:   PetscOptionsGetReal(NULL,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;
230:   PetscErrorCode         ierr;
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: {
298:   Field          rLocal[3];
299:   PetscScalar    phi[3]    = {0.0, 0.0, 0.0};
300:   PetscReal      xI = i*hx, yI = j*hy, hxhy = hx*hy;
301:   Field          res;
302:   PetscInt       q, k;

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

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

347:   PetscMemzero(&rLocal,sizeof(rLocal));
348:   for (q = 0; q < 4; q++) {
349:     phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
350:     phi[1] = quadPoints[q*2];
351:     phi[2] = quadPoints[q*2+1];
352:     res.u  = quadWeights[q]*PetscExpScalar(u[0].u*phi[0] + u[1].u*phi[1] + u[2].u*phi[2]);
353:     res.v  = quadWeights[q]*PetscExpScalar(u[0].v*phi[0] + u[1].v*phi[1] + u[2].v*phi[2]);

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

373: /*
374:    FormFunctionLocal - Evaluates nonlinear function, F(x).

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

388:   /* Naive Jacobian calculation:

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

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

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

414:        2 (1)    (0)
415:      i,j+1 --- i+1,j+1
416:        |  \      |
417:        |   \     |
418:        |    \    |
419:        |     \   |
420:        |      \  |
421:       i,j  --- i+1,j
422:        0         1 (2)

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

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

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

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

546:         f[j+1][i].u = x[j+1][i].u - uExact.u;
547:         f[j+1][i].v = x[j+1][i].v - uExact.v;
548:         f[j+1][i].p = x[j+1][i].p - uExact.p;
549:       }
550:     }
551:   }

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

565: PetscErrorCode nonlinearJacobian(PetscReal lambda, Field u[], PetscScalar J[])
566: {
568:   return(0);
569: }

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

591:   alpha   = user->alpha;
592:   lambda  = user->lambda;
593:   hx      = 1.0/(PetscReal)(info->mx-1);
594:   hy      = 1.0/(PetscReal)(info->my-1);
595:   sc      = hx*hy*lambda;
596:   hxhy    = hx*hy;
597:   detJInv = hxhy;

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

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

633:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
634:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);

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

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

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

642:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
643: #if defined(PRES_BC)
644:         ident.i = i; ident.j = j; ident.c = 2;

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

676:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
677:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);

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

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

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

685:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
686: #if defined(PRES_BC)
687:         ident.i = i+1; ident.j = j; ident.c = 2;

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

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

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

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

727:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
728: #if defined(PRES_BC)
729:         ident.i = i+1; ident.j = j+1; ident.c = 2;

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

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

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

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

769:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
770: #if defined(PRES_BC)
771:         ident.i = i; ident.j = j+1; ident.c = 2;

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

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

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

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

874:   /*
875:      Tell the matrix we will never add a new nonzero location to the
876:      matrix. If we do, it will generate an error.
877:   */
878:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);

880:   CreateNullSpace(info->da,&N);
881:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&N,&nullspace);
882:   VecDestroy(&N);
883:   MatSetNullSpace(A,nullspace);
884:   MatNullSpaceDestroy(&nullspace);
885:   return(0);
886: }

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

904:   DMDAGetLocalInfo(da, &info);
905:   DMGetLocalVector(da, &fLocalVec);
906:   DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
907:   DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
908:   DMDAVecGetArray(da, fLocalVec, &f);

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

937:         / x_U \ = / -1  0 \ / x_L \ + / hx \
938:         \ y_U /   \  0 -1 / \ y_L /   \ hy /
939:        */
940:       for (q = 0; q < 4; q++) {
941:         phi[0]  = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
942:         phi[1]  = quadPoints[q*2];
943:         phi[2]  = quadPoints[q*2+1];
944:         u.u     = uLocal[2].u*phi[0]+ uLocal[3].u*phi[1] + uLocal[1].u*phi[2];
945:         u.v     = uLocal[2].v*phi[0]+ uLocal[3].v*phi[1] + uLocal[1].v*phi[2];
946:         u.p     = uLocal[0].p*phi[0]+ uLocal[1].p*phi[1] + uLocal[3].p*phi[2];
947:         x       = (1.0 - quadPoints[q*2] + (PetscReal)i)*hx;
948:         y       = (1.0 - quadPoints[q*2+1] + (PetscReal)j)*hy;
949:         ExactSolution(PetscAbsScalar(x), PetscAbsScalar(y), &uExact);
950:         *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)));
951:       }
952:     }
953:   }

955:   DMDAVecRestoreArray(da, fLocalVec, &f);
956:   DMRestoreLocalVector(da, &fLocalVec);
957:   return(0);
958: }