Actual source code: tsimpl.h
petsc-3.7.7 2017-09-25
Report Typos and Errors

      1: #ifndef __TSIMPL_H

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

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

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

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

     24: PETSC_EXTERN PetscErrorCode TSRKRegisterAll(void);
     25: PETSC_EXTERN PetscErrorCode TSARKIMEXRegisterAll(void);
     26: PETSC_EXTERN PetscErrorCode TSRosWRegisterAll(void);
     27: PETSC_EXTERN PetscErrorCode TSGLRegisterAll(void);
     28: PETSC_EXTERN PetscErrorCode TSGLAdaptRegisterAll(void);

     30: typedef struct _TSOps *TSOps;

     32: struct _TSOps {
     33:   PetscErrorCode (*snesfunction)(SNES,Vec,Vec,TS);
     34:   PetscErrorCode (*snesjacobian)(SNES,Vec,Mat,Mat,TS);
     35:   PetscErrorCode (*setup)(TS);
     36:   PetscErrorCode (*step)(TS);
     37:   PetscErrorCode (*solve)(TS);
     38:   PetscErrorCode (*interpolate)(TS,PetscReal,Vec);
     39:   PetscErrorCode (*evaluatewlte)(TS,NormType,PetscInt*,PetscReal*);
     40:   PetscErrorCode (*evaluatestep)(TS,PetscInt,Vec,PetscBool*);
     41:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,TS);
     42:   PetscErrorCode (*destroy)(TS);
     43:   PetscErrorCode (*view)(TS,PetscViewer);
     44:   PetscErrorCode (*reset)(TS);
     45:   PetscErrorCode (*linearstability)(TS,PetscReal,PetscReal,PetscReal*,PetscReal*);
     46:   PetscErrorCode (*load)(TS,PetscViewer);
     47:   PetscErrorCode (*rollback)(TS);
     48:   PetscErrorCode (*getstages)(TS,PetscInt*,Vec**);
     49:   PetscErrorCode (*adjointstep)(TS);
     50:   PetscErrorCode (*adjointsetup)(TS);
     51:   PetscErrorCode (*adjointintegral)(TS);
     52:   PetscErrorCode (*forwardintegral)(TS);
     53: };

     55: /*
     56:    TSEvent - Abstract object to handle event monitoring
     57: */
     58: typedef struct _n_TSEvent *TSEvent;

     60: typedef struct _TSTrajectoryOps *TSTrajectoryOps;

     62: struct _TSTrajectoryOps {
     63:   PetscErrorCode (*view)(TSTrajectory,PetscViewer);
     64:   PetscErrorCode (*destroy)(TSTrajectory);
     65:   PetscErrorCode (*set)(TSTrajectory,TS,PetscInt,PetscReal,Vec);
     66:   PetscErrorCode (*get)(TSTrajectory,TS,PetscInt,PetscReal*);
     67:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,TSTrajectory);
     68:   PetscErrorCode (*setup)(TSTrajectory,TS);
     69: };

     71: struct _p_TSTrajectory {
     72:   PETSCHEADER(struct _TSTrajectoryOps);
     73:   PetscViewer monitor;
     74:   PetscInt    setupcalled;             /* true if setup has been called */
     75:   PetscInt    recomps;                 /* counter for recomputations in the adjoint run */
     76:   PetscInt    diskreads,diskwrites;    /* counters for disk checkpoint reads and writes */
     77:   void        *data;
     78: };

     80: struct _p_TS {
     81:   PETSCHEADER(struct _TSOps);
     82:   TSProblemType  problem_type;
     83:   TSEquationType equation_type;

     85:   DM             dm;
     86:   Vec            vec_sol; /* solution vector in first and second order equations */
     87:   Vec            vec_dot; /* time derivative vector in second order equations */
     88:   TSAdapt        adapt;
     89:   TSEvent        event;

     91:   /* ---------------- User (or PETSc) Provided stuff ---------------------*/
     92:   PetscErrorCode (*monitor[MAXTSMONITORS])(TS,PetscInt,PetscReal,Vec,void*);
     93:   PetscErrorCode (*monitordestroy[MAXTSMONITORS])(void**);
     94:   void *monitorcontext[MAXTSMONITORS];
     95:   PetscInt  numbermonitors;
     96:   PetscErrorCode (*adjointmonitor[MAXTSMONITORS])(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*);
     97:   PetscErrorCode (*adjointmonitordestroy[MAXTSMONITORS])(void**);
     98:   void *adjointmonitorcontext[MAXTSMONITORS];
     99:   PetscInt  numberadjointmonitors;

    101:   PetscErrorCode (*prestep)(TS);
    102:   PetscErrorCode (*prestage)(TS,PetscReal);
    103:   PetscErrorCode (*poststage)(TS,PetscReal,PetscInt,Vec*);
    104:   PetscErrorCode (*poststep)(TS);
    105:   PetscErrorCode (*functiondomainerror)(TS,PetscReal,Vec,PetscBool*);

    107:   /* ---------------------- Sensitivity Analysis support -----------------*/
    108:   TSTrajectory trajectory;          /* All solutions are kept here for the entire time integration process */
    109:   Vec       *vecs_sensi;            /* one vector for each cost function */
    110:   Vec       *vecs_sensip;
    111:   PetscInt  numcost;                /* number of cost functions */
    112:   Vec       vec_costintegral;
    113:   PetscInt  adjointsetupcalled;
    114:   PetscInt  adjoint_max_steps;
    115:   PetscBool adjoint_solve;          /* immediately call TSAdjointSolve() after TSSolve() is complete */
    116:   PetscBool costintegralfwd;        /* cost integral is evaluated in the forward run if true */
    117:   Vec       vec_costintegrand;      /* workspace for Adjoint computations */
    118:   Mat       Jacp;
    119:   void      *rhsjacobianpctx;
    120:   void      *costintegrandctx;
    121:   Vec       *vecs_drdy;
    122:   Vec       *vecs_drdp;

    124:   PetscErrorCode (*rhsjacobianp)(TS,PetscReal,Vec,Mat,void*);
    125:   PetscErrorCode (*costintegrand)(TS,PetscReal,Vec,Vec,void*);
    126:   PetscErrorCode (*drdyfunction)(TS,PetscReal,Vec,Vec*,void*);
    127:   PetscErrorCode (*drdpfunction)(TS,PetscReal,Vec,Vec*,void*);

    129:   /* ---------------------- IMEX support ---------------------------------*/
    130:   /* These extra slots are only used when the user provides both Implicit and RHS */
    131:   Mat Arhs;     /* Right hand side matrix */
    132:   Mat Brhs;     /* Right hand side preconditioning matrix */
    133:   Vec Frhs;     /* Right hand side function value */

    135:   /* This is a general caching scheme to avoid recomputing the Jacobian at a place that has been previously been evaluated.
    136:    * The present use case is that TSComputeRHSFunctionLinear() evaluates the Jacobian once and we don't want it to be immeditely re-evaluated.
    137:    */
    138:   struct {
    139:     PetscReal time;             /* The time at which the matrices were last evaluated */
    140:     Vec X;                      /* Solution vector at which the Jacobian was last evaluated */
    141:     PetscObjectState Xstate;    /* State of the solution vector */
    142:     MatStructure mstructure;    /* The structure returned */
    143:     /* Flag to unshift Jacobian before calling the IJacobian or RHSJacobian functions.  This is useful
    144:      * if the user would like to reuse (part of) the Jacobian from the last evaluation. */
    145:     PetscBool reuse;
    146:     PetscReal scale,shift;
    147:   } rhsjacobian;

    149:   struct {
    150:     PetscReal shift;            /* The derivative of the lhs wrt to Xdot */
    151:   } ijacobian;

    153:   /* --------------------Nonlinear Iteration------------------------------*/
    154:   SNES     snes;
    155:   PetscInt ksp_its;                /* total number of linear solver iterations */
    156:   PetscInt snes_its;               /* total number of nonlinear solver iterations */
    157:   PetscInt num_snes_failures;
    158:   PetscInt max_snes_failures;

    160:   /* --- Data that is unique to each particular solver --- */
    161:   PetscInt setupcalled;             /* true if setup has been called */
    162:   void     *data;                   /* implementationspecific data */
    163:   void     *user;                   /* user context */

    165:   /* ------------------  Parameters -------------------------------------- */
    166:   PetscInt  max_steps;              /* max number of steps */
    167:   PetscReal max_time;               /* max time allowed */

    169:   /* --------------------------------------------------------------------- */

    171:   PetscBool steprollback;           /* flag to indicate that the step was rolled back */
    172:   PetscBool steprestart;            /* flag to indicate that the timestepper has to discard any history and restart */
    173:   PetscInt  steps;                  /* steps taken so far in latest call to TSSolve() */
    174:   PetscInt  total_steps;            /* steps taken in all calls to TSSolve() since the TS was created or since TSSetUp() was called */
    175:   PetscReal ptime;                  /* time at the start of the current step (stage time is internal if it exists) */
    176:   PetscReal time_step;              /* current time increment */
    177:   PetscReal ptime_prev;             /* time at the start of the previous step */
    178:   PetscReal ptime_prev_rollback;    /* time at the start of the 2nd previous step to recover from rollback */
    179:   PetscReal solvetime;              /* time at the conclusion of TSSolve() */

    181:   TSConvergedReason reason;
    182:   PetscBool errorifstepfailed;
    183:   PetscInt  reject,max_reject;
    184:   TSExactFinalTimeOption exact_final_time;

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

    190:   /* ------------------- Default work-area management ------------------ */
    191:   PetscInt nwork;
    192:   Vec      *work;
    193: };

    195: struct _TSAdaptOps {
    196:   PetscErrorCode (*choose)(TSAdapt,TS,PetscReal,PetscInt*,PetscReal*,PetscBool*,PetscReal*);
    197:   PetscErrorCode (*destroy)(TSAdapt);
    198:   PetscErrorCode (*reset)(TSAdapt);
    199:   PetscErrorCode (*view)(TSAdapt,PetscViewer);
    200:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,TSAdapt);
    201:   PetscErrorCode (*load)(TSAdapt,PetscViewer);
    202: };

    204: struct _p_TSAdapt {
    205:   PETSCHEADER(struct _TSAdaptOps);
    206:   void *data;
    207:   PetscErrorCode (*checkstage)(TSAdapt,TS,PetscReal,Vec,PetscBool*);
    208:   struct {
    209:     PetscInt   n;                /* number of candidate schemes, including the one currently in use */
    210:     PetscBool  inuse_set;        /* the current scheme has been set */
    211:     const char *name[16];        /* name of the scheme */
    212:     PetscInt   order[16];        /* classical order of each scheme */
    213:     PetscInt   stageorder[16];   /* stage order of each scheme */
    214:     PetscReal  ccfl[16];         /* stability limit relative to explicit Euler */
    215:     PetscReal  cost[16];         /* relative measure of the amount of work required for each scheme */
    216:   } candidates;
    217:   PetscReal   dt_min,dt_max;
    218:   PetscReal   scale_solve_failed; /* Scale step by this factor if solver (linear or nonlinear) fails. */
    219:   PetscViewer monitor;
    220:   NormType    wnormtype;
    221: };

    223: typedef struct _p_DMTS *DMTS;
    224: typedef struct _DMTSOps *DMTSOps;
    225: struct _DMTSOps {
    226:   TSRHSFunction rhsfunction;
    227:   TSRHSJacobian rhsjacobian;

    229:   TSIFunction ifunction;
    230:   PetscErrorCode (*ifunctionview)(void*,PetscViewer);
    231:   PetscErrorCode (*ifunctionload)(void**,PetscViewer);

    233:   TSIJacobian ijacobian;
    234:   PetscErrorCode (*ijacobianview)(void*,PetscViewer);
    235:   PetscErrorCode (*ijacobianload)(void**,PetscViewer);

    237:   TSI2Function i2function;
    238:   TSI2Jacobian i2jacobian;

    240:   TSSolutionFunction solution;
    241:   TSForcingFunction  forcing;

    243:   PetscErrorCode (*destroy)(DMTS);
    244:   PetscErrorCode (*duplicate)(DMTS,DMTS);
    245: };

    247: struct _p_DMTS {
    248:   PETSCHEADER(struct _DMTSOps);
    249:   void *rhsfunctionctx;
    250:   void *rhsjacobianctx;

    252:   void *ifunctionctx;
    253:   void *ijacobianctx;

    255:   void *i2functionctx;
    256:   void *i2jacobianctx;

    258:   void *solutionctx;
    259:   void *forcingctx;

    261:   void *data;

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

    273: PETSC_EXTERN PetscErrorCode DMGetDMTS(DM,DMTS*);
    274: PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM,DMTS*);
    275: PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM,DM);
    276: PETSC_EXTERN PetscErrorCode DMTSView(DMTS,PetscViewer);
    277: PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS,PetscViewer);
    278: PETSC_EXTERN PetscErrorCode DMTSCopy(DMTS,DMTS);

    280: typedef enum {TSEVENT_NONE,TSEVENT_LOCATED_INTERVAL,TSEVENT_PROCESSING,TSEVENT_ZERO,TSEVENT_RESET_NEXTSTEP} TSEventStatus;

    282: struct _n_TSEvent {
    283:   PetscScalar    *fvalue;          /* value of event function at the end of the step*/
    284:   PetscScalar    *fvalue_prev;     /* value of event function at start of the step (left end-point of event interval) */
    285:   PetscReal       ptime_prev;      /* time at step start (left end-point of event interval) */
    286:   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) */
    287:   PetscReal       ptime_right;     /* time on the right end-point of the event interval */
    288:   PetscScalar    *fvalue_right;    /* value of event function at the right end-point of the event interval */
    289:   PetscInt       *side;            /* Used for detecting repetition of end-point, -1 => left, +1 => right */
    290:   PetscReal       timestep_prev;   /* previous time step */
    291:   PetscReal       timestep_orig;   /* initial time step */
    292:   PetscBool      *zerocrossing;    /* Flag to signal zero crossing detection */
    293:   PetscErrorCode  (*eventhandler)(TS,PetscReal,Vec,PetscScalar*,void*); /* User event handler function */
    294:   PetscErrorCode  (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,PetscBool,void*); /* User post event function */
    295:   void           *ctx;              /* User context for event handler and post even functions */
    296:   PetscInt       *direction;        /* Zero crossing direction: 1 -> Going positive, -1 -> Going negative, 0 -> Any */
    297:   PetscBool      *terminate;        /* 1 -> Terminate time stepping, 0 -> continue */
    298:   PetscInt        nevents;          /* Number of events to handle */
    299:   PetscInt        nevents_zero;     /* Number of event zero detected */
    300:   PetscInt       *events_zero;      /* List of events that have reached zero */
    301:   PetscReal      *vtol;             /* Vector tolerances for event zero check */
    302:   TSEventStatus   status;           /* Event status */
    303:   PetscInt        iterctr;          /* Iteration counter */
    304:   PetscViewer     monitor;
    305:   /* Struct to record the events */
    306:   struct {
    307:     PetscInt  ctr;        /* recorder counter */
    308:     PetscReal *time;      /* Event times */
    309:     PetscInt  *stepnum;   /* Step numbers */
    310:     PetscInt  *nevents;   /* Number of events occuring at the event times */
    311:     PetscInt  **eventidx; /* Local indices of the events in the event list */
    312:   } recorder;
    313:   PetscInt  recsize; /* Size of recorder stack */
    314: };

    316: PETSC_EXTERN PetscErrorCode TSEventInitialize(TSEvent,TS,PetscReal,Vec);
    317: PETSC_EXTERN PetscErrorCode TSEventDestroy(TSEvent*);
    318: PETSC_EXTERN PetscErrorCode TSEventHandler(TS);
    319: PETSC_EXTERN PetscErrorCode TSAdjointEventHandler(TS);

    321: PETSC_EXTERN PetscLogEvent TS_AdjointStep, TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;

    323: typedef enum {TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
    324:               TS_STEP_PENDING,    /* vec_sol advanced, but step has not been accepted yet */
    325:               TS_STEP_COMPLETE    /* step accepted and ptime, steps, etc have been advanced */
    326: } TSStepStatus;

    328: struct _n_TSMonitorLGCtx {
    329:   PetscDrawLG    lg;
    330:   PetscInt       howoften;  /* when > 0 uses step % howoften, when negative only final solution plotted */
    331:   PetscInt       ksp_its,snes_its;
    332:   char           **names;
    333:   char           **displaynames;
    334:   PetscInt       ndisplayvariables;
    335:   PetscInt       *displayvariables;
    336:   PetscReal      *displayvalues;
    337:   PetscErrorCode (*transform)(void*,Vec,Vec*);
    338:   PetscErrorCode (*transformdestroy)(void*);
    339:   void           *transformctx;
    340: };

    342: struct _n_TSMonitorEnvelopeCtx {
    343:   Vec max,min;
    344: };

    346: /*
    347:     Checks if the user provide a TSSetIFunction() but an explicit method is called; generate an error in that case
    348: */
    349: PETSC_STATIC_INLINE PetscErrorCode TSCheckImplicitTerm(TS ts)
    350: {
    351:   TSIFunction      ifunction;
    352:   DM               dm;
    353:   PetscErrorCode   ierr;

    356:   TSGetDM(ts,&dm);
    357:   DMTSGetIFunction(dm,&ifunction,NULL);
    358:   if (ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"You are attempting to use an explicit ODE integrator but provided an implicit function definition with TSSetIFunction()");
    359:   return(0);
    360: }

    362: PETSC_EXTERN PetscLogEvent TSTrajectory_Set, TSTrajectory_Get, TSTrajectory_DiskWrite, TSTrajectory_DiskRead;

    364: #endif