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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscInitialize(&argc,&argv,(char*)0,help);
99: MPI_Comm_size(PETSC_COMM_WORLD,&size);
102: m = 60;
103: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
104: PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug);
105: flg_string = PETSC_FALSE;
106: PetscOptionsGetBool(NULL,NULL,"-test_string_viewer",&flg_string,NULL);
108: appctx.m = m;
109: appctx.h = 1.0/(m-1.0);
110: appctx.norm_2 = 0.0;
111: appctx.norm_max = 0.0;
113: PetscPrintf(PETSC_COMM_SELF,"Solving a linear TS problem on 1 processor\n");
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Create vector data structures
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: /*
120: Create vector data structures for approximate and exact solutions
121: */
122: VecCreateSeq(PETSC_COMM_SELF,m,&u);
123: VecDuplicate(u,&appctx.solution);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Set up displays to show graphs of the solution and error
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,380,400,160,&appctx.viewer1);
130: PetscViewerDrawGetDraw(appctx.viewer1,0,&draw);
131: PetscDrawSetDoubleBuffer(draw);
132: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,0,400,160,&appctx.viewer2);
133: PetscViewerDrawGetDraw(appctx.viewer2,0,&draw);
134: PetscDrawSetDoubleBuffer(draw);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Create timestepping solver context
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: TSCreate(PETSC_COMM_SELF,&ts);
141: TSSetProblemType(ts,TS_LINEAR);
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Set optional user-defined monitoring routine
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: if (!flg_string) {
148: TSMonitorSet(ts,Monitor,&appctx,NULL);
149: }
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Create matrix data structure; set matrix evaluation routine.
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: MatCreate(PETSC_COMM_SELF,&A);
157: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m);
158: MatSetFromOptions(A);
159: MatSetUp(A);
161: flg = PETSC_FALSE;
162: PetscOptionsGetBool(NULL,NULL,"-use_ifunc",&flg,NULL);
163: if (!flg) {
164: appctx.A = NULL;
165: PetscOptionsGetBool(NULL,NULL,"-time_dependent_rhs",&flg,NULL);
166: if (flg) {
167: /*
168: For linear problems with a time-dependent f(u,t) in the equation
169: u_t = f(u,t), the user provides the discretized right-hand-side
170: as a time-dependent matrix.
171: */
172: TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);
173: TSSetRHSJacobian(ts,A,A,RHSMatrixHeat,&appctx);
174: } else {
175: /*
176: For linear problems with a time-independent f(u) in the equation
177: u_t = f(u), the user provides the discretized right-hand-side
178: as a matrix only once, and then sets the special Jacobian evaluation
179: routine TSComputeRHSJacobianConstant() which will NOT recompute the Jacobian.
180: */
181: RHSMatrixHeat(ts,0.0,u,A,A,&appctx);
182: TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);
183: TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx);
184: }
185: } else {
186: Mat J;
188: RHSMatrixHeat(ts,0.0,u,A,A,&appctx);
189: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&J);
190: TSSetIFunction(ts,NULL,IFunctionHeat,&appctx);
191: TSSetIJacobian(ts,J,J,IJacobianHeat,&appctx);
192: MatDestroy(&J);
194: PetscObjectReference((PetscObject)A);
195: appctx.A = A;
196: appctx.oshift = PETSC_MIN_REAL;
197: }
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Set solution vector and initial timestep
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
202: dt = appctx.h*appctx.h/2.0;
203: TSSetTimeStep(ts,dt);
205: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206: Customize timestepping solver:
207: - Set the solution method to be the Backward Euler method.
208: - Set timestepping duration info
209: Then set runtime options, which can override these defaults.
210: For example,
211: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
212: to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
213: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
215: TSSetMaxSteps(ts,time_steps_max);
216: TSSetMaxTime(ts,time_total_max);
217: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
218: TSSetFromOptions(ts);
220: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221: Solve the problem
222: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224: /*
225: Evaluate initial conditions
226: */
227: InitialConditions(u,&appctx);
229: /*
230: Run the timestepping solver
231: */
232: TSSolve(ts,u);
233: TSGetStepNumber(ts,&steps);
235: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: View timestepping solver info
237: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
239: 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));
240: if (!flg_string) {
241: TSView(ts,PETSC_VIEWER_STDOUT_SELF);
242: } else {
243: PetscViewer stringviewer;
244: char string[512];
245: const char *outstring;
247: PetscViewerStringOpen(PETSC_COMM_WORLD,string,sizeof(string),&stringviewer);
248: TSView(ts,stringviewer);
249: PetscViewerStringGetStringRead(stringviewer,&outstring,NULL);
251: PetscPrintf(PETSC_COMM_WORLD,"Output from string viewer:%s\n",outstring);
252: PetscViewerDestroy(&stringviewer);
253: }
255: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
256: Free work space. All PETSc objects should be destroyed when they
257: are no longer needed.
258: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260: TSDestroy(&ts);
261: MatDestroy(&A);
262: VecDestroy(&u);
263: PetscViewerDestroy(&appctx.viewer1);
264: PetscViewerDestroy(&appctx.viewer2);
265: VecDestroy(&appctx.solution);
266: MatDestroy(&appctx.A);
268: /*
269: Always call PetscFinalize() before exiting a program. This routine
270: - finalizes the PETSc libraries as well as MPI
271: - provides summary and diagnostic information if certain runtime
272: options are chosen (e.g., -log_view).
273: */
274: PetscFinalize();
275: return 0;
276: }
277: /* --------------------------------------------------------------------- */
278: /*
279: InitialConditions - Computes the solution at the initial time.
281: Input Parameter:
282: u - uninitialized solution vector (global)
283: appctx - user-defined application context
285: Output Parameter:
286: u - vector with solution at initial time (global)
287: */
288: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
289: {
290: PetscScalar *u_localptr,h = appctx->h;
291: PetscInt i;
293: /*
294: Get a pointer to vector data.
295: - For default PETSc vectors, VecGetArray() returns a pointer to
296: the data array. Otherwise, the routine is implementation dependent.
297: - You MUST call VecRestoreArray() when you no longer need access to
298: the array.
299: - Note that the Fortran interface to VecGetArray() differs from the
300: C version. See the users manual for details.
301: */
302: VecGetArrayWrite(u,&u_localptr);
304: /*
305: We initialize the solution array by simply writing the solution
306: directly into the array locations. Alternatively, we could use
307: VecSetValues() or VecSetValuesLocal().
308: */
309: for (i=0; i<appctx->m; i++) u_localptr[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
311: /*
312: Restore vector
313: */
314: VecRestoreArrayWrite(u,&u_localptr);
316: /*
317: Print debugging information if desired
318: */
319: if (appctx->debug) {
320: PetscPrintf(PETSC_COMM_WORLD,"Initial guess vector\n");
321: VecView(u,PETSC_VIEWER_STDOUT_SELF);
322: }
324: return 0;
325: }
326: /* --------------------------------------------------------------------- */
327: /*
328: ExactSolution - Computes the exact solution at a given time.
330: Input Parameters:
331: t - current time
332: solution - vector in which exact solution will be computed
333: appctx - user-defined application context
335: Output Parameter:
336: solution - vector with the newly computed exact solution
337: */
338: PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
339: {
340: PetscScalar *s_localptr,h = appctx->h,ex1,ex2,sc1,sc2,tc = t;
341: PetscInt i;
343: /*
344: Get a pointer to vector data.
345: */
346: VecGetArrayWrite(solution,&s_localptr);
348: /*
349: Simply write the solution directly into the array locations.
350: Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
351: */
352: ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*tc);
353: ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*tc);
354: sc1 = PETSC_PI*6.*h; 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 %3D: 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; v[1] = stwo; v[2] = sone;
493: for (i=mstart; i<mend; i++) {
494: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
495: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
496: }
498: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
499: Complete the matrix assembly process and set some options
500: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
501: /*
502: Assemble matrix, using the 2-step process:
503: MatAssemblyBegin(), MatAssemblyEnd()
504: Computations can be done while messages are in transition
505: by placing code between these two statements.
506: */
507: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
508: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
510: /*
511: Set and option to indicate that we will never add a new nonzero location
512: to the matrix. If we do, it will generate an error.
513: */
514: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
516: return 0;
517: }
519: PetscErrorCode IFunctionHeat(TS ts,PetscReal t,Vec X,Vec Xdot,Vec r,void *ctx)
520: {
521: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
523: MatMult(appctx->A,X,r);
524: VecAYPX(r,-1.0,Xdot);
525: return 0;
526: }
528: PetscErrorCode IJacobianHeat(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal s,Mat A,Mat B,void *ctx)
529: {
530: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
532: if (appctx->oshift == s) return 0;
533: MatCopy(appctx->A,A,SAME_NONZERO_PATTERN);
534: MatScale(A,-1);
535: MatShift(A,s);
536: MatCopy(A,B,SAME_NONZERO_PATTERN);
537: appctx->oshift = s;
538: return 0;
539: }
541: /*TEST
543: test:
544: args: -nox -ts_type ssp -ts_dt 0.0005
546: test:
547: suffix: 2
548: args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1
550: test:
551: suffix: 3
552: args: -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason
553: filter: sed "s/ATOL/RTOL/g"
554: requires: !single
556: test:
557: suffix: 4
558: args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason
559: filter: sed "s/ATOL/RTOL/g"
561: test:
562: suffix: 5
563: args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs
564: filter: sed "s/ATOL/RTOL/g"
566: test:
567: requires: !single
568: suffix: pod_guess
569: args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason
571: test:
572: requires: !single
573: suffix: pod_guess_Ainner
574: 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
576: test:
577: requires: !single
578: suffix: fischer_guess
579: args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason
581: test:
582: requires: !single
583: suffix: fischer_guess_2
584: 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
586: test:
587: requires: !single
588: suffix: fischer_guess_3
589: 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
591: test:
592: requires: !single
593: suffix: stringview
594: args: -nox -ts_type rosw -test_string_viewer
596: test:
597: requires: !single
598: suffix: stringview_euler
599: args: -nox -ts_type euler -test_string_viewer
601: TEST*/