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