Actual source code: ex3.c

  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";

  9: /* ------------------------------------------------------------------------

 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.

 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>

 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)

 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)

 37:     The parallel version of this code is ts/tutorials/ex4.c

 39:   ------------------------------------------------------------------------- */

 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: */

 51: #include <petscts.h>
 52: #include <petscdraw.h>

 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;

 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 *);

 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;

 93:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:      Initialize program and set problem parameters
 95:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 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!");

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));

108:   appctx.m        = m;
109:   appctx.h        = 1.0 / (m - 1.0);
110:   appctx.norm_2   = 0.0;
111:   appctx.norm_max = 0.0;

113:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n"));

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:      Create vector data structures
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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));

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:      Set up displays to show graphs of the solution and error
127:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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));

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:      Create timestepping solver context
138:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

140:   PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
141:   PetscCall(TSSetProblemType(ts, TS_LINEAR));

143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:      Set optional user-defined monitoring routine
145:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

147:   if (!flg_string) PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

151:      Create matrix data structure; set matrix evaluation routine.
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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));

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;

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));

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:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

200:   dt = appctx.h * appctx.h / 2.0;
201:   PetscCall(TSSetTimeStep(ts, dt));

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:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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));

218:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219:      Solve the problem
220:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

222:   /*
223:      Evaluate initial conditions
224:   */
225:   PetscCall(InitialConditions(u, &appctx));

227:   /*
228:      Run the timestepping solver
229:   */
230:   PetscCall(TSSolve(ts, u));
231:   PetscCall(TSGetStepNumber(ts, &steps));

233:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234:      View timestepping solver info
235:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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;

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:   }

253:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254:      Free work space.  All PETSc objects should be destroyed when they
255:      are no longer needed.
256:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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));

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.

279:    Input Parameter:
280:    u - uninitialized solution vector (global)
281:    appctx - user-defined application context

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;

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));

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);

310:   /*
311:      Restore vector
312:   */
313:   PetscCall(VecRestoreArrayWrite(u, &u_localptr));

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:   }

323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }
325: /* --------------------------------------------------------------------- */
326: /*
327:    ExactSolution - Computes the exact solution at a given time.

329:    Input Parameters:
330:    t - current time
331:    solution - vector in which exact solution will be computed
332:    appctx - user-defined application context

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;

342:   PetscFunctionBeginUser;
343:   /*
344:      Get a pointer to vector data.
345:   */
346:   PetscCall(VecGetArrayWrite(solution, &s_localptr));

348:   /*
349:      Simply write the solution directly into the array locations.
350:      Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
351:   */
352:   ex1 = PetscExpScalar(-36. * PETSC_PI * PETSC_PI * tc);
353:   ex2 = PetscExpScalar(-4. * PETSC_PI * PETSC_PI * tc);
354:   sc1 = PETSC_PI * 6. * h;
355:   sc2 = PETSC_PI * 2. * h;
356:   for (i = 0; i < appctx->m; i++) s_localptr[i] = PetscSinScalar(sc1 * (PetscReal)i) * ex1 + 3. * PetscSinScalar(sc2 * (PetscReal)i) * ex2;

358:   /*
359:      Restore vector
360:   */
361:   PetscCall(VecRestoreArrayWrite(solution, &s_localptr));
362:   PetscFunctionReturn(PETSC_SUCCESS);
363: }
364: /* --------------------------------------------------------------------- */
365: /*
366:    Monitor - User-provided routine to monitor the solution computed at
367:    each timestep.  This example plots the solution and computes the
368:    error in two different norms.

370:    This example also demonstrates changing the timestep via TSSetTimeStep().

372:    Input Parameters:
373:    ts     - the timestep context
374:    step   - the count of the current step (with 0 meaning the
375:              initial condition)
376:    time   - the current time
377:    u      - the solution at this timestep
378:    ctx    - the user-provided context for this monitoring routine.
379:             In this case we use the application context which contains
380:             information about the problem size, workspace and the exact
381:             solution.
382: */
383: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx)
384: {
385:   AppCtx   *appctx = (AppCtx *)ctx; /* user-defined application context */
386:   PetscReal norm_2, norm_max, dt, dttol;

388:   PetscFunctionBeginUser;
389:   /*
390:      View a graph of the current iterate
391:   */
392:   PetscCall(VecView(u, appctx->viewer2));

394:   /*
395:      Compute the exact solution
396:   */
397:   PetscCall(ExactSolution(time, appctx->solution, appctx));

399:   /*
400:      Print debugging information if desired
401:   */
402:   if (appctx->debug) {
403:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n"));
404:     PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
405:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
406:     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
407:   }

409:   /*
410:      Compute the 2-norm and max-norm of the error
411:   */
412:   PetscCall(VecAXPY(appctx->solution, -1.0, u));
413:   PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
414:   norm_2 = PetscSqrtReal(appctx->h) * norm_2;
415:   PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));

417:   PetscCall(TSGetTimeStep(ts, &dt));
418:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep %3" PetscInt_FMT ": step size = %g, time = %g, 2-norm error = %g, max norm error = %g\n", step, (double)dt, (double)time, (double)norm_2, (double)norm_max));

420:   appctx->norm_2 += norm_2;
421:   appctx->norm_max += norm_max;

423:   dttol = .0001;
424:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, NULL));
425:   if (dt < dttol) {
426:     dt *= .999;
427:     PetscCall(TSSetTimeStep(ts, dt));
428:   }

430:   /*
431:      View a graph of the error
432:   */
433:   PetscCall(VecView(appctx->solution, appctx->viewer1));

435:   /*
436:      Print debugging information if desired
437:   */
438:   if (appctx->debug) {
439:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
440:     PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
441:   }

443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }
445: /* --------------------------------------------------------------------- */
446: /*
447:    RHSMatrixHeat - User-provided routine to compute the right-hand-side
448:    matrix for the heat equation.

450:    Input Parameters:
451:    ts - the TS context
452:    t - current time
453:    global_in - global input vector
454:    dummy - optional user-defined context, as set by TSetRHSJacobian()

456:    Output Parameters:
457:    AA - Jacobian matrix
458:    BB - optionally different preconditioning matrix
459:    str - flag indicating matrix structure

461:    Notes:
462:    Recall that MatSetValues() uses 0-based row and column numbers
463:    in Fortran as well as in C.
464: */
465: PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)
466: {
467:   Mat         A      = AA;            /* Jacobian matrix */
468:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
469:   PetscInt    mstart = 0;
470:   PetscInt    mend   = appctx->m;
471:   PetscInt    i, idx[3];
472:   PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;

474:   PetscFunctionBeginUser;
475:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
476:      Compute entries for the locally owned part of the matrix
477:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
478:   /*
479:      Set matrix rows corresponding to boundary data
480:   */

482:   mstart = 0;
483:   v[0]   = 1.0;
484:   PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
485:   mstart++;

487:   mend--;
488:   v[0] = 1.0;
489:   PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));

491:   /*
492:      Set matrix rows corresponding to interior data.  We construct the
493:      matrix one row at a time.
494:   */
495:   v[0] = sone;
496:   v[1] = stwo;
497:   v[2] = sone;
498:   for (i = mstart; i < mend; i++) {
499:     idx[0] = i - 1;
500:     idx[1] = i;
501:     idx[2] = i + 1;
502:     PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
503:   }

505:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
506:      Complete the matrix assembly process and set some options
507:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
508:   /*
509:      Assemble matrix, using the 2-step process:
510:        MatAssemblyBegin(), MatAssemblyEnd()
511:      Computations can be done while messages are in transition
512:      by placing code between these two statements.
513:   */
514:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
515:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

517:   /*
518:      Set and option to indicate that we will never add a new nonzero location
519:      to the matrix. If we do, it will generate an error.
520:   */
521:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));

523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

526: PetscErrorCode IFunctionHeat(TS ts, PetscReal t, Vec X, Vec Xdot, Vec r, void *ctx)
527: {
528:   AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */

530:   PetscFunctionBeginUser;
531:   PetscCall(MatMult(appctx->A, X, r));
532:   PetscCall(VecAYPX(r, -1.0, Xdot));
533:   PetscFunctionReturn(PETSC_SUCCESS);
534: }

536: PetscErrorCode IJacobianHeat(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal s, Mat A, Mat B, void *ctx)
537: {
538:   AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */

540:   PetscFunctionBeginUser;
541:   if (appctx->oshift == s) PetscFunctionReturn(PETSC_SUCCESS);
542:   PetscCall(MatCopy(appctx->A, A, SAME_NONZERO_PATTERN));
543:   PetscCall(MatScale(A, -1));
544:   PetscCall(MatShift(A, s));
545:   PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
546:   appctx->oshift = s;
547:   PetscFunctionReturn(PETSC_SUCCESS);
548: }

550: /*TEST

552:     test:
553:       args: -nox -ts_type ssp -ts_dt 0.0005

555:     test:
556:       suffix: 2
557:       args: -nox -ts_type ssp -ts_dt 0.0005 -time_dependent_rhs 1

559:     test:
560:       suffix: 3
561:       args:  -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason
562:       filter: sed "s/ATOL/RTOL/g"
563:       requires: !single

565:     test:
566:       suffix: 4
567:       args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason
568:       filter: sed "s/ATOL/RTOL/g"

570:     test:
571:       suffix: 5
572:       args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs
573:       filter: sed "s/ATOL/RTOL/g"

575:     test:
576:       requires: !single
577:       suffix: pod_guess
578:       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason

580:     test:
581:       requires: !single
582:       suffix: pod_guess_Ainner
583:       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type pod -ksp_guess_pod_Ainner -pc_type none -ksp_converged_reason

585:     test:
586:       requires: !single
587:       suffix: fischer_guess
588:       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason

590:     test:
591:       requires: !single
592:       suffix: fischer_guess_2
593:       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 2,10 -pc_type none -ksp_converged_reason

595:     test:
596:       requires: !single
597:       suffix: fischer_guess_3
598:       args: -nox -ts_type beuler -use_ifunc -ts_dt 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 3,10 -pc_type none -ksp_converged_reason

600:     test:
601:       requires: !single
602:       suffix: stringview
603:       args: -nox -ts_type rosw -test_string_viewer

605:     test:
606:       requires: !single
607:       suffix: stringview_euler
608:       args: -nox -ts_type euler -test_string_viewer

610: TEST*/