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