1 #pragma once
2
3 #include <petscts.h>
4 #include <petsc/private/petscimpl.h>
5
6 /* SUBMANSEC = TS */
7
8 /*
9 Timesteping context.
10 General DAE: F(t,U,U_t) = 0, required Jacobian is G'(U) where G(U) = F(t,U,U0+a*U)
11 General ODE: U_t = F(t,U) <-- the right-hand-side function
12 Linear ODE: U_t = A(t) U <-- the right-hand-side matrix
13 Linear (no time) ODE: U_t = A U <-- the right-hand-side matrix
14 */
15
16 /*
17 Maximum number of monitors you can run with a single TS
18 */
19 #define MAXTSMONITORS 10
20
21 PETSC_EXTERN PetscBool TSRegisterAllCalled;
22 PETSC_EXTERN PetscErrorCode TSRegisterAll(void);
23 PETSC_EXTERN PetscErrorCode TSAdaptRegisterAll(void);
24
25 PETSC_EXTERN PetscErrorCode TSRKRegisterAll(void);
26 PETSC_EXTERN PetscErrorCode TSMPRKRegisterAll(void);
27 PETSC_EXTERN PetscErrorCode TSARKIMEXRegisterAll(void);
28 PETSC_EXTERN PetscErrorCode TSRosWRegisterAll(void);
29 PETSC_EXTERN PetscErrorCode TSGLLERegisterAll(void);
30 PETSC_EXTERN PetscErrorCode TSGLLEAdaptRegisterAll(void);
31 PETSC_EXTERN PetscErrorCode TSIRKRegisterAll(void);
32
33 typedef struct _TSOps *TSOps;
34
35 struct _TSOps {
36 PetscErrorCode (*snesfunction)(SNES, Vec, Vec, TS);
37 PetscErrorCode (*snesjacobian)(SNES, Vec, Mat, Mat, TS);
38 PetscErrorCode (*setup)(TS);
39 PetscErrorCode (*step)(TS);
40 PetscErrorCode (*solve)(TS);
41 PetscErrorCode (*interpolate)(TS, PetscReal, Vec);
42 PetscErrorCode (*evaluatewlte)(TS, NormType, PetscInt *, PetscReal *);
43 PetscErrorCode (*evaluatestep)(TS, PetscInt, Vec, PetscBool *);
44 PetscErrorCode (*setfromoptions)(TS, PetscOptionItems);
45 PetscErrorCode (*destroy)(TS);
46 PetscErrorCode (*view)(TS, PetscViewer);
47 PetscErrorCode (*reset)(TS);
48 PetscErrorCode (*linearstability)(TS, PetscReal, PetscReal, PetscReal *, PetscReal *);
49 PetscErrorCode (*load)(TS, PetscViewer);
50 PetscErrorCode (*rollback)(TS);
51 PetscErrorCode (*getstages)(TS, PetscInt *, Vec *[]);
52 PetscErrorCode (*adjointstep)(TS);
53 PetscErrorCode (*adjointsetup)(TS);
54 PetscErrorCode (*adjointreset)(TS);
55 PetscErrorCode (*adjointintegral)(TS);
56 PetscErrorCode (*forwardsetup)(TS);
57 PetscErrorCode (*forwardreset)(TS);
58 PetscErrorCode (*forwardstep)(TS);
59 PetscErrorCode (*forwardintegral)(TS);
60 PetscErrorCode (*forwardgetstages)(TS, PetscInt *, Mat *[]);
61 PetscErrorCode (*getsolutioncomponents)(TS, PetscInt *, Vec *);
62 PetscErrorCode (*getauxsolution)(TS, Vec *);
63 PetscErrorCode (*gettimeerror)(TS, PetscInt, Vec *);
64 PetscErrorCode (*settimeerror)(TS, Vec);
65 PetscErrorCode (*startingmethod)(TS);
66 PetscErrorCode (*initcondition)(TS, Vec);
67 PetscErrorCode (*exacterror)(TS, Vec, Vec);
68 PetscErrorCode (*resizeregister)(TS, PetscBool);
69 };
70
71 /*
72 TSEvent - Abstract object to handle event monitoring
73 */
74 typedef struct _n_TSEvent *TSEvent;
75
76 typedef struct _TSTrajectoryOps *TSTrajectoryOps;
77
78 struct _TSTrajectoryOps {
79 PetscErrorCode (*view)(TSTrajectory, PetscViewer);
80 PetscErrorCode (*reset)(TSTrajectory);
81 PetscErrorCode (*destroy)(TSTrajectory);
82 PetscErrorCode (*set)(TSTrajectory, TS, PetscInt, PetscReal, Vec);
83 PetscErrorCode (*get)(TSTrajectory, TS, PetscInt, PetscReal *);
84 PetscErrorCode (*setfromoptions)(TSTrajectory, PetscOptionItems);
85 PetscErrorCode (*setup)(TSTrajectory, TS);
86 };
87
88 /* TSHistory is an helper object that allows inquiring
89 the TSTrajectory by time and not by the step number only */
90 typedef struct _n_TSHistory *TSHistory;
91
92 struct _p_TSTrajectory {
93 PETSCHEADER(struct _TSTrajectoryOps);
94 TSHistory tsh; /* associates times to unique step ids */
95 /* stores necessary data to reconstruct states and derivatives via Lagrangian interpolation */
96 struct {
97 PetscInt order; /* interpolation order */
98 Vec *W; /* work vectors */
99 PetscScalar *L; /* workspace for Lagrange basis */
100 PetscReal *T; /* Lagrange times (stored) */
101 Vec *WW; /* just an array of pointers */
102 PetscBool *TT; /* workspace for Lagrange */
103 PetscReal *TW; /* Lagrange times (workspace) */
104
105 /* caching */
106 PetscBool caching;
107 struct {
108 PetscObjectId id;
109 PetscObjectState state;
110 PetscReal time;
111 PetscInt step;
112 } Ucached;
113 struct {
114 PetscObjectId id;
115 PetscObjectState state;
116 PetscReal time;
117 PetscInt step;
118 } Udotcached;
119 } lag;
120 Vec U, Udot; /* used by TSTrajectory{Get|Restore}UpdatedHistoryVecs */
121 PetscBool usehistory; /* whether to use TSHistory */
122 PetscBool solution_only; /* whether we dump just the solution or also the stages */
123 PetscBool adjoint_solve_mode; /* whether we will use the Trajectory inside a TSAdjointSolve() or not */
124 PetscViewer monitor;
125 PetscBool setupcalled; /* true if setup has been called */
126 PetscInt recomps; /* counter for recomputations in the adjoint run */
127 PetscInt diskreads, diskwrites; /* counters for disk checkpoint reads and writes */
128 char **names; /* the name of each variable; each process has only the local names */
129 PetscBool keepfiles; /* keep the files generated during the run after the run is complete */
130 char *dirname, *filetemplate; /* directory name and file name template for disk checkpoints */
131 char *dirfiletemplate; /* complete directory and file name template for disk checkpoints */
132 PetscErrorCode (*transform)(void *, Vec, Vec *);
133 PetscCtxDestroyFn *transformdestroy;
134 void *transformctx;
135 void *data;
136 };
137
138 typedef struct _TS_RHSSplitLink *TS_RHSSplitLink;
139 struct _TS_RHSSplitLink {
140 TS ts;
141 char *splitname;
142 IS is;
143 TS_RHSSplitLink next;
144 PetscLogEvent event;
145 };
146
147 typedef struct _TS_EvaluationTimes *TSEvaluationTimes;
148 struct _TS_EvaluationTimes {
149 PetscInt num_time_points; /* number of time points */
150 PetscReal *time_points; /* array of the time span */
151 PetscReal reltol; /* relative tolerance for span point detection */
152 PetscReal abstol; /* absolute tolerance for span point detection */
153 PetscReal worktol; /* the ultimate tolerance (variable), maintained within a single TS time step for consistency */
154 PetscInt time_point_idx; /* index of the time_point to be reached next */
155 PetscInt sol_idx; /* index into sol_vecs and sol_times */
156 Vec *sol_vecs; /* array of the solutions at the specified time points */
157 PetscReal *sol_times; /* array of times that sol_vecs was taken at */
158 };
159
160 struct _p_TS {
161 PETSCHEADER(struct _TSOps);
162 TSProblemType problem_type;
163 TSEquationType equation_type;
164
165 DM dm;
166 Vec vec_sol; /* solution vector in first and second order equations */
167 Vec vec_sol0; /* solution vector at the beginning of the step */
168 Vec vec_dot; /* time derivative vector in second order equations */
169 TSAdapt adapt;
170 TSAdaptType default_adapt_type;
171 TSEvent event;
172
173 /* ---------------- Resize ---------------------*/
174 PetscBool resizerollback;
175 PetscObjectList resizetransferobjs;
176
177 /* ---------------- User (or PETSc) Provided stuff ---------------------*/
178 PetscErrorCode (*monitor[MAXTSMONITORS])(TS, PetscInt, PetscReal, Vec, void *);
179 PetscCtxDestroyFn *monitordestroy[MAXTSMONITORS];
180 void *monitorcontext[MAXTSMONITORS];
181 PetscInt numbermonitors;
182 PetscErrorCode (*adjointmonitor[MAXTSMONITORS])(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *);
183 PetscCtxDestroyFn *adjointmonitordestroy[MAXTSMONITORS];
184 void *adjointmonitorcontext[MAXTSMONITORS];
185 PetscInt numberadjointmonitors;
186
187 PetscErrorCode (*prestep)(TS);
188 PetscErrorCode (*prestage)(TS, PetscReal);
189 PetscErrorCode (*poststage)(TS, PetscReal, PetscInt, Vec *);
190 PetscErrorCode (*postevaluate)(TS);
191 PetscErrorCode (*poststep)(TS);
192 PetscErrorCode (*functiondomainerror)(TS, PetscReal, Vec, PetscBool *);
193 PetscErrorCode (*resizesetup)(TS, PetscInt, PetscReal, Vec, PetscBool *, void *);
194 PetscErrorCode (*resizetransfer)(TS, PetscInt, Vec[], Vec[], void *);
195 void *resizectx;
196
197 /* ---------------------- Sensitivity Analysis support -----------------*/
198 TSTrajectory trajectory; /* All solutions are kept here for the entire time integration process */
199 Vec *vecs_sensi; /* one vector for each cost function */
200 Vec *vecs_sensip;
201 PetscInt numcost; /* number of cost functions */
202 Vec vec_costintegral;
203 PetscBool adjointsetupcalled;
204 PetscInt adjoint_steps;
205 PetscInt adjoint_max_steps;
206 PetscBool adjoint_solve; /* immediately call TSAdjointSolve() after TSSolve() is complete */
207 PetscBool costintegralfwd; /* cost integral is evaluated in the forward run if true */
208 Vec vec_costintegrand; /* workspace for Adjoint computations */
209 Mat Jacp, Jacprhs;
210 void *ijacobianpctx, *rhsjacobianpctx;
211 void *costintegrandctx;
212 Vec *vecs_drdu;
213 Vec *vecs_drdp;
214 Vec vec_drdu_col, vec_drdp_col;
215
216 /* first-order adjoint */
217 PetscErrorCode (*rhsjacobianp)(TS, PetscReal, Vec, Mat, void *);
218 PetscErrorCode (*ijacobianp)(TS, PetscReal, Vec, Vec, PetscReal, Mat, void *);
219 PetscErrorCode (*costintegrand)(TS, PetscReal, Vec, Vec, void *);
220 PetscErrorCode (*drdufunction)(TS, PetscReal, Vec, Vec *, void *);
221 PetscErrorCode (*drdpfunction)(TS, PetscReal, Vec, Vec *, void *);
222
223 /* second-order adjoint */
224 Vec *vecs_sensi2;
225 Vec *vecs_sensi2p;
226 Vec vec_dir; /* directional vector for optimization */
227 Vec *vecs_fuu, *vecs_guu;
228 Vec *vecs_fup, *vecs_gup;
229 Vec *vecs_fpu, *vecs_gpu;
230 Vec *vecs_fpp, *vecs_gpp;
231 void *ihessianproductctx, *rhshessianproductctx;
232 PetscErrorCode (*ihessianproduct_fuu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
233 PetscErrorCode (*ihessianproduct_fup)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
234 PetscErrorCode (*ihessianproduct_fpu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
235 PetscErrorCode (*ihessianproduct_fpp)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
236 PetscErrorCode (*rhshessianproduct_guu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
237 PetscErrorCode (*rhshessianproduct_gup)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
238 PetscErrorCode (*rhshessianproduct_gpu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
239 PetscErrorCode (*rhshessianproduct_gpp)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
240
241 /* specific to forward sensitivity analysis */
242 Mat mat_sensip; /* matrix storing forward sensitivities */
243 Vec vec_sensip_col; /* space for a column of the sensip matrix */
244 Vec *vecs_integral_sensip; /* one vector for each integral */
245 PetscInt num_parameters;
246 PetscInt num_initialvalues;
247 void *vecsrhsjacobianpctx;
248 PetscBool forwardsetupcalled;
249 PetscBool forward_solve;
250 PetscErrorCode (*vecsrhsjacobianp)(TS, PetscReal, Vec, Vec *, void *);
251
252 /* ---------------------- IMEX support ---------------------------------*/
253 /* These extra slots are only used when the user provides both Implicit and RHS */
254 Mat Arhs; /* Right hand side matrix */
255 Mat Brhs; /* Right hand side matrix used to construct the preconditioner */
256 Vec Frhs; /* Right hand side function value */
257
258 /* This is a general caching scheme to avoid recomputing the Jacobian at a place that has been previously been evaluated.
259 * The present use case is that TSComputeRHSFunctionLinear() evaluates the Jacobian once and we don't want it to be immeditely re-evaluated.
260 */
261 struct {
262 PetscReal time; /* The time at which the matrices were last evaluated */
263 PetscObjectId Xid; /* Unique ID of solution vector at which the Jacobian was last evaluated */
264 PetscObjectState Xstate; /* State of the solution vector */
265 MatStructure mstructure; /* The structure returned */
266 /* Flag to unshift Jacobian before calling the IJacobian or RHSJacobian functions. This is useful
267 * if the user would like to reuse (part of) the Jacobian from the last evaluation. */
268 PetscBool reuse;
269 PetscReal scale, shift;
270 } rhsjacobian;
271
272 struct {
273 PetscReal shift; /* The derivative of the lhs wrt to Xdot */
274 } ijacobian;
275
276 MatStructure axpy_pattern; /* information about the nonzero pattern of the RHS Jacobian in reference to the implicit Jacobian */
277 /* --------------------Nonlinear Iteration------------------------------*/
278 SNES snes;
279 PetscBool usessnes; /* Flag set by each TSType to indicate if the type actually uses a SNES;
280 this works around the design flaw that a SNES is ALWAYS created with TS even when it is not needed.*/
281 PetscInt ksp_its; /* total number of linear solver iterations */
282 PetscInt snes_its; /* total number of nonlinear solver iterations */
283 PetscInt num_snes_failures;
284 PetscInt max_snes_failures;
285
286 /* --- Logging --- */
287 PetscInt ifuncs, rhsfuncs, ijacs, rhsjacs;
288
289 /* --- Data that is unique to each particular solver --- */
290 PetscBool setupcalled; /* true if setup has been called */
291 void *data; /* implementationspecific data */
292 void *ctx; /* user context */
293
294 PetscBool steprollback; /* flag to indicate that the step was rolled back */
295 PetscBool steprestart; /* flag to indicate that the timestepper has to discard any history and restart */
296 PetscBool stepresize; /* flag to indicate that the discretization was resized */
297 PetscInt steps; /* steps taken so far in all successive calls to TSSolve() */
298 PetscReal ptime; /* time at the start of the current step (stage time is internal if it exists) */
299 PetscReal time_step; /* current time increment */
300 PetscReal time_step0; /* proposed time increment at the beginning of the step */
301 PetscReal initial_time_step; /* proposed time step at start of TSSolve(), actual time step used may differ */
302 PetscReal ptime_prev; /* time at the start of the previous step */
303 PetscReal ptime_prev_rollback; /* time at the start of the 2nd previous step to recover from rollback */
304 PetscReal solvetime; /* time at the conclusion of TSSolve() */
305 PetscBool stifflyaccurate; /* flag to indicate that the method is stiffly accurate */
306
307 TSConvergedReason reason;
308 PetscBool errorifstepfailed;
309 PetscInt reject, max_reject;
310 TSExactFinalTimeOption exact_final_time;
311
312 PetscObjectParameterDeclare(PetscReal, rtol); /* Relative and absolute tolerance for local truncation error */
313 PetscObjectParameterDeclare(PetscReal, atol);
314 PetscObjectParameterDeclare(PetscReal, max_time); /* max time allowed */
315 PetscObjectParameterDeclare(PetscInt, max_steps); /* maximum time-step number to execute until (possibly with nonzero starting value) */
316 PetscObjectParameterDeclare(PetscInt, run_steps); /* maximum number of time steps for TSSolve to take on each call */
317 Vec vatol, vrtol; /* Relative and absolute tolerance in vector form */
318 PetscReal cfltime, cfltime_local;
319 PetscInt start_step; /* step number at start of current run */
320
321 PetscBool testjacobian;
322 PetscBool testjacobiantranspose;
323 /* ------------------- Default work-area management ------------------ */
324 PetscInt nwork;
325 Vec *work;
326
327 /* ---------------------- RHS splitting support ---------------------------------*/
328 PetscInt num_rhs_splits;
329 TS_RHSSplitLink tsrhssplit;
330 PetscBool use_splitrhsfunction;
331 SNES snesrhssplit;
332
333 /* ---------------------- Quadrature integration support ---------------------------------*/
334 TS quadraturets;
335
336 /* ---------------------- Time span support ---------------------------------*/
337 TSEvaluationTimes eval_times;
338 };
339
340 struct _TSAdaptOps {
341 PetscErrorCode (*choose)(TSAdapt, TS, PetscReal, PetscInt *, PetscReal *, PetscBool *, PetscReal *, PetscReal *, PetscReal *);
342 PetscErrorCode (*destroy)(TSAdapt);
343 PetscErrorCode (*reset)(TSAdapt);
344 PetscErrorCode (*view)(TSAdapt, PetscViewer);
345 PetscErrorCode (*setfromoptions)(TSAdapt, PetscOptionItems);
346 PetscErrorCode (*load)(TSAdapt, PetscViewer);
347 };
348
349 struct _p_TSAdapt {
350 PETSCHEADER(struct _TSAdaptOps);
351 void *data;
352 PetscErrorCode (*checkstage)(TSAdapt, TS, PetscReal, Vec, PetscBool *);
353 struct {
354 PetscInt n; /* number of candidate schemes, including the one currently in use */
355 PetscBool inuse_set; /* the current scheme has been set */
356 const char *name[16]; /* name of the scheme */
357 PetscInt order[16]; /* classical order of each scheme */
358 PetscInt stageorder[16]; /* stage order of each scheme */
359 PetscReal ccfl[16]; /* stability limit relative to explicit Euler */
360 PetscReal cost[16]; /* relative measure of the amount of work required for each scheme */
361 } candidates;
362 PetscBool always_accept;
363 PetscReal safety; /* safety factor relative to target error/stability goal */
364 PetscReal reject_safety; /* extra safety factor if the last step was rejected */
365 PetscReal clip[2]; /* admissible time step decrease/increase factors */
366 PetscReal dt_min, dt_max; /* admissible minimum and maximum time step */
367 PetscReal ignore_max; /* minimum value of the solution to be considered by the adaptor */
368 PetscBool glee_use_local; /* GLEE adaptor uses global or local error */
369 PetscReal scale_solve_failed; /* scale step by this factor if solver (linear or nonlinear) fails. */
370 PetscReal matchstepfac[2]; /* factors to control the behaviour of matchstep */
371 NormType wnormtype;
372 PetscViewer monitor;
373 PetscInt timestepjustdecreased_delay; /* number of timesteps after a decrease in the timestep before the timestep can be increased */
374 PetscInt timestepjustdecreased;
375 PetscReal dt_eval_times_cached; /* time step before hitting a TS evaluation time point */
376 };
377
378 typedef struct _p_DMTS *DMTS;
379 typedef struct _DMTSOps *DMTSOps;
380 struct _DMTSOps {
381 TSRHSFunctionFn *rhsfunction;
382 TSRHSJacobianFn *rhsjacobian;
383
384 TSIFunctionFn *ifunction;
385 PetscErrorCode (*ifunctionview)(void *, PetscViewer);
386 PetscErrorCode (*ifunctionload)(void **, PetscViewer);
387
388 TSIJacobianFn *ijacobian;
389 PetscErrorCode (*ijacobianview)(void *, PetscViewer);
390 PetscErrorCode (*ijacobianload)(void **, PetscViewer);
391
392 TSI2FunctionFn *i2function;
393 TSI2JacobianFn *i2jacobian;
394
395 TSTransientVariableFn *transientvar;
396
397 TSSolutionFn *solution;
398 TSForcingFn *forcing;
399
400 PetscErrorCode (*destroy)(DMTS);
401 PetscErrorCode (*duplicate)(DMTS, DMTS);
402 };
403
404 /*S
405 DMTS - Object held by a `DM` that contains all the callback functions and their contexts needed by a `TS`
406
407 Level: developer
408
409 Notes:
410 Users provide callback functions and their contexts to `TS` using, for example, `TSSetIFunction()`. These values are stored
411 in a `DMTS` that is contained in the `DM` associated with the `TS`. If no `DM` was provided by
412 the user with `TSSetDM()` it is automatically created by `TSGetDM()` with `DMShellCreate()`.
413
414 Users very rarely need to worked directly with the `DMTS` object, rather they work with the `TS` and the `DM` they created
415
416 Multiple `DM` can share a single `DMTS`, often each `DM` is associated with
417 a grid refinement level. `DMGetDMTS()` returns the `DMTS` associated with a `DM`. `DMGetDMTSWrite()` returns a unique
418 `DMTS` that is only associated with the current `DM`, making a copy of the shared `DMTS` if needed (copy-on-write).
419
420 See `DMKSP` for details on why there is a needed for `DMTS` instead of simply storing the user callbacks directly in the `DM` or the `TS`
421
422 Developer Note:
423 The original `dm` inside the `DMTS` is NOT reference counted (to prevent a reference count loop between a `DM` and a `DMSNES`).
424 The `DM` on which this context was first created is cached here to implement one-way
425 copy-on-write. When `DMGetDMTSWrite()` sees a request using a different `DM`, it makes a copy of the `DMTS`.
426
427 .seealso: [](ch_ts), `TSCreate()`, `DM`, `DMGetDMTSWrite()`, `DMGetDMTS()`, `TSSetIFunction()`, `DMTSSetRHSFunctionContextDestroy()`,
428 `DMTSSetRHSJacobian()`, `DMTSGetRHSJacobian()`, `DMTSSetRHSJacobianContextDestroy()`, `DMTSSetIFunction()`, `DMTSGetIFunction()`,
429 `DMTSSetIFunctionContextDestroy()`, `DMTSSetIJacobian()`, `DMTSGetIJacobian()`, `DMTSSetIJacobianContextDestroy()`,
430 `DMTSSetI2Function()`, `DMTSGetI2Function()`, `DMTSSetI2FunctionContextDestroy()`, `DMTSSetI2Jacobian()`,
431 `DMTSGetI2Jacobian()`, `DMTSSetI2JacobianContextDestroy()`, `DMKSP`, `DMSNES`
432 S*/
433 struct _p_DMTS {
434 PETSCHEADER(struct _DMTSOps);
435 PetscContainer rhsfunctionctxcontainer;
436 PetscContainer rhsjacobianctxcontainer;
437
438 PetscContainer ifunctionctxcontainer;
439 PetscContainer ijacobianctxcontainer;
440
441 PetscContainer i2functionctxcontainer;
442 PetscContainer i2jacobianctxcontainer;
443
444 void *transientvarctx;
445
446 void *solutionctx;
447 void *forcingctx;
448
449 void *data;
450
451 /* See the developer note for DMTS above */
452 DM originaldm;
453 };
454
455 PETSC_INTERN PetscErrorCode DMTSUnsetRHSFunctionContext_Internal(DM);
456 PETSC_INTERN PetscErrorCode DMTSUnsetRHSJacobianContext_Internal(DM);
457 PETSC_INTERN PetscErrorCode DMTSUnsetIFunctionContext_Internal(DM);
458 PETSC_INTERN PetscErrorCode DMTSUnsetIJacobianContext_Internal(DM);
459 PETSC_INTERN PetscErrorCode DMTSUnsetI2FunctionContext_Internal(DM);
460 PETSC_INTERN PetscErrorCode DMTSUnsetI2JacobianContext_Internal(DM);
461
462 PETSC_EXTERN PetscErrorCode DMGetDMTS(DM, DMTS *);
463 PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM, DMTS *);
464 PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM, DM);
465 PETSC_EXTERN PetscErrorCode DMTSView(DMTS, PetscViewer);
466 PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS, PetscViewer);
467 PETSC_EXTERN PetscErrorCode DMTSCopy(DMTS, DMTS);
468
469 struct _n_TSEvent {
470 PetscReal *fvalue_prev; /* value of indicator function at the left end-point of the event interval */
471 PetscReal *fvalue; /* value of indicator function at the current point */
472 PetscReal *fvalue_right; /* value of indicator function at the right end-point of the event interval */
473 PetscInt *fsign_prev; /* sign of indicator function at the left end-point of the event interval */
474 PetscInt *fsign; /* sign of indicator function at the current point */
475 PetscInt *fsign_right; /* sign of indicator function at the right end-point of the event interval */
476 PetscReal ptime_prev; /* time at the previous point (left end-point of the event interval) */
477 PetscReal ptime_right; /* time at the right end-point of the event interval */
478 PetscReal ptime_cache; /* the point visited by the TS before the event interval was detected; cached - to reuse if necessary */
479 PetscReal timestep_cache; /* time step considered by the TS before the event interval was detected; cached - to reuse if necessary */
480 PetscInt *side; /* upon bracket subdivision, indicates which sub-bracket is taken further, -1 -> left one, +1 -> right one, +2 -> neither, 0 -> zero-crossing located */
481 PetscInt *side_prev; /* counts the repeating previous side's (with values: -n <=> '-1'*n; +n <=> '+1'*n); used in the Anderson-Bjorck iteration */
482 PetscReal timestep_postevent; /* first time step after the event; can be PETSC_DECIDE */
483 PetscReal timestep_2nd_postevent; /* second time step after the event; can be PETSC_DECIDE */
484 PetscReal timestep_min; /* minimum time step */
485 PetscBool *justrefined_AB; /* this flag shows if the given indicator function i = [0..nevents) participated in Anderson-Bjorck process in the last iteration of TSEventHandler() */
486 PetscReal *gamma_AB; /* cumulative scaling factor for the Anderson-Bjorck iteration */
487 PetscErrorCode (*indicator)(TS, PetscReal, Vec, PetscReal *, void *); /* this callback defines the user function(s) whose sign changes indicate events */
488 PetscErrorCode (*postevent)(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *); /* user post-event callback */
489 void *ctx; /* user context for indicator and postevent callbacks */
490 PetscInt *direction; /* zero crossing direction to trigger the event: +1 -> going positive, -1 -> going negative, 0 -> any */
491 PetscBool *terminate; /* 1 -> terminate time stepping on event location, 0 -> continue */
492 PetscInt nevents; /* number of events (indicator functions) to handle on the current MPI process */
493 PetscInt nevents_zero; /* number of events triggered */
494 PetscInt *events_zero; /* list of the events triggered */
495 PetscReal *vtol; /* array of tolerances for the indicator function zero check */
496 PetscInt iterctr; /* iteration counter: used both for reporting and as a status indicator */
497 PetscBool processing; /* this flag shows if the event-resolving iterations are in progress, or the post-event dt handling is in progress */
498 PetscBool revisit_right; /* [sync] "revisit the bracket's right end", if true, then fvalue(s) are not calculated, but are taken from fvalue_right(s) */
499 PetscViewer monitor;
500 /* Struct to record the events */
501 struct {
502 PetscInt ctr; /* Recorder counter */
503 PetscReal *time; /* Event times */
504 PetscInt *stepnum; /* Step numbers */
505 PetscInt *nevents; /* Number of events occurring at the event times */
506 PetscInt **eventidx; /* Local indices of the events in the event list */
507 } recorder;
508 PetscInt recsize; /* Size of recorder stack */
509 PetscInt refct; /* Reference count */
510 };
511
512 PETSC_EXTERN PetscErrorCode TSEventInitialize(TSEvent, TS, PetscReal, Vec);
513 PETSC_EXTERN PetscErrorCode TSEventDestroy(TSEvent *);
514 PETSC_EXTERN PetscErrorCode TSEventHandler(TS);
515 PETSC_EXTERN PetscErrorCode TSAdjointEventHandler(TS);
516
517 PETSC_EXTERN PetscLogEvent TS_AdjointStep;
518 PETSC_EXTERN PetscLogEvent TS_Step;
519 PETSC_EXTERN PetscLogEvent TS_PseudoComputeTimeStep;
520 PETSC_EXTERN PetscLogEvent TS_FunctionEval;
521 PETSC_EXTERN PetscLogEvent TS_JacobianEval;
522 PETSC_EXTERN PetscLogEvent TS_ForwardStep;
523
524 typedef enum {
525 TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
526 TS_STEP_PENDING, /* vec_sol advanced, but step has not been accepted yet */
527 TS_STEP_COMPLETE /* step accepted and ptime, steps, etc have been advanced */
528 } TSStepStatus;
529
530 struct _n_TSMonitorLGCtx {
531 PetscDrawLG lg;
532 PetscBool semilogy;
533 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
534 PetscInt ksp_its, snes_its;
535 char **names;
536 char **displaynames;
537 PetscInt ndisplayvariables;
538 PetscInt *displayvariables;
539 PetscReal *displayvalues;
540 PetscErrorCode (*transform)(void *, Vec, Vec *);
541 PetscCtxDestroyFn *transformdestroy;
542 void *transformctx;
543 };
544
545 struct _n_TSMonitorSPCtx {
546 PetscDrawSP sp;
547 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
548 PetscInt retain; /* Retain n points plotted to show trajectories, or -1 for all points */
549 PetscBool phase; /* Plot in phase space rather than coordinate space */
550 PetscBool multispecies; /* Change scatter point color based on species */
551 PetscInt ksp_its, snes_its;
552 };
553
554 struct _n_TSMonitorHGCtx {
555 PetscDrawHG *hg;
556 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
557 PetscInt Ns; /* The number of species to histogram */
558 PetscBool velocity; /* Plot in velocity space rather than coordinate space */
559 };
560
561 struct _n_TSMonitorEnvelopeCtx {
562 Vec max, min;
563 };
564
565 /*
566 Checks if the user provide a TSSetIFunction() but an explicit method is called; generate an error in that case
567 */
TSCheckImplicitTerm(TS ts)568 static inline PetscErrorCode TSCheckImplicitTerm(TS ts)
569 {
570 TSIFunctionFn *ifunction;
571 DM dm;
572
573 PetscFunctionBegin;
574 PetscCall(TSGetDM(ts, &dm));
575 PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
576 PetscCheck(!ifunction, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_INCOMP, "You are attempting to use an explicit ODE integrator but provided an implicit function definition with TSSetIFunction()");
577 PetscFunctionReturn(PETSC_SUCCESS);
578 }
579
580 PETSC_INTERN PetscErrorCode TSGetRHSMats_Private(TS, Mat *, Mat *);
581 /* this is declared here as TSHistory is not public */
582 PETSC_EXTERN PetscErrorCode TSAdaptHistorySetTSHistory(TSAdapt, TSHistory, PetscBool);
583
584 PETSC_INTERN PetscErrorCode TSTrajectoryReconstruct_Private(TSTrajectory, TS, PetscReal, Vec, Vec);
585 PETSC_INTERN PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory, TS);
586
587 PETSC_EXTERN PetscLogEvent TSTrajectory_Set;
588 PETSC_EXTERN PetscLogEvent TSTrajectory_Get;
589 PETSC_EXTERN PetscLogEvent TSTrajectory_GetVecs;
590 PETSC_EXTERN PetscLogEvent TSTrajectory_SetUp;
591 PETSC_EXTERN PetscLogEvent TSTrajectory_DiskWrite;
592 PETSC_EXTERN PetscLogEvent TSTrajectory_DiskRead;
593
594 struct _n_TSMonitorDrawCtx {
595 PetscViewer viewer;
596 Vec initialsolution;
597 PetscBool showinitial;
598 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
599 PetscBool showtimestepandtime;
600 };
601
602 struct _n_TSMonitorVTKCtx {
603 char *filenametemplate;
604 PetscInt interval; /* when > 0 uses step % interval, when negative only final solution plotted */
605 };
606
607 struct _n_TSMonitorSolutionCtx {
608 PetscBool skip_initial; // Skip the viewer the first time TSMonitorSolution is run (within a single call to `TSSolve()`)
609 };
610