1 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"; 8 9 /* ------------------------------------------------------------------------ 10 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) = 0, u(t,1) = 0, 16 and the initial condition 17 u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x). 18 This is a linear, second-order, parabolic equation. 19 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> 26 27 We compare the approximate solution with the exact solution, given by 28 u_exact(x,t) = exp(-36*pi*pi*t) * sin(6*pi*x) + 29 3*exp(-4*pi*pi*t) * sin(2*pi*x) 30 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 f(u) 36 37 The parallel version of this code is ts/tutorials/ex4.c 38 39 ------------------------------------------------------------------------- */ 40 41 /* 42 Include "ts.h" so that we can use TS solvers. Note that this file 43 automatically includes: 44 petscsys.h - base PETSc routines vec.h - vectors 45 sys.h - system routines mat.h - matrices 46 is.h - index sets ksp.h - Krylov subspace methods 47 viewer.h - viewers pc.h - preconditioners 48 snes.h - nonlinear solvers 49 */ 50 51 #include <petscts.h> 52 #include <petscdraw.h> 53 54 /* 55 User-defined application context - contains data needed by the 56 application-provided call-back routines. 57 */ 58 typedef struct { 59 Vec solution; /* global exact solution vector */ 60 PetscInt m; /* total number of grid points */ 61 PetscReal h; /* mesh width h = 1/(m-1) */ 62 PetscBool debug; /* flag (1 indicates activation of debugging printouts) */ 63 PetscViewer viewer1, viewer2; /* viewers for the solution and error */ 64 PetscReal norm_2, norm_max; /* error norms */ 65 } AppCtx; 66 67 /* 68 User-defined routines 69 */ 70 extern PetscErrorCode InitialConditions(Vec, AppCtx *); 71 extern PetscErrorCode RHSMatrixHeat(TS, PetscReal, Vec, Mat, Mat, void *); 72 extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *); 73 extern PetscErrorCode ExactSolution(PetscReal, Vec, AppCtx *); 74 extern PetscErrorCode MyBCRoutine(TS, PetscReal, Vec, void *); 75 76 int main(int argc, char **argv) 77 { 78 AppCtx appctx; /* user-defined application context */ 79 TS ts; /* timestepping context */ 80 Mat A; /* matrix data structure */ 81 Vec u; /* approximate solution vector */ 82 PetscReal time_total_max = 100.0; /* default max total time */ 83 PetscInt time_steps_max = 100; /* default max timesteps */ 84 PetscDraw draw; /* drawing context */ 85 PetscInt steps, m; 86 PetscMPIInt size; 87 PetscReal dt; 88 PetscReal ftime; 89 PetscBool flg; 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Initialize program and set problem parameters 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 93 94 PetscFunctionBeginUser; 95 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 96 MPI_Comm_size(PETSC_COMM_WORLD, &size); 97 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 98 99 m = 60; 100 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 101 PetscCall(PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug)); 102 103 appctx.m = m; 104 appctx.h = 1.0 / (m - 1.0); 105 appctx.norm_2 = 0.0; 106 appctx.norm_max = 0.0; 107 108 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n")); 109 110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 111 Create vector data structures 112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 113 114 /* 115 Create vector data structures for approximate and exact solutions 116 */ 117 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &u)); 118 PetscCall(VecDuplicate(u, &appctx.solution)); 119 120 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 121 Set up displays to show graphs of the solution and error 122 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 123 124 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 380, 400, 160, &appctx.viewer1)); 125 PetscCall(PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw)); 126 PetscCall(PetscDrawSetDoubleBuffer(draw)); 127 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 0, 400, 160, &appctx.viewer2)); 128 PetscCall(PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw)); 129 PetscCall(PetscDrawSetDoubleBuffer(draw)); 130 131 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 132 Create timestepping solver context 133 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 134 135 PetscCall(TSCreate(PETSC_COMM_SELF, &ts)); 136 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 137 138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139 Set optional user-defined monitoring routine 140 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 141 142 PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL)); 143 144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 146 Create matrix data structure; set matrix evaluation routine. 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 149 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 150 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m)); 151 PetscCall(MatSetFromOptions(A)); 152 PetscCall(MatSetUp(A)); 153 154 PetscCall(PetscOptionsHasName(NULL, NULL, "-time_dependent_rhs", &flg)); 155 if (flg) { 156 /* 157 For linear problems with a time-dependent f(u,t) in the equation 158 u_t = f(u,t), the user provides the discretized right-hand-side 159 as a time-dependent matrix. 160 */ 161 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx)); 162 PetscCall(TSSetRHSJacobian(ts, A, A, RHSMatrixHeat, &appctx)); 163 } else { 164 /* 165 For linear problems with a time-independent f(u) in the equation 166 u_t = f(u), the user provides the discretized right-hand-side 167 as a matrix only once, and then sets a null matrix evaluation 168 routine. 169 */ 170 PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx)); 171 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx)); 172 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &appctx)); 173 } 174 175 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 176 Set solution vector and initial timestep 177 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 178 179 dt = appctx.h * appctx.h / 2.0; 180 PetscCall(TSSetTimeStep(ts, dt)); 181 PetscCall(TSSetSolution(ts, u)); 182 183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 Customize timestepping solver: 185 - Set the solution method to be the Backward Euler method. 186 - Set timestepping duration info 187 Then set runtime options, which can override these defaults. 188 For example, 189 -ts_max_steps <maxsteps> -ts_max_time <maxtime> 190 to override the defaults set by TSSetMaxSteps()/TSSetMaxTime(). 191 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 192 193 PetscCall(TSSetMaxSteps(ts, time_steps_max)); 194 PetscCall(TSSetMaxTime(ts, time_total_max)); 195 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 196 PetscCall(TSSetFromOptions(ts)); 197 198 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 199 Solve the problem 200 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 201 202 /* 203 Evaluate initial conditions 204 */ 205 PetscCall(InitialConditions(u, &appctx)); 206 207 /* 208 Run the timestepping solver 209 */ 210 PetscCall(TSSolve(ts, u)); 211 PetscCall(TSGetSolveTime(ts, &ftime)); 212 PetscCall(TSGetStepNumber(ts, &steps)); 213 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 View timestepping solver info 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 218 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))); 219 PetscCall(TSView(ts, PETSC_VIEWER_STDOUT_SELF)); 220 221 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 222 Free work space. All PETSc objects should be destroyed when they 223 are no longer needed. 224 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 225 226 PetscCall(TSDestroy(&ts)); 227 PetscCall(MatDestroy(&A)); 228 PetscCall(VecDestroy(&u)); 229 PetscCall(PetscViewerDestroy(&appctx.viewer1)); 230 PetscCall(PetscViewerDestroy(&appctx.viewer2)); 231 PetscCall(VecDestroy(&appctx.solution)); 232 233 /* 234 Always call PetscFinalize() before exiting a program. This routine 235 - finalizes the PETSc libraries as well as MPI 236 - provides summary and diagnostic information if certain runtime 237 options are chosen (e.g., -log_view). 238 */ 239 PetscCall(PetscFinalize()); 240 return 0; 241 } 242 /* --------------------------------------------------------------------- */ 243 /* 244 InitialConditions - Computes the solution at the initial time. 245 246 Input Parameter: 247 u - uninitialized solution vector (global) 248 appctx - user-defined application context 249 250 Output Parameter: 251 u - vector with solution at initial time (global) 252 */ 253 PetscErrorCode InitialConditions(Vec u, AppCtx *appctx) 254 { 255 PetscScalar *u_localptr; 256 PetscInt i; 257 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 286 return 0; 287 } 288 /* --------------------------------------------------------------------- */ 289 /* 290 ExactSolution - Computes the exact solution at a given time. 291 292 Input Parameters: 293 t - current time 294 solution - vector in which exact solution will be computed 295 appctx - user-defined application context 296 297 Output Parameter: 298 solution - vector with the newly computed exact solution 299 */ 300 PetscErrorCode ExactSolution(PetscReal t, Vec solution, AppCtx *appctx) 301 { 302 PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2; 303 PetscInt i; 304 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 return 0; 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 */ 345 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) 346 { 347 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */ 348 PetscReal norm_2, norm_max, dt, dttol; 349 PetscBool flg; 350 351 /* 352 View a graph of the current iterate 353 */ 354 PetscCall(VecView(u, appctx->viewer2)); 355 356 /* 357 Compute the exact solution 358 */ 359 PetscCall(ExactSolution(crtime, appctx->solution, appctx)); 360 361 /* 362 Print debugging information if desired 363 */ 364 if (appctx->debug) { 365 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n")); 366 PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF)); 367 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n")); 368 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF)); 369 } 370 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 379 PetscCall(TSGetTimeStep(ts, &dt)); 380 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)); 381 appctx->norm_2 += norm_2; 382 appctx->norm_max += norm_max; 383 384 dttol = .0001; 385 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, &flg)); 386 if (dt < dttol) { 387 dt *= .999; 388 PetscCall(TSSetTimeStep(ts, dt)); 389 } 390 391 /* 392 View a graph of the error 393 */ 394 PetscCall(VecView(appctx->solution, appctx->viewer1)); 395 396 /* 397 Print debugging information if desired 398 */ 399 if (appctx->debug) { 400 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n")); 401 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF)); 402 } 403 404 return 0; 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 preconditioning matrix 420 str - flag indicating matrix structure 421 422 Notes: 423 Recall that MatSetValues() uses 0-based row and column numbers 424 in Fortran as well as in C. 425 */ 426 PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx) 427 { 428 Mat A = AA; /* Jacobian matrix */ 429 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */ 430 PetscInt mstart = 0; 431 PetscInt mend = appctx->m; 432 PetscInt i, idx[3]; 433 PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo; 434 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 483 return 0; 484 } 485 /* --------------------------------------------------------------------- */ 486 /* 487 Input Parameters: 488 ts - the TS context 489 t - current time 490 f - function 491 ctx - optional user-defined context, as set by TSetBCFunction() 492 */ 493 PetscErrorCode MyBCRoutine(TS ts, PetscReal t, Vec f, void *ctx) 494 { 495 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */ 496 PetscInt m = appctx->m; 497 PetscScalar *fa; 498 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 505 return 0; 506 } 507 508 /*TEST 509 510 test: 511 args: -nox -ts_max_steps 4 512 513 TEST*/ 514