xref: /petsc/src/ts/tutorials/ex3.c (revision d5b43468fb8780a8feea140ccd6fa3e6a50411cc)
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   /*
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 
323   return 0;
324 }
325 /* --------------------------------------------------------------------- */
326 /*
327    ExactSolution - Computes the exact solution at a given time.
328 
329    Input Parameters:
330    t - current time
331    solution - vector in which exact solution will be computed
332    appctx - user-defined application context
333 
334    Output Parameter:
335    solution - vector with the newly computed exact solution
336 */
337 PetscErrorCode ExactSolution(PetscReal t, Vec solution, AppCtx *appctx)
338 {
339   PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2, tc = t;
340   PetscInt     i;
341 
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   return 0;
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   /*
388      View a graph of the current iterate
389   */
390   PetscCall(VecView(u, appctx->viewer2));
391 
392   /*
393      Compute the exact solution
394   */
395   PetscCall(ExactSolution(time, appctx->solution, appctx));
396 
397   /*
398      Print debugging information if desired
399   */
400   if (appctx->debug) {
401     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n"));
402     PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
403     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
404     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
405   }
406 
407   /*
408      Compute the 2-norm and max-norm of the error
409   */
410   PetscCall(VecAXPY(appctx->solution, -1.0, u));
411   PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
412   norm_2 = PetscSqrtReal(appctx->h) * norm_2;
413   PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
414 
415   PetscCall(TSGetTimeStep(ts, &dt));
416   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));
417 
418   appctx->norm_2 += norm_2;
419   appctx->norm_max += norm_max;
420 
421   dttol = .0001;
422   PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, NULL));
423   if (dt < dttol) {
424     dt *= .999;
425     PetscCall(TSSetTimeStep(ts, dt));
426   }
427 
428   /*
429      View a graph of the error
430   */
431   PetscCall(VecView(appctx->solution, appctx->viewer1));
432 
433   /*
434      Print debugging information if desired
435   */
436   if (appctx->debug) {
437     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
438     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
439   }
440 
441   return 0;
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   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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 
520   return 0;
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   PetscCall(MatMult(appctx->A, X, r));
528   PetscCall(VecAYPX(r, -1.0, Xdot));
529   return 0;
530 }
531 
532 PetscErrorCode IJacobianHeat(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal s, Mat A, Mat B, void *ctx)
533 {
534   AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
535 
536   if (appctx->oshift == s) return 0;
537   PetscCall(MatCopy(appctx->A, A, SAME_NONZERO_PATTERN));
538   PetscCall(MatScale(A, -1));
539   PetscCall(MatShift(A, s));
540   PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
541   appctx->oshift = s;
542   return 0;
543 }
544 
545 /*TEST
546 
547     test:
548       args: -nox -ts_type ssp -ts_dt 0.0005
549 
550     test:
551       suffix: 2
552       args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1
553 
554     test:
555       suffix: 3
556       args:  -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason
557       filter: sed "s/ATOL/RTOL/g"
558       requires: !single
559 
560     test:
561       suffix: 4
562       args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason
563       filter: sed "s/ATOL/RTOL/g"
564 
565     test:
566       suffix: 5
567       args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs
568       filter: sed "s/ATOL/RTOL/g"
569 
570     test:
571       requires: !single
572       suffix: pod_guess
573       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason
574 
575     test:
576       requires: !single
577       suffix: pod_guess_Ainner
578       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
579 
580     test:
581       requires: !single
582       suffix: fischer_guess
583       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason
584 
585     test:
586       requires: !single
587       suffix: fischer_guess_2
588       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
589 
590     test:
591       requires: !single
592       suffix: fischer_guess_3
593       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
594 
595     test:
596       requires: !single
597       suffix: stringview
598       args: -nox -ts_type rosw -test_string_viewer
599 
600     test:
601       requires: !single
602       suffix: stringview_euler
603       args: -nox -ts_type euler -test_string_viewer
604 
605 TEST*/
606