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 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 */ 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 */ 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 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 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 PetscFunctionBeginUser; 436 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 437 Compute entries for the locally owned part of the matrix 438 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 439 /* 440 Set matrix rows corresponding to boundary data 441 */ 442 443 mstart = 0; 444 v[0] = 1.0; 445 PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES)); 446 mstart++; 447 448 mend--; 449 v[0] = 1.0; 450 PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES)); 451 452 /* 453 Set matrix rows corresponding to interior data. We construct the 454 matrix one row at a time. 455 */ 456 v[0] = sone; 457 v[1] = stwo; 458 v[2] = sone; 459 for (i = mstart; i < mend; i++) { 460 idx[0] = i - 1; 461 idx[1] = i; 462 idx[2] = i + 1; 463 PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES)); 464 } 465 466 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 467 Complete the matrix assembly process and set some options 468 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 469 /* 470 Assemble matrix, using the 2-step process: 471 MatAssemblyBegin(), MatAssemblyEnd() 472 Computations can be done while messages are in transition 473 by placing code between these two statements. 474 */ 475 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 476 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 477 478 /* 479 Set and option to indicate that we will never add a new nonzero location 480 to the matrix. If we do, it will generate an error. 481 */ 482 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 483 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 500 PetscCall(VecGetArray(f, &fa)); 501 fa[0] = 0.0; 502 fa[m - 1] = 1.0; 503 PetscCall(VecRestoreArray(f, &fa)); 504 PetscCall(PetscPrintf(PETSC_COMM_SELF, "t=%g\n", (double)t)); 505 PetscFunctionReturn(PETSC_SUCCESS); 506 } 507 508 /*TEST 509 510 test: 511 args: -nox -ts_max_steps 4 512 513 TEST*/ 514