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: /*T
 12:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
 13:    Concepts: DA^using distributed arrays;
 14:    Concepts: multicomponent
 15:    Processors: n
 16: T*/

 18: /* ------------------------------------------------------------------------

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

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

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

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

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

 52:   ------------------------------------------------------------------------- */

 54: /* 
 55:    Include "petscda.h" so that we can use distributed arrays (DAs).
 56:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 57:    file automatically includes:
 58:      petsc.h       - base PETSc routines   petscvec.h - vectors
 59:      petscsys.h    - system routines       petscmat.h - matrices
 60:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 61:      petscviewer.h - viewers               petscpc.h  - preconditioners
 62:      petscksp.h   - linear solvers 
 63: */
 64:  #include petscsnes.h
 65:  #include petscda.h
 66:  #include petscdmmg.h

 68: /* 
 69:    User-defined routines and data structures
 70: */
 71: typedef struct {
 72:   PetscScalar u,v,omega,temp;
 73: } Field;


 80: typedef struct {
 81:    PassiveReal  lidvelocity,prandtl,grashof;  /* physical parameters */
 82:    PetscTruth     draw_contours;                /* flag - 1 indicates drawing contours */
 83: } AppCtx;

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

 97:   PetscInitialize(&argc,&argv,(char *)0,help);
 98:   comm = PETSC_COMM_WORLD;

100:   PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,PETSC_NULL);
101:   PreLoadBegin(PETSC_TRUE,"SetUp");
102:     DMMGCreate(comm,nlevels,&user,&dmmg);

104:     /*
105:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
106:       for principal unknowns (x) and governing residuals (f)
107:     */
108:     DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
109:     DMMGSetDM(dmmg,(DM)da);
110:     DADestroy(da);

112:     DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
113:                      PETSC_IGNORE,PETSC_IGNORE);
114:     /* 
115:      Problem parameters (velocity of lid, prandtl, and grashof numbers)
116:     */
117:     user.lidvelocity = 1.0/(mx*my);
118:     user.prandtl     = 1.0;
119:     user.grashof     = 1.0;
120:     PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
121:     PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
122:     PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);
123:     PetscOptionsHasName(PETSC_NULL,"-contours",&user.draw_contours);

125:     DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
126:     DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
127:     DASetFieldName(DMMGGetDA(dmmg),2,"Omega");
128:     DASetFieldName(DMMGGetDA(dmmg),3,"temperature");

130:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131:        Create user context, set problem data, create vector data structures.
132:        Also, compute the initial guess.
133:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

135:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:        Create nonlinear solver context

138:        Process adiC(36): FormFunctionLocal FormFunctionLocali
139:        Process blockadiC(4): FormFunctionLocali4
140:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
142:     DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,admf_FormFunctionLocali);
143:     DMMGSetSNESLocalib(dmmg,FormFunctionLocali4,0,admfb_FormFunctionLocali4);

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


149:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:        Solve the nonlinear system
151:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152:     DMMGSetInitialGuess(dmmg,FormInitialGuess);

154:   PreLoadStage("Solve");
155:     DMMGSolve(dmmg);

157:     snes = DMMGGetSNES(dmmg);
158:     SNESGetIterationNumber(snes,&its);
159:     PetscPrintf(comm,"Number of Newton iterations = %D\n", its);

161:     /*
162:       Visualize solution
163:     */

165:     if (user.draw_contours) {
166:       VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
167:     }

169:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:        Free work space.  All PETSc objects should be destroyed when they
171:        are no longer needed.
172:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

174:     DMMGDestroy(dmmg);
175:   PreLoadEnd();

177:   PetscFinalize();
178:   return 0;
179: }

181: /* ------------------------------------------------------------------- */


186: /* 
187:    FormInitialGuess - Forms initial approximation.

189:    Input Parameters:
190:    user - user-defined application context
191:    X - vector

193:    Output Parameter:
194:    X - vector
195:  */
196: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
197: {
198:   AppCtx         *user = (AppCtx*)dmmg->user;
199:   DA             da = (DA)dmmg->dm;
200:   PetscInt       i,j,mx,xs,ys,xm,ym;
202:   PetscReal      grashof,dx;
203:   Field          **x;

205:   grashof = user->grashof;

207:   DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
208:   dx  = 1.0/(mx-1);

210:   /*
211:      Get local grid boundaries (for 2-dimensional DA):
212:        xs, ys   - starting grid indices (no ghost points)
213:        xm, ym   - widths of local grid (no ghost points)
214:   */
215:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

217:   /*
218:      Get a pointer to vector data.
219:        - For default PETSc vectors, VecGetArray() returns a pointer to
220:          the data array.  Otherwise, the routine is implementation dependent.
221:        - You MUST call VecRestoreArray() when you no longer need access to
222:          the array.
223:   */
224:   DAVecGetArray(da,X,&x);

226:   /*
227:      Compute initial guess over the locally owned part of the grid
228:      Initial condition is motionless fluid and equilibrium temperature
229:   */
230:   for (j=ys; j<ys+ym; j++) {
231:     for (i=xs; i<xs+xm; i++) {
232:       x[j][i].u     = 0.0;
233:       x[j][i].v     = 0.0;
234:       x[j][i].omega = 0.0;
235:       x[j][i].temp  = (grashof>0)*i*dx;
236:     }
237:   }

239:   /*
240:      Restore vector
241:   */
242:   DAVecRestoreArray(da,X,&x);
243:   return 0;
244: }
245: 
246: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
247:  {
248:   AppCtx         *user = (AppCtx*)ptr;
250:   PetscInt       xints,xinte,yints,yinte,i,j;
251:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
252:   PetscReal      grashof,prandtl,lid;
253:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

256:   grashof = user->grashof;
257:   prandtl = user->prandtl;
258:   lid     = user->lidvelocity;

260:   /* 
261:      Define mesh intervals ratios for uniform grid.

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

266:      
267:   */
268:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
269:   hx = 1.0/dhx;                   hy = 1.0/dhy;
270:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

274:   /* Test whether we are on the bottom edge of the global array */
275:   if (yints == 0) {
276:     j = 0;
277:     yints = yints + 1;
278:     /* bottom edge */
279:     for (i=info->xs; i<info->xs+info->xm; i++) {
280:       f[j][i].u     = x[j][i].u;
281:       f[j][i].v     = x[j][i].v;
282:       f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
283:       f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
284:     }
285:   }

287:   /* Test whether we are on the top edge of the global array */
288:   if (yinte == info->my) {
289:     j = info->my - 1;
290:     yinte = yinte - 1;
291:     /* top edge */
292:     for (i=info->xs; i<info->xs+info->xm; i++) {
293:         f[j][i].u     = x[j][i].u - lid;
294:         f[j][i].v     = x[j][i].v;
295:         f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
296:         f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
297:     }
298:   }

300:   /* Test whether we are on the left edge of the global array */
301:   if (xints == 0) {
302:     i = 0;
303:     xints = xints + 1;
304:     /* left edge */
305:     for (j=info->ys; j<info->ys+info->ym; j++) {
306:       f[j][i].u     = x[j][i].u;
307:       f[j][i].v     = x[j][i].v;
308:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
309:       f[j][i].temp  = x[j][i].temp;
310:     }
311:   }

313:   /* Test whether we are on the right edge of the global array */
314:   if (xinte == info->mx) {
315:     i = info->mx - 1;
316:     xinte = xinte - 1;
317:     /* right edge */
318:     for (j=info->ys; j<info->ys+info->ym; j++) {
319:       f[j][i].u     = x[j][i].u;
320:       f[j][i].v     = x[j][i].v;
321:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
322:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
323:     }
324:   }

326:   /* Compute over the interior points */
327:   for (j=yints; j<yinte; j++) {
328:     for (i=xints; i<xinte; i++) {

330:         /*
331:           convective coefficients for upwinding
332:         */
333:         vx = x[j][i].u; avx = PetscAbsScalar(vx);
334:         vxp = .5*(vx+avx); vxm = .5*(vx-avx);
335:         vy = x[j][i].v; avy = PetscAbsScalar(vy);
336:         vyp = .5*(vy+avy); vym = .5*(vy-avy);

338:         /* U velocity */
339:         u          = x[j][i].u;
340:         uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
341:         uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
342:         f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

344:         /* V velocity */
345:         u          = x[j][i].v;
346:         uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
347:         uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
348:         f[j][i].v  = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

350:         /* Omega */
351:         u          = x[j][i].omega;
352:         uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
353:         uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
354:         f[j][i].omega = uxx + uyy +
355:                         (vxp*(u - x[j][i-1].omega) +
356:                           vxm*(x[j][i+1].omega - u)) * hy +
357:                         (vyp*(u - x[j-1][i].omega) +
358:                           vym*(x[j+1][i].omega - u)) * hx -
359:                         .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

361:         /* Temperature */
362:         u             = x[j][i].temp;
363:         uxx           = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
364:         uyy           = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
365:         f[j][i].temp =  uxx + uyy  + prandtl * (
366:                         (vxp*(u - x[j][i-1].temp) +
367:                           vxm*(x[j][i+1].temp - u)) * hy +
368:                         (vyp*(u - x[j-1][i].temp) +
369:                                  vym*(x[j+1][i].temp - u)) * hx);
370:     }
371:   }

373:   /*
374:      Flop count (multiply-adds are counted as 2 operations)
375:   */
376:   PetscLogFlops(84*info->ym*info->xm);
377:   return(0);
378: }

380: /*
381:     This function that evaluates the function for a single 
382:     degree of freedom. It is used by the -dmmg_fas solver
383: */
384: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
385: {
386:   AppCtx      *user = (AppCtx*)ptr;
387:   PetscInt    i,j,c;
388:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
389:   PassiveReal grashof,prandtl,lid;
390:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

393:   grashof = user->grashof;
394:   prandtl = user->prandtl;
395:   lid     = user->lidvelocity;

397:   /* 
398:      Define mesh intervals ratios for uniform grid.
399:      [Note: FD formulae below are normalized by multiplying through by
400:      local volume element to obtain coefficients O(1) in two dimensions.]
401:   */
402:   dhx = (PetscReal)(info->mx-1);     dhy = (PetscReal)(info->my-1);
403:   hx = 1.0/dhx;                   hy = 1.0/dhy;
404:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

408:   /* Test whether we are on the right edge of the global array */
409:   if (i == info->mx-1) {
410:     if (c == 0)      *f = x[j][i].u;
411:     else if (c == 1) *f = x[j][i].v;
412:     else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
413:     else             *f = x[j][i].temp - (PetscReal)(grashof>0);

415:   /* Test whether we are on the left edge of the global array */
416:   } else if (i == 0) {
417:     if (c == 0)      *f = x[j][i].u;
418:     else if (c == 1) *f = x[j][i].v;
419:     else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
420:     else             *f = x[j][i].temp;

422:   /* Test whether we are on the top edge of the global array */
423:   } else if (j == info->my-1) {
424:     if (c == 0)      *f = x[j][i].u - lid;
425:     else if (c == 1) *f = x[j][i].v;
426:     else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
427:     else             *f = x[j][i].temp-x[j-1][i].temp;

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

436:   /* Compute over the interior points */
437:   } else {
438:     /*
439:       convective coefficients for upwinding
440:     */
441:     vx = x[j][i].u; avx = PetscAbsScalar(vx);
442:     vxp = .5*(vx+avx); vxm = .5*(vx-avx);
443:     vy = x[j][i].v; avy = PetscAbsScalar(vy);
444:     vyp = .5*(vy+avy); vym = .5*(vy-avy);

446:     /* U velocity */
447:     if (c == 0) {
448:       u          = x[j][i].u;
449:       uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
450:       uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
451:       *f         = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

453:     /* V velocity */
454:     } else if (c == 1) {
455:       u          = x[j][i].v;
456:       uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
457:       uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
458:       *f         = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
459: 
460:     /* Omega */
461:     } else if (c == 2) {
462:       u          = x[j][i].omega;
463:       uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
464:       uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
465:       *f         = uxx + uyy +
466:         (vxp*(u - x[j][i-1].omega) +
467:          vxm*(x[j][i+1].omega - u)) * hy +
468:         (vyp*(u - x[j-1][i].omega) +
469:          vym*(x[j+1][i].omega - u)) * hx -
470:          .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
471: 
472:     /* Temperature */
473:     } else {
474:       u           = x[j][i].temp;
475:       uxx         = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
476:       uyy         = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
477:       *f          =  uxx + uyy  + prandtl * (
478:                                              (vxp*(u - x[j][i-1].temp) +
479:                                               vxm*(x[j][i+1].temp - u)) * hy +
480:                                              (vyp*(u - x[j-1][i].temp) +
481:                                              vym*(x[j+1][i].temp - u)) * hx);
482:     }
483:   }

485:   return(0);
486: }

488: /*
489:     This function that evaluates the function for a single 
490:     grid point. It is used by the -dmmg_fas -dmmg_fas_block solver
491: */
492: PetscErrorCode FormFunctionLocali4(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *ff,void *ptr)
493: {
494:   Field       *f = (Field*)ff;
495:   AppCtx      *user = (AppCtx*)ptr;
496:   PetscInt    i,j;
497:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
498:   PassiveReal grashof,prandtl,lid;
499:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

502:   grashof = user->grashof;
503:   prandtl = user->prandtl;
504:   lid     = user->lidvelocity;

506:   /* 
507:      Define mesh intervals ratios for uniform grid.
508:      [Note: FD formulae below are normalized by multiplying through by
509:      local volume element to obtain coefficients O(1) in two dimensions.]
510:   */
511:   dhx = (PetscReal)(info->mx-1);     dhy = (PetscReal)(info->my-1);
512:   hx = 1.0/dhx;                   hy = 1.0/dhy;
513:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

517:   /* Test whether we are on the right edge of the global array */
518:   if (i == info->mx-1) {
519:     f->u = x[j][i].u;
520:     f->v = x[j][i].v;
521:     f->omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
522:     f->temp = x[j][i].temp - (PetscReal)(grashof>0);

524:     /* Test whether we are on the left edge of the global array */
525:   } else if (i == 0) {
526:     f->u = x[j][i].u;
527:     f->v = x[j][i].v;
528:     f->omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
529:     f->temp = x[j][i].temp;
530: 
531:     /* Test whether we are on the top edge of the global array */
532:   } else if (j == info->my-1) {
533:     f->u = x[j][i].u - lid;
534:     f->v = x[j][i].v;
535:     f->omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
536:     f->temp = x[j][i].temp-x[j-1][i].temp;

538:     /* Test whether we are on the bottom edge of the global array */
539:   } else if (j == 0) {
540:     f->u = x[j][i].u;
541:     f->v = x[j][i].v;
542:     f->omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
543:     f->temp = x[j][i].temp - x[j+1][i].temp;

545:     /* Compute over the interior points */
546:   } else {
547:     /*
548:       convective coefficients for upwinding
549:     */
550:     vx = x[j][i].u; avx = PetscAbsScalar(vx);
551:     vxp = .5*(vx+avx); vxm = .5*(vx-avx);
552:     vy = x[j][i].v; avy = PetscAbsScalar(vy);
553:     vyp = .5*(vy+avy); vym = .5*(vy-avy);

555:     /* U velocity */
556:     u          = x[j][i].u;
557:     uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
558:     uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
559:     f->u        = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

561:     /* V velocity */
562:     u          = x[j][i].v;
563:     uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
564:     uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
565:     f->v        = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
566: 
567:     /* Omega */
568:     u          = x[j][i].omega;
569:     uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
570:     uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
571:     f->omega    = uxx + uyy +
572:       (vxp*(u - x[j][i-1].omega) +
573:        vxm*(x[j][i+1].omega - u)) * hy +
574:       (vyp*(u - x[j-1][i].omega) +
575:        vym*(x[j+1][i].omega - u)) * hx -
576:        .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
577: 
578:     /* Temperature */
579: 
580:     u           = x[j][i].temp;
581:     uxx         = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
582:     uyy         = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
583:     f->temp     =  uxx + uyy  + prandtl * (
584:                                            (vxp*(u - x[j][i-1].temp) +
585:                                             vxm*(x[j][i+1].temp - u)) * hy +
586:                                            (vyp*(u - x[j-1][i].temp) +
587:                                             vym*(x[j+1][i].temp - u)) * hx);
588:   }
589:   return(0);
590: }