Actual source code: ex27.c
2: static char help[] = "Nonlinear driven cavity with multigrid and pseudo timestepping 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: dU/dt - Lap(U) - Grad_y(Omega) = 0
26: dV/dt - Lap(V) + Grad_x(Omega) = 0
27: d(omega)/dt - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
28: dT/dt - 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).
53: This example is used in the paper
55: Todd S. Coffey and C. T. Kelley and David E. Keyes,
56: Pseudotransient Continuation and Differential-Algebraic Equations, 2003.
58: however Figure 3.1 was generated with a slightly different algorithm (see
59: targets runex27 and runex27_p) than described in the paper. In particular,
60: the described algorithm is linearly implicit, advancing to the next step
61: after one Newton step, so that the steady state residual is always used, but
62: the figure was generated by converging each step to a relative tolerance of
63: 1.e-3. On the example problem, setting -snes_max_it 1 has only minor impact
64: on number of steps, but significantly reduces the required number of linear
65: solves.
67: ------------------------------------------------------------------------- */
69: /*
70: Include "petscda.h" so that we can use distributed arrays (DAs).
71: Include "petscsnes.h" so that we can use SNES solvers. Note that this
72: file automatically includes:
73: petscsys.h - base PETSc routines petscvec.h - vectors
74: petscmat.h - matrices
75: petscis.h - index sets petscksp.h - Krylov subspace methods
76: petscviewer.h - viewers petscpc.h - preconditioners
77: petscksp.h - linear solvers
78: */
79: #include petscsnes.h
80: #include petscda.h
81: #include petscdmmg.h
83: /*
84: User-defined routines and data structures
85: */
87: typedef struct {
88: PassiveScalar fnorm_ini,dt_ini,cfl_ini;
89: PassiveScalar fnorm,dt,cfl;
90: PassiveScalar ptime;
91: PassiveScalar cfl_max,max_time;
92: PassiveScalar fnorm_ratio;
93: PetscInt ires,iramp,itstep;
94: PetscInt max_steps,print_freq;
95: PetscInt LocalTimeStepping;
96: PetscTruth use_parabolic;
97: } TstepCtx;
99: /*
100: The next two structures are essentially the same. The difference is that
101: the first does not contain derivative information (as used by ADIC) while the
102: second does. The first is used only to contain the solution at the previous time-step
103: which is a constant in the computation of the current time step and hence passive
104: (meaning not active in terms of derivatives).
105: */
106: typedef struct {
107: PassiveScalar u,v,omega,temp;
108: } PassiveField;
110: typedef struct {
111: PetscScalar u,v,omega,temp;
112: } Field;
115: typedef struct {
116: PetscInt mglevels;
117: PetscInt cycles; /* numbers of time steps for integration */
118: PassiveReal lidvelocity,prandtl,grashof; /* physical parameters */
119: PetscTruth draw_contours; /* indicates drawing contours of solution */
120: PetscTruth PreLoading;
121: } Parameter;
123: typedef struct {
124: Vec Xold,func;
125: TstepCtx *tsCtx;
126: Parameter *param;
127: } AppCtx;
140: int main(int argc,char **argv)
141: {
142: DMMG *dmmg; /* multilevel grid structure */
143: AppCtx *user; /* user-defined work context */
144: TstepCtx tsCtx;
145: Parameter param;
146: PetscInt mx,my,i;
148: MPI_Comm comm;
149: DA da;
151: PetscInitialize(&argc,&argv,(char *)0,help);
152: comm = PETSC_COMM_WORLD;
155: PreLoadBegin(PETSC_TRUE,"SetUp");
157: param.PreLoading = PreLoading;
158: DMMGCreate(comm,2,&user,&dmmg);
159: param.mglevels = DMMGGetLevels(dmmg);
162: /*
163: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
164: for principal unknowns (x) and governing residuals (f)
165: */
166: DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
167: DMMGSetDM(dmmg,(DM)da);
168: DADestroy(da);
170: DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
171: PETSC_IGNORE,PETSC_IGNORE);
172: /*
173: Problem parameters (velocity of lid, prandtl, and grashof numbers)
174: */
175: param.lidvelocity = 1.0/(mx*my);
176: param.prandtl = 1.0;
177: param.grashof = 1.0;
178: PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",¶m.lidvelocity,PETSC_NULL);
179: PetscOptionsGetReal(PETSC_NULL,"-prandtl",¶m.prandtl,PETSC_NULL);
180: PetscOptionsGetReal(PETSC_NULL,"-grashof",¶m.grashof,PETSC_NULL);
181: PetscOptionsHasName(PETSC_NULL,"-contours",¶m.draw_contours);
183: DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
184: DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
185: DASetFieldName(DMMGGetDA(dmmg),2,"Omega");
186: DASetFieldName(DMMGGetDA(dmmg),3,"temperature");
188: /*======================================================================*/
189: /* Initilize stuff related to time stepping */
190: /*======================================================================*/
191: tsCtx.fnorm_ini = 0.0; tsCtx.cfl_ini = 50.0; tsCtx.cfl_max = 1.0e+06;
192: tsCtx.max_steps = 50; tsCtx.max_time = 1.0e+12; tsCtx.iramp = -50;
193: tsCtx.dt = 0.01; tsCtx.fnorm_ratio = 1.0e+10;
194: tsCtx.LocalTimeStepping = 0;
195: tsCtx.use_parabolic = PETSC_FALSE;
196: PetscOptionsGetTruth(PETSC_NULL,"-use_parabolic",&tsCtx.use_parabolic,PETSC_IGNORE);
197: PetscOptionsGetInt(PETSC_NULL,"-max_st",&tsCtx.max_steps,PETSC_NULL);
198: PetscOptionsGetReal(PETSC_NULL,"-ts_rtol",&tsCtx.fnorm_ratio,PETSC_NULL);
199: PetscOptionsGetReal(PETSC_NULL,"-cfl_ini",&tsCtx.cfl_ini,PETSC_NULL);
200: PetscOptionsGetReal(PETSC_NULL,"-cfl_max",&tsCtx.cfl_max,PETSC_NULL);
201: tsCtx.print_freq = tsCtx.max_steps;
202: PetscOptionsGetInt(PETSC_NULL,"-print_freq",&tsCtx.print_freq,PETSC_NULL);
203:
204: PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
205: for (i=0; i<param.mglevels; i++) {
206: VecDuplicate(dmmg[i]->x, &(user[i].Xold));
207: VecDuplicate(dmmg[i]->x, &(user[i].func));
208: user[i].tsCtx = &tsCtx;
209: user[i].param = ¶m;
210: dmmg[i]->user = &user[i];
211: }
212: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213: Create user context, set problem data, create vector data structures.
214: Also, compute the initial guess.
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216:
217: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218: Create nonlinear solver context
219:
220: Process adiC(20): AddTSTermLocal FormFunctionLocal FormFunctionLocali
221: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
223: DMMGSetFromOptions(dmmg);
224: DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
225:
226: PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, grashof # = %G\n",
227: param.lidvelocity,param.prandtl,param.grashof);
228:
229:
230: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: Solve the nonlinear system
232: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233: DMMGSetInitialGuess(dmmg,FormInitialGuess);
234:
235: PreLoadStage("Solve");
236: Initialize(dmmg);
237: Update(dmmg);
238: /*
239: Visualize solution
240: */
241:
242: if (param.draw_contours) {
243: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
244: }
245: /*VecView(DMMGGetx(dmmg),PETSC_VIEWER_STDOUT_WORLD);*/
246: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247: Free work space. All PETSc objects should be destroyed when they
248: are no longer needed.
249: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
250:
251: for (i=0; i<param.mglevels; i++) {
252: VecDestroy(user[i].Xold);
253: VecDestroy(user[i].func);
254: }
255: PetscFree(user);
256: DMMGDestroy(dmmg);
257: PreLoadEnd();
258:
259: PetscFinalize();
260: return 0;
261: }
263: /* ------------------------------------------------------------------- */
266: /* ------------------------------------------------------------------- */
267: PetscErrorCode Initialize(DMMG *dmmg)
268: {
269: AppCtx *user = (AppCtx*)dmmg[0]->user;
270: DA da;
271: Parameter *param = user->param;
272: PetscInt i,j,mx,xs,ys,xm,ym,mglevel;
274: PetscReal grashof,dx;
275: Field **x;
277: da = (DA)(dmmg[param->mglevels-1]->dm);
278: grashof = user->param->grashof;
280: DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
281: dx = 1.0/(mx-1);
283: /*
284: Get local grid boundaries (for 2-dimensional DA):
285: xs, ys - starting grid indices (no ghost points)
286: xm, ym - widths of local grid (no ghost points)
287: */
288: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
290: /*
291: Get a pointer to vector data.
292: - For default PETSc vectors, VecGetArray() returns a pointer to
293: the data array. Otherwise, the routine is implementation dependent.
294: - You MUST call VecRestoreArray() when you no longer need access to
295: the array.
296: */
297: DAVecGetArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
299: /*
300: Compute initial guess over the locally owned part of the grid
301: Initial condition is motionless fluid and equilibrium temperature
302: */
303: for (j=ys; j<ys+ym; j++) {
304: for (i=xs; i<xs+xm; i++) {
305: x[j][i].u = 0.0;
306: x[j][i].v = 0.0;
307: x[j][i].omega = 0.0;
308: x[j][i].temp = (grashof>0)*i*dx;
309: }
310: }
312: /*
313: Restore vector
314: */
315: DAVecRestoreArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
316: /* Restrict Xold to coarser levels */
317: for (mglevel=param->mglevels-1; mglevel>0; mglevel--) {
318: MatRestrict(dmmg[mglevel]->R, ((AppCtx*)dmmg[mglevel]->user)->Xold, ((AppCtx*)dmmg[mglevel-1]->user)->Xold);
319: VecPointwiseMult(((AppCtx*)dmmg[mglevel-1]->user)->Xold,((AppCtx*)dmmg[mglevel-1]->user)->Xold,dmmg[mglevel]->Rscale);
320: }
321:
322: /* Store X in the qold for time stepping */
323: /*VecDuplicate(X,&tsCtx->qold);
324: VecDuplicate(X,&tsCtx->func);
325: VecCopy(X,tsCtx->Xold);
326: tsCtx->ires = 0;
327: SNESComputeFunction(snes,tsCtx->Xold,tsCtx->func);
328: VecNorm(tsCtx->func,NORM_2,&tsCtx->fnorm_ini);
329: tsCtx->ires = 1;
330: PetscPrintf(PETSC_COMM_WORLD,"Initial function norm is %G\n",tsCtx->fnorm_ini);*/
331: return 0;
332: }
333: /* ------------------------------------------------------------------- */
336: /*
337: FormInitialGuess - Forms initial approximation.
339: Input Parameters:
340: user - user-defined application context
341: X - vector
343: Output Parameter:
344: X - vector
345: */
346: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
347: {
348: AppCtx *user = (AppCtx*)dmmg->user;
349: TstepCtx *tsCtx = user->tsCtx;
352: VecCopy(user->Xold, X);
353: /* calculate the residual on fine mesh */
354: if (user->tsCtx->fnorm_ini == 0.0) {
355: tsCtx->ires = 0;
356: SNESComputeFunction(dmmg->snes,user->Xold,user->func);
357: VecNorm(user->func,NORM_2,&tsCtx->fnorm_ini);
358: tsCtx->ires = 1;
359: tsCtx->cfl = tsCtx->cfl_ini;
360: }
361: return 0;
362: }
364: /*---------------------------------------------------------------------*/
367: PetscErrorCode AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,void *ptr)
368: /*---------------------------------------------------------------------*/
369: {
370: AppCtx *user = (AppCtx*)ptr;
371: TstepCtx *tsCtx = user->tsCtx;
372: DA da = info->da;
374: PetscInt i,j, xints,xinte,yints,yinte;
375: PetscReal hx,hy,dhx,dhy,hxhy;
376: PassiveScalar dtinv;
377: PassiveField **xold;
378: PetscTruth use_parab = tsCtx->use_parabolic;
381: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
382: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
383: hx = 1.0/dhx; hy = 1.0/dhy;
384: hxhy = hx*hy;
385: DAVecGetArray(da,user->Xold,&xold);
386: dtinv = hxhy/(tsCtx->cfl*tsCtx->dt);
387: /*
388: use_parab = PETSC_TRUE for parabolic equations; all the four equations have temporal term.
389: = PETSC_FALSE for differential algebraic equtions (DAE);
390: velocity equations do not have temporal term.
391: */
392: for (j=yints; j<yinte; j++) {
393: for (i=xints; i<xinte; i++) {
394: if (use_parab) {
395: f[j][i].u += dtinv*(x[j][i].u-xold[j][i].u);
396: f[j][i].v += dtinv*(x[j][i].v-xold[j][i].v);
397: }
398: f[j][i].omega += dtinv*(x[j][i].omega-xold[j][i].omega);
399: f[j][i].temp += dtinv*(x[j][i].temp-xold[j][i].temp);
400: }
401: }
402: DAVecRestoreArray(da,user->Xold,&xold);
403: return(0);
404: }
408: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
409: {
410: AppCtx *user = (AppCtx*)ptr;
411: TstepCtx *tsCtx = user->tsCtx;
413: PetscInt xints,xinte,yints,yinte,i,j;
414: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
415: PetscReal grashof,prandtl,lid;
416: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
419: grashof = user->param->grashof;
420: prandtl = user->param->prandtl;
421: lid = user->param->lidvelocity;
423: /*
424: Define mesh intervals ratios for uniform grid.
425: [Note: FD formulae below are normalized by multiplying through by
426: local volume element to obtain coefficients O(1) in two dimensions.]
427: */
428: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
429: hx = 1.0/dhx; hy = 1.0/dhy;
430: hxdhy = hx*dhy; hydhx = hy*dhx;
432: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
434: /* Test whether we are on the bottom edge of the global array */
435: if (yints == 0) {
436: j = 0;
437: yints = yints + 1;
438: /* bottom edge */
439: for (i=info->xs; i<info->xs+info->xm; i++) {
440: f[j][i].u = x[j][i].u;
441: f[j][i].v = x[j][i].v;
442: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][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 top edge of the global array */
448: if (yinte == info->my) {
449: j = info->my - 1;
450: yinte = yinte - 1;
451: /* top edge */
452: for (i=info->xs; i<info->xs+info->xm; i++) {
453: f[j][i].u = x[j][i].u - lid;
454: f[j][i].v = x[j][i].v;
455: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
456: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
457: }
458: }
460: /* Test whether we are on the left edge of the global array */
461: if (xints == 0) {
462: i = 0;
463: xints = xints + 1;
464: /* left 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+1].v - x[j][i].v)*dhx;
469: f[j][i].temp = x[j][i].temp;
470: }
471: }
473: /* Test whether we are on the right edge of the global array */
474: if (xinte == info->mx) {
475: i = info->mx - 1;
476: xinte = xinte - 1;
477: /* right edge */
478: for (j=info->ys; j<info->ys+info->ym; j++) {
479: f[j][i].u = x[j][i].u;
480: f[j][i].v = x[j][i].v;
481: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
482: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
483: }
484: }
486: /* Compute over the interior points */
487: for (j=yints; j<yinte; j++) {
488: for (i=xints; i<xinte; i++) {
490: /*
491: convective coefficients for upwinding
492: */
493: vx = x[j][i].u; avx = PetscAbsScalar(vx);
494: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
495: vy = x[j][i].v; avy = PetscAbsScalar(vy);
496: vyp = .5*(vy+avy); vym = .5*(vy-avy);
498: /* U velocity */
499: u = x[j][i].u;
500: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
501: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
502: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
504: /* V velocity */
505: u = x[j][i].v;
506: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
507: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
508: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
510: /* Omega */
511: u = x[j][i].omega;
512: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
513: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
514: f[j][i].omega = uxx + uyy +
515: (vxp*(u - x[j][i-1].omega) +
516: vxm*(x[j][i+1].omega - u)) * hy +
517: (vyp*(u - x[j-1][i].omega) +
518: vym*(x[j+1][i].omega - u)) * hx -
519: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
521: /* Temperature */
522: u = x[j][i].temp;
523: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
524: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
525: f[j][i].temp = uxx + uyy + prandtl * (
526: (vxp*(u - x[j][i-1].temp) +
527: vxm*(x[j][i+1].temp - u)) * hy +
528: (vyp*(u - x[j-1][i].temp) +
529: vym*(x[j+1][i].temp - u)) * hx);
530: }
531: }
533: /* Add time step contribution */
534: if (tsCtx->ires) {
535: AddTSTermLocal(info,x,f,ptr);
536: }
537: /*
538: Flop count (multiply-adds are counted as 2 operations)
539: */
540: PetscLogFlops(84.0*info->ym*info->xm);
541: return(0);
542: }
546: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
547: {
548: AppCtx *user = (AppCtx*)ptr;
549: PetscInt i,j,c;
550: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
551: PassiveReal grashof,prandtl,lid;
552: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
555: grashof = user->param->grashof;
556: prandtl = user->param->prandtl;
557: lid = user->param->lidvelocity;
559: /*
560: Define mesh intervals ratios for uniform grid.
561: [Note: FD formulae below are normalized by multiplying through by
562: local volume element to obtain coefficients O(1) in two dimensions.]
563: */
564: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
565: hx = 1.0/dhx; hy = 1.0/dhy;
566: hxdhy = hx*dhy; hydhx = hy*dhx;
568: i = st->i; j = st->j; c = st->c;
570: /* Test whether we are on the right edge of the global array */
571: if (i == info->mx-1) {
572: if (c == 0) *f = x[j][i].u;
573: else if (c == 1) *f = x[j][i].v;
574: else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
575: else *f = x[j][i].temp - (PetscReal)(grashof>0);
577: /* Test whether we are on the left edge of the global array */
578: } else if (i == 0) {
579: if (c == 0) *f = x[j][i].u;
580: else if (c == 1) *f = x[j][i].v;
581: else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
582: else *f = x[j][i].temp;
584: /* Test whether we are on the top edge of the global array */
585: } else if (j == info->my-1) {
586: if (c == 0) *f = x[j][i].u - lid;
587: else if (c == 1) *f = x[j][i].v;
588: else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
589: else *f = x[j][i].temp-x[j-1][i].temp;
591: /* Test whether we are on the bottom edge of the global array */
592: } else if (j == 0) {
593: if (c == 0) *f = x[j][i].u;
594: else if (c == 1) *f = x[j][i].v;
595: else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
596: else *f = x[j][i].temp-x[j+1][i].temp;
598: /* Compute over the interior points */
599: } else {
600: /*
601: convective coefficients for upwinding
602: */
603: vx = x[j][i].u; avx = PetscAbsScalar(vx);
604: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
605: vy = x[j][i].v; avy = PetscAbsScalar(vy);
606: vyp = .5*(vy+avy); vym = .5*(vy-avy);
608: /* U velocity */
609: if (c == 0) {
610: u = x[j][i].u;
611: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
612: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
613: *f = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
615: /* V velocity */
616: } else if (c == 1) {
617: u = x[j][i].v;
618: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
619: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
620: *f = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
621:
622: /* Omega */
623: } else if (c == 2) {
624: u = x[j][i].omega;
625: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
626: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
627: *f = uxx + uyy +
628: (vxp*(u - x[j][i-1].omega) +
629: vxm*(x[j][i+1].omega - u)) * hy +
630: (vyp*(u - x[j-1][i].omega) +
631: vym*(x[j+1][i].omega - u)) * hx -
632: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
633:
634: /* Temperature */
635: } else {
636: u = x[j][i].temp;
637: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
638: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
639: *f = uxx + uyy + prandtl * (
640: (vxp*(u - x[j][i-1].temp) +
641: vxm*(x[j][i+1].temp - u)) * hy +
642: (vyp*(u - x[j-1][i].temp) +
643: vym*(x[j+1][i].temp - u)) * hx);
644: }
645: }
647: return(0);
648: }
651: /*---------------------------------------------------------------------*/
654: PetscErrorCode Update(DMMG *dmmg)
655: /*---------------------------------------------------------------------*/
656: {
657: AppCtx *user = (AppCtx *) ((dmmg[0])->user);
658: TstepCtx *tsCtx = user->tsCtx;
659: Parameter *param = user->param;
660: SNES snes;
662: PetscScalar fratio;
663: PetscScalar time1,time2,cpuloc;
664: PetscInt max_steps,its;
665: PetscTruth print_flag = PETSC_FALSE;
666: PetscInt nfailsCum = 0,nfails = 0;
669: /* Note: print_flag displays diagnostic info, not convergence behavior. Use '-snes_monitor' for converges info. */
670: PetscOptionsHasName(PETSC_NULL,"-print",&print_flag);
671: if (user->param->PreLoading)
672: max_steps = 1;
673: else
674: max_steps = tsCtx->max_steps;
675: fratio = 1.0;
676:
677: PetscGetTime(&time1);
678: for (tsCtx->itstep = 0; (tsCtx->itstep < max_steps) &&
679: (fratio <= tsCtx->fnorm_ratio); tsCtx->itstep++) {
680: DMMGSolve(dmmg);
681: snes = DMMGGetSNES(dmmg);
682: SNESGetIterationNumber(snes,&its);
683: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n", its);
684: SNESGetNonlinearStepFailures(snes,&nfails);
685: nfailsCum += nfails; nfails = 0;
686: if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
687: /*tsCtx->qcur = DMMGGetx(dmmg);
688: VecCopy(tsCtx->qcur,tsCtx->qold);*/
689:
690: VecCopy(dmmg[param->mglevels-1]->x, ((AppCtx*)dmmg[param->mglevels-1]->user)->Xold);
691: for (its=param->mglevels-1; its>0 ;its--) {
692: MatRestrict(dmmg[its]->R, ((AppCtx*)dmmg[its]->user)->Xold, ((AppCtx*)dmmg[its-1]->user)->Xold);
693: VecPointwiseMult(((AppCtx*)dmmg[its-1]->user)->Xold,((AppCtx*)dmmg[its-1]->user)->Xold,dmmg[its]->Rscale);
694: }
696:
697: ComputeTimeStep(snes,((AppCtx*)dmmg[param->mglevels-1]->user));
698: if (print_flag) {
699: PetscPrintf(PETSC_COMM_WORLD,"At Time Step %D cfl = %G and fnorm = %G\n",
700: tsCtx->itstep,tsCtx->cfl,tsCtx->fnorm);
701: PetscPrintf(PETSC_COMM_WORLD,"Wall clock time needed %G seconds for %D time steps\n",
702: cpuloc,tsCtx->itstep);
703: }
704: fratio = tsCtx->fnorm_ini/tsCtx->fnorm;
705: PetscGetTime(&time2);
706: cpuloc = time2-time1;
707: MPI_Barrier(PETSC_COMM_WORLD);
708: } /* End of time step loop */
710: if (print_flag) {
711: PetscPrintf(PETSC_COMM_WORLD,"Total wall clock time needed %G seconds for %D time steps\n",
712: cpuloc,tsCtx->itstep);
713: PetscPrintf(PETSC_COMM_WORLD,"cfl = %G fnorm = %G\n",tsCtx->cfl,tsCtx->fnorm);
714: }
715: if (user->param->PreLoading) {
716: tsCtx->fnorm_ini = 0.0;
717: PetscPrintf(PETSC_COMM_WORLD,"Preloading done ...\n");
718: }
719: /*
720: {
721: Vec xx,yy;
722: PetscScalar fnorm,fnorm1;
723: SNESGetFunctionNorm(snes,&fnorm);
724: xx = DMMGGetx(dmmg);
725: VecDuplicate(xx,&yy);
726: SNESComputeFunction(snes,xx,yy);
727: VecNorm(yy,NORM_2,&fnorm1);
728: PetscPrintf(PETSC_COMM_WORLD,"fnorm = %G, fnorm1 = %G\n",fnorm,fnorm1);
729:
730: }
731: */
733: return(0);
734: }
736: /*---------------------------------------------------------------------*/
739: PetscErrorCode ComputeTimeStep(SNES snes,void *ptr)
740: /*---------------------------------------------------------------------*/
741: {
742: AppCtx *user = (AppCtx*)ptr;
743: TstepCtx *tsCtx = user->tsCtx;
744: Vec func = user->func;
745: PetscScalar inc = 1.1, newcfl;
747: /*int iramp = tsCtx->iramp;*/
748:
750: tsCtx->dt = 0.01;
751: PetscOptionsGetScalar(PETSC_NULL,"-deltat",&tsCtx->dt,PETSC_NULL);
752: tsCtx->ires = 0;
753: SNESComputeFunction(snes,user->Xold,user->func);
754: tsCtx->ires = 1;
755: VecNorm(func,NORM_2,&tsCtx->fnorm);
756: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
757: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
758: /* first time through so compute initial function norm */
759: /*if (tsCtx->fnorm_ini == 0.0) {
760: tsCtx->fnorm_ini = tsCtx->fnorm;
761: tsCtx->cfl = tsCtx->cfl_ini;
762: } else {
763: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
764: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
765: }*/
766: return(0);
767: }