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