Actual source code: tsimpl.h
petsc-3.6.2 2015-10-02
Report Typos and Errors

      2: #ifndef __TSIMPL_H

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

      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);

     24: typedef struct _TSOps *TSOps;

     26: struct _TSOps {
     27:   PetscErrorCode (*snesfunction)(SNES,Vec,Vec,TS);
     28:   PetscErrorCode (*snesjacobian)(SNES,Vec,Mat,Mat,TS);
     29:   PetscErrorCode (*setup)(TS);
     30:   PetscErrorCode (*step)(TS);
     31:   PetscErrorCode (*solve)(TS);
     32:   PetscErrorCode (*interpolate)(TS,PetscReal,Vec);
     33:   PetscErrorCode (*evaluatestep)(TS,PetscInt,Vec,PetscBool*);
     34:   PetscErrorCode (*setfromoptions)(PetscOptions*,TS);
     35:   PetscErrorCode (*destroy)(TS);
     36:   PetscErrorCode (*view)(TS,PetscViewer);
     37:   PetscErrorCode (*reset)(TS);
     38:   PetscErrorCode (*linearstability)(TS,PetscReal,PetscReal,PetscReal*,PetscReal*);
     39:   PetscErrorCode (*load)(TS,PetscViewer);
     40:   PetscErrorCode (*rollback)(TS);
     41:   PetscErrorCode (*getstages)(TS,PetscInt*,Vec**);
     42:   PetscErrorCode (*adjointstep)(TS);
     43:   PetscErrorCode (*adjointsetup)(TS);
     44: };

     46: /*
     47:    TSEvent - Abstract object to handle event monitoring
     48: */
     49: typedef struct _p_TSEvent *TSEvent;

     51: typedef struct _TSTrajectoryOps *TSTrajectoryOps;

     53: struct _TSTrajectoryOps {
     54:   PetscErrorCode (*view)(TSTrajectory,PetscViewer);
     55:   PetscErrorCode (*destroy)(TSTrajectory);
     56:   PetscErrorCode (*set)(TSTrajectory,TS,PetscInt,PetscReal,Vec);
     57:   PetscErrorCode (*get)(TSTrajectory,TS,PetscInt,PetscReal*);
     58: };

     60: struct _p_TSTrajectory {
     61:   PETSCHEADER(struct _TSTrajectoryOps);
     62:   void *data;
     63: };

     65: struct _p_TS {
     66:   PETSCHEADER(struct _TSOps);
     67:   DM            dm;
     68:   TSProblemType problem_type;
     69:   Vec           vec_sol;
     70:   TSAdapt       adapt;
     71:   TSEvent       event;

     73:   /* ---------------- User (or PETSc) Provided stuff ---------------------*/
     74:   PetscErrorCode (*monitor[MAXTSMONITORS])(TS,PetscInt,PetscReal,Vec,void*); /* returns control to user after */
     75:   PetscErrorCode (*monitordestroy[MAXTSMONITORS])(void**);
     76:   void *monitorcontext[MAXTSMONITORS];                 /* residual calculation, allows user */
     77:   PetscInt  numbermonitors;                                 /* to, for instance, print residual norm, etc. */

     79:   PetscErrorCode (*prestep)(TS);
     80:   PetscErrorCode (*prestage)(TS,PetscReal);
     81:   PetscErrorCode (*poststage)(TS,PetscReal,PetscInt,Vec*);
     82:   PetscErrorCode (*poststep)(TS);

     84:   /* ---------------------- Sensitivity Analysis support -----------------*/
     85:   TSTrajectory trajectory;   /* All solutions are kept here for the entire time integration process */
     86:   Vec       *vecs_sensi;             /* one vector for each cost function */
     87:   Vec       *vecs_sensip;
     88:   PetscInt  numcost;                 /* number of cost functions */
     89:   Vec       vec_costintegral;
     90:   PetscInt  adjointsetupcalled;
     91:   PetscInt  adjoint_max_steps;
     92:   PetscBool adjoint_solve;          /* immediately call TSAdjointSolve() after TSSolve() is complete */
     93:   PetscBool costintegralfwd;        /* cost integral is evaluated in the forward run if true */
     94:   /* workspace for Adjoint computations */
     95:   Vec       vec_costintegrand;
     96:   Mat       Jacp;
     97:   void      *rhsjacobianpctx;
     98:   void      *costintegrandctx;
     99:   Vec       *vecs_drdy;
    100:   Vec       *vecs_drdp;

    102:   PetscErrorCode (*rhsjacobianp)(TS,PetscReal,Vec,Mat,void*);
    103:   PetscErrorCode (*costintegrand)(TS,PetscReal,Vec,Vec,void*);
    104:   PetscErrorCode (*drdyfunction)(TS,PetscReal,Vec,Vec*,void*);
    105:   PetscErrorCode (*drdpfunction)(TS,PetscReal,Vec,Vec*,void*);

    107:   /* ---------------------- IMEX support ---------------------------------*/
    108:   /* These extra slots are only used when the user provides both Implicit and RHS */
    109:   Mat Arhs;     /* Right hand side matrix */
    110:   Mat Brhs;     /* Right hand side preconditioning matrix */
    111:   Vec Frhs;     /* Right hand side function value */

    113:   /* This is a general caching scheme to avoid recomputing the Jacobian at a place that has been previously been evaluated.
    114:    * The present use case is that TSComputeRHSFunctionLinear() evaluates the Jacobian once and we don't want it to be immeditely re-evaluated.
    115:    */
    116:   struct {
    117:     PetscReal time;             /* The time at which the matrices were last evaluated */
    118:     Vec X;                      /* Solution vector at which the Jacobian was last evaluated */
    119:     PetscObjectState Xstate;    /* State of the solution vector */
    120:     MatStructure mstructure;    /* The structure returned */
    121:     /* Flag to unshift Jacobian before calling the IJacobian or RHSJacobian functions.  This is useful
    122:      * if the user would like to reuse (part of) the Jacobian from the last evaluation. */
    123:     PetscBool reuse;
    124:     PetscReal scale,shift;
    125:   } rhsjacobian;

    127:   struct {
    128:     PetscReal shift;            /* The derivative of the lhs wrt to Xdot */
    129:   } ijacobian;

    131:   /* ---------------------Nonlinear Iteration------------------------------*/
    132:   SNES  snes;

    134:   /* --- Data that is unique to each particular solver --- */
    135:   PetscInt setupcalled;             /* true if setup has been called */
    136:   void     *data;                   /* implementationspecific data */
    137:   void     *user;                   /* user context */

    139:   /* ------------------  Parameters -------------------------------------- */
    140:   PetscInt  max_steps;              /* max number of steps */
    141:   PetscReal max_time;               /* max time allowed */
    142:   PetscReal time_step;              /* current/completed time increment */
    143:   PetscReal time_step_prev;         /* previous time step  */

    145:   /*
    146:      these are temporary to support increasing the time step if nonlinear solver convergence remains good
    147:      and time_step was previously cut due to failed nonlinear solver
    148:   */
    149:   PetscReal time_step_orig;            /* original time step requested by user */
    150:   PetscInt  time_steps_since_decrease; /* number of timesteps since timestep was decreased due to lack of convergence */
    151:   /* ----------------------------------------------------------------------------------------------------------------*/

    153:   PetscBool steprollback;           /* Is the current step rolled back? */
    154:   PetscInt  steps;                  /* steps taken so far in latest call to TSSolve() */
    155:   PetscInt  total_steps;            /* steps taken in all calls to TSSolve() since the TS was created or since TSSetUp() was called */
    156:   PetscReal ptime;                  /* time at the start of the current step (stage time is internal if it exists) */
    157:   PetscReal ptime_prev;             /* time at the start of the previous step */
    158:   PetscReal solvetime;              /* time at the conclusion of TSSolve() */
    159:   PetscInt  ksp_its;                /* total number of linear solver iterations */
    160:   PetscInt  snes_its;               /* total number of nonlinear solver iterations */

    162:   PetscInt num_snes_failures;
    163:   PetscInt max_snes_failures;
    164:   TSConvergedReason reason;
    165:   TSEquationType equation_type;
    166:   PetscBool errorifstepfailed;
    167:   TSExactFinalTimeOption  exact_final_time;
    168:   PetscBool retain_stages;
    169:   PetscInt reject,max_reject;

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

    175:   /* ------------------- Default work-area management ------------------ */
    176:   PetscInt nwork;
    177:   Vec      *work;
    178: };

    180: struct _TSAdaptOps {
    181:   PetscErrorCode (*choose)(TSAdapt,TS,PetscReal,PetscInt*,PetscReal*,PetscBool*,PetscReal*);
    182:   PetscErrorCode (*destroy)(TSAdapt);
    183:   PetscErrorCode (*reset)(TSAdapt);
    184:   PetscErrorCode (*view)(TSAdapt,PetscViewer);
    185:   PetscErrorCode (*setfromoptions)(PetscOptions*,TSAdapt);
    186:   PetscErrorCode (*load)(TSAdapt,PetscViewer);
    187: };

    189: struct _p_TSAdapt {
    190:   PETSCHEADER(struct _TSAdaptOps);
    191:   void *data;
    192:   PetscErrorCode (*checkstage)(TSAdapt,TS,PetscBool*);
    193:   struct {
    194:     PetscInt   n;                /* number of candidate schemes, including the one currently in use */
    195:     PetscBool  inuse_set;        /* the current scheme has been set */
    196:     const char *name[16];        /* name of the scheme */
    197:     PetscInt   order[16];        /* classical order of each scheme */
    198:     PetscInt   stageorder[16];   /* stage order of each scheme */
    199:     PetscReal  ccfl[16];         /* stability limit relative to explicit Euler */
    200:     PetscReal  cost[16];         /* relative measure of the amount of work required for each scheme */
    201:   } candidates;
    202:   PetscReal   dt_min,dt_max;
    203:   PetscReal   scale_solve_failed; /* Scale step by this factor if solver (linear or nonlinear) fails. */
    204:   PetscViewer monitor;
    205:   NormType    wnormtype;
    206: };

    208: typedef struct _p_DMTS *DMTS;
    209: typedef struct _DMTSOps *DMTSOps;
    210: struct _DMTSOps {
    211:   TSRHSFunction rhsfunction;
    212:   TSRHSJacobian rhsjacobian;

    214:   TSIFunction ifunction;
    215:   PetscErrorCode (*ifunctionview)(void*,PetscViewer);
    216:   PetscErrorCode (*ifunctionload)(void**,PetscViewer);

    218:   TSIJacobian ijacobian;
    219:   PetscErrorCode (*ijacobianview)(void*,PetscViewer);
    220:   PetscErrorCode (*ijacobianload)(void**,PetscViewer);

    222:   TSSolutionFunction solution;
    223:   PetscErrorCode (*forcing)(TS,PetscReal,Vec,void*);

    225:   PetscErrorCode (*destroy)(DMTS);
    226:   PetscErrorCode (*duplicate)(DMTS,DMTS);
    227: };

    229: struct _p_DMTS {
    230:   PETSCHEADER(struct _DMTSOps);
    231:   void *rhsfunctionctx;
    232:   void *rhsjacobianctx;

    234:   void *ifunctionctx;
    235:   void *ijacobianctx;

    237:   void *solutionctx;
    238:   void *forcingctx;

    240:   void *data;

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

    252: PETSC_EXTERN PetscErrorCode DMGetDMTS(DM,DMTS*);
    253: PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM,DMTS*);
    254: PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM,DM);
    255: PETSC_EXTERN PetscErrorCode DMTSView(DMTS,PetscViewer);
    256: PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS,PetscViewer);

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

    260: /* Maximum number of event times that can be recorded */
    261: #define MAXEVENTRECORDERS 24

    263: struct _p_TSEvent {
    264:   PetscScalar    *fvalue;          /* value of event function at the end of the step*/
    265:   PetscScalar    *fvalue_prev;     /* value of event function at start of the step */
    266:   PetscReal       ptime;           /* time at step end */
    267:   PetscReal       ptime_prev;      /* time at step start */
    268:   PetscErrorCode  (*monitor)(TS,PetscReal,Vec,PetscScalar*,void*); /* User event monitor function */
    269:   PetscErrorCode  (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,PetscBool,void*); /* User post event function */
    270:   PetscBool      *terminate;        /* 1 -> Terminate time stepping, 0 -> continue */
    271:   PetscInt       *direction;        /* Zero crossing direction: 1 -> Going positive, -1 -> Going negative, 0 -> Any */
    272:   PetscInt        nevents;          /* Number of events to handle */
    273:   PetscInt        nevents_zero;     /* Number of event zero detected */
    274:   PetscInt        *events_zero;      /* List of events that have reached zero */
    275:   void           *monitorcontext;
    276:   PetscReal      *vtol;             /* Vector tolerances for event zero check */
    277:   TSEventStatus   status;           /* Event status */
    278:   PetscReal       tstepend;         /* End time of step */
    279:   PetscReal       initial_timestep; /* Initial time step */
    280:   PetscViewer     mon;
    281:   /* Struct to record the events */
    282:   struct {
    283:     PetscInt  ctr;                          /* recorder counter */
    284:     PetscReal time[MAXEVENTRECORDERS];      /* Event times */
    285:     PetscInt  stepnum[MAXEVENTRECORDERS];   /* Step numbers */
    286:     PetscInt  nevents[MAXEVENTRECORDERS];   /* Number of events occuring at the event times */
    287:     PetscInt  *eventidx[MAXEVENTRECORDERS]; /* Local indices of the events in the event list */
    288:   } recorder;
    289: };

    291: PETSC_EXTERN PetscErrorCode TSEventMonitor(TS);
    292: PETSC_EXTERN PetscErrorCode TSEventMonitorDestroy(TSEvent*);
    293: PETSC_EXTERN PetscErrorCode TSAdjointEventMonitor(TS);
    294: PETSC_EXTERN PetscErrorCode TSEventMonitorInitialize(TS);

    296: PETSC_EXTERN PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval;

    298: typedef enum {TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
    299:               TS_STEP_PENDING,    /* vec_sol advanced, but step has not been accepted yet */
    300:               TS_STEP_COMPLETE    /* step accepted and ptime, steps, etc have been advanced */
    301: } TSStepStatus;

    303: struct _n_TSMonitorLGCtx {
    304:   PetscDrawLG    lg;
    305:   PetscInt       howoften;  /* when > 0 uses step % howoften, when negative only final solution plotted */
    306:   PetscInt       ksp_its,snes_its;
    307:   char           **names;
    308:   char           **displaynames;
    309:   PetscInt       ndisplayvariables;
    310:   PetscInt       *displayvariables;
    311:   PetscReal      *displayvalues;
    312:   PetscErrorCode (*transform)(void*,Vec,Vec*);
    313:   PetscErrorCode (*transformdestroy)(void*);
    314:   void           *transformctx;
    315: };

    317: struct _n_TSMonitorEnvelopeCtx {
    318:   Vec max,min;
    319: };

    322: #endif