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 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, (char *)0, 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 */ 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 286 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 306 /* 307 Get a pointer to vector data. 308 */ 309 PetscCall(VecGetArray(solution, &s_localptr)); 310 311 /* 312 Simply write the solution directly into the array locations. 313 Alternatively, we culd use VecSetValues() or VecSetValuesLocal(). 314 */ 315 ex1 = PetscExpReal(-36. * PETSC_PI * PETSC_PI * t); 316 ex2 = PetscExpReal(-4. * PETSC_PI * PETSC_PI * t); 317 sc1 = PETSC_PI * 6. * h; 318 sc2 = PETSC_PI * 2. * h; 319 for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscSinReal(PetscRealPart(sc1) * (PetscReal)i) * ex1 + 3. * PetscSinReal(PetscRealPart(sc2) * (PetscReal)i) * ex2; 320 321 /* 322 Restore vector 323 */ 324 PetscCall(VecRestoreArray(solution, &s_localptr)); 325 PetscFunctionReturn(PETSC_SUCCESS); 326 } 327 /* --------------------------------------------------------------------- */ 328 /* 329 Monitor - User-provided routine to monitor the solution computed at 330 each timestep. This example plots the solution and computes the 331 error in two different norms. 332 333 This example also demonstrates changing the timestep via TSSetTimeStep(). 334 335 Input Parameters: 336 ts - the timestep context 337 step - the count of the current step (with 0 meaning the 338 initial condition) 339 crtime - the current time 340 u - the solution at this timestep 341 ctx - the user-provided context for this monitoring routine. 342 In this case we use the application context which contains 343 information about the problem size, workspace and the exact 344 solution. 345 */ 346 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) 347 { 348 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */ 349 PetscReal norm_2, norm_max, dt, dttol; 350 PetscBool flg; 351 352 PetscFunctionBeginUser; 353 /* 354 View a graph of the current iterate 355 */ 356 PetscCall(VecView(u, appctx->viewer2)); 357 358 /* 359 Compute the exact solution 360 */ 361 PetscCall(ExactSolution(crtime, appctx->solution, appctx)); 362 363 /* 364 Print debugging information if desired 365 */ 366 if (appctx->debug) { 367 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n")); 368 PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF)); 369 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n")); 370 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF)); 371 } 372 373 /* 374 Compute the 2-norm and max-norm of the error 375 */ 376 PetscCall(VecAXPY(appctx->solution, -1.0, u)); 377 PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2)); 378 norm_2 = PetscSqrtReal(appctx->h) * norm_2; 379 PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max)); 380 381 PetscCall(TSGetTimeStep(ts, &dt)); 382 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)); 383 appctx->norm_2 += norm_2; 384 appctx->norm_max += norm_max; 385 386 dttol = .0001; 387 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, &flg)); 388 if (dt < dttol) { 389 dt *= .999; 390 PetscCall(TSSetTimeStep(ts, dt)); 391 } 392 393 /* 394 View a graph of the error 395 */ 396 PetscCall(VecView(appctx->solution, appctx->viewer1)); 397 398 /* 399 Print debugging information if desired 400 */ 401 if (appctx->debug) { 402 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n")); 403 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF)); 404 } 405 406 PetscFunctionReturn(PETSC_SUCCESS); 407 } 408 /* --------------------------------------------------------------------- */ 409 /* 410 RHSMatrixHeat - User-provided routine to compute the right-hand-side 411 matrix for the heat equation. 412 413 Input Parameters: 414 ts - the TS context 415 t - current time 416 global_in - global input vector 417 dummy - optional user-defined context, as set by TSetRHSJacobian() 418 419 Output Parameters: 420 AA - Jacobian matrix 421 BB - optionally different preconditioning matrix 422 str - flag indicating matrix structure 423 424 Notes: 425 Recall that MatSetValues() uses 0-based row and column numbers 426 in Fortran as well as in C. 427 */ 428 PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx) 429 { 430 Mat A = AA; /* Jacobian matrix */ 431 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */ 432 PetscInt mstart = 0; 433 PetscInt mend = appctx->m; 434 PetscInt i, idx[3]; 435 PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo; 436 437 PetscFunctionBeginUser; 438 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 439 Compute entries for the locally owned part of the matrix 440 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 441 /* 442 Set matrix rows corresponding to boundary data 443 */ 444 445 mstart = 0; 446 v[0] = 1.0; 447 PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES)); 448 mstart++; 449 450 mend--; 451 v[0] = 1.0; 452 PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES)); 453 454 /* 455 Set matrix rows corresponding to interior data. We construct the 456 matrix one row at a time. 457 */ 458 v[0] = sone; 459 v[1] = stwo; 460 v[2] = sone; 461 for (i = mstart; i < mend; i++) { 462 idx[0] = i - 1; 463 idx[1] = i; 464 idx[2] = i + 1; 465 PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES)); 466 } 467 468 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 469 Complete the matrix assembly process and set some options 470 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 471 /* 472 Assemble matrix, using the 2-step process: 473 MatAssemblyBegin(), MatAssemblyEnd() 474 Computations can be done while messages are in transition 475 by placing code between these two statements. 476 */ 477 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 478 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 479 480 /* 481 Set and option to indicate that we will never add a new nonzero location 482 to the matrix. If we do, it will generate an error. 483 */ 484 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 485 486 PetscFunctionReturn(PETSC_SUCCESS); 487 } 488 /* --------------------------------------------------------------------- */ 489 /* 490 Input Parameters: 491 ts - the TS context 492 t - current time 493 f - function 494 ctx - optional user-defined context, as set by TSetBCFunction() 495 */ 496 PetscErrorCode MyBCRoutine(TS ts, PetscReal t, Vec f, void *ctx) 497 { 498 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */ 499 PetscInt m = appctx->m; 500 PetscScalar *fa; 501 502 PetscFunctionBeginUser; 503 PetscCall(VecGetArray(f, &fa)); 504 fa[0] = 0.0; 505 fa[m - 1] = 1.0; 506 PetscCall(VecRestoreArray(f, &fa)); 507 PetscCall(PetscPrintf(PETSC_COMM_SELF, "t=%g\n", (double)t)); 508 509 PetscFunctionReturn(PETSC_SUCCESS); 510 } 511 512 /*TEST 513 514 test: 515 args: -nox -ts_max_steps 4 516 517 TEST*/ 518