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