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