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*/
 17: /* ------------------------------------------------------------------------

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

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

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

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

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

 51:   ------------------------------------------------------------------------- */

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

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

 74: PetscErrorCode FormInitialGuess(DMMG,Vec);
 75: PetscErrorCode FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
 76: PetscErrorCode FormFunctionLocali(DALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);
 77: PetscErrorCode FormFunctionLocali4(DALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);

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

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

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

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

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

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

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

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

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

137:        Process adiC(36): FormFunctionLocal FormFunctionLocali
138:        Process blockadiC(4): FormFunctionLocali4
139:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
141:     DMMGSetFromOptions(dmmg);
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: /********  PetscDraw draw;
178:   PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),0,&draw);
179:   PetscDrawSetPause(draw,-1);
180:   PetscDrawPause(draw); */
181:   PetscFinalize();
182:   return 0;
183: }

185: /* ------------------------------------------------------------------- */


190: /* 
191:    FormInitialGuess - Forms initial approximation.

193:    Input Parameters:
194:    user - user-defined application context
195:    X - vector

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

209:   grashof = user->grashof;

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

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

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

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

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

260:   grashof = user->grashof;
261:   prandtl = user->prandtl;
262:   lid     = user->lidvelocity;

264:   /* 
265:      Define mesh intervals ratios for uniform grid.

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

270:      
271:   */
272:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
273:   hx = 1.0/dhx;                   hy = 1.0/dhy;
274:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

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

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

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

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

330:   /* Compute over the interior points */
331:   for (j=yints; j<yinte; j++) {
332:     for (i=xints; i<xinte; i++) {

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

342:         /* U velocity */
343:         u          = x[j][i].u;
344:         uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
345:         uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
346:         f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

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

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

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

377:   /*
378:      Flop count (multiply-adds are counted as 2 operations)
379:   */
380:   PetscLogFlops(84.0*info->ym*info->xm);
381:   return(0);
382: }

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

397:   grashof = user->grashof;
398:   prandtl = user->prandtl;
399:   lid     = user->lidvelocity;

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

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

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

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

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

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

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

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

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

489:   return(0);
490: }

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

506:   grashof = user->grashof;
507:   prandtl = user->prandtl;
508:   lid     = user->lidvelocity;

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

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

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

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

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

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

559:     /* U velocity */
560:     u          = x[j][i].u;
561:     uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
562:     uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
563:     f->u        = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

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