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 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 108 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 109 110 m = 60; 111 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 112 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-debug",&appctx.debug)); 113 flg_string = PETSC_FALSE; 114 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_string_viewer",&flg_string,NULL)); 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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Solving a linear TS problem on 1 processor\n")); 122 123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124 Create vector data structures 125 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 126 127 /* 128 Create vector data structures for approximate and exact solutions 129 */ 130 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,m,&u)); 131 CHKERRQ(VecDuplicate(u,&appctx.solution)); 132 133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 134 Set up displays to show graphs of the solution and error 135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 136 137 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,380,400,160,&appctx.viewer1)); 138 CHKERRQ(PetscViewerDrawGetDraw(appctx.viewer1,0,&draw)); 139 CHKERRQ(PetscDrawSetDoubleBuffer(draw)); 140 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_SELF,0,"",80,0,400,160,&appctx.viewer2)); 141 CHKERRQ(PetscViewerDrawGetDraw(appctx.viewer2,0,&draw)); 142 CHKERRQ(PetscDrawSetDoubleBuffer(draw)); 143 144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 Create timestepping solver context 146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 147 148 CHKERRQ(TSCreate(PETSC_COMM_SELF,&ts)); 149 CHKERRQ(TSSetProblemType(ts,TS_LINEAR)); 150 151 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152 Set optional user-defined monitoring routine 153 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 154 155 if (!flg_string) { 156 CHKERRQ(TSMonitorSet(ts,Monitor,&appctx,NULL)); 157 } 158 159 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160 161 Create matrix data structure; set matrix evaluation routine. 162 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 163 164 CHKERRQ(MatCreate(PETSC_COMM_SELF,&A)); 165 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m)); 166 CHKERRQ(MatSetFromOptions(A)); 167 CHKERRQ(MatSetUp(A)); 168 169 flg = PETSC_FALSE; 170 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-use_ifunc",&flg,NULL)); 171 if (!flg) { 172 appctx.A = NULL; 173 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-time_dependent_rhs",&flg,NULL)); 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 CHKERRQ(TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx)); 181 CHKERRQ(TSSetRHSJacobian(ts,A,A,RHSMatrixHeat,&appctx)); 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 CHKERRQ(RHSMatrixHeat(ts,0.0,u,A,A,&appctx)); 190 CHKERRQ(TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx)); 191 CHKERRQ(TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx)); 192 } 193 } else { 194 Mat J; 195 196 CHKERRQ(RHSMatrixHeat(ts,0.0,u,A,A,&appctx)); 197 CHKERRQ(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&J)); 198 CHKERRQ(TSSetIFunction(ts,NULL,IFunctionHeat,&appctx)); 199 CHKERRQ(TSSetIJacobian(ts,J,J,IJacobianHeat,&appctx)); 200 CHKERRQ(MatDestroy(&J)); 201 202 CHKERRQ(PetscObjectReference((PetscObject)A)); 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 CHKERRQ(TSSetTimeStep(ts,dt)); 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 CHKERRQ(TSSetMaxSteps(ts,time_steps_max)); 224 CHKERRQ(TSSetMaxTime(ts,time_total_max)); 225 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 226 CHKERRQ(TSSetFromOptions(ts)); 227 228 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 229 Solve the problem 230 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 231 232 /* 233 Evaluate initial conditions 234 */ 235 CHKERRQ(InitialConditions(u,&appctx)); 236 237 /* 238 Run the timestepping solver 239 */ 240 CHKERRQ(TSSolve(ts,u)); 241 CHKERRQ(TSGetStepNumber(ts,&steps)); 242 243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244 View timestepping solver info 245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 246 247 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))); 248 if (!flg_string) { 249 CHKERRQ(TSView(ts,PETSC_VIEWER_STDOUT_SELF)); 250 } else { 251 PetscViewer stringviewer; 252 char string[512]; 253 const char *outstring; 254 255 CHKERRQ(PetscViewerStringOpen(PETSC_COMM_WORLD,string,sizeof(string),&stringviewer)); 256 CHKERRQ(TSView(ts,stringviewer)); 257 CHKERRQ(PetscViewerStringGetStringRead(stringviewer,&outstring,NULL)); 258 PetscCheck((char*)outstring == (char*)string,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"String returned from viewer does not equal original string"); 259 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Output from string viewer:%s\n",outstring)); 260 CHKERRQ(PetscViewerDestroy(&stringviewer)); 261 } 262 263 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 264 Free work space. All PETSc objects should be destroyed when they 265 are no longer needed. 266 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 267 268 CHKERRQ(TSDestroy(&ts)); 269 CHKERRQ(MatDestroy(&A)); 270 CHKERRQ(VecDestroy(&u)); 271 CHKERRQ(PetscViewerDestroy(&appctx.viewer1)); 272 CHKERRQ(PetscViewerDestroy(&appctx.viewer2)); 273 CHKERRQ(VecDestroy(&appctx.solution)); 274 CHKERRQ(MatDestroy(&appctx.A)); 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 PetscInt i; 300 301 /* 302 Get a pointer to vector data. 303 - For default PETSc vectors, VecGetArray() returns a pointer to 304 the data array. Otherwise, the routine is implementation dependent. 305 - You MUST call VecRestoreArray() when you no longer need access to 306 the array. 307 - Note that the Fortran interface to VecGetArray() differs from the 308 C version. See the users manual for details. 309 */ 310 CHKERRQ(VecGetArrayWrite(u,&u_localptr)); 311 312 /* 313 We initialize the solution array by simply writing the solution 314 directly into the array locations. Alternatively, we could use 315 VecSetValues() or VecSetValuesLocal(). 316 */ 317 for (i=0; i<appctx->m; i++) u_localptr[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h); 318 319 /* 320 Restore vector 321 */ 322 CHKERRQ(VecRestoreArrayWrite(u,&u_localptr)); 323 324 /* 325 Print debugging information if desired 326 */ 327 if (appctx->debug) { 328 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Initial guess vector\n")); 329 CHKERRQ(VecView(u,PETSC_VIEWER_STDOUT_SELF)); 330 } 331 332 return 0; 333 } 334 /* --------------------------------------------------------------------- */ 335 /* 336 ExactSolution - Computes the exact solution at a given time. 337 338 Input Parameters: 339 t - current time 340 solution - vector in which exact solution will be computed 341 appctx - user-defined application context 342 343 Output Parameter: 344 solution - vector with the newly computed exact solution 345 */ 346 PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx) 347 { 348 PetscScalar *s_localptr,h = appctx->h,ex1,ex2,sc1,sc2,tc = t; 349 PetscInt i; 350 351 /* 352 Get a pointer to vector data. 353 */ 354 CHKERRQ(VecGetArrayWrite(solution,&s_localptr)); 355 356 /* 357 Simply write the solution directly into the array locations. 358 Alternatively, we culd use VecSetValues() or VecSetValuesLocal(). 359 */ 360 ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*tc); 361 ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*tc); 362 sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h; 363 for (i=0; i<appctx->m; i++) s_localptr[i] = PetscSinScalar(sc1*(PetscReal)i)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i)*ex2; 364 365 /* 366 Restore vector 367 */ 368 CHKERRQ(VecRestoreArrayWrite(solution,&s_localptr)); 369 return 0; 370 } 371 /* --------------------------------------------------------------------- */ 372 /* 373 Monitor - User-provided routine to monitor the solution computed at 374 each timestep. This example plots the solution and computes the 375 error in two different norms. 376 377 This example also demonstrates changing the timestep via TSSetTimeStep(). 378 379 Input Parameters: 380 ts - the timestep context 381 step - the count of the current step (with 0 meaning the 382 initial condition) 383 time - the current time 384 u - the solution at this timestep 385 ctx - the user-provided context for this monitoring routine. 386 In this case we use the application context which contains 387 information about the problem size, workspace and the exact 388 solution. 389 */ 390 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx) 391 { 392 AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */ 393 PetscReal norm_2,norm_max,dt,dttol; 394 395 /* 396 View a graph of the current iterate 397 */ 398 CHKERRQ(VecView(u,appctx->viewer2)); 399 400 /* 401 Compute the exact solution 402 */ 403 CHKERRQ(ExactSolution(time,appctx->solution,appctx)); 404 405 /* 406 Print debugging information if desired 407 */ 408 if (appctx->debug) { 409 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Computed solution vector\n")); 410 CHKERRQ(VecView(u,PETSC_VIEWER_STDOUT_SELF)); 411 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n")); 412 CHKERRQ(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 413 } 414 415 /* 416 Compute the 2-norm and max-norm of the error 417 */ 418 CHKERRQ(VecAXPY(appctx->solution,-1.0,u)); 419 CHKERRQ(VecNorm(appctx->solution,NORM_2,&norm_2)); 420 norm_2 = PetscSqrtReal(appctx->h)*norm_2; 421 CHKERRQ(VecNorm(appctx->solution,NORM_MAX,&norm_max)); 422 423 CHKERRQ(TSGetTimeStep(ts,&dt)); 424 CHKERRQ(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)); 425 426 appctx->norm_2 += norm_2; 427 appctx->norm_max += norm_max; 428 429 dttol = .0001; 430 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-dttol",&dttol,NULL)); 431 if (dt < dttol) { 432 dt *= .999; 433 CHKERRQ(TSSetTimeStep(ts,dt)); 434 } 435 436 /* 437 View a graph of the error 438 */ 439 CHKERRQ(VecView(appctx->solution,appctx->viewer1)); 440 441 /* 442 Print debugging information if desired 443 */ 444 if (appctx->debug) { 445 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error vector\n")); 446 CHKERRQ(VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF)); 447 } 448 449 return 0; 450 } 451 /* --------------------------------------------------------------------- */ 452 /* 453 RHSMatrixHeat - User-provided routine to compute the right-hand-side 454 matrix for the heat equation. 455 456 Input Parameters: 457 ts - the TS context 458 t - current time 459 global_in - global input vector 460 dummy - optional user-defined context, as set by TSetRHSJacobian() 461 462 Output Parameters: 463 AA - Jacobian matrix 464 BB - optionally different preconditioning matrix 465 str - flag indicating matrix structure 466 467 Notes: 468 Recall that MatSetValues() uses 0-based row and column numbers 469 in Fortran as well as in C. 470 */ 471 PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,void *ctx) 472 { 473 Mat A = AA; /* Jacobian matrix */ 474 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 475 PetscInt mstart = 0; 476 PetscInt mend = appctx->m; 477 PetscInt i,idx[3]; 478 PetscScalar v[3],stwo = -2./(appctx->h*appctx->h),sone = -.5*stwo; 479 480 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 481 Compute entries for the locally owned part of the matrix 482 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 483 /* 484 Set matrix rows corresponding to boundary data 485 */ 486 487 mstart = 0; 488 v[0] = 1.0; 489 CHKERRQ(MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES)); 490 mstart++; 491 492 mend--; 493 v[0] = 1.0; 494 CHKERRQ(MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES)); 495 496 /* 497 Set matrix rows corresponding to interior data. We construct the 498 matrix one row at a time. 499 */ 500 v[0] = sone; v[1] = stwo; v[2] = sone; 501 for (i=mstart; i<mend; i++) { 502 idx[0] = i-1; idx[1] = i; idx[2] = i+1; 503 CHKERRQ(MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES)); 504 } 505 506 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 507 Complete the matrix assembly process and set some options 508 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 509 /* 510 Assemble matrix, using the 2-step process: 511 MatAssemblyBegin(), MatAssemblyEnd() 512 Computations can be done while messages are in transition 513 by placing code between these two statements. 514 */ 515 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 516 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 517 518 /* 519 Set and option to indicate that we will never add a new nonzero location 520 to the matrix. If we do, it will generate an error. 521 */ 522 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 523 524 return 0; 525 } 526 527 PetscErrorCode IFunctionHeat(TS ts,PetscReal t,Vec X,Vec Xdot,Vec r,void *ctx) 528 { 529 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 530 531 CHKERRQ(MatMult(appctx->A,X,r)); 532 CHKERRQ(VecAYPX(r,-1.0,Xdot)); 533 return 0; 534 } 535 536 PetscErrorCode IJacobianHeat(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal s,Mat A,Mat B,void *ctx) 537 { 538 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 539 540 if (appctx->oshift == s) return 0; 541 CHKERRQ(MatCopy(appctx->A,A,SAME_NONZERO_PATTERN)); 542 CHKERRQ(MatScale(A,-1)); 543 CHKERRQ(MatShift(A,s)); 544 CHKERRQ(MatCopy(A,B,SAME_NONZERO_PATTERN)); 545 appctx->oshift = s; 546 return 0; 547 } 548 549 /*TEST 550 551 test: 552 args: -nox -ts_type ssp -ts_dt 0.0005 553 554 test: 555 suffix: 2 556 args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1 557 558 test: 559 suffix: 3 560 args: -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason 561 filter: sed "s/ATOL/RTOL/g" 562 requires: !single 563 564 test: 565 suffix: 4 566 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason 567 filter: sed "s/ATOL/RTOL/g" 568 569 test: 570 suffix: 5 571 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs 572 filter: sed "s/ATOL/RTOL/g" 573 574 test: 575 requires: !single 576 suffix: pod_guess 577 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason 578 579 test: 580 requires: !single 581 suffix: pod_guess_Ainner 582 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 583 584 test: 585 requires: !single 586 suffix: fischer_guess 587 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason 588 589 test: 590 requires: !single 591 suffix: fischer_guess_2 592 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 593 594 test: 595 requires: !single 596 suffix: fischer_guess_3 597 args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 3,10 -pc_type none -ksp_converged_reason 598 599 test: 600 requires: !single 601 suffix: stringview 602 args: -nox -ts_type rosw -test_string_viewer 603 604 test: 605 requires: !single 606 suffix: stringview_euler 607 args: -nox -ts_type euler -test_string_viewer 608 609 TEST*/ 610