xref: /petsc/src/ts/tutorials/ex4.c (revision a69119a591a03a9d906b29c0a4e9802e4d7c9795)
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   -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        mpiexec -n <procs> 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 uniprocessor version of this code is ts/tutorials/ex3.c
38 
39   ------------------------------------------------------------------------- */
40 
41 /*
42    Include "petscdmda.h" so that we can use distributed arrays (DMDAs) to manage
43    the parallel grid.  Include "petscts.h" so that we can use TS solvers.
44    Note that this file 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 <petscdm.h>
53 #include <petscdmda.h>
54 #include <petscts.h>
55 #include <petscdraw.h>
56 
57 /*
58    User-defined application context - contains data needed by the
59    application-provided call-back routines.
60 */
61 typedef struct {
62   MPI_Comm    comm;             /* communicator */
63   DM          da;               /* distributed array data structure */
64   Vec         localwork;        /* local ghosted work vector */
65   Vec         u_local;          /* local ghosted approximate solution vector */
66   Vec         solution;         /* global exact solution vector */
67   PetscInt    m;                /* total number of grid points */
68   PetscReal   h;                /* mesh width h = 1/(m-1) */
69   PetscBool   debug;            /* flag (1 indicates activation of debugging printouts) */
70   PetscViewer viewer1, viewer2; /* viewers for the solution and error */
71   PetscReal   norm_2, norm_max; /* error norms */
72 } AppCtx;
73 
74 /*
75    User-defined routines
76 */
77 extern PetscErrorCode InitialConditions(Vec, AppCtx *);
78 extern PetscErrorCode RHSMatrixHeat(TS, PetscReal, Vec, Mat, Mat, void *);
79 extern PetscErrorCode RHSFunctionHeat(TS, PetscReal, Vec, Vec, void *);
80 extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
81 extern PetscErrorCode ExactSolution(PetscReal, Vec, AppCtx *);
82 
83 int main(int argc, char **argv) {
84   AppCtx        appctx;               /* user-defined application context */
85   TS            ts;                   /* timestepping context */
86   Mat           A;                    /* matrix data structure */
87   Vec           u;                    /* approximate solution vector */
88   PetscReal     time_total_max = 1.0; /* default max total time */
89   PetscInt      time_steps_max = 100; /* default max timesteps */
90   PetscDraw     draw;                 /* drawing context */
91   PetscInt      steps, m;
92   PetscMPIInt   size;
93   PetscReal     dt, ftime;
94   PetscBool     flg;
95   TSProblemType tsproblem = TS_LINEAR;
96 
97   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98      Initialize program and set problem parameters
99      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100 
101   PetscFunctionBeginUser;
102   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
103   appctx.comm = PETSC_COMM_WORLD;
104 
105   m = 60;
106   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
107   PetscCall(PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug));
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   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
114   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solving a linear TS problem, number of processors = %d\n", size));
115 
116   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117      Create vector data structures
118      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119   /*
120      Create distributed array (DMDA) to manage parallel grid and vectors
121      and to set up the ghost point communication pattern.  There are M
122      total grid values spread equally among all the processors.
123   */
124 
125   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m, 1, 1, NULL, &appctx.da));
126   PetscCall(DMSetFromOptions(appctx.da));
127   PetscCall(DMSetUp(appctx.da));
128 
129   /*
130      Extract global and local vectors from DMDA; we use these to store the
131      approximate solution.  Then duplicate these for remaining vectors that
132      have the same types.
133   */
134   PetscCall(DMCreateGlobalVector(appctx.da, &u));
135   PetscCall(DMCreateLocalVector(appctx.da, &appctx.u_local));
136 
137   /*
138      Create local work vector for use in evaluating right-hand-side function;
139      create global work vector for storing exact solution.
140   */
141   PetscCall(VecDuplicate(appctx.u_local, &appctx.localwork));
142   PetscCall(VecDuplicate(u, &appctx.solution));
143 
144   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145      Set up displays to show graphs of the solution and error
146      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147 
148   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "", 80, 380, 400, 160, &appctx.viewer1));
149   PetscCall(PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw));
150   PetscCall(PetscDrawSetDoubleBuffer(draw));
151   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, "", 80, 0, 400, 160, &appctx.viewer2));
152   PetscCall(PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw));
153   PetscCall(PetscDrawSetDoubleBuffer(draw));
154 
155   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156      Create timestepping solver context
157      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158 
159   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
160 
161   flg = PETSC_FALSE;
162   PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonlinear", &flg, NULL));
163   PetscCall(TSSetProblemType(ts, flg ? TS_NONLINEAR : TS_LINEAR));
164 
165   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166      Set optional user-defined monitoring routine
167      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168   PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
169 
170   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 
172      Create matrix data structure; set matrix evaluation routine.
173      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174 
175   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
176   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
177   PetscCall(MatSetFromOptions(A));
178   PetscCall(MatSetUp(A));
179 
180   flg = PETSC_FALSE;
181   PetscCall(PetscOptionsGetBool(NULL, NULL, "-time_dependent_rhs", &flg, NULL));
182   if (flg) {
183     /*
184        For linear problems with a time-dependent f(u,t) in the equation
185        u_t = f(u,t), the user provides the discretized right-hand-side
186        as a time-dependent matrix.
187     */
188     PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
189     PetscCall(TSSetRHSJacobian(ts, A, A, RHSMatrixHeat, &appctx));
190   } else {
191     /*
192        For linear problems with a time-independent f(u) in the equation
193        u_t = f(u), the user provides the discretized right-hand-side
194        as a matrix only once, and then sets a null matrix evaluation
195        routine.
196     */
197     PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx));
198     PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
199     PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &appctx));
200   }
201 
202   if (tsproblem == TS_NONLINEAR) {
203     SNES snes;
204     PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionHeat, &appctx));
205     PetscCall(TSGetSNES(ts, &snes));
206     PetscCall(SNESSetJacobian(snes, NULL, NULL, SNESComputeJacobianDefault, NULL));
207   }
208 
209   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210      Set solution vector and initial timestep
211      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212 
213   dt = appctx.h * appctx.h / 2.0;
214   PetscCall(TSSetTimeStep(ts, dt));
215   PetscCall(TSSetSolution(ts, u));
216 
217   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218      Customize timestepping solver:
219        - Set the solution method to be the Backward Euler method.
220        - Set timestepping duration info
221      Then set runtime options, which can override these defaults.
222      For example,
223           -ts_max_steps <maxsteps> -ts_max_time <maxtime>
224      to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
225      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
226 
227   PetscCall(TSSetMaxSteps(ts, time_steps_max));
228   PetscCall(TSSetMaxTime(ts, time_total_max));
229   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
230   PetscCall(TSSetFromOptions(ts));
231 
232   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233      Solve the problem
234      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235 
236   /*
237      Evaluate initial conditions
238   */
239   PetscCall(InitialConditions(u, &appctx));
240 
241   /*
242      Run the timestepping solver
243   */
244   PetscCall(TSSolve(ts, u));
245   PetscCall(TSGetSolveTime(ts, &ftime));
246   PetscCall(TSGetStepNumber(ts, &steps));
247 
248   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249      View timestepping solver info
250      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
251   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Total timesteps %" PetscInt_FMT ", Final time %g\n", steps, (double)ftime));
252   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Avg. error (2 norm) = %g Avg. error (max norm) = %g\n", (double)(appctx.norm_2 / steps), (double)(appctx.norm_max / steps)));
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.localwork));
265   PetscCall(VecDestroy(&appctx.solution));
266   PetscCall(VecDestroy(&appctx.u_local));
267   PetscCall(DMDestroy(&appctx.da));
268 
269   /*
270      Always call PetscFinalize() before exiting a program.  This routine
271        - finalizes the PETSc libraries as well as MPI
272        - provides summary and diagnostic information if certain runtime
273          options are chosen (e.g., -log_view).
274   */
275   PetscCall(PetscFinalize());
276   return 0;
277 }
278 /* --------------------------------------------------------------------- */
279 /*
280    InitialConditions - Computes the solution at the initial time.
281 
282    Input Parameter:
283    u - uninitialized solution vector (global)
284    appctx - user-defined application context
285 
286    Output Parameter:
287    u - vector with solution at initial time (global)
288 */
289 PetscErrorCode InitialConditions(Vec u, AppCtx *appctx) {
290   PetscScalar *u_localptr, h = appctx->h;
291   PetscInt     i, mybase, myend;
292 
293   /*
294      Determine starting point of each processor's range of
295      grid values.
296   */
297   PetscCall(VecGetOwnershipRange(u, &mybase, &myend));
298 
299   /*
300     Get a pointer to vector data.
301     - For default PETSc vectors, VecGetArray() returns a pointer to
302       the data array.  Otherwise, the routine is implementation dependent.
303     - You MUST call VecRestoreArray() when you no longer need access to
304       the array.
305     - Note that the Fortran interface to VecGetArray() differs from the
306       C version.  See the users manual for details.
307   */
308   PetscCall(VecGetArray(u, &u_localptr));
309 
310   /*
311      We initialize the solution array by simply writing the solution
312      directly into the array locations.  Alternatively, we could use
313      VecSetValues() or VecSetValuesLocal().
314   */
315   for (i = mybase; i < myend; i++) u_localptr[i - mybase] = PetscSinScalar(PETSC_PI * i * 6. * h) + 3. * PetscSinScalar(PETSC_PI * i * 2. * h);
316 
317   /*
318      Restore vector
319   */
320   PetscCall(VecRestoreArray(u, &u_localptr));
321 
322   /*
323      Print debugging information if desired
324   */
325   if (appctx->debug) {
326     PetscCall(PetscPrintf(appctx->comm, "initial guess vector\n"));
327     PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
328   }
329 
330   return 0;
331 }
332 /* --------------------------------------------------------------------- */
333 /*
334    ExactSolution - Computes the exact solution at a given time.
335 
336    Input Parameters:
337    t - current time
338    solution - vector in which exact solution will be computed
339    appctx - user-defined application context
340 
341    Output Parameter:
342    solution - vector with the newly computed exact solution
343 */
344 PetscErrorCode ExactSolution(PetscReal t, Vec solution, AppCtx *appctx) {
345   PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2;
346   PetscInt     i, mybase, myend;
347 
348   /*
349      Determine starting and ending points of each processor's
350      range of grid values
351   */
352   PetscCall(VecGetOwnershipRange(solution, &mybase, &myend));
353 
354   /*
355      Get a pointer to vector data.
356   */
357   PetscCall(VecGetArray(solution, &s_localptr));
358 
359   /*
360      Simply write the solution directly into the array locations.
361      Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
362   */
363   ex1 = PetscExpReal(-36. * PETSC_PI * PETSC_PI * t);
364   ex2 = PetscExpReal(-4. * PETSC_PI * PETSC_PI * t);
365   sc1 = PETSC_PI * 6. * h;
366   sc2 = PETSC_PI * 2. * h;
367   for (i = mybase; i < myend; i++) s_localptr[i - mybase] = PetscSinScalar(sc1 * (PetscReal)i) * ex1 + 3. * PetscSinScalar(sc2 * (PetscReal)i) * ex2;
368 
369   /*
370      Restore vector
371   */
372   PetscCall(VecRestoreArray(solution, &s_localptr));
373   return 0;
374 }
375 /* --------------------------------------------------------------------- */
376 /*
377    Monitor - User-provided routine to monitor the solution computed at
378    each timestep.  This example plots the solution and computes the
379    error in two different norms.
380 
381    Input Parameters:
382    ts     - the timestep context
383    step   - the count of the current step (with 0 meaning the
384              initial condition)
385    time   - the current time
386    u      - the solution at this timestep
387    ctx    - the user-provided context for this monitoring routine.
388             In this case we use the application context which contains
389             information about the problem size, workspace and the exact
390             solution.
391 */
392 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx) {
393   AppCtx   *appctx = (AppCtx *)ctx; /* user-defined application context */
394   PetscReal norm_2, norm_max;
395 
396   /*
397      View a graph of the current iterate
398   */
399   PetscCall(VecView(u, appctx->viewer2));
400 
401   /*
402      Compute the exact solution
403   */
404   PetscCall(ExactSolution(time, appctx->solution, appctx));
405 
406   /*
407      Print debugging information if desired
408   */
409   if (appctx->debug) {
410     PetscCall(PetscPrintf(appctx->comm, "Computed solution vector\n"));
411     PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
412     PetscCall(PetscPrintf(appctx->comm, "Exact solution vector\n"));
413     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_WORLD));
414   }
415 
416   /*
417      Compute the 2-norm and max-norm of the error
418   */
419   PetscCall(VecAXPY(appctx->solution, -1.0, u));
420   PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
421   norm_2 = PetscSqrtReal(appctx->h) * norm_2;
422   PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
423   if (norm_2 < 1e-14) norm_2 = 0;
424   if (norm_max < 1e-14) norm_max = 0;
425 
426   /*
427      PetscPrintf() causes only the first processor in this
428      communicator to print the timestep information.
429   */
430   PetscCall(PetscPrintf(appctx->comm, "Timestep %" PetscInt_FMT ": time = %g 2-norm error = %g max norm error = %g\n", step, (double)time, (double)norm_2, (double)norm_max));
431   appctx->norm_2 += norm_2;
432   appctx->norm_max += norm_max;
433 
434   /*
435      View a graph of the error
436   */
437   PetscCall(VecView(appctx->solution, appctx->viewer1));
438 
439   /*
440      Print debugging information if desired
441   */
442   if (appctx->debug) {
443     PetscCall(PetscPrintf(appctx->comm, "Error vector\n"));
444     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_WORLD));
445   }
446 
447   return 0;
448 }
449 
450 /* --------------------------------------------------------------------- */
451 /*
452    RHSMatrixHeat - User-provided routine to compute the right-hand-side
453    matrix for the heat equation.
454 
455    Input Parameters:
456    ts - the TS context
457    t - current time
458    global_in - global input vector
459    dummy - optional user-defined context, as set by TSetRHSJacobian()
460 
461    Output Parameters:
462    AA - Jacobian matrix
463    BB - optionally different preconditioning matrix
464    str - flag indicating matrix structure
465 
466   Notes:
467   RHSMatrixHeat computes entries for the locally owned part of the system.
468    - Currently, all PETSc parallel matrix formats are partitioned by
469      contiguous chunks of rows across the processors.
470    - Each processor needs to insert only elements that it owns
471      locally (but any non-local elements will be sent to the
472      appropriate processor during matrix assembly).
473    - Always specify global row and columns of matrix entries when
474      using MatSetValues(); we could alternatively use MatSetValuesLocal().
475    - Here, we set all entries for a particular row at once.
476    - Note that MatSetValues() uses 0-based row and column numbers
477      in Fortran as well as in C.
478 */
479 PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx) {
480   Mat         A      = AA;            /* Jacobian matrix */
481   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
482   PetscInt    i, mstart, mend, idx[3];
483   PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
484 
485   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
486      Compute entries for the locally owned part of the matrix
487      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
488 
489   PetscCall(MatGetOwnershipRange(A, &mstart, &mend));
490 
491   /*
492      Set matrix rows corresponding to boundary data
493   */
494 
495   if (mstart == 0) { /* first processor only */
496     v[0] = 1.0;
497     PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
498     mstart++;
499   }
500 
501   if (mend == appctx->m) { /* last processor only */
502     mend--;
503     v[0] = 1.0;
504     PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
505   }
506 
507   /*
508      Set matrix rows corresponding to interior data.  We construct the
509      matrix one row at a time.
510   */
511   v[0] = sone;
512   v[1] = stwo;
513   v[2] = sone;
514   for (i = mstart; i < mend; i++) {
515     idx[0] = i - 1;
516     idx[1] = i;
517     idx[2] = i + 1;
518     PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
519   }
520 
521   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
522      Complete the matrix assembly process and set some options
523      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
524   /*
525      Assemble matrix, using the 2-step process:
526        MatAssemblyBegin(), MatAssemblyEnd()
527      Computations can be done while messages are in transition
528      by placing code between these two statements.
529   */
530   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
531   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
532 
533   /*
534      Set and option to indicate that we will never add a new nonzero location
535      to the matrix. If we do, it will generate an error.
536   */
537   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
538 
539   return 0;
540 }
541 
542 PetscErrorCode RHSFunctionHeat(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx) {
543   Mat A;
544 
545   PetscFunctionBeginUser;
546   PetscCall(TSGetRHSJacobian(ts, &A, NULL, NULL, &ctx));
547   PetscCall(RHSMatrixHeat(ts, t, globalin, A, NULL, ctx));
548   /* PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); */
549   PetscCall(MatMult(A, globalin, globalout));
550   PetscFunctionReturn(0);
551 }
552 
553 /*TEST
554 
555     test:
556       args: -ts_view -nox
557 
558     test:
559       suffix: 2
560       args: -ts_view -nox
561       nsize: 3
562 
563     test:
564       suffix: 3
565       args: -ts_view -nox -nonlinear
566 
567     test:
568       suffix: 4
569       args: -ts_view -nox -nonlinear
570       nsize: 3
571       timeoutfactor: 3
572 
573     test:
574       suffix: sundials
575       requires: sundials2
576       args: -nox -ts_type sundials -ts_max_steps 5 -nonlinear
577       nsize: 4
578 
579     test:
580       suffix: sundials_dense
581       requires: sundials2
582       args: -nox -ts_type sundials -ts_sundials_use_dense -ts_max_steps 5 -nonlinear
583       nsize: 1
584 
585 TEST*/
586