Actual source code: ex27.c

  2: static char help[] = "Nonlinear driven cavity with multigrid and pseudo timestepping 2d.\n\
  3:   \n\
  4: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
  5: The flow can be driven with the lid or with bouyancy or both:\n\
  6:   -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
  7:   -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
  8:   -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
  9:   -contours : draw contour plots of solution\n\n";

 11: /*T
 12:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
 13:    Concepts: DA^using distributed arrays;
 14:    Concepts: multicomponent
 15:    Processors: n
 16: T*/

 18: /* ------------------------------------------------------------------------

 20:     We thank David E. Keyes for contributing the driven cavity discretization
 21:     within this example code.

 23:     This problem is modeled by the partial differential equation system
 24:   
 25:         dU/dt          - Lap(U)     - Grad_y(Omega) = 0
 26:         dV/dt          - Lap(V)     + Grad_x(Omega) = 0
 27:         d(omega)/dt    - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
 28:         dT/dt          - Lap(T)     + PR*Div([U*T,V*T]) = 0

 30:     in the unit square, which is uniformly discretized in each of x and
 31:     y in this simple encoding.

 33:     No-slip, rigid-wall Dirichlet conditions are used for [U,V].
 34:     Dirichlet conditions are used for Omega, based on the definition of
 35:     vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
 36:     constant coordinate boundary, the tangential derivative is zero.
 37:     Dirichlet conditions are used for T on the left and right walls,
 38:     and insulation homogeneous Neumann conditions are used for T on
 39:     the top and bottom walls. 

 41:     A finite difference approximation with the usual 5-point stencil 
 42:     is used to discretize the boundary value problem to obtain a 
 43:     nonlinear system of equations.  Upwinding is used for the divergence
 44:     (convective) terms and central for the gradient (source) terms.
 45:     
 46:     The Jacobian can be either
 47:       * formed via finite differencing using coloring (the default), or
 48:       * applied matrix-free via the option -snes_mf 
 49:         (for larger grid problems this variant may not converge 
 50:         without a preconditioner due to ill-conditioning).


 53:     This example is used in the paper

 55:         Todd S. Coffey and C. T. Kelley and David E. Keyes,
 56:         Pseudotransient Continuation and Differential-Algebraic Equations, 2003.

 58:     however Figure 3.1 was generated with a slightly different algorithm (see
 59:     targets runex27 and runex27_p) than described in the paper.  In particular,
 60:     the described algorithm is linearly implicit, advancing to the next step
 61:     after one Newton step, so that the steady state residual is always used, but
 62:     the figure was generated by converging each step to a relative tolerance of
 63:     1.e-3.  On the example problem, setting -snes_max_it 1 has only minor impact
 64:     on number of steps, but significantly reduces the required number of linear
 65:     solves.

 67:   ------------------------------------------------------------------------- */

 69: /* 
 70:    Include "petscda.h" so that we can use distributed arrays (DAs).
 71:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 72:    file automatically includes:
 73:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 74:      petscmat.h - matrices
 75:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 76:      petscviewer.h - viewers               petscpc.h  - preconditioners
 77:      petscksp.h   - linear solvers 
 78: */
 79:  #include petscsnes.h
 80:  #include petscda.h
 81:  #include petscdmmg.h

 83: /* 
 84:    User-defined routines and data structures
 85: */

 87: typedef struct {
 88:   PassiveScalar  fnorm_ini,dt_ini,cfl_ini;
 89:   PassiveScalar  fnorm,dt,cfl;
 90:   PassiveScalar  ptime;
 91:   PassiveScalar  cfl_max,max_time;
 92:   PassiveScalar  fnorm_ratio;
 93:   PetscInt       ires,iramp,itstep;
 94:   PetscInt       max_steps,print_freq;
 95:   PetscInt       LocalTimeStepping;
 96:   PetscTruth     use_parabolic;
 97: } TstepCtx;

 99: /*
100:     The next two structures are essentially the same. The difference is that
101:   the first does not contain derivative information (as used by ADIC) while the
102:   second does. The first is used only to contain the solution at the previous time-step
103:   which is a constant in the computation of the current time step and hence passive 
104:   (meaning not active in terms of derivatives).
105: */
106: typedef struct {
107:   PassiveScalar u,v,omega,temp;
108: } PassiveField;

110: typedef struct {
111:   PetscScalar u,v,omega,temp;
112: } Field;


115: typedef struct {
116:   PetscInt     mglevels;
117:   PetscInt     cycles;                       /* numbers of time steps for integration */
118:   PassiveReal  lidvelocity,prandtl,grashof;  /* physical parameters */
119:   PetscTruth   draw_contours;                /* indicates drawing contours of solution */
120:   PetscTruth   PreLoading;
121: } Parameter;

123: typedef struct {
124:   Vec          Xold,func;
125:   TstepCtx     *tsCtx;
126:   Parameter    *param;
127: } AppCtx;



140: int main(int argc,char **argv)
141: {
142:   DMMG           *dmmg;               /* multilevel grid structure */
143:   AppCtx         *user;                /* user-defined work context */
144:   TstepCtx       tsCtx;
145:   Parameter      param;
146:   PetscInt       mx,my,i;
148:   MPI_Comm       comm;
149:   DA             da;

151:   PetscInitialize(&argc,&argv,(char *)0,help);
152:   comm = PETSC_COMM_WORLD;


155:   PreLoadBegin(PETSC_TRUE,"SetUp");

157:     param.PreLoading = PreLoading;
158:     DMMGCreate(comm,2,&user,&dmmg);
159:     param.mglevels = DMMGGetLevels(dmmg);


162:     /*
163:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
164:       for principal unknowns (x) and governing residuals (f)
165:     */
166:     DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
167:     DMMGSetDM(dmmg,(DM)da);
168:     DADestroy(da);

170:     DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
171:                      PETSC_IGNORE,PETSC_IGNORE);
172:     /* 
173:      Problem parameters (velocity of lid, prandtl, and grashof numbers)
174:     */
175:     param.lidvelocity = 1.0/(mx*my);
176:     param.prandtl     = 1.0;
177:     param.grashof     = 1.0;
178:     PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&param.lidvelocity,PETSC_NULL);
179:     PetscOptionsGetReal(PETSC_NULL,"-prandtl",&param.prandtl,PETSC_NULL);
180:     PetscOptionsGetReal(PETSC_NULL,"-grashof",&param.grashof,PETSC_NULL);
181:     PetscOptionsHasName(PETSC_NULL,"-contours",&param.draw_contours);

183:     DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
184:     DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
185:     DASetFieldName(DMMGGetDA(dmmg),2,"Omega");
186:     DASetFieldName(DMMGGetDA(dmmg),3,"temperature");

188:     /*======================================================================*/
189:     /* Initilize stuff related to time stepping */
190:     /*======================================================================*/
191:     tsCtx.fnorm_ini = 0.0;  tsCtx.cfl_ini     = 50.0;    tsCtx.cfl_max = 1.0e+06;
192:     tsCtx.max_steps = 50;   tsCtx.max_time    = 1.0e+12; tsCtx.iramp   = -50;
193:     tsCtx.dt        = 0.01; tsCtx.fnorm_ratio = 1.0e+10;
194:     tsCtx.LocalTimeStepping = 0;
195:     tsCtx.use_parabolic     = PETSC_FALSE;
196:     PetscOptionsGetTruth(PETSC_NULL,"-use_parabolic",&tsCtx.use_parabolic,PETSC_IGNORE);
197:     PetscOptionsGetInt(PETSC_NULL,"-max_st",&tsCtx.max_steps,PETSC_NULL);
198:     PetscOptionsGetReal(PETSC_NULL,"-ts_rtol",&tsCtx.fnorm_ratio,PETSC_NULL);
199:     PetscOptionsGetReal(PETSC_NULL,"-cfl_ini",&tsCtx.cfl_ini,PETSC_NULL);
200:     PetscOptionsGetReal(PETSC_NULL,"-cfl_max",&tsCtx.cfl_max,PETSC_NULL);
201:     tsCtx.print_freq = tsCtx.max_steps;
202:     PetscOptionsGetInt(PETSC_NULL,"-print_freq",&tsCtx.print_freq,PETSC_NULL);
203: 
204:     PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
205:     for (i=0; i<param.mglevels; i++) {
206:       VecDuplicate(dmmg[i]->x, &(user[i].Xold));
207:       VecDuplicate(dmmg[i]->x, &(user[i].func));
208:       user[i].tsCtx = &tsCtx;
209:       user[i].param = &param;
210:       dmmg[i]->user = &user[i];
211:     }
212:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213:        Create user context, set problem data, create vector data structures.
214:        Also, compute the initial guess.
215:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: 
217:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218:        Create nonlinear solver context
219:        
220:        Process adiC(20):  AddTSTermLocal FormFunctionLocal FormFunctionLocali
221:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
223:     DMMGSetFromOptions(dmmg);
224:     DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
225: 
226:     PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, grashof # = %G\n",
227:                        param.lidvelocity,param.prandtl,param.grashof);
228: 
229: 
230:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231:        Solve the nonlinear system
232:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233:     DMMGSetInitialGuess(dmmg,FormInitialGuess);
234: 
235:     PreLoadStage("Solve");
236:     Initialize(dmmg);
237:     Update(dmmg);
238:     /*
239:       Visualize solution
240:     */
241: 
242:     if (param.draw_contours) {
243:       VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
244:     }
245:     /*VecView(DMMGGetx(dmmg),PETSC_VIEWER_STDOUT_WORLD);*/
246:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247:        Free work space.  All PETSc objects should be destroyed when they
248:        are no longer needed.
249:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
250: 
251:     for (i=0; i<param.mglevels; i++) {
252:       VecDestroy(user[i].Xold);
253:       VecDestroy(user[i].func);
254:     }
255:     PetscFree(user);
256:     DMMGDestroy(dmmg);
257:     PreLoadEnd();
258: 
259:   PetscFinalize();
260:   return 0;
261: }

263: /* ------------------------------------------------------------------- */
266: /* ------------------------------------------------------------------- */
267: PetscErrorCode Initialize(DMMG *dmmg)
268: {
269:   AppCtx         *user = (AppCtx*)dmmg[0]->user;
270:   DA             da;
271:   Parameter      *param = user->param;
272:   PetscInt       i,j,mx,xs,ys,xm,ym,mglevel;
274:   PetscReal      grashof,dx;
275:   Field          **x;

277:   da = (DA)(dmmg[param->mglevels-1]->dm);
278:   grashof = user->param->grashof;

280:   DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
281:   dx  = 1.0/(mx-1);

283:   /*
284:      Get local grid boundaries (for 2-dimensional DA):
285:        xs, ys   - starting grid indices (no ghost points)
286:        xm, ym   - widths of local grid (no ghost points)
287:   */
288:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

290:   /*
291:      Get a pointer to vector data.
292:        - For default PETSc vectors, VecGetArray() returns a pointer to
293:          the data array.  Otherwise, the routine is implementation dependent.
294:        - You MUST call VecRestoreArray() when you no longer need access to
295:          the array.
296:   */
297:   DAVecGetArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);

299:   /*
300:      Compute initial guess over the locally owned part of the grid
301:      Initial condition is motionless fluid and equilibrium temperature
302:   */
303:   for (j=ys; j<ys+ym; j++) {
304:     for (i=xs; i<xs+xm; i++) {
305:       x[j][i].u     = 0.0;
306:       x[j][i].v     = 0.0;
307:       x[j][i].omega = 0.0;
308:       x[j][i].temp  = (grashof>0)*i*dx;
309:     }
310:   }

312:   /*
313:      Restore vector
314:   */
315:   DAVecRestoreArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,&x);
316:   /* Restrict Xold to coarser levels */
317:   for (mglevel=param->mglevels-1; mglevel>0; mglevel--) {
318:     MatRestrict(dmmg[mglevel]->R, ((AppCtx*)dmmg[mglevel]->user)->Xold, ((AppCtx*)dmmg[mglevel-1]->user)->Xold);
319:     VecPointwiseMult(((AppCtx*)dmmg[mglevel-1]->user)->Xold,((AppCtx*)dmmg[mglevel-1]->user)->Xold,dmmg[mglevel]->Rscale);
320:   }
321: 
322:   /* Store X in the qold for time stepping */
323:   /*VecDuplicate(X,&tsCtx->qold);
324:   VecDuplicate(X,&tsCtx->func);
325:   VecCopy(X,tsCtx->Xold);
326:   tsCtx->ires = 0;
327:   SNESComputeFunction(snes,tsCtx->Xold,tsCtx->func);
328:   VecNorm(tsCtx->func,NORM_2,&tsCtx->fnorm_ini);
329:   tsCtx->ires = 1;
330:   PetscPrintf(PETSC_COMM_WORLD,"Initial function norm is %G\n",tsCtx->fnorm_ini);*/
331:   return 0;
332: }
333: /* ------------------------------------------------------------------- */
336: /* 
337:    FormInitialGuess - Forms initial approximation.

339:    Input Parameters:
340:    user - user-defined application context
341:    X - vector

343:    Output Parameter:
344:    X - vector
345:  */
346: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
347: {
348:   AppCtx         *user = (AppCtx*)dmmg->user;
349:   TstepCtx       *tsCtx = user->tsCtx;

352:   VecCopy(user->Xold, X);
353:   /* calculate the residual on fine mesh */
354:   if (user->tsCtx->fnorm_ini == 0.0) {
355:     tsCtx->ires = 0;
356:     SNESComputeFunction(dmmg->snes,user->Xold,user->func);
357:     VecNorm(user->func,NORM_2,&tsCtx->fnorm_ini);
358:     tsCtx->ires = 1;
359:     tsCtx->cfl = tsCtx->cfl_ini;
360:   }
361:   return 0;
362: }

364: /*---------------------------------------------------------------------*/
367: PetscErrorCode AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,void *ptr)
368: /*---------------------------------------------------------------------*/
369: {
370:   AppCtx         *user = (AppCtx*)ptr;
371:   TstepCtx       *tsCtx = user->tsCtx;
372:   DA             da = info->da;
374:   PetscInt       i,j, xints,xinte,yints,yinte;
375:   PetscReal      hx,hy,dhx,dhy,hxhy;
376:   PassiveScalar  dtinv;
377:   PassiveField   **xold;
378:   PetscTruth     use_parab = tsCtx->use_parabolic;

381:   xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
382:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
383:   hx = 1.0/dhx;                   hy = 1.0/dhy;
384:   hxhy = hx*hy;
385:   DAVecGetArray(da,user->Xold,&xold);
386:   dtinv = hxhy/(tsCtx->cfl*tsCtx->dt);
387:   /* 
388:      use_parab = PETSC_TRUE for parabolic equations; all the four equations have temporal term.
389:                = PETSC_FALSE for differential algebraic equtions (DAE); 
390:                  velocity equations do not have temporal term.
391:   */
392:   for (j=yints; j<yinte; j++) {
393:     for (i=xints; i<xinte; i++) {
394:       if (use_parab) {
395:         f[j][i].u     += dtinv*(x[j][i].u-xold[j][i].u);
396:         f[j][i].v     += dtinv*(x[j][i].v-xold[j][i].v);
397:       }
398:       f[j][i].omega += dtinv*(x[j][i].omega-xold[j][i].omega);
399:       f[j][i].temp  += dtinv*(x[j][i].temp-xold[j][i].temp);
400:     }
401:   }
402:   DAVecRestoreArray(da,user->Xold,&xold);
403:   return(0);
404: }

408: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
409:  {
410:   AppCtx         *user = (AppCtx*)ptr;
411:   TstepCtx       *tsCtx = user->tsCtx;
413:   PetscInt       xints,xinte,yints,yinte,i,j;
414:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
415:   PetscReal      grashof,prandtl,lid;
416:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

419:   grashof = user->param->grashof;
420:   prandtl = user->param->prandtl;
421:   lid     = user->param->lidvelocity;

423:   /* 
424:      Define mesh intervals ratios for uniform grid.
425:      [Note: FD formulae below are normalized by multiplying through by
426:      local volume element to obtain coefficients O(1) in two dimensions.]
427:   */
428:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
429:   hx = 1.0/dhx;                   hy = 1.0/dhy;
430:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

432:   xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;

434:   /* Test whether we are on the bottom edge of the global array */
435:   if (yints == 0) {
436:     j = 0;
437:     yints = yints + 1;
438:     /* bottom edge */
439:     for (i=info->xs; i<info->xs+info->xm; i++) {
440:         f[j][i].u     = x[j][i].u;
441:         f[j][i].v     = x[j][i].v;
442:         f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
443:         f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
444:     }
445:   }

447:   /* Test whether we are on the top edge of the global array */
448:   if (yinte == info->my) {
449:     j = info->my - 1;
450:     yinte = yinte - 1;
451:     /* top edge */
452:     for (i=info->xs; i<info->xs+info->xm; i++) {
453:         f[j][i].u     = x[j][i].u - lid;
454:         f[j][i].v     = x[j][i].v;
455:         f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
456:         f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
457:     }
458:   }

460:   /* Test whether we are on the left edge of the global array */
461:   if (xints == 0) {
462:     i = 0;
463:     xints = xints + 1;
464:     /* left edge */
465:     for (j=info->ys; j<info->ys+info->ym; j++) {
466:       f[j][i].u     = x[j][i].u;
467:       f[j][i].v     = x[j][i].v;
468:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
469:       f[j][i].temp  = x[j][i].temp;
470:     }
471:   }

473:   /* Test whether we are on the right edge of the global array */
474:   if (xinte == info->mx) {
475:     i = info->mx - 1;
476:     xinte = xinte - 1;
477:     /* right edge */
478:     for (j=info->ys; j<info->ys+info->ym; j++) {
479:       f[j][i].u     = x[j][i].u;
480:       f[j][i].v     = x[j][i].v;
481:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
482:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
483:     }
484:   }

486:   /* Compute over the interior points */
487:   for (j=yints; j<yinte; j++) {
488:     for (i=xints; i<xinte; i++) {

490:         /*
491:           convective coefficients for upwinding
492:         */
493:         vx = x[j][i].u; avx = PetscAbsScalar(vx);
494:         vxp = .5*(vx+avx); vxm = .5*(vx-avx);
495:         vy = x[j][i].v; avy = PetscAbsScalar(vy);
496:         vyp = .5*(vy+avy); vym = .5*(vy-avy);

498:         /* U velocity */
499:         u          = x[j][i].u;
500:         uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
501:         uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
502:         f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

504:         /* V velocity */
505:         u          = x[j][i].v;
506:         uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
507:         uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
508:         f[j][i].v  = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

510:         /* Omega */
511:         u          = x[j][i].omega;
512:         uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
513:         uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
514:         f[j][i].omega = uxx + uyy +
515:                         (vxp*(u - x[j][i-1].omega) +
516:                           vxm*(x[j][i+1].omega - u)) * hy +
517:                         (vyp*(u - x[j-1][i].omega) +
518:                           vym*(x[j+1][i].omega - u)) * hx -
519:                         .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

521:         /* Temperature */
522:         u             = x[j][i].temp;
523:         uxx           = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
524:         uyy           = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
525:         f[j][i].temp =  uxx + uyy  + prandtl * (
526:                         (vxp*(u - x[j][i-1].temp) +
527:                           vxm*(x[j][i+1].temp - u)) * hy +
528:                         (vyp*(u - x[j-1][i].temp) +
529:                                  vym*(x[j+1][i].temp - u)) * hx);
530:     }
531:   }

533:   /* Add time step contribution */
534:   if (tsCtx->ires) {
535:     AddTSTermLocal(info,x,f,ptr);
536:   }
537:   /*
538:      Flop count (multiply-adds are counted as 2 operations)
539:   */
540:   PetscLogFlops(84.0*info->ym*info->xm);
541:   return(0);
542: }

546: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
547:  {
548:   AppCtx      *user = (AppCtx*)ptr;
549:   PetscInt    i,j,c;
550:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
551:   PassiveReal grashof,prandtl,lid;
552:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

555:   grashof = user->param->grashof;
556:   prandtl = user->param->prandtl;
557:   lid     = user->param->lidvelocity;

559:   /* 
560:      Define mesh intervals ratios for uniform grid.
561:      [Note: FD formulae below are normalized by multiplying through by
562:      local volume element to obtain coefficients O(1) in two dimensions.]
563:   */
564:   dhx = (PetscReal)(info->mx-1);     dhy = (PetscReal)(info->my-1);
565:   hx = 1.0/dhx;                   hy = 1.0/dhy;
566:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

568:   i = st->i; j = st->j; c = st->c;

570:   /* Test whether we are on the right edge of the global array */
571:   if (i == info->mx-1) {
572:     if (c == 0) *f     = x[j][i].u;
573:     else if (c == 1) *f     = x[j][i].v;
574:     else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
575:     else *f  = x[j][i].temp - (PetscReal)(grashof>0);

577:   /* Test whether we are on the left edge of the global array */
578:   } else if (i == 0) {
579:     if (c == 0) *f     = x[j][i].u;
580:     else if (c == 1) *f     = x[j][i].v;
581:     else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
582:     else *f  = x[j][i].temp;

584:   /* Test whether we are on the top edge of the global array */
585:   } else if (j == info->my-1) {
586:     if (c == 0) *f     = x[j][i].u - lid;
587:     else if (c == 1) *f     = x[j][i].v;
588:     else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
589:     else *f  = x[j][i].temp-x[j-1][i].temp;

591:   /* Test whether we are on the bottom edge of the global array */
592:   } else if (j == 0) {
593:     if (c == 0) *f     = x[j][i].u;
594:     else if (c == 1) *f     = x[j][i].v;
595:     else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
596:     else *f  = x[j][i].temp-x[j+1][i].temp;

598:   /* Compute over the interior points */
599:   } else {
600:     /*
601:       convective coefficients for upwinding
602:     */
603:     vx = x[j][i].u; avx = PetscAbsScalar(vx);
604:     vxp = .5*(vx+avx); vxm = .5*(vx-avx);
605:     vy = x[j][i].v; avy = PetscAbsScalar(vy);
606:     vyp = .5*(vy+avy); vym = .5*(vy-avy);

608:     /* U velocity */
609:     if (c == 0) {
610:       u          = x[j][i].u;
611:       uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
612:       uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
613:       *f         = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

615:     /* V velocity */
616:     } else if (c == 1) {
617:       u          = x[j][i].v;
618:       uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
619:       uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
620:       *f         = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
621: 
622:     /* Omega */
623:     } else if (c == 2) {
624:       u          = x[j][i].omega;
625:       uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
626:       uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
627:       *f         = uxx + uyy +
628:         (vxp*(u - x[j][i-1].omega) +
629:          vxm*(x[j][i+1].omega - u)) * hy +
630:         (vyp*(u - x[j-1][i].omega) +
631:          vym*(x[j+1][i].omega - u)) * hx -
632:         .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
633: 
634:     /* Temperature */
635:     } else {
636:       u           = x[j][i].temp;
637:       uxx         = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
638:       uyy         = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
639:       *f          =  uxx + uyy  + prandtl * (
640:                                              (vxp*(u - x[j][i-1].temp) +
641:                                               vxm*(x[j][i+1].temp - u)) * hy +
642:                                              (vyp*(u - x[j-1][i].temp) +
643:                                               vym*(x[j+1][i].temp - u)) * hx);
644:     }
645:   }

647:   return(0);
648: }


651: /*---------------------------------------------------------------------*/
654: PetscErrorCode Update(DMMG *dmmg)
655: /*---------------------------------------------------------------------*/
656: {
657:  AppCtx         *user = (AppCtx *) ((dmmg[0])->user);
658:  TstepCtx         *tsCtx = user->tsCtx;
659:  Parameter      *param = user->param;
660:  SNES           snes;
662:  PetscScalar         fratio;
663:  PetscScalar         time1,time2,cpuloc;
664:  PetscInt       max_steps,its;
665:  PetscTruth     print_flag = PETSC_FALSE;
666:  PetscInt       nfailsCum = 0,nfails = 0;

669:   /* Note: print_flag displays diagnostic info, not convergence behavior. Use '-snes_monitor' for converges info. */
670:   PetscOptionsHasName(PETSC_NULL,"-print",&print_flag);
671:   if (user->param->PreLoading)
672:    max_steps = 1;
673:   else
674:    max_steps = tsCtx->max_steps;
675:   fratio = 1.0;
676: 
677:   PetscGetTime(&time1);
678:   for (tsCtx->itstep = 0; (tsCtx->itstep < max_steps) &&
679:          (fratio <= tsCtx->fnorm_ratio); tsCtx->itstep++) {
680:     DMMGSolve(dmmg);
681:     snes = DMMGGetSNES(dmmg);
682:     SNESGetIterationNumber(snes,&its);
683:     PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n", its);
684:     SNESGetNonlinearStepFailures(snes,&nfails);
685:     nfailsCum += nfails; nfails = 0;
686:     if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
687:     /*tsCtx->qcur = DMMGGetx(dmmg);
688:       VecCopy(tsCtx->qcur,tsCtx->qold);*/
689: 
690:     VecCopy(dmmg[param->mglevels-1]->x, ((AppCtx*)dmmg[param->mglevels-1]->user)->Xold);
691:     for (its=param->mglevels-1; its>0 ;its--) {
692:       MatRestrict(dmmg[its]->R, ((AppCtx*)dmmg[its]->user)->Xold, ((AppCtx*)dmmg[its-1]->user)->Xold);
693:       VecPointwiseMult(((AppCtx*)dmmg[its-1]->user)->Xold,((AppCtx*)dmmg[its-1]->user)->Xold,dmmg[its]->Rscale);
694:     }

696: 
697:     ComputeTimeStep(snes,((AppCtx*)dmmg[param->mglevels-1]->user));
698:     if (print_flag) {
699:       PetscPrintf(PETSC_COMM_WORLD,"At Time Step %D cfl = %G and fnorm = %G\n",
700:                          tsCtx->itstep,tsCtx->cfl,tsCtx->fnorm);
701:       PetscPrintf(PETSC_COMM_WORLD,"Wall clock time needed %G seconds for %D time steps\n",
702:                          cpuloc,tsCtx->itstep);
703:     }
704:     fratio = tsCtx->fnorm_ini/tsCtx->fnorm;
705:     PetscGetTime(&time2);
706:     cpuloc = time2-time1;
707:     MPI_Barrier(PETSC_COMM_WORLD);
708:   } /* End of time step loop */

710:   if (print_flag) {
711:     PetscPrintf(PETSC_COMM_WORLD,"Total wall clock time needed %G seconds for %D time steps\n",
712:                      cpuloc,tsCtx->itstep);
713:     PetscPrintf(PETSC_COMM_WORLD,"cfl = %G fnorm = %G\n",tsCtx->cfl,tsCtx->fnorm);
714:   }
715:   if (user->param->PreLoading) {
716:     tsCtx->fnorm_ini = 0.0;
717:     PetscPrintf(PETSC_COMM_WORLD,"Preloading done ...\n");
718:   }
719:   /*
720:   {
721:     Vec xx,yy;
722:     PetscScalar fnorm,fnorm1;
723:     SNESGetFunctionNorm(snes,&fnorm);
724:     xx = DMMGGetx(dmmg);
725:     VecDuplicate(xx,&yy);
726:     SNESComputeFunction(snes,xx,yy);
727:     VecNorm(yy,NORM_2,&fnorm1);
728:     PetscPrintf(PETSC_COMM_WORLD,"fnorm = %G, fnorm1 = %G\n",fnorm,fnorm1);
729:     
730:   }
731:   */

733:   return(0);
734: }

736: /*---------------------------------------------------------------------*/
739: PetscErrorCode ComputeTimeStep(SNES snes,void *ptr)
740: /*---------------------------------------------------------------------*/
741: {
742:   AppCtx         *user = (AppCtx*)ptr;
743:   TstepCtx       *tsCtx = user->tsCtx;
744:   Vec                 func = user->func;
745:   PetscScalar    inc = 1.1,  newcfl;
747:   /*int               iramp = tsCtx->iramp;*/
748: 
750:   tsCtx->dt = 0.01;
751:   PetscOptionsGetScalar(PETSC_NULL,"-deltat",&tsCtx->dt,PETSC_NULL);
752:   tsCtx->ires = 0;
753:   SNESComputeFunction(snes,user->Xold,user->func);
754:   tsCtx->ires = 1;
755:   VecNorm(func,NORM_2,&tsCtx->fnorm);
756:   newcfl     = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
757:   tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
758:   /* first time through so compute initial function norm */
759:   /*if (tsCtx->fnorm_ini == 0.0) {
760:     tsCtx->fnorm_ini = tsCtx->fnorm;
761:     tsCtx->cfl       = tsCtx->cfl_ini;
762:   } else {
763:     newcfl     = inc*tsCtx->cfl_ini*tsCtx->fnorm_ini/tsCtx->fnorm;
764:     tsCtx->cfl = PetscMin(newcfl,tsCtx->cfl_max);
765:     }*/
766:   return(0);
767: }