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