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";
7
8 /* ------------------------------------------------------------------------
9
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.
18
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>
25
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)
29
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)
35
36 The parallel version of this code is ts/tutorials/ex4.c
37
38 ------------------------------------------------------------------------- */
39
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>
51
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;
64
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 *);
72
main(int argc,char ** argv)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;
86
87 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88 Initialize program and set problem parameters
89 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90
91 PetscFunctionBeginUser;
92 PetscCall(PetscInitialize(&argc, &argv, NULL, 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!");
95
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;
103
104 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n"));
105
106 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107 Create vector data structures
108 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109
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));
115
116 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117 Set up displays to show graphs of the solution and error
118 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119
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));
126
127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 Create timestepping solver context
129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130
131 PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
132 PetscCall(TSSetProblemType(ts, TS_LINEAR));
133
134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135 Set optional user-defined monitoring routine
136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137
138 PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
139
140 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141
142 Create matrix data structure; set matrix evaluation routine.
143 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144
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));
149
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 }
170
171 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172 Set solution vector and initial timestep
173 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174
175 dt = appctx.h * appctx.h / 2.0;
176 PetscCall(TSSetTimeStep(ts, dt));
177 PetscCall(TSSetSolution(ts, u));
178
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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188
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));
193
194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195 Solve the problem
196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197
198 /*
199 Evaluate initial conditions
200 */
201 PetscCall(InitialConditions(u, &appctx));
202
203 /*
204 Run the timestepping solver
205 */
206 PetscCall(TSSolve(ts, u));
207 PetscCall(TSGetSolveTime(ts, &ftime));
208 PetscCall(TSGetStepNumber(ts, &steps));
209
210 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211 View timestepping solver info
212 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213
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));
216
217 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218 Free work space. All PETSc objects should be destroyed when they
219 are no longer needed.
220 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221
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));
228
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.
241
242 Input Parameter:
243 u - uninitialized solution vector (global)
244 appctx - user-defined application context
245
246 Output Parameter:
247 u - vector with solution at initial time (global)
248 */
InitialConditions(Vec u,AppCtx * appctx)249 PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
250 {
251 PetscScalar *u_localptr, h = appctx->h;
252 PetscInt i;
253
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));
265
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);
272
273 /*
274 Restore vector
275 */
276 PetscCall(VecRestoreArray(u, &u_localptr));
277
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 }
285 PetscFunctionReturn(PETSC_SUCCESS);
286 }
287 /* --------------------------------------------------------------------- */
288 /*
289 ExactSolution - Computes the exact solution at a given time.
290
291 Input Parameters:
292 t - current time
293 solution - vector in which exact solution will be computed
294 appctx - user-defined application context
295
296 Output Parameter:
297 solution - vector with the newly computed exact solution
298 */
ExactSolution(PetscReal t,Vec solution,AppCtx * appctx)299 PetscErrorCode ExactSolution(PetscReal t, Vec solution, AppCtx *appctx)
300 {
301 PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2, tc = t;
302 PetscInt i;
303
304 PetscFunctionBeginUser;
305 /*
306 Get a pointer to vector data.
307 */
308 PetscCall(VecGetArray(solution, &s_localptr));
309
310 /*
311 Simply write the solution directly into the array locations.
312 Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
313 */
314 ex1 = PetscExpScalar(-36. * PETSC_PI * PETSC_PI * tc);
315 ex2 = PetscExpScalar(-4. * PETSC_PI * PETSC_PI * tc);
316 sc1 = PETSC_PI * 6. * h;
317 sc2 = PETSC_PI * 2. * h;
318 for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscCosScalar(sc1 * (PetscReal)i) * ex1 + 3. * PetscCosScalar(sc2 * (PetscReal)i) * ex2;
319
320 /*
321 Restore vector
322 */
323 PetscCall(VecRestoreArray(solution, &s_localptr));
324 PetscFunctionReturn(PETSC_SUCCESS);
325 }
326 /* --------------------------------------------------------------------- */
327 /*
328 Monitor - User-provided routine to monitor the solution computed at
329 each timestep. This example plots the solution and computes the
330 error in two different norms.
331
332 Input Parameters:
333 ts - the timestep context
334 step - the count of the current step (with 0 meaning the
335 initial condition)
336 time - the current time
337 u - the solution at this timestep
338 ctx - the user-provided context for this monitoring routine.
339 In this case we use the application context which contains
340 information about the problem size, workspace and the exact
341 solution.
342 */
Monitor(TS ts,PetscInt step,PetscReal time,Vec u,PetscCtx ctx)343 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, PetscCtx ctx)
344 {
345 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
346 PetscReal norm_2, norm_max;
347
348 PetscFunctionBeginUser;
349 /*
350 View a graph of the current iterate
351 */
352 PetscCall(VecView(u, appctx->viewer2));
353
354 /*
355 Compute the exact solution
356 */
357 PetscCall(ExactSolution(time, appctx->solution, appctx));
358
359 /*
360 Print debugging information if desired
361 */
362 if (appctx->debug) {
363 printf("Computed solution vector\n");
364 PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
365 printf("Exact solution vector\n");
366 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
367 }
368
369 /*
370 Compute the 2-norm and max-norm of the error
371 */
372 PetscCall(VecAXPY(appctx->solution, -1.0, u));
373 PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
374 norm_2 = PetscSqrtReal(appctx->h) * norm_2;
375 PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
376 if (norm_2 < 1e-14) norm_2 = 0;
377 if (norm_max < 1e-14) norm_max = 0;
378
379 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));
380 appctx->norm_2 += norm_2;
381 appctx->norm_max += norm_max;
382
383 /*
384 View a graph of the error
385 */
386 PetscCall(VecView(appctx->solution, appctx->viewer1));
387
388 /*
389 Print debugging information if desired
390 */
391 if (appctx->debug) {
392 printf("Error vector\n");
393 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
394 }
395 PetscFunctionReturn(PETSC_SUCCESS);
396 }
397 /* --------------------------------------------------------------------- */
398 /*
399 RHSMatrixHeat - User-provided routine to compute the right-hand-side
400 matrix for the heat equation.
401
402 Input Parameters:
403 ts - the TS context
404 t - current time
405 global_in - global input vector
406 dummy - optional user-defined context, as set by TSetRHSJacobian()
407
408 Output Parameters:
409 AA - Jacobian matrix
410 BB - optionally different matrix used to construct the preconditioner
411
412 Notes:
413 Recall that MatSetValues() uses 0-based row and column numbers
414 in Fortran as well as in C.
415 */
RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,PetscCtx ctx)416 PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, PetscCtx ctx)
417 {
418 Mat A = AA; /* Jacobian matrix */
419 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
420 PetscInt mstart = 0;
421 PetscInt mend = appctx->m;
422 PetscInt i, idx[3];
423 PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
424
425 PetscFunctionBeginUser;
426 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
427 Compute entries for the locally owned part of the matrix
428 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
429 /*
430 Set matrix rows corresponding to boundary data
431 */
432
433 mstart = 0;
434 v[0] = 1.0;
435 PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
436 mstart++;
437
438 mend--;
439 v[0] = 1.0;
440 PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
441
442 /*
443 Set matrix rows corresponding to interior data. We construct the
444 matrix one row at a time.
445 */
446 v[0] = sone;
447 v[1] = stwo;
448 v[2] = sone;
449 for (i = mstart; i < mend; i++) {
450 idx[0] = i - 1;
451 idx[1] = i;
452 idx[2] = i + 1;
453 PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
454 }
455
456 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
457 Complete the matrix assembly process and set some options
458 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
459 /*
460 Assemble matrix, using the 2-step process:
461 MatAssemblyBegin(), MatAssemblyEnd()
462 Computations can be done while messages are in transition
463 by placing code between these two statements.
464 */
465 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
466 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
467
468 /*
469 Set and option to indicate that we will never add a new nonzero location
470 to the matrix. If we do, it will generate an error.
471 */
472 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
473 PetscFunctionReturn(PETSC_SUCCESS);
474 }
475
476 /*TEST
477
478 test:
479 requires: x
480
481 test:
482 suffix: nox
483 args: -nox
484
485 TEST*/
486