Actual source code: ex26.c

      1: static char help[] = "Transient nonlinear driven cavity in 2d.\n\
      2:   \n\
      3: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
      4: The flow can be driven with the lid or with buoyancy or both:\n\
      5:   -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
      6:   -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
      7:   -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
      8:   -contours : draw contour plots of solution\n\n";
      9: /*
     10:       See src/snes/tutorials/ex19.c for the steady-state version.

     12:       There used to be a SNES example (src/snes/tutorials/ex27.c) that
     13:       implemented this algorithm without using TS and was used for the numerical
     14:       results in the paper

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

     19:       That example was removed because it used obsolete interfaces, but the
     20:       algorithms from the paper can be reproduced using this example.

     22:       Note: The paper describes the algorithm as being linearly implicit but the
     23:       numerical results were created using nonlinearly implicit Euler.  The
     24:       algorithm as described (linearly implicit) is more efficient and is the
     25:       default when using TSPSEUDO.  If you want to reproduce the numerical
     26:       results from the paper, you'll have to change the SNES to converge the
     27:       nonlinear solve (e.g., -snes_type newtonls).  The DAE versus ODE variants
     28:       are controlled using the -parabolic option.

     30:       Comment preserved from snes/tutorials/ex27.c, since removed:

     32:         [H]owever Figure 3.1 was generated with a slightly different algorithm
     33:         (see targets runex27 and runex27_p) than described in the paper.  In
     34:         particular, the described algorithm is linearly implicit, advancing to
     35:         the next step after one Newton step, so that the steady state residual
     36:         is always used, but the figure was generated by converging each step to
     37:         a relative tolerance of 1.e-3.  On the example problem, setting
     38:         -snes_type ksponly has only minor impact on number of steps, but
     39:         significantly reduces the required number of linear solves.

     41:       See also https://lists.mcs.anl.gov/pipermail/petsc-dev/2010-March/002362.html
     42: */

     44: /* ------------------------------------------------------------------------

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

     49:     This problem is modeled by the partial differential equation system

     51:         - Lap(U) - Grad_y(Omega) = 0
     52:         - Lap(V) + Grad_x(Omega) = 0
     53:         Omega_t - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
     54:         T_t - Lap(T) + PR*Div([U*T,V*T]) = 0

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

     59:     No-slip, rigid-wall Dirichlet conditions are used for [U,V].
     60:     Dirichlet conditions are used for Omega, based on the definition of
     61:     vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
     62:     constant coordinate boundary, the tangential derivative is zero.
     63:     Dirichlet conditions are used for T on the left and right walls,
     64:     and insulation homogeneous Neumann conditions are used for T on
     65:     the top and bottom walls.

     67:     A finite difference approximation with the usual 5-point stencil
     68:     is used to discretize the boundary value problem to obtain a
     69:     nonlinear system of equations.  Upwinding is used for the divergence
     70:     (convective) terms and central for the gradient (source) terms.

     72:     The Jacobian can be either
     73:       * formed via finite differencing using coloring (the default), or
     74:       * applied matrix-free via the option -snes_mf
     75:         (for larger grid problems this variant may not converge
     76:         without a preconditioner due to ill-conditioning).

     78:   ------------------------------------------------------------------------- */

     80: /*
     81:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
     82:    Include "petscts.h" so that we can use TS solvers.  Note that this
     83:    file automatically includes:
     84:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     85:      petscmat.h - matrices
     86:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     87:      petscviewer.h - viewers               petscpc.h  - preconditioners
     88:      petscksp.h   - linear solvers         petscsnes.h - nonlinear solvers
     89: */
     90: #include <petscts.h>
     91: #include <petscdm.h>
     92: #include <petscdmda.h>

     94: /*
     95:    User-defined routines and data structures
     96: */
     97: typedef struct {
     98:   PetscScalar u, v, omega, temp;
     99: } Field;

    101: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, void *);

    103: typedef struct {
    104:   PetscReal lidvelocity, prandtl, grashof; /* physical parameters */
    105:   PetscBool parabolic;                     /* allow a transient term corresponding roughly to artificial compressibility */
    106:   PetscReal cfl_initial;                   /* CFL for first time step */
    107: } AppCtx;

    109: PetscErrorCode FormInitialSolution(TS, Vec, AppCtx *);

    111: int main(int argc, char **argv)
    112: {
    113:   AppCtx            user; /* user-defined work context */
    114:   PetscInt          mx, my, steps;
    115:   TS                ts;
    116:   DM                da;
    117:   Vec               X;
    118:   PetscReal         ftime;
    119:   TSConvergedReason reason;

    121:   PetscFunctionBeginUser;
    122:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    123:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    124:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da));
    125:   PetscCall(DMSetFromOptions(da));
    126:   PetscCall(DMSetUp(da));
    127:   PetscCall(TSSetDM(ts, (DM)da));

    129:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
    130:   /*
    131:      Problem parameters (velocity of lid, prandtl, and grashof numbers)
    132:   */
    133:   user.lidvelocity = 1.0 / (mx * my);
    134:   user.prandtl     = 1.0;
    135:   user.grashof     = 1.0;
    136:   user.parabolic   = PETSC_FALSE;
    137:   user.cfl_initial = 50.;

    139:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Driven cavity/natural convection options", "");
    140:   PetscCall(PetscOptionsReal("-lidvelocity", "Lid velocity, related to Reynolds number", "", user.lidvelocity, &user.lidvelocity, NULL));
    141:   PetscCall(PetscOptionsReal("-prandtl", "Ratio of viscous to thermal diffusivity", "", user.prandtl, &user.prandtl, NULL));
    142:   PetscCall(PetscOptionsReal("-grashof", "Ratio of buoyant to viscous forces", "", user.grashof, &user.grashof, NULL));
    143:   PetscCall(PetscOptionsBool("-parabolic", "Relax incompressibility to make the system parabolic instead of differential-algebraic", "", user.parabolic, &user.parabolic, NULL));
    144:   PetscCall(PetscOptionsReal("-cfl_initial", "Advective CFL for the first time step", "", user.cfl_initial, &user.cfl_initial, NULL));
    145:   PetscOptionsEnd();

    147:   PetscCall(DMDASetFieldName(da, 0, "x-velocity"));
    148:   PetscCall(DMDASetFieldName(da, 1, "y-velocity"));
    149:   PetscCall(DMDASetFieldName(da, 2, "Omega"));
    150:   PetscCall(DMDASetFieldName(da, 3, "temperature"));

    152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    153:      Create user context, set problem data, create vector data structures.
    154:      Also, compute the initial guess.
    155:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    158:      Create time integration context
    159:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    160:   PetscCall(DMSetApplicationContext(da, &user));
    161:   PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (DMDATSIFunctionLocalFn *)FormIFunctionLocal, &user));
    162:   PetscCall(TSSetMaxSteps(ts, 10000));
    163:   PetscCall(TSSetMaxTime(ts, 1e12));
    164:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    165:   PetscCall(TSSetTimeStep(ts, user.cfl_initial / (user.lidvelocity * mx)));
    166:   PetscCall(TSSetFromOptions(ts));

    168:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "x%" PetscInt_FMT " grid, lid velocity = %g, prandtl # = %g, grashof # = %g\n", mx, my, (double)user.lidvelocity, (double)user.prandtl, (double)user.grashof));

    170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    171:      Solve the nonlinear system
    172:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    174:   PetscCall(DMCreateGlobalVector(da, &X));
    175:   PetscCall(FormInitialSolution(ts, X, &user));

    177:   PetscCall(TSSolve(ts, X));
    178:   PetscCall(TSGetSolveTime(ts, &ftime));
    179:   PetscCall(TSGetStepNumber(ts, &steps));
    180:   PetscCall(TSGetConvergedReason(ts, &reason));

    182:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

    184:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    185:      Free work space.  All PETSc objects should be destroyed when they
    186:      are no longer needed.
    187:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    188:   PetscCall(VecDestroy(&X));
    189:   PetscCall(DMDestroy(&da));
    190:   PetscCall(TSDestroy(&ts));

    192:   PetscCall(PetscFinalize());
    193:   return 0;
    194: }

    196: /* ------------------------------------------------------------------- */

    198: /*
    199:    FormInitialSolution - Forms initial approximation.

    201:    Input Parameters:
    202:    user - user-defined application context
    203:    X - vector

    205:    Output Parameter:
    206:    X - vector
    207:  */
    208: PetscErrorCode FormInitialSolution(TS ts, Vec X, AppCtx *user)
    209: {
    210:   DM        da;
    211:   PetscInt  i, j, mx, xs, ys, xm, ym;
    212:   PetscReal grashof, dx;
    213:   Field   **x;

    215:   PetscFunctionBeginUser;
    216:   grashof = user->grashof;
    217:   PetscCall(TSGetDM(ts, &da));
    218:   PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
    219:   dx = 1.0 / (mx - 1);

    221:   /*
    222:      Get local grid boundaries (for 2-dimensional DMDA):
    223:        xs, ys   - starting grid indices (no ghost points)
    224:        xm, ym   - widths of local grid (no ghost points)
    225:   */
    226:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

    228:   /*
    229:      Get a pointer to vector data.
    230:        - For default PETSc vectors, VecGetArray() returns a pointer to
    231:          the data array.  Otherwise, the routine is implementation dependent.
    232:        - You MUST call VecRestoreArray() when you no longer need access to
    233:          the array.
    234:   */
    235:   PetscCall(DMDAVecGetArray(da, X, &x));

    237:   /*
    238:      Compute initial guess over the locally owned part of the grid
    239:      Initial condition is motionless fluid and equilibrium temperature
    240:   */
    241:   for (j = ys; j < ys + ym; j++) {
    242:     for (i = xs; i < xs + xm; i++) {
    243:       x[j][i].u     = 0.0;
    244:       x[j][i].v     = 0.0;
    245:       x[j][i].omega = 0.0;
    246:       x[j][i].temp  = (grashof > 0) * i * dx;
    247:     }
    248:   }

    250:   /*
    251:      Restore vector
    252:   */
    253:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    254:   PetscFunctionReturn(PETSC_SUCCESS);
    255: }

    257: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal ptime, Field **x, Field **xdot, Field **f, void *ptr)
    258: {
    259:   AppCtx     *user = (AppCtx *)ptr;
    260:   PetscInt    xints, xinte, yints, yinte, i, j;
    261:   PetscReal   hx, hy, dhx, dhy, hxdhy, hydhx;
    262:   PetscReal   grashof, prandtl, lid;
    263:   PetscScalar u, udot, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym;

    265:   PetscFunctionBeginUser;
    266:   grashof = user->grashof;
    267:   prandtl = user->prandtl;
    268:   lid     = user->lidvelocity;

    270:   /*
    271:      Define mesh intervals ratios for uniform grid.

    273:      Note: FD formulae below are normalized by multiplying through by
    274:      local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.

    276:   */
    277:   dhx   = (PetscReal)(info->mx - 1);
    278:   dhy   = (PetscReal)(info->my - 1);
    279:   hx    = 1.0 / dhx;
    280:   hy    = 1.0 / dhy;
    281:   hxdhy = hx * dhy;
    282:   hydhx = hy * dhx;

    284:   xints = info->xs;
    285:   xinte = info->xs + info->xm;
    286:   yints = info->ys;
    287:   yinte = info->ys + info->ym;

    289:   /* Test whether we are on the bottom edge of the global array */
    290:   if (yints == 0) {
    291:     j     = 0;
    292:     yints = yints + 1;
    293:     /* bottom edge */
    294:     for (i = info->xs; i < info->xs + info->xm; i++) {
    295:       f[j][i].u     = x[j][i].u;
    296:       f[j][i].v     = x[j][i].v;
    297:       f[j][i].omega = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy;
    298:       f[j][i].temp  = x[j][i].temp - x[j + 1][i].temp;
    299:     }
    300:   }

    302:   /* Test whether we are on the top edge of the global array */
    303:   if (yinte == info->my) {
    304:     j     = info->my - 1;
    305:     yinte = yinte - 1;
    306:     /* top edge */
    307:     for (i = info->xs; i < info->xs + info->xm; i++) {
    308:       f[j][i].u     = x[j][i].u - lid;
    309:       f[j][i].v     = x[j][i].v;
    310:       f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy;
    311:       f[j][i].temp  = x[j][i].temp - x[j - 1][i].temp;
    312:     }
    313:   }

    315:   /* Test whether we are on the left edge of the global array */
    316:   if (xints == 0) {
    317:     i     = 0;
    318:     xints = xints + 1;
    319:     /* left edge */
    320:     for (j = info->ys; j < info->ys + info->ym; j++) {
    321:       f[j][i].u     = x[j][i].u;
    322:       f[j][i].v     = x[j][i].v;
    323:       f[j][i].omega = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx;
    324:       f[j][i].temp  = x[j][i].temp;
    325:     }
    326:   }

    328:   /* Test whether we are on the right edge of the global array */
    329:   if (xinte == info->mx) {
    330:     i     = info->mx - 1;
    331:     xinte = xinte - 1;
    332:     /* right edge */
    333:     for (j = info->ys; j < info->ys + info->ym; j++) {
    334:       f[j][i].u     = x[j][i].u;
    335:       f[j][i].v     = x[j][i].v;
    336:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx;
    337:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof > 0);
    338:     }
    339:   }

    341:   /* Compute over the interior points */
    342:   for (j = yints; j < yinte; j++) {
    343:     for (i = xints; i < xinte; i++) {
    344:       /*
    345:         convective coefficients for upwinding
    346:       */
    347:       vx  = x[j][i].u;
    348:       avx = PetscAbsScalar(vx);
    349:       vxp = .5 * (vx + avx);
    350:       vxm = .5 * (vx - avx);
    351:       vy  = x[j][i].v;
    352:       avy = PetscAbsScalar(vy);
    353:       vyp = .5 * (vy + avy);
    354:       vym = .5 * (vy - avy);

    356:       /* U velocity */
    357:       u         = x[j][i].u;
    358:       udot      = user->parabolic ? xdot[j][i].u : 0.;
    359:       uxx       = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
    360:       uyy       = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
    361:       f[j][i].u = udot + uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx;

    363:       /* V velocity */
    364:       u         = x[j][i].v;
    365:       udot      = user->parabolic ? xdot[j][i].v : 0.;
    366:       uxx       = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx;
    367:       uyy       = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy;
    368:       f[j][i].v = udot + uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy;

    370:       /* Omega */
    371:       u   = x[j][i].omega;
    372:       uxx = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx;
    373:       uyy = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy;
    374:       f[j][i].omega = (xdot[j][i].omega + uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (vyp * (u - x[j - 1][i].omega) + vym * (x[j + 1][i].omega - u)) * hx - .5 * grashof * (x[j][i + 1].temp - x[j][i - 1].temp) * hy);

    376:       /* Temperature */
    377:       u            = x[j][i].temp;
    378:       uxx          = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx;
    379:       uyy          = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy;
    380:       f[j][i].temp = (xdot[j][i].temp + uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) * hy + (vyp * (u - x[j - 1][i].temp) + vym * (x[j + 1][i].temp - u)) * hx));
    381:     }
    382:   }

    384:   /*
    385:      Flop count (multiply-adds are counted as 2 operations)
    386:   */
    387:   PetscCall(PetscLogFlops(84.0 * info->ym * info->xm));
    388:   PetscFunctionReturn(PETSC_SUCCESS);
    389: }

    391: /*TEST

    393:     test:
    394:       args: -da_grid_x 20 -da_grid_y 20 -lidvelocity 100 -grashof 1e3 -ts_max_steps 100 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type rosw -ts_rosw_type ra3pw -ts_monitor -ts_monitor_solution_vtk 'foo-%03d.vts'
    395:       requires: !complex !single

    397:     test:
    398:       suffix: 2
    399:       nsize: 4
    400:       args: -da_grid_x 20 -da_grid_y 20 -lidvelocity 100 -grashof 1e3 -ts_max_steps 100 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type rosw -ts_rosw_type ra3pw -ts_monitor -ts_monitor_solution_vtk 'foo-%03d.vts' -ts_monitor_solution_vtk_interval -1
    401:       requires: !complex !single

    403:     test:
    404:       suffix: 3
    405:       nsize: 4
    406:       args: -da_refine 2 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 -pc_type none -ts_type beuler -ts_monitor -snes_monitor_short -snes_type aspin -da_overlap 4 -npc_sub_ksp_type preonly -npc_sub_pc_type lu
    407:       requires: !complex !single

    409:     test:
    410:       suffix: 4
    411:       nsize: 2
    412:       args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3
    413:       requires: !complex !single

    415:     test:
    416:       suffix: asm
    417:       nsize: 4
    418:       args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3
    419:       requires: !complex !single

    421: TEST*/