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