Actual source code: tsimpl.h
petsc-3.6.4 2016-04-12
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