1*26bd1501SBarry Smith #if !defined(PETSCGL_H) 2*26bd1501SBarry Smith #define PETSCGL_H 33d177a5cSEmil Constantinescu 43d177a5cSEmil Constantinescu #include <petsc/private/tsimpl.h> 53d177a5cSEmil Constantinescu 63d177a5cSEmil Constantinescu typedef enum {TSGLLEERROR_FORWARD,TSGLLEERROR_BACKWARD} TSGLLEErrorDirection; 73d177a5cSEmil Constantinescu 83d177a5cSEmil Constantinescu typedef struct _TSGLLEScheme *TSGLLEScheme; 93d177a5cSEmil Constantinescu struct _TSGLLEScheme { 103d177a5cSEmil Constantinescu PetscInt p; /* order of the method */ 113d177a5cSEmil Constantinescu PetscInt q; /* stage-order of the method */ 123d177a5cSEmil Constantinescu PetscInt r; /* number of items carried between stages */ 133d177a5cSEmil Constantinescu PetscInt s; /* number of stages */ 143d177a5cSEmil Constantinescu PetscScalar *c; /* location of the stages */ 153d177a5cSEmil Constantinescu PetscScalar *a,*b,*u,*v; /* tableau for the method */ 163d177a5cSEmil Constantinescu 173d177a5cSEmil Constantinescu /* For use in rescale & modify */ 183d177a5cSEmil Constantinescu PetscScalar *alpha; /* X_n(t_n) - X_{n-1}(t_n) = - alpha^T h^{p+1} x^{(p+1)}(t_n) */ 193d177a5cSEmil Constantinescu PetscScalar *beta; /* - beta^T h^{p+2} x^{(p+2)}(t_n) */ 203d177a5cSEmil Constantinescu PetscScalar *gamma; /* - gamma^T h^{p+2} f' x^{(p+1)}(t_n) + O(h^{p+3}) */ 213d177a5cSEmil Constantinescu 223d177a5cSEmil Constantinescu /* Error estimates */ 233d177a5cSEmil Constantinescu /* h^{p+1}x^{(p+1)}(t_n) ~= phi[0]*h*Ydot + psi[0]*X[1:] */ 243d177a5cSEmil Constantinescu /* h^{p+2}x^{(p+2)}(t_n) ~= phi[1]*h*Ydot + psi[1]*X[1:] */ 253d177a5cSEmil Constantinescu /* h^{p+2}f' x^{(p+1)}(t_n) ~= phi[2]*h*Ydot + psi[2]*X[1:] */ 263d177a5cSEmil Constantinescu PetscScalar *phi; /* dim=[3][s] for estimating higher moments, see B,J,W 2007 */ 273d177a5cSEmil Constantinescu PetscScalar *psi; /* dim=[3][r-1], [0 psi^T] of B,J,W 2007 */ 283d177a5cSEmil Constantinescu PetscScalar *stage_error; 293d177a5cSEmil Constantinescu 303d177a5cSEmil Constantinescu /* Desirable properties which enable extra optimizations */ 313d177a5cSEmil Constantinescu PetscBool stiffly_accurate; /* Last row of [A U] is equal t first row of [B V]? */ 323d177a5cSEmil Constantinescu PetscBool fsal; /* First Same As Last: X[1] = h*Ydot[s-1] (and stiffly accurate) */ 333d177a5cSEmil Constantinescu }; 343d177a5cSEmil Constantinescu 353d177a5cSEmil Constantinescu typedef struct TS_GLLE { 363d177a5cSEmil Constantinescu TSGLLEAcceptFunction Accept; /* Decides whether to accept a given time step, given estimates of local truncation error */ 373d177a5cSEmil Constantinescu TSGLLEAdapt adapt; 383d177a5cSEmil Constantinescu 393d177a5cSEmil Constantinescu /* These names are only stored so that they can be printed in TSView_GLLE() without making these schemes full-blown 403d177a5cSEmil Constantinescu objects (the implementations I'm thinking of do not have state and I'm lazy). */ 413d177a5cSEmil Constantinescu char accept_name[256]; 423d177a5cSEmil Constantinescu 433d177a5cSEmil Constantinescu /* specific to the family of GL method */ 443d177a5cSEmil Constantinescu PetscErrorCode (*EstimateHigherMoments)(TSGLLEScheme,PetscReal,Vec*,Vec*,Vec*); /* Provide local error estimates */ 453d177a5cSEmil Constantinescu PetscErrorCode (*CompleteStep)(TSGLLEScheme,PetscReal,TSGLLEScheme,PetscReal,Vec*,Vec*,Vec*); 463d177a5cSEmil Constantinescu PetscErrorCode (*Destroy)(struct TS_GLLE*); 473d177a5cSEmil Constantinescu PetscErrorCode (*View)(struct TS_GLLE*,PetscViewer); 483d177a5cSEmil Constantinescu char type_name[256]; 493d177a5cSEmil Constantinescu PetscInt nschemes; 503d177a5cSEmil Constantinescu TSGLLEScheme *schemes; 513d177a5cSEmil Constantinescu 523d177a5cSEmil Constantinescu Vec *X; /* Items to carry between steps */ 533d177a5cSEmil Constantinescu Vec *Xold; /* Values of these items at the last step */ 543d177a5cSEmil Constantinescu Vec W; /* = 1/(atol+rtol*|X0|), used for WRMS norm */ 553d177a5cSEmil 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)} */ 563d177a5cSEmil Constantinescu PetscReal wrms_atol,wrms_rtol; 573d177a5cSEmil Constantinescu 583d177a5cSEmil Constantinescu /* Stages (Y,Ydot) are computed sequentially */ 593d177a5cSEmil Constantinescu Vec *Ydot; /* Derivatives of stage vectors, must be stored */ 603d177a5cSEmil Constantinescu Vec Y; /* Stage vector, only used while solving the stage so we don't need to store it */ 613d177a5cSEmil Constantinescu Vec Z; /* Affine vector */ 623d177a5cSEmil Constantinescu PetscReal scoeff; /* Ydot = Z + shift*Y; shift = scoeff/ts->time_step */ 633d177a5cSEmil Constantinescu PetscReal stage_time; /* time at current stage */ 643d177a5cSEmil Constantinescu PetscInt stage; /* index of the stage we are currently solving for */ 653d177a5cSEmil Constantinescu 663d177a5cSEmil Constantinescu /* Runtime options */ 673d177a5cSEmil Constantinescu PetscInt current_scheme; 683d177a5cSEmil Constantinescu PetscInt max_order,min_order,start_order; 693d177a5cSEmil Constantinescu PetscBool extrapolate; /* use extrapolation to produce initial Newton iterate? */ 703d177a5cSEmil Constantinescu TSGLLEErrorDirection error_direction; /* TSGLLEERROR_FORWARD or TSGLLEERROR_BACKWARD */ 713d177a5cSEmil Constantinescu 723d177a5cSEmil Constantinescu PetscInt max_step_rejections; 733d177a5cSEmil Constantinescu 743d177a5cSEmil Constantinescu PetscBool setupcalled; 753d177a5cSEmil Constantinescu void *data; 763d177a5cSEmil Constantinescu } TS_GLLE; 773d177a5cSEmil Constantinescu 783d177a5cSEmil Constantinescu #endif 79