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: DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
209:
210: PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, grashof # = %G\n",
211: param.lidvelocity,param.prandtl,param.grashof);
212:
213:
214: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: Solve the nonlinear system
216: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217: DMMGSetInitialGuess(dmmg,FormInitialGuess);
218:
219: PreLoadStage("Solve");
220: Initialize(dmmg);
221: Update(dmmg);
222: /*
223: Visualize solution
224: */
225:
226: if (param.draw_contours) {
227: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
228: }
229: /*VecView(DMMGGetx(dmmg),PETSC_VIEWER_STDOUT_WORLD);*/
230: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: Free work space. All PETSc objects should be destroyed when they
232: are no longer needed.
233: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234:
235: for (i=0; i<param.mglevels; i++) {
236: VecDestroy(user[i].Xold);
237: VecDestroy(user[i].func);
238: }
239: PetscFree(user);
240: DMMGDestroy(dmmg);
241: PreLoadEnd();
242:
243: PetscFinalize();
244: return 0;
245: }
247: /* ------------------------------------------------------------------- */
250: /* ------------------------------------------------------------------- */
251: PetscErrorCode Initialize(DMMG *dmmg)
252: {
253: AppCtx *user = (AppCtx*)dmmg[0]->user;
254: DA da;
255: Parameter *param = user->param;
256: PetscInt i,j,mx,xs,ys,xm,ym,mglevel;
258: PetscReal grashof,dx;
259: Field **x;
261: da = (DA)(dmmg[param->mglevels-1]->dm);
262: grashof = user->param->grashof;
264: DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
265: dx = 1.0/(mx-1);
267: /*
268: Get local grid boundaries (for 2-dimensional DA):
269: xs, ys - starting grid indices (no ghost points)
270: xm, ym - widths of local grid (no ghost points)
271: */
272: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
274: /*
275: Get a pointer to vector data.
276: - For default PETSc vectors, VecGetArray() returns a pointer to
277: the data array. Otherwise, the routine is implementation dependent.
278: - You MUST call VecRestoreArray() when you no longer need access to
279: the array.
280: */
281: DAVecGetArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
283: /*
284: Compute initial guess over the locally owned part of the grid
285: Initial condition is motionless fluid and equilibrium temperature
286: */
287: for (j=ys; j<ys+ym; j++) {
288: for (i=xs; i<xs+xm; i++) {
289: x[j][i].u = 0.0;
290: x[j][i].v = 0.0;
291: x[j][i].omega = 0.0;
292: x[j][i].temp = (grashof>0)*i*dx;
293: }
294: }
296: /*
297: Restore vector
298: */
299: DAVecRestoreArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
300: /* Restrict Xold to coarser levels */
301: for (mglevel=param->mglevels-1; mglevel>0; mglevel--) {
302: MatRestrict(dmmg[mglevel]->R, ((AppCtx*)dmmg[mglevel]->user)->Xold, ((AppCtx*)dmmg[mglevel-1]->user)->Xold);
303: VecPointwiseMult(((AppCtx*)dmmg[mglevel-1]->user)->Xold,((AppCtx*)dmmg[mglevel-1]->user)->Xold,dmmg[mglevel]->Rscale);
304: }
305:
306: /* Store X in the qold for time stepping */
307: /*VecDuplicate(X,&tsCtx->qold);
308: VecDuplicate(X,&tsCtx->func);
309: VecCopy(X,tsCtx->Xold);
310: tsCtx->ires = 0;
311: SNESComputeFunction(snes,tsCtx->Xold,tsCtx->func);
312: VecNorm(tsCtx->func,NORM_2,&tsCtx->fnorm_ini);
313: tsCtx->ires = 1;
314: PetscPrintf(PETSC_COMM_WORLD,"Initial function norm is %G\n",tsCtx->fnorm_ini);*/
315: return 0;
316: }
317: /* ------------------------------------------------------------------- */
320: /*
321: FormInitialGuess - Forms initial approximation.
323: Input Parameters:
324: user - user-defined application context
325: X - vector
327: Output Parameter:
328: X - vector
329: */
330: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
331: {
332: AppCtx *user = (AppCtx*)dmmg->user;
333: TstepCtx *tsCtx = user->tsCtx;
336: VecCopy(user->Xold, X);
337: /* calculate the residual on fine mesh */
338: if (user->tsCtx->fnorm_ini == 0.0) {
339: tsCtx->ires = 0;
340: SNESComputeFunction(dmmg->snes,user->Xold,user->func);
341: VecNorm(user->func,NORM_2,&tsCtx->fnorm_ini);
342: tsCtx->ires = 1;
343: tsCtx->cfl = tsCtx->cfl_ini;
344: }
345: return 0;
346: }
348: /*---------------------------------------------------------------------*/
351: PetscErrorCode AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,void *ptr)
352: /*---------------------------------------------------------------------*/
353: {
354: AppCtx *user = (AppCtx*)ptr;
355: TstepCtx *tsCtx = user->tsCtx;
356: DA da = info->da;
358: PetscInt i,j, xints,xinte,yints,yinte;
359: PetscReal hx,hy,dhx,dhy,hxhy;
360: PassiveScalar dtinv;
361: PassiveField **xold;
362: PetscTruth use_parab = tsCtx->use_parabolic;
365: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
366: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
367: hx = 1.0/dhx; hy = 1.0/dhy;
368: hxhy = hx*hy;
369: DAVecGetArray(da,user->Xold,&xold);
370: dtinv = hxhy/(tsCtx->cfl*tsCtx->dt);
371: /*
372: use_parab = PETSC_TRUE for parabolic equations; all the four equations have temporal term.
373: = PETSC_FALSE for differential algebraic equtions (DAE);
374: velocity equations do not have temporal term.
375: */
376: for (j=yints; j<yinte; j++) {
377: for (i=xints; i<xinte; i++) {
378: if (use_parab) {
379: f[j][i].u += dtinv*(x[j][i].u-xold[j][i].u);
380: f[j][i].v += dtinv*(x[j][i].v-xold[j][i].v);
381: }
382: f[j][i].omega += dtinv*(x[j][i].omega-xold[j][i].omega);
383: f[j][i].temp += dtinv*(x[j][i].temp-xold[j][i].temp);
384: }
385: }
386: DAVecRestoreArray(da,user->Xold,&xold);
387: return(0);
388: }
392: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
393: {
394: AppCtx *user = (AppCtx*)ptr;
395: TstepCtx *tsCtx = user->tsCtx;
397: PetscInt xints,xinte,yints,yinte,i,j;
398: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
399: PetscReal grashof,prandtl,lid;
400: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
403: grashof = user->param->grashof;
404: prandtl = user->param->prandtl;
405: lid = user->param->lidvelocity;
407: /*
408: Define mesh intervals ratios for uniform grid.
409: [Note: FD formulae below are normalized by multiplying through by
410: local volume element to obtain coefficients O(1) in two dimensions.]
411: */
412: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
413: hx = 1.0/dhx; hy = 1.0/dhy;
414: hxdhy = hx*dhy; hydhx = hy*dhx;
416: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
418: /* Test whether we are on the bottom edge of the global array */
419: if (yints == 0) {
420: j = 0;
421: yints = yints + 1;
422: /* bottom edge */
423: for (i=info->xs; i<info->xs+info->xm; i++) {
424: f[j][i].u = x[j][i].u;
425: f[j][i].v = x[j][i].v;
426: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
427: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
428: }
429: }
431: /* Test whether we are on the top edge of the global array */
432: if (yinte == info->my) {
433: j = info->my - 1;
434: yinte = yinte - 1;
435: /* top edge */
436: for (i=info->xs; i<info->xs+info->xm; i++) {
437: f[j][i].u = x[j][i].u - lid;
438: f[j][i].v = x[j][i].v;
439: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
440: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
441: }
442: }
444: /* Test whether we are on the left edge of the global array */
445: if (xints == 0) {
446: i = 0;
447: xints = xints + 1;
448: /* left edge */
449: for (j=info->ys; j<info->ys+info->ym; j++) {
450: f[j][i].u = x[j][i].u;
451: f[j][i].v = x[j][i].v;
452: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
453: f[j][i].temp = x[j][i].temp;
454: }
455: }
457: /* Test whether we are on the right edge of the global array */
458: if (xinte == info->mx) {
459: i = info->mx - 1;
460: xinte = xinte - 1;
461: /* right edge */
462: for (j=info->ys; j<info->ys+info->ym; j++) {
463: f[j][i].u = x[j][i].u;
464: f[j][i].v = x[j][i].v;
465: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
466: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
467: }
468: }
470: /* Compute over the interior points */
471: for (j=yints; j<yinte; j++) {
472: for (i=xints; i<xinte; i++) {
474: /*
475: convective coefficients for upwinding
476: */
477: vx = x[j][i].u; avx = PetscAbsScalar(vx);
478: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
479: vy = x[j][i].v; avy = PetscAbsScalar(vy);
480: vyp = .5*(vy+avy); vym = .5*(vy-avy);
482: /* U velocity */
483: u = x[j][i].u;
484: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
485: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
486: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
488: /* V velocity */
489: u = x[j][i].v;
490: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
491: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
492: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
494: /* Omega */
495: u = x[j][i].omega;
496: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
497: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
498: f[j][i].omega = uxx + uyy +
499: (vxp*(u - x[j][i-1].omega) +
500: vxm*(x[j][i+1].omega - u)) * hy +
501: (vyp*(u - x[j-1][i].omega) +
502: vym*(x[j+1][i].omega - u)) * hx -
503: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
505: /* Temperature */
506: u = x[j][i].temp;
507: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
508: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
509: f[j][i].temp = uxx + uyy + prandtl * (
510: (vxp*(u - x[j][i-1].temp) +
511: vxm*(x[j][i+1].temp - u)) * hy +
512: (vyp*(u - x[j-1][i].temp) +
513: vym*(x[j+1][i].temp - u)) * hx);
514: }
515: }
517: /* Add time step contribution */
518: if (tsCtx->ires) {
519: AddTSTermLocal(info,x,f,ptr);
520: }
521: /*
522: Flop count (multiply-adds are counted as 2 operations)
523: */
524: PetscLogFlops(84*info->ym*info->xm);
525: return(0);
526: }
530: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
531: {
532: AppCtx *user = (AppCtx*)ptr;
533: PetscInt i,j,c;
534: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
535: PassiveReal grashof,prandtl,lid;
536: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
539: grashof = user->param->grashof;
540: prandtl = user->param->prandtl;
541: lid = user->param->lidvelocity;
543: /*
544: Define mesh intervals ratios for uniform grid.
545: [Note: FD formulae below are normalized by multiplying through by
546: local volume element to obtain coefficients O(1) in two dimensions.]
547: */
548: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
549: hx = 1.0/dhx; hy = 1.0/dhy;
550: hxdhy = hx*dhy; hydhx = hy*dhx;
552: i = st->i; j = st->j; c = st->c;
554: /* Test whether we are on the right edge of the global array */
555: if (i == info->mx-1) {
556: if (c == 0) *f = x[j][i].u;
557: else if (c == 1) *f = x[j][i].v;
558: else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
559: else *f = x[j][i].temp - (PetscReal)(grashof>0);
561: /* Test whether we are on the left edge of the global array */
562: } else if (i == 0) {
563: if (c == 0) *f = x[j][i].u;
564: else if (c == 1) *f = x[j][i].v;
565: else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
566: else *f = x[j][i].temp;
568: /* Test whether we are on the top edge of the global array */
569: } else if (j == info->my-1) {
570: if (c == 0) *f = x[j][i].u - lid;
571: else if (c == 1) *f = x[j][i].v;
572: else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
573: else *f = x[j][i].temp-x[j-1][i].temp;
575: /* Test whether we are on the bottom edge of the global array */
576: } else if (j == 0) {
577: if (c == 0) *f = x[j][i].u;
578: else if (c == 1) *f = x[j][i].v;
579: else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
580: else *f = x[j][i].temp-x[j+1][i].temp;
582: /* Compute over the interior points */
583: } else {
584: /*
585: convective coefficients for upwinding
586: */
587: vx = x[j][i].u; avx = PetscAbsScalar(vx);
588: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
589: vy = x[j][i].v; avy = PetscAbsScalar(vy);
590: vyp = .5*(vy+avy); vym = .5*(vy-avy);
592: /* U velocity */
593: if (c == 0) {
594: u = x[j][i].u;
595: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
596: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
597: *f = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
599: /* V velocity */
600: } else if (c == 1) {
601: u = x[j][i].v;
602: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
603: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
604: *f = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
605:
606: /* Omega */
607: } else if (c == 2) {
608: u = x[j][i].omega;
609: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
610: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
611: *f = uxx + uyy +
612: (vxp*(u - x[j][i-1].omega) +
613: vxm*(x[j][i+1].omega - u)) * hy +
614: (vyp*(u - x[j-1][i].omega) +
615: vym*(x[j+1][i].omega - u)) * hx -
616: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
617:
618: /* Temperature */
619: } else {
620: u = x[j][i].temp;
621: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
622: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
623: *f = uxx + uyy + prandtl * (
624: (vxp*(u - x[j][i-1].temp) +
625: vxm*(x[j][i+1].temp - u)) * hy +
626: (vyp*(u - x[j-1][i].temp) +
627: vym*(x[j+1][i].temp - u)) * hx);
628: }
629: }
631: return(0);
632: }
635: /*---------------------------------------------------------------------*/
638: PetscErrorCode Update(DMMG *dmmg)
639: /*---------------------------------------------------------------------*/
640: {
641: AppCtx *user = (AppCtx *) ((dmmg[0])->user);
642: TstepCtx *tsCtx = user->tsCtx;
643: Parameter *param = user->param;
644: SNES snes;
646: PetscScalar fratio;
647: PetscScalar time1,time2,cpuloc;
648: PetscInt max_steps,its;
649: PetscTruth print_flag = PETSC_FALSE;
650: PetscInt nfailsCum = 0,nfails = 0;
654: PetscOptionsHasName(PETSC_NULL,"-print",&print_flag);
655: if (user->param->PreLoading)
656: max_steps = 1;
657: else
658: max_steps = tsCtx->max_steps;
659: fratio = 1.0;
660:
661: PetscGetTime(&time1);
662: for (tsCtx->itstep = 0; (tsCtx->itstep < max_steps) &&
663: (fratio <= tsCtx->fnorm_ratio); tsCtx->itstep++) {
664: DMMGSolve(dmmg);
665: snes = DMMGGetSNES(dmmg);
666: SNESGetIterationNumber(snes,&its);
667: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n", its);
668: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
669: nfailsCum += nfails; nfails = 0;
670: if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
671: /*tsCtx->qcur = DMMGGetx(dmmg);
672: VecCopy(tsCtx->qcur,tsCtx->qold);*/
674: VecCopy(dmmg[param->mglevels-1]->x, ((AppCtx*)dmmg[param->mglevels-1]->user)->Xold);
675: for (its=param->mglevels-1; its>0 ;its--) {
676: MatRestrict(dmmg[its]->R, ((AppCtx*)dmmg[its]->user)->Xold, ((AppCtx*)dmmg[its-1]->user)->Xold);
677: VecPointwiseMult(((AppCtx*)dmmg[its-1]->user)->Xold,((AppCtx*)dmmg[its-1]->user)->Xold,dmmg[its]->Rscale);
678: }
680:
681: ComputeTimeStep(snes,((AppCtx*)dmmg[param->mglevels-1]->user));
682: if (print_flag) {
683: PetscPrintf(PETSC_COMM_WORLD,"At Time Step %D cfl = %G and fnorm = %G\n",
684: tsCtx->itstep,tsCtx->cfl,tsCtx->fnorm);
685: PetscPrintf(PETSC_COMM_WORLD,"Wall clock time needed %G seconds for %D time steps\n",
686: cpuloc,tsCtx->itstep);
687: }
688: fratio = tsCtx->fnorm_ini/tsCtx->fnorm;
689: PetscGetTime(&time2);
690: cpuloc = time2-time1;
691: MPI_Barrier(PETSC_COMM_WORLD);
692: } /* End of time step loop */
693:
694: PetscPrintf(PETSC_COMM_WORLD,"Total wall clock time needed %G seconds for %D time steps\n",
695: cpuloc,tsCtx->itstep);
696: PetscPrintf(PETSC_COMM_WORLD,"cfl = %G fnorm = %G\n",tsCtx->cfl,tsCtx->fnorm);
697: if (user->param->PreLoading) {
698: tsCtx->fnorm_ini = 0.0;
699: PetscPrintf(PETSC_COMM_WORLD,"Preloading done ...\n");
700: }
701: /*
702: {
703: Vec xx,yy;
704: PetscScalar fnorm,fnorm1;
705: SNESGetFunctionNorm(snes,&fnorm);
706: xx = DMMGGetx(dmmg);
707: VecDuplicate(xx,&yy);
708: SNESComputeFunction(snes,xx,yy);
709: VecNorm(yy,NORM_2,&fnorm1);
710: PetscPrintf(PETSC_COMM_WORLD,"fnorm = %G, fnorm1 = %G\n",fnorm,fnorm1);
711:
712: }
713: */
715: return(0);
716: }
718: /*---------------------------------------------------------------------*/
721: PetscErrorCode ComputeTimeStep(SNES snes,void *ptr)
722: /*---------------------------------------------------------------------*/
723: {
724: AppCtx *user = (AppCtx*)ptr;
725: TstepCtx *tsCtx = user->tsCtx;
726: Vec func = user->func;
727: PetscScalar inc = 1.1, newcfl;
729: /*int iramp = tsCtx->iramp;*/
730:
732: tsCtx->dt = 0.01;
733: PetscOptionsGetScalar(PETSC_NULL,"-deltat",&tsCtx->dt,PETSC_NULL);
734: tsCtx->ires = 0;
735: SNESComputeFunction(snes,user->Xold,user->func);
736: tsCtx->ires = 1;
737: VecNorm(func,NORM_2,&tsCtx->fnorm);
738: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
739: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
740: /* first time through so compute initial function norm */
741: /*if (tsCtx->fnorm_ini == 0.0) {
742: tsCtx->fnorm_ini = tsCtx->fnorm;
743: tsCtx->cfl = tsCtx->cfl_ini;
744: } else {
745: newcfl = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
746: tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
747: }*/
748: return(0);
749: }