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",&param.lidvelocity,PETSC_NULL);
164:     PetscOptionsGetReal(PETSC_NULL,"-prandtl",&param.prandtl,PETSC_NULL);
165:     PetscOptionsGetReal(PETSC_NULL,"-grashof",&param.grashof,PETSC_NULL);
166:     PetscOptionsHasName(PETSC_NULL,"-contours",&param.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 = &param;
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: }