Actual source code: petsc-tsimpl.h
petsc-3.4.2 2013-07-02

      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 5

     21: typedef struct _TSOps *TSOps;

     23: struct _TSOps {
     24:   PetscErrorCode (*snesfunction)(SNES,Vec,Vec,TS);
     25:   PetscErrorCode (*snesjacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,TS);
     26:   PetscErrorCode (*setup)(TS);
     27:   PetscErrorCode (*step)(TS);
     28:   PetscErrorCode (*solve)(TS);
     29:   PetscErrorCode (*interpolate)(TS,PetscReal,Vec);
     30:   PetscErrorCode (*evaluatestep)(TS,PetscInt,Vec,PetscBool*);
     31:   PetscErrorCode (*setfromoptions)(TS);
     32:   PetscErrorCode (*destroy)(TS);
     33:   PetscErrorCode (*view)(TS,PetscViewer);
     34:   PetscErrorCode (*reset)(TS);
     35:   PetscErrorCode (*linearstability)(TS,PetscReal,PetscReal,PetscReal*,PetscReal*);
     36:   PetscErrorCode (*load)(TS,PetscViewer);
     37: };

     39: struct _p_TS {
     40:   PETSCHEADER(struct _TSOps);
     41:   DM            dm;
     42:   TSProblemType problem_type;
     43:   Vec           vec_sol;
     44:   TSAdapt       adapt;

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

     52:   PetscErrorCode (*prestep)(TS);
     53:   PetscErrorCode (*prestage)(TS,PetscReal);
     54:   PetscErrorCode (*poststep)(TS);

     56:   /* ---------------------- IMEX support ---------------------------------*/
     57:   /* These extra slots are only used when the user provides both Implicit and RHS */
     58:   Mat Arhs;     /* Right hand side matrix */
     59:   Mat Brhs;     /* Right hand side preconditioning matrix */
     60:   Vec Frhs;     /* Right hand side function value */

     62:   /* This is a general caching scheme to avoid recomputing the Jacobian at a place that has been previously been evaluated.
     63:    * The present use case is that TSComputeRHSFunctionLinear() evaluates the Jacobian once and we don't want it to be immeditely re-evaluated.
     64:    */
     65:   struct {
     66:     PetscReal time;             /* The time at which the matrices were last evaluated */
     67:     Vec X;                      /* Solution vector at which the Jacobian was last evaluated */
     68:     PetscInt Xstate;            /* State of the solution vector */
     69:     MatStructure mstructure;    /* The structure returned */
     70:   } rhsjacobian;

     72:   struct {
     73:     PetscReal time;             /* The time at which the matrices were last evaluated */
     74:     Vec X;                      /* Solution vector at which the Jacobian was last evaluated */
     75:     Vec Xdot;                   /* Time derivative of the state vector at which the Jacobian was last evaluated */
     76:     PetscInt Xstate;            /* State of the solution vector */
     77:     PetscInt Xdotstate;         /* State of the solution vector */
     78:     MatStructure mstructure;    /* The structure returned */
     79:     PetscReal shift;            /* The derivative of the lhs wrt to Xdot */
     80:     PetscBool imex;             /* Flag of the method if it was started as an imex method */
     81:   } ijacobian;

     83:   /* ---------------------Nonlinear Iteration------------------------------*/
     84:   SNES  snes;

     86:   /* --- Data that is unique to each particular solver --- */
     87:   PetscInt setupcalled;             /* true if setup has been called */
     88:   void     *data;                   /* implementationspecific data */
     89:   void     *user;                   /* user context */

     91:   /* ------------------  Parameters -------------------------------------- */
     92:   PetscInt  max_steps;              /* max number of steps */
     93:   PetscReal max_time;               /* max time allowed */
     94:   PetscReal time_step;              /* current/completed time increment */
     95:   PetscReal time_step_prev;         /* previous time step  */

     97:   /*
     98:      these are temporary to support increasing the time step if nonlinear solver convergence remains good
     99:      and time_step was previously cut due to failed nonlinear solver
    100:   */
    101:   PetscReal time_step_orig;            /* original time step requested by user */
    102:   PetscInt  time_steps_since_decrease; /* number of timesteps since timestep was decreased due to lack of convergence */
    103:   /* ----------------------------------------------------------------------------------------------------------------*/

    105:   PetscInt  steps;                  /* steps taken so far */
    106:   PetscReal ptime;                  /* time at the start of the current step (stage time is internal if it exists) */
    107:   PetscReal solvetime;              /* time at the conclusion of TSSolve() */
    108:   PetscInt  ksp_its;                /* total number of linear solver iterations */
    109:   PetscInt  snes_its;               /* total number of nonlinear solver iterations */

    111:   PetscInt num_snes_failures;
    112:   PetscInt max_snes_failures;
    113:   TSConvergedReason reason;
    114:   TSEquationType equation_type;
    115:   PetscBool errorifstepfailed;
    116:   TSExactFinalTimeOption  exact_final_time;
    117:   PetscBool retain_stages;
    118:   PetscInt reject,max_reject;

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

    124:   /* ------------------- Default work-area management ------------------ */
    125:   PetscInt nwork;
    126:   Vec      *work;
    127: };

    129: struct _TSAdaptOps {
    130:   PetscErrorCode (*choose)(TSAdapt,TS,PetscReal,PetscInt*,PetscReal*,PetscBool*,PetscReal*);
    131:   PetscErrorCode (*checkstage)(TSAdapt,TS,PetscBool*);
    132:   PetscErrorCode (*destroy)(TSAdapt);
    133:   PetscErrorCode (*view)(TSAdapt,PetscViewer);
    134:   PetscErrorCode (*setfromoptions)(TSAdapt);
    135:   PetscErrorCode (*load)(TSAdapt,PetscViewer);
    136: };

    138: struct _p_TSAdapt {
    139:   PETSCHEADER(struct _TSAdaptOps);
    140:   void *data;
    141:   struct {
    142:     PetscInt   n;                /* number of candidate schemes, including the one currently in use */
    143:     PetscBool  inuse_set;        /* the current scheme has been set */
    144:     const char *name[16];        /* name of the scheme */
    145:     PetscInt   order[16];        /* classical order of each scheme */
    146:     PetscInt   stageorder[16];   /* stage order of each scheme */
    147:     PetscReal  ccfl[16];         /* stability limit relative to explicit Euler */
    148:     PetscReal  cost[16];         /* relative measure of the amount of work required for each scheme */
    149:   } candidates;
    150:   PetscReal   dt_min,dt_max;
    151:   PetscReal   scale_solve_failed; /* Scale step by this factor if solver (linear or nonlinear) fails. */
    152:   PetscViewer monitor;
    153: };

    155: typedef struct _p_DMTS *DMTS;
    156: typedef struct _DMTSOps *DMTSOps;
    157: struct _DMTSOps {
    158:   TSRHSFunction rhsfunction;
    159:   TSRHSJacobian rhsjacobian;

    161:   TSIFunction ifunction;
    162:   PetscErrorCode (*ifunctionview)(void*,PetscViewer);
    163:   PetscErrorCode (*ifunctionload)(void**,PetscViewer);

    165:   TSIJacobian ijacobian;
    166:   PetscErrorCode (*ijacobianview)(void*,PetscViewer);
    167:   PetscErrorCode (*ijacobianload)(void**,PetscViewer);

    169:   TSSolutionFunction solution;
    170:   PetscErrorCode (*forcing)(TS,PetscReal,Vec,void*);

    172:   PetscErrorCode (*destroy)(DMTS);
    173:   PetscErrorCode (*duplicate)(DMTS,DMTS);
    174: };

    176: struct _p_DMTS {
    177:   PETSCHEADER(struct _DMTSOps);
    178:   void *rhsfunctionctx;
    179:   void *rhsjacobianctx;

    181:   void *ifunctionctx;
    182:   void *ijacobianctx;

    184:   void *solutionctx;
    185:   void *forcingctx;

    187:   void *data;

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

    199: PETSC_EXTERN PetscErrorCode DMGetDMTS(DM,DMTS*);
    200: PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM,DMTS*);
    201: PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM,DM);
    202: PETSC_EXTERN PetscErrorCode DMTSView(DMTS,PetscViewer);
    203: PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS,PetscViewer);

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

    208: typedef enum {TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
    209:               TS_STEP_PENDING,    /* vec_sol advanced, but step has not been accepted yet */
    210:               TS_STEP_COMPLETE    /* step accepted and ptime, steps, etc have been advanced */
    211: } TSStepStatus;

    213: #endif