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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: PetscFunctionBeginUser;
93: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
94: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
95: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
97: m = 60;
98: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
99: PetscCall(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: PetscCall(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: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &u));
115: PetscCall(VecDuplicate(u, &appctx.solution));
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Set up displays to show graphs of the solution and error
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 380, 400, 160, &appctx.viewer1));
122: PetscCall(PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw));
123: PetscCall(PetscDrawSetDoubleBuffer(draw));
124: PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 0, 400, 160, &appctx.viewer2));
125: PetscCall(PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw));
126: PetscCall(PetscDrawSetDoubleBuffer(draw));
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Create timestepping solver context
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
133: PetscCall(TSSetProblemType(ts, TS_LINEAR));
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Set optional user-defined monitoring routine
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Create matrix data structure; set matrix evaluation routine.
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: PetscCall(MatCreate(PETSC_COMM_SELF, &A));
147: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
148: PetscCall(MatSetFromOptions(A));
149: PetscCall(MatSetUp(A));
151: PetscCall(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: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
159: PetscCall(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: PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx));
168: PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
169: PetscCall(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: PetscCall(TSSetTimeStep(ts, dt));
178: PetscCall(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: PetscCall(TSSetMaxSteps(ts, time_steps_max));
191: PetscCall(TSSetMaxTime(ts, time_total_max));
192: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
193: PetscCall(TSSetFromOptions(ts));
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Solve the problem
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199: /*
200: Evaluate initial conditions
201: */
202: PetscCall(InitialConditions(u, &appctx));
204: /*
205: Run the timestepping solver
206: */
207: PetscCall(TSSolve(ts, u));
208: PetscCall(TSGetSolveTime(ts, &ftime));
209: PetscCall(TSGetStepNumber(ts, &steps));
211: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212: View timestepping solver info
213: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
215: 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)));
216: PetscCall(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: PetscCall(TSDestroy(&ts));
224: PetscCall(MatDestroy(&A));
225: PetscCall(VecDestroy(&u));
226: PetscCall(PetscViewerDestroy(&appctx.viewer1));
227: PetscCall(PetscViewerDestroy(&appctx.viewer2));
228: PetscCall(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: PetscCall(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: PetscFunctionBeginUser;
256: /*
257: Get a pointer to vector data.
258: - For default PETSc vectors, VecGetArray() returns a pointer to
259: the data array. Otherwise, the routine is implementation dependent.
260: - You MUST call VecRestoreArray() when you no longer need access to
261: the array.
262: - Note that the Fortran interface to VecGetArray() differs from the
263: C version. See the users manual for details.
264: */
265: PetscCall(VecGetArray(u, &u_localptr));
267: /*
268: We initialize the solution array by simply writing the solution
269: directly into the array locations. Alternatively, we could use
270: VecSetValues() or VecSetValuesLocal().
271: */
272: for (i = 0; i < appctx->m; i++) u_localptr[i] = PetscCosScalar(PETSC_PI * i * 6. * h) + 3. * PetscCosScalar(PETSC_PI * i * 2. * h);
274: /*
275: Restore vector
276: */
277: PetscCall(VecRestoreArray(u, &u_localptr));
279: /*
280: Print debugging information if desired
281: */
282: if (appctx->debug) {
283: printf("initial guess vector\n");
284: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
285: }
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, tc = t;
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 = PetscExpScalar(-36. * PETSC_PI * PETSC_PI * tc);
317: ex2 = PetscExpScalar(-4. * PETSC_PI * PETSC_PI * tc);
318: sc1 = PETSC_PI * 6. * h;
319: sc2 = PETSC_PI * 2. * h;
320: for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscCosScalar(sc1 * (PetscReal)i) * ex1 + 3. * PetscCosScalar(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: Input Parameters:
335: ts - the timestep context
336: step - the count of the current step (with 0 meaning the
337: initial condition)
338: time - the current time
339: u - the solution at this timestep
340: ctx - the user-provided context for this monitoring routine.
341: In this case we use the application context which contains
342: information about the problem size, workspace and the exact
343: solution.
344: */
345: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx)
346: {
347: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
348: PetscReal norm_2, norm_max;
350: PetscFunctionBeginUser;
351: /*
352: View a graph of the current iterate
353: */
354: PetscCall(VecView(u, appctx->viewer2));
356: /*
357: Compute the exact solution
358: */
359: PetscCall(ExactSolution(time, appctx->solution, appctx));
361: /*
362: Print debugging information if desired
363: */
364: if (appctx->debug) {
365: printf("Computed solution vector\n");
366: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
367: printf("Exact solution vector\n");
368: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
369: }
371: /*
372: Compute the 2-norm and max-norm of the error
373: */
374: PetscCall(VecAXPY(appctx->solution, -1.0, u));
375: PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
376: norm_2 = PetscSqrtReal(appctx->h) * norm_2;
377: PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
378: if (norm_2 < 1e-14) norm_2 = 0;
379: if (norm_max < 1e-14) norm_max = 0;
381: 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));
382: appctx->norm_2 += norm_2;
383: appctx->norm_max += norm_max;
385: /*
386: View a graph of the error
387: */
388: PetscCall(VecView(appctx->solution, appctx->viewer1));
390: /*
391: Print debugging information if desired
392: */
393: if (appctx->debug) {
394: printf("Error vector\n");
395: PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
396: }
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
400: /* --------------------------------------------------------------------- */
401: /*
402: RHSMatrixHeat - User-provided routine to compute the right-hand-side
403: matrix for the heat equation.
405: Input Parameters:
406: ts - the TS context
407: t - current time
408: global_in - global input vector
409: dummy - optional user-defined context, as set by TSetRHSJacobian()
411: Output Parameters:
412: AA - Jacobian matrix
413: BB - optionally different preconditioning matrix
414: str - flag indicating matrix structure
416: Notes:
417: Recall that MatSetValues() uses 0-based row and column numbers
418: in Fortran as well as in C.
419: */
420: PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)
421: {
422: Mat A = AA; /* Jacobian matrix */
423: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
424: PetscInt mstart = 0;
425: PetscInt mend = appctx->m;
426: PetscInt i, idx[3];
427: PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
429: PetscFunctionBeginUser;
430: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
431: Compute entries for the locally owned part of the matrix
432: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
433: /*
434: Set matrix rows corresponding to boundary data
435: */
437: mstart = 0;
438: v[0] = 1.0;
439: PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
440: mstart++;
442: mend--;
443: v[0] = 1.0;
444: PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
446: /*
447: Set matrix rows corresponding to interior data. We construct the
448: matrix one row at a time.
449: */
450: v[0] = sone;
451: v[1] = stwo;
452: v[2] = sone;
453: for (i = mstart; i < mend; i++) {
454: idx[0] = i - 1;
455: idx[1] = i;
456: idx[2] = i + 1;
457: PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
458: }
460: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
461: Complete the matrix assembly process and set some options
462: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
463: /*
464: Assemble matrix, using the 2-step process:
465: MatAssemblyBegin(), MatAssemblyEnd()
466: Computations can be done while messages are in transition
467: by placing code between these two statements.
468: */
469: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
470: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
472: /*
473: Set and option to indicate that we will never add a new nonzero location
474: to the matrix. If we do, it will generate an error.
475: */
476: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*TEST
483: test:
484: requires: x
486: test:
487: suffix: nox
488: args: -nox
490: TEST*/