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