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: /*
     72:    TSEvent - Abstract object to handle event monitoring
     73: */
     74: typedef struct _n_TSEvent *TSEvent;

     76: typedef struct _TSTrajectoryOps *TSTrajectoryOps;

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

     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;

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

    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:   PetscInt    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:   PetscErrorCode (*transformdestroy)(void *);
    134:   void *transformctx;
    135:   void *data;
    136: };

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

    147: typedef struct _TS_TimeSpan *TSTimeSpan;
    148: struct _TS_TimeSpan {
    149:   PetscInt   num_span_times; /* number of time points */
    150:   PetscReal *span_times;     /* 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   spanctr;        /* counter of the time points that have been reached */
    155:   Vec       *vecs_sol;       /* array of the solutions at the specified time points */
    156: };

    158: struct _p_TS {
    159:   PETSCHEADER(struct _TSOps);
    160:   TSProblemType  problem_type;
    161:   TSEquationType equation_type;

    163:   DM          dm;
    164:   Vec         vec_sol;  /* solution vector in first and second order equations */
    165:   Vec         vec_sol0; /* solution vector at the beginning of the step */
    166:   Vec         vec_dot;  /* time derivative vector in second order equations */
    167:   TSAdapt     adapt;
    168:   TSAdaptType default_adapt_type;
    169:   TSEvent     event;

    171:   /* ---------------- Resize ---------------------*/
    172:   PetscBool       resizerollback;
    173:   PetscObjectList resizetransferobjs;

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

    186:   PetscErrorCode (*prestep)(TS);
    187:   PetscErrorCode (*prestage)(TS, PetscReal);
    188:   PetscErrorCode (*poststage)(TS, PetscReal, PetscInt, Vec *);
    189:   PetscErrorCode (*postevaluate)(TS);
    190:   PetscErrorCode (*poststep)(TS);
    191:   PetscErrorCode (*functiondomainerror)(TS, PetscReal, Vec, PetscBool *);
    192:   PetscErrorCode (*resizesetup)(TS, PetscInt, PetscReal, Vec, PetscBool *, void *);
    193:   PetscErrorCode (*resizetransfer)(TS, PetscInt, Vec[], Vec[], void *);
    194:   void *resizectx;

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

    215:   /* first-order adjoint */
    216:   PetscErrorCode (*rhsjacobianp)(TS, PetscReal, Vec, Mat, void *);
    217:   PetscErrorCode (*ijacobianp)(TS, PetscReal, Vec, Vec, PetscReal, Mat, void *);
    218:   PetscErrorCode (*costintegrand)(TS, PetscReal, Vec, Vec, void *);
    219:   PetscErrorCode (*drdufunction)(TS, PetscReal, Vec, Vec *, void *);
    220:   PetscErrorCode (*drdpfunction)(TS, PetscReal, Vec, Vec *, void *);

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

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

    251:   /* ---------------------- IMEX support ---------------------------------*/
    252:   /* These extra slots are only used when the user provides both Implicit and RHS */
    253:   Mat Arhs; /* Right hand side matrix */
    254:   Mat Brhs; /* Right hand side preconditioning matrix */
    255:   Vec Frhs; /* Right hand side function value */

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

    271:   struct {
    272:     PetscReal shift; /* The derivative of the lhs wrt to Xdot */
    273:   } ijacobian;

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

    285:   /* --- Logging --- */
    286:   PetscInt ifuncs, rhsfuncs, ijacs, rhsjacs;

    288:   /* --- Data that is unique to each particular solver --- */
    289:   PetscInt setupcalled; /* true if setup has been called */
    290:   void    *data;        /* implementationspecific data */
    291:   void    *user;        /* user context */

    293:   PetscBool steprollback;        /* flag to indicate that the step was rolled back */
    294:   PetscBool steprestart;         /* flag to indicate that the timestepper has to discard any history and restart */
    295:   PetscBool stepresize;          /* flag to indicate that the discretization was resized */
    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 time_step0;          /* proposed time increment at the beginning of the step */
    300:   PetscReal ptime_prev;          /* time at the start of the previous step */
    301:   PetscReal ptime_prev_rollback; /* time at the start of the 2nd previous step to recover from rollback */
    302:   PetscReal solvetime;           /* time at the conclusion of TSSolve() */
    303:   PetscBool stifflyaccurate;     /* flag to indicate that the method is stiffly accurate */

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

    310:   PetscObjectParameterDeclare(PetscReal, rtol); /* Relative and absolute tolerance for local truncation error */
    311:   PetscObjectParameterDeclare(PetscReal, atol);
    312:   PetscObjectParameterDeclare(PetscReal, max_time); /* max time allowed */
    313:   PetscObjectParameterDeclare(PetscInt, max_steps); /* max number of steps */
    314:   Vec       vatol, vrtol;                           /* Relative and absolute tolerance in vector form */
    315:   PetscReal cfltime, cfltime_local;

    317:   PetscBool testjacobian;
    318:   PetscBool testjacobiantranspose;
    319:   /* ------------------- Default work-area management ------------------ */
    320:   PetscInt nwork;
    321:   Vec     *work;

    323:   /* ---------------------- RHS splitting support ---------------------------------*/
    324:   PetscInt        num_rhs_splits;
    325:   TS_RHSSplitLink tsrhssplit;
    326:   PetscBool       use_splitrhsfunction;
    327:   SNES            snesrhssplit;

    329:   /* ---------------------- Quadrature integration support ---------------------------------*/
    330:   TS quadraturets;

    332:   /* ---------------------- Time span support ---------------------------------*/
    333:   TSTimeSpan tspan;
    334: };

    336: struct _TSAdaptOps {
    337:   PetscErrorCode (*choose)(TSAdapt, TS, PetscReal, PetscInt *, PetscReal *, PetscBool *, PetscReal *, PetscReal *, PetscReal *);
    338:   PetscErrorCode (*destroy)(TSAdapt);
    339:   PetscErrorCode (*reset)(TSAdapt);
    340:   PetscErrorCode (*view)(TSAdapt, PetscViewer);
    341:   PetscErrorCode (*setfromoptions)(TSAdapt, PetscOptionItems *);
    342:   PetscErrorCode (*load)(TSAdapt, PetscViewer);
    343: };

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

    374: typedef struct _p_DMTS  *DMTS;
    375: typedef struct _DMTSOps *DMTSOps;
    376: struct _DMTSOps {
    377:   TSRHSFunctionFn *rhsfunction;
    378:   TSRHSJacobianFn *rhsjacobian;

    380:   TSIFunctionFn *ifunction;
    381:   PetscErrorCode (*ifunctionview)(void *, PetscViewer);
    382:   PetscErrorCode (*ifunctionload)(void **, PetscViewer);

    384:   TSIJacobianFn *ijacobian;
    385:   PetscErrorCode (*ijacobianview)(void *, PetscViewer);
    386:   PetscErrorCode (*ijacobianload)(void **, PetscViewer);

    388:   TSI2FunctionFn *i2function;
    389:   TSI2JacobianFn *i2jacobian;

    391:   TSTransientVariableFn *transientvar;

    393:   TSSolutionFn *solution;
    394:   TSForcingFn  *forcing;

    396:   PetscErrorCode (*destroy)(DMTS);
    397:   PetscErrorCode (*duplicate)(DMTS, DMTS);
    398: };

    400: /*S
    401:    DMTS - Object held by a `DM` that contains all the callback functions and their contexts needed by a `TS`

    403:    Level: developer

    405:    Notes:
    406:    Users provide callback functions and their contexts to `TS` using, for example, `TSSetIFunction()`. These values are stored
    407:    in a `DMTS` that is contained in the `DM` associated with the `TS`. If no `DM` was provided by
    408:    the user with `TSSetDM()` it is automatically created by `TSGetDM()` with `DMShellCreate()`.

    410:    Users very rarely need to worked directly with the `DMTS` object, rather they work with the `TS` and the `DM` they created

    412:    Multiple `DM` can share a single `DMTS`, often each `DM` is associated with
    413:    a grid refinement level. `DMGetDMTS()` returns the `DMTS` associated with a `DM`. `DMGetDMTSWrite()` returns a unique
    414:    `DMTS` that is only associated with the current `DM`, making a copy of the shared `DMTS` if needed (copy-on-write).

    416:    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`

    418:    Developer Note:
    419:    The original `dm` inside the `DMTS` is NOT reference counted  (to prevent a reference count loop between a `DM` and a `DMSNES`).
    420:    The `DM` on which this context was first created is cached here to implement one-way
    421:    copy-on-write. When `DMGetDMTSWrite()` sees a request using a different `DM`, it makes a copy of the `DMTS`.

    423: .seealso: [](ch_ts), `TSCreate()`, `DM`, `DMGetDMTSWrite()`, `DMGetDMTS()`, `TSSetIFunction()`, `DMTSSetRHSFunctionContextDestroy()`,
    424:           `DMTSSetRHSJacobian()`, `DMTSGetRHSJacobian()`, `DMTSSetRHSJacobianContextDestroy()`, `DMTSSetIFunction()`, `DMTSGetIFunction()`,
    425:           `DMTSSetIFunctionContextDestroy()`, `DMTSSetIJacobian()`, `DMTSGetIJacobian()`, `DMTSSetIJacobianContextDestroy()`,
    426:           `DMTSSetI2Function()`, `DMTSGetI2Function()`, `DMTSSetI2FunctionContextDestroy()`, `DMTSSetI2Jacobian()`,
    427:           `DMTSGetI2Jacobian()`, `DMTSSetI2JacobianContextDestroy()`, `DMKSP`, `DMSNES`
    428: S*/
    429: struct _p_DMTS {
    430:   PETSCHEADER(struct _DMTSOps);
    431:   PetscContainer rhsfunctionctxcontainer;
    432:   PetscContainer rhsjacobianctxcontainer;

    434:   PetscContainer ifunctionctxcontainer;
    435:   PetscContainer ijacobianctxcontainer;

    437:   PetscContainer i2functionctxcontainer;
    438:   PetscContainer i2jacobianctxcontainer;

    440:   void *transientvarctx;

    442:   void *solutionctx;
    443:   void *forcingctx;

    445:   void *data;

    447:   /* See the developer note for DMTS above */
    448:   DM originaldm;
    449: };

    451: PETSC_EXTERN PetscErrorCode DMTSUnsetRHSFunctionContext_Internal(DM);
    452: PETSC_EXTERN PetscErrorCode DMTSUnsetRHSJacobianContext_Internal(DM);
    453: PETSC_EXTERN PetscErrorCode DMTSUnsetIFunctionContext_Internal(DM);
    454: PETSC_EXTERN PetscErrorCode DMTSUnsetIJacobianContext_Internal(DM);
    455: PETSC_EXTERN PetscErrorCode DMTSUnsetI2FunctionContext_Internal(DM);
    456: PETSC_EXTERN PetscErrorCode DMTSUnsetI2JacobianContext_Internal(DM);

    458: PETSC_EXTERN PetscErrorCode DMGetDMTS(DM, DMTS *);
    459: PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM, DMTS *);
    460: PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM, DM);
    461: PETSC_EXTERN PetscErrorCode DMTSView(DMTS, PetscViewer);
    462: PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS, PetscViewer);
    463: PETSC_EXTERN PetscErrorCode DMTSCopy(DMTS, DMTS);

    465: struct _n_TSEvent {
    466:   PetscReal *fvalue_prev;                                                                   /* value of indicator function at the left end-point of the event interval */
    467:   PetscReal *fvalue;                                                                        /* value of indicator function at the current point */
    468:   PetscReal *fvalue_right;                                                                  /* value of indicator function at the right end-point of the event interval */
    469:   PetscInt  *fsign_prev;                                                                    /* sign of indicator function at the left end-point of the event interval */
    470:   PetscInt  *fsign;                                                                         /* sign of indicator function at the current point */
    471:   PetscInt  *fsign_right;                                                                   /* sign of indicator function at the right end-point of the event interval */
    472:   PetscReal  ptime_prev;                                                                    /* time at the previous point (left end-point of the event interval) */
    473:   PetscReal  ptime_right;                                                                   /* time at the right end-point of the event interval */
    474:   PetscReal  ptime_cache;                                                                   /* the point visited by the TS before the event interval was detected; cached - to reuse if necessary */
    475:   PetscReal  timestep_cache;                                                                /* time step considered by the TS before the event interval was detected; cached - to reuse if necessary */
    476:   PetscInt  *side;                                                                          /* upon bracket subdivision, indicates which sub-bracket is taken further, -1 -> left one, +1 -> right one, +2 -> neither, 0 -> zero-crossing located */
    477:   PetscInt  *side_prev;                                                                     /* counts the repeating previous side's (with values: -n <=> '-1'*n; +n <=> '+1'*n); used in the Anderson-Bjorck iteration */
    478:   PetscReal  timestep_postevent;                                                            /* first time step after the event; can be PETSC_DECIDE */
    479:   PetscReal  timestep_2nd_postevent;                                                        /* second time step after the event; can be PETSC_DECIDE */
    480:   PetscReal  timestep_min;                                                                  /* minimum time step */
    481:   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() */
    482:   PetscReal *gamma_AB;                                                                      /* cumulative scaling factor for the Anderson-Bjorck iteration */
    483:   PetscErrorCode (*indicator)(TS, PetscReal, Vec, PetscReal *, void *);                     /* this callback defines the user function(s) whose sign changes indicate events */
    484:   PetscErrorCode (*postevent)(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *); /* user post-event callback */
    485:   void       *ctx;                                                                          /* user context for indicator and postevent callbacks */
    486:   PetscInt   *direction;                                                                    /* zero crossing direction to trigger the event: +1 -> going positive, -1 -> going negative, 0 -> any */
    487:   PetscBool  *terminate;                                                                    /* 1 -> terminate time stepping on event location, 0 -> continue */
    488:   PetscInt    nevents;                                                                      /* number of events (indicator functions) to handle on the current MPI process */
    489:   PetscInt    nevents_zero;                                                                 /* number of events triggered */
    490:   PetscInt   *events_zero;                                                                  /* list of the events triggered */
    491:   PetscReal  *vtol;                                                                         /* array of tolerances for the indicator function zero check */
    492:   PetscInt    iterctr;                                                                      /* iteration counter: used both for reporting and as a status indicator */
    493:   PetscBool   processing;                                                                   /* this flag shows if the event-resolving iterations are in progress, or the post-event dt handling is in progress */
    494:   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) */
    495:   PetscViewer monitor;
    496:   /* Struct to record the events */
    497:   struct {
    498:     PetscInt   ctr;      /* Recorder counter */
    499:     PetscReal *time;     /* Event times */
    500:     PetscInt  *stepnum;  /* Step numbers */
    501:     PetscInt  *nevents;  /* Number of events occurring at the event times */
    502:     PetscInt **eventidx; /* Local indices of the events in the event list */
    503:   } recorder;
    504:   PetscInt recsize; /* Size of recorder stack */
    505:   PetscInt refct;   /* Reference count */
    506: };

    508: PETSC_EXTERN PetscErrorCode TSEventInitialize(TSEvent, TS, PetscReal, Vec);
    509: PETSC_EXTERN PetscErrorCode TSEventDestroy(TSEvent *);
    510: PETSC_EXTERN PetscErrorCode TSEventHandler(TS);
    511: PETSC_EXTERN PetscErrorCode TSAdjointEventHandler(TS);

    513: PETSC_EXTERN PetscLogEvent TS_AdjointStep;
    514: PETSC_EXTERN PetscLogEvent TS_Step;
    515: PETSC_EXTERN PetscLogEvent TS_PseudoComputeTimeStep;
    516: PETSC_EXTERN PetscLogEvent TS_FunctionEval;
    517: PETSC_EXTERN PetscLogEvent TS_JacobianEval;
    518: PETSC_EXTERN PetscLogEvent TS_ForwardStep;

    520: typedef enum {
    521:   TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
    522:   TS_STEP_PENDING,    /* vec_sol advanced, but step has not been accepted yet */
    523:   TS_STEP_COMPLETE    /* step accepted and ptime, steps, etc have been advanced */
    524: } TSStepStatus;

    526: struct _n_TSMonitorLGCtx {
    527:   PetscDrawLG lg;
    528:   PetscBool   semilogy;
    529:   PetscInt    howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
    530:   PetscInt    ksp_its, snes_its;
    531:   char      **names;
    532:   char      **displaynames;
    533:   PetscInt    ndisplayvariables;
    534:   PetscInt   *displayvariables;
    535:   PetscReal  *displayvalues;
    536:   PetscErrorCode (*transform)(void *, Vec, Vec *);
    537:   PetscErrorCode (*transformdestroy)(void *);
    538:   void *transformctx;
    539: };

    541: struct _n_TSMonitorSPCtx {
    542:   PetscDrawSP sp;
    543:   PetscInt    howoften;     /* when > 0 uses step % howoften, when negative only final solution plotted */
    544:   PetscInt    retain;       /* Retain n points plotted to show trajectories, or -1 for all points */
    545:   PetscBool   phase;        /* Plot in phase space rather than coordinate space */
    546:   PetscBool   multispecies; /* Change scatter point color based on species */
    547:   PetscInt    ksp_its, snes_its;
    548: };

    550: struct _n_TSMonitorHGCtx {
    551:   PetscDrawHG *hg;
    552:   PetscInt     howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
    553:   PetscInt     Ns;       /* The number of species to histogram */
    554:   PetscBool    velocity; /* Plot in velocity space rather than coordinate space */
    555: };

    557: struct _n_TSMonitorEnvelopeCtx {
    558:   Vec max, min;
    559: };

    561: /*
    562:     Checks if the user provide a TSSetIFunction() but an explicit method is called; generate an error in that case
    563: */
    564: static inline PetscErrorCode TSCheckImplicitTerm(TS ts)
    565: {
    566:   TSIFunctionFn *ifunction;
    567:   DM             dm;

    569:   PetscFunctionBegin;
    570:   PetscCall(TSGetDM(ts, &dm));
    571:   PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
    572:   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()");
    573:   PetscFunctionReturn(PETSC_SUCCESS);
    574: }

    576: PETSC_EXTERN PetscErrorCode TSGetRHSMats_Private(TS, Mat *, Mat *);
    577: /* this is declared here as TSHistory is not public */
    578: PETSC_EXTERN PetscErrorCode TSAdaptHistorySetTSHistory(TSAdapt, TSHistory, PetscBool);

    580: PETSC_INTERN PetscErrorCode TSTrajectoryReconstruct_Private(TSTrajectory, TS, PetscReal, Vec, Vec);
    581: PETSC_INTERN PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory, TS);

    583: PETSC_EXTERN PetscLogEvent TSTrajectory_Set;
    584: PETSC_EXTERN PetscLogEvent TSTrajectory_Get;
    585: PETSC_EXTERN PetscLogEvent TSTrajectory_GetVecs;
    586: PETSC_EXTERN PetscLogEvent TSTrajectory_SetUp;
    587: PETSC_EXTERN PetscLogEvent TSTrajectory_DiskWrite;
    588: PETSC_EXTERN PetscLogEvent TSTrajectory_DiskRead;

    590: struct _n_TSMonitorDrawCtx {
    591:   PetscViewer viewer;
    592:   Vec         initialsolution;
    593:   PetscBool   showinitial;
    594:   PetscInt    howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
    595:   PetscBool   showtimestepandtime;
    596: };

    598: struct _n_TSMonitorVTKCtx {
    599:   char    *filenametemplate;
    600:   PetscInt interval; /* when > 0 uses step % interval, when negative only final solution plotted */
    601: };