Actual source code: tsimpl.h

  1: #pragma once

  3: #include <petscts.h>
  4: #include <petsc/private/petscimpl.h>

  6: /*
  7:     Timesteping context.
  8:       General DAE: F(t,U,U_t) = 0, required Jacobian is G'(U) where G(U) = F(t,U,U0+a*U)
  9:       General ODE: U_t = F(t,U) <-- the right-hand-side function
 10:       Linear  ODE: U_t = A(t) U <-- the right-hand-side matrix
 11:       Linear (no time) ODE: U_t = A U <-- the right-hand-side matrix
 12: */

 14: /*
 15:      Maximum number of monitors you can run with a single TS
 16: */
 17: #define MAXTSMONITORS 10

 19: PETSC_EXTERN PetscBool      TSRegisterAllCalled;
 20: PETSC_EXTERN PetscErrorCode TSRegisterAll(void);
 21: PETSC_EXTERN PetscErrorCode TSAdaptRegisterAll(void);

 23: PETSC_EXTERN PetscErrorCode TSRKRegisterAll(void);
 24: PETSC_EXTERN PetscErrorCode TSMPRKRegisterAll(void);
 25: PETSC_EXTERN PetscErrorCode TSARKIMEXRegisterAll(void);
 26: PETSC_EXTERN PetscErrorCode TSRosWRegisterAll(void);
 27: PETSC_EXTERN PetscErrorCode TSGLLERegisterAll(void);
 28: PETSC_EXTERN PetscErrorCode TSGLLEAdaptRegisterAll(void);
 29: PETSC_EXTERN PetscErrorCode TSIRKRegisterAll(void);

 31: typedef struct _TSOps *TSOps;

 33: struct _TSOps {
 34:   PetscErrorCode (*snesfunction)(SNES, Vec, Vec, TS);
 35:   PetscErrorCode (*snesjacobian)(SNES, Vec, Mat, Mat, TS);
 36:   PetscErrorCode (*setup)(TS);
 37:   PetscErrorCode (*step)(TS);
 38:   PetscErrorCode (*solve)(TS);
 39:   PetscErrorCode (*interpolate)(TS, PetscReal, Vec);
 40:   PetscErrorCode (*evaluatewlte)(TS, NormType, PetscInt *, PetscReal *);
 41:   PetscErrorCode (*evaluatestep)(TS, PetscInt, Vec, PetscBool *);
 42:   PetscErrorCode (*setfromoptions)(TS, PetscOptionItems *);
 43:   PetscErrorCode (*destroy)(TS);
 44:   PetscErrorCode (*view)(TS, PetscViewer);
 45:   PetscErrorCode (*reset)(TS);
 46:   PetscErrorCode (*linearstability)(TS, PetscReal, PetscReal, PetscReal *, PetscReal *);
 47:   PetscErrorCode (*load)(TS, PetscViewer);
 48:   PetscErrorCode (*rollback)(TS);
 49:   PetscErrorCode (*getstages)(TS, PetscInt *, Vec *[]);
 50:   PetscErrorCode (*adjointstep)(TS);
 51:   PetscErrorCode (*adjointsetup)(TS);
 52:   PetscErrorCode (*adjointreset)(TS);
 53:   PetscErrorCode (*adjointintegral)(TS);
 54:   PetscErrorCode (*forwardsetup)(TS);
 55:   PetscErrorCode (*forwardreset)(TS);
 56:   PetscErrorCode (*forwardstep)(TS);
 57:   PetscErrorCode (*forwardintegral)(TS);
 58:   PetscErrorCode (*forwardgetstages)(TS, PetscInt *, Mat *[]);
 59:   PetscErrorCode (*getsolutioncomponents)(TS, PetscInt *, Vec *);
 60:   PetscErrorCode (*getauxsolution)(TS, Vec *);
 61:   PetscErrorCode (*gettimeerror)(TS, PetscInt, Vec *);
 62:   PetscErrorCode (*settimeerror)(TS, Vec);
 63:   PetscErrorCode (*startingmethod)(TS);
 64:   PetscErrorCode (*initcondition)(TS, Vec);
 65:   PetscErrorCode (*exacterror)(TS, Vec, Vec);
 66:   PetscErrorCode (*resizeregister)(TS, PetscBool);
 67: };

 69: /*
 70:    TSEvent - Abstract object to handle event monitoring
 71: */
 72: typedef struct _n_TSEvent *TSEvent;

 74: typedef struct _TSTrajectoryOps *TSTrajectoryOps;

 76: struct _TSTrajectoryOps {
 77:   PetscErrorCode (*view)(TSTrajectory, PetscViewer);
 78:   PetscErrorCode (*reset)(TSTrajectory);
 79:   PetscErrorCode (*destroy)(TSTrajectory);
 80:   PetscErrorCode (*set)(TSTrajectory, TS, PetscInt, PetscReal, Vec);
 81:   PetscErrorCode (*get)(TSTrajectory, TS, PetscInt, PetscReal *);
 82:   PetscErrorCode (*setfromoptions)(TSTrajectory, PetscOptionItems *);
 83:   PetscErrorCode (*setup)(TSTrajectory, TS);
 84: };

 86: /* TSHistory is an helper object that allows inquiring
 87:    the TSTrajectory by time and not by the step number only */
 88: typedef struct _n_TSHistory *TSHistory;

 90: struct _p_TSTrajectory {
 91:   PETSCHEADER(struct _TSTrajectoryOps);
 92:   TSHistory tsh; /* associates times to unique step ids */
 93:   /* stores necessary data to reconstruct states and derivatives via Lagrangian interpolation */
 94:   struct {
 95:     PetscInt     order; /* interpolation order */
 96:     Vec         *W;     /* work vectors */
 97:     PetscScalar *L;     /* workspace for Lagrange basis */
 98:     PetscReal   *T;     /* Lagrange times (stored) */
 99:     Vec         *WW;    /* just an array of pointers */
100:     PetscBool   *TT;    /* workspace for Lagrange */
101:     PetscReal   *TW;    /* Lagrange times (workspace) */

103:     /* caching */
104:     PetscBool caching;
105:     struct {
106:       PetscObjectId    id;
107:       PetscObjectState state;
108:       PetscReal        time;
109:       PetscInt         step;
110:     } Ucached;
111:     struct {
112:       PetscObjectId    id;
113:       PetscObjectState state;
114:       PetscReal        time;
115:       PetscInt         step;
116:     } Udotcached;
117:   } lag;
118:   Vec         U, Udot;            /* used by TSTrajectory{Get|Restore}UpdatedHistoryVecs */
119:   PetscBool   usehistory;         /* whether to use TSHistory */
120:   PetscBool   solution_only;      /* whether we dump just the solution or also the stages */
121:   PetscBool   adjoint_solve_mode; /* whether we will use the Trajectory inside a TSAdjointSolve() or not */
122:   PetscViewer monitor;
123:   PetscInt    setupcalled;            /* true if setup has been called */
124:   PetscInt    recomps;                /* counter for recomputations in the adjoint run */
125:   PetscInt    diskreads, diskwrites;  /* counters for disk checkpoint reads and writes */
126:   char      **names;                  /* the name of each variable; each process has only the local names */
127:   PetscBool   keepfiles;              /* keep the files generated during the run after the run is complete */
128:   char       *dirname, *filetemplate; /* directory name and file name template for disk checkpoints */
129:   char       *dirfiletemplate;        /* complete directory and file name template for disk checkpoints */
130:   PetscErrorCode (*transform)(void *, Vec, Vec *);
131:   PetscErrorCode (*transformdestroy)(void *);
132:   void *transformctx;
133:   void *data;
134: };

136: typedef struct _TS_RHSSplitLink *TS_RHSSplitLink;
137: struct _TS_RHSSplitLink {
138:   TS              ts;
139:   char           *splitname;
140:   IS              is;
141:   TS_RHSSplitLink next;
142:   PetscLogEvent   event;
143: };

145: typedef struct _TS_TimeSpan *TSTimeSpan;
146: struct _TS_TimeSpan {
147:   PetscInt   num_span_times; /* number of time points */
148:   PetscReal *span_times;     /* array of the time span */
149:   PetscReal  reltol;         /* relative tolerance for span point detection */
150:   PetscReal  abstol;         /* absolute tolerance for span point detection */
151:   PetscInt   spanctr;        /* counter of the time points that have been reached */
152:   Vec       *vecs_sol;       /* array of the solutions at the specified time points */
153: };

155: struct _p_TS {
156:   PETSCHEADER(struct _TSOps);
157:   TSProblemType  problem_type;
158:   TSEquationType equation_type;

160:   DM          dm;
161:   Vec         vec_sol; /* solution vector in first and second order equations */
162:   Vec         vec_dot; /* time derivative vector in second order equations */
163:   TSAdapt     adapt;
164:   TSAdaptType default_adapt_type;
165:   TSEvent     event;

167:   /* ---------------- Resize ---------------------*/
168:   PetscObjectList resizetransferobjs;

170:   /* ---------------- User (or PETSc) Provided stuff ---------------------*/
171:   PetscErrorCode (*monitor[MAXTSMONITORS])(TS, PetscInt, PetscReal, Vec, void *);
172:   PetscErrorCode (*monitordestroy[MAXTSMONITORS])(void **);
173:   void    *monitorcontext[MAXTSMONITORS];
174:   PetscInt numbermonitors;
175:   PetscErrorCode (*adjointmonitor[MAXTSMONITORS])(TS, PetscInt, PetscReal, Vec, PetscInt, Vec *, Vec *, void *);
176:   PetscErrorCode (*adjointmonitordestroy[MAXTSMONITORS])(void **);
177:   void    *adjointmonitorcontext[MAXTSMONITORS];
178:   PetscInt numberadjointmonitors;
179:   PetscInt monitorFrequency; /* Number of timesteps between monitor output */

181:   PetscErrorCode (*prestep)(TS);
182:   PetscErrorCode (*prestage)(TS, PetscReal);
183:   PetscErrorCode (*poststage)(TS, PetscReal, PetscInt, Vec *);
184:   PetscErrorCode (*postevaluate)(TS);
185:   PetscErrorCode (*poststep)(TS);
186:   PetscErrorCode (*functiondomainerror)(TS, PetscReal, Vec, PetscBool *);
187:   PetscErrorCode (*resizesetup)(TS, PetscInt, PetscReal, Vec, PetscBool *, void *);
188:   PetscErrorCode (*resizetransfer)(TS, PetscInt, Vec[], Vec[], void *);
189:   void *resizectx;

191:   /* ---------------------- Sensitivity Analysis support -----------------*/
192:   TSTrajectory trajectory; /* All solutions are kept here for the entire time integration process */
193:   Vec         *vecs_sensi; /* one vector for each cost function */
194:   Vec         *vecs_sensip;
195:   PetscInt     numcost; /* number of cost functions */
196:   Vec          vec_costintegral;
197:   PetscInt     adjointsetupcalled;
198:   PetscInt     adjoint_steps;
199:   PetscInt     adjoint_max_steps;
200:   PetscBool    adjoint_solve;     /* immediately call TSAdjointSolve() after TSSolve() is complete */
201:   PetscBool    costintegralfwd;   /* cost integral is evaluated in the forward run if true */
202:   Vec          vec_costintegrand; /* workspace for Adjoint computations */
203:   Mat          Jacp, Jacprhs;
204:   void        *ijacobianpctx, *rhsjacobianpctx;
205:   void        *costintegrandctx;
206:   Vec         *vecs_drdu;
207:   Vec         *vecs_drdp;
208:   Vec          vec_drdu_col, vec_drdp_col;

210:   /* first-order adjoint */
211:   PetscErrorCode (*rhsjacobianp)(TS, PetscReal, Vec, Mat, void *);
212:   PetscErrorCode (*ijacobianp)(TS, PetscReal, Vec, Vec, PetscReal, Mat, void *);
213:   PetscErrorCode (*costintegrand)(TS, PetscReal, Vec, Vec, void *);
214:   PetscErrorCode (*drdufunction)(TS, PetscReal, Vec, Vec *, void *);
215:   PetscErrorCode (*drdpfunction)(TS, PetscReal, Vec, Vec *, void *);

217:   /* second-order adjoint */
218:   Vec  *vecs_sensi2;
219:   Vec  *vecs_sensi2p;
220:   Vec   vec_dir; /* directional vector for optimization */
221:   Vec  *vecs_fuu, *vecs_guu;
222:   Vec  *vecs_fup, *vecs_gup;
223:   Vec  *vecs_fpu, *vecs_gpu;
224:   Vec  *vecs_fpp, *vecs_gpp;
225:   void *ihessianproductctx, *rhshessianproductctx;
226:   PetscErrorCode (*ihessianproduct_fuu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
227:   PetscErrorCode (*ihessianproduct_fup)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
228:   PetscErrorCode (*ihessianproduct_fpu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
229:   PetscErrorCode (*ihessianproduct_fpp)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
230:   PetscErrorCode (*rhshessianproduct_guu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
231:   PetscErrorCode (*rhshessianproduct_gup)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
232:   PetscErrorCode (*rhshessianproduct_gpu)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);
233:   PetscErrorCode (*rhshessianproduct_gpp)(TS, PetscReal, Vec, Vec *, Vec, Vec *, void *);

235:   /* specific to forward sensitivity analysis */
236:   Mat       mat_sensip;           /* matrix storing forward sensitivities */
237:   Vec       vec_sensip_col;       /* space for a column of the sensip matrix */
238:   Vec      *vecs_integral_sensip; /* one vector for each integral */
239:   PetscInt  num_parameters;
240:   PetscInt  num_initialvalues;
241:   void     *vecsrhsjacobianpctx;
242:   PetscInt  forwardsetupcalled;
243:   PetscBool forward_solve;
244:   PetscErrorCode (*vecsrhsjacobianp)(TS, PetscReal, Vec, Vec *, void *);

246:   /* ---------------------- IMEX support ---------------------------------*/
247:   /* These extra slots are only used when the user provides both Implicit and RHS */
248:   Mat Arhs; /* Right hand side matrix */
249:   Mat Brhs; /* Right hand side preconditioning matrix */
250:   Vec Frhs; /* Right hand side function value */

252:   /* This is a general caching scheme to avoid recomputing the Jacobian at a place that has been previously been evaluated.
253:    * The present use case is that TSComputeRHSFunctionLinear() evaluates the Jacobian once and we don't want it to be immeditely re-evaluated.
254:    */
255:   struct {
256:     PetscReal        time;       /* The time at which the matrices were last evaluated */
257:     PetscObjectId    Xid;        /* Unique ID of solution vector at which the Jacobian was last evaluated */
258:     PetscObjectState Xstate;     /* State of the solution vector */
259:     MatStructure     mstructure; /* The structure returned */
260:     /* Flag to unshift Jacobian before calling the IJacobian or RHSJacobian functions.  This is useful
261:      * if the user would like to reuse (part of) the Jacobian from the last evaluation. */
262:     PetscBool reuse;
263:     PetscReal scale, shift;
264:   } rhsjacobian;

266:   struct {
267:     PetscReal shift; /* The derivative of the lhs wrt to Xdot */
268:   } ijacobian;

270:   MatStructure axpy_pattern; /* information about the nonzero pattern of the RHS Jacobian in reference to the implicit Jacobian */
271:   /* --------------------Nonlinear Iteration------------------------------*/
272:   SNES      snes;
273:   PetscBool usessnes; /* Flag set by each TSType to indicate if the type actually uses a SNES;
274:                            this works around the design flaw that a SNES is ALWAYS created with TS even when it is not needed.*/
275:   PetscInt  ksp_its;  /* total number of linear solver iterations */
276:   PetscInt  snes_its; /* total number of nonlinear solver iterations */
277:   PetscInt  num_snes_failures;
278:   PetscInt  max_snes_failures;

280:   /* --- Logging --- */
281:   PetscInt ifuncs, rhsfuncs, ijacs, rhsjacs;

283:   /* --- Data that is unique to each particular solver --- */
284:   PetscInt setupcalled; /* true if setup has been called */
285:   void    *data;        /* implementationspecific data */
286:   void    *user;        /* user context */

288:   /* ------------------  Parameters -------------------------------------- */
289:   PetscInt  max_steps; /* max number of steps */
290:   PetscReal max_time;  /* max time allowed */

292:   /* --------------------------------------------------------------------- */

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:   PetscInt  steps;               /* steps taken so far in all successive calls to TSSolve() */
297:   PetscReal ptime;               /* time at the start of the current step (stage time is internal if it exists) */
298:   PetscReal time_step;           /* current time increment */
299:   PetscReal ptime_prev;          /* time at the start of the previous step */
300:   PetscReal ptime_prev_rollback; /* time at the start of the 2nd previous step to recover from rollback */
301:   PetscReal solvetime;           /* time at the conclusion of TSSolve() */
302:   PetscBool stifflyaccurate;     /* flag to indicate that the method is stiffly accurate */

304:   TSConvergedReason      reason;
305:   PetscBool              errorifstepfailed;
306:   PetscInt               reject, max_reject;
307:   TSExactFinalTimeOption exact_final_time;

309:   PetscReal atol, rtol;   /* Relative and absolute tolerance for local truncation error */
310:   Vec       vatol, vrtol; /* Relative and absolute tolerance in vector form */
311:   PetscReal cfltime, cfltime_local;

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;

324:   /* ---------------------- Quadrature integration support ---------------------------------*/
325:   TS quadraturets;

327:   /* ---------------------- Time span support ---------------------------------*/
328:   TSTimeSpan tspan;
329: };

331: struct _TSAdaptOps {
332:   PetscErrorCode (*choose)(TSAdapt, TS, PetscReal, PetscInt *, PetscReal *, PetscBool *, PetscReal *, PetscReal *, PetscReal *);
333:   PetscErrorCode (*destroy)(TSAdapt);
334:   PetscErrorCode (*reset)(TSAdapt);
335:   PetscErrorCode (*view)(TSAdapt, PetscViewer);
336:   PetscErrorCode (*setfromoptions)(TSAdapt, PetscOptionItems *);
337:   PetscErrorCode (*load)(TSAdapt, PetscViewer);
338: };

340: struct _p_TSAdapt {
341:   PETSCHEADER(struct _TSAdaptOps);
342:   void *data;
343:   PetscErrorCode (*checkstage)(TSAdapt, TS, PetscReal, Vec, PetscBool *);
344:   struct {
345:     PetscInt    n;              /* number of candidate schemes, including the one currently in use */
346:     PetscBool   inuse_set;      /* the current scheme has been set */
347:     const char *name[16];       /* name of the scheme */
348:     PetscInt    order[16];      /* classical order of each scheme */
349:     PetscInt    stageorder[16]; /* stage order of each scheme */
350:     PetscReal   ccfl[16];       /* stability limit relative to explicit Euler */
351:     PetscReal   cost[16];       /* relative measure of the amount of work required for each scheme */
352:   } candidates;
353:   PetscBool   always_accept;
354:   PetscReal   safety;             /* safety factor relative to target error/stability goal */
355:   PetscReal   reject_safety;      /* extra safety factor if the last step was rejected */
356:   PetscReal   clip[2];            /* admissible time step decrease/increase factors */
357:   PetscReal   dt_min, dt_max;     /* admissible minimum and maximum time step */
358:   PetscReal   ignore_max;         /* minimum value of the solution to be considered by the adaptor */
359:   PetscBool   glee_use_local;     /* GLEE adaptor uses global or local error */
360:   PetscReal   scale_solve_failed; /* scale step by this factor if solver (linear or nonlinear) fails. */
361:   PetscReal   matchstepfac[2];    /* factors to control the behaviour of matchstep */
362:   NormType    wnormtype;
363:   PetscViewer monitor;
364:   PetscInt    timestepjustdecreased_delay; /* number of timesteps after a decrease in the timestep before the timestep can be increased */
365:   PetscInt    timestepjustdecreased;
366:   PetscReal   dt_span_cached; /* time step before hitting a TS span time point */
367: };

369: typedef struct _p_DMTS  *DMTS;
370: typedef struct _DMTSOps *DMTSOps;
371: struct _DMTSOps {
372:   TSRHSFunction rhsfunction;
373:   TSRHSJacobian rhsjacobian;

375:   TSIFunction ifunction;
376:   PetscErrorCode (*ifunctionview)(void *, PetscViewer);
377:   PetscErrorCode (*ifunctionload)(void **, PetscViewer);

379:   TSIJacobian ijacobian;
380:   PetscErrorCode (*ijacobianview)(void *, PetscViewer);
381:   PetscErrorCode (*ijacobianload)(void **, PetscViewer);

383:   TSI2Function i2function;
384:   TSI2Jacobian i2jacobian;

386:   TSTransientVariable transientvar;

388:   TSSolutionFunction solution;
389:   TSForcingFunction  forcing;

391:   PetscErrorCode (*destroy)(DMTS);
392:   PetscErrorCode (*duplicate)(DMTS, DMTS);
393: };

395: struct _p_DMTS {
396:   PETSCHEADER(struct _DMTSOps);
397:   PetscContainer rhsfunctionctxcontainer;
398:   PetscContainer rhsjacobianctxcontainer;

400:   PetscContainer ifunctionctxcontainer;
401:   PetscContainer ijacobianctxcontainer;

403:   PetscContainer i2functionctxcontainer;
404:   PetscContainer i2jacobianctxcontainer;

406:   void *transientvarctx;

408:   void *solutionctx;
409:   void *forcingctx;

411:   void *data;

413:   /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
414:    * copy-on-write. When DMGetDMTSWrite() sees a request using a different DM, it makes a copy. Thus, if a user
415:    * only interacts directly with one level, e.g., using TSSetIFunction(), then coarse levels of a multilevel item
416:    * integrator are built, then the user changes the routine with another call to TSSetIFunction(), it automatically
417:    * propagates to all the levels. If instead, they get out a specific level and set the function on that level,
418:    * subsequent changes to the original level will no longer propagate to that level.
419:    */
420:   DM originaldm;
421: };

423: PETSC_EXTERN PetscErrorCode DMTSUnsetRHSFunctionContext_Internal(DM);
424: PETSC_EXTERN PetscErrorCode DMTSUnsetRHSJacobianContext_Internal(DM);
425: PETSC_EXTERN PetscErrorCode DMTSUnsetIFunctionContext_Internal(DM);
426: PETSC_EXTERN PetscErrorCode DMTSUnsetIJacobianContext_Internal(DM);
427: PETSC_EXTERN PetscErrorCode DMTSUnsetI2FunctionContext_Internal(DM);
428: PETSC_EXTERN PetscErrorCode DMTSUnsetI2JacobianContext_Internal(DM);

430: PETSC_EXTERN PetscErrorCode DMGetDMTS(DM, DMTS *);
431: PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM, DMTS *);
432: PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM, DM);
433: PETSC_EXTERN PetscErrorCode DMTSView(DMTS, PetscViewer);
434: PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS, PetscViewer);
435: PETSC_EXTERN PetscErrorCode DMTSCopy(DMTS, DMTS);

437: typedef enum {
438:   TSEVENT_NONE,
439:   TSEVENT_LOCATED_INTERVAL,
440:   TSEVENT_PROCESSING,
441:   TSEVENT_ZERO,
442:   TSEVENT_RESET_NEXTSTEP
443: } TSEventStatus;

445: struct _n_TSEvent {
446:   PetscScalar *fvalue;                                                                      /* value of event function at the end of the step*/
447:   PetscScalar *fvalue_prev;                                                                 /* value of event function at start of the step (left end-point of event interval) */
448:   PetscReal    ptime_prev;                                                                  /* time at step start (left end-point of event interval) */
449:   PetscReal    ptime_end;                                                                   /* end time of step (when an event interval is detected, ptime_end is fixed to the time at step end during event processing) */
450:   PetscReal    ptime_right;                                                                 /* time on the right end-point of the event interval */
451:   PetscScalar *fvalue_right;                                                                /* value of event function at the right end-point of the event interval */
452:   PetscInt    *side;                                                                        /* Used for detecting repetition of end-point, -1 => left, +1 => right */
453:   PetscReal    timestep_prev;                                                               /* previous time step */
454:   PetscReal    timestep_posteventinterval;                                                  /* time step immediately after the event interval */
455:   PetscReal    timestep_postevent;                                                          /* time step immediately after the event */
456:   PetscReal    timestep_min;                                                                /* Minimum time step */
457:   PetscBool   *zerocrossing;                                                                /* Flag to signal zero crossing detection */
458:   PetscErrorCode (*eventhandler)(TS, PetscReal, Vec, PetscScalar *, void *);                /* User event handler function */
459:   PetscErrorCode (*postevent)(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *); /* User post event function */
460:   void         *ctx;                                                                        /* User context for event handler and post even functions */
461:   PetscInt     *direction;                                                                  /* Zero crossing direction: 1 -> Going positive, -1 -> Going negative, 0 -> Any */
462:   PetscBool    *terminate;                                                                  /* 1 -> Terminate time stepping, 0 -> continue */
463:   PetscInt      nevents;                                                                    /* Number of events to handle */
464:   PetscInt      nevents_zero;                                                               /* Number of event zero detected */
465:   PetscInt     *events_zero;                                                                /* List of events that have reached zero */
466:   PetscReal    *vtol;                                                                       /* Vector tolerances for event zero check */
467:   TSEventStatus status;                                                                     /* Event status */
468:   PetscInt      iterctr;                                                                    /* Iteration counter */
469:   PetscViewer   monitor;
470:   /* Struct to record the events */
471:   struct {
472:     PetscInt   ctr;      /* recorder counter */
473:     PetscReal *time;     /* Event times */
474:     PetscInt  *stepnum;  /* Step numbers */
475:     PetscInt  *nevents;  /* Number of events occurring at the event times */
476:     PetscInt **eventidx; /* Local indices of the events in the event list */
477:   } recorder;
478:   PetscInt recsize; /* Size of recorder stack */
479:   PetscInt refct;   /* reference count */
480: };

482: PETSC_EXTERN PetscErrorCode TSEventInitialize(TSEvent, TS, PetscReal, Vec);
483: PETSC_EXTERN PetscErrorCode TSEventDestroy(TSEvent *);
484: PETSC_EXTERN PetscErrorCode TSEventHandler(TS);
485: PETSC_EXTERN PetscErrorCode TSAdjointEventHandler(TS);

487: PETSC_EXTERN PetscLogEvent TS_AdjointStep;
488: PETSC_EXTERN PetscLogEvent TS_Step;
489: PETSC_EXTERN PetscLogEvent TS_PseudoComputeTimeStep;
490: PETSC_EXTERN PetscLogEvent TS_FunctionEval;
491: PETSC_EXTERN PetscLogEvent TS_JacobianEval;
492: PETSC_EXTERN PetscLogEvent TS_ForwardStep;

494: typedef enum {
495:   TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
496:   TS_STEP_PENDING,    /* vec_sol advanced, but step has not been accepted yet */
497:   TS_STEP_COMPLETE    /* step accepted and ptime, steps, etc have been advanced */
498: } TSStepStatus;

500: struct _n_TSMonitorLGCtx {
501:   PetscDrawLG lg;
502:   PetscBool   semilogy;
503:   PetscInt    howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
504:   PetscInt    ksp_its, snes_its;
505:   char      **names;
506:   char      **displaynames;
507:   PetscInt    ndisplayvariables;
508:   PetscInt   *displayvariables;
509:   PetscReal  *displayvalues;
510:   PetscErrorCode (*transform)(void *, Vec, Vec *);
511:   PetscErrorCode (*transformdestroy)(void *);
512:   void *transformctx;
513: };

515: struct _n_TSMonitorSPCtx {
516:   PetscDrawSP sp;
517:   PetscInt    howoften;     /* when > 0 uses step % howoften, when negative only final solution plotted */
518:   PetscInt    retain;       /* Retain n points plotted to show trajectories, or -1 for all points */
519:   PetscBool   phase;        /* Plot in phase space rather than coordinate space */
520:   PetscBool   multispecies; /* Change scatter point color based on species */
521:   PetscInt    ksp_its, snes_its;
522: };

524: struct _n_TSMonitorHGCtx {
525:   PetscDrawHG *hg;
526:   PetscInt     howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
527:   PetscInt     Ns;       /* The number of species to histogram */
528:   PetscBool    velocity; /* Plot in velocity space rather than coordinate space */
529: };

531: struct _n_TSMonitorEnvelopeCtx {
532:   Vec max, min;
533: };

535: /*
536:     Checks if the user provide a TSSetIFunction() but an explicit method is called; generate an error in that case
537: */
538: static inline PetscErrorCode TSCheckImplicitTerm(TS ts)
539: {
540:   TSIFunction ifunction;
541:   DM          dm;

543:   PetscFunctionBegin;
544:   PetscCall(TSGetDM(ts, &dm));
545:   PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
546:   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()");
547:   PetscFunctionReturn(PETSC_SUCCESS);
548: }

550: PETSC_EXTERN PetscErrorCode TSGetRHSMats_Private(TS, Mat *, Mat *);
551: /* this is declared here as TSHistory is not public */
552: PETSC_EXTERN PetscErrorCode TSAdaptHistorySetTSHistory(TSAdapt, TSHistory, PetscBool);

554: PETSC_INTERN PetscErrorCode TSTrajectoryReconstruct_Private(TSTrajectory, TS, PetscReal, Vec, Vec);
555: PETSC_INTERN PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory, TS);

557: PETSC_EXTERN PetscLogEvent TSTrajectory_Set;
558: PETSC_EXTERN PetscLogEvent TSTrajectory_Get;
559: PETSC_EXTERN PetscLogEvent TSTrajectory_GetVecs;
560: PETSC_EXTERN PetscLogEvent TSTrajectory_SetUp;
561: PETSC_EXTERN PetscLogEvent TSTrajectory_DiskWrite;
562: PETSC_EXTERN PetscLogEvent TSTrajectory_DiskRead;

564: struct _n_TSMonitorDrawCtx {
565:   PetscViewer viewer;
566:   Vec         initialsolution;
567:   PetscBool   showinitial;
568:   PetscInt    howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
569:   PetscBool   showtimestepandtime;
570: };