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: /*F-----------------------------------------------------------------------
We thank David E. Keyes for contributing the driven cavity discretization within this example code.
This problem is modeled by the partial differential equation system
in the unit square, which is uniformly discretized in each of x and y in this simple encoding.
No-slip, rigid-wall Dirichlet conditions are used for $ [U,V]$.
Dirichlet conditions are used for Omega, based on the definition of
vorticity: $ \Omega = - \nabla_y U + \nabla_x V$, where along each
constant coordinate boundary, the tangential derivative is zero.
Dirichlet conditions are used for T on the left and right walls,
and insulation homogeneous Neumann conditions are used for T on
the top and bottom walls.
A finite difference approximation with the usual 5-point stencil
is used to discretize the boundary value problem to obtain a
nonlinear system of equations. Upwinding is used for the divergence
(convective) terms and central for the gradient (source) terms.
The Jacobian can be either
* formed via finite differencing using coloring (the default), or
* applied matrix-free via the option -snes_mf
(for larger grid problems this variant may not converge
without a preconditioner due to ill-conditioning).
56: ------------------------------------------------------------------------F*/
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 SNES 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: }