4: #include <petsc/private/tsimpl.h>
6: typedef enum {TSGLERROR_FORWARD,TSGLERROR_BACKWARD} TSGLErrorDirection;
8: typedef struct _TSGLScheme *TSGLScheme;
9: struct _TSGLScheme {
10: PetscInt p; /* order of the method */
11: PetscInt q; /* stage-order of the method */
12: PetscInt r; /* number of items carried between stages */
13: PetscInt s; /* number of stages */
14: PetscScalar *c; /* location of the stages */
15: PetscScalar *a,*b,*u,*v; /* tableau for the method */
17: /* For use in rescale & modify */
18: PetscScalar *alpha; /* X_n(t_n) - X_{n-1}(t_n) = - alpha^T h^{p+1} x^{(p+1)}(t_n) */
19: PetscScalar *beta; /* - beta^T h^{p+2} x^{(p+2)}(t_n) */
20: PetscScalar *gamma; /* - gamma^T h^{p+2} f' x^{(p+1)}(t_n) + O(h^{p+3}) */
22: /* Error estimates */
23: /* h^{p+1}x^{(p+1)}(t_n) ~= phi[0]*h*Ydot + psi[0]*X[1:] */
24: /* h^{p+2}x^{(p+2)}(t_n) ~= phi[1]*h*Ydot + psi[1]*X[1:] */
25: /* h^{p+2}f' x^{(p+1)}(t_n) ~= phi[2]*h*Ydot + psi[2]*X[1:] */
26: PetscScalar *phi; /* dim=[3][s] for estimating higher moments, see B,J,W 2007 */
27: PetscScalar *psi; /* dim=[3][r-1], [0 psi^T] of B,J,W 2007 */
28: PetscScalar *stage_error;
30: /* Desirable properties which enable extra optimizations */
31: PetscBool stiffly_accurate; /* Last row of [A U] is equal t first row of [B V]? */
32: PetscBool fsal; /* First Same As Last: X[1] = h*Ydot[s-1] (and stiffly accurate) */
33: };
35: typedef struct TS_GL {
36: TSGLAcceptFunction Accept; /* Decides whether to accept a given time step, given estimates of local truncation error */
37: TSGLAdapt adapt;
39: /* These names are only stored so that they can be printed in TSView_GL() without making these schemes full-blown
40: objects (the implementations I'm thinking of do not have state and I'm lazy). */
41: char accept_name[256];
43: /* specific to the family of GL method */
44: PetscErrorCode (*EstimateHigherMoments)(TSGLScheme,PetscReal,Vec*,Vec*,Vec*); /* Provide local error estimates */
45: PetscErrorCode (*CompleteStep)(TSGLScheme,PetscReal,TSGLScheme,PetscReal,Vec*,Vec*,Vec*);
46: PetscErrorCode (*Destroy)(struct TS_GL*);
47: PetscErrorCode (*View)(struct TS_GL*,PetscViewer);
48: char type_name[256];
49: PetscInt nschemes;
50: TSGLScheme *schemes;
52: Vec *X; /* Items to carry between steps */
53: Vec *Xold; /* Values of these items at the last step */
54: Vec W; /* = 1/(atol+rtol*|X0|), used for WRMS norm */
55: Vec *himom; /* len=3, Estimates of h^{p+1}x^{(p+1)}, h^{p+2}x^{(p+2)}, h^{p+2}(df/dx) x^{(p+1)} */
56: PetscReal wrms_atol,wrms_rtol;
58: /* Stages (Y,Ydot) are computed sequentially */
59: Vec *Ydot; /* Derivatives of stage vectors, must be stored */
60: Vec Y; /* Stage vector, only used while solving the stage so we don't need to store it */
61: Vec Z; /* Affine vector */
62: PetscReal scoeff; /* Ydot = Z + shift*Y; shift = scoeff/ts->time_step */
63: PetscReal stage_time; /* time at current stage */
64: PetscInt stage; /* index of the stage we are currently solving for */
66: /* Runtime options */
67: PetscInt current_scheme;
68: PetscInt max_order,min_order,start_order;
69: PetscBool extrapolate; /* use extrapolation to produce initial Newton iterate? */
70: TSGLErrorDirection error_direction; /* TSGLERROR_FORWARD or TSGLERROR_BACKWARD */
72: PetscInt max_step_rejections;
74: PetscBool setupcalled;
75: void *data;
76: } TS_GL;
78: #endif