Actual source code: ex19.c

      2: static char help[] = "Nonlinear driven cavity with multigrid in 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:         - Lap(U) - Grad_y(Omega) = 0
     26:         - Lap(V) + Grad_x(Omega) = 0
     27:         - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
     28:         - 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: */
     71: typedef struct {
     72:   PetscScalar u,v,omega,temp;
     73: } Field;

     79: typedef struct {
     80:    PassiveReal  lidvelocity,prandtl,grashof;  /* physical parameters */
     81:    PetscTruth     draw_contours;                /* flag - 1 indicates drawing contours */
     82: } AppCtx;

     86: int main(int argc,char **argv)
     87: {
     88:   DMMG           *dmmg;               /* multilevel grid structure */
     89:   AppCtx         user;                /* user-defined work context */
     90:   PetscInt       mx,my,its;
     92:   MPI_Comm       comm;
     93:   SNES           snes;
     94:   DA             da;

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

    100:   PreLoadBegin(PETSC_TRUE,"SetUp");
    101:     DMMGCreate(comm,2,&user,&dmmg);

    104:     /*
    105:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
    106:       for principal unknowns (x) and governing residuals (f)
    107:     */
    108:     DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
    109:     DMMGSetDM(dmmg,(DM)da);
    110:     DADestroy(da);

    112:     DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
    113:                      PETSC_IGNORE,PETSC_IGNORE);
    114:     /*
    115:      Problem parameters (velocity of lid, prandtl, and grashof numbers)
    116:     */
    117:     user.lidvelocity = 1.0/(mx*my);
    118:     user.prandtl     = 1.0;
    119:     user.grashof     = 1.0;
    120:     PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
    121:     PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
    122:     PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);
    123:     PetscOptionsHasName(PETSC_NULL,"-contours",&user.draw_contours);

    125:     DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
    126:     DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
    127:     DASetFieldName(DMMGGetDA(dmmg),2,"Omega");
    128:     DASetFieldName(DMMGGetDA(dmmg),3,"temperature");

    130:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    131:        Create user context, set problem data, create vector data structures.
    132:        Also, compute the initial guess.
    133:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    135:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    136:        Create nonlinear solver context

    138:        Process adiC(36): FormFunctionLocal FormFunctionLocali
    139:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    140:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
    141:     DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);

    143:     PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g\n",
    144:                        user.lidvelocity,user.prandtl,user.grashof);

    147:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    148:        Solve the nonlinear system
    149:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    150:     DMMGSetInitialGuess(dmmg,FormInitialGuess);

    152:   PreLoadStage("Solve");
    153:     DMMGSolve(dmmg);

    155:     snes = DMMGGetSNES(dmmg);
    156:     SNESGetIterationNumber(snes,&its);
    157:     PetscPrintf(comm,"Number of Newton iterations = %D\n", its);

    159:     /*
    160:       Visualize solution
    161:     */

    163:     if (user.draw_contours) {
    164:       VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
    165:     }

    167:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    168:        Free work space.  All PETSc objects should be destroyed when they
    169:        are no longer needed.
    170:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    172:     DMMGDestroy(dmmg);
    173:   PreLoadEnd();

    175:   PetscFinalize();
    176:   return 0;
    177: }

    179: /* ------------------------------------------------------------------- */

    184: /*
    185:    FormInitialGuess - Forms initial approximation.

    187:    Input Parameters:
    188:    user - user-defined application context
    189:    X - vector

    191:    Output Parameter:
    192:    X - vector
    193:  */
    194: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
    195: {
    196:   AppCtx         *user = (AppCtx*)dmmg->user;
    197:   DA             da = (DA)dmmg->dm;
    198:   PetscInt       i,j,mx,xs,ys,xm,ym;
    200:   PetscReal      grashof,dx;
    201:   Field          **x;

    203:   grashof = user->grashof;

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

    208:   /*
    209:      Get local grid boundaries (for 2-dimensional DA):
    210:        xs, ys   - starting grid indices (no ghost points)
    211:        xm, ym   - widths of local grid (no ghost points)
    212:   */
    213:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

    215:   /*
    216:      Get a pointer to vector data.
    217:        - For default PETSc vectors, VecGetArray() returns a pointer to
    218:          the data array.  Otherwise, the routine is implementation dependent.
    219:        - You MUST call VecRestoreArray() when you no longer need access to
    220:          the array.
    221:   */
    222:   DAVecGetArray(da,X,&x);

    224:   /*
    225:      Compute initial guess over the locally owned part of the grid
    226:      Initial condition is motionless fluid and equilibrium temperature
    227:   */
    228:   for (j=ys; j<ys+ym; j++) {
    229:     for (i=xs; i<xs+xm; i++) {
    230:       x[j][i].u     = 0.0;
    231:       x[j][i].v     = 0.0;
    232:       x[j][i].omega = 0.0;
    233:       x[j][i].temp  = (grashof>0)*i*dx;
    234:     }
    235:   }

    237:   /*
    238:      Restore vector
    239:   */
    240:   DAVecRestoreArray(da,X,&x);
    241:   return 0;
    242: }
    243: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
    244:  {
    245:   AppCtx         *user = (AppCtx*)ptr;
    247:   PetscInt       xints,xinte,yints,yinte,i,j;
    248:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
    249:   PetscReal      grashof,prandtl,lid;
    250:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

    253:   grashof = user->grashof;
    254:   prandtl = user->prandtl;
    255:   lid     = user->lidvelocity;

    257:   /*
    258:      Define mesh intervals ratios for uniform grid.
    259:      [Note: FD formulae below are normalized by multiplying through by
    260:      local volume element to obtain coefficients O(1) in two dimensions.]
    261:   */
    262:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
    263:   hx = 1.0/dhx;                   hy = 1.0/dhy;
    264:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

    268:   /* Test whether we are on the bottom edge of the global array */
    269:   if (yints == 0) {
    270:     j = 0;
    271:     yints = yints + 1;
    272:     /* bottom edge */
    273:     for (i=info->xs; i<info->xs+info->xm; i++) {
    274:         f[j][i].u     = x[j][i].u;
    275:         f[j][i].v     = x[j][i].v;
    276:         f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
    277:         f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
    278:     }
    279:   }

    281:   /* Test whether we are on the top edge of the global array */
    282:   if (yinte == info->my) {
    283:     j = info->my - 1;
    284:     yinte = yinte - 1;
    285:     /* top edge */
    286:     for (i=info->xs; i<info->xs+info->xm; i++) {
    287:         f[j][i].u     = x[j][i].u - lid;
    288:         f[j][i].v     = x[j][i].v;
    289:         f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
    290:         f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
    291:     }
    292:   }

    294:   /* Test whether we are on the left edge of the global array */
    295:   if (xints == 0) {
    296:     i = 0;
    297:     xints = xints + 1;
    298:     /* left edge */
    299:     for (j=info->ys; j<info->ys+info->ym; j++) {
    300:       f[j][i].u     = x[j][i].u;
    301:       f[j][i].v     = x[j][i].v;
    302:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
    303:       f[j][i].temp  = x[j][i].temp;
    304:     }
    305:   }

    307:   /* Test whether we are on the right edge of the global array */
    308:   if (xinte == info->mx) {
    309:     i = info->mx - 1;
    310:     xinte = xinte - 1;
    311:     /* right edge */
    312:     for (j=info->ys; j<info->ys+info->ym; j++) {
    313:       f[j][i].u     = x[j][i].u;
    314:       f[j][i].v     = x[j][i].v;
    315:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
    316:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
    317:     }
    318:   }

    320:   /* Compute over the interior points */
    321:   for (j=yints; j<yinte; j++) {
    322:     for (i=xints; i<xinte; i++) {

    324:         /*
    325:           convective coefficients for upwinding
    326:         */
    327:         vx = x[j][i].u; avx = PetscAbsScalar(vx);
    328:         vxp = .5*(vx+avx); vxm = .5*(vx-avx);
    329:         vy = x[j][i].v; avy = PetscAbsScalar(vy);
    330:         vyp = .5*(vy+avy); vym = .5*(vy-avy);

    332:         /* U velocity */
    333:         u          = x[j][i].u;
    334:         uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
    335:         uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
    336:         f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

    338:         /* V velocity */
    339:         u          = x[j][i].v;
    340:         uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
    341:         uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
    342:         f[j][i].v  = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

    344:         /* Omega */
    345:         u          = x[j][i].omega;
    346:         uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
    347:         uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
    348:         f[j][i].omega = uxx + uyy +
    349:                         (vxp*(u - x[j][i-1].omega) +
    350:                           vxm*(x[j][i+1].omega - u)) * hy +
    351:                         (vyp*(u - x[j-1][i].omega) +
    352:                           vym*(x[j+1][i].omega - u)) * hx -
    353:                         .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

    355:         /* Temperature */
    356:         u             = x[j][i].temp;
    357:         uxx           = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
    358:         uyy           = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
    359:         f[j][i].temp =  uxx + uyy  + prandtl * (
    360:                         (vxp*(u - x[j][i-1].temp) +
    361:                           vxm*(x[j][i+1].temp - u)) * hy +
    362:                         (vyp*(u - x[j-1][i].temp) +
    363:                                  vym*(x[j+1][i].temp - u)) * hx);
    364:     }
    365:   }

    367:   /*
    368:      Flop count (multiply-adds are counted as 2 operations)
    369:   */
    370:   PetscLogFlops(84*info->ym*info->xm);
    371:   return(0);
    372: }

    374: /*
    375:     This is an experimental function and can be safely ignored.
    376: */
    377: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
    378:  {
    379:   AppCtx      *user = (AppCtx*)ptr;
    380:   PetscInt    i,j,c;
    381:   PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
    382:   PassiveReal grashof,prandtl,lid;
    383:   PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

    386:   grashof = user->grashof;
    387:   prandtl = user->prandtl;
    388:   lid     = user->lidvelocity;

    390:   /*
    391:      Define mesh intervals ratios for uniform grid.
    392:      [Note: FD formulae below are normalized by multiplying through by
    393:      local volume element to obtain coefficients O(1) in two dimensions.]
    394:   */
    395:   dhx = (PetscReal)(info->mx-1);     dhy = (PetscReal)(info->my-1);
    396:   hx = 1.0/dhx;                   hy = 1.0/dhy;
    397:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

    401:   /* Test whether we are on the right edge of the global array */
    402:   if (i == info->mx-1) {
    403:     if (c == 0)      *f = x[j][i].u;
    404:     else if (c == 1) *f = x[j][i].v;
    405:     else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
    406:     else             *f = x[j][i].temp - (PetscReal)(grashof>0);

    408:   /* Test whether we are on the left edge of the global array */
    409:   } else if (i == 0) {
    410:     if (c == 0)      *f = x[j][i].u;
    411:     else if (c == 1) *f = x[j][i].v;
    412:     else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
    413:     else             *f = x[j][i].temp;

    415:   /* Test whether we are on the top edge of the global array */
    416:   } else if (j == info->my-1) {
    417:     if (c == 0)      *f = x[j][i].u - lid;
    418:     else if (c == 1) *f = x[j][i].v;
    419:     else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
    420:     else             *f = x[j][i].temp-x[j-1][i].temp;

    422:   /* Test whether we are on the bottom edge of the global array */
    423:   } else if (j == 0) {
    424:     if (c == 0)      *f = x[j][i].u;
    425:     else if (c == 1) *f = x[j][i].v;
    426:     else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
    427:     else             *f = x[j][i].temp - x[j+1][i].temp;

    429:   /* Compute over the interior points */
    430:   } else {
    431:     /*
    432:       convective coefficients for upwinding
    433:     */
    434:     vx = x[j][i].u; avx = PetscAbsScalar(vx);
    435:     vxp = .5*(vx+avx); vxm = .5*(vx-avx);
    436:     vy = x[j][i].v; avy = PetscAbsScalar(vy);
    437:     vyp = .5*(vy+avy); vym = .5*(vy-avy);

    439:     /* U velocity */
    440:     if (c == 0) {
    441:       u          = x[j][i].u;
    442:       uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
    443:       uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
    444:       *f         = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

    446:     /* V velocity */
    447:     } else if (c == 1) {
    448:       u          = x[j][i].v;
    449:       uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
    450:       uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
    451:       *f         = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
    452:
    453:     /* Omega */
    454:     } else if (c == 2) {
    455:       u          = x[j][i].omega;
    456:       uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
    457:       uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
    458:       *f         = uxx + uyy +
    459:         (vxp*(u - x[j][i-1].omega) +
    460:          vxm*(x[j][i+1].omega - u)) * hy +
    461:         (vyp*(u - x[j-1][i].omega) +
    462:          vym*(x[j+1][i].omega - u)) * hx -
    463:         .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
    464:
    465:     /* Temperature */
    466:     } else {
    467:       u           = x[j][i].temp;
    468:       uxx         = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
    469:       uyy         = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
    470:       *f          =  uxx + uyy  + prandtl * (
    471:                                              (vxp*(u - x[j][i-1].temp) +
    472:                                               vxm*(x[j][i+1].temp - u)) * hy +
    473:                                              (vyp*(u - x[j-1][i].temp) +
    474:                                               vym*(x[j+1][i].temp - u)) * hx);
    475:     }
    476:   }

    478:   return(0);
    479: }