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