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 { 78 AppCtx appctx; /* user-defined application context */ 79 TS ts; /* timestepping context */ 80 Mat A; /* matrix data structure */ 81 Vec u; /* approximate solution vector */ 82 PetscReal time_total_max = 100.0; /* default max total time */ 83 PetscInt time_steps_max = 100; /* default max timesteps */ 84 PetscDraw draw; /* drawing context */ 85 PetscInt steps, m; 86 PetscMPIInt size; 87 PetscReal dt; 88 PetscReal ftime; 89 PetscBool flg; 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Initialize program and set problem parameters 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 93 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 /* 258 Get a pointer to vector data. 259 - For default PETSc vectors, VecGetArray() returns a pointer to 260 the data array. Otherwise, the routine is implementation dependent. 261 - You MUST call VecRestoreArray() when you no longer need access to 262 the array. 263 - Note that the Fortran interface to VecGetArray() differs from the 264 C version. See the users manual for details. 265 */ 266 PetscCall(VecGetArray(u,&u_localptr)); 267 268 /* 269 We initialize the solution array by simply writing the solution 270 directly into the array locations. Alternatively, we could use 271 VecSetValues() or VecSetValuesLocal(). 272 */ 273 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); 274 275 /* 276 Restore vector 277 */ 278 PetscCall(VecRestoreArray(u,&u_localptr)); 279 280 /* 281 Print debugging information if desired 282 */ 283 if (appctx->debug) PetscCall(VecView(u,PETSC_VIEWER_STDOUT_SELF)); 284 285 return 0; 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 /* 305 Get a pointer to vector data. 306 */ 307 PetscCall(VecGetArray(solution,&s_localptr)); 308 309 /* 310 Simply write the solution directly into the array locations. 311 Alternatively, we culd use VecSetValues() or VecSetValuesLocal(). 312 */ 313 ex1 = PetscExpReal(-36.*PETSC_PI*PETSC_PI*t); ex2 = PetscExpReal(-4.*PETSC_PI*PETSC_PI*t); 314 sc1 = PETSC_PI*6.*h; 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 { 344 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 345 PetscReal norm_2, norm_max, dt, dttol; 346 PetscBool flg; 347 348 /* 349 View a graph of the current iterate 350 */ 351 PetscCall(VecView(u,appctx->viewer2)); 352 353 /* 354 Compute the exact solution 355 */ 356 PetscCall(ExactSolution(crtime,appctx->solution,appctx)); 357 358 /* 359 Print debugging information if desired 360 */ 361 if (appctx->debug) { 362 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Computed solution vector\n")); 363 PetscCall(VecView(u,PETSC_VIEWER_STDOUT_SELF)); 364 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n")); 365 PetscCall(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 366 } 367 368 /* 369 Compute the 2-norm and max-norm of the error 370 */ 371 PetscCall(VecAXPY(appctx->solution,-1.0,u)); 372 PetscCall(VecNorm(appctx->solution,NORM_2,&norm_2)); 373 norm_2 = PetscSqrtReal(appctx->h)*norm_2; 374 PetscCall(VecNorm(appctx->solution,NORM_MAX,&norm_max)); 375 376 PetscCall(TSGetTimeStep(ts,&dt)); 377 if (norm_2 > 1.e-2) { 378 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)); 379 } 380 appctx->norm_2 += norm_2; 381 appctx->norm_max += norm_max; 382 383 dttol = .0001; 384 PetscCall(PetscOptionsGetReal(NULL,NULL,"-dttol",&dttol,&flg)); 385 if (dt < dttol) { 386 dt *= .999; 387 PetscCall(TSSetTimeStep(ts,dt)); 388 } 389 390 /* 391 View a graph of the error 392 */ 393 PetscCall(VecView(appctx->solution,appctx->viewer1)); 394 395 /* 396 Print debugging information if desired 397 */ 398 if (appctx->debug) { 399 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error vector\n")); 400 PetscCall(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 401 } 402 403 return 0; 404 } 405 /* --------------------------------------------------------------------- */ 406 /* 407 RHSMatrixHeat - User-provided routine to compute the right-hand-side 408 matrix for the heat equation. 409 410 Input Parameters: 411 ts - the TS context 412 t - current time 413 global_in - global input vector 414 dummy - optional user-defined context, as set by TSetRHSJacobian() 415 416 Output Parameters: 417 AA - Jacobian matrix 418 BB - optionally different preconditioning matrix 419 str - flag indicating matrix structure 420 421 Notes: 422 Recall that MatSetValues() uses 0-based row and column numbers 423 in Fortran as well as in C. 424 */ 425 PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx) 426 { 427 Mat A = AA; /* Jacobian matrix */ 428 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 429 PetscInt mstart = 0; 430 PetscInt mend = appctx->m; 431 PetscInt i, idx[3]; 432 PetscScalar v[3], stwo = -2./(appctx->h*appctx->h), sone = -.5*stwo; 433 434 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 435 Compute entries for the locally owned part of the matrix 436 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 437 /* 438 Set matrix rows corresponding to boundary data 439 */ 440 441 mstart = 0; 442 v[0] = 1.0; 443 PetscCall(MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES)); 444 mstart++; 445 446 mend--; 447 v[0] = 1.0; 448 PetscCall(MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES)); 449 450 /* 451 Set matrix rows corresponding to interior data. We construct the 452 matrix one row at a time. 453 */ 454 v[0] = sone; v[1] = stwo; v[2] = sone; 455 for (i=mstart; i<mend; i++) { 456 idx[0] = i-1; idx[1] = i; 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 { 490 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 491 PetscInt m = appctx->m; 492 PetscScalar *fa; 493 494 PetscCall(VecGetArray(f,&fa)); 495 fa[0] = 0.0; 496 fa[m-1] = 1.0; 497 PetscCall(VecRestoreArray(f,&fa)); 498 PetscCall(PetscPrintf(PETSC_COMM_SELF,"t=%g\n",(double)t)); 499 500 return 0; 501 } 502 503 /*TEST 504 505 test: 506 args: -nox -ts_max_steps 4 507 508 TEST*/ 509