Actual source code: ex7.c

  1: /* Program usage:  mpiexec -n <procs> ex7 [-help] [all PETSc options] */

  3: static char help[] = "Nonlinear PDE in 2d.\n\
  4: We solve the Stokes equation in a 2D rectangular\n\
  5: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";

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

 13: /* ------------------------------------------------------------------------

 15:     The Stokes equation is given by the partial differential equation
 16:   
 17:         -alpha*Laplacian u - \nabla p = f,  0 < x,y < 1,

 19:           \nabla \cdot u              = 0
 20:   
 21:     with boundary conditions
 22:    
 23:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 24:   
 25:     A P2/P1 finite element approximation is used to discretize the boundary
 26:     value problem on the two triangles which make up each rectangle in the DMDA
 27:     to obtain a nonlinear system of equations.

 29:   ------------------------------------------------------------------------- */

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

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

 57: typedef struct {
 58:   PetscScalar u;
 59:   PetscScalar v;
 60:   PetscScalar p;
 61: } Field;

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

 70: static PetscScalar Identity[9] = {0.0833333, 0.0416667, 0.0416667,
 71:                                   0.0416667, 0.0833333, 0.0416667,
 72:                                   0.0416667, 0.0416667, 0.0833333};

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

 81: static PetscScalar Divergence[18] = {-0.1666667, 0.1666667, 0.0,
 82:                                      -0.1666667, 0.0,       0.1666667,
 83:                                      -0.1666667, 0.1666667, 0.0,
 84:                                      -0.1666667, 0.0,       0.1666667,
 85:                                      -0.1666667, 0.1666667, 0.0,
 86:                                      -0.1666667, 0.0,       0.1666667};

 88: /* These are */
 89: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
 90:                                     0.07503111, 0.64494897,
 91:                                     0.66639025, 0.15505103,
 92:                                     0.28001992, 0.64494897};
 93: static PetscScalar quadWeights[4] = {0.15902069,  0.09097931,  0.15902069,  0.09097931};

 95: /* 
 96:    User-defined routines
 97: */

107: int main(int argc,char **argv)
108: {
109:   DMMG                  *dmmg;                 /* hierarchy manager */
110:   DM                     da;
111:   SNES                   snes;                 /* nonlinear solver */
112:   AppCtx                *user;                 /* user-defined work context */
113:   PetscBag               bag;
114:   PetscInt               its;                  /* iterations for convergence */
115:   SNESConvergedReason    reason;
116:   PetscBool              drawContours;         /* flag for drawing contours */
117:   PetscErrorCode         ierr;
118:   PetscReal              lambda_max = 6.81, lambda_min = 0.0, error;

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Initialize program
122:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123:   PetscInitialize(&argc,&argv,(char *)0,help);

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:      Initialize problem parameters
127:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128:   PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
129:   PetscBagGetData(bag, (void **) &user);
130:   PetscBagSetName(bag, "params", "Parameters for SNES example 4");
131:   PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
132:   PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
133:   PetscBagSetFromOptions(bag);
134:   PetscOptionsGetReal(PETSC_NULL,"-alpha",&user->alpha,PETSC_NULL);
135:   PetscOptionsGetReal(PETSC_NULL,"-lambda",&user->lambda,PETSC_NULL);
136:   if (user->lambda > lambda_max || user->lambda < lambda_min) {
137:     SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
138:   }

140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:      Create multilevel DM data structure (DMMG) to manage hierarchical solvers
142:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143:   DMMGCreate(PETSC_COMM_WORLD,1,user,&dmmg);

145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:      Create distributed array (DMDA) to manage parallel grid and vectors
147:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,
149:                     3,1,PETSC_NULL,PETSC_NULL,&da);
150:   DMDASetFieldName(da, 0, "ooblek");
151:   DMMGSetDM(dmmg, (DM) da);
152:   DMDestroy(&da);

154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:      Set the discretization functions
156:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157:   DMMGSetSNESLocal(dmmg, FormFunctionLocal, FormJacobianLocal, 0, 0);
158:   DMMGSetFromOptions(dmmg);

160:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161:      Set the problem null space
162:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163:   DMMGSetNullSpace(dmmg, PETSC_FALSE, 1, CreateNullSpace);

165:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166:      Evaluate initial guess
167:      Note: The user should initialize the vector, x, with the initial guess
168:      for the nonlinear solver prior to calling SNESSolve().  In particular,
169:      to employ an initial guess of zero, the user should explicitly set
170:      this vector to zero by calling VecSet().
171:   */
172:   DMMGSetInitialGuess(dmmg, FormInitialGuess);

174:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175:      Solve nonlinear system
176:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177:   DMMGSolve(dmmg);
178:   snes = DMMGGetSNES(dmmg);
179:   SNESGetIterationNumber(snes,&its);
180:   SNESGetConvergedReason(snes, &reason);

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
184:   L_2Error(DMMGGetDM(dmmg), DMMGGetx(dmmg), &error, user);
185:   PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %G\n", error);
186:   PrintVector(dmmg[0], DMMGGetx(dmmg));

188:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189:      Free work space.  All PETSc objects should be destroyed when they
190:      are no longer needed.
191:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192:   DMMGDestroy(dmmg);
193:   PetscBagDestroy(&bag);
194:   PetscFinalize();
195:   return(0);
196: }

200: PetscErrorCode PrintVector(DMMG dmmg, Vec U)
201: {
202:   DM             da =  dmmg->dm;
203:   Field        **u;
204:   PetscInt       i,j,xs,ys,xm,ym;

208:   DMDAVecGetArray(da,U,&u);
209:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
210:   for(j = ys+ym-1; j >= ys; j--) {
211:     for(i = xs; i < xs+xm; i++) {
212:       printf("u[%d][%d] = (%g, %g, %g) ", j, i, u[j][i].u, u[j][i].v, u[j][i].p);
213:     }
214:     printf("\n");
215:   }
216:   DMDAVecRestoreArray(da,U,&u);
217:   return(0);
218: }

222: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, Field *u)
223: {
225:   (*u).u = x;
226:   (*u).v = -y;
227:   (*u).p = 0.0;
228:   return(0);
229: }

233: PetscErrorCode CreateNullSpace(DMMG dmmg, Vec *nulls)
234: {
235:   DM             da =  dmmg->dm;
236:   Vec            X = nulls[0];
237:   Field        **x;
238:   PetscInt       xs,ys,xm,ym,i,j;

242:   DMDAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL);
243:   DMDAVecGetArray(da, X, &x);
244:   for(j = ys; j < ys+ym; j++) {
245:     for(i = xs; i < xs+xm; i++) {
246:       x[j][i].u = 0.0;
247:       x[j][i].v = 0.0;
248:       x[j][i].p = 1.0;
249:     }
250:   }
251:   DMDAVecRestoreArray(da, X, &x);
252:   PrintVector(dmmg, X);
253:   return(0);
254: }

258: /* 
259:    FormInitialGuess - Forms initial approximation.

261:    Input Parameters:
262:    dmmg - The DMMG context
263:    X - vector

265:    Output Parameter:
266:    X - vector
267: */
268: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
269: {
270:   AppCtx        *user = (AppCtx *) dmmg->user;
271:   DM             da =  dmmg->dm;
272:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
274:   PetscReal      lambda,temp1,temp,hx,hy;
275:   Field        **x;

278:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
279:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

281:   lambda = user->lambda;
282:   hx     = 1.0/(PetscReal)(Mx-1);
283:   hy     = 1.0/(PetscReal)(My-1);
284:   if (lambda == 0.0) {
285:     temp1  = 0.0;
286:   } else {
287:     temp1  = lambda/(lambda + 1.0);
288:   }

290:   /*
291:      Get a pointer to vector data.
292:        - For default PETSc vectors, VecGetArray() returns a pointer to
293:          the data array.  Otherwise, the routine is implementation dependent.
294:        - You MUST call VecRestoreArray() when you no longer need access to
295:          the array.
296:   */
297:   DMDAVecGetArray(da,X,&x);

299:   /*
300:      Get local grid boundaries (for 2-dimensional DMDA):
301:        xs, ys   - starting grid indices (no ghost points)
302:        xm, ym   - widths of local grid (no ghost points)

304:   */
305:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

307:   /*
308:      Compute initial guess over the locally owned part of the grid
309:   */
310:   for (j=ys; j<ys+ym; j++) {
311:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
312:     for (i=xs; i<xs+xm; i++) {
313: #define CHECK_SOLUTION
314: #ifdef CHECK_SOLUTION
315:         ExactSolution(i*hx, j*hy, &x[j][i]);
316: #else
317:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
318:         /* Boundary conditions are usually zero Dirichlet */
319:         ExactSolution(i*hx, j*hy, &x[j][i]);
320:       } else {
321:         x[j][i].u = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
322:         x[j][i].v = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
323:         x[j][i].p = 1.0;
324:       }
325: #endif
326:     }
327:   }

329:   /*
330:      Restore vector
331:   */
332:   DMDAVecRestoreArray(da,X,&x);
333:   PrintVector(dmmg, X);
334:   return(0);
335: }

339: PetscErrorCode constantResidual(PetscReal lambda, PetscBool  isLower, int i, int j, PetscReal hx, PetscReal hy, Field r[])
340: {
341:   Field       rLocal[3] = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
342:   PetscScalar phi[3] = {0.0, 0.0, 0.0};
343:   PetscReal   xI = i*hx, yI = j*hy, hxhy = hx*hy, x, y;
344:   Field       res;
345:   PetscInt    q, k;

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:     if (isLower) {
353:       x    = xI + quadPoints[q*2]*hx;
354:       y    = yI + quadPoints[q*2+1]*hy;
355:     } else {
356:       x    = xI + 1.0 - quadPoints[q*2]*hx;
357:       y    = yI + 1.0 - quadPoints[q*2+1]*hy;
358:     }
359:     res.u  = quadWeights[q]*(0.0);
360:     res.v  = quadWeights[q]*(0.0);
361:     res.p  = quadWeights[q]*(0.0);
362:     for(k = 0; k < 3; k++) {
363:       rLocal[k].u += phi[k]*res.u;
364:       rLocal[k].v += phi[k]*res.v;
365:       rLocal[k].p += phi[k]*res.p;
366:     }
367:   }
368:   for(k = 0; k < 3; k++) {
369:     printf("  constLocal[%d] = (%g, %g, %g)\n", k, lambda*hxhy*rLocal[k].u, lambda*hxhy*rLocal[k].v, hxhy*rLocal[k].p);
370:     r[k].u += lambda*hxhy*rLocal[k].u;
371:     r[k].v += lambda*hxhy*rLocal[k].v;
372:     r[k].p += hxhy*rLocal[k].p;
373:   }
374:   return(0);
375: }

377: PetscErrorCode nonlinearResidual(PetscReal lambda, Field u[], Field r[]) {
378:   Field       rLocal[3] = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
379:   PetscScalar phi[3] = {0.0, 0.0, 0.0};
380:   Field       res;
381:   PetscInt q;

384:   for(q = 0; q < 4; q++) {
385:     phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
386:     phi[1] = quadPoints[q*2];
387:     phi[2] = quadPoints[q*2+1];
388:     res.u  = quadWeights[q]*PetscExpScalar(u[0].u*phi[0] + u[1].u*phi[1] + u[2].u*phi[2]);
389:     res.v  = quadWeights[q]*PetscExpScalar(u[0].v*phi[0] + u[1].v*phi[1] + u[2].v*phi[2]);
390:     rLocal[0].u += phi[0]*res.u;
391:     rLocal[0].v += phi[0]*res.v;
392:     rLocal[1].u += phi[1]*res.u;
393:     rLocal[1].v += phi[1]*res.v;
394:     rLocal[2].u += phi[2]*res.u;
395:     rLocal[2].v += phi[2]*res.v;
396:   }
397:   r[0].u += lambda*rLocal[0].u;
398:   r[0].v += lambda*rLocal[0].v;
399:   r[1].u += lambda*rLocal[1].u;
400:   r[1].v += lambda*rLocal[1].v;
401:   r[2].u += lambda*rLocal[2].u;
402:   r[2].v += lambda*rLocal[2].v;
403:   return(0);
404: }

408: /* 
409:    FormFunctionLocal - Evaluates nonlinear function, F(x).

411:        Process adiC(36): FormFunctionLocal

413:  */
414: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, AppCtx *user)
415: {
416:   Field          uLocal[3];
417:   Field          rLocal[3];
418:   PetscScalar    G[4];
419:   Field          uExact;
420:   PetscReal      alpha,lambda,hx,hy,hxhy,sc,detJInv;
421:   PetscInt       i,j,k,l;


426:   /* Naive Jacobian calculation:

428:      J = / 1/hx  0   \ J^{-1} = / hx   0 \  1/|J| = hx*hy = |J^{-1}|
429:          \  0   1/hy /          \  0  hy /
430:    */
431:   alpha   = user->alpha;
432:   lambda  = user->lambda;
433:   hx      = 1.0/(PetscReal)(info->mx-1);
434:   hy      = 1.0/(PetscReal)(info->my-1);
435:   sc      = hx*hy*lambda;
436:   hxhy    = hx*hy;
437:   detJInv = hxhy;
438:   G[0] = (1.0/(hx*hx)) * detJInv;
439:   G[1] = 0.0;
440:   G[2] = G[1];
441:   G[3] = (1.0/(hy*hy)) * detJInv;
442:   for(k = 0; k < 4; k++) {
443:     printf("G[%d] = %g\n", k, G[k]);
444:   }

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

451:        2 (1)    (0)
452:      i,j+1 --- i+1,j+1
453:        |  \      |
454:        |   \     |
455:        |    \    |
456:        |     \   |
457:        |      \  |
458:       i,j  --- i+1,j
459:        0         1 (2)

461:      and therefore we do not loop over the last vertex in each dimension.
462:   */
463:   for(j = info->ys; j < info->ys+info->ym-1; j++) {
464:     for(i = info->xs; i < info->xs+info->xm-1; i++) {
465:       /* Lower element */
466:       uLocal[0] = x[j][i];
467:       uLocal[1] = x[j][i+1];
468:       uLocal[2] = x[j+1][i];
469:       printf("Lower Solution ElementVector for (%d, %d)\n", i, j);
470:       for(k = 0; k < 3; k++) {
471:         printf("  uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);
472:       }
473:       for(k = 0; k < 3; k++) {
474:         rLocal[k].u = 0.0;
475:         rLocal[k].v = 0.0;
476:         rLocal[k].p = 0.0;
477:         for(l = 0; l < 3; l++) {
478:           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;
479:           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;
480:           /* rLocal[k].u += hxhy*Identity[k*3+l]*uLocal[l].u; */
481:           /* rLocal[k].p += hxhy*Identity[k*3+l]*uLocal[l].p; */
482:           /* Gradient */
483:           rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
484:           rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
485:           /* Divergence */
486:           rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
487:         }
488:       }
489:       printf("Lower Laplacian ElementVector for (%d, %d)\n", i, j);
490:       for(k = 0; k < 3; k++) {
491:         printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
492:       }
493:       constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
494:       printf("Lower Laplacian+Constant 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:       nonlinearResidual(0.0*sc, uLocal, rLocal);
499:       printf("Lower Full nonlinear 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:       f[j][i].u   += rLocal[0].u;
504:       f[j][i].v   += rLocal[0].v;
505:       f[j][i].p   += rLocal[0].p;
506:       f[j][i+1].u += rLocal[1].u;
507:       f[j][i+1].v += rLocal[1].v;
508:       f[j][i+1].p += rLocal[1].p;
509:       f[j+1][i].u += rLocal[2].u;
510:       f[j+1][i].v += rLocal[2].v;
511:       f[j+1][i].p += rLocal[2].p;
512:       /* Upper element */
513:       uLocal[0] = x[j+1][i+1];
514:       uLocal[1] = x[j+1][i];
515:       uLocal[2] = x[j][i+1];
516:       printf("Upper Solution ElementVector for (%d, %d)\n", i, j);
517:       for(k = 0; k < 3; k++) {
518:         printf("  uLocal[%d] = (%g, %g, %g)\n", k, uLocal[k].u, uLocal[k].v, uLocal[k].p);
519:       }
520:       for(k = 0; k < 3; k++) {
521:         rLocal[k].u = 0.0;
522:         rLocal[k].v = 0.0;
523:         rLocal[k].p = 0.0;
524:         for(l = 0; l < 3; l++) {
525:           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;
526:           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;
527:           /* rLocal[k].p += Identity[k*3+l]*uLocal[l].p; */
528:           /* Gradient */
529:           rLocal[k].u += hx*Gradient[(k*2+0)*3 + l]*uLocal[l].p;
530:           rLocal[k].v += hy*Gradient[(k*2+1)*3 + l]*uLocal[l].p;
531:           /* Divergence */
532:           rLocal[k].p += hx*Divergence[(k*2+0)*3 + l]*uLocal[l].u + hy*Divergence[(k*2+1)*3 + l]*uLocal[l].v;
533:         }
534:       }
535:       printf("Upper Laplacian ElementVector for (%d, %d)\n", i, j);
536:       for(k = 0; k < 3; k++) {
537:         printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
538:       }
539:       constantResidual(1.0, PETSC_BOOL, i, j, hx, hy, rLocal);
540:       printf("Upper Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
541:       for(k = 0; k < 3; k++) {
542:         printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
543:       }
544:       nonlinearResidual(0.0*sc, uLocal, rLocal);
545:       printf("Upper Full nonlinear ElementVector for (%d, %d)\n", i, j);
546:       for(k = 0; k < 3; k++) {
547:         printf("  rLocal[%d] = (%g, %g, %g)\n", k, rLocal[k].u, rLocal[k].v, rLocal[k].p);
548:       }
549:       f[j+1][i+1].u += rLocal[0].u;
550:       f[j+1][i+1].v += rLocal[0].v;
551:       f[j+1][i+1].p += rLocal[0].p;
552:       f[j+1][i].u   += rLocal[1].u;
553:       f[j+1][i].v   += rLocal[1].v;
554:       f[j+1][i].p   += rLocal[1].p;
555:       f[j][i+1].u   += rLocal[2].u;
556:       f[j][i+1].v   += rLocal[2].v;
557:       f[j][i+1].p   += rLocal[2].p;
558:       /* Boundary conditions */
559:       if (i == 0 || j == 0) {
560:         ExactSolution(i*hx, j*hy, &uExact);
561:         f[j][i].u = x[j][i].u - uExact.u;
562:         f[j][i].v = x[j][i].v - uExact.v;
563:         f[j][i].p = x[j][i].p - uExact.p;
564:       }
565:       if ((i == info->mx-2) || (j == 0)) {
566:         ExactSolution((i+1)*hx, j*hy, &uExact);
567:         f[j][i+1].u = x[j][i+1].u - uExact.u;
568:         f[j][i+1].v = x[j][i+1].v - uExact.v;
569:         f[j][i+1].p = x[j][i+1].p - uExact.p;
570:       }
571:       if ((i == info->mx-2) || (j == info->my-2)) {
572:         ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
573:         f[j+1][i+1].u = x[j+1][i+1].u - uExact.u;
574:         f[j+1][i+1].v = x[j+1][i+1].v - uExact.v;
575:         f[j+1][i+1].p = x[j+1][i+1].p - uExact.p;
576:       }
577:       if ((i == 0) || (j == info->my-2)) {
578:         ExactSolution(i*hx, (j+1)*hy, &uExact);
579:         f[j+1][i].u = x[j+1][i].u - uExact.u;
580:         f[j+1][i].v = x[j+1][i].v - uExact.v;
581:         f[j+1][i].p = x[j+1][i].p - uExact.p;
582:       }
583:     }
584:   }

586:   for(j = info->ys+info->ym-1; j >= info->ys; j--) {
587:     for(i = info->xs; i < info->xs+info->xm; i++) {
588:       printf("f[%d][%d] = (%g, %g, %g) ", j, i, f[j][i].u, f[j][i].v, f[j][i].p);
589:     }
590:     printf("\n");
591:   }
592:   PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
593:   return(0);
594: }

596: PetscErrorCode nonlinearJacobian(PetscReal lambda, Field u[], PetscScalar J[]) {
598:   return(0);
599: }

603: /*
604:    FormJacobianLocal - Evaluates Jacobian matrix.
605: */
606: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field **x, Mat jac, AppCtx *user)
607: {
608:   Field          uLocal[4];
609:   PetscScalar    JLocal[144];
610:   MatStencil     rows[4*3], cols[4*3], ident;
611:   PetscInt       lowerRow[3] = {0, 1, 3};
612:   PetscInt       upperRow[3] = {2, 3, 1};
613:   PetscInt       hasLower[3], hasUpper[3], velocityRows[4], pressureRows[4];
614:   PetscScalar    alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
615:   PetscInt       i,j,k,l,numRows,dof = info->dof;

619:   alpha  = user->alpha;
620:   lambda = user->lambda;
621:   hx     = 1.0/(PetscReal)(info->mx-1);
622:   hy     = 1.0/(PetscReal)(info->my-1);
623:   sc     = hx*hy*lambda;
624:   hxhy   = hx*hy;
625:   detJInv = hxhy;
626:   G[0] = (1.0/(hx*hx)) * detJInv;
627:   G[1] = 0.0;
628:   G[2] = G[1];
629:   G[3] = (1.0/(hy*hy)) * detJInv;
630:   for(k = 0; k < 4; k++) {
631:     printf("G[%d] = %g\n", k, G[k]);
632:   }

634:   MatZeroEntries(jac);
635:   /* 
636:      Compute entries for the locally owned part of the Jacobian.
637:       - Currently, all PETSc parallel matrix formats are partitioned by
638:         contiguous chunks of rows across the processors. 
639:       - Each processor needs to insert only elements that it owns
640:         locally (but any non-local elements will be sent to the
641:         appropriate processor during matrix assembly). 
642:       - Here, we set all entries for a particular row at once.
643:       - We can set matrix entries either using either
644:         MatSetValuesLocal() or MatSetValues(), as discussed above.
645:   */
646: #define NOT_PRES_BC 1
647:   for (j=info->ys; j<info->ys+info->ym-1; j++) {
648:     for (i=info->xs; i<info->xs+info->xm-1; i++) {
649:       PetscMemzero(JLocal, 144 * sizeof(PetscScalar));
650:       numRows = 0;
651:       /* Lower element */
652:       uLocal[0] = x[j][i];
653:       uLocal[1] = x[j][i+1];
654:       uLocal[2] = x[j+1][i+1];
655:       uLocal[3] = x[j+1][i];
656:       /* i,j */
657:       if (i == 0 || j == 0) {
658:         hasLower[0] = 0;
659:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
660:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
661:         ident.i = i; ident.j = j; ident.c = 0;
662:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
663:         ident.i = i; ident.j = j; ident.c = 1;
664:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
665: #ifdef PRES_BC
666:         ident.i = i; ident.j = j; ident.c = 2;
667:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
668: #endif
669:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
670:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
671:       } else {
672:         hasLower[0] = 1;
673:         velocityRows[0] = numRows;
674:         rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 0;
675:         numRows++;
676:         rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 1;
677:         numRows++;
678: #ifdef PRES_BC
679:         pressureRows[0] = numRows;
680:         rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
681:         numRows++;
682: #endif
683:       }
684: #ifndef PRES_BC
685:       pressureRows[0] = numRows;
686:       rows[numRows].i = i; rows[numRows].j = j; rows[numRows].c = 2;
687:       numRows++;
688: #endif
689:       cols[0*dof+0].i = i; cols[0*dof+0].j = j; cols[0*dof+0].c = 0;
690:       cols[0*dof+1].i = i; cols[0*dof+1].j = j; cols[0*dof+1].c = 1;
691:       cols[0*dof+2].i = i; cols[0*dof+2].j = j; cols[0*dof+2].c = 2;
692:       /* i+1,j */
693:       if ((i == info->mx-2) || (j == 0)) {
694:         hasLower[1] = 0;
695:         hasUpper[2] = 0;
696:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
697:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
698:         ident.i = i+1; ident.j = j; ident.c = 0;
699:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
700:         ident.i = i+1; ident.j = j; ident.c = 1;
701:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
702: #ifdef PRES_BC
703:         ident.i = i+1; ident.j = j; ident.c = 2;
704:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
705: #endif
706:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
707:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
708:       } else {
709:         hasLower[1] = 1;
710:         hasUpper[2] = 1;
711:         velocityRows[1] = numRows;
712:         rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 0;
713:         numRows++;
714:         rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 1;
715:         numRows++;
716: #ifdef PRES_BC
717:         pressureRows[1] = numRows;
718:         rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
719:         numRows++;
720: #endif
721:       }
722: #ifndef PRES_BC
723:       pressureRows[1] = numRows;
724:       rows[numRows].i = i+1; rows[numRows].j = j; rows[numRows].c = 2;
725:       numRows++;
726: #endif
727:       cols[1*dof+0].i = i+1; cols[1*dof+0].j = j; cols[1*dof+0].c = 0;
728:       cols[1*dof+1].i = i+1; cols[1*dof+1].j = j; cols[1*dof+1].c = 1;
729:       cols[1*dof+2].i = i+1; cols[1*dof+2].j = j; cols[1*dof+2].c = 2;
730:       /* i+1,j+1 */
731:       if ((i == info->mx-2) || (j == info->my-2)) {
732:         hasUpper[0] = 0;
733:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
734:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
735:         ident.i = i+1; ident.j = j+1; ident.c = 0;
736:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
737:         ident.i = i+1; ident.j = j+1; ident.c = 1;
738:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
739: #ifdef PRES_BC
740:         ident.i = i+1; ident.j = j+1; ident.c = 2;
741:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
742: #endif
743:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
744:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
745:       } else {
746:         hasUpper[0] = 1;
747:         velocityRows[2] = numRows;
748:         rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 0;
749:         numRows++;
750:         rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 1;
751:         numRows++;
752: #ifdef PRES_BC
753:         pressureRows[2] = numRows;
754:         rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
755:         numRows++;
756: #endif
757:       }
758: #ifndef PRES_BC
759:       pressureRows[2] = numRows;
760:       rows[numRows].i = i+1; rows[numRows].j = j+1; rows[numRows].c = 2;
761:       numRows++;
762: #endif
763:       cols[2*dof+0].i = i+1; cols[2*dof+0].j = j+1; cols[2*dof+0].c = 0;
764:       cols[2*dof+1].i = i+1; cols[2*dof+1].j = j+1; cols[2*dof+1].c = 1;
765:       cols[2*dof+2].i = i+1; cols[2*dof+2].j = j+1; cols[2*dof+2].c = 2;
766:       /* i,j+1 */
767:       if ((i == 0) || (j == info->my-2)) {
768:         hasLower[2] = 0;
769:         hasUpper[1] = 0;
770:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
771:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
772:         ident.i = i; ident.j = j+1; ident.c = 0;
773:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
774:         ident.i = i; ident.j = j+1; ident.c = 1;
775:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
776: #ifdef PRES_BC
777:         ident.i = i; ident.j = j+1; ident.c = 2;
778:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
779: #endif
780:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
781:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
782:       } else {
783:         hasLower[2] = 1;
784:         hasUpper[1] = 1;
785:         velocityRows[3] = numRows;
786:         rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 0;
787:         numRows++;
788:         rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 1;
789:         numRows++;
790: #ifdef PRES_BC
791:         pressureRows[3] = numRows;
792:         rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
793:         numRows++;
794: #endif
795:       }
796: #ifndef PRES_BC
797:       pressureRows[3] = numRows;
798:       rows[numRows].i = i; rows[numRows].j = j+1; rows[numRows].c = 2;
799:       numRows++;
800: #endif
801:       cols[3*dof+0].i = i; cols[3*dof+0].j = j+1; cols[3*dof+0].c = 0;
802:       cols[3*dof+1].i = i; cols[3*dof+1].j = j+1; cols[3*dof+1].c = 1;
803:       cols[3*dof+2].i = i; cols[3*dof+2].j = j+1; cols[3*dof+2].c = 2;

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

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

868:   /* 
869:      Assemble matrix, using the 2-step process:
870:        MatAssemblyBegin(), MatAssemblyEnd().
871:   */
872:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
873:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
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);
879:   return(0);
880: }

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

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

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

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

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