Actual source code: ex6.c


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

  9: /* ------------------------------------------------------------------------

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

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

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

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

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

 39:   ------------------------------------------------------------------------- */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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

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

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

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

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

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

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

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

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

218:   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)));
219:   PetscCall(TSView(ts, PETSC_VIEWER_STDOUT_SELF));

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

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

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

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

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

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

270:   /*
271:      We initialize the solution array by simply writing the solution
272:      directly into the array locations.  Alternatively, we could use
273:      VecSetValues() or VecSetValuesLocal().
274:   */
275:   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);

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

282:   /*
283:      Print debugging information if desired
284:   */
285:   if (appctx->debug) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));

287:   PetscFunctionReturn(PETSC_SUCCESS);
288: }
289: /* --------------------------------------------------------------------- */
290: /*
291:    ExactSolution - Computes the exact solution at a given time.

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

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

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

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

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

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

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

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

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

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

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

382:   PetscCall(TSGetTimeStep(ts, &dt));
383:   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));
384:   appctx->norm_2 += norm_2;
385:   appctx->norm_max += norm_max;

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

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

399:   /*
400:      Print debugging information if desired
401:   */
402:   if (appctx->debug) {
403:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
404:     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
405:   }

407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }
409: /* --------------------------------------------------------------------- */
410: /*
411:    RHSMatrixHeat - User-provided routine to compute the right-hand-side
412:    matrix for the heat equation.

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

420:    Output Parameters:
421:    AA - Jacobian matrix
422:    BB - optionally different preconditioning matrix
423:    str - flag indicating matrix structure

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

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

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

451:   mend--;
452:   v[0] = 1.0;
453:   PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));

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

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

481:   /*
482:      Set and option to indicate that we will never add a new nonzero location
483:      to the matrix. If we do, it will generate an error.
484:   */
485:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));

487:   PetscFunctionReturn(PETSC_SUCCESS);
488: }
489: /* --------------------------------------------------------------------- */
490: /*
491:    Input Parameters:
492:    ts - the TS context
493:    t - current time
494:    f - function
495:    ctx - optional user-defined context, as set by TSetBCFunction()
496:  */
497: PetscErrorCode MyBCRoutine(TS ts, PetscReal t, Vec f, void *ctx)
498: {
499:   AppCtx      *appctx = (AppCtx *)ctx; /* user-defined application context */
500:   PetscInt     m      = appctx->m;
501:   PetscScalar *fa;

503:   PetscFunctionBeginUser;
504:   PetscCall(VecGetArray(f, &fa));
505:   fa[0]     = 0.0;
506:   fa[m - 1] = 1.0;
507:   PetscCall(VecRestoreArray(f, &fa));
508:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "t=%g\n", (double)t));

510:   PetscFunctionReturn(PETSC_SUCCESS);
511: }

513: /*TEST

515:     test:
516:       args: -nox -ts_max_steps 4

518: TEST*/