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