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