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: }