xref: /petsc/src/ts/impls/implicit/glle/glle.h (revision d5b43468fb8780a8feea140ccd6fa3e6a50411cc)
1 #ifndef PETSCGL_H
2 #define PETSCGL_H
3 
4 #include <petsc/private/tsimpl.h>
5 
6 typedef enum {
7   TSGLLEERROR_FORWARD,
8   TSGLLEERROR_BACKWARD
9 } TSGLLEErrorDirection;
10 
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 */
19 
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}) */
24 
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;
32 
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 };
37 
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;
41 
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];
45 
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;
54 
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;
60 
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 */
68 
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 */
74 
75   PetscInt max_step_rejections;
76 
77   PetscBool setupcalled;
78   void     *data;
79 } TS_GLLE;
80 
81 #endif
82