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).
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: */
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: PetscInt ires,iramp,itstep;
79: PetscInt max_steps,print_freq;
80: PetscInt 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: PetscInt mglevels;
102: PetscInt 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;
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: PetscInt mx,my,i;
133: MPI_Comm comm;
134: DA da;
136: PetscInitialize(&argc,&argv,(char *)0,help);
137: comm = PETSC_COMM_WORLD;
140: PreLoadBegin(PETSC_TRUE,"SetUp");
142: param.PreLoading = PreLoading;
143: DMMGCreate(comm,2,&user,&dmmg);
144: param.mglevels = DMMGGetLevels(dmmg);
147: /*
148: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
149: for principal unknowns (x) and governing residuals (f)
150: */
151: DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
152: DMMGSetDM(dmmg,(DM)da);
153: DADestroy(da);
155: DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
156: PETSC_IGNORE,PETSC_IGNORE);
157: /*
158: Problem parameters (velocity of lid, prandtl, and grashof numbers)
159: */
160: param.lidvelocity = 1.0/(mx*my);
161: param.prandtl = 1.0;
162: param.grashof = 1.0;
163: PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",¶m.lidvelocity,PETSC_NULL);
164: PetscOptionsGetReal(PETSC_NULL,"-prandtl",¶m.prandtl,PETSC_NULL);
165: PetscOptionsGetReal(PETSC_NULL,"-grashof",¶m.grashof,PETSC_NULL);
166: PetscOptionsHasName(PETSC_NULL,"-contours",¶m.draw_contours);
168: DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
169: DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
170: DASetFieldName(DMMGGetDA(dmmg),2,"Omega");
171: DASetFieldName(DMMGGetDA(dmmg),3,"temperature");
173: /*======================================================================*/
174: /* Initilize stuff related to time stepping */
175: /*======================================================================*/
176: tsCtx.fnorm_ini = 0.0; tsCtx.cfl_ini = 50.0; tsCtx.cfl_max = 1.0e+06;
177: tsCtx.max_steps = 50; tsCtx.max_time = 1.0e+12; tsCtx.iramp = -50;
178: tsCtx.dt = 0.01; tsCtx.fnorm_ratio = 1.0e+10;
179: tsCtx.LocalTimeStepping = 0;
180: tsCtx.use_parabolic = PETSC_FALSE;
181: PetscOptionsGetTruth(PETSC_NULL,"-use_parabolic",&tsCtx.use_parabolic,PETSC_IGNORE);
182: PetscOptionsGetInt(PETSC_NULL,"-max_st",&tsCtx.max_steps,PETSC_NULL);
183: PetscOptionsGetReal(PETSC_NULL,"-ts_rtol",&tsCtx.fnorm_ratio,PETSC_NULL);
184: PetscOptionsGetReal(PETSC_NULL,"-cfl_ini",&tsCtx.cfl_ini,PETSC_NULL);
185: PetscOptionsGetReal(PETSC_NULL,"-cfl_max",&tsCtx.cfl_max,PETSC_NULL);
186: tsCtx.print_freq = tsCtx.max_steps;
187: PetscOptionsGetInt(PETSC_NULL,"-print_freq",&tsCtx.print_freq,PETSC_NULL);
188:
189: PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
190: for (i=0; i<param.mglevels; i++) {
191: VecDuplicate(dmmg[i]->x, &(user[i].Xold));
192: VecDuplicate(dmmg[i]->x, &(user[i].func));
193: user[i].tsCtx = &tsCtx;
194: user[i].param = ¶m;
195: dmmg[i]->user = &user[i];
196: }
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Create user context, set problem data, create vector data structures.
199: Also, compute the initial guess.
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201:
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Create nonlinear solver context
204:
205: Process adiC(20): AddTSTermLocal FormFunctionLocal FormFunctionLocali
206: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
208: DMMGSetFromOptions(dmmg);
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: PetscErrorCode Initialize(DMMG *dmmg)
253: {
254: AppCtx *user = (AppCtx*)dmmg[0]->user;
255: DA da;
256: Parameter *param = user->param;
257: PetscInt i,j,mx,xs,ys,xm,ym,mglevel;
259: PetscReal grashof,dx;
260: Field **x;
262: da = (DA)(dmmg[param->mglevels-1]->dm);
263: grashof = user->param->grashof;
265: DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
266: dx = 1.0/(mx-1);
268: /*
269: Get local grid boundaries (for 2-dimensional DA):
270: xs, ys - starting grid indices (no ghost points)
271: xm, ym - widths of local grid (no ghost points)
272: */
273: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
275: /*
276: Get a pointer to vector data.
277: - For default PETSc vectors, VecGetArray() returns a pointer to
278: the data array. Otherwise, the routine is implementation dependent.
279: - You MUST call VecRestoreArray() when you no longer need access to
280: the array.
281: */
282: DAVecGetArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
284: /*
285: Compute initial guess over the locally owned part of the grid
286: Initial condition is motionless fluid and equilibrium temperature
287: */
288: for (j=ys; j<ys+ym; j++) {
289: for (i=xs; i<xs+xm; i++) {
290: x[j][i].u = 0.0;
291: x[j][i].v = 0.0;
292: x[j][i].omega = 0.0;
293: x[j][i].temp = (grashof>0)*i*dx;
294: }
295: }
297: /*
298: Restore vector
299: */
300: DAVecRestoreArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
301: /* Restrict Xold to coarser levels */
302: for (mglevel=param->mglevels-1; mglevel>0; mglevel--) {
303: MatRestrict(dmmg[mglevel]->R, ((AppCtx*)dmmg[mglevel]->user)->Xold, ((AppCtx*)dmmg[mglevel-1]->user)->Xold);
304: VecPointwiseMult(((AppCtx*)dmmg[mglevel-1]->user)->Xold,((AppCtx*)dmmg[mglevel-1]->user)->Xold,dmmg[mglevel]->Rscale);
305: }
306:
307: /* Store X in the qold for time stepping */
308: /*VecDuplicate(X,&tsCtx->qold);
309: VecDuplicate(X,&tsCtx->func);
310: VecCopy(X,tsCtx->Xold);
311: tsCtx->ires = 0;
312: SNESComputeFunction(snes,tsCtx->Xold,tsCtx->func);
313: VecNorm(tsCtx->func,NORM_2,&tsCtx->fnorm_ini);
314: tsCtx->ires = 1;
315: PetscPrintf(PETSC_COMM_WORLD,"Initial function norm is %G\n",tsCtx->fnorm_ini);*/
316: return 0;
317: }
318: /* ------------------------------------------------------------------- */
321: /*
322: FormInitialGuess - Forms initial approximation.
324: Input Parameters:
325: user - user-defined application context
326: X - vector
328: Output Parameter:
329: X - vector
330: */
331: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
332: {
333: AppCtx *user = (AppCtx*)dmmg->user;
334: TstepCtx *tsCtx = user->tsCtx;
337: VecCopy(user->Xold, X);
338: /* calculate the residual on fine mesh */
339: if (user->tsCtx->fnorm_ini == 0.0) {
340: tsCtx->ires = 0;
341: SNESComputeFunction(dmmg->snes,user->Xold,user->func);
342: VecNorm(user->func,NORM_2,&tsCtx->fnorm_ini);
343: tsCtx->ires = 1;
344: tsCtx->cfl = tsCtx->cfl_ini;
345: }
346: return 0;
347: }
349: /*---------------------------------------------------------------------*/
352: PetscErrorCode AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,void *ptr)
353: /*---------------------------------------------------------------------*/
354: {
355: AppCtx *user = (AppCtx*)ptr;
356: TstepCtx *tsCtx = user->tsCtx;
357: DA da = info->da;
359: PetscInt i,j, xints,xinte,yints,yinte;
360: PetscReal hx,hy,dhx,dhy,hxhy;
361: PassiveScalar dtinv;
362: PassiveField **xold;
363: PetscTruth use_parab = tsCtx->use_parabolic;
366: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
367: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
368: hx = 1.0/dhx; hy = 1.0/dhy;
369: hxhy = hx*hy;
370: DAVecGetArray(da,user->Xold,&xold);
371: dtinv = hxhy/(tsCtx->cfl*tsCtx->dt);
372: /*
373: use_parab = PETSC_TRUE for parabolic equations; all the four equations have temporal term.
374: = PETSC_FALSE for differential algebraic equtions (DAE);
375: velocity equations do not have temporal term.
376: */
377: for (j=yints; j<yinte; j++) {
378: for (i=xints; i<xinte; i++) {
379: if (use_parab) {
380: f[j][i].u += dtinv*(x[j][i].u-xold[j][i].u);
381: f[j][i].v += dtinv*(x[j][i].v-xold[j][i].v);
382: }
383: f[j][i].omega += dtinv*(x[j][i].omega-xold[j][i].omega);
384: f[j][i].temp += dtinv*(x[j][i].temp-xold[j][i].temp);
385: }
386: }
387: DAVecRestoreArray(da,user->Xold,&xold);
388: return(0);
389: }
393: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
394: {
395: AppCtx *user = (AppCtx*)ptr;
396: TstepCtx *tsCtx = user->tsCtx;
398: PetscInt xints,xinte,yints,yinte,i,j;
399: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
400: PetscReal grashof,prandtl,lid;
401: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
404: grashof = user->param->grashof;
405: prandtl = user->param->prandtl;
406: lid = user->param->lidvelocity;
408: /*
409: Define mesh intervals ratios for uniform grid.
410: [Note: FD formulae below are normalized by multiplying through by
411: local volume element to obtain coefficients O(1) in two dimensions.]
412: */
413: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
414: hx = 1.0/dhx; hy = 1.0/dhy;
415: hxdhy = hx*dhy; hydhx = hy*dhx;
417: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
419: /* Test whether we are on the bottom edge of the global array */
420: if (yints == 0) {
421: j = 0;
422: yints = yints + 1;
423: /* bottom edge */
424: for (i=info->xs; i<info->xs+info->xm; i++) {
425: f[j][i].u = x[j][i].u;
426: f[j][i].v = x[j][i].v;
427: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
428: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
429: }
430: }
432: /* Test whether we are on the top edge of the global array */
433: if (yinte == info->my) {
434: j = info->my - 1;
435: yinte = yinte - 1;
436: /* top edge */
437: for (i=info->xs; i<info->xs+info->xm; i++) {
438: f[j][i].u = x[j][i].u - lid;
439: f[j][i].v = x[j][i].v;
440: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
441: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
442: }
443: }
445: /* Test whether we are on the left edge of the global array */
446: if (xints == 0) {
447: i = 0;
448: xints = xints + 1;
449: /* left edge */
450: for (j=info->ys; j<info->ys+info->ym; j++) {
451: f[j][i].u = x[j][i].u;
452: f[j][i].v = x[j][i].v;
453: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
454: f[j][i].temp = x[j][i].temp;
455: }
456: }
458: /* Test whether we are on the right edge of the global array */
459: if (xinte == info->mx) {
460: i = info->mx - 1;
461: xinte = xinte - 1;
462: /* right edge */
463: for (j=info->ys; j<info->ys+info->ym; j++) {
464: f[j][i].u = x[j][i].u;
465: f[j][i].v = x[j][i].v;
466: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
467: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
468: }
469: }
471: /* Compute over the interior points */
472: for (j=yints; j<yinte; j++) {
473: for (i=xints; i<xinte; i++) {
475: /*
476: convective coefficients for upwinding
477: */
478: vx = x[j][i].u; avx = PetscAbsScalar(vx);
479: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
480: vy = x[j][i].v; avy = PetscAbsScalar(vy);
481: vyp = .5*(vy+avy); vym = .5*(vy-avy);
483: /* U velocity */
484: u = x[j][i].u;
485: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
486: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
487: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
489: /* V velocity */
490: u = x[j][i].v;
491: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
492: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
493: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
495: /* Omega */
496: u = x[j][i].omega;
497: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
498: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
499: f[j][i].omega = uxx + uyy +
500: (vxp*(u - x[j][i-1].omega) +
501: vxm*(x[j][i+1].omega - u)) * hy +
502: (vyp*(u - x[j-1][i].omega) +
503: vym*(x[j+1][i].omega - u)) * hx -
504: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
506: /* Temperature */
507: u = x[j][i].temp;
508: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
509: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
510: f[j][i].temp = uxx + uyy + prandtl * (
511: (vxp*(u - x[j][i-1].temp) +
512: vxm*(x[j][i+1].temp - u)) * hy +
513: (vyp*(u - x[j-1][i].temp) +
514: vym*(x[j+1][i].temp - u)) * hx);
515: }
516: }
518: /* Add time step contribution */
519: if (tsCtx->ires) {
520: AddTSTermLocal(info,x,f,ptr);
521: }
522: /*
523: Flop count (multiply-adds are counted as 2 operations)
524: */
525: PetscLogFlops(84*info->ym*info->xm);
526: return(0);
527: }
531: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
532: {
533: AppCtx *user = (AppCtx*)ptr;
534: PetscInt i,j,c;
535: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
536: PassiveReal grashof,prandtl,lid;
537: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
540: grashof = user->param->grashof;
541: prandtl = user->param->prandtl;
542: lid = user->param->lidvelocity;
544: /*
545: Define mesh intervals ratios for uniform grid.
546: [Note: FD formulae below are normalized by multiplying through by
547: local volume element to obtain coefficients O(1) in two dimensions.]
548: */
549: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
550: hx = 1.0/dhx; hy = 1.0/dhy;
551: hxdhy = hx*dhy; hydhx = hy*dhx;
553: i = st->i; j = st->j; c = st->c;
555: /* Test whether we are on the right edge of the global array */
556: if (i == info->mx-1) {
557: if (c == 0) *f = x[j][i].u;
558: else if (c == 1) *f = x[j][i].v;
559: else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
560: else *f = x[j][i].temp - (PetscReal)(grashof>0);
562: /* Test whether we are on the left edge of the global array */
563: } else if (i == 0) {
564: if (c == 0) *f = x[j][i].u;
565: else if (c == 1) *f = x[j][i].v;
566: else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
567: else *f = x[j][i].temp;
569: /* Test whether we are on the top edge of the global array */
570: } else if (j == info->my-1) {
571: if (c == 0) *f = x[j][i].u - lid;
572: else if (c == 1) *f = x[j][i].v;
573: else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
574: else *f = x[j][i].temp-x[j-1][i].temp;
576: /* Test whether we are on the bottom edge of the global array */
577: } else if (j == 0) {
578: if (c == 0) *f = x[j][i].u;
579: else if (c == 1) *f = x[j][i].v;
580: else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
581: else *f = x[j][i].temp-x[j+1][i].temp;
583: /* Compute over the interior points */
584: } else {
585: /*
586: convective coefficients for upwinding
587: */
588: vx = x[j][i].u; avx = PetscAbsScalar(vx);
589: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
590: vy = x[j][i].v; avy = PetscAbsScalar(vy);
591: vyp = .5*(vy+avy); vym = .5*(vy-avy);
593: /* U velocity */
594: if (c == 0) {
595: u = x[j][i].u;
596: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
597: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
598: *f = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
600: /* V velocity */
601: } else if (c == 1) {
602: u = x[j][i].v;
603: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
604: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
605: *f = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
606:
607: /* Omega */
608: } else if (c == 2) {
609: u = x[j][i].omega;
610: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
611: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
612: *f = uxx + uyy +
613: (vxp*(u - x[j][i-1].omega) +
614: vxm*(x[j][i+1].omega - u)) * hy +
615: (vyp*(u - x[j-1][i].omega) +
616: vym*(x[j+1][i].omega - u)) * hx -
617: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
618:
619: /* Temperature */
620: } else {
621: u = x[j][i].temp;
622: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
623: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
624: *f = uxx + uyy + prandtl * (
625: (vxp*(u - x[j][i-1].temp) +
626: vxm*(x[j][i+1].temp - u)) * hy +
627: (vyp*(u - x[j-1][i].temp) +
628: vym*(x[j+1][i].temp - u)) * hx);
629: }
630: }
632: return(0);
633: }
636: /*---------------------------------------------------------------------*/
639: PetscErrorCode Update(DMMG *dmmg)
640: /*---------------------------------------------------------------------*/
641: {
642: AppCtx *user = (AppCtx *) ((dmmg[0])->user);
643: TstepCtx *tsCtx = user->tsCtx;
644: Parameter *param = user->param;
645: SNES snes;
647: PetscScalar fratio;
648: PetscScalar time1,time2,cpuloc;
649: PetscInt max_steps,its;
650: PetscTruth print_flag = PETSC_FALSE;
651: PetscInt nfailsCum = 0,nfails = 0;
654: /* Note: print_flag displays diagnostic info, not convergence behavior. Use '-snes_monitor' for converges info. */
655: PetscOptionsHasName(PETSC_NULL,"-print",&print_flag);
656: if (user->param->PreLoading)
657: max_steps = 1;
658: else
659: max_steps = tsCtx->max_steps;
660: fratio = 1.0;
661:
662: PetscGetTime(&time1);
663: for (tsCtx->itstep = 0; (tsCtx->itstep < max_steps) &&
664: (fratio <= tsCtx->fnorm_ratio); tsCtx->itstep++) {
665: DMMGSolve(dmmg);
666: snes = DMMGGetSNES(dmmg);
667: SNESGetIterationNumber(snes,&its);
668: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n", its);
669: SNESGetNonlinearStepFailures(snes,&nfails);
670: nfailsCum += nfails; nfails = 0;
671: if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
672: /*tsCtx->qcur = DMMGGetx(dmmg);
673: VecCopy(tsCtx->qcur,tsCtx->qold);*/
674:
675: VecCopy(dmmg[param->mglevels-1]->x, ((AppCtx*)dmmg[param->mglevels-1]->user)->Xold);
676: for (its=param->mglevels-1; its>0 ;its--) {
677: MatRestrict(dmmg[its]->R, ((AppCtx*)dmmg[its]->user)->Xold, ((AppCtx*)dmmg[its-1]->user)->Xold);
678: VecPointwiseMult(((AppCtx*)dmmg[its-1]->user)->Xold,((AppCtx*)dmmg[its-1]->user)->Xold,dmmg[its]->Rscale);
679: }
681:
682: ComputeTimeStep(snes,((AppCtx*)dmmg[param->mglevels-1]->user));
683: if (print_flag) {
684: PetscPrintf(PETSC_COMM_WORLD,"At Time Step %D cfl = %G and fnorm = %G\n",
685: tsCtx->itstep,tsCtx->cfl,tsCtx->fnorm);
686: PetscPrintf(PETSC_COMM_WORLD,"Wall clock time needed %G seconds for %D time steps\n",
687: cpuloc,tsCtx->itstep);
688: }
689: fratio = tsCtx->fnorm_ini/tsCtx->fnorm;
690: PetscGetTime(&time2);
691: cpuloc = time2-time1;
692: MPI_Barrier(PETSC_COMM_WORLD);
693: } /* End of time step loop */
695: if (print_flag) {
696: PetscPrintf(PETSC_COMM_WORLD,"Total wall clock time needed %G seconds for %D time steps\n",
697: cpuloc,tsCtx->itstep);
698: PetscPrintf(PETSC_COMM_WORLD,"cfl = %G fnorm = %G\n",tsCtx->cfl,tsCtx->fnorm);
699: }
700: if (user->param->PreLoading) {
701: tsCtx->fnorm_ini = 0.0;
702: PetscPrintf(PETSC_COMM_WORLD,"Preloading done ...\n");
703: }
704: /*
705: {
706: Vec xx,yy;
707: PetscScalar fnorm,fnorm1;
708: SNESGetFunctionNorm(snes,&fnorm);
709: xx = DMMGGetx(dmmg);
710: VecDuplicate(xx,&yy);
711: SNESComputeFunction(snes,xx,yy);
712: VecNorm(yy,NORM_2,&fnorm1);
713: PetscPrintf(PETSC_COMM_WORLD,"fnorm = %G, fnorm1 = %G\n",fnorm,fnorm1);
714:
715: }
716: */
718: return(0);
719: }
721: /*---------------------------------------------------------------------*/
724: PetscErrorCode ComputeTimeStep(SNES snes,void *ptr)
725: /*---------------------------------------------------------------------*/
726: {
727: AppCtx *user = (AppCtx*)ptr;
728: TstepCtx *tsCtx = user->tsCtx;
729: Vec func = user->func;
730: PetscScalar inc = 1.1, newcfl;
732: /*int iramp = tsCtx->iramp;*/
733:
735: tsCtx->dt = 0.01;
736: PetscOptionsGetScalar(PETSC_NULL,"-deltat",&tsCtx->dt,PETSC_NULL);
737: tsCtx->ires = 0;
738: SNESComputeFunction(snes,user->Xold,user->func);
739: tsCtx->ires = 1;
740: VecNorm(func,NORM_2,&tsCtx->fnorm);
741: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
742: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
743: /* first time through so compute initial function norm */
744: /*if (tsCtx->fnorm_ini == 0.0) {
745: tsCtx->fnorm_ini = tsCtx->fnorm;
746: tsCtx->cfl = tsCtx->cfl_ini;
747: } else {
748: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
749: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
750: }*/
751: return(0);
752: }