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 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 109 110 m = 60; 111 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 112 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug)); 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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Solving a linear TS problem on 1 processor\n")); 120 121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 122 Create vector data structures 123 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 124 125 /* 126 Create vector data structures for approximate and exact solutions 127 */ 128 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,m,&u)); 129 CHKERRQ(VecDuplicate(u,&appctx.solution)); 130 131 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 132 Set up displays to show graphs of the solution and error 133 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 134 135 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,380,400,160,&appctx.viewer1)); 136 CHKERRQ(PetscViewerDrawGetDraw(appctx.viewer1,0,&draw)); 137 CHKERRQ(PetscDrawSetDoubleBuffer(draw)); 138 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,0,400,160,&appctx.viewer2)); 139 CHKERRQ(PetscViewerDrawGetDraw(appctx.viewer2,0,&draw)); 140 CHKERRQ(PetscDrawSetDoubleBuffer(draw)); 141 142 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 Create timestepping solver context 144 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 145 146 CHKERRQ(TSCreate(PETSC_COMM_SELF,&ts)); 147 CHKERRQ(TSSetProblemType(ts,TS_LINEAR)); 148 149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 Set optional user-defined monitoring routine 151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152 153 CHKERRQ(TSMonitorSet(ts,Monitor,&appctx,NULL)); 154 155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156 157 Create matrix data structure; set matrix evaluation routine. 158 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 159 160 CHKERRQ(MatCreate(PETSC_COMM_SELF,&A)); 161 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m)); 162 CHKERRQ(MatSetFromOptions(A)); 163 CHKERRQ(MatSetUp(A)); 164 165 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-time_dependent_rhs",&flg)); 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 CHKERRQ(TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx)); 173 CHKERRQ(TSSetRHSJacobian(ts,A,A,RHSMatrixHeat,&appctx)); 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 CHKERRQ(RHSMatrixHeat(ts,0.0,u,A,A,&appctx)); 182 CHKERRQ(TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx)); 183 CHKERRQ(TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx)); 184 } 185 186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 187 Set solution vector and initial timestep 188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 189 190 dt = appctx.h*appctx.h/2.0; 191 CHKERRQ(TSSetTimeStep(ts,dt)); 192 CHKERRQ(TSSetSolution(ts,u)); 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 CHKERRQ(TSSetMaxSteps(ts,time_steps_max)); 205 CHKERRQ(TSSetMaxTime(ts,time_total_max)); 206 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 207 CHKERRQ(TSSetFromOptions(ts)); 208 209 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 210 Solve the problem 211 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 212 213 /* 214 Evaluate initial conditions 215 */ 216 CHKERRQ(InitialConditions(u,&appctx)); 217 218 /* 219 Run the timestepping solver 220 */ 221 CHKERRQ(TSSolve(ts,u)); 222 CHKERRQ(TSGetSolveTime(ts,&ftime)); 223 CHKERRQ(TSGetStepNumber(ts,&steps)); 224 225 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 226 View timestepping solver info 227 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 228 229 CHKERRQ(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))); 230 CHKERRQ(TSView(ts,PETSC_VIEWER_STDOUT_SELF)); 231 232 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233 Free work space. All PETSc objects should be destroyed when they 234 are no longer needed. 235 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 236 237 CHKERRQ(TSDestroy(&ts)); 238 CHKERRQ(MatDestroy(&A)); 239 CHKERRQ(VecDestroy(&u)); 240 CHKERRQ(PetscViewerDestroy(&appctx.viewer1)); 241 CHKERRQ(PetscViewerDestroy(&appctx.viewer2)); 242 CHKERRQ(VecDestroy(&appctx.solution)); 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 269 /* 270 Get a pointer to vector data. 271 - For default PETSc vectors, VecGetArray() returns a pointer to 272 the data array. Otherwise, the routine is implementation dependent. 273 - You MUST call VecRestoreArray() when you no longer need access to 274 the array. 275 - Note that the Fortran interface to VecGetArray() differs from the 276 C version. See the users manual for details. 277 */ 278 CHKERRQ(VecGetArray(u,&u_localptr)); 279 280 /* 281 We initialize the solution array by simply writing the solution 282 directly into the array locations. Alternatively, we could use 283 VecSetValues() or VecSetValuesLocal(). 284 */ 285 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); 286 287 /* 288 Restore vector 289 */ 290 CHKERRQ(VecRestoreArray(u,&u_localptr)); 291 292 /* 293 Print debugging information if desired 294 */ 295 if (appctx->debug) { 296 CHKERRQ(VecView(u,PETSC_VIEWER_STDOUT_SELF)); 297 } 298 299 return 0; 300 } 301 /* --------------------------------------------------------------------- */ 302 /* 303 ExactSolution - Computes the exact solution at a given time. 304 305 Input Parameters: 306 t - current time 307 solution - vector in which exact solution will be computed 308 appctx - user-defined application context 309 310 Output Parameter: 311 solution - vector with the newly computed exact solution 312 */ 313 PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx) 314 { 315 PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2; 316 PetscInt i; 317 318 /* 319 Get a pointer to vector data. 320 */ 321 CHKERRQ(VecGetArray(solution,&s_localptr)); 322 323 /* 324 Simply write the solution directly into the array locations. 325 Alternatively, we culd use VecSetValues() or VecSetValuesLocal(). 326 */ 327 ex1 = PetscExpReal(-36.*PETSC_PI*PETSC_PI*t); ex2 = PetscExpReal(-4.*PETSC_PI*PETSC_PI*t); 328 sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h; 329 for (i=0; i<appctx->m; i++) s_localptr[i] = PetscSinReal(PetscRealPart(sc1)*(PetscReal)i)*ex1 + 3.*PetscSinReal(PetscRealPart(sc2)*(PetscReal)i)*ex2; 330 331 /* 332 Restore vector 333 */ 334 CHKERRQ(VecRestoreArray(solution,&s_localptr)); 335 return 0; 336 } 337 /* --------------------------------------------------------------------- */ 338 /* 339 Monitor - User-provided routine to monitor the solution computed at 340 each timestep. This example plots the solution and computes the 341 error in two different norms. 342 343 This example also demonstrates changing the timestep via TSSetTimeStep(). 344 345 Input Parameters: 346 ts - the timestep context 347 step - the count of the current step (with 0 meaning the 348 initial condition) 349 crtime - the current time 350 u - the solution at this timestep 351 ctx - the user-provided context for this monitoring routine. 352 In this case we use the application context which contains 353 information about the problem size, workspace and the exact 354 solution. 355 */ 356 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal crtime,Vec u,void *ctx) 357 { 358 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 359 PetscReal norm_2, norm_max, dt, dttol; 360 PetscBool flg; 361 362 /* 363 View a graph of the current iterate 364 */ 365 CHKERRQ(VecView(u,appctx->viewer2)); 366 367 /* 368 Compute the exact solution 369 */ 370 CHKERRQ(ExactSolution(crtime,appctx->solution,appctx)); 371 372 /* 373 Print debugging information if desired 374 */ 375 if (appctx->debug) { 376 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Computed solution vector\n")); 377 CHKERRQ(VecView(u,PETSC_VIEWER_STDOUT_SELF)); 378 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n")); 379 CHKERRQ(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 380 } 381 382 /* 383 Compute the 2-norm and max-norm of the error 384 */ 385 CHKERRQ(VecAXPY(appctx->solution,-1.0,u)); 386 CHKERRQ(VecNorm(appctx->solution,NORM_2,&norm_2)); 387 norm_2 = PetscSqrtReal(appctx->h)*norm_2; 388 CHKERRQ(VecNorm(appctx->solution,NORM_MAX,&norm_max)); 389 390 CHKERRQ(TSGetTimeStep(ts,&dt)); 391 if (norm_2 > 1.e-2) { 392 CHKERRQ(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)); 393 } 394 appctx->norm_2 += norm_2; 395 appctx->norm_max += norm_max; 396 397 dttol = .0001; 398 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-dttol",&dttol,&flg)); 399 if (dt < dttol) { 400 dt *= .999; 401 CHKERRQ(TSSetTimeStep(ts,dt)); 402 } 403 404 /* 405 View a graph of the error 406 */ 407 CHKERRQ(VecView(appctx->solution,appctx->viewer1)); 408 409 /* 410 Print debugging information if desired 411 */ 412 if (appctx->debug) { 413 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error vector\n")); 414 CHKERRQ(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 415 } 416 417 return 0; 418 } 419 /* --------------------------------------------------------------------- */ 420 /* 421 RHSMatrixHeat - User-provided routine to compute the right-hand-side 422 matrix for the heat equation. 423 424 Input Parameters: 425 ts - the TS context 426 t - current time 427 global_in - global input vector 428 dummy - optional user-defined context, as set by TSetRHSJacobian() 429 430 Output Parameters: 431 AA - Jacobian matrix 432 BB - optionally different preconditioning matrix 433 str - flag indicating matrix structure 434 435 Notes: 436 Recall that MatSetValues() uses 0-based row and column numbers 437 in Fortran as well as in C. 438 */ 439 PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx) 440 { 441 Mat A = AA; /* Jacobian matrix */ 442 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 443 PetscInt mstart = 0; 444 PetscInt mend = appctx->m; 445 PetscInt i, idx[3]; 446 PetscScalar v[3], stwo = -2./(appctx->h*appctx->h), sone = -.5*stwo; 447 448 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 449 Compute entries for the locally owned part of the matrix 450 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 451 /* 452 Set matrix rows corresponding to boundary data 453 */ 454 455 mstart = 0; 456 v[0] = 1.0; 457 CHKERRQ(MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES)); 458 mstart++; 459 460 mend--; 461 v[0] = 1.0; 462 CHKERRQ(MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES)); 463 464 /* 465 Set matrix rows corresponding to interior data. We construct the 466 matrix one row at a time. 467 */ 468 v[0] = sone; v[1] = stwo; v[2] = sone; 469 for (i=mstart; i<mend; i++) { 470 idx[0] = i-1; idx[1] = i; idx[2] = i+1; 471 CHKERRQ(MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES)); 472 } 473 474 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 475 Complete the matrix assembly process and set some options 476 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 477 /* 478 Assemble matrix, using the 2-step process: 479 MatAssemblyBegin(), MatAssemblyEnd() 480 Computations can be done while messages are in transition 481 by placing code between these two statements. 482 */ 483 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 484 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 485 486 /* 487 Set and option to indicate that we will never add a new nonzero location 488 to the matrix. If we do, it will generate an error. 489 */ 490 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 491 492 return 0; 493 } 494 /* --------------------------------------------------------------------- */ 495 /* 496 Input Parameters: 497 ts - the TS context 498 t - current time 499 f - function 500 ctx - optional user-defined context, as set by TSetBCFunction() 501 */ 502 PetscErrorCode MyBCRoutine(TS ts,PetscReal t,Vec f,void *ctx) 503 { 504 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 505 PetscInt m = appctx->m; 506 PetscScalar *fa; 507 508 CHKERRQ(VecGetArray(f,&fa)); 509 fa[0] = 0.0; 510 fa[m-1] = 1.0; 511 CHKERRQ(VecRestoreArray(f,&fa)); 512 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"t=%g\n",(double)t)); 513 514 return 0; 515 } 516 517 /*TEST 518 519 test: 520 args: -nox -ts_max_steps 4 521 522 TEST*/ 523