Actual source code: tsimpl.h
petsc-3.10.5 2019-03-28
1: #ifndef __TSIMPL_H
4: #include <petscts.h>
5: #include <petsc/private/petscimpl.h>
7: /*
8: Timesteping context.
9: General DAE: F(t,U,U_t) = 0, required Jacobian is G'(U) where G(U) = F(t,U,U0+a*U)
10: General ODE: U_t = F(t,U) <-- the right-hand-side function
11: Linear ODE: U_t = A(t) U <-- the right-hand-side matrix
12: Linear (no time) ODE: U_t = A U <-- the right-hand-side matrix
13: */
15: /*
16: Maximum number of monitors you can run with a single TS
17: */
18: #define MAXTSMONITORS 10
20: PETSC_EXTERN PetscBool TSRegisterAllCalled;
21: PETSC_EXTERN PetscErrorCode TSRegisterAll(void);
22: PETSC_EXTERN PetscErrorCode TSAdaptRegisterAll(void);
24: PETSC_EXTERN PetscErrorCode TSRKRegisterAll(void);
25: PETSC_EXTERN PetscErrorCode TSARKIMEXRegisterAll(void);
26: PETSC_EXTERN PetscErrorCode TSRosWRegisterAll(void);
27: PETSC_EXTERN PetscErrorCode TSGLLERegisterAll(void);
28: PETSC_EXTERN PetscErrorCode TSGLLEAdaptRegisterAll(void);
30: typedef struct _TSOps *TSOps;
32: struct _TSOps {
33: PetscErrorCode (*snesfunction)(SNES,Vec,Vec,TS);
34: PetscErrorCode (*snesjacobian)(SNES,Vec,Mat,Mat,TS);
35: PetscErrorCode (*setup)(TS);
36: PetscErrorCode (*step)(TS);
37: PetscErrorCode (*solve)(TS);
38: PetscErrorCode (*interpolate)(TS,PetscReal,Vec);
39: PetscErrorCode (*evaluatewlte)(TS,NormType,PetscInt*,PetscReal*);
40: PetscErrorCode (*evaluatestep)(TS,PetscInt,Vec,PetscBool*);
41: PetscErrorCode (*setfromoptions)(PetscOptionItems*,TS);
42: PetscErrorCode (*destroy)(TS);
43: PetscErrorCode (*view)(TS,PetscViewer);
44: PetscErrorCode (*reset)(TS);
45: PetscErrorCode (*linearstability)(TS,PetscReal,PetscReal,PetscReal*,PetscReal*);
46: PetscErrorCode (*load)(TS,PetscViewer);
47: PetscErrorCode (*rollback)(TS);
48: PetscErrorCode (*getstages)(TS,PetscInt*,Vec**);
49: PetscErrorCode (*adjointstep)(TS);
50: PetscErrorCode (*adjointsetup)(TS);
51: PetscErrorCode (*adjointintegral)(TS);
52: PetscErrorCode (*forwardsetup)(TS);
53: PetscErrorCode (*forwardstep)(TS);
54: PetscErrorCode (*forwardintegral)(TS);
55: PetscErrorCode (*getsolutioncomponents)(TS,PetscInt*,Vec*);
56: PetscErrorCode (*getauxsolution)(TS,Vec*);
57: PetscErrorCode (*gettimeerror)(TS,PetscInt,Vec*);
58: PetscErrorCode (*settimeerror)(TS,Vec);
59: PetscErrorCode (*startingmethod) (TS);
60: };
62: /*
63: TSEvent - Abstract object to handle event monitoring
64: */
65: typedef struct _n_TSEvent *TSEvent;
67: typedef struct _TSTrajectoryOps *TSTrajectoryOps;
69: struct _TSTrajectoryOps {
70: PetscErrorCode (*view)(TSTrajectory,PetscViewer);
71: PetscErrorCode (*reset)(TSTrajectory);
72: PetscErrorCode (*destroy)(TSTrajectory);
73: PetscErrorCode (*set)(TSTrajectory,TS,PetscInt,PetscReal,Vec);
74: PetscErrorCode (*get)(TSTrajectory,TS,PetscInt,PetscReal*);
75: PetscErrorCode (*setfromoptions)(PetscOptionItems*,TSTrajectory);
76: PetscErrorCode (*setup)(TSTrajectory,TS);
77: };
79: struct _p_TSTrajectory {
80: PETSCHEADER(struct _TSTrajectoryOps);
81: PetscViewer monitor;
82: PetscInt setupcalled; /* true if setup has been called */
83: PetscInt recomps; /* counter for recomputations in the adjoint run */
84: PetscInt diskreads,diskwrites; /* counters for disk checkpoint reads and writes */
85: char **names; /* the name of each variable; each process has only the local names */
86: PetscBool keepfiles; /* keep the files generated during the run after the run is complete */
87: char *dirname,*filetemplate; /* directory name and file name template for disk checkpoints */
88: char *dirfiletemplate; /* complete directory and file name template for disk checkpoints */
89: PetscErrorCode (*transform)(void*,Vec,Vec*);
90: PetscErrorCode (*transformdestroy)(void*);
91: void* transformctx;
92: void *data;
93: };
95: typedef struct _TS_RHSSplitLink *TS_RHSSplitLink;
96: struct _TS_RHSSplitLink {
97: TS ts;
98: char *splitname;
99: IS is;
100: TS_RHSSplitLink next;
101: PetscLogEvent event;
102: };
104: struct _p_TS {
105: PETSCHEADER(struct _TSOps);
106: TSProblemType problem_type;
107: TSEquationType equation_type;
109: DM dm;
110: Vec vec_sol; /* solution vector in first and second order equations */
111: Vec vec_dot; /* time derivative vector in second order equations */
112: TSAdapt adapt;
113: TSAdaptType default_adapt_type;
114: TSEvent event;
116: /* ---------------- User (or PETSc) Provided stuff ---------------------*/
117: PetscErrorCode (*monitor[MAXTSMONITORS])(TS,PetscInt,PetscReal,Vec,void*);
118: PetscErrorCode (*monitordestroy[MAXTSMONITORS])(void**);
119: void *monitorcontext[MAXTSMONITORS];
120: PetscInt numbermonitors;
121: PetscErrorCode (*adjointmonitor[MAXTSMONITORS])(TS,PetscInt,PetscReal,Vec,PetscInt,Vec*,Vec*,void*);
122: PetscErrorCode (*adjointmonitordestroy[MAXTSMONITORS])(void**);
123: void *adjointmonitorcontext[MAXTSMONITORS];
124: PetscInt numberadjointmonitors;
126: PetscErrorCode (*prestep)(TS);
127: PetscErrorCode (*prestage)(TS,PetscReal);
128: PetscErrorCode (*poststage)(TS,PetscReal,PetscInt,Vec*);
129: PetscErrorCode (*postevaluate)(TS);
130: PetscErrorCode (*poststep)(TS);
131: PetscErrorCode (*functiondomainerror)(TS,PetscReal,Vec,PetscBool*);
133: /* ---------------------- Sensitivity Analysis support -----------------*/
134: TSTrajectory trajectory; /* All solutions are kept here for the entire time integration process */
135: Vec *vecs_sensi; /* one vector for each cost function */
136: Vec *vecs_sensip;
137: PetscInt numcost; /* number of cost functions */
138: Vec vec_costintegral;
139: PetscInt adjointsetupcalled;
140: PetscInt adjoint_steps;
141: PetscInt adjoint_max_steps;
142: PetscBool adjoint_solve; /* immediately call TSAdjointSolve() after TSSolve() is complete */
143: PetscBool costintegralfwd; /* cost integral is evaluated in the forward run if true */
144: Vec vec_costintegrand; /* workspace for Adjoint computations */
145: Mat Jacp;
146: void *rhsjacobianpctx;
147: void *costintegrandctx;
148: Vec *vecs_drdy;
149: Vec *vecs_drdp;
151: PetscErrorCode (*rhsjacobianp)(TS,PetscReal,Vec,Mat,void*);
152: PetscErrorCode (*costintegrand)(TS,PetscReal,Vec,Vec,void*);
153: PetscErrorCode (*drdyfunction)(TS,PetscReal,Vec,Vec*,void*);
154: PetscErrorCode (*drdpfunction)(TS,PetscReal,Vec,Vec*,void*);
156: /* specific to forward sensitivity analysis */
157: Mat mat_sensip; /* matrix storing forward sensitivities */
158: Vec *vecs_integral_sensip; /* one vector for each integral */
159: PetscInt num_parameters;
160: PetscInt num_initialvalues;
161: void *vecsrhsjacobianpctx;
162: PetscInt forwardsetupcalled;
163: PetscBool forward_solve;
164: PetscErrorCode (*vecsrhsjacobianp)(TS,PetscReal,Vec,Vec*,void*);
166: /* ---------------------- IMEX support ---------------------------------*/
167: /* These extra slots are only used when the user provides both Implicit and RHS */
168: Mat Arhs; /* Right hand side matrix */
169: Mat Brhs; /* Right hand side preconditioning matrix */
170: Vec Frhs; /* Right hand side function value */
172: /* This is a general caching scheme to avoid recomputing the Jacobian at a place that has been previously been evaluated.
173: * The present use case is that TSComputeRHSFunctionLinear() evaluates the Jacobian once and we don't want it to be immeditely re-evaluated.
174: */
175: struct {
176: PetscReal time; /* The time at which the matrices were last evaluated */
177: PetscObjectId Xid; /* Unique ID of solution vector at which the Jacobian was last evaluated */
178: PetscObjectState Xstate; /* State of the solution vector */
179: MatStructure mstructure; /* The structure returned */
180: /* Flag to unshift Jacobian before calling the IJacobian or RHSJacobian functions. This is useful
181: * if the user would like to reuse (part of) the Jacobian from the last evaluation. */
182: PetscBool reuse;
183: PetscReal scale,shift;
184: } rhsjacobian;
186: struct {
187: PetscReal shift; /* The derivative of the lhs wrt to Xdot */
188: } ijacobian;
190: /* --------------------Nonlinear Iteration------------------------------*/
191: SNES snes;
192: PetscBool usessnes; /* Flag set by each TSType to indicate if the type actually uses a SNES;
193: this works around the design flaw that a SNES is ALWAYS created with TS even when it is not needed.*/
194: PetscInt ksp_its; /* total number of linear solver iterations */
195: PetscInt snes_its; /* total number of nonlinear solver iterations */
196: PetscInt num_snes_failures;
197: PetscInt max_snes_failures;
199: /* --- Data that is unique to each particular solver --- */
200: PetscInt setupcalled; /* true if setup has been called */
201: void *data; /* implementationspecific data */
202: void *user; /* user context */
204: /* ------------------ Parameters -------------------------------------- */
205: PetscInt max_steps; /* max number of steps */
206: PetscReal max_time; /* max time allowed */
208: /* --------------------------------------------------------------------- */
210: PetscBool steprollback; /* flag to indicate that the step was rolled back */
211: PetscBool steprestart; /* flag to indicate that the timestepper has to discard any history and restart */
212: PetscInt steps; /* steps taken so far in all successive calls to TSSolve() */
213: PetscReal ptime; /* time at the start of the current step (stage time is internal if it exists) */
214: PetscReal time_step; /* current time increment */
215: PetscReal ptime_prev; /* time at the start of the previous step */
216: PetscReal ptime_prev_rollback; /* time at the start of the 2nd previous step to recover from rollback */
217: PetscReal solvetime; /* time at the conclusion of TSSolve() */
219: TSConvergedReason reason;
220: PetscBool errorifstepfailed;
221: PetscInt reject,max_reject;
222: TSExactFinalTimeOption exact_final_time;
224: PetscReal atol,rtol; /* Relative and absolute tolerance for local truncation error */
225: Vec vatol,vrtol; /* Relative and absolute tolerance in vector form */
226: PetscReal cfltime,cfltime_local;
228: PetscBool testjacobian;
229: PetscBool testjacobiantranspose;
230: /* ------------------- Default work-area management ------------------ */
231: PetscInt nwork;
232: Vec *work;
234: /* ---------------------- RHS splitting support ---------------------------------*/
235: PetscInt num_rhs_splits;
236: TS_RHSSplitLink tsrhssplit;
237: };
239: struct _TSAdaptOps {
240: PetscErrorCode (*choose)(TSAdapt,TS,PetscReal,PetscInt*,PetscReal*,PetscBool*,PetscReal*,PetscReal*,PetscReal*);
241: PetscErrorCode (*destroy)(TSAdapt);
242: PetscErrorCode (*reset)(TSAdapt);
243: PetscErrorCode (*view)(TSAdapt,PetscViewer);
244: PetscErrorCode (*setfromoptions)(PetscOptionItems*,TSAdapt);
245: PetscErrorCode (*load)(TSAdapt,PetscViewer);
246: };
248: struct _p_TSAdapt {
249: PETSCHEADER(struct _TSAdaptOps);
250: void *data;
251: PetscErrorCode (*checkstage)(TSAdapt,TS,PetscReal,Vec,PetscBool*);
252: struct {
253: PetscInt n; /* number of candidate schemes, including the one currently in use */
254: PetscBool inuse_set; /* the current scheme has been set */
255: const char *name[16]; /* name of the scheme */
256: PetscInt order[16]; /* classical order of each scheme */
257: PetscInt stageorder[16]; /* stage order of each scheme */
258: PetscReal ccfl[16]; /* stability limit relative to explicit Euler */
259: PetscReal cost[16]; /* relative measure of the amount of work required for each scheme */
260: } candidates;
261: PetscBool always_accept;
262: PetscReal safety; /* safety factor relative to target error/stability goal */
263: PetscReal reject_safety; /* extra safety factor if the last step was rejected */
264: PetscReal clip[2]; /* admissible time step decrease/increase factors */
265: PetscReal dt_min,dt_max; /* admissible minimum and maximum time step */
266: PetscReal scale_solve_failed; /* scale step by this factor if solver (linear or nonlinear) fails. */
267: NormType wnormtype;
268: PetscViewer monitor;
269: PetscInt timestepjustincreased;
270: };
272: typedef struct _p_DMTS *DMTS;
273: typedef struct _DMTSOps *DMTSOps;
274: struct _DMTSOps {
275: TSRHSFunction rhsfunction;
276: TSRHSJacobian rhsjacobian;
278: TSIFunction ifunction;
279: PetscErrorCode (*ifunctionview)(void*,PetscViewer);
280: PetscErrorCode (*ifunctionload)(void**,PetscViewer);
282: TSIJacobian ijacobian;
283: PetscErrorCode (*ijacobianview)(void*,PetscViewer);
284: PetscErrorCode (*ijacobianload)(void**,PetscViewer);
286: TSI2Function i2function;
287: TSI2Jacobian i2jacobian;
289: TSSolutionFunction solution;
290: TSForcingFunction forcing;
292: PetscErrorCode (*destroy)(DMTS);
293: PetscErrorCode (*duplicate)(DMTS,DMTS);
294: };
296: struct _p_DMTS {
297: PETSCHEADER(struct _DMTSOps);
298: void *rhsfunctionctx;
299: void *rhsjacobianctx;
301: void *ifunctionctx;
302: void *ijacobianctx;
304: void *i2functionctx;
305: void *i2jacobianctx;
307: void *solutionctx;
308: void *forcingctx;
310: void *data;
312: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
313: * copy-on-write. When DMGetDMTSWrite() sees a request using a different DM, it makes a copy. Thus, if a user
314: * only interacts directly with one level, e.g., using TSSetIFunction(), then coarse levels of a multilevel item
315: * integrator are built, then the user changes the routine with another call to TSSetIFunction(), it automatically
316: * propagates to all the levels. If instead, they get out a specific level and set the function on that level,
317: * subsequent changes to the original level will no longer propagate to that level.
318: */
319: DM originaldm;
320: };
322: PETSC_EXTERN PetscErrorCode DMGetDMTS(DM,DMTS*);
323: PETSC_EXTERN PetscErrorCode DMGetDMTSWrite(DM,DMTS*);
324: PETSC_EXTERN PetscErrorCode DMCopyDMTS(DM,DM);
325: PETSC_EXTERN PetscErrorCode DMTSView(DMTS,PetscViewer);
326: PETSC_EXTERN PetscErrorCode DMTSLoad(DMTS,PetscViewer);
327: PETSC_EXTERN PetscErrorCode DMTSCopy(DMTS,DMTS);
329: typedef enum {TSEVENT_NONE,TSEVENT_LOCATED_INTERVAL,TSEVENT_PROCESSING,TSEVENT_ZERO,TSEVENT_RESET_NEXTSTEP} TSEventStatus;
331: struct _n_TSEvent {
332: PetscScalar *fvalue; /* value of event function at the end of the step*/
333: PetscScalar *fvalue_prev; /* value of event function at start of the step (left end-point of event interval) */
334: PetscReal ptime_prev; /* time at step start (left end-point of event interval) */
335: PetscReal ptime_end; /* end time of step (when an event interval is detected, ptime_end is fixed to the time at step end during event processing) */
336: PetscReal ptime_right; /* time on the right end-point of the event interval */
337: PetscScalar *fvalue_right; /* value of event function at the right end-point of the event interval */
338: PetscInt *side; /* Used for detecting repetition of end-point, -1 => left, +1 => right */
339: PetscReal timestep_prev; /* previous time step */
340: PetscReal timestep_orig; /* initial time step */
341: PetscBool *zerocrossing; /* Flag to signal zero crossing detection */
342: PetscErrorCode (*eventhandler)(TS,PetscReal,Vec,PetscScalar*,void*); /* User event handler function */
343: PetscErrorCode (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,PetscBool,void*); /* User post event function */
344: void *ctx; /* User context for event handler and post even functions */
345: PetscInt *direction; /* Zero crossing direction: 1 -> Going positive, -1 -> Going negative, 0 -> Any */
346: PetscBool *terminate; /* 1 -> Terminate time stepping, 0 -> continue */
347: PetscInt nevents; /* Number of events to handle */
348: PetscInt nevents_zero; /* Number of event zero detected */
349: PetscInt *events_zero; /* List of events that have reached zero */
350: PetscReal *vtol; /* Vector tolerances for event zero check */
351: TSEventStatus status; /* Event status */
352: PetscInt iterctr; /* Iteration counter */
353: PetscViewer monitor;
354: /* Struct to record the events */
355: struct {
356: PetscInt ctr; /* recorder counter */
357: PetscReal *time; /* Event times */
358: PetscInt *stepnum; /* Step numbers */
359: PetscInt *nevents; /* Number of events occuring at the event times */
360: PetscInt **eventidx; /* Local indices of the events in the event list */
361: } recorder;
362: PetscInt recsize; /* Size of recorder stack */
363: PetscInt refct; /* reference count */
364: };
366: PETSC_EXTERN PetscErrorCode TSEventInitialize(TSEvent,TS,PetscReal,Vec);
367: PETSC_EXTERN PetscErrorCode TSEventDestroy(TSEvent*);
368: PETSC_EXTERN PetscErrorCode TSEventHandler(TS);
369: PETSC_EXTERN PetscErrorCode TSAdjointEventHandler(TS);
371: PETSC_EXTERN PetscLogEvent TS_AdjointStep;
372: PETSC_EXTERN PetscLogEvent TS_Step;
373: PETSC_EXTERN PetscLogEvent TS_PseudoComputeTimeStep;
374: PETSC_EXTERN PetscLogEvent TS_FunctionEval;
375: PETSC_EXTERN PetscLogEvent TS_JacobianEval;
376: PETSC_EXTERN PetscLogEvent TS_ForwardStep;
378: typedef enum {TS_STEP_INCOMPLETE, /* vec_sol, ptime, etc point to beginning of step */
379: TS_STEP_PENDING, /* vec_sol advanced, but step has not been accepted yet */
380: TS_STEP_COMPLETE /* step accepted and ptime, steps, etc have been advanced */
381: } TSStepStatus;
383: struct _n_TSMonitorLGCtx {
384: PetscDrawLG lg;
385: PetscBool semilogy;
386: PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
387: PetscInt ksp_its,snes_its;
388: char **names;
389: char **displaynames;
390: PetscInt ndisplayvariables;
391: PetscInt *displayvariables;
392: PetscReal *displayvalues;
393: PetscErrorCode (*transform)(void*,Vec,Vec*);
394: PetscErrorCode (*transformdestroy)(void*);
395: void *transformctx;
396: };
398: struct _n_TSMonitorEnvelopeCtx {
399: Vec max,min;
400: };
402: /*
403: Checks if the user provide a TSSetIFunction() but an explicit method is called; generate an error in that case
404: */
405: PETSC_STATIC_INLINE PetscErrorCode TSCheckImplicitTerm(TS ts)
406: {
407: TSIFunction ifunction;
408: DM dm;
409: PetscErrorCode ierr;
412: TSGetDM(ts,&dm);
413: DMTSGetIFunction(dm,&ifunction,NULL);
414: if (ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"You are attempting to use an explicit ODE integrator but provided an implicit function definition with TSSetIFunction()");
415: return(0);
416: }
418: PETSC_EXTERN PetscLogEvent TSTrajectory_Set;
419: PETSC_EXTERN PetscLogEvent TSTrajectory_Get;
420: PETSC_EXTERN PetscLogEvent TSTrajectory_DiskWrite;
421: PETSC_EXTERN PetscLogEvent TSTrajectory_DiskRead;
423: struct _n_TSMonitorDrawCtx {
424: PetscViewer viewer;
425: Vec initialsolution;
426: PetscBool showinitial;
427: PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */
428: PetscBool showtimestepandtime;
429: };
430: #endif