Actual source code: ex3.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: -use_ifunc : Use IFunction/IJacobian interface\n\
7: -debug : Activate debugging printouts\n\
8: -nox : Deactivate x-window graphics\n\n";
10: /* ------------------------------------------------------------------------
12: This program solves the one-dimensional heat equation (also called the
13: diffusion equation),
14: u_t = u_xx,
15: on the domain 0 <= x <= 1, with the boundary conditions
16: u(t,0) = 0, u(t,1) = 0,
17: and the initial condition
18: u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x).
19: This is a linear, second-order, parabolic equation.
21: We discretize the right-hand side using finite differences with
22: uniform grid spacing h:
23: u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2)
24: We then demonstrate time evolution using the various TS methods by
25: running the program via
26: ex3 -ts_type <timestepping solver>
28: We compare the approximate solution with the exact solution, given by
29: u_exact(x,t) = exp(-36*pi*pi*t) * sin(6*pi*x) +
30: 3*exp(-4*pi*pi*t) * sin(2*pi*x)
32: Notes:
33: This code demonstrates the TS solver interface to two variants of
34: linear problems, u_t = f(u,t), namely
35: - time-dependent f: f(u,t) is a function of t
36: - time-independent f: f(u,t) is simply f(u)
38: The parallel version of this code is ts/tutorials/ex4.c
40: ------------------------------------------------------------------------- */
42: /*
43: Include "petscts.h" so that we can use TS solvers. Note that this file
44: automatically includes:
45: petscsys.h - base PETSc routines petscvec.h - vectors
46: petscmat.h - matrices
47: petscis.h - index sets petscksp.h - Krylov subspace methods
48: petscviewer.h - viewers petscpc.h - preconditioners
49: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
50: */
52: #include <petscts.h>
53: #include <petscdraw.h>
55: /*
56: User-defined application context - contains data needed by the
57: application-provided call-back routines.
58: */
59: typedef struct {
60: Vec solution; /* global exact solution vector */
61: PetscInt m; /* total number of grid points */
62: PetscReal h; /* mesh width h = 1/(m-1) */
63: PetscBool debug; /* flag (1 indicates activation of debugging printouts) */
64: PetscViewer viewer1, viewer2; /* viewers for the solution and error */
65: PetscReal norm_2, norm_max; /* error norms */
66: Mat A; /* RHS mat, used with IFunction interface */
67: PetscReal oshift; /* old shift applied, prevent to recompute the IJacobian */
68: } AppCtx;
70: /*
71: User-defined routines
72: */
73: extern PetscErrorCode InitialConditions(Vec, AppCtx *);
74: extern PetscErrorCode RHSMatrixHeat(TS, PetscReal, Vec, Mat, Mat, void *);
75: extern PetscErrorCode IFunctionHeat(TS, PetscReal, Vec, Vec, Vec, void *);
76: extern PetscErrorCode IJacobianHeat(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
77: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
78: extern PetscErrorCode ExactSolution(PetscReal, Vec, AppCtx *);
80: int main(int argc, char **argv)
81: {
82: AppCtx appctx; /* user-defined application context */
83: TS ts; /* timestepping context */
84: Mat A; /* matrix data structure */
85: Vec u; /* approximate solution vector */
86: PetscReal time_total_max = 100.0; /* default max total time */
87: PetscInt time_steps_max = 100; /* default max timesteps */
88: PetscDraw draw; /* drawing context */
89: PetscInt steps, m;
90: PetscMPIInt size;
91: PetscReal dt;
92: PetscBool flg, flg_string;
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Initialize program and set problem parameters
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PetscInitialize(&argc, &argv, (char *)0, help);
100: MPI_Comm_size(PETSC_COMM_WORLD, &size);
103: m = 60;
104: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
105: PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug);
106: flg_string = PETSC_FALSE;
107: PetscOptionsGetBool(NULL, NULL, "-test_string_viewer", &flg_string, NULL);
109: appctx.m = m;
110: appctx.h = 1.0 / (m - 1.0);
111: appctx.norm_2 = 0.0;
112: appctx.norm_max = 0.0;
114: PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n");
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Create vector data structures
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: /*
121: Create vector data structures for approximate and exact solutions
122: */
123: VecCreateSeq(PETSC_COMM_SELF, m, &u);
124: VecDuplicate(u, &appctx.solution);
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Set up displays to show graphs of the solution and error
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 380, 400, 160, &appctx.viewer1);
131: PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw);
132: PetscDrawSetDoubleBuffer(draw);
133: PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 0, 400, 160, &appctx.viewer2);
134: PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw);
135: PetscDrawSetDoubleBuffer(draw);
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Create timestepping solver context
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: TSCreate(PETSC_COMM_SELF, &ts);
142: TSSetProblemType(ts, TS_LINEAR);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Set optional user-defined monitoring routine
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: if (!flg_string) TSMonitorSet(ts, Monitor, &appctx, NULL);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Create matrix data structure; set matrix evaluation routine.
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: MatCreate(PETSC_COMM_SELF, &A);
156: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m);
157: MatSetFromOptions(A);
158: MatSetUp(A);
160: flg = PETSC_FALSE;
161: PetscOptionsGetBool(NULL, NULL, "-use_ifunc", &flg, NULL);
162: if (!flg) {
163: appctx.A = NULL;
164: PetscOptionsGetBool(NULL, NULL, "-time_dependent_rhs", &flg, NULL);
165: if (flg) {
166: /*
167: For linear problems with a time-dependent f(u,t) in the equation
168: u_t = f(u,t), the user provides the discretized right-hand-side
169: as a time-dependent matrix.
170: */
171: TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx);
172: TSSetRHSJacobian(ts, A, A, RHSMatrixHeat, &appctx);
173: } else {
174: /*
175: For linear problems with a time-independent f(u) in the equation
176: u_t = f(u), the user provides the discretized right-hand-side
177: as a matrix only once, and then sets the special Jacobian evaluation
178: routine TSComputeRHSJacobianConstant() which will NOT recompute the Jacobian.
179: */
180: RHSMatrixHeat(ts, 0.0, u, A, A, &appctx);
181: TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx);
182: TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &appctx);
183: }
184: } else {
185: Mat J;
187: RHSMatrixHeat(ts, 0.0, u, A, A, &appctx);
188: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &J);
189: TSSetIFunction(ts, NULL, IFunctionHeat, &appctx);
190: TSSetIJacobian(ts, J, J, IJacobianHeat, &appctx);
191: MatDestroy(&J);
193: PetscObjectReference((PetscObject)A);
194: appctx.A = A;
195: appctx.oshift = PETSC_MIN_REAL;
196: }
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Set solution vector and initial timestep
199: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: dt = appctx.h * appctx.h / 2.0;
202: TSSetTimeStep(ts, dt);
204: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: Customize timestepping solver:
206: - Set the solution method to be the Backward Euler method.
207: - Set timestepping duration info
208: Then set runtime options, which can override these defaults.
209: For example,
210: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
211: to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
212: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
214: TSSetMaxSteps(ts, time_steps_max);
215: TSSetMaxTime(ts, time_total_max);
216: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
217: TSSetFromOptions(ts);
219: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: Solve the problem
221: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
223: /*
224: Evaluate initial conditions
225: */
226: InitialConditions(u, &appctx);
228: /*
229: Run the timestepping solver
230: */
231: TSSolve(ts, u);
232: TSGetStepNumber(ts, &steps);
234: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: View timestepping solver info
236: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
238: 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));
239: if (!flg_string) {
240: TSView(ts, PETSC_VIEWER_STDOUT_SELF);
241: } else {
242: PetscViewer stringviewer;
243: char string[512];
244: const char *outstring;
246: PetscViewerStringOpen(PETSC_COMM_WORLD, string, sizeof(string), &stringviewer);
247: TSView(ts, stringviewer);
248: PetscViewerStringGetStringRead(stringviewer, &outstring, NULL);
250: PetscPrintf(PETSC_COMM_WORLD, "Output from string viewer:%s\n", outstring);
251: PetscViewerDestroy(&stringviewer);
252: }
254: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: Free work space. All PETSc objects should be destroyed when they
256: are no longer needed.
257: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
259: TSDestroy(&ts);
260: MatDestroy(&A);
261: VecDestroy(&u);
262: PetscViewerDestroy(&appctx.viewer1);
263: PetscViewerDestroy(&appctx.viewer2);
264: VecDestroy(&appctx.solution);
265: MatDestroy(&appctx.A);
267: /*
268: Always call PetscFinalize() before exiting a program. This routine
269: - finalizes the PETSc libraries as well as MPI
270: - provides summary and diagnostic information if certain runtime
271: options are chosen (e.g., -log_view).
272: */
273: PetscFinalize();
274: return 0;
275: }
276: /* --------------------------------------------------------------------- */
277: /*
278: InitialConditions - Computes the solution at the initial time.
280: Input Parameter:
281: u - uninitialized solution vector (global)
282: appctx - user-defined application context
284: Output Parameter:
285: u - vector with solution at initial time (global)
286: */
287: PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
288: {
289: PetscScalar *u_localptr, h = appctx->h;
290: PetscInt i;
292: /*
293: Get a pointer to vector data.
294: - For default PETSc vectors, VecGetArray() returns a pointer to
295: the data array. Otherwise, the routine is implementation dependent.
296: - You MUST call VecRestoreArray() when you no longer need access to
297: the array.
298: - Note that the Fortran interface to VecGetArray() differs from the
299: C version. See the users manual for details.
300: */
301: VecGetArrayWrite(u, &u_localptr);
303: /*
304: We initialize the solution array by simply writing the solution
305: directly into the array locations. Alternatively, we could use
306: VecSetValues() or VecSetValuesLocal().
307: */
308: for (i = 0; i < appctx->m; i++) u_localptr[i] = PetscSinScalar(PETSC_PI * i * 6. * h) + 3. * PetscSinScalar(PETSC_PI * i * 2. * h);
310: /*
311: Restore vector
312: */
313: VecRestoreArrayWrite(u, &u_localptr);
315: /*
316: Print debugging information if desired
317: */
318: if (appctx->debug) {
319: PetscPrintf(PETSC_COMM_WORLD, "Initial guess vector\n");
320: VecView(u, PETSC_VIEWER_STDOUT_SELF);
321: }
323: return 0;
324: }
325: /* --------------------------------------------------------------------- */
326: /*
327: ExactSolution - Computes the exact solution at a given time.
329: Input Parameters:
330: t - current time
331: solution - vector in which exact solution will be computed
332: appctx - user-defined application context
334: Output Parameter:
335: solution - vector with the newly computed exact solution
336: */
337: PetscErrorCode ExactSolution(PetscReal t, Vec solution, AppCtx *appctx)
338: {
339: PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2, tc = t;
340: PetscInt i;
342: /*
343: Get a pointer to vector data.
344: */
345: VecGetArrayWrite(solution, &s_localptr);
347: /*
348: Simply write the solution directly into the array locations.
349: Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
350: */
351: ex1 = PetscExpScalar(-36. * PETSC_PI * PETSC_PI * tc);
352: ex2 = PetscExpScalar(-4. * PETSC_PI * PETSC_PI * tc);
353: sc1 = PETSC_PI * 6. * h;
354: sc2 = PETSC_PI * 2. * h;
355: for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscSinScalar(sc1 * (PetscReal)i) * ex1 + 3. * PetscSinScalar(sc2 * (PetscReal)i) * ex2;
357: /*
358: Restore vector
359: */
360: VecRestoreArrayWrite(solution, &s_localptr);
361: return 0;
362: }
363: /* --------------------------------------------------------------------- */
364: /*
365: Monitor - User-provided routine to monitor the solution computed at
366: each timestep. This example plots the solution and computes the
367: error in two different norms.
369: This example also demonstrates changing the timestep via TSSetTimeStep().
371: Input Parameters:
372: ts - the timestep context
373: step - the count of the current step (with 0 meaning the
374: initial condition)
375: time - the current time
376: u - the solution at this timestep
377: ctx - the user-provided context for this monitoring routine.
378: In this case we use the application context which contains
379: information about the problem size, workspace and the exact
380: solution.
381: */
382: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx)
383: {
384: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
385: PetscReal norm_2, norm_max, dt, dttol;
387: /*
388: View a graph of the current iterate
389: */
390: VecView(u, appctx->viewer2);
392: /*
393: Compute the exact solution
394: */
395: ExactSolution(time, appctx->solution, appctx);
397: /*
398: Print debugging information if desired
399: */
400: if (appctx->debug) {
401: PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n");
402: VecView(u, PETSC_VIEWER_STDOUT_SELF);
403: PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n");
404: VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF);
405: }
407: /*
408: Compute the 2-norm and max-norm of the error
409: */
410: VecAXPY(appctx->solution, -1.0, u);
411: VecNorm(appctx->solution, NORM_2, &norm_2);
412: norm_2 = PetscSqrtReal(appctx->h) * norm_2;
413: VecNorm(appctx->solution, NORM_MAX, &norm_max);
415: TSGetTimeStep(ts, &dt);
416: PetscPrintf(PETSC_COMM_WORLD, "Timestep %3" PetscInt_FMT ": step size = %g, time = %g, 2-norm error = %g, max norm error = %g\n", step, (double)dt, (double)time, (double)norm_2, (double)norm_max);
418: appctx->norm_2 += norm_2;
419: appctx->norm_max += norm_max;
421: dttol = .0001;
422: PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, NULL);
423: if (dt < dttol) {
424: dt *= .999;
425: TSSetTimeStep(ts, dt);
426: }
428: /*
429: View a graph of the error
430: */
431: VecView(appctx->solution, appctx->viewer1);
433: /*
434: Print debugging information if desired
435: */
436: if (appctx->debug) {
437: PetscPrintf(PETSC_COMM_SELF, "Error vector\n");
438: VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF);
439: }
441: return 0;
442: }
443: /* --------------------------------------------------------------------- */
444: /*
445: RHSMatrixHeat - User-provided routine to compute the right-hand-side
446: matrix for the heat equation.
448: Input Parameters:
449: ts - the TS context
450: t - current time
451: global_in - global input vector
452: dummy - optional user-defined context, as set by TSetRHSJacobian()
454: Output Parameters:
455: AA - Jacobian matrix
456: BB - optionally different preconditioning matrix
457: str - flag indicating matrix structure
459: Notes:
460: Recall that MatSetValues() uses 0-based row and column numbers
461: in Fortran as well as in C.
462: */
463: PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)
464: {
465: Mat A = AA; /* Jacobian matrix */
466: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
467: PetscInt mstart = 0;
468: PetscInt mend = appctx->m;
469: PetscInt i, idx[3];
470: PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
472: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
473: Compute entries for the locally owned part of the matrix
474: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
475: /*
476: Set matrix rows corresponding to boundary data
477: */
479: mstart = 0;
480: v[0] = 1.0;
481: MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES);
482: mstart++;
484: mend--;
485: v[0] = 1.0;
486: MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES);
488: /*
489: Set matrix rows corresponding to interior data. We construct the
490: matrix one row at a time.
491: */
492: v[0] = sone;
493: v[1] = stwo;
494: v[2] = sone;
495: for (i = mstart; i < mend; i++) {
496: idx[0] = i - 1;
497: idx[1] = i;
498: idx[2] = i + 1;
499: MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES);
500: }
502: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
503: Complete the matrix assembly process and set some options
504: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
505: /*
506: Assemble matrix, using the 2-step process:
507: MatAssemblyBegin(), MatAssemblyEnd()
508: Computations can be done while messages are in transition
509: by placing code between these two statements.
510: */
511: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
512: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
514: /*
515: Set and option to indicate that we will never add a new nonzero location
516: to the matrix. If we do, it will generate an error.
517: */
518: MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
520: return 0;
521: }
523: PetscErrorCode IFunctionHeat(TS ts, PetscReal t, Vec X, Vec Xdot, Vec r, void *ctx)
524: {
525: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
527: MatMult(appctx->A, X, r);
528: VecAYPX(r, -1.0, Xdot);
529: return 0;
530: }
532: PetscErrorCode IJacobianHeat(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal s, Mat A, Mat B, void *ctx)
533: {
534: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
536: if (appctx->oshift == s) return 0;
537: MatCopy(appctx->A, A, SAME_NONZERO_PATTERN);
538: MatScale(A, -1);
539: MatShift(A, s);
540: MatCopy(A, B, SAME_NONZERO_PATTERN);
541: appctx->oshift = s;
542: return 0;
543: }
545: /*TEST
547: test:
548: args: -nox -ts_type ssp -ts_dt 0.0005
550: test:
551: suffix: 2
552: args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1
554: test:
555: suffix: 3
556: args: -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason
557: filter: sed "s/ATOL/RTOL/g"
558: requires: !single
560: test:
561: suffix: 4
562: args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason
563: filter: sed "s/ATOL/RTOL/g"
565: test:
566: suffix: 5
567: args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs
568: filter: sed "s/ATOL/RTOL/g"
570: test:
571: requires: !single
572: suffix: pod_guess
573: args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason
575: test:
576: requires: !single
577: suffix: pod_guess_Ainner
578: args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -ksp_guess_pod_Ainner -pc_type none -ksp_converged_reason
580: test:
581: requires: !single
582: suffix: fischer_guess
583: args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason
585: test:
586: requires: !single
587: suffix: fischer_guess_2
588: args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 2,10 -pc_type none -ksp_converged_reason
590: test:
591: requires: !single
592: suffix: fischer_guess_3
593: args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 3,10 -pc_type none -ksp_converged_reason
595: test:
596: requires: !single
597: suffix: stringview
598: args: -nox -ts_type rosw -test_string_viewer
600: test:
601: requires: !single
602: suffix: stringview_euler
603: args: -nox -ts_type euler -test_string_viewer
605: TEST*/