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