xref: /petsc/src/ts/tutorials/ex3.c (revision 4ad8454beace47809662cdae21ee081016eaa39a)
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, (char *)0, 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, void *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 preconditioning matrix
457    str - flag indicating matrix structure
458 
459    Notes:
460    Recall that MatSetValues() uses 0-based row and column numbers
461    in Fortran as well as in C.
462 */
463 PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)
464 {
465   Mat         A      = AA;            /* Jacobian matrix */
466   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
467   PetscInt    mstart = 0;
468   PetscInt    mend   = appctx->m;
469   PetscInt    i, idx[3];
470   PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
471 
472   PetscFunctionBeginUser;
473   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
474      Compute entries for the locally owned part of the matrix
475      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
476   /*
477      Set matrix rows corresponding to boundary data
478   */
479 
480   mstart = 0;
481   v[0]   = 1.0;
482   PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
483   mstart++;
484 
485   mend--;
486   v[0] = 1.0;
487   PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
488 
489   /*
490      Set matrix rows corresponding to interior data.  We construct the
491      matrix one row at a time.
492   */
493   v[0] = sone;
494   v[1] = stwo;
495   v[2] = sone;
496   for (i = mstart; i < mend; i++) {
497     idx[0] = i - 1;
498     idx[1] = i;
499     idx[2] = i + 1;
500     PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
501   }
502 
503   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
504      Complete the matrix assembly process and set some options
505      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
506   /*
507      Assemble matrix, using the 2-step process:
508        MatAssemblyBegin(), MatAssemblyEnd()
509      Computations can be done while messages are in transition
510      by placing code between these two statements.
511   */
512   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
513   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
514 
515   /*
516      Set and option to indicate that we will never add a new nonzero location
517      to the matrix. If we do, it will generate an error.
518   */
519   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
520   PetscFunctionReturn(PETSC_SUCCESS);
521 }
522 
523 PetscErrorCode IFunctionHeat(TS ts, PetscReal t, Vec X, Vec Xdot, Vec r, void *ctx)
524 {
525   AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
526 
527   PetscFunctionBeginUser;
528   PetscCall(MatMult(appctx->A, X, r));
529   PetscCall(VecAYPX(r, -1.0, Xdot));
530   PetscFunctionReturn(PETSC_SUCCESS);
531 }
532 
533 PetscErrorCode IJacobianHeat(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal s, Mat A, Mat B, void *ctx)
534 {
535   AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
536 
537   PetscFunctionBeginUser;
538   if (appctx->oshift == s) PetscFunctionReturn(PETSC_SUCCESS);
539   PetscCall(MatCopy(appctx->A, A, SAME_NONZERO_PATTERN));
540   PetscCall(MatScale(A, -1));
541   PetscCall(MatShift(A, s));
542   PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
543   appctx->oshift = s;
544   PetscFunctionReturn(PETSC_SUCCESS);
545 }
546 
547 /*TEST
548 
549     test:
550       args: -nox -ts_type ssp -ts_dt 0.0005
551 
552     test:
553       suffix: 2
554       args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1
555 
556     test:
557       suffix: 3
558       args: -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason
559       filter: sed "s/ATOL/RTOL/g"
560       requires: !single
561 
562     test:
563       suffix: 4
564       args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason
565       filter: sed "s/ATOL/RTOL/g"
566 
567     test:
568       suffix: 5
569       args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs
570       filter: sed "s/ATOL/RTOL/g"
571 
572     test:
573       requires: !single
574       suffix: pod_guess
575       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason
576 
577     test:
578       requires: !single
579       suffix: pod_guess_Ainner
580       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
581 
582     test:
583       requires: !single
584       suffix: fischer_guess
585       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason
586 
587     test:
588       requires: !single
589       suffix: fischer_guess_2
590       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
591 
592     test:
593       requires: !single
594       suffix: fischer_guess_3
595       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
596 
597     test:
598       requires: !single
599       suffix: stringview
600       args: -nox -ts_type rosw -test_string_viewer
601 
602     test:
603       requires: !single
604       suffix: stringview_euler
605       args: -nox -ts_type euler -test_string_viewer
606 
607 TEST*/
608