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
main(int argc,char ** argv)79 int main(int argc, char **argv)
80 {
81 AppCtx appctx; /* user-defined application context */
82 TS ts; /* timestepping context */
83 Mat A; /* matrix data structure */
84 Vec u; /* approximate solution vector */
85 PetscReal time_total_max = 100.0; /* default max total time */
86 PetscInt time_steps_max = 100; /* default max timesteps */
87 PetscDraw draw; /* drawing context */
88 PetscInt steps, m;
89 PetscMPIInt size;
90 PetscReal dt;
91 PetscBool flg, flg_string;
92
93 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94 Initialize program and set problem parameters
95 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96
97 PetscFunctionBeginUser;
98 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
99 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
100 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
101
102 m = 60;
103 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
104 PetscCall(PetscOptionsHasName(NULL, NULL, "-debug", &appctx.debug));
105 flg_string = PETSC_FALSE;
106 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_string_viewer", &flg_string, NULL));
107
108 appctx.m = m;
109 appctx.h = 1.0 / (m - 1.0);
110 appctx.norm_2 = 0.0;
111 appctx.norm_max = 0.0;
112
113 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving a linear TS problem on 1 processor\n"));
114
115 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116 Create vector data structures
117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118
119 /*
120 Create vector data structures for approximate and exact solutions
121 */
122 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &u));
123 PetscCall(VecDuplicate(u, &appctx.solution));
124
125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126 Set up displays to show graphs of the solution and error
127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128
129 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 380, 400, 160, &appctx.viewer1));
130 PetscCall(PetscViewerDrawGetDraw(appctx.viewer1, 0, &draw));
131 PetscCall(PetscDrawSetDoubleBuffer(draw));
132 PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, "", 80, 0, 400, 160, &appctx.viewer2));
133 PetscCall(PetscViewerDrawGetDraw(appctx.viewer2, 0, &draw));
134 PetscCall(PetscDrawSetDoubleBuffer(draw));
135
136 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137 Create timestepping solver context
138 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139
140 PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
141 PetscCall(TSSetProblemType(ts, TS_LINEAR));
142
143 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 Set optional user-defined monitoring routine
145 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146
147 if (!flg_string) PetscCall(TSMonitorSet(ts, Monitor, &appctx, NULL));
148
149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150
151 Create matrix data structure; set matrix evaluation routine.
152 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153
154 PetscCall(MatCreate(PETSC_COMM_SELF, &A));
155 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
156 PetscCall(MatSetFromOptions(A));
157 PetscCall(MatSetUp(A));
158
159 flg = PETSC_FALSE;
160 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_ifunc", &flg, NULL));
161 if (!flg) {
162 appctx.A = NULL;
163 PetscCall(PetscOptionsGetBool(NULL, NULL, "-time_dependent_rhs", &flg, NULL));
164 if (flg) {
165 /*
166 For linear problems with a time-dependent f(u,t) in the equation
167 u_t = f(u,t), the user provides the discretized right-hand side
168 as a time-dependent matrix.
169 */
170 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
171 PetscCall(TSSetRHSJacobian(ts, A, A, RHSMatrixHeat, &appctx));
172 } else {
173 /*
174 For linear problems with a time-independent f(u) in the equation
175 u_t = f(u), the user provides the discretized right-hand side
176 as a matrix only once, and then sets the special Jacobian evaluation
177 routine TSComputeRHSJacobianConstant() which will NOT recompute the Jacobian.
178 */
179 PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx));
180 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &appctx));
181 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &appctx));
182 }
183 } else {
184 Mat J;
185
186 PetscCall(RHSMatrixHeat(ts, 0.0, u, A, A, &appctx));
187 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &J));
188 PetscCall(TSSetIFunction(ts, NULL, IFunctionHeat, &appctx));
189 PetscCall(TSSetIJacobian(ts, J, J, IJacobianHeat, &appctx));
190 PetscCall(MatDestroy(&J));
191
192 PetscCall(PetscObjectReference((PetscObject)A));
193 appctx.A = A;
194 appctx.oshift = PETSC_MIN_REAL;
195 }
196 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197 Set solution vector and initial timestep
198 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199
200 dt = appctx.h * appctx.h / 2.0;
201 PetscCall(TSSetTimeStep(ts, dt));
202
203 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204 Customize timestepping solver:
205 - Set the solution method to be the Backward Euler method.
206 - Set timestepping duration info
207 Then set runtime options, which can override these defaults.
208 For example,
209 -ts_max_steps <maxsteps> -ts_max_time <maxtime>
210 to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
211 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212
213 PetscCall(TSSetMaxSteps(ts, time_steps_max));
214 PetscCall(TSSetMaxTime(ts, time_total_max));
215 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
216 PetscCall(TSSetFromOptions(ts));
217
218 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219 Solve the problem
220 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221
222 /*
223 Evaluate initial conditions
224 */
225 PetscCall(InitialConditions(u, &appctx));
226
227 /*
228 Run the timestepping solver
229 */
230 PetscCall(TSSolve(ts, u));
231 PetscCall(TSGetStepNumber(ts, &steps));
232
233 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234 View timestepping solver info
235 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
236
237 PetscCall(PetscPrintf(PETSC_COMM_SELF, "avg. error (2 norm) = %g, avg. error (max norm) = %g\n", (double)(appctx.norm_2 / steps), (double)(appctx.norm_max / steps)));
238 if (!flg_string) {
239 PetscCall(TSView(ts, PETSC_VIEWER_STDOUT_SELF));
240 } else {
241 PetscViewer stringviewer;
242 char string[512];
243 const char *outstring;
244
245 PetscCall(PetscViewerStringOpen(PETSC_COMM_WORLD, string, sizeof(string), &stringviewer));
246 PetscCall(TSView(ts, stringviewer));
247 PetscCall(PetscViewerStringGetStringRead(stringviewer, &outstring, NULL));
248 PetscCheck((char *)outstring == (char *)string, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "String returned from viewer does not equal original string");
249 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output from string viewer:%s\n", outstring));
250 PetscCall(PetscViewerDestroy(&stringviewer));
251 }
252
253 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254 Free work space. All PETSc objects should be destroyed when they
255 are no longer needed.
256 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
257
258 PetscCall(TSDestroy(&ts));
259 PetscCall(MatDestroy(&A));
260 PetscCall(VecDestroy(&u));
261 PetscCall(PetscViewerDestroy(&appctx.viewer1));
262 PetscCall(PetscViewerDestroy(&appctx.viewer2));
263 PetscCall(VecDestroy(&appctx.solution));
264 PetscCall(MatDestroy(&appctx.A));
265
266 /*
267 Always call PetscFinalize() before exiting a program. This routine
268 - finalizes the PETSc libraries as well as MPI
269 - provides summary and diagnostic information if certain runtime
270 options are chosen (e.g., -log_view).
271 */
272 PetscCall(PetscFinalize());
273 return 0;
274 }
275 /* --------------------------------------------------------------------- */
276 /*
277 InitialConditions - Computes the solution at the initial time.
278
279 Input Parameter:
280 u - uninitialized solution vector (global)
281 appctx - user-defined application context
282
283 Output Parameter:
284 u - vector with solution at initial time (global)
285 */
InitialConditions(Vec u,AppCtx * appctx)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 */
ExactSolution(PetscReal t,Vec solution,AppCtx * appctx)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 */
Monitor(TS ts,PetscInt step,PetscReal time,Vec u,PetscCtx ctx)382 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec u, PetscCtx ctx)
383 {
384 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
385 PetscReal norm_2, norm_max, dt, dttol;
386
387 PetscFunctionBeginUser;
388 /*
389 View a graph of the current iterate
390 */
391 PetscCall(VecView(u, appctx->viewer2));
392
393 /*
394 Compute the exact solution
395 */
396 PetscCall(ExactSolution(time, appctx->solution, appctx));
397
398 /*
399 Print debugging information if desired
400 */
401 if (appctx->debug) {
402 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Computed solution vector\n"));
403 PetscCall(VecView(u, PETSC_VIEWER_STDOUT_SELF));
404 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Exact solution vector\n"));
405 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
406 }
407
408 /*
409 Compute the 2-norm and max-norm of the error
410 */
411 PetscCall(VecAXPY(appctx->solution, -1.0, u));
412 PetscCall(VecNorm(appctx->solution, NORM_2, &norm_2));
413 norm_2 = PetscSqrtReal(appctx->h) * norm_2;
414 PetscCall(VecNorm(appctx->solution, NORM_MAX, &norm_max));
415
416 PetscCall(TSGetTimeStep(ts, &dt));
417 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep %3" PetscInt_FMT ": step size = %g, time = %g, 2-norm error = %g, max norm error = %g\n", step, (double)dt, (double)time, (double)norm_2, (double)norm_max));
418
419 appctx->norm_2 += norm_2;
420 appctx->norm_max += norm_max;
421
422 dttol = .0001;
423 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dttol", &dttol, NULL));
424 if (dt < dttol) {
425 dt *= .999;
426 PetscCall(TSSetTimeStep(ts, dt));
427 }
428
429 /*
430 View a graph of the error
431 */
432 PetscCall(VecView(appctx->solution, appctx->viewer1));
433
434 /*
435 Print debugging information if desired
436 */
437 if (appctx->debug) {
438 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error vector\n"));
439 PetscCall(VecView(appctx->solution, PETSC_VIEWER_STDOUT_SELF));
440 }
441 PetscFunctionReturn(PETSC_SUCCESS);
442 }
443 /* --------------------------------------------------------------------- */
444 /*
445 RHSMatrixHeat - User-provided routine to compute the right-hand-side
446 matrix for the heat equation.
447
448 Input Parameters:
449 ts - the TS context
450 t - current time
451 global_in - global input vector
452 dummy - optional user-defined context, as set by TSetRHSJacobian()
453
454 Output Parameters:
455 AA - Jacobian matrix
456 BB - optionally different matrix used to construct the preconditioner
457
458 Notes:
459 Recall that MatSetValues() uses 0-based row and column numbers
460 in Fortran as well as in C.
461 */
RHSMatrixHeat(TS ts,PetscReal t,Vec X,Mat AA,Mat BB,PetscCtx ctx)462 PetscErrorCode RHSMatrixHeat(TS ts, PetscReal t, Vec X, Mat AA, Mat BB, PetscCtx ctx)
463 {
464 Mat A = AA; /* Jacobian matrix */
465 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
466 PetscInt mstart = 0;
467 PetscInt mend = appctx->m;
468 PetscInt i, idx[3];
469 PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;
470
471 PetscFunctionBeginUser;
472 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
473 Compute entries for the locally owned part of the matrix
474 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
475 /*
476 Set matrix rows corresponding to boundary data
477 */
478
479 mstart = 0;
480 v[0] = 1.0;
481 PetscCall(MatSetValues(A, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
482 mstart++;
483
484 mend--;
485 v[0] = 1.0;
486 PetscCall(MatSetValues(A, 1, &mend, 1, &mend, v, INSERT_VALUES));
487
488 /*
489 Set matrix rows corresponding to interior data. We construct the
490 matrix one row at a time.
491 */
492 v[0] = sone;
493 v[1] = stwo;
494 v[2] = sone;
495 for (i = mstart; i < mend; i++) {
496 idx[0] = i - 1;
497 idx[1] = i;
498 idx[2] = i + 1;
499 PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
500 }
501
502 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
503 Complete the matrix assembly process and set some options
504 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
505 /*
506 Assemble matrix, using the 2-step process:
507 MatAssemblyBegin(), MatAssemblyEnd()
508 Computations can be done while messages are in transition
509 by placing code between these two statements.
510 */
511 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
512 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
513
514 /*
515 Set and option to indicate that we will never add a new nonzero location
516 to the matrix. If we do, it will generate an error.
517 */
518 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
519 PetscFunctionReturn(PETSC_SUCCESS);
520 }
521
IFunctionHeat(TS ts,PetscReal t,Vec X,Vec Xdot,Vec r,PetscCtx ctx)522 PetscErrorCode IFunctionHeat(TS ts, PetscReal t, Vec X, Vec Xdot, Vec r, PetscCtx ctx)
523 {
524 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
525
526 PetscFunctionBeginUser;
527 PetscCall(MatMult(appctx->A, X, r));
528 PetscCall(VecAYPX(r, -1.0, Xdot));
529 PetscFunctionReturn(PETSC_SUCCESS);
530 }
531
IJacobianHeat(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal s,Mat A,Mat B,PetscCtx ctx)532 PetscErrorCode IJacobianHeat(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal s, Mat A, Mat B, PetscCtx ctx)
533 {
534 AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
535
536 PetscFunctionBeginUser;
537 if (appctx->oshift == s) PetscFunctionReturn(PETSC_SUCCESS);
538 PetscCall(MatCopy(appctx->A, A, SAME_NONZERO_PATTERN));
539 PetscCall(MatScale(A, -1));
540 PetscCall(MatShift(A, s));
541 PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
542 appctx->oshift = s;
543 PetscFunctionReturn(PETSC_SUCCESS);
544 }
545
546 /*TEST
547
548 test:
549 args: -nox -ts_type ssp -ts_time_step 0.0005
550
551 test:
552 suffix: 2
553 args: -nox -ts_type ssp -ts_time_step 0.0005 -time_dependent_rhs 1
554
555 test:
556 suffix: 3
557 args: -nox -ts_type rosw -ts_max_steps 3 -ksp_converged_reason
558 filter: sed "s/ATOL/RTOL/g"
559 requires: !single
560
561 test:
562 suffix: 4
563 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason
564 filter: sed "s/ATOL/RTOL/g"
565
566 test:
567 suffix: 5
568 args: -nox -ts_type beuler -ts_max_steps 3 -ksp_converged_reason -time_dependent_rhs
569 filter: sed "s/ATOL/RTOL/g"
570
571 test:
572 requires: !single
573 suffix: pod_guess
574 args: -nox -ts_type beuler -use_ifunc -ts_time_step 0.0005 -ksp_guess_type pod -pc_type none -ksp_converged_reason
575
576 test:
577 requires: !single
578 suffix: pod_guess_Ainner
579 args: -nox -ts_type beuler -use_ifunc -ts_time_step 0.0005 -ksp_guess_type pod -ksp_guess_pod_Ainner -pc_type none -ksp_converged_reason
580
581 test:
582 requires: !single
583 suffix: fischer_guess
584 args: -nox -ts_type beuler -use_ifunc -ts_time_step 0.0005 -ksp_guess_type fischer -pc_type none -ksp_converged_reason
585
586 test:
587 requires: !single
588 suffix: fischer_guess_2
589 args: -nox -ts_type beuler -use_ifunc -ts_time_step 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 2,10 -pc_type none -ksp_converged_reason
590
591 test:
592 requires: !single
593 suffix: fischer_guess_3
594 args: -nox -ts_type beuler -use_ifunc -ts_time_step 0.0005 -ksp_guess_type fischer -ksp_guess_fischer_model 3,10 -pc_type none -ksp_converged_reason
595
596 test:
597 requires: !single
598 suffix: stringview
599 args: -nox -ts_type rosw -test_string_viewer
600
601 test:
602 requires: !single
603 suffix: stringview_euler
604 args: -nox -ts_type euler -test_string_viewer
605
606 TEST*/
607