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