Actual source code: petsc-tsimpl.h

petsc-3.5.4 2015-05-23
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 5

 21: typedef struct _TSOps *TSOps;

 23: struct _TSOps {
 24:   PetscErrorCode (*snesfunction)(SNES,Vec,Vec,TS);
 25:   PetscErrorCode (*snesjacobian)(SNES,Vec,Mat,Mat,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:   PetscErrorCode (*rollback)(TS);
 38: };

 40: /* 
 41:    TSEvent - Abstract object to handle event monitoring
 42: */
 43: typedef struct _p_TSEvent *TSEvent;

 45: struct _p_TS {
 46:   PETSCHEADER(struct _TSOps);
 47:   DM            dm;
 48:   TSProblemType problem_type;
 49:   Vec           vec_sol;
 50:   TSAdapt       adapt;
 51:   TSEvent       event;

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

 59:   PetscErrorCode (*prestep)(TS);
 60:   PetscErrorCode (*prestage)(TS,PetscReal);
 61:   PetscErrorCode (*poststage)(TS,PetscReal,PetscInt,Vec*);
 62:   PetscErrorCode (*poststep)(TS);

 64:   /* ---------------------- IMEX support ---------------------------------*/
 65:   /* These extra slots are only used when the user provides both Implicit and RHS */
 66:   Mat Arhs;     /* Right hand side matrix */
 67:   Mat Brhs;     /* Right hand side preconditioning matrix */
 68:   Vec Frhs;     /* Right hand side function value */

 70:   /* This is a general caching scheme to avoid recomputing the Jacobian at a place that has been previously been evaluated.
 71:    * The present use case is that TSComputeRHSFunctionLinear() evaluates the Jacobian once and we don't want it to be immeditely re-evaluated.
 72:    */
 73:   struct {
 74:     PetscReal time;             /* The time at which the matrices were last evaluated */
 75:     Vec X;                      /* Solution vector at which the Jacobian was last evaluated */
 76:     PetscObjectState Xstate;    /* State of the solution vector */
 77:     MatStructure mstructure;    /* The structure returned */
 78:     /* Flag to unshift Jacobian before calling the IJacobian or RHSJacobian functions.  This is useful
 79:      * if the user would like to reuse (part of) the Jacobian from the last evaluation. */
 80:     PetscBool reuse;
 81:     PetscReal scale,shift;
 82:   } rhsjacobian;

 84:   struct {
 85:     PetscReal shift;            /* The derivative of the lhs wrt to Xdot */
 86:   } ijacobian;

 88:   /* ---------------------Nonlinear Iteration------------------------------*/
 89:   SNES  snes;

 91:   /* --- Data that is unique to each particular solver --- */
 92:   PetscInt setupcalled;             /* true if setup has been called */
 93:   void     *data;                   /* implementationspecific data */
 94:   void     *user;                   /* user context */

 96:   /* ------------------  Parameters -------------------------------------- */
 97:   PetscInt  max_steps;              /* max number of steps */
 98:   PetscReal max_time;               /* max time allowed */
 99:   PetscReal time_step;              /* current/completed time increment */
100:   PetscReal time_step_prev;         /* previous time step  */

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

110:   PetscInt  steps;                  /* steps taken so far */
111:   PetscReal ptime;                  /* time at the start of the current step (stage time is internal if it exists) */
112:   PetscReal ptime_prev;             /* time at the start of the previous step */
113:   PetscReal solvetime;              /* time at the conclusion of TSSolve() */
114:   PetscInt  ksp_its;                /* total number of linear solver iterations */
115:   PetscInt  snes_its;               /* total number of nonlinear solver iterations */

117:   PetscInt num_snes_failures;
118:   PetscInt max_snes_failures;
119:   TSConvergedReason reason;
120:   TSEquationType equation_type;
121:   PetscBool errorifstepfailed;
122:   TSExactFinalTimeOption  exact_final_time;
123:   PetscBool retain_stages;
124:   PetscInt reject,max_reject;

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

130:   /* ------------------- Default work-area management ------------------ */
131:   PetscInt nwork;
132:   Vec      *work;
133: };

135: struct _TSAdaptOps {
136:   PetscErrorCode (*choose)(TSAdapt,TS,PetscReal,PetscInt*,PetscReal*,PetscBool*,PetscReal*);
137:   PetscErrorCode (*checkstage)(TSAdapt,TS,PetscBool*);
138:   PetscErrorCode (*destroy)(TSAdapt);
139:   PetscErrorCode (*view)(TSAdapt,PetscViewer);
140:   PetscErrorCode (*setfromoptions)(TSAdapt);
141:   PetscErrorCode (*load)(TSAdapt,PetscViewer);
142: };

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

161: typedef struct _p_DMTS *DMTS;
162: typedef struct _DMTSOps *DMTSOps;
163: struct _DMTSOps {
164:   TSRHSFunction rhsfunction;
165:   TSRHSJacobian rhsjacobian;

167:   TSIFunction ifunction;
168:   PetscErrorCode (*ifunctionview)(void*,PetscViewer);
169:   PetscErrorCode (*ifunctionload)(void**,PetscViewer);

171:   TSIJacobian ijacobian;
172:   PetscErrorCode (*ijacobianview)(void*,PetscViewer);
173:   PetscErrorCode (*ijacobianload)(void**,PetscViewer);

175:   TSSolutionFunction solution;
176:   PetscErrorCode (*forcing)(TS,PetscReal,Vec,void*);

178:   PetscErrorCode (*destroy)(DMTS);
179:   PetscErrorCode (*duplicate)(DMTS,DMTS);
180: };

182: struct _p_DMTS {
183:   PETSCHEADER(struct _DMTSOps);
184:   void *rhsfunctionctx;
185:   void *rhsjacobianctx;

187:   void *ifunctionctx;
188:   void *ijacobianctx;

190:   void *solutionctx;
191:   void *forcingctx;

193:   void *data;

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

205: PETSC_EXTERN PetscErrorCode DMGetDMTS(DM,DMTS*);
206: PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM,DMTS*);
207: PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM,DM);
208: PETSC_EXTERN PetscErrorCode DMTSView(DMTS,PetscViewer);
209: PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS,PetscViewer);

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

213: struct _p_TSEvent {
214:   PetscScalar    *fvalue;          /* value of event function at the end of the step*/
215:   PetscScalar    *fvalue_prev;     /* value of event function at start of the step */
216:   PetscReal       ptime;           /* time at step end */
217:   PetscReal       ptime_prev;      /* time at step start */
218:   PetscErrorCode  (*monitor)(TS,PetscReal,Vec,PetscScalar*,void*); /* User event monitor function */
219:   PetscErrorCode  (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,void*); /* User post event function */
220:   PetscBool      *terminate;        /* 1 -> Terminate time stepping, 0 -> continue */
221:   PetscInt       *direction;        /* Zero crossing direction: 1 -> Going positive, -1 -> Going negative, 0 -> Any */
222:   PetscInt        nevents;          /* Number of events to handle */
223:   PetscInt        nevents_zero;     /* Number of event zero detected */
224:   PetscInt        *events_zero;      /* List of events that have reached zero */
225:   void           *monitorcontext;
226:   PetscReal       tol;              /* Tolerance for event zero check */
227:   TSEventStatus   status;           /* Event status */
228:   PetscReal       tstepend;         /* End time of step */
229:   PetscReal       initial_timestep; /* Initial time step */
230: };

232: PETSC_EXTERN PetscErrorCode TSEventMonitor(TS);
233: PETSC_EXTERN PetscErrorCode TSEventMonitorDestroy(TSEvent*);

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

237: typedef enum {TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
238:               TS_STEP_PENDING,    /* vec_sol advanced, but step has not been accepted yet */
239:               TS_STEP_COMPLETE    /* step accepted and ptime, steps, etc have been advanced */
240: } TSStepStatus;

242: struct _n_TSMonitorLGCtx {
243:   PetscDrawLG lg;
244:   PetscInt    howoften;  /* when > 0 uses step % howoften, when negative only final solution plotted */
245:   PetscInt    ksp_its,snes_its;
246: };

248: #endif