Actual source code: ex27.c
1: /*$Id: ex19.c,v 1.30 2001/08/07 21:31:17 bsmith Exp $*/
3: static char help[] = "Nonlinear driven cavity with multigrid and pusedo timestepping 2d.\n\
4: \n\
5: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
6: The flow can be driven with the lid or with bouyancy or both:\n\
7: -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
8: -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
9: -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
10: -contours : draw contour plots of solution\n\n";
12: /*T
13: Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
14: Concepts: DA^using distributed arrays;
15: Concepts: multicomponent
16: Processors: n
17: T*/
19: /* ------------------------------------------------------------------------
21: We thank David E. Keyes for contributing the driven cavity discretization
22: within this example code.
24: This problem is modeled by the partial differential equation system
25:
26: dU/dt - Lap(U) - Grad_y(Omega) = 0
27: dV/dt - Lap(V) + Grad_x(Omega) = 0
28: d(omega)/dt - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
29: dT/dt - Lap(T) + PR*Div([U*T,V*T]) = 0
31: in the unit square, which is uniformly discretized in each of x and
32: y in this simple encoding.
34: No-slip, rigid-wall Dirichlet conditions are used for [U,V].
35: Dirichlet conditions are used for Omega, based on the definition of
36: vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
37: constant coordinate boundary, the tangential derivative is zero.
38: Dirichlet conditions are used for T on the left and right walls,
39: and insulation homogeneous Neumann conditions are used for T on
40: the top and bottom walls.
42: A finite difference approximation with the usual 5-point stencil
43: is used to discretize the boundary value problem to obtain a
44: nonlinear system of equations. Upwinding is used for the divergence
45: (convective) terms and central for the gradient (source) terms.
46:
47: The Jacobian can be either
48: * formed via finite differencing using coloring (the default), or
49: * applied matrix-free via the option -snes_mf
50: (for larger grid problems this variant may not converge
51: without a preconditioner due to ill-conditioning).
53: ------------------------------------------------------------------------- */
55: /*
56: Include "petscda.h" so that we can use distributed arrays (DAs).
57: Include "petscsnes.h" so that we can use SNES solvers. Note that this
58: file automatically includes:
59: petsc.h - base PETSc routines petscvec.h - vectors
60: petscsys.h - system routines petscmat.h - matrices
61: petscis.h - index sets petscksp.h - Krylov subspace methods
62: petscviewer.h - viewers petscpc.h - preconditioners
63: petscksp.h - linear solvers
64: */
65: #include petscsnes.h
66: #include petscda.h
68: /*
69: User-defined routines and data structures
70: */
72: typedef struct {
73: PassiveScalar fnorm_ini,dt_ini,cfl_ini;
74: PassiveScalar fnorm,dt,cfl;
75: PassiveScalar ptime;
76: PassiveScalar cfl_max,max_time;
77: PassiveScalar fnorm_ratio;
78: int ires,iramp,itstep;
79: int max_steps,print_freq;
80: int LocalTimeStepping;
81: PetscTruth use_parabolic;
82: } TstepCtx;
84: /*
85: The next two structures are essentially the same. The difference is that
86: the first does not contain derivative information (as used by ADIC) while the
87: second does. The first is used only to contain the solution at the previous time-step
88: which is a constant in the computation of the current time step and hence passive
89: (meaning not active in terms of derivatives).
90: */
91: typedef struct {
92: PassiveScalar u,v,omega,temp;
93: } PassiveField;
95: typedef struct {
96: PetscScalar u,v,omega,temp;
97: } Field;
100: typedef struct {
101: int mglevels;
102: int cycles; /* numbers of time steps for integration */
103: PassiveReal lidvelocity,prandtl,grashof; /* physical parameters */
104: PetscTruth draw_contours; /* indicates drawing contours of solution */
105: PetscTruth PreLoading;
106: } Parameter;
108: typedef struct {
109: Vec Xold,func;
110: TstepCtx *tsCtx;
111: Parameter *param;
112: } AppCtx;
114: extern int FormInitialGuess(SNES,Vec,void*);
115: extern int FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
116: extern int FormFunctionLocali(DALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);
117: extern int Update(DMMG *);
118: extern int Initialize(DMMG *);
119: extern int ComputeTimeStep(SNES,void*);
120: extern int AddTSTermLocal(DALocalInfo*,Field**,Field**,void*);
125: int main(int argc,char **argv)
126: {
127: DMMG *dmmg; /* multilevel grid structure */
128: AppCtx *user; /* user-defined work context */
129: TstepCtx tsCtx;
130: Parameter param;
131: int mx,my,its;
132: int i,ierr;
133: MPI_Comm comm;
134: SNES snes;
135: DA da;
137: PetscInitialize(&argc,&argv,(char *)0,help);
138: comm = PETSC_COMM_WORLD;
141: PreLoadBegin(PETSC_TRUE,"SetUp");
143: param.PreLoading = PreLoading;
144: DMMGCreate(comm,2,&user,&dmmg);
145: param.mglevels = DMMGGetLevels(dmmg);
148: /*
149: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
150: for principal unknowns (x) and governing residuals (f)
151: */
152: DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
153: DMMGSetDM(dmmg,(DM)da);
154: DADestroy(da);
156: DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
157: PETSC_IGNORE,PETSC_IGNORE);
158: /*
159: Problem parameters (velocity of lid, prandtl, and grashof numbers)
160: */
161: param.lidvelocity = 1.0/(mx*my);
162: param.prandtl = 1.0;
163: param.grashof = 1.0;
164: PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",¶m.lidvelocity,PETSC_NULL);
165: PetscOptionsGetReal(PETSC_NULL,"-prandtl",¶m.prandtl,PETSC_NULL);
166: PetscOptionsGetReal(PETSC_NULL,"-grashof",¶m.grashof,PETSC_NULL);
167: PetscOptionsHasName(PETSC_NULL,"-contours",¶m.draw_contours);
169: DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
170: DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
171: DASetFieldName(DMMGGetDA(dmmg),2,"Omega");
172: DASetFieldName(DMMGGetDA(dmmg),3,"temperature");
174: /*======================================================================*/
175: /* Initilize stuff related to time stepping */
176: /*======================================================================*/
177: tsCtx.fnorm_ini = 0.0; tsCtx.cfl_ini = 50.0; tsCtx.cfl_max = 1.0e+06;
178: tsCtx.max_steps = 50; tsCtx.max_time = 1.0e+12; tsCtx.iramp = -50;
179: tsCtx.dt = 0.01; tsCtx.fnorm_ratio = 1.0e+10;
180: tsCtx.LocalTimeStepping = 0;
181: tsCtx.use_parabolic = PETSC_FALSE;
182: PetscOptionsGetLogical(PETSC_NULL,"-use_parabolic",&tsCtx.use_parabolic,PETSC_IGNORE);
183: PetscOptionsGetInt(PETSC_NULL,"-max_st",&tsCtx.max_steps,PETSC_NULL);
184: PetscOptionsGetReal(PETSC_NULL,"-ts_rtol",&tsCtx.fnorm_ratio,PETSC_NULL);
185: PetscOptionsGetReal(PETSC_NULL,"-cfl_ini",&tsCtx.cfl_ini,PETSC_NULL);
186: PetscOptionsGetReal(PETSC_NULL,"-cfl_max",&tsCtx.cfl_max,PETSC_NULL);
187: tsCtx.print_freq = tsCtx.max_steps;
188: PetscOptionsGetInt(PETSC_NULL,"-print_freq",&tsCtx.print_freq,PETSC_NULL);
189:
190: PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
191: for (i=0; i<param.mglevels; i++) {
192: VecDuplicate(dmmg[i]->x, &(user[i].Xold));
193: VecDuplicate(dmmg[i]->x, &(user[i].func));
194: user[i].tsCtx = &tsCtx;
195: user[i].param = ¶m;
196: dmmg[i]->user = &user[i];
197: }
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Create user context, set problem data, create vector data structures.
200: Also, compute the initial guess.
201: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
202:
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Create nonlinear solver context
205:
206: Process adiC(20): AddTSTermLocal FormFunctionLocal FormFunctionLocali
207: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
209: DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
210:
211: PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g\n",
212: param.lidvelocity,param.prandtl,param.grashof);
213:
214:
215: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: Solve the nonlinear system
217: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218: DMMGSetInitialGuess(dmmg,FormInitialGuess);
219:
220: PreLoadStage("Solve");
221: Initialize(dmmg);
222: Update(dmmg);
223: /*
224: Visualize solution
225: */
226:
227: if (param.draw_contours) {
228: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
229: }
230: /*VecView(DMMGGetx(dmmg),PETSC_VIEWER_STDOUT_WORLD);*/
231: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: Free work space. All PETSc objects should be destroyed when they
233: are no longer needed.
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235:
236: for (i=0; i<param.mglevels; i++) {
237: VecDestroy(user[i].Xold);
238: VecDestroy(user[i].func);
239: }
240: PetscFree(user);
241: DMMGDestroy(dmmg);
242: PreLoadEnd();
243:
244: PetscFinalize();
245: return 0;
246: }
248: /* ------------------------------------------------------------------- */
251: /* ------------------------------------------------------------------- */
252: int Initialize(DMMG *dmmg)
253: {
254: AppCtx *user = (AppCtx*)dmmg[0]->user;
255: DA da;
256: TstepCtx *tsCtx = user->tsCtx;
257: Parameter *param = user->param;
258: int i,j,mx,ierr,xs,ys,xm,ym;
259: int mglevel;
260: PetscReal grashof,dx;
261: Field **x;
263: da = (DA)(dmmg[param->mglevels-1]->dm);
264: grashof = user->param->grashof;
266: DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
267: dx = 1.0/(mx-1);
269: /*
270: Get local grid boundaries (for 2-dimensional DA):
271: xs, ys - starting grid indices (no ghost points)
272: xm, ym - widths of local grid (no ghost points)
273: */
274: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
276: /*
277: Get a pointer to vector data.
278: - For default PETSc vectors, VecGetArray() returns a pointer to
279: the data array. Otherwise, the routine is implementation dependent.
280: - You MUST call VecRestoreArray() when you no longer need access to
281: the array.
282: */
283: DAVecGetArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
285: /*
286: Compute initial guess over the locally owned part of the grid
287: Initial condition is motionless fluid and equilibrium temperature
288: */
289: for (j=ys; j<ys+ym; j++) {
290: for (i=xs; i<xs+xm; i++) {
291: x[j][i].u = 0.0;
292: x[j][i].v = 0.0;
293: x[j][i].omega = 0.0;
294: x[j][i].temp = (grashof>0)*i*dx;
295: }
296: }
298: /*
299: Restore vector
300: */
301: DAVecRestoreArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
302: /* Restrict Xold to coarser levels */
303: for (mglevel=param->mglevels-1; mglevel>0; mglevel--) {
304: MatRestrict(dmmg[mglevel]->R, ((AppCtx*)dmmg[mglevel]->user)->Xold, ((AppCtx*)dmmg[mglevel-1]->user)->Xold);
305: VecPointwiseMult(dmmg[mglevel]->Rscale,((AppCtx*)dmmg[mglevel-1]->user)->Xold,((AppCtx*)dmmg[mglevel-1]->user)->Xold);
306: }
307:
308: /* Store X in the qold for time stepping */
309: /*VecDuplicate(X,&tsCtx->qold);
310: VecDuplicate(X,&tsCtx->func);
311: VecCopy(X,tsCtx->Xold);
312: tsCtx->ires = 0;
313: SNESComputeFunction(snes,tsCtx->Xold,tsCtx->func);
314: VecNorm(tsCtx->func,NORM_2,&tsCtx->fnorm_ini);
315: tsCtx->ires = 1;
316: PetscPrintf(PETSC_COMM_WORLD,"Initial function norm is %g\n",tsCtx->fnorm_ini);*/
317: return 0;
318: }
319: /* ------------------------------------------------------------------- */
322: /*
323: FormInitialGuess - Forms initial approximation.
325: Input Parameters:
326: user - user-defined application context
327: X - vector
329: Output Parameter:
330: X - vector
331: */
332: int FormInitialGuess(SNES snes,Vec X,void *ptr)
333: {
334: DMMG dmmg = (DMMG)ptr;
335: AppCtx *user = (AppCtx*)dmmg->user;
336: TstepCtx *tsCtx = user->tsCtx;
337: int ierr;
338: VecCopy(user->Xold, X);
340: /* calculate the residual on fine mesh */
341: if (user->tsCtx->fnorm_ini == 0.0) {
342: tsCtx->ires = 0;
343: SNESComputeFunction(snes,user->Xold,user->func);
344: VecNorm(user->func,NORM_2,&tsCtx->fnorm_ini);
345: tsCtx->ires = 1;
346: tsCtx->cfl = tsCtx->cfl_ini;
347: }
349: return 0;
350: }
352: /*---------------------------------------------------------------------*/
355: int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,void *ptr)
356: /*---------------------------------------------------------------------*/
357: {
358: AppCtx *user = (AppCtx*)ptr;
359: TstepCtx *tsCtx = user->tsCtx;
360: DA da = info->da;
361: int ierr,i,j, xints,xinte,yints,yinte;
362: PetscReal hx,hy,dhx,dhy,hxhy;
363: PassiveScalar dtinv;
364: PassiveField **xold;
365: PetscTruth use_parab = tsCtx->use_parabolic;
368: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
369: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
370: hx = 1.0/dhx; hy = 1.0/dhy;
371: hxhy = hx*hy;
372: DAVecGetArray(da,user->Xold,&xold);
373: dtinv = hxhy/(tsCtx->cfl*tsCtx->dt);
374: /*
375: use_parab = PETSC_TRUE for parabolic equations; all the four equations have temporal term.
376: = PETSC_FALSE for differential algebraic equtions (DAE);
377: velocity equations do not have temporal term.
378: */
379: for (j=yints; j<yinte; j++) {
380: for (i=xints; i<xinte; i++) {
381: if (use_parab) {
382: f[j][i].u += dtinv*(x[j][i].u-xold[j][i].u);
383: f[j][i].v += dtinv*(x[j][i].v-xold[j][i].v);
384: }
385: f[j][i].omega += dtinv*(x[j][i].omega-xold[j][i].omega);
386: f[j][i].temp += dtinv*(x[j][i].temp-xold[j][i].temp);
387: }
388: }
389: DAVecRestoreArray(da,user->Xold,&xold);
390: return(0);
391: }
395: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
396: {
397: AppCtx *user = (AppCtx*)ptr;
398: TstepCtx *tsCtx = user->tsCtx;
399: int ierr,i,j;
400: int xints,xinte,yints,yinte;
401: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
402: PetscReal grashof,prandtl,lid;
403: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
406: grashof = user->param->grashof;
407: prandtl = user->param->prandtl;
408: lid = user->param->lidvelocity;
410: /*
411: Define mesh intervals ratios for uniform grid.
412: [Note: FD formulae below are normalized by multiplying through by
413: local volume element to obtain coefficients O(1) in two dimensions.]
414: */
415: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
416: hx = 1.0/dhx; hy = 1.0/dhy;
417: hxdhy = hx*dhy; hydhx = hy*dhx;
419: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
421: /* Test whether we are on the bottom edge of the global array */
422: if (yints == 0) {
423: j = 0;
424: yints = yints + 1;
425: /* bottom edge */
426: for (i=info->xs; i<info->xs+info->xm; i++) {
427: f[j][i].u = x[j][i].u;
428: f[j][i].v = x[j][i].v;
429: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
430: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
431: }
432: }
434: /* Test whether we are on the top edge of the global array */
435: if (yinte == info->my) {
436: j = info->my - 1;
437: yinte = yinte - 1;
438: /* top edge */
439: for (i=info->xs; i<info->xs+info->xm; i++) {
440: f[j][i].u = x[j][i].u - lid;
441: f[j][i].v = x[j][i].v;
442: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
443: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
444: }
445: }
447: /* Test whether we are on the left edge of the global array */
448: if (xints == 0) {
449: i = 0;
450: xints = xints + 1;
451: /* left edge */
452: for (j=info->ys; j<info->ys+info->ym; j++) {
453: f[j][i].u = x[j][i].u;
454: f[j][i].v = x[j][i].v;
455: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
456: f[j][i].temp = x[j][i].temp;
457: }
458: }
460: /* Test whether we are on the right edge of the global array */
461: if (xinte == info->mx) {
462: i = info->mx - 1;
463: xinte = xinte - 1;
464: /* right edge */
465: for (j=info->ys; j<info->ys+info->ym; j++) {
466: f[j][i].u = x[j][i].u;
467: f[j][i].v = x[j][i].v;
468: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
469: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
470: }
471: }
473: /* Compute over the interior points */
474: for (j=yints; j<yinte; j++) {
475: for (i=xints; i<xinte; i++) {
477: /*
478: convective coefficients for upwinding
479: */
480: vx = x[j][i].u; avx = PetscAbsScalar(vx);
481: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
482: vy = x[j][i].v; avy = PetscAbsScalar(vy);
483: vyp = .5*(vy+avy); vym = .5*(vy-avy);
485: /* U velocity */
486: u = x[j][i].u;
487: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
488: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
489: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
491: /* V velocity */
492: u = x[j][i].v;
493: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
494: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
495: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
497: /* Omega */
498: u = x[j][i].omega;
499: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
500: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
501: f[j][i].omega = uxx + uyy +
502: (vxp*(u - x[j][i-1].omega) +
503: vxm*(x[j][i+1].omega - u)) * hy +
504: (vyp*(u - x[j-1][i].omega) +
505: vym*(x[j+1][i].omega - u)) * hx -
506: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
508: /* Temperature */
509: u = x[j][i].temp;
510: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
511: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
512: f[j][i].temp = uxx + uyy + prandtl * (
513: (vxp*(u - x[j][i-1].temp) +
514: vxm*(x[j][i+1].temp - u)) * hy +
515: (vyp*(u - x[j-1][i].temp) +
516: vym*(x[j+1][i].temp - u)) * hx);
517: }
518: }
520: /* Add time step contribution */
521: if (tsCtx->ires) {
522: AddTSTermLocal(info,x,f,ptr);
523: }
524: /*
525: Flop count (multiply-adds are counted as 2 operations)
526: */
527: PetscLogFlops(84*info->ym*info->xm);
528: return(0);
529: }
533: int FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
534: {
535: AppCtx *user = (AppCtx*)ptr;
536: int i,j,c;
537: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
538: PassiveReal grashof,prandtl,lid;
539: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
542: grashof = user->param->grashof;
543: prandtl = user->param->prandtl;
544: lid = user->param->lidvelocity;
546: /*
547: Define mesh intervals ratios for uniform grid.
548: [Note: FD formulae below are normalized by multiplying through by
549: local volume element to obtain coefficients O(1) in two dimensions.]
550: */
551: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
552: hx = 1.0/dhx; hy = 1.0/dhy;
553: hxdhy = hx*dhy; hydhx = hy*dhx;
555: i = st->i; j = st->j; c = st->c;
557: /* Test whether we are on the right edge of the global array */
558: if (i == info->mx-1) {
559: if (c == 0) *f = x[j][i].u;
560: else if (c == 1) *f = x[j][i].v;
561: else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
562: else *f = x[j][i].temp - (PetscReal)(grashof>0);
564: /* Test whether we are on the left edge of the global array */
565: } else if (i == 0) {
566: if (c == 0) *f = x[j][i].u;
567: else if (c == 1) *f = x[j][i].v;
568: else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
569: else *f = x[j][i].temp;
571: /* Test whether we are on the top edge of the global array */
572: } else if (j == info->my-1) {
573: if (c == 0) *f = x[j][i].u - lid;
574: else if (c == 1) *f = x[j][i].v;
575: else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
576: else *f = x[j][i].temp-x[j-1][i].temp;
578: /* Test whether we are on the bottom edge of the global array */
579: } else if (j == 0) {
580: if (c == 0) *f = x[j][i].u;
581: else if (c == 1) *f = x[j][i].v;
582: else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
583: else *f = x[j][i].temp-x[j+1][i].temp;
585: /* Compute over the interior points */
586: } else {
587: /*
588: convective coefficients for upwinding
589: */
590: vx = x[j][i].u; avx = PetscAbsScalar(vx);
591: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
592: vy = x[j][i].v; avy = PetscAbsScalar(vy);
593: vyp = .5*(vy+avy); vym = .5*(vy-avy);
595: /* U velocity */
596: if (c == 0) {
597: u = x[j][i].u;
598: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
599: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
600: *f = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
602: /* V velocity */
603: } else if (c == 1) {
604: u = x[j][i].v;
605: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
606: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
607: *f = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
608:
609: /* Omega */
610: } else if (c == 2) {
611: u = x[j][i].omega;
612: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
613: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
614: *f = uxx + uyy +
615: (vxp*(u - x[j][i-1].omega) +
616: vxm*(x[j][i+1].omega - u)) * hy +
617: (vyp*(u - x[j-1][i].omega) +
618: vym*(x[j+1][i].omega - u)) * hx -
619: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
620:
621: /* Temperature */
622: } else {
623: u = x[j][i].temp;
624: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
625: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
626: *f = uxx + uyy + prandtl * (
627: (vxp*(u - x[j][i-1].temp) +
628: vxm*(x[j][i+1].temp - u)) * hy +
629: (vyp*(u - x[j-1][i].temp) +
630: vym*(x[j+1][i].temp - u)) * hx);
631: }
632: }
634: return(0);
635: }
638: /*---------------------------------------------------------------------*/
641: int Update(DMMG *dmmg)
642: /*---------------------------------------------------------------------*/
643: {
644:
645: AppCtx *user = (AppCtx *) ((dmmg[0])->user);
646: TstepCtx *tsCtx = user->tsCtx;
647: Parameter *param = user->param;
648: SNES snes;
649: int ierr,its;
650: PetscScalar fratio;
651: PetscScalar time1,time2,cpuloc;
652: int max_steps;
653: PetscTruth print_flag = PETSC_FALSE;
654: FILE *fptr = 0;
655: int nfailsCum = 0,nfails = 0;
659: PetscOptionsHasName(PETSC_NULL,"-print",&print_flag);
660: if (user->param->PreLoading)
661: max_steps = 1;
662: else
663: max_steps = tsCtx->max_steps;
664: fratio = 1.0;
665:
666: PetscGetTime(&time1);
667: for (tsCtx->itstep = 0; (tsCtx->itstep < max_steps) &&
668: (fratio <= tsCtx->fnorm_ratio); tsCtx->itstep++) {
669: DMMGSolve(dmmg);
670: snes = DMMGGetSNES(dmmg);
671: SNESGetIterationNumber(snes,&its);
672: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n", its);
673: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
674: nfailsCum += nfails; nfails = 0;
675: if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
676: /*tsCtx->qcur = DMMGGetx(dmmg);
677: VecCopy(tsCtx->qcur,tsCtx->qold);*/
679: VecCopy(dmmg[param->mglevels-1]->x, ((AppCtx*)dmmg[param->mglevels-1]->user)->Xold);
680: for (its=param->mglevels-1; its>0 ;its--) {
681: MatRestrict(dmmg[its]->R, ((AppCtx*)dmmg[its]->user)->Xold, ((AppCtx*)dmmg[its-1]->user)->Xold);
682: VecPointwiseMult(dmmg[its]->Rscale,((AppCtx*)dmmg[its-1]->user)->Xold,((AppCtx*)dmmg[its-1]->user)->Xold);
683: }
685:
686: ComputeTimeStep(snes,((AppCtx*)dmmg[param->mglevels-1]->user));
687: if (print_flag) {
688: PetscPrintf(PETSC_COMM_WORLD,"At Time Step %d cfl = %g and fnorm = %g\n",
689: tsCtx->itstep,tsCtx->cfl,tsCtx->fnorm);
690: PetscPrintf(PETSC_COMM_WORLD,"Wall clock time needed %g seconds for %d time steps\n",
691: cpuloc,tsCtx->itstep);
692: }
693: fratio = tsCtx->fnorm_ini/tsCtx->fnorm;
694: PetscGetTime(&time2);
695: cpuloc = time2-time1;
696: MPI_Barrier(PETSC_COMM_WORLD);
697: } /* End of time step loop */
698:
699: PetscPrintf(PETSC_COMM_WORLD,"Total wall clock time needed %g seconds for %d time steps\n",
700: cpuloc,tsCtx->itstep);
701: PetscPrintf(PETSC_COMM_WORLD,"cfl = %g fnorm = %g\n",tsCtx->cfl,tsCtx->fnorm);
702: if (user->param->PreLoading) {
703: tsCtx->fnorm_ini = 0.0;
704: PetscPrintf(PETSC_COMM_WORLD,"Preloading done ...\n");
705: }
706: /*
707: {
708: Vec xx,yy;
709: PetscScalar fnorm,fnorm1;
710: SNESGetFunctionNorm(snes,&fnorm);
711: xx = DMMGGetx(dmmg);
712: VecDuplicate(xx,&yy);
713: SNESComputeFunction(snes,xx,yy);
714: VecNorm(yy,NORM_2,&fnorm1);
715: PetscPrintf(PETSC_COMM_WORLD,"fnorm = %g, fnorm1 = %g\n",fnorm,fnorm1);
716:
717: }
718: */
720: return(0);
721: }
723: /*---------------------------------------------------------------------*/
726: int ComputeTimeStep(SNES snes,void *ptr)
727: /*---------------------------------------------------------------------*/
728: {
729: AppCtx *user = (AppCtx*)ptr;
730: TstepCtx *tsCtx = user->tsCtx;
731: Vec func = user->func;
732: PetscScalar inc = 1.1, newcfl;
733: int ierr;
734: /*int iramp = tsCtx->iramp;*/
735:
737: tsCtx->dt = 0.01;
738: PetscOptionsGetScalar(PETSC_NULL,"-deltat",&tsCtx->dt,PETSC_NULL);
739: tsCtx->ires = 0;
740: SNESComputeFunction(snes,user->Xold,user->func);
741: tsCtx->ires = 1;
742: VecNorm(func,NORM_2,&tsCtx->fnorm);
743: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
744: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
745: /* first time through so compute initial function norm */
746: /*if (tsCtx->fnorm_ini == 0.0) {
747: tsCtx->fnorm_ini = tsCtx->fnorm;
748: tsCtx->cfl = tsCtx->cfl_ini;
749: } else {
750: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
751: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
752: }*/
753: return(0);
754: }