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, (char *)0, 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 preconditioning matrix
420:    str - flag indicating matrix structure

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

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

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

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

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

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

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

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

508: /*TEST

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

513: TEST*/