Actual source code: ex19tu.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: DMDA^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 "petscdmda.h" so that we can use distributed arrays (DMDAs).
 56:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 57:    file automatically includes:
 58:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 59:      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 <petscdmda.h>
 66: #include <petscdmmg.h>
 67: #include <../src/snes/impls/ls/lsimpl.h>
 68: /* 
 69:    User-defined routines and data structures
 70: */
 71: typedef struct {
 72:   PetscScalar u,v,omega;
 73: } Field;


 80: typedef struct {
 81:   PassiveReal  lidvelocity,prandtl,grashof,re;  /* physical parameters */
 82:    PetscBool      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;
 93:   MPI_Comm       comm;
 94:   SNES           snes;
 95:   DM             da;

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


101:   // in order only run once, I block it: PetscPreLoadBegin(PETSC_TRUE,"SetUp");
102:     DMMGCreate(comm,2,&user,&dmmg);


105:     /*
106:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
107:       for principal unknowns (x) and governing residuals (f)
108:     */
109:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
110:     DMMGSetDM(dmmg,(DM)da);
111:     DMDestroy(&da);
112:     PetscPrintf(comm,"mx = %d, my= %d\n",
113:                        mx,my);

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

132:     DMDASetFieldName(DMMGGetDM(dmmg),0,"x-velocity");
133:     DMDASetFieldName(DMMGGetDM(dmmg),1,"y-velocity");
134:     DMDASetFieldName(DMMGGetDM(dmmg),2,"Omega");
135:     //  DMDASetFieldName(DMMGGetDM(dmmg),3,"temperature");

137:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138:        Create user context, set problem data, create vector data structures.
139:        Also, compute the initial guess.
140:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

142:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:        Create nonlinear solver context

145:        Process adiC(36): FormFunctionLocal FormFunctionLocali FormFunctionLocali4
146:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
148:     DMMGSetFromOptions(dmmg);
149:     DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
150:     DMMGSetSNESLocalib(dmmg,FormFunctionLocali4,ad_FormFunctionLocali4,admf_FormFunctionLocali4);

152:     PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g, re=%g \n",
153:                        user.lidvelocity,user.prandtl,user.grashof,user.re);


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

161:     //I block it:  PetscPreLoadStage("Solve");
162:     DMMGSolve(dmmg);

164:     snes = DMMGGetSNES(dmmg);
165:     SNESGetIterationNumber(snes,&its);
166:     PetscPrintf(comm,"Number of Newton iterations = %D\n", its);
167: 
168:     /*
169:       Visualize solution
170:     */

172:     if (user.draw_contours) {
173:       VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
174:       // VecView(DMMGGetx(dmmg),PETSC_VIEWER_STDOUT_WORLD);
175:     }

177:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:        Free work space.  All PETSc objects should be destroyed when they
179:        are no longer needed.
180:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

182:     DMMGDestroy(dmmg);
183:     //PetscPreLoadEnd();

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: }
253: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
254:  {
255:   AppCtx         *user = (AppCtx*)ptr;
257:   PetscInt       xints,xinte,yints,yinte,i,j;
258:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
259:   PetscReal      grashof,prandtl,lid,re;
260:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

263:   grashof = user->grashof;
264:   prandtl = user->prandtl;
265:   lid     = user->lidvelocity;
266:   re      = user->re;
267:   /* 
268:      Define mesh intervals ratios for uniform grid.
269:      [Note: FD formulae below are normalized by multiplying through by
270:      local volume element to obtain coefficients O(1) in two dimensions.]
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)/re +
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 is an experimental function and can be safely ignored.
386: */
387: PetscErrorCode FormFunctionLocali(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
388:  {
389:   AppCtx      *user = (AppCtx*)ptr;
390:   PetscInt    i,j,c;
391:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
392:   PassiveReal grashof,prandtl,lid,re;
393:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

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

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

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

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

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

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

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

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

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

488:   return(0);
489: }

491: /*
492:     This is an experimental function and can be safely ignored.
493: */
494: PetscErrorCode FormFunctionLocali4(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *ff,void *ptr)
495:  {
496:    Field *f = (Field*)ff;
497:   AppCtx      *user = (AppCtx*)ptr;
498:   PetscInt    i,j;
499:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
500:   PassiveReal grashof,prandtl,lid,re;
501:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

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

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

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

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

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

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

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

563:     /* V velocity */
564:     u          = x[j][i].v;
565:     uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
566:     uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
567:     f->v        = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
568: 
569:     /* Omega */
570:     u          = x[j][i].omega;
571:     uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
572:     uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
573:     f->omega    = uxx + uyy +
574:                     (vx*(u - x[j][i-1].omega)) * hy +
575:                     (vy*(u - x[j-1][i].omega)) * hx ;-
576:                          .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

578:     u           = x[j][i].temp;
579:     uxx         = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
580:     uyy         = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
581:     f->temp     =  uxx + uyy  + prandtl * (
582:                                            (vxp*(u - x[j][i-1].temp) +
583:                                             vxm*(x[j][i+1].temp - u)) * hy +
584:                                            (vyp*(u - x[j-1][i].temp) +
585:                                            vym*(x[j+1][i].temp - u)) * hx);
586:   }

588:   return(0);
589: }