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*/