Actual source code: ex5.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) = 1, u(t,1) = 1,
 16:    and the initial condition
 17:        u(0,x) = cos(6*pi*x) + 3*cos(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) * cos(6*pi*x) +
 29:                       3*exp(-4*pi*pi*t) * cos(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 just f(u)

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

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

 41: /*
 42:    Include "petscts.h" so that we can use TS solvers.  Note that this file
 43:    automatically includes:
 44:      petscsys.h       - base PETSc routines   petscvec.h  - vectors
 45:      petscmat.h  - matrices
 46:      petscis.h     - index sets            petscksp.h  - Krylov subspace methods
 47:      petscviewer.h - viewers               petscpc.h   - preconditioners
 48:      petscksp.h   - linear solvers        petscsnes.h - nonlinear solvers
 49: */
 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 *);

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

 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:      Initialize program and set problem parameters
 90:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 93:   PetscInitialize(&argc, &argv, (char *)0, help);
 94:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 97:   m = 60;
 98:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 99:   PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug);
100:   appctx.m        = m;
101:   appctx.h        = 1.0 / (m - 1.0);
102:   appctx.norm_2   = 0.0;
103:   appctx.norm_max = 0.0;

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

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Create vector data structures
109:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:      Create timestepping solver context
130:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

132:   TSCreate(PETSC_COMM_SELF, &ts);
133:   TSSetProblemType(ts, TS_LINEAR);

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Set optional user-defined monitoring routine
137:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

139:   TSMonitorSet(ts, Monitor, &appctx, NULL);

141:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

143:      Create matrix data structure; set matrix evaluation routine.
144:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

146:   MatCreate(PETSC_COMM_SELF, &A);
147:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m);
148:   MatSetFromOptions(A);
149:   MatSetUp(A);

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

172:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173:      Set solution vector and initial timestep
174:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

176:   dt = appctx.h * appctx.h / 2.0;
177:   TSSetTimeStep(ts, dt);
178:   TSSetSolution(ts, u);

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

190:   TSSetMaxSteps(ts, time_steps_max);
191:   TSSetMaxTime(ts, time_total_max);
192:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
193:   TSSetFromOptions(ts);

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196:      Solve the problem
197:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

199:   /*
200:      Evaluate initial conditions
201:   */
202:   InitialConditions(u, &appctx);

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

211:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212:      View timestepping solver info
213:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

215:   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));
216:   TSView(ts, PETSC_VIEWER_STDOUT_SELF);

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

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

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

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

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

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

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

273:   /*
274:      Restore vector
275:   */
276:   VecRestoreArray(u, &u_localptr);

278:   /*
279:      Print debugging information if desired
280:   */
281:   if (appctx->debug) {
282:     printf("initial guess vector\n");
283:     VecView(u, PETSC_VIEWER_STDOUT_SELF);
284:   }

286:   return 0;
287: }
288: /* --------------------------------------------------------------------- */
289: /*
290:    ExactSolution - Computes the exact solution at a given time.

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

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

305:   /*
306:      Get a pointer to vector data.
307:   */
308:   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 = PetscExpScalar(-36. * PETSC_PI * PETSC_PI * tc);
315:   ex2 = PetscExpScalar(-4. * PETSC_PI * PETSC_PI * tc);
316:   sc1 = PETSC_PI * 6. * h;
317:   sc2 = PETSC_PI * 2. * h;
318:   for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscCosScalar(sc1 * (PetscReal)i) * ex1 + 3. * PetscCosScalar(sc2 * (PetscReal)i) * ex2;

320:   /*
321:      Restore vector
322:   */
323:   VecRestoreArray(solution, &s_localptr);
324:   return 0;
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:    Input Parameters:
333:    ts     - the timestep context
334:    step   - the count of the current step (with 0 meaning the
335:              initial condition)
336:    time   - the current time
337:    u      - the solution at this timestep
338:    ctx    - the user-provided context for this monitoring routine.
339:             In this case we use the application context which contains
340:             information about the problem size, workspace and the exact
341:             solution.
342: */
343: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx)
344: {
345:   AppCtx   *appctx = (AppCtx *)ctx; /* user-defined application context */
346:   PetscReal norm_2, norm_max;

348:   /*
349:      View a graph of the current iterate
350:   */
351:   VecView(u, appctx->viewer2);

353:   /*
354:      Compute the exact solution
355:   */
356:   ExactSolution(time, appctx->solution, appctx);

358:   /*
359:      Print debugging information if desired
360:   */
361:   if (appctx->debug) {
362:     printf("Computed solution vector\n");
363:     VecView(u, PETSC_VIEWER_STDOUT_SELF);
364:     printf("Exact solution vector\n");
365:     VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF);
366:   }

368:   /*
369:      Compute the 2-norm and max-norm of the error
370:   */
371:   VecAXPY(appctx->solution, -1.0, u);
372:   VecNorm(appctx->solution, NORM_2, &norm_2);
373:   norm_2 = PetscSqrtReal(appctx->h) * norm_2;
374:   VecNorm(appctx->solution, NORM_MAX, &norm_max);
375:   if (norm_2 < 1e-14) norm_2 = 0;
376:   if (norm_max < 1e-14) norm_max = 0;

378:   PetscPrintf(PETSC_COMM_WORLD, "Timestep %" PetscInt_FMT ": time = %g, 2-norm error = %g, max norm error = %g\n", step, (double)time, (double)norm_2, (double)norm_max);
379:   appctx->norm_2 += norm_2;
380:   appctx->norm_max += norm_max;

382:   /*
383:      View a graph of the error
384:   */
385:   VecView(appctx->solution, appctx->viewer1);

387:   /*
388:      Print debugging information if desired
389:   */
390:   if (appctx->debug) {
391:     printf("Error vector\n");
392:     VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF);
393:   }

395:   return 0;
396: }
397: /* --------------------------------------------------------------------- */
398: /*
399:    RHSMatrixHeat - User-provided routine to compute the right-hand-side
400:    matrix for the heat equation.

402:    Input Parameters:
403:    ts - the TS context
404:    t - current time
405:    global_in - global input vector
406:    dummy - optional user-defined context, as set by TSetRHSJacobian()

408:    Output Parameters:
409:    AA - Jacobian matrix
410:    BB - optionally different preconditioning matrix
411:    str - flag indicating matrix structure

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

426:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
427:      Compute entries for the locally owned part of the matrix
428:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
429:   /*
430:      Set matrix rows corresponding to boundary data
431:   */

433:   mstart = 0;
434:   v[0]   = 1.0;
435:   MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES);
436:   mstart++;

438:   mend--;
439:   v[0] = 1.0;
440:   MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES);

442:   /*
443:      Set matrix rows corresponding to interior data.  We construct the
444:      matrix one row at a time.
445:   */
446:   v[0] = sone;
447:   v[1] = stwo;
448:   v[2] = sone;
449:   for (i = mstart; i < mend; i++) {
450:     idx[0] = i - 1;
451:     idx[1] = i;
452:     idx[2] = i + 1;
453:     MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES);
454:   }

456:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
457:      Complete the matrix assembly process and set some options
458:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
459:   /*
460:      Assemble matrix, using the 2-step process:
461:        MatAssemblyBegin(), MatAssemblyEnd()
462:      Computations can be done while messages are in transition
463:      by placing code between these two statements.
464:   */
465:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
466:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

468:   /*
469:      Set and option to indicate that we will never add a new nonzero location
470:      to the matrix. If we do, it will generate an error.
471:   */
472:   MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);

474:   return 0;
475: }

477: /*TEST

479:     test:
480:       requires: x

482:     test:
483:       suffix: nox
484:       args: -nox

486: TEST*/