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