Actual source code: ex6.c

      1: static char help[] = "Solves a simple time-dependent linear PDE (the heat equation).\n\
      2: Input parameters include:\n\
      3:   -m <points>, where <points> = number of grid points\n\
      4:   -time_dependent_rhs : Treat the problem as having a time-dependent right-hand side\n\
      5:   -debug              : Activate debugging printouts\n\
      6:   -nox                : Deactivate x-window graphics\n\n";

      8: /* ------------------------------------------------------------------------

     10:    This program solves the one-dimensional heat equation (also called the
     11:    diffusion equation),
     12:        u_t = u_xx,
     13:    on the domain 0 <= x <= 1, with the boundary conditions
     14:        u(t,0) = 0, u(t,1) = 0,
     15:    and the initial condition
     16:        u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x).
     17:    This is a linear, second-order, parabolic equation.

     19:    We discretize the right-hand side using finite differences with
     20:    uniform grid spacing h:
     21:        u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2)
     22:    We then demonstrate time evolution using the various TS methods by
     23:    running the program via
     24:        ex3 -ts_type <timestepping solver>

     26:    We compare the approximate solution with the exact solution, given by
     27:        u_exact(x,t) = exp(-36*pi*pi*t) * sin(6*pi*x) +
     28:                       3*exp(-4*pi*pi*t) * sin(2*pi*x)

     30:    Notes:
     31:    This code demonstrates the TS solver interface to two variants of
     32:    linear problems, u_t = f(u,t), namely
     33:      - time-dependent f:   f(u,t) is a function of t
     34:      - time-independent f: f(u,t) is simply f(u)

     36:     The parallel version of this code is ts/tutorials/ex4.c

     38:   ------------------------------------------------------------------------- */

     40: /*
     41:    Include "ts.h" so that we can use TS solvers.  Note that this file
     42:    automatically includes:
     43:      petscsys.h  - base PETSc routines   vec.h  - vectors
     44:      sys.h    - system routines       mat.h  - matrices
     45:      is.h     - index sets            ksp.h  - Krylov subspace methods
     46:      viewer.h - viewers               pc.h   - preconditioners
     47:      snes.h - nonlinear solvers
     48: */

     50: #include <petscts.h>
     51: #include <petscdraw.h>

     53: /*
     54:    User-defined application context - contains data needed by the
     55:    application-provided call-back routines.
     56: */
     57: typedef struct {
     58:   Vec         solution;         /* global exact solution vector */
     59:   PetscInt    m;                /* total number of grid points */
     60:   PetscReal   h;                /* mesh width h = 1/(m-1) */
     61:   PetscBool   debug;            /* flag (1 indicates activation of debugging printouts) */
     62:   PetscViewer viewer1, viewer2; /* viewers for the solution and error */
     63:   PetscReal   norm_2, norm_max; /* error norms */
     64: } AppCtx;

     66: /*
     67:    User-defined routines
     68: */
     69: extern PetscErrorCode InitialConditions(Vec, AppCtx *);
     70: extern PetscErrorCode RHSMatrixHeat(TS, PetscReal, Vec, Mat, Mat, void *);
     71: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
     72: extern PetscErrorCode ExactSolution(PetscReal, Vec, AppCtx *);
     73: extern PetscErrorCode MyBCRoutine(TS, PetscReal, Vec, void *);

     75: int main(int argc, char **argv)
     76: {
     77:   AppCtx      appctx;                 /* user-defined application context */
     78:   TS          ts;                     /* timestepping context */
     79:   Mat         A;                      /* matrix data structure */
     80:   Vec         u;                      /* approximate solution vector */
     81:   PetscReal   time_total_max = 100.0; /* default max total time */
     82:   PetscInt    time_steps_max = 100;   /* default max timesteps */
     83:   PetscDraw   draw;                   /* drawing context */
     84:   PetscInt    steps, m;
     85:   PetscMPIInt size;
     86:   PetscReal   dt;
     87:   PetscReal   ftime;
     88:   PetscBool   flg;
     89:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     90:      Initialize program and set problem parameters
     91:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     93:   PetscFunctionBeginUser;
     94:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     95:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
     96:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

     98:   m = 60;
     99:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
    100:   PetscCall(PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug));

    102:   appctx.m        = m;
    103:   appctx.h        = 1.0 / (m - 1.0);
    104:   appctx.norm_2   = 0.0;
    105:   appctx.norm_max = 0.0;

    107:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n"));

    109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    110:      Create vector data structures
    111:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    113:   /*
    114:      Create vector data structures for approximate and exact solutions
    115:   */
    116:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &u));
    117:   PetscCall(VecDuplicate(u, &appctx.solution));

    119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    120:      Set up displays to show graphs of the solution and error
    121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    123:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 380, 400, 160, &appctx.viewer1));
    124:   PetscCall(PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw));
    125:   PetscCall(PetscDrawSetDoubleBuffer(draw));
    126:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 0, 400, 160, &appctx.viewer2));
    127:   PetscCall(PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw));
    128:   PetscCall(PetscDrawSetDoubleBuffer(draw));

    130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    131:      Create timestepping solver context
    132:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    134:   PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
    135:   PetscCall(TSSetProblemType(ts, TS_LINEAR));

    137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    138:      Set optional user-defined monitoring routine
    139:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    141:   PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));

    143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    145:      Create matrix data structure; set matrix evaluation routine.
    146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    148:   PetscCall(MatCreate(PETSC_COMM_SELF, &A));
    149:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
    150:   PetscCall(MatSetFromOptions(A));
    151:   PetscCall(MatSetUp(A));

    153:   PetscCall(PetscOptionsHasName(NULL, NULL, "-time_dependent_rhs", &flg));
    154:   if (flg) {
    155:     /*
    156:        For linear problems with a time-dependent f(u,t) in the equation
    157:        u_t = f(u,t), the user provides the discretized right-hand-side
    158:        as a time-dependent matrix.
    159:     */
    160:     PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
    161:     PetscCall(TSSetRHSJacobian(ts, A, A, RHSMatrixHeat, &appctx));
    162:   } else {
    163:     /*
    164:        For linear problems with a time-independent f(u) in the equation
    165:        u_t = f(u), the user provides the discretized right-hand side
    166:        as a matrix only once, and then sets a null matrix evaluation
    167:        routine.
    168:     */
    169:     PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx));
    170:     PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
    171:     PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &appctx));
    172:   }

    174:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    175:      Set solution vector and initial timestep
    176:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    178:   dt = appctx.h * appctx.h / 2.0;
    179:   PetscCall(TSSetTimeStep(ts, dt));
    180:   PetscCall(TSSetSolution(ts, u));

    182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    183:      Customize timestepping solver:
    184:        - Set the solution method to be the Backward Euler method.
    185:        - Set timestepping duration info
    186:      Then set runtime options, which can override these defaults.
    187:      For example,
    188:           -ts_max_steps <maxsteps> -ts_max_time <maxtime>
    189:      to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
    190:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    192:   PetscCall(TSSetMaxSteps(ts, time_steps_max));
    193:   PetscCall(TSSetMaxTime(ts, time_total_max));
    194:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    195:   PetscCall(TSSetFromOptions(ts));

    197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    198:      Solve the problem
    199:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    201:   /*
    202:      Evaluate initial conditions
    203:   */
    204:   PetscCall(InitialConditions(u, &appctx));

    206:   /*
    207:      Run the timestepping solver
    208:   */
    209:   PetscCall(TSSolve(ts, u));
    210:   PetscCall(TSGetSolveTime(ts, &ftime));
    211:   PetscCall(TSGetStepNumber(ts, &steps));

    213:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    214:      View timestepping solver info
    215:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    217:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "avg. error (2 norm) = %g, avg. error (max norm) = %g\n", (double)(appctx.norm_2 / steps), (double)(appctx.norm_max / steps)));
    218:   PetscCall(TSView(ts, PETSC_VIEWER_STDOUT_SELF));

    220:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    221:      Free work space.  All PETSc objects should be destroyed when they
    222:      are no longer needed.
    223:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    225:   PetscCall(TSDestroy(&ts));
    226:   PetscCall(MatDestroy(&A));
    227:   PetscCall(VecDestroy(&u));
    228:   PetscCall(PetscViewerDestroy(&appctx.viewer1));
    229:   PetscCall(PetscViewerDestroy(&appctx.viewer2));
    230:   PetscCall(VecDestroy(&appctx.solution));

    232:   /*
    233:      Always call PetscFinalize() before exiting a program.  This routine
    234:        - finalizes the PETSc libraries as well as MPI
    235:        - provides summary and diagnostic information if certain runtime
    236:          options are chosen (e.g., -log_view).
    237:   */
    238:   PetscCall(PetscFinalize());
    239:   return 0;
    240: }
    241: /* --------------------------------------------------------------------- */
    242: /*
    243:    InitialConditions - Computes the solution at the initial time.

    245:    Input Parameter:
    246:    u - uninitialized solution vector (global)
    247:    appctx - user-defined application context

    249:    Output Parameter:
    250:    u - vector with solution at initial time (global)
    251: */
    252: PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
    253: {
    254:   PetscScalar *u_localptr;
    255:   PetscInt     i;

    257:   PetscFunctionBeginUser;
    258:   /*
    259:     Get a pointer to vector data.
    260:     - For default PETSc vectors, VecGetArray() returns a pointer to
    261:       the data array.  Otherwise, the routine is implementation dependent.
    262:     - You MUST call VecRestoreArray() when you no longer need access to
    263:       the array.
    264:     - Note that the Fortran interface to VecGetArray() differs from the
    265:       C version.  See the users manual for details.
    266:   */
    267:   PetscCall(VecGetArray(u, &u_localptr));

    269:   /*
    270:      We initialize the solution array by simply writing the solution
    271:      directly into the array locations.  Alternatively, we could use
    272:      VecSetValues() or VecSetValuesLocal().
    273:   */
    274:   for (i = 0; i < appctx->m; i++) u_localptr[i] = PetscSinReal(PETSC_PI * i * 6. * appctx->h) + 3. * PetscSinReal(PETSC_PI * i * 2. * appctx->h);

    276:   /*
    277:      Restore vector
    278:   */
    279:   PetscCall(VecRestoreArray(u, &u_localptr));

    281:   /*
    282:      Print debugging information if desired
    283:   */
    284:   if (appctx->debug) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }
    287: /* --------------------------------------------------------------------- */
    288: /*
    289:    ExactSolution - Computes the exact solution at a given time.

    291:    Input Parameters:
    292:    t - current time
    293:    solution - vector in which exact solution will be computed
    294:    appctx - user-defined application context

    296:    Output Parameter:
    297:    solution - vector with the newly computed exact solution
    298: */
    299: PetscErrorCode ExactSolution(PetscReal t, Vec solution, AppCtx *appctx)
    300: {
    301:   PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2;
    302:   PetscInt     i;

    304:   PetscFunctionBeginUser;
    305:   /*
    306:      Get a pointer to vector data.
    307:   */
    308:   PetscCall(VecGetArray(solution, &s_localptr));

    310:   /*
    311:      Simply write the solution directly into the array locations.
    312:      Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
    313:   */
    314:   ex1 = PetscExpReal(-36. * PETSC_PI * PETSC_PI * t);
    315:   ex2 = PetscExpReal(-4. * PETSC_PI * PETSC_PI * t);
    316:   sc1 = PETSC_PI * 6. * h;
    317:   sc2 = PETSC_PI * 2. * h;
    318:   for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscSinReal(PetscRealPart(sc1) * (PetscReal)i) * ex1 + 3. * PetscSinReal(PetscRealPart(sc2) * (PetscReal)i) * ex2;

    320:   /*
    321:      Restore vector
    322:   */
    323:   PetscCall(VecRestoreArray(solution, &s_localptr));
    324:   PetscFunctionReturn(PETSC_SUCCESS);
    325: }
    326: /* --------------------------------------------------------------------- */
    327: /*
    328:    Monitor - User-provided routine to monitor the solution computed at
    329:    each timestep.  This example plots the solution and computes the
    330:    error in two different norms.

    332:    This example also demonstrates changing the timestep via TSSetTimeStep().

    334:    Input Parameters:
    335:    ts     - the timestep context
    336:    step   - the count of the current step (with 0 meaning the
    337:              initial condition)
    338:    crtime  - the current time
    339:    u      - the solution at this timestep
    340:    ctx    - the user-provided context for this monitoring routine.
    341:             In this case we use the application context which contains
    342:             information about the problem size, workspace and the exact
    343:             solution.
    344: */
    345: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
    346: {
    347:   AppCtx   *appctx = (AppCtx *)ctx; /* user-defined application context */
    348:   PetscReal norm_2, norm_max, dt, dttol;
    349:   PetscBool flg;

    351:   PetscFunctionBeginUser;
    352:   /*
    353:      View a graph of the current iterate
    354:   */
    355:   PetscCall(VecView(u, appctx->viewer2));

    357:   /*
    358:      Compute the exact solution
    359:   */
    360:   PetscCall(ExactSolution(crtime, appctx->solution, appctx));

    362:   /*
    363:      Print debugging information if desired
    364:   */
    365:   if (appctx->debug) {
    366:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n"));
    367:     PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
    368:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
    369:     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
    370:   }

    372:   /*
    373:      Compute the 2-norm and max-norm of the error
    374:   */
    375:   PetscCall(VecAXPY(appctx->solution, -1.0, u));
    376:   PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
    377:   norm_2 = PetscSqrtReal(appctx->h) * norm_2;
    378:   PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));

    380:   PetscCall(TSGetTimeStep(ts, &dt));
    381:   if (norm_2 > 1.e-2) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Timestep %" PetscInt_FMT ": step size = %g, time = %g, 2-norm error = %g, max norm error = %g\n", step, (double)dt, (double)crtime, (double)norm_2, (double)norm_max));
    382:   appctx->norm_2 += norm_2;
    383:   appctx->norm_max += norm_max;

    385:   dttol = .0001;
    386:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, &flg));
    387:   if (dt < dttol) {
    388:     dt *= .999;
    389:     PetscCall(TSSetTimeStep(ts, dt));
    390:   }

    392:   /*
    393:      View a graph of the error
    394:   */
    395:   PetscCall(VecView(appctx->solution, appctx->viewer1));

    397:   /*
    398:      Print debugging information if desired
    399:   */
    400:   if (appctx->debug) {
    401:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
    402:     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
    403:   }
    404:   PetscFunctionReturn(PETSC_SUCCESS);
    405: }
    406: /* --------------------------------------------------------------------- */
    407: /*
    408:    RHSMatrixHeat - User-provided routine to compute the right-hand-side
    409:    matrix for the heat equation.

    411:    Input Parameters:
    412:    ts - the TS context
    413:    t - current time
    414:    global_in - global input vector
    415:    dummy - optional user-defined context, as set by TSetRHSJacobian()

    417:    Output Parameters:
    418:    AA - Jacobian matrix
    419:    BB - optionally different matrix used to construct the preconditioner

    421:    Notes:
    422:    Recall that MatSetValues() uses 0-based row and column numbers
    423:    in Fortran as well as in C.
    424: */
    425: PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)
    426: {
    427:   Mat         A      = AA;            /* Jacobian matrix */
    428:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
    429:   PetscInt    mstart = 0;
    430:   PetscInt    mend   = appctx->m;
    431:   PetscInt    i, idx[3];
    432:   PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;

    434:   PetscFunctionBeginUser;
    435:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    436:      Compute entries for the locally owned part of the matrix
    437:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    438:   /*
    439:      Set matrix rows corresponding to boundary data
    440:   */

    442:   mstart = 0;
    443:   v[0]   = 1.0;
    444:   PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
    445:   mstart++;

    447:   mend--;
    448:   v[0] = 1.0;
    449:   PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));

    451:   /*
    452:      Set matrix rows corresponding to interior data.  We construct the
    453:      matrix one row at a time.
    454:   */
    455:   v[0] = sone;
    456:   v[1] = stwo;
    457:   v[2] = sone;
    458:   for (i = mstart; i < mend; i++) {
    459:     idx[0] = i - 1;
    460:     idx[1] = i;
    461:     idx[2] = i + 1;
    462:     PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
    463:   }

    465:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    466:      Complete the matrix assembly process and set some options
    467:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    468:   /*
    469:      Assemble matrix, using the 2-step process:
    470:        MatAssemblyBegin(), MatAssemblyEnd()
    471:      Computations can be done while messages are in transition
    472:      by placing code between these two statements.
    473:   */
    474:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    475:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    477:   /*
    478:      Set and option to indicate that we will never add a new nonzero location
    479:      to the matrix. If we do, it will generate an error.
    480:   */
    481:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
    482:   PetscFunctionReturn(PETSC_SUCCESS);
    483: }
    484: /* --------------------------------------------------------------------- */
    485: /*
    486:    Input Parameters:
    487:    ts - the TS context
    488:    t - current time
    489:    f - function
    490:    ctx - optional user-defined context, as set by TSetBCFunction()
    491:  */
    492: PetscErrorCode MyBCRoutine(TS ts, PetscReal t, Vec f, void *ctx)
    493: {
    494:   AppCtx      *appctx = (AppCtx *)ctx; /* user-defined application context */
    495:   PetscInt     m      = appctx->m;
    496:   PetscScalar *fa;

    498:   PetscFunctionBeginUser;
    499:   PetscCall(VecGetArray(f, &fa));
    500:   fa[0]     = 0.0;
    501:   fa[m - 1] = 1.0;
    502:   PetscCall(VecRestoreArray(f, &fa));
    503:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "t=%g\n", (double)t));
    504:   PetscFunctionReturn(PETSC_SUCCESS);
    505: }

    507: /*TEST

    509:     test:
    510:       args: -nox -ts_max_steps 4

    512: TEST*/