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));
286: PetscFunctionReturn(PETSC_SUCCESS);
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;
303: PetscInt i;
305: PetscFunctionBeginUser;
306: /*
307: Get a pointer to vector data.
308: */
309: PetscCall(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);
316: ex2 = PetscExpReal(-4. * PETSC_PI * PETSC_PI * t);
317: sc1 = PETSC_PI * 6. * h;
318: sc2 = PETSC_PI * 2. * h;
319: for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscSinReal(PetscRealPart(sc1) * (PetscReal)i) * ex1 + 3. * PetscSinReal(PetscRealPart(sc2) * (PetscReal)i) * ex2;
321: /*
322: Restore vector
323: */
324: PetscCall(VecRestoreArray(solution, &s_localptr));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
327: /* --------------------------------------------------------------------- */
328: /*
329: Monitor - User-provided routine to monitor the solution computed at
330: each timestep. This example plots the solution and computes the
331: error in two different norms.
333: This example also demonstrates changing the timestep via TSSetTimeStep().
335: Input Parameters:
336: ts - the timestep context
337: step - the count of the current step (with 0 meaning the
338: initial condition)
339: crtime - the current time
340: u - the solution at this timestep
341: ctx - the user-provided context for this monitoring routine.
342: In this case we use the application context which contains
343: information about the problem size, workspace and the exact
344: solution.
345: */
346: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
347: {
348: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
349: PetscReal norm_2, norm_max, dt, dttol;
350: PetscBool flg;
352: PetscFunctionBeginUser;
353: /*
354: View a graph of the current iterate
355: */
356: PetscCall(VecView(u, appctx->viewer2));
358: /*
359: Compute the exact solution
360: */
361: PetscCall(ExactSolution(crtime, appctx->solution, appctx));
363: /*
364: Print debugging information if desired
365: */
366: if (appctx->debug) {
367: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n"));
368: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
369: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
370: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
371: }
373: /*
374: Compute the 2-norm and max-norm of the error
375: */
376: PetscCall(VecAXPY(appctx->solution, -1.0, u));
377: PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
378: norm_2 = PetscSqrtReal(appctx->h) * norm_2;
379: PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
381: PetscCall(TSGetTimeStep(ts, &dt));
382: 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));
383: appctx->norm_2 += norm_2;
384: appctx->norm_max += norm_max;
386: dttol = .0001;
387: PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, &flg));
388: if (dt < dttol) {
389: dt *= .999;
390: PetscCall(TSSetTimeStep(ts, dt));
391: }
393: /*
394: View a graph of the error
395: */
396: PetscCall(VecView(appctx->solution, appctx->viewer1));
398: /*
399: Print debugging information if desired
400: */
401: if (appctx->debug) {
402: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
403: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
404: }
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
408: /* --------------------------------------------------------------------- */
409: /*
410: RHSMatrixHeat - User-provided routine to compute the right-hand-side
411: matrix for the heat equation.
413: Input Parameters:
414: ts - the TS context
415: t - current time
416: global_in - global input vector
417: dummy - optional user-defined context, as set by TSetRHSJacobian()
419: Output Parameters:
420: AA - Jacobian matrix
421: BB - optionally different preconditioning matrix
422: str - flag indicating matrix structure
424: Notes:
425: Recall that MatSetValues() uses 0-based row and column numbers
426: in Fortran as well as in C.
427: */
428: PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)
429: {
430: Mat A = AA; /* Jacobian matrix */
431: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
432: PetscInt mstart = 0;
433: PetscInt mend = appctx->m;
434: PetscInt i, idx[3];
435: PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
437: PetscFunctionBeginUser;
438: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
439: Compute entries for the locally owned part of the matrix
440: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
441: /*
442: Set matrix rows corresponding to boundary data
443: */
445: mstart = 0;
446: v[0] = 1.0;
447: PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
448: mstart++;
450: mend--;
451: v[0] = 1.0;
452: PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
454: /*
455: Set matrix rows corresponding to interior data. We construct the
456: matrix one row at a time.
457: */
458: v[0] = sone;
459: v[1] = stwo;
460: v[2] = sone;
461: for (i = mstart; i < mend; i++) {
462: idx[0] = i - 1;
463: idx[1] = i;
464: idx[2] = i + 1;
465: PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
466: }
468: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
469: Complete the matrix assembly process and set some options
470: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
471: /*
472: Assemble matrix, using the 2-step process:
473: MatAssemblyBegin(), MatAssemblyEnd()
474: Computations can be done while messages are in transition
475: by placing code between these two statements.
476: */
477: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
478: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
480: /*
481: Set and option to indicate that we will never add a new nonzero location
482: to the matrix. If we do, it will generate an error.
483: */
484: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
488: /* --------------------------------------------------------------------- */
489: /*
490: Input Parameters:
491: ts - the TS context
492: t - current time
493: f - function
494: ctx - optional user-defined context, as set by TSetBCFunction()
495: */
496: PetscErrorCode MyBCRoutine(TS ts, PetscReal t, Vec f, void *ctx)
497: {
498: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
499: PetscInt m = appctx->m;
500: PetscScalar *fa;
502: PetscFunctionBeginUser;
503: PetscCall(VecGetArray(f, &fa));
504: fa[0] = 0.0;
505: fa[m - 1] = 1.0;
506: PetscCall(VecRestoreArray(f, &fa));
507: PetscCall(PetscPrintf(PETSC_COMM_SELF, "t=%g\n", (double)t));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: /*TEST
514: test:
515: args: -nox -ts_max_steps 4
517: TEST*/