Actual source code: ex19.c

  2: static char help[] = "Nonlinear driven cavity with multigrid in 2d.\n\
  3:   \n\
  4: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
  5: The flow can be driven with the lid or with bouyancy or both:\n\
  6:   -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
  7:   -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
  8:   -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
  9:   -contours : draw contour plots of solution\n\n";

 11: /*
 12:       THIS EXAMPLE IS DEPRECATED, PLEASE SEE ex50.c FOR THE CURRENT APPROACH
 13:       DMMG operations are deprecated 
 14: */

 16: /*T
 17:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
 18:    Concepts: DMDA^using distributed arrays;
 19:    Concepts: multicomponent
 20:    Processors: n
 21: T*/
 22: /* ------------------------------------------------------------------------

 24:     We thank David E. Keyes for contributing the driven cavity discretization
 25:     within this example code.

 27:     This problem is modeled by the partial differential equation system
 28:   
 29:         - Lap(U) - Grad_y(Omega) = 0
 30:         - Lap(V) + Grad_x(Omega) = 0
 31:         - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
 32:         - Lap(T) + PR*Div([U*T,V*T]) = 0

 34:     in the unit square, which is uniformly discretized in each of x and
 35:     y in this simple encoding.

 37:     No-slip, rigid-wall Dirichlet conditions are used for [U,V].
 38:     Dirichlet conditions are used for Omega, based on the definition of
 39:     vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
 40:     constant coordinate boundary, the tangential derivative is zero.
 41:     Dirichlet conditions are used for T on the left and right walls,
 42:     and insulation homogeneous Neumann conditions are used for T on
 43:     the top and bottom walls. 

 45:     A finite difference approximation with the usual 5-point stencil 
 46:     is used to discretize the boundary value problem to obtain a 
 47:     nonlinear system of equations.  Upwinding is used for the divergence
 48:     (convective) terms and central for the gradient (source) terms.
 49:     
 50:     The Jacobian can be either
 51:       * formed via finite differencing using coloring (the default), or
 52:       * applied matrix-free via the option -snes_mf 
 53:         (for larger grid problems this variant may not converge 
 54:         without a preconditioner due to ill-conditioning).

 56:   ------------------------------------------------------------------------- */

 58: /* 
 59:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 60:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 61:    file automatically includes:
 62:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 63:      petscmat.h - matrices
 64:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 65:      petscviewer.h - viewers               petscpc.h  - preconditioners
 66:      petscksp.h   - linear solvers 
 67: */
 68: #include <petscsnes.h>
 69: #include <petscdmda.h>
 70: #include <petscdmmg.h>

 72: /* 
 73:    User-defined routines and data structures
 74: */
 75: typedef struct {
 76:   PetscScalar u,v,omega,temp;
 77: } Field;

 79: PetscErrorCode FormInitialGuess(DMMG,Vec);
 80: PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);
 81: PetscErrorCode FormFunctionLocali(DMDALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);
 82: PetscErrorCode FormFunctionLocali4(DMDALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);

 84: typedef struct {
 85:    PassiveReal  lidvelocity,prandtl,grashof;  /* physical parameters */
 86:    PetscBool      draw_contours;                /* flag - 1 indicates drawing contours */
 87: } AppCtx;

 91: int main(int argc,char **argv)
 92: {
 93:   DMMG           *dmmg;               /* multilevel grid structure */
 94:   AppCtx         user;                /* user-defined work context */
 95:   PetscInt       mx,my,its,nlevels=2;
 97:   MPI_Comm       comm;
 98:   SNES           snes;
 99:   DM             da;

101:   PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return(1);
102:   comm = PETSC_COMM_WORLD;

104:   PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,PETSC_NULL);
105:   PetscPreLoadBegin(PETSC_TRUE,"SetUp");
106:     DMMGCreate(comm,nlevels,&user,&dmmg);

108:     /*
109:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
110:       for principal unknowns (x) and governing residuals (f)
111:     */
112:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
113:     DMMGSetDM(dmmg,(DM)da);
114:     DMDestroy(&da);

116:     DMDAGetInfo(DMMGGetDM(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
117:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
118:     /* 
119:      Problem parameters (velocity of lid, prandtl, and grashof numbers)
120:     */
121:     user.lidvelocity = 1.0/(mx*my);
122:     user.prandtl     = 1.0;
123:     user.grashof     = 1.0;
124:     PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
125:     PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
126:     PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);
127:     PetscOptionsHasName(PETSC_NULL,"-contours",&user.draw_contours);

129:     DMDASetFieldName(DMMGGetDM(dmmg),0,"x-velocity");
130:     DMDASetFieldName(DMMGGetDM(dmmg),1,"y-velocity");
131:     DMDASetFieldName(DMMGGetDM(dmmg),2,"Omega");
132:     DMDASetFieldName(DMMGGetDM(dmmg),3,"temperature");

134:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:        Create user context, set problem data, create vector data structures.
136:        Also, compute the initial guess.
137:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

139:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140:        Create nonlinear solver context

142:        Process adiC(36): FormFunctionLocal FormFunctionLocali
143:        Process blockadiC(4): FormFunctionLocali4
144:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
146:     DMMGSetFromOptions(dmmg);
147:     /*DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,admf_FormFunctionLocali);
148:       DMMGSetSNESLocalib(dmmg,FormFunctionLocali4,0,admfb_FormFunctionLocali4); */

150:     PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, grashof # = %G\n",
151:                        user.lidvelocity,user.prandtl,user.grashof);


154:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:        Solve the nonlinear system
156:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157:     DMMGSetInitialGuess(dmmg,FormInitialGuess);

159:   PetscPreLoadStage("Solve");
160:     DMMGSolve(dmmg);
161:     snes = DMMGGetSNES(dmmg);
162:     SNESGetIterationNumber(snes,&its);
163:     PetscPrintf(comm,"Number of Newton iterations = %D\n", its);

165:     /*
166:       Visualize solution
167:     */

169:     if (user.draw_contours) {
170:       VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
171:     }

173:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174:        Free work space.  All PETSc objects should be destroyed when they
175:        are no longer needed.
176:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

178:     DMMGDestroy(dmmg);
179:   PetscPreLoadEnd();

181: /********  PetscDraw draw;
182:   PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),0,&draw);
183:   PetscDrawSetPause(draw,-1);
184:   PetscDrawPause(draw); */
185:   PetscFinalize();
186:   return 0;
187: }

189: /* ------------------------------------------------------------------- */


194: /* 
195:    FormInitialGuess - Forms initial approximation.

197:    Input Parameters:
198:    user - user-defined application context
199:    X - vector

201:    Output Parameter:
202:    X - vector
203:  */
204: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
205: {
206:   AppCtx         *user = (AppCtx*)dmmg->user;
207:   DM             da = dmmg->dm;
208:   PetscInt       i,j,mx,xs,ys,xm,ym;
210:   PetscReal      grashof,dx;
211:   Field          **x;

213:   grashof = user->grashof;

215:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
216:   dx  = 1.0/(mx-1);

218:   /*
219:      Get local grid boundaries (for 2-dimensional DMDA):
220:        xs, ys   - starting grid indices (no ghost points)
221:        xm, ym   - widths of local grid (no ghost points)
222:   */
223:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

225:   /*
226:      Get a pointer to vector data.
227:        - For default PETSc vectors, VecGetArray() returns a pointer to
228:          the data array.  Otherwise, the routine is implementation dependent.
229:        - You MUST call VecRestoreArray() when you no longer need access to
230:          the array.
231:   */
232:   DMDAVecGetArray(da,X,&x);

234:   /*
235:      Compute initial guess over the locally owned part of the grid
236:      Initial condition is motionless fluid and equilibrium temperature
237:   */
238:   for (j=ys; j<ys+ym; j++) {
239:     for (i=xs; i<xs+xm; i++) {
240:       x[j][i].u     = 0.0;
241:       x[j][i].v     = 0.0;
242:       x[j][i].omega = 0.0;
243:       x[j][i].temp  = (grashof>0)*i*dx;
244:     }
245:   }

247:   /*
248:      Restore vector
249:   */
250:   DMDAVecRestoreArray(da,X,&x);
251:   return 0;
252: }

256: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
257:  {
258:   AppCtx         *user = (AppCtx*)ptr;
260:   PetscInt       xints,xinte,yints,yinte,i,j;
261:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
262:   PetscReal      grashof,prandtl,lid;
263:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

266:   grashof = user->grashof;
267:   prandtl = user->prandtl;
268:   lid     = user->lidvelocity;

270:   /* 
271:      Define mesh intervals ratios for uniform grid.

273:      Note: FD formulae below are normalized by multiplying through by
274:      local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.

276:      
277:   */
278:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
279:   hx = 1.0/dhx;                   hy = 1.0/dhy;
280:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

282:   xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;

284:   /* Test whether we are on the bottom edge of the global array */
285:   if (yints == 0) {
286:     j = 0;
287:     yints = yints + 1;
288:     /* bottom edge */
289:     for (i=info->xs; i<info->xs+info->xm; i++) {
290:       f[j][i].u     = x[j][i].u;
291:       f[j][i].v     = x[j][i].v;
292:       f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
293:       f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
294:     }
295:   }

297:   /* Test whether we are on the top edge of the global array */
298:   if (yinte == info->my) {
299:     j = info->my - 1;
300:     yinte = yinte - 1;
301:     /* top edge */
302:     for (i=info->xs; i<info->xs+info->xm; i++) {
303:         f[j][i].u     = x[j][i].u - lid;
304:         f[j][i].v     = x[j][i].v;
305:         f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
306:         f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
307:     }
308:   }

310:   /* Test whether we are on the left edge of the global array */
311:   if (xints == 0) {
312:     i = 0;
313:     xints = xints + 1;
314:     /* left edge */
315:     for (j=info->ys; j<info->ys+info->ym; j++) {
316:       f[j][i].u     = x[j][i].u;
317:       f[j][i].v     = x[j][i].v;
318:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
319:       f[j][i].temp  = x[j][i].temp;
320:     }
321:   }

323:   /* Test whether we are on the right edge of the global array */
324:   if (xinte == info->mx) {
325:     i = info->mx - 1;
326:     xinte = xinte - 1;
327:     /* right edge */
328:     for (j=info->ys; j<info->ys+info->ym; j++) {
329:       f[j][i].u     = x[j][i].u;
330:       f[j][i].v     = x[j][i].v;
331:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
332:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
333:     }
334:   }

336:   /* Compute over the interior points */
337:   for (j=yints; j<yinte; j++) {
338:     for (i=xints; i<xinte; i++) {

340:         /*
341:           convective coefficients for upwinding
342:         */
343:         vx = x[j][i].u; avx = PetscAbsScalar(vx);
344:         vxp = .5*(vx+avx); vxm = .5*(vx-avx);
345:         vy = x[j][i].v; avy = PetscAbsScalar(vy);
346:         vyp = .5*(vy+avy); vym = .5*(vy-avy);

348:         /* U velocity */
349:         u          = x[j][i].u;
350:         uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
351:         uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
352:         f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

354:         /* V velocity */
355:         u          = x[j][i].v;
356:         uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
357:         uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
358:         f[j][i].v  = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

360:         /* Omega */
361:         u          = x[j][i].omega;
362:         uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
363:         uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
364:         f[j][i].omega = uxx + uyy +
365:                         (vxp*(u - x[j][i-1].omega) +
366:                           vxm*(x[j][i+1].omega - u)) * hy +
367:                         (vyp*(u - x[j-1][i].omega) +
368:                           vym*(x[j+1][i].omega - u)) * hx -
369:                         .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

371:         /* Temperature */
372:         u             = x[j][i].temp;
373:         uxx           = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
374:         uyy           = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
375:         f[j][i].temp =  uxx + uyy  + prandtl * (
376:                         (vxp*(u - x[j][i-1].temp) +
377:                           vxm*(x[j][i+1].temp - u)) * hy +
378:                         (vyp*(u - x[j-1][i].temp) +
379:                                  vym*(x[j+1][i].temp - u)) * hx);
380:     }
381:   }

383:   /*
384:      Flop count (multiply-adds are counted as 2 operations)
385:   */
386:   PetscLogFlops(84.0*info->ym*info->xm);
387:   return(0);
388: }

392: /*
393:     This function that evaluates the function for a single 
394:     degree of freedom. It is used by the -dmmg_fas solver
395: */
396: PetscErrorCode FormFunctionLocali(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
397: {
398:   AppCtx      *user = (AppCtx*)ptr;
399:   PetscInt    i,j,c;
400:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
401:   PassiveReal grashof,prandtl,lid;
402:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

405:   grashof = user->grashof;
406:   prandtl = user->prandtl;
407:   lid     = user->lidvelocity;

409:   /* 
410:      Define mesh intervals ratios for uniform grid.
411:      [Note: FD formulae below are normalized by multiplying through by
412:      local volume element to obtain coefficients O(1) in two dimensions.]
413:   */
414:   dhx = (PetscReal)(info->mx-1);     dhy = (PetscReal)(info->my-1);
415:   hx = 1.0/dhx;                   hy = 1.0/dhy;
416:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

418:   i = st->i; j = st->j; c = st->c;

420:   /* Test whether we are on the right edge of the global array */
421:   if (i == info->mx-1) {
422:     if (c == 0)      *f = x[j][i].u;
423:     else if (c == 1) *f = x[j][i].v;
424:     else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
425:     else             *f = x[j][i].temp - (PetscReal)(grashof>0);

427:   /* Test whether we are on the left edge of the global array */
428:   } else if (i == 0) {
429:     if (c == 0)      *f = x[j][i].u;
430:     else if (c == 1) *f = x[j][i].v;
431:     else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
432:     else             *f = x[j][i].temp;

434:   /* Test whether we are on the top edge of the global array */
435:   } else if (j == info->my-1) {
436:     if (c == 0)      *f = x[j][i].u - lid;
437:     else if (c == 1) *f = x[j][i].v;
438:     else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
439:     else             *f = x[j][i].temp-x[j-1][i].temp;

441:   /* Test whether we are on the bottom edge of the global array */
442:   } else if (j == 0) {
443:     if (c == 0)      *f = x[j][i].u;
444:     else if (c == 1) *f = x[j][i].v;
445:     else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
446:     else             *f = x[j][i].temp - x[j+1][i].temp;

448:   /* Compute over the interior points */
449:   } else {
450:     /*
451:       convective coefficients for upwinding
452:     */
453:     vx = x[j][i].u; avx = PetscAbsScalar(vx);
454:     vxp = .5*(vx+avx); vxm = .5*(vx-avx);
455:     vy = x[j][i].v; avy = PetscAbsScalar(vy);
456:     vyp = .5*(vy+avy); vym = .5*(vy-avy);

458:     /* U velocity */
459:     if (c == 0) {
460:       u          = x[j][i].u;
461:       uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
462:       uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
463:       *f         = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

465:     /* V velocity */
466:     } else if (c == 1) {
467:       u          = x[j][i].v;
468:       uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
469:       uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
470:       *f         = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
471: 
472:     /* Omega */
473:     } else if (c == 2) {
474:       u          = x[j][i].omega;
475:       uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
476:       uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
477:       *f         = uxx + uyy +
478:         (vxp*(u - x[j][i-1].omega) +
479:          vxm*(x[j][i+1].omega - u)) * hy +
480:         (vyp*(u - x[j-1][i].omega) +
481:          vym*(x[j+1][i].omega - u)) * hx -
482:          .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
483: 
484:     /* Temperature */
485:     } else {
486:       u           = x[j][i].temp;
487:       uxx         = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
488:       uyy         = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
489:       *f          =  uxx + uyy  + prandtl * (
490:                                              (vxp*(u - x[j][i-1].temp) +
491:                                               vxm*(x[j][i+1].temp - u)) * hy +
492:                                              (vyp*(u - x[j-1][i].temp) +
493:                                              vym*(x[j+1][i].temp - u)) * hx);
494:     }
495:   }

497:   return(0);
498: }

502: /*
503:     This function that evaluates the function for a single 
504:     grid point. It is used by the -dmmg_fas -dmmg_fas_block solver
505: */
506: PetscErrorCode FormFunctionLocali4(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *ff,void *ptr)
507: {
508:   Field       *f = (Field*)ff;
509:   AppCtx      *user = (AppCtx*)ptr;
510:   PetscInt    i,j;
511:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
512:   PassiveReal grashof,prandtl,lid;
513:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

516:   grashof = user->grashof;
517:   prandtl = user->prandtl;
518:   lid     = user->lidvelocity;

520:   /* 
521:      Define mesh intervals ratios for uniform grid.
522:      [Note: FD formulae below are normalized by multiplying through by
523:      local volume element to obtain coefficients O(1) in two dimensions.]
524:   */
525:   dhx = (PetscReal)(info->mx-1);     dhy = (PetscReal)(info->my-1);
526:   hx = 1.0/dhx;                   hy = 1.0/dhy;
527:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

529:   i = st->i; j = st->j;

531:   /* Test whether we are on the right edge of the global array */
532:   if (i == info->mx-1) {
533:     f->u = x[j][i].u;
534:     f->v = x[j][i].v;
535:     f->omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
536:     f->temp = x[j][i].temp - (PetscReal)(grashof>0);

538:     /* Test whether we are on the left edge of the global array */
539:   } else if (i == 0) {
540:     f->u = x[j][i].u;
541:     f->v = x[j][i].v;
542:     f->omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
543:     f->temp = x[j][i].temp;
544: 
545:     /* Test whether we are on the top edge of the global array */
546:   } else if (j == info->my-1) {
547:     f->u = x[j][i].u - lid;
548:     f->v = x[j][i].v;
549:     f->omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
550:     f->temp = x[j][i].temp-x[j-1][i].temp;

552:     /* Test whether we are on the bottom edge of the global array */
553:   } else if (j == 0) {
554:     f->u = x[j][i].u;
555:     f->v = x[j][i].v;
556:     f->omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
557:     f->temp = x[j][i].temp - x[j+1][i].temp;

559:     /* Compute over the interior points */
560:   } else {
561:     /*
562:       convective coefficients for upwinding
563:     */
564:     vx = x[j][i].u; avx = PetscAbsScalar(vx);
565:     vxp = .5*(vx+avx); vxm = .5*(vx-avx);
566:     vy = x[j][i].v; avy = PetscAbsScalar(vy);
567:     vyp = .5*(vy+avy); vym = .5*(vy-avy);

569:     /* U velocity */
570:     u          = x[j][i].u;
571:     uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
572:     uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
573:     f->u        = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

575:     /* V velocity */
576:     u          = x[j][i].v;
577:     uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
578:     uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
579:     f->v        = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
580: 
581:     /* Omega */
582:     u          = x[j][i].omega;
583:     uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
584:     uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
585:     f->omega    = uxx + uyy +
586:       (vxp*(u - x[j][i-1].omega) +
587:        vxm*(x[j][i+1].omega - u)) * hy +
588:       (vyp*(u - x[j-1][i].omega) +
589:        vym*(x[j+1][i].omega - u)) * hx -
590:        .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
591: 
592:     /* Temperature */
593: 
594:     u           = x[j][i].temp;
595:     uxx         = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
596:     uyy         = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
597:     f->temp     =  uxx + uyy  + prandtl * (
598:                                            (vxp*(u - x[j][i-1].temp) +
599:                                             vxm*(x[j][i+1].temp - u)) * hy +
600:                                            (vyp*(u - x[j-1][i].temp) +
601:                                             vym*(x[j+1][i].temp - u)) * hx);
602:   }
603:   return(0);
604: }