Actual source code: ex5.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) = 1, u(t,1) = 1,
15: and the initial condition
16: u(0,x) = cos(6*pi*x) + 3*cos(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) * cos(6*pi*x) +
28: 3*exp(-4*pi*pi*t) * cos(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 just f(u)
36: The parallel version of this code is ts/tutorials/ex4.c
38: ------------------------------------------------------------------------- */
40: /*
41: Include "petscts.h" so that we can use TS solvers. Note that this file
42: automatically includes:
43: petscsys.h - base PETSc routines petscvec.h - vectors
44: petscmat.h - matrices
45: petscis.h - index sets petscksp.h - Krylov subspace methods
46: petscviewer.h - viewers petscpc.h - preconditioners
47: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
48: */
49: #include <petscts.h>
50: #include <petscdraw.h>
52: /*
53: User-defined application context - contains data needed by the
54: application-provided call-back routines.
55: */
56: typedef struct {
57: Vec solution; /* global exact solution vector */
58: PetscInt m; /* total number of grid points */
59: PetscReal h; /* mesh width h = 1/(m-1) */
60: PetscBool debug; /* flag (1 indicates activation of debugging printouts) */
61: PetscViewer viewer1, viewer2; /* viewers for the solution and error */
62: PetscReal norm_2, norm_max; /* error norms */
63: } AppCtx;
65: /*
66: User-defined routines
67: */
68: extern PetscErrorCode InitialConditions(Vec, AppCtx *);
69: extern PetscErrorCode RHSMatrixHeat(TS, PetscReal, Vec, Mat, Mat, void *);
70: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
71: extern PetscErrorCode ExactSolution(PetscReal, Vec, AppCtx *);
73: int main(int argc, char **argv)
74: {
75: AppCtx appctx; /* user-defined application context */
76: TS ts; /* timestepping context */
77: Mat A; /* matrix data structure */
78: Vec u; /* approximate solution vector */
79: PetscReal time_total_max = 100.0; /* default max total time */
80: PetscInt time_steps_max = 100; /* default max timesteps */
81: PetscDraw draw; /* drawing context */
82: PetscInt steps, m;
83: PetscMPIInt size;
84: PetscBool flg;
85: PetscReal dt, ftime;
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Initialize program and set problem parameters
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: PetscFunctionBeginUser;
92: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
93: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
94: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
96: m = 60;
97: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
98: PetscCall(PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug));
99: appctx.m = m;
100: appctx.h = 1.0 / (m - 1.0);
101: appctx.norm_2 = 0.0;
102: appctx.norm_max = 0.0;
104: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n"));
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create vector data structures
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Create vector data structures for approximate and exact solutions
112: */
113: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &u));
114: PetscCall(VecDuplicate(u, &appctx.solution));
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Set up displays to show graphs of the solution and error
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 380, 400, 160, &appctx.viewer1));
121: PetscCall(PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw));
122: PetscCall(PetscDrawSetDoubleBuffer(draw));
123: PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 0, 400, 160, &appctx.viewer2));
124: PetscCall(PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw));
125: PetscCall(PetscDrawSetDoubleBuffer(draw));
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Create timestepping solver context
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
132: PetscCall(TSSetProblemType(ts, TS_LINEAR));
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Set optional user-defined monitoring routine
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Create matrix data structure; set matrix evaluation routine.
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: PetscCall(MatCreate(PETSC_COMM_SELF, &A));
146: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
147: PetscCall(MatSetFromOptions(A));
148: PetscCall(MatSetUp(A));
150: PetscCall(PetscOptionsHasName(NULL, NULL, "-time_dependent_rhs", &flg));
151: if (flg) {
152: /*
153: For linear problems with a time-dependent f(u,t) in the equation
154: u_t = f(u,t), the user provides the discretized right-hand-side
155: as a time-dependent matrix.
156: */
157: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
158: PetscCall(TSSetRHSJacobian(ts, A, A, RHSMatrixHeat, &appctx));
159: } else {
160: /*
161: For linear problems with a time-independent f(u) in the equation
162: u_t = f(u), the user provides the discretized right-hand-side
163: as a matrix only once, and then sets a null matrix evaluation
164: routine.
165: */
166: PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx));
167: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
168: PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &appctx));
169: }
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Set solution vector and initial timestep
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: dt = appctx.h * appctx.h / 2.0;
176: PetscCall(TSSetTimeStep(ts, dt));
177: PetscCall(TSSetSolution(ts, u));
179: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: Customize timestepping solver:
181: - Set the solution method to be the Backward Euler method.
182: - Set timestepping duration info
183: Then set runtime options, which can override these defaults.
184: For example,
185: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
186: to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: PetscCall(TSSetMaxSteps(ts, time_steps_max));
190: PetscCall(TSSetMaxTime(ts, time_total_max));
191: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
192: PetscCall(TSSetFromOptions(ts));
194: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: Solve the problem
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: /*
199: Evaluate initial conditions
200: */
201: PetscCall(InitialConditions(u, &appctx));
203: /*
204: Run the timestepping solver
205: */
206: PetscCall(TSSolve(ts, u));
207: PetscCall(TSGetSolveTime(ts, &ftime));
208: PetscCall(TSGetStepNumber(ts, &steps));
210: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: View timestepping solver info
212: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
214: 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)));
215: PetscCall(TSView(ts, PETSC_VIEWER_STDOUT_SELF));
217: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218: Free work space. All PETSc objects should be destroyed when they
219: are no longer needed.
220: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222: PetscCall(TSDestroy(&ts));
223: PetscCall(MatDestroy(&A));
224: PetscCall(VecDestroy(&u));
225: PetscCall(PetscViewerDestroy(&appctx.viewer1));
226: PetscCall(PetscViewerDestroy(&appctx.viewer2));
227: PetscCall(VecDestroy(&appctx.solution));
229: /*
230: Always call PetscFinalize() before exiting a program. This routine
231: - finalizes the PETSc libraries as well as MPI
232: - provides summary and diagnostic information if certain runtime
233: options are chosen (e.g., -log_view).
234: */
235: PetscCall(PetscFinalize());
236: return 0;
237: }
238: /* --------------------------------------------------------------------- */
239: /*
240: InitialConditions - Computes the solution at the initial time.
242: Input Parameter:
243: u - uninitialized solution vector (global)
244: appctx - user-defined application context
246: Output Parameter:
247: u - vector with solution at initial time (global)
248: */
249: PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
250: {
251: PetscScalar *u_localptr, h = appctx->h;
252: PetscInt i;
254: PetscFunctionBeginUser;
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: PetscCall(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: PetscCall(VecRestoreArray(u, &u_localptr));
278: /*
279: Print debugging information if desired
280: */
281: if (appctx->debug) {
282: printf("initial guess vector\n");
283: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
284: }
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, tc = t;
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 = PetscExpScalar(-36. * PETSC_PI * PETSC_PI * tc);
316: ex2 = PetscExpScalar(-4. * PETSC_PI * PETSC_PI * tc);
317: sc1 = PETSC_PI * 6. * h;
318: sc2 = PETSC_PI * 2. * h;
319: for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscCosScalar(sc1 * (PetscReal)i) * ex1 + 3. * PetscCosScalar(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: Input Parameters:
334: ts - the timestep context
335: step - the count of the current step (with 0 meaning the
336: initial condition)
337: time - 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 time, Vec u, void *ctx)
345: {
346: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
347: PetscReal norm_2, norm_max;
349: PetscFunctionBeginUser;
350: /*
351: View a graph of the current iterate
352: */
353: PetscCall(VecView(u, appctx->viewer2));
355: /*
356: Compute the exact solution
357: */
358: PetscCall(ExactSolution(time, appctx->solution, appctx));
360: /*
361: Print debugging information if desired
362: */
363: if (appctx->debug) {
364: printf("Computed solution vector\n");
365: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
366: printf("Exact solution vector\n");
367: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
368: }
370: /*
371: Compute the 2-norm and max-norm of the error
372: */
373: PetscCall(VecAXPY(appctx->solution, -1.0, u));
374: PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
375: norm_2 = PetscSqrtReal(appctx->h) * norm_2;
376: PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
377: if (norm_2 < 1e-14) norm_2 = 0;
378: if (norm_max < 1e-14) norm_max = 0;
380: PetscCall(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));
381: appctx->norm_2 += norm_2;
382: appctx->norm_max += norm_max;
384: /*
385: View a graph of the error
386: */
387: PetscCall(VecView(appctx->solution, appctx->viewer1));
389: /*
390: Print debugging information if desired
391: */
392: if (appctx->debug) {
393: printf("Error vector\n");
394: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
395: }
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
399: /* --------------------------------------------------------------------- */
400: /*
401: RHSMatrixHeat - User-provided routine to compute the right-hand-side
402: matrix for the heat equation.
404: Input Parameters:
405: ts - the TS context
406: t - current time
407: global_in - global input vector
408: dummy - optional user-defined context, as set by TSetRHSJacobian()
410: Output Parameters:
411: AA - Jacobian matrix
412: BB - optionally different preconditioning matrix
413: str - flag indicating matrix structure
415: Notes:
416: Recall that MatSetValues() uses 0-based row and column numbers
417: in Fortran as well as in C.
418: */
419: PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)
420: {
421: Mat A = AA; /* Jacobian matrix */
422: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
423: PetscInt mstart = 0;
424: PetscInt mend = appctx->m;
425: PetscInt i, idx[3];
426: PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
428: PetscFunctionBeginUser;
429: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
430: Compute entries for the locally owned part of the matrix
431: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
432: /*
433: Set matrix rows corresponding to boundary data
434: */
436: mstart = 0;
437: v[0] = 1.0;
438: PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
439: mstart++;
441: mend--;
442: v[0] = 1.0;
443: PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
445: /*
446: Set matrix rows corresponding to interior data. We construct the
447: matrix one row at a time.
448: */
449: v[0] = sone;
450: v[1] = stwo;
451: v[2] = sone;
452: for (i = mstart; i < mend; i++) {
453: idx[0] = i - 1;
454: idx[1] = i;
455: idx[2] = i + 1;
456: PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
457: }
459: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460: Complete the matrix assembly process and set some options
461: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
462: /*
463: Assemble matrix, using the 2-step process:
464: MatAssemblyBegin(), MatAssemblyEnd()
465: Computations can be done while messages are in transition
466: by placing code between these two statements.
467: */
468: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
469: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
471: /*
472: Set and option to indicate that we will never add a new nonzero location
473: to the matrix. If we do, it will generate an error.
474: */
475: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: /*TEST
482: test:
483: requires: x
485: test:
486: suffix: nox
487: args: -nox
489: TEST*/