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