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:     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: }