Actual source code: eptorsion1.c

tao-2.1-p0 2012-07-24
  1: /* Program usage: mpirun -np 1 eptorsion1 [-help] [all TAO options] */

  3: /* ----------------------------------------------------------------------

  5:   Elastic-plastic torsion problem.  

  7:   The elastic plastic torsion problem arises from the determination 
  8:   of the stress field on an infinitely long cylindrical bar, which is
  9:   equivalent to the solution of the following problem:

 11:   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
 12:   
 13:   where C is the torsion angle per unit length.

 15:   The multiprocessor version of this code is eptorsion2.c.

 17: ---------------------------------------------------------------------- */

 19: /*
 20:   Include "taosolver.h" so that we can use TAO solvers.  Note that this 
 21:   file automatically includes files for lower-level support, such as those
 22:   provided by the PETSc library:
 23:      petsc.h       - base PETSc routines   petscvec.h - vectors
 24:      petscsys.h    - sysem routines        petscmat.h - matrices
 25:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 26:      petscviewer.h - viewers               petscpc.h  - preconditioners
 27: */

 29: #include "taosolver.h"


 32: static  char help[]=
 33: "Demonstrates use of the TAO package to solve \n\
 34: unconstrained minimization problems on a single processor.  This example \n\
 35: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
 36: test suite.\n\
 37: The command line options are:\n\
 38:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 39:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 40:   -par <param>, where <param> = angle of twist per unit length\n\n";

 42: /*T 
 43:    Concepts: TAO - Solving an unconstrained minimization problem
 44:    Routines: TaoInitialize(); TaoFinalize(); 
 45:    Routines: TaoCreate(); TaoSetType();
 46:    Routines: TaoSetInitialVector(); 
 47:    Routines: TaoSetObjectiveAndGradientRoutine();
 48:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 49:    Routines: TaoGetKSP(); TaoSolve();
 50:    Routines: TaoGetTerminationReason(); TaoDestroy();
 51:    Processors: 1
 52: T*/ 

 54: /* 
 55:    User-defined application context - contains data needed by the 
 56:    application-provided call-back routines, FormFunction(),
 57:    FormGradient(), and FormHessian().
 58: */

 60: typedef struct {
 61:    PetscReal  param;      /* nonlinearity parameter */
 62:    PetscInt   mx, my;     /* discretization in x- and y-directions */
 63:    PetscInt   ndim;       /* problem dimension */
 64:    Vec     s, y, xvec; /* work space for computing Hessian */
 65:    PetscReal  hx, hy;     /* mesh spacing in x- and y-directions */
 66: } AppCtx;

 68: /* -------- User-defined Routines --------- */

 70: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
 71: PetscErrorCode FormFunction(TaoSolver,Vec,PetscReal*,void*);
 72: PetscErrorCode FormGradient(TaoSolver,Vec,Vec,void*);
 73: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*, MatStructure *,void*);
 74: PetscErrorCode HessianProductMat(Mat,Vec,Vec);
 75: PetscErrorCode HessianProduct(void*,Vec,Vec);
 76: PetscErrorCode MatrixFreeHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure*,void*);
 77: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal *,Vec,void *);

 81: PetscErrorCode main(int argc,char **argv)
 82: {
 84:   PetscInt      mx=10;               /* discretization in x-direction */
 85:   PetscInt      my=10;               /* discretization in y-direction */
 86:   Vec         x;                   /* solution, gradient vectors */
 87:   PetscBool   flg;                 /* A return value when checking for use options */
 88:   TaoSolver   tao;                 /* TaoSolver solver context */
 89:   Mat         H;                   /* Hessian matrix */
 90:   TaoSolverTerminationReason reason;        
 91:   KSP         ksp;                 /* PETSc Krylov subspace solver */
 92:   AppCtx      user;                /* application context */
 93:   PetscMPIInt size;                /* number of processes */
 94:   PetscReal one=1.0;


 97:   /* Initialize TAO,PETSc */
 98:   PetscInitialize(&argc,&argv,(char *)0,help);
 99:   TaoInitialize(&argc,&argv,(char *)0,help);

101:   /* Optional:  Check  that only one processor is being used. */
102:   MPI_Comm_size(MPI_COMM_WORLD,&size); 
103:   if (size >1) {
104:     PetscPrintf(PETSC_COMM_SELF,"This example is intended for single processor use!\n");
105:     PetscPrintf(PETSC_COMM_SELF,"Try the example eptorsion2!\n");
106:     SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");
107:   }

109:   /* Specify default parameters for the problem, check for command-line overrides */
110:   user.param = 5.0;
111:   PetscOptionsGetInt(PETSC_NULL,"-my",&my,&flg); 
112:   PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,&flg); 
113:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,&flg); 


116:   PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
117:   PetscPrintf(PETSC_COMM_SELF,"mx: %D     my: %D   \n\n",mx,my);  
118:   user.ndim = mx * my; user.mx = mx; user.my = my;

120:   user.hx = one/(mx+1); user.hy = one/(my+1);


123:   /* Allocate vectors */
124:   VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y); 
125:   VecDuplicate(user.y,&user.s); 
126:   VecDuplicate(user.y,&x); 

128:   /* The TAO code begins here */

130:   /* Create TAO solver and set desired solution method */
131:   TaoCreate(PETSC_COMM_SELF,&tao); 
132:   TaoSetType(tao,"tao_lmvm"); 

134:   /* Set solution vector with an initial guess */
135:   FormInitialGuess(&user,x); 
136:   TaoSetInitialVector(tao,x); 

138:   /* Set routine for function and gradient evaluation */
139:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user); 

141:   /* From command line options, determine if using matrix-free hessian */
142:   PetscOptionsHasName(PETSC_NULL,"-my_tao_mf",&flg); 
143:   if (flg) {
144:     MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,
145:                           user.ndim,(void*)&user,&H); 
146:     MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat); CHKERRQ
147: (ierr);
148:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE); 

150:     TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user); 


153:     /* Set null preconditioner.  Alternatively, set user-provided 
154:        preconditioner or explicitly form preconditioning matrix */
155:     PetscOptionsSetValue("-tao_pc_type","none"); 

157:   } else {

159:     MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,PETSC_NULL,&H); 
160:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE); 

162:     TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user); 

164:   }



168:   /* Modify the PETSc KSP structure */
169:   PetscOptionsSetValue("-tao_ksp_type","cg"); 

171:   /* Check for any TAO command line options */
172:   TaoSetFromOptions(tao); 


175:   /* SOLVE THE APPLICATION */
176:   TaoSolve(tao);  
177:   TaoGetKSP(tao,&ksp); 
178:   if (ksp) {
179:     KSPView(ksp,PETSC_VIEWER_STDOUT_SELF); 
180:   }

182:   /* 
183:      To View TAO solver information use
184:       TaoView(tao); 
185:   */

187:   /* Get information on termination */
188:   TaoGetTerminationReason(tao,&reason); 
189:   if (reason <= 0){
190:     PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
191:   }

193:   /* Free TAO data structures */
194:   TaoDestroy(&tao); 

196:   /* Free PETSc data structures */
197:   VecDestroy(&user.s); 
198:   VecDestroy(&user.y); 
199:   VecDestroy(&x); 
200:   MatDestroy(&H); 

202:   /* Finalize TAO, PETSc */
203:   TaoFinalize();
204:   PetscFinalize();

206:   return 0;
207: }

209: /* ------------------------------------------------------------------- */
212: /* 
213:     FormInitialGuess - Computes an initial approximation to the solution.

215:     Input Parameters:
216: .   user - user-defined application context
217: .   X    - vector

219:     Output Parameters:
220: .   X    - vector
221: */
222: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
223: {
224:   PetscReal hx = user->hx, hy = user->hy, temp;
225:   PetscReal val;
227:   PetscInt i, j, k, nx = user->mx, ny = user->my;

229:   /* Compute initial guess */
230:   for (j=0; j<ny; j++) {
231:     temp = PetscMin(j+1,ny-j)*hy;
232:     for (i=0; i<nx; i++) {
233:       k   = nx*j + i;
234:       val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
235:       VecSetValues(X,1,&k,&val,ADD_VALUES); 
236:     }
237:   }
238:   VecAssemblyBegin(X); 
239:   VecAssemblyEnd(X); 
240:   return 0;
241: }
242: /* ------------------------------------------------------------------- */
245: /* 
246:    FormFunctionGradient - Evaluates the function and corresponding gradient.
247:     
248:    Input Parameters:
249:    tao - the TaoSolver context
250:    X   - the input vector 
251:    ptr - optional user-defined context, as set by TaoSetFunction()

253:    Output Parameters:
254:    f   - the newly evaluated function
255:    G   - the newly evaluated gradient
256: */
257: PetscErrorCode FormFunctionGradient(TaoSolver tao,Vec X,PetscReal *f,Vec G,void *ptr)
258: {
260:   FormFunction(tao,X,f,ptr);
261:   FormGradient(tao,X,G,ptr);
262:   return 0;
263: }
264: /* ------------------------------------------------------------------- */
267: /* 
268:    FormFunction - Evaluates the function, f(X).

270:    Input Parameters:
271: .  tao - the TaoSolver context
272: .  X   - the input vector 
273: .  ptr - optional user-defined context, as set by TaoSetFunction()

275:    Output Parameters:
276: .  f    - the newly evaluated function
277: */
278: PetscErrorCode FormFunction(TaoSolver tao,Vec X,PetscReal *f,void *ptr)
279: {
280:   AppCtx *user = (AppCtx *) ptr;
281:   PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
282:   PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
283:   PetscReal v, cdiv3 = user->param/three;
284:   PetscReal *x;
286:   PetscInt  nx = user->mx, ny = user->my, i, j, k;

288:   /* Get pointer to vector data */
289:   VecGetArray(X,&x); 

291:   /* Compute function contributions over the lower triangular elements */
292:   for (j=-1; j<ny; j++) {
293:     for (i=-1; i<nx; i++) {
294:       k = nx*j + i;
295:       v = zero;
296:       vr = zero;
297:       vt = zero;
298:       if (i >= 0 && j >= 0) v = x[k];
299:       if (i < nx-1 && j > -1) vr = x[k+1];
300:       if (i > -1 && j < ny-1) vt = x[k+nx];
301:       dvdx = (vr-v)/hx;
302:       dvdy = (vt-v)/hy;
303:       fquad += dvdx*dvdx + dvdy*dvdy;
304:       flin -= cdiv3*(v+vr+vt);
305:     }
306:   }

308:   /* Compute function contributions over the upper triangular elements */
309:   for (j=0; j<=ny; j++) {
310:     for (i=0; i<=nx; i++) {
311:       k = nx*j + i;
312:       vb = zero;
313:       vl = zero;
314:       v = zero;
315:       if (i < nx && j > 0) vb = x[k-nx];
316:       if (i > 0 && j < ny) vl = x[k-1];
317:       if (i < nx && j < ny) v = x[k];
318:       dvdx = (v-vl)/hx;
319:       dvdy = (v-vb)/hy;
320:       fquad = fquad + dvdx*dvdx + dvdy*dvdy;
321:       flin = flin - cdiv3*(vb+vl+v);
322:     }
323:   }

325:   /* Restore vector */
326:   VecRestoreArray(X,&x); 

328:   /* Scale the function */
329:   area = p5*hx*hy;
330:   *f = area*(p5*fquad+flin);

332:   PetscLogFlops(nx*ny*24); 
333:   return 0;
334: }
335: /* ------------------------------------------------------------------- */
338: /*  
339:     FormGradient - Evaluates the gradient, G(X).              

341:     Input Parameters:
342: .   tao  - the TaoSolver context
343: .   X    - input vector
344: .   ptr  - optional user-defined context
345:     
346:     Output Parameters:
347: .   G - vector containing the newly evaluated gradient
348: */
349: PetscErrorCode FormGradient(TaoSolver tao,Vec X,Vec G,void *ptr)
350: {
351:   AppCtx *user = (AppCtx *) ptr;
352:   PetscReal zero=0.0, p5=0.5, three = 3.0, area, val, *x;
354:   PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
355:   PetscReal hx = user->hx, hy = user->hy;
356:   PetscReal vb, vl, vr, vt, dvdx, dvdy;
357:   PetscReal v, cdiv3 = user->param/three;

359:   /* Initialize gradient to zero */
360:   VecSet(G, zero); 

362:   /* Get pointer to vector data */
363:   VecGetArray(X,&x); 

365:   /* Compute gradient contributions over the lower triangular elements */
366:   for (j=-1; j<ny; j++) {
367:     for (i=-1; i<nx; i++) {
368:       k  = nx*j + i;
369:       v  = zero;
370:       vr = zero;
371:       vt = zero;
372:       if (i >= 0 && j >= 0)    v = x[k];
373:       if (i < nx-1 && j > -1) vr = x[k+1];
374:       if (i > -1 && j < ny-1) vt = x[k+nx];
375:       dvdx = (vr-v)/hx;
376:       dvdy = (vt-v)/hy;
377:       if (i != -1 && j != -1) {
378:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
379:         VecSetValues(G,1,&ind,&val,ADD_VALUES); 
380:       }
381:       if (i != nx-1 && j != -1) {
382:         ind = k+1; val =  dvdx/hx - cdiv3;
383:         VecSetValues(G,1,&ind,&val,ADD_VALUES); 
384:       }
385:       if (i != -1 && j != ny-1) {
386:         ind = k+nx; val = dvdy/hy - cdiv3;
387:         VecSetValues(G,1,&ind,&val,ADD_VALUES); 
388:       }
389:     }
390:   }

392:   /* Compute gradient contributions over the upper triangular elements */
393:   for (j=0; j<=ny; j++) {
394:     for (i=0; i<=nx; i++) {
395:       k = nx*j + i;
396:       vb = zero;
397:       vl = zero;
398:       v  = zero;
399:       if (i < nx && j > 0) vb = x[k-nx];
400:       if (i > 0 && j < ny) vl = x[k-1];
401:       if (i < nx && j < ny) v = x[k];
402:       dvdx = (v-vl)/hx;
403:       dvdy = (v-vb)/hy;
404:       if (i != nx && j != 0) {
405:         ind = k-nx; val = - dvdy/hy - cdiv3;
406:         VecSetValues(G,1,&ind,&val,ADD_VALUES); 
407:       }
408:       if (i != 0 && j != ny) {
409:         ind = k-1; val =  - dvdx/hx - cdiv3;
410:         VecSetValues(G,1,&ind,&val,ADD_VALUES); 
411:       }
412:       if (i != nx && j != ny) {
413:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
414:         VecSetValues(G,1,&ind,&val,ADD_VALUES); 
415:       }
416:     }
417:   }

419:   /* Restore vector */
420:   VecRestoreArray(X,&x); 

422:   /* Assemble gradient vector */
423:   VecAssemblyBegin(G); 
424:   VecAssemblyEnd(G); 

426:   /* Scale the gradient */
427:   area = p5*hx*hy;
428:   VecScale(G, area); 
429:   
430:   PetscLogFlops(nx*ny*24); 
431:   return 0;
432: }

434: /* ------------------------------------------------------------------- */
437: /* 
438:    FormHessian - Forms the Hessian matrix.

440:    Input Parameters:
441: .  tao - the TaoSolver context
442: .  X    - the input vector
443: .  ptr  - optional user-defined context, as set by TaoSetHessian()
444:    
445:    Output Parameters:
446: .  H     - Hessian matrix
447: .  PrecH - optionally different preconditioning Hessian
448: .  flag  - flag indicating matrix structure

450:    Notes:
451:    This routine is intended simply as an example of the interface
452:    to a Hessian evaluation routine.  Since this example compute the
453:    Hessian a column at a time, it is not particularly efficient and
454:    is not recommended.
455: */
456: PetscErrorCode FormHessian(TaoSolver tao,Vec X,Mat *HH,Mat *Hpre, MatStructure *flg, void *ptr)
457: {
458:   AppCtx     *user = (AppCtx *) ptr;
460:   PetscInt   i,j, ndim = user->ndim;
461:   PetscReal  *y, zero = 0.0, one = 1.0;
462:   Mat H=*HH;
463:   PetscBool assembled;
464:   *Hpre = H;

466:   /* Set location of vector */
467:   user->xvec = X;

469:   /* Initialize Hessian entries and work vector to zero */
470:   MatAssembled(H,&assembled); 
471:   if (assembled){MatZeroEntries(H);  }

473:   VecSet(user->s, zero); 

475:   /* Loop over matrix columns to compute entries of the Hessian */
476:   for (j=0; j<ndim; j++) {

478:     VecSetValues(user->s,1,&j,&one,INSERT_VALUES); 
479:     VecAssemblyBegin(user->s); 
480:     VecAssemblyEnd(user->s); 

482:     HessianProduct(ptr,user->s,user->y); 

484:     VecSetValues(user->s,1,&j,&zero,INSERT_VALUES); 
485:     VecAssemblyBegin(user->s); 
486:     VecAssemblyEnd(user->s); 

488:     VecGetArray(user->y,&y); 
489:     for (i=0; i<ndim; i++) {
490:       if (y[i] != zero) {
491:         MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES); 
492:       }
493:     }
494:     VecRestoreArray(user->y,&y); 

496:   }

498:   *flg=SAME_NONZERO_PATTERN;

500:   /* Assemble matrix  */
501:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); 
502:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); 
503:   return 0;
504: }


507: /* ------------------------------------------------------------------- */
510: /* 
511:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
512:    products.
513:     
514:    Input Parameters:
515: .  tao - the TaoSolver context
516: .  X    - the input vector
517: .  ptr  - optional user-defined context, as set by TaoSetHessian()
518:    
519:    Output Parameters:
520: .  H     - Hessian matrix
521: .  PrecH - optionally different preconditioning Hessian
522: .  flag  - flag indicating matrix structure
523: */
524: PetscErrorCode MatrixFreeHessian(TaoSolver tao,Vec X,Mat *H,Mat *PrecH,
525:                       MatStructure *flag,void *ptr)
526: {
527:   AppCtx     *user = (AppCtx *) ptr;

529:   /* Sets location of vector for use in computing matrix-vector products
530:      of the form H(X)*y  */

532:   user->xvec = X;   
533:   return 0;
534: }

536: /* ------------------------------------------------------------------- */
539: /* 
540:    HessianProductMat - Computes the matrix-vector product
541:    y = mat*svec.

543:    Input Parameters:
544: .  mat  - input matrix
545: .  svec - input vector

547:    Output Parameters:
548: .  y    - solution vector
549: */
550: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
551: {
552:   void *ptr;
554:   MatShellGetContext(mat,&ptr); 
555:   HessianProduct(ptr,svec,y); 
556:   

558:   return 0;
559: }

561: /* ------------------------------------------------------------------- */
564: /* 
565:    Hessian Product - Computes the matrix-vector product: 
566:    y = f''(x)*svec.

568:    Input Parameters
569: .  ptr  - pointer to the user-defined context
570: .  svec - input vector

572:    Output Parameters:
573: .  y    - product vector
574: */
575: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
576: {
577:   AppCtx *user = (AppCtx *)ptr;
578:   PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
579:   PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
581:   PetscInt nx, ny, i, j, k, ind;

583:   nx   = user->mx;
584:   ny   = user->my;
585:   hx   = user->hx;
586:   hy   = user->hy;
587:   hxhx = one/(hx*hx);
588:   hyhy = one/(hy*hy);

590:   /* Get pointers to vector data */
591:   VecGetArray(user->xvec,&x); 
592:   VecGetArray(svec,&s); 

594:   /* Initialize product vector to zero */
595:   VecSet(y, zero); 

597:   /* Compute f''(x)*s over the lower triangular elements */
598:   for (j=-1; j<ny; j++) {
599:     for (i=-1; i<nx; i++) {
600:        k = nx*j + i;
601:        v = zero;
602:        vr = zero;
603:        vt = zero;
604:        if (i != -1 && j != -1) v = s[k];
605:        if (i != nx-1 && j != -1) {
606:          vr = s[k+1];
607:          ind = k+1; val = hxhx*(vr-v);
608:          VecSetValues(y,1,&ind,&val,ADD_VALUES); 
609:        }
610:        if (i != -1 && j != ny-1) {
611:          vt = s[k+nx];
612:          ind = k+nx; val = hyhy*(vt-v);
613:          VecSetValues(y,1,&ind,&val,ADD_VALUES); 
614:        }
615:        if (i != -1 && j != -1) {
616:          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
617:          VecSetValues(y,1,&ind,&val,ADD_VALUES); 
618:        }
619:      }
620:    }
621:   
622:   /* Compute f''(x)*s over the upper triangular elements */
623:   for (j=0; j<=ny; j++) {
624:     for (i=0; i<=nx; i++) {
625:        k = nx*j + i;
626:        v = zero;
627:        vl = zero;
628:        vb = zero;
629:        if (i != nx && j != ny) v = s[k];
630:        if (i != nx && j != 0) {
631:          vb = s[k-nx];
632:          ind = k-nx; val = hyhy*(vb-v);
633:          VecSetValues(y,1,&ind,&val,ADD_VALUES); 
634:        }
635:        if (i != 0 && j != ny) {
636:          vl = s[k-1];
637:          ind = k-1; val = hxhx*(vl-v);
638:          VecSetValues(y,1,&ind,&val,ADD_VALUES); 
639:        }
640:        if (i != nx && j != ny) {
641:          ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
642:          VecSetValues(y,1,&ind,&val,ADD_VALUES); 
643:        }
644:     }
645:   }
646:   /* Restore vector data */
647:   VecRestoreArray(svec,&s); 
648:   VecRestoreArray(user->xvec,&x); 

650:   /* Assemble vector */
651:   VecAssemblyBegin(y); 
652:   VecAssemblyEnd(y); 
653:  
654:   /* Scale resulting vector by area */
655:   area = p5*hx*hy;
656:   VecScale(y, area); 

658:   PetscLogFlops(nx*ny*18); 
659:   
660:   return 0;
661: }