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