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) = 0, u(t,1) = 0,
15 and the initial condition
16 u(0,x) = sin(6*pi*x) + 3*sin(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) * sin(6*pi*x) +
28 3*exp(-4*pi*pi*t) * sin(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 f(u)
35
36 The parallel version of this code is ts/tutorials/ex4.c
37
38 ------------------------------------------------------------------------- */
39
40 /*
41 Include "ts.h" so that we can use TS solvers. Note that this file
42 automatically includes:
43 petscsys.h - base PETSc routines vec.h - vectors
44 sys.h - system routines mat.h - matrices
45 is.h - index sets ksp.h - Krylov subspace methods
46 viewer.h - viewers pc.h - preconditioners
47 snes.h - nonlinear solvers
48 */
49
50 #include <petscts.h>
51 #include <petscdraw.h>
52
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;
65
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 *);
73 extern PetscErrorCode MyBCRoutine(TS, PetscReal, Vec, void *);
74
main(int argc,char ** argv)75 int main(int argc, char **argv)
76 {
77 AppCtx appctx; /* user-defined application context */
78 TS ts; /* timestepping context */
79 Mat A; /* matrix data structure */
80 Vec u; /* approximate solution vector */
81 PetscReal time_total_max = 100.0; /* default max total time */
82 PetscInt time_steps_max = 100; /* default max timesteps */
83 PetscDraw draw; /* drawing context */
84 PetscInt steps, m;
85 PetscMPIInt size;
86 PetscReal dt;
87 PetscReal ftime;
88 PetscBool flg;
89 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90 Initialize program and set problem parameters
91 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92
93 PetscFunctionBeginUser;
94 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
95 MPI_Comm_size(PETSC_COMM_WORLD, &size);
96 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
97
98 m = 60;
99 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
100 PetscCall(PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug));
101
102 appctx.m = m;
103 appctx.h = 1.0 / (m - 1.0);
104 appctx.norm_2 = 0.0;
105 appctx.norm_max = 0.0;
106
107 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n"));
108
109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110 Create vector data structures
111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112
113 /*
114 Create vector data structures for approximate and exact solutions
115 */
116 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &u));
117 PetscCall(VecDuplicate(u, &appctx.solution));
118
119 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120 Set up displays to show graphs of the solution and error
121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122
123 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 380, 400, 160, &appctx.viewer1));
124 PetscCall(PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw));
125 PetscCall(PetscDrawSetDoubleBuffer(draw));
126 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 0, 400, 160, &appctx.viewer2));
127 PetscCall(PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw));
128 PetscCall(PetscDrawSetDoubleBuffer(draw));
129
130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131 Create timestepping solver context
132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133
134 PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
135 PetscCall(TSSetProblemType(ts, TS_LINEAR));
136
137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 Set optional user-defined monitoring routine
139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140
141 PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
142
143 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144
145 Create matrix data structure; set matrix evaluation routine.
146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147
148 PetscCall(MatCreate(PETSC_COMM_SELF, &A));
149 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
150 PetscCall(MatSetFromOptions(A));
151 PetscCall(MatSetUp(A));
152
153 PetscCall(PetscOptionsHasName(NULL, NULL, "-time_dependent_rhs", &flg));
154 if (flg) {
155 /*
156 For linear problems with a time-dependent f(u,t) in the equation
157 u_t = f(u,t), the user provides the discretized right-hand-side
158 as a time-dependent matrix.
159 */
160 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
161 PetscCall(TSSetRHSJacobian(ts, A, A, RHSMatrixHeat, &appctx));
162 } else {
163 /*
164 For linear problems with a time-independent f(u) in the equation
165 u_t = f(u), the user provides the discretized right-hand side
166 as a matrix only once, and then sets a null matrix evaluation
167 routine.
168 */
169 PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx));
170 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
171 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &appctx));
172 }
173
174 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175 Set solution vector and initial timestep
176 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177
178 dt = appctx.h * appctx.h / 2.0;
179 PetscCall(TSSetTimeStep(ts, dt));
180 PetscCall(TSSetSolution(ts, u));
181
182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183 Customize timestepping solver:
184 - Set the solution method to be the Backward Euler method.
185 - Set timestepping duration info
186 Then set runtime options, which can override these defaults.
187 For example,
188 -ts_max_steps <maxsteps> -ts_max_time <maxtime>
189 to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
190 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191
192 PetscCall(TSSetMaxSteps(ts, time_steps_max));
193 PetscCall(TSSetMaxTime(ts, time_total_max));
194 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
195 PetscCall(TSSetFromOptions(ts));
196
197 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198 Solve the problem
199 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200
201 /*
202 Evaluate initial conditions
203 */
204 PetscCall(InitialConditions(u, &appctx));
205
206 /*
207 Run the timestepping solver
208 */
209 PetscCall(TSSolve(ts, u));
210 PetscCall(TSGetSolveTime(ts, &ftime));
211 PetscCall(TSGetStepNumber(ts, &steps));
212
213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214 View timestepping solver info
215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216
217 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)));
218 PetscCall(TSView(ts, PETSC_VIEWER_STDOUT_SELF));
219
220 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221 Free work space. All PETSc objects should be destroyed when they
222 are no longer needed.
223 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224
225 PetscCall(TSDestroy(&ts));
226 PetscCall(MatDestroy(&A));
227 PetscCall(VecDestroy(&u));
228 PetscCall(PetscViewerDestroy(&appctx.viewer1));
229 PetscCall(PetscViewerDestroy(&appctx.viewer2));
230 PetscCall(VecDestroy(&appctx.solution));
231
232 /*
233 Always call PetscFinalize() before exiting a program. This routine
234 - finalizes the PETSc libraries as well as MPI
235 - provides summary and diagnostic information if certain runtime
236 options are chosen (e.g., -log_view).
237 */
238 PetscCall(PetscFinalize());
239 return 0;
240 }
241 /* --------------------------------------------------------------------- */
242 /*
243 InitialConditions - Computes the solution at the initial time.
244
245 Input Parameter:
246 u - uninitialized solution vector (global)
247 appctx - user-defined application context
248
249 Output Parameter:
250 u - vector with solution at initial time (global)
251 */
InitialConditions(Vec u,AppCtx * appctx)252 PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
253 {
254 PetscScalar *u_localptr;
255 PetscInt i;
256
257 PetscFunctionBeginUser;
258 /*
259 Get a pointer to vector data.
260 - For default PETSc vectors, VecGetArray() returns a pointer to
261 the data array. Otherwise, the routine is implementation dependent.
262 - You MUST call VecRestoreArray() when you no longer need access to
263 the array.
264 - Note that the Fortran interface to VecGetArray() differs from the
265 C version. See the users manual for details.
266 */
267 PetscCall(VecGetArray(u, &u_localptr));
268
269 /*
270 We initialize the solution array by simply writing the solution
271 directly into the array locations. Alternatively, we could use
272 VecSetValues() or VecSetValuesLocal().
273 */
274 for (i = 0; i < appctx->m; i++) u_localptr[i] = PetscSinReal(PETSC_PI * i * 6. * appctx->h) + 3. * PetscSinReal(PETSC_PI * i * 2. * appctx->h);
275
276 /*
277 Restore vector
278 */
279 PetscCall(VecRestoreArray(u, &u_localptr));
280
281 /*
282 Print debugging information if desired
283 */
284 if (appctx->debug) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
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;
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 = PetscExpReal(-36. * PETSC_PI * PETSC_PI * t);
315 ex2 = PetscExpReal(-4. * PETSC_PI * PETSC_PI * t);
316 sc1 = PETSC_PI * 6. * h;
317 sc2 = PETSC_PI * 2. * h;
318 for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscSinReal(PetscRealPart(sc1) * (PetscReal)i) * ex1 + 3. * PetscSinReal(PetscRealPart(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 This example also demonstrates changing the timestep via TSSetTimeStep().
333
334 Input Parameters:
335 ts - the timestep context
336 step - the count of the current step (with 0 meaning the
337 initial condition)
338 crtime - 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 */
Monitor(TS ts,PetscInt step,PetscReal crtime,Vec u,PetscCtx ctx)345 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, PetscCtx ctx)
346 {
347 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
348 PetscReal norm_2, norm_max, dt, dttol;
349 PetscBool flg;
350
351 PetscFunctionBeginUser;
352 /*
353 View a graph of the current iterate
354 */
355 PetscCall(VecView(u, appctx->viewer2));
356
357 /*
358 Compute the exact solution
359 */
360 PetscCall(ExactSolution(crtime, appctx->solution, appctx));
361
362 /*
363 Print debugging information if desired
364 */
365 if (appctx->debug) {
366 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n"));
367 PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
368 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
369 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
370 }
371
372 /*
373 Compute the 2-norm and max-norm of the error
374 */
375 PetscCall(VecAXPY(appctx->solution, -1.0, u));
376 PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
377 norm_2 = PetscSqrtReal(appctx->h) * norm_2;
378 PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
379
380 PetscCall(TSGetTimeStep(ts, &dt));
381 if (norm_2 > 1.e-2) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Timestep %" PetscInt_FMT ": step size = %g, time = %g, 2-norm error = %g, max norm error = %g\n", step, (double)dt, (double)crtime, (double)norm_2, (double)norm_max));
382 appctx->norm_2 += norm_2;
383 appctx->norm_max += norm_max;
384
385 dttol = .0001;
386 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, &flg));
387 if (dt < dttol) {
388 dt *= .999;
389 PetscCall(TSSetTimeStep(ts, dt));
390 }
391
392 /*
393 View a graph of the error
394 */
395 PetscCall(VecView(appctx->solution, appctx->viewer1));
396
397 /*
398 Print debugging information if desired
399 */
400 if (appctx->debug) {
401 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
402 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
403 }
404 PetscFunctionReturn(PETSC_SUCCESS);
405 }
406 /* --------------------------------------------------------------------- */
407 /*
408 RHSMatrixHeat - User-provided routine to compute the right-hand-side
409 matrix for the heat equation.
410
411 Input Parameters:
412 ts - the TS context
413 t - current time
414 global_in - global input vector
415 dummy - optional user-defined context, as set by TSetRHSJacobian()
416
417 Output Parameters:
418 AA - Jacobian matrix
419 BB - optionally different matrix used to construct the preconditioner
420
421 Notes:
422 Recall that MatSetValues() uses 0-based row and column numbers
423 in Fortran as well as in C.
424 */
RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,PetscCtx ctx)425 PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, PetscCtx ctx)
426 {
427 Mat A = AA; /* Jacobian matrix */
428 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
429 PetscInt mstart = 0;
430 PetscInt mend = appctx->m;
431 PetscInt i, idx[3];
432 PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
433
434 PetscFunctionBeginUser;
435 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436 Compute entries for the locally owned part of the matrix
437 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
438 /*
439 Set matrix rows corresponding to boundary data
440 */
441
442 mstart = 0;
443 v[0] = 1.0;
444 PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
445 mstart++;
446
447 mend--;
448 v[0] = 1.0;
449 PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
450
451 /*
452 Set matrix rows corresponding to interior data. We construct the
453 matrix one row at a time.
454 */
455 v[0] = sone;
456 v[1] = stwo;
457 v[2] = sone;
458 for (i = mstart; i < mend; i++) {
459 idx[0] = i - 1;
460 idx[1] = i;
461 idx[2] = i + 1;
462 PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
463 }
464
465 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
466 Complete the matrix assembly process and set some options
467 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
468 /*
469 Assemble matrix, using the 2-step process:
470 MatAssemblyBegin(), MatAssemblyEnd()
471 Computations can be done while messages are in transition
472 by placing code between these two statements.
473 */
474 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
475 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
476
477 /*
478 Set and option to indicate that we will never add a new nonzero location
479 to the matrix. If we do, it will generate an error.
480 */
481 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
482 PetscFunctionReturn(PETSC_SUCCESS);
483 }
484 /* --------------------------------------------------------------------- */
485 /*
486 Input Parameters:
487 ts - the TS context
488 t - current time
489 f - function
490 ctx - optional user-defined context, as set by TSetBCFunction()
491 */
MyBCRoutine(TS ts,PetscReal t,Vec f,PetscCtx ctx)492 PetscErrorCode MyBCRoutine(TS ts, PetscReal t, Vec f, PetscCtx ctx)
493 {
494 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
495 PetscInt m = appctx->m;
496 PetscScalar *fa;
497
498 PetscFunctionBeginUser;
499 PetscCall(VecGetArray(f, &fa));
500 fa[0] = 0.0;
501 fa[m - 1] = 1.0;
502 PetscCall(VecRestoreArray(f, &fa));
503 PetscCall(PetscPrintf(PETSC_COMM_SELF, "t=%g\n", (double)t));
504 PetscFunctionReturn(PETSC_SUCCESS);
505 }
506
507 /*TEST
508
509 test:
510 args: -nox -ts_max_steps 4
511
512 TEST*/
513