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:   PetscInitialize(&argc,&argv,(char*)0,help);
 95:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 98:   m    = 60;
 99:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
100:   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:   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:   VecCreateSeq(PETSC_COMM_SELF,m,&u);
117:   VecDuplicate(u,&appctx.solution);

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

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

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

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

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

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

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

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

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

153:   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:     TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);
161:     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:     RHSMatrixHeat(ts,0.0,u,A,A,&appctx);
170:     TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);
171:     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:   TSSetTimeStep(ts,dt);
180:   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:   TSSetMaxSteps(ts,time_steps_max);
193:   TSSetMaxTime(ts,time_total_max);
194:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
195:   TSSetFromOptions(ts);

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

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

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

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

217:   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:   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:   TSDestroy(&ts);
226:   MatDestroy(&A);
227:   VecDestroy(&u);
228:   PetscViewerDestroy(&appctx.viewer1);
229:   PetscViewerDestroy(&appctx.viewer2);
230:   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:   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:   /*
258:     Get a pointer to vector data.
259:     - For default PETSc vectors, VecGetArray() returns a pointer to
260:       the data array.  Otherwise, the routine is implementation dependent.
261:     - You MUST call VecRestoreArray() when you no longer need access to
262:       the array.
263:     - Note that the Fortran interface to VecGetArray() differs from the
264:       C version.  See the users manual for details.
265:   */
266:   VecGetArray(u,&u_localptr);

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

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

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

287:   return 0;
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:   /*
307:      Get a pointer to vector data.
308:   */
309:   VecGetArray(solution,&s_localptr);

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

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

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

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

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

355:   /*
356:      Compute the exact solution
357:   */
358:   ExactSolution(crtime,appctx->solution,appctx);

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

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

378:   TSGetTimeStep(ts,&dt);
379:   if (norm_2 > 1.e-2) {
380:     PetscPrintf(PETSC_COMM_SELF,"Timestep %D: 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);
381:   }
382:   appctx->norm_2   += norm_2;
383:   appctx->norm_max += norm_max;

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

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

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

405:   return 0;
406: }
407: /* --------------------------------------------------------------------- */
408: /*
409:    RHSMatrixHeat - User-provided routine to compute the right-hand-side
410:    matrix for the heat equation.

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

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

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

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:   MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
446:   mstart++;

448:   mend--;
449:   v[0] = 1.0;
450:   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; v[1] = stwo; v[2] = sone;
457:   for (i=mstart; i<mend; i++) {
458:     idx[0] = i-1; idx[1] = i; idx[2] = i+1;
459:     MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
460:   }

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

474:   /*
475:      Set and option to indicate that we will never add a new nonzero location
476:      to the matrix. If we do, it will generate an error.
477:   */
478:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);

480:   return 0;
481: }
482: /* --------------------------------------------------------------------- */
483: /*
484:    Input Parameters:
485:    ts - the TS context
486:    t - current time
487:    f - function
488:    ctx - optional user-defined context, as set by TSetBCFunction()
489:  */
490: PetscErrorCode MyBCRoutine(TS ts,PetscReal t,Vec f,void *ctx)
491: {
492:   AppCtx         *appctx = (AppCtx*) ctx;      /* user-defined application context */
493:   PetscInt       m = appctx->m;
494:   PetscScalar    *fa;

496:   VecGetArray(f,&fa);
497:   fa[0]   = 0.0;
498:   fa[m-1] = 1.0;
499:   VecRestoreArray(f,&fa);
500:   PetscPrintf(PETSC_COMM_SELF,"t=%g\n",(double)t);

502:   return 0;
503: }

505: /*TEST

507:     test:
508:       args: -nox -ts_max_steps 4

510: TEST*/