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 -use_ifunc : Use IFunction/IJacobian interface\n\ 7 -debug : Activate debugging printouts\n\ 8 -nox : Deactivate x-window graphics\n\n"; 9 10 /* 11 Concepts: TS^time-dependent linear problems 12 Concepts: TS^heat equation 13 Concepts: TS^diffusion equation 14 Processors: 1 15 */ 16 17 /* ------------------------------------------------------------------------ 18 19 This program solves the one-dimensional heat equation (also called the 20 diffusion equation), 21 u_t = u_xx, 22 on the domain 0 <= x <= 1, with the boundary conditions 23 u(t,0) = 0, u(t,1) = 0, 24 and the initial condition 25 u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x). 26 This is a linear, second-order, parabolic equation. 27 28 We discretize the right-hand side using finite differences with 29 uniform grid spacing h: 30 u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2) 31 We then demonstrate time evolution using the various TS methods by 32 running the program via 33 ex3 -ts_type <timestepping solver> 34 35 We compare the approximate solution with the exact solution, given by 36 u_exact(x,t) = exp(-36*pi*pi*t) * sin(6*pi*x) + 37 3*exp(-4*pi*pi*t) * sin(2*pi*x) 38 39 Notes: 40 This code demonstrates the TS solver interface to two variants of 41 linear problems, u_t = f(u,t), namely 42 - time-dependent f: f(u,t) is a function of t 43 - time-independent f: f(u,t) is simply f(u) 44 45 The parallel version of this code is ts/tutorials/ex4.c 46 47 ------------------------------------------------------------------------- */ 48 49 /* 50 Include "petscts.h" so that we can use TS solvers. Note that this file 51 automatically includes: 52 petscsys.h - base PETSc routines petscvec.h - vectors 53 petscmat.h - matrices 54 petscis.h - index sets petscksp.h - Krylov subspace methods 55 petscviewer.h - viewers petscpc.h - preconditioners 56 petscksp.h - linear solvers petscsnes.h - nonlinear solvers 57 */ 58 59 #include <petscts.h> 60 #include <petscdraw.h> 61 62 /* 63 User-defined application context - contains data needed by the 64 application-provided call-back routines. 65 */ 66 typedef struct { 67 Vec solution; /* global exact solution vector */ 68 PetscInt m; /* total number of grid points */ 69 PetscReal h; /* mesh width h = 1/(m-1) */ 70 PetscBool debug; /* flag (1 indicates activation of debugging printouts) */ 71 PetscViewer viewer1,viewer2; /* viewers for the solution and error */ 72 PetscReal norm_2,norm_max; /* error norms */ 73 Mat A; /* RHS mat, used with IFunction interface */ 74 PetscReal oshift; /* old shift applied, prevent to recompute the IJacobian */ 75 } AppCtx; 76 77 /* 78 User-defined routines 79 */ 80 extern PetscErrorCode InitialConditions(Vec,AppCtx*); 81 extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Vec,Mat,Mat,void*); 82 extern PetscErrorCode IFunctionHeat(TS,PetscReal,Vec,Vec,Vec,void*); 83 extern PetscErrorCode IJacobianHeat(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 84 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); 85 extern PetscErrorCode ExactSolution(PetscReal,Vec,AppCtx*); 86 87 int main(int argc,char **argv) 88 { 89 AppCtx appctx; /* user-defined application context */ 90 TS ts; /* timestepping context */ 91 Mat A; /* matrix data structure */ 92 Vec u; /* approximate solution vector */ 93 PetscReal time_total_max = 100.0; /* default max total time */ 94 PetscInt time_steps_max = 100; /* default max timesteps */ 95 PetscDraw draw; /* drawing context */ 96 PetscErrorCode ierr; 97 PetscInt steps,m; 98 PetscMPIInt size; 99 PetscReal dt; 100 PetscBool flg,flg_string; 101 102 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 103 Initialize program and set problem parameters 104 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 105 106 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 107 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 108 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 109 110 m = 60; 111 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 112 ierr = PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug);CHKERRQ(ierr); 113 flg_string = PETSC_FALSE; 114 ierr = PetscOptionsGetBool(NULL,NULL,"-test_string_viewer",&flg_string,NULL);CHKERRQ(ierr); 115 116 appctx.m = m; 117 appctx.h = 1.0/(m-1.0); 118 appctx.norm_2 = 0.0; 119 appctx.norm_max = 0.0; 120 121 ierr = PetscPrintf(PETSC_COMM_SELF,"Solving a linear TS problem on 1 processor\n");CHKERRQ(ierr); 122 123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124 Create vector data structures 125 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 126 127 /* 128 Create vector data structures for approximate and exact solutions 129 */ 130 ierr = VecCreateSeq(PETSC_COMM_SELF,m,&u);CHKERRQ(ierr); 131 ierr = VecDuplicate(u,&appctx.solution);CHKERRQ(ierr); 132 133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 134 Set up displays to show graphs of the solution and error 135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 136 137 ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,380,400,160,&appctx.viewer1);CHKERRQ(ierr); 138 ierr = PetscViewerDrawGetDraw(appctx.viewer1,0,&draw);CHKERRQ(ierr); 139 ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr); 140 ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,0,400,160,&appctx.viewer2);CHKERRQ(ierr); 141 ierr = PetscViewerDrawGetDraw(appctx.viewer2,0,&draw);CHKERRQ(ierr); 142 ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr); 143 144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 Create timestepping solver context 146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 147 148 ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr); 149 ierr = TSSetProblemType(ts,TS_LINEAR);CHKERRQ(ierr); 150 151 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152 Set optional user-defined monitoring routine 153 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 154 155 if (!flg_string) { 156 ierr = TSMonitorSet(ts,Monitor,&appctx,NULL);CHKERRQ(ierr); 157 } 158 159 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160 161 Create matrix data structure; set matrix evaluation routine. 162 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 163 164 ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr); 165 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m);CHKERRQ(ierr); 166 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 167 ierr = MatSetUp(A);CHKERRQ(ierr); 168 169 flg = PETSC_FALSE; 170 ierr = PetscOptionsGetBool(NULL,NULL,"-use_ifunc",&flg,NULL);CHKERRQ(ierr); 171 if (!flg) { 172 appctx.A = NULL; 173 ierr = PetscOptionsGetBool(NULL,NULL,"-time_dependent_rhs",&flg,NULL);CHKERRQ(ierr); 174 if (flg) { 175 /* 176 For linear problems with a time-dependent f(u,t) in the equation 177 u_t = f(u,t), the user provides the discretized right-hand-side 178 as a time-dependent matrix. 179 */ 180 ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);CHKERRQ(ierr); 181 ierr = TSSetRHSJacobian(ts,A,A,RHSMatrixHeat,&appctx);CHKERRQ(ierr); 182 } else { 183 /* 184 For linear problems with a time-independent f(u) in the equation 185 u_t = f(u), the user provides the discretized right-hand-side 186 as a matrix only once, and then sets the special Jacobian evaluation 187 routine TSComputeRHSJacobianConstant() which will NOT recompute the Jacobian. 188 */ 189 ierr = RHSMatrixHeat(ts,0.0,u,A,A,&appctx);CHKERRQ(ierr); 190 ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);CHKERRQ(ierr); 191 ierr = TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx);CHKERRQ(ierr); 192 } 193 } else { 194 Mat J; 195 196 ierr = RHSMatrixHeat(ts,0.0,u,A,A,&appctx);CHKERRQ(ierr); 197 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&J);CHKERRQ(ierr); 198 ierr = TSSetIFunction(ts,NULL,IFunctionHeat,&appctx);CHKERRQ(ierr); 199 ierr = TSSetIJacobian(ts,J,J,IJacobianHeat,&appctx);CHKERRQ(ierr); 200 ierr = MatDestroy(&J);CHKERRQ(ierr); 201 202 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 203 appctx.A = A; 204 appctx.oshift = PETSC_MIN_REAL; 205 } 206 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207 Set solution vector and initial timestep 208 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 209 210 dt = appctx.h*appctx.h/2.0; 211 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 212 213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 214 Customize timestepping solver: 215 - Set the solution method to be the Backward Euler method. 216 - Set timestepping duration info 217 Then set runtime options, which can override these defaults. 218 For example, 219 -ts_max_steps <maxsteps> -ts_max_time <maxtime> 220 to override the defaults set by TSSetMaxSteps()/TSSetMaxTime(). 221 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 222 223 ierr = TSSetMaxSteps(ts,time_steps_max);CHKERRQ(ierr); 224 ierr = TSSetMaxTime(ts,time_total_max);CHKERRQ(ierr); 225 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 226 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 227 228 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 229 Solve the problem 230 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 231 232 /* 233 Evaluate initial conditions 234 */ 235 ierr = InitialConditions(u,&appctx);CHKERRQ(ierr); 236 237 /* 238 Run the timestepping solver 239 */ 240 ierr = TSSolve(ts,u);CHKERRQ(ierr); 241 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 242 243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244 View timestepping solver info 245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 246 247 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); 248 if (!flg_string) { 249 ierr = TSView(ts,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 250 } else { 251 PetscViewer stringviewer; 252 char string[512]; 253 const char *outstring; 254 255 ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,string,sizeof(string),&stringviewer);CHKERRQ(ierr); 256 ierr = TSView(ts,stringviewer);CHKERRQ(ierr); 257 ierr = PetscViewerStringGetStringRead(stringviewer,&outstring,NULL);CHKERRQ(ierr); 258 if ((char*)outstring != (char*)string) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"String returned from viewer does not equal original string"); 259 ierr = PetscPrintf(PETSC_COMM_WORLD,"Output from string viewer:%s\n",outstring);CHKERRQ(ierr); 260 ierr = PetscViewerDestroy(&stringviewer);CHKERRQ(ierr); 261 } 262 263 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 264 Free work space. All PETSc objects should be destroyed when they 265 are no longer needed. 266 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 267 268 ierr = TSDestroy(&ts);CHKERRQ(ierr); 269 ierr = MatDestroy(&A);CHKERRQ(ierr); 270 ierr = VecDestroy(&u);CHKERRQ(ierr); 271 ierr = PetscViewerDestroy(&appctx.viewer1);CHKERRQ(ierr); 272 ierr = PetscViewerDestroy(&appctx.viewer2);CHKERRQ(ierr); 273 ierr = VecDestroy(&appctx.solution);CHKERRQ(ierr); 274 ierr = MatDestroy(&appctx.A);CHKERRQ(ierr); 275 276 /* 277 Always call PetscFinalize() before exiting a program. This routine 278 - finalizes the PETSc libraries as well as MPI 279 - provides summary and diagnostic information if certain runtime 280 options are chosen (e.g., -log_view). 281 */ 282 ierr = PetscFinalize(); 283 return ierr; 284 } 285 /* --------------------------------------------------------------------- */ 286 /* 287 InitialConditions - Computes the solution at the initial time. 288 289 Input Parameter: 290 u - uninitialized solution vector (global) 291 appctx - user-defined application context 292 293 Output Parameter: 294 u - vector with solution at initial time (global) 295 */ 296 PetscErrorCode InitialConditions(Vec u,AppCtx *appctx) 297 { 298 PetscScalar *u_localptr,h = appctx->h; 299 PetscErrorCode ierr; 300 PetscInt i; 301 302 /* 303 Get a pointer to vector data. 304 - For default PETSc vectors, VecGetArray() returns a pointer to 305 the data array. Otherwise, the routine is implementation dependent. 306 - You MUST call VecRestoreArray() when you no longer need access to 307 the array. 308 - Note that the Fortran interface to VecGetArray() differs from the 309 C version. See the users manual for details. 310 */ 311 ierr = VecGetArray(u,&u_localptr);CHKERRQ(ierr); 312 313 /* 314 We initialize the solution array by simply writing the solution 315 directly into the array locations. Alternatively, we could use 316 VecSetValues() or VecSetValuesLocal(). 317 */ 318 for (i=0; i<appctx->m; i++) u_localptr[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h); 319 320 /* 321 Restore vector 322 */ 323 ierr = VecRestoreArray(u,&u_localptr);CHKERRQ(ierr); 324 325 /* 326 Print debugging information if desired 327 */ 328 if (appctx->debug) { 329 ierr = PetscPrintf(PETSC_COMM_WORLD,"Initial guess vector\n");CHKERRQ(ierr); 330 ierr = VecView(u,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 331 } 332 333 return 0; 334 } 335 /* --------------------------------------------------------------------- */ 336 /* 337 ExactSolution - Computes the exact solution at a given time. 338 339 Input Parameters: 340 t - current time 341 solution - vector in which exact solution will be computed 342 appctx - user-defined application context 343 344 Output Parameter: 345 solution - vector with the newly computed exact solution 346 */ 347 PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx) 348 { 349 PetscScalar *s_localptr,h = appctx->h,ex1,ex2,sc1,sc2,tc = t; 350 PetscErrorCode ierr; 351 PetscInt i; 352 353 /* 354 Get a pointer to vector data. 355 */ 356 ierr = VecGetArray(solution,&s_localptr);CHKERRQ(ierr); 357 358 /* 359 Simply write the solution directly into the array locations. 360 Alternatively, we culd use VecSetValues() or VecSetValuesLocal(). 361 */ 362 ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*tc); 363 ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*tc); 364 sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h; 365 for (i=0; i<appctx->m; i++) s_localptr[i] = PetscSinScalar(sc1*(PetscReal)i)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i)*ex2; 366 367 /* 368 Restore vector 369 */ 370 ierr = VecRestoreArray(solution,&s_localptr);CHKERRQ(ierr); 371 return 0; 372 } 373 /* --------------------------------------------------------------------- */ 374 /* 375 Monitor - User-provided routine to monitor the solution computed at 376 each timestep. This example plots the solution and computes the 377 error in two different norms. 378 379 This example also demonstrates changing the timestep via TSSetTimeStep(). 380 381 Input Parameters: 382 ts - the timestep context 383 step - the count of the current step (with 0 meaning the 384 initial condition) 385 time - the current time 386 u - the solution at this timestep 387 ctx - the user-provided context for this monitoring routine. 388 In this case we use the application context which contains 389 information about the problem size, workspace and the exact 390 solution. 391 */ 392 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx) 393 { 394 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 395 PetscErrorCode ierr; 396 PetscReal norm_2,norm_max,dt,dttol; 397 398 /* 399 View a graph of the current iterate 400 */ 401 ierr = VecView(u,appctx->viewer2);CHKERRQ(ierr); 402 403 /* 404 Compute the exact solution 405 */ 406 ierr = ExactSolution(time,appctx->solution,appctx);CHKERRQ(ierr); 407 408 /* 409 Print debugging information if desired 410 */ 411 if (appctx->debug) { 412 ierr = PetscPrintf(PETSC_COMM_SELF,"Computed solution vector\n");CHKERRQ(ierr); 413 ierr = VecView(u,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 414 ierr = PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n");CHKERRQ(ierr); 415 ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 416 } 417 418 /* 419 Compute the 2-norm and max-norm of the error 420 */ 421 ierr = VecAXPY(appctx->solution,-1.0,u);CHKERRQ(ierr); 422 ierr = VecNorm(appctx->solution,NORM_2,&norm_2);CHKERRQ(ierr); 423 norm_2 = PetscSqrtReal(appctx->h)*norm_2; 424 ierr = VecNorm(appctx->solution,NORM_MAX,&norm_max);CHKERRQ(ierr); 425 if (norm_2 < 1e-14) norm_2 = 0; 426 if (norm_max < 1e-14) norm_max = 0; 427 428 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 429 ierr = PetscPrintf(PETSC_COMM_WORLD,"Timestep %3D: step size = %g, time = %g, 2-norm error = %g, max norm error = %g\n",step,(double)dt,(double)time,(double)norm_2,(double)norm_max);CHKERRQ(ierr); 430 431 appctx->norm_2 += norm_2; 432 appctx->norm_max += norm_max; 433 434 dttol = .0001; 435 ierr = PetscOptionsGetReal(NULL,NULL,"-dttol",&dttol,NULL);CHKERRQ(ierr); 436 if (dt < dttol) { 437 dt *= .999; 438 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 439 } 440 441 /* 442 View a graph of the error 443 */ 444 ierr = VecView(appctx->solution,appctx->viewer1);CHKERRQ(ierr); 445 446 /* 447 Print debugging information if desired 448 */ 449 if (appctx->debug) { 450 ierr = PetscPrintf(PETSC_COMM_SELF,"Error vector\n");CHKERRQ(ierr); 451 ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 452 } 453 454 return 0; 455 } 456 /* --------------------------------------------------------------------- */ 457 /* 458 RHSMatrixHeat - User-provided routine to compute the right-hand-side 459 matrix for the heat equation. 460 461 Input Parameters: 462 ts - the TS context 463 t - current time 464 global_in - global input vector 465 dummy - optional user-defined context, as set by TSetRHSJacobian() 466 467 Output Parameters: 468 AA - Jacobian matrix 469 BB - optionally different preconditioning matrix 470 str - flag indicating matrix structure 471 472 Notes: 473 Recall that MatSetValues() uses 0-based row and column numbers 474 in Fortran as well as in C. 475 */ 476 PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx) 477 { 478 Mat A = AA; /* Jacobian matrix */ 479 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 480 PetscInt mstart = 0; 481 PetscInt mend = appctx->m; 482 PetscErrorCode ierr; 483 PetscInt i,idx[3]; 484 PetscScalar v[3],stwo = -2./(appctx->h*appctx->h),sone = -.5*stwo; 485 486 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 487 Compute entries for the locally owned part of the matrix 488 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 489 /* 490 Set matrix rows corresponding to boundary data 491 */ 492 493 mstart = 0; 494 v[0] = 1.0; 495 ierr = MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);CHKERRQ(ierr); 496 mstart++; 497 498 mend--; 499 v[0] = 1.0; 500 ierr = MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);CHKERRQ(ierr); 501 502 /* 503 Set matrix rows corresponding to interior data. We construct the 504 matrix one row at a time. 505 */ 506 v[0] = sone; v[1] = stwo; v[2] = sone; 507 for (i=mstart; i<mend; i++) { 508 idx[0] = i-1; idx[1] = i; idx[2] = i+1; 509 ierr = MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 510 } 511 512 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 513 Complete the matrix assembly process and set some options 514 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 515 /* 516 Assemble matrix, using the 2-step process: 517 MatAssemblyBegin(), MatAssemblyEnd() 518 Computations can be done while messages are in transition 519 by placing code between these two statements. 520 */ 521 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 522 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 523 524 /* 525 Set and option to indicate that we will never add a new nonzero location 526 to the matrix. If we do, it will generate an error. 527 */ 528 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 529 530 return 0; 531 } 532 533 PetscErrorCode IFunctionHeat(TS ts,PetscReal t,Vec X,Vec Xdot,Vec r,void *ctx) 534 { 535 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 536 PetscErrorCode ierr; 537 538 ierr = MatMult(appctx->A,X,r);CHKERRQ(ierr); 539 ierr = VecAYPX(r,-1.0,Xdot);CHKERRQ(ierr); 540 return 0; 541 } 542 543 PetscErrorCode IJacobianHeat(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal s,Mat A,Mat B,void *ctx) 544 { 545 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 546 PetscErrorCode ierr; 547 548 if (appctx->oshift == s) return 0; 549 ierr = MatCopy(appctx->A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 550 ierr = MatScale(A,-1);CHKERRQ(ierr); 551 ierr = MatShift(A,s);CHKERRQ(ierr); 552 ierr = MatCopy(A,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 553 appctx->oshift = s; 554 return 0; 555 } 556 557 /*TEST 558 559 test: 560 args: -nox -ts_type ssp -ts_dt 0.0005 561 562 test: 563 suffix: 2 564 args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1 565 566 test: 567 suffix: 3 568 args: -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason 569 filter: sed "s/ATOL/RTOL/g" 570 requires: !single 571 572 test: 573 suffix: 4 574 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason 575 filter: sed "s/ATOL/RTOL/g" 576 577 test: 578 suffix: 5 579 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs 580 filter: sed "s/ATOL/RTOL/g" 581 582 test: 583 requires: !single 584 suffix: pod_guess 585 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason 586 587 test: 588 requires: !single 589 suffix: pod_guess_Ainner 590 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -ksp_guess_pod_Ainner -pc_type none -ksp_converged_reason 591 592 test: 593 requires: !single 594 suffix: fischer_guess 595 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason 596 597 test: 598 requires: !single 599 suffix: fischer_guess_2 600 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 2,10 -pc_type none -ksp_converged_reason 601 602 test: 603 requires: !single 604 suffix: stringview 605 args: -nox -ts_type rosw -test_string_viewer 606 607 test: 608 requires: !single 609 suffix: stringview_euler 610 args: -nox -ts_type euler -test_string_viewer 611 612 TEST*/ 613