1e27a552bSJed Brown /* 261692a83SJed Brown Code for timestepping with Rosenbrock W methods 3e27a552bSJed Brown 4e27a552bSJed Brown Notes: 5e27a552bSJed Brown The general system is written as 6e27a552bSJed Brown 7f9c1d6abSBarry Smith F(t,U,Udot) = G(t,U) 8e27a552bSJed Brown 9f9c1d6abSBarry Smith where F represents the stiff part of the physics and G represents the non-stiff part. 10f9c1d6abSBarry Smith This method is designed to be linearly implicit on F and can use an approximate and lagged Jacobian. 11e27a552bSJed Brown 12e27a552bSJed Brown */ 13b45d2f2cSJed Brown #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 14e27a552bSJed Brown 1561692a83SJed Brown #include <../src/mat/blockinvert.h> 1661692a83SJed Brown 1719fd82e9SBarry Smith static TSRosWType TSRosWDefault = TSROSWRA34PW2; 18e27a552bSJed Brown static PetscBool TSRosWRegisterAllCalled; 19e27a552bSJed Brown static PetscBool TSRosWPackageInitialized; 20e27a552bSJed Brown 2161692a83SJed Brown typedef struct _RosWTableau *RosWTableau; 2261692a83SJed Brown struct _RosWTableau { 23e27a552bSJed Brown char *name; 24e27a552bSJed Brown PetscInt order; /* Classical approximation order of the method */ 25e27a552bSJed Brown PetscInt s; /* Number of stages */ 26f4aed992SEmil Constantinescu PetscInt pinterp; /* Interpolation order */ 2761692a83SJed Brown PetscReal *A; /* Propagation table, strictly lower triangular */ 2861692a83SJed Brown PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */ 29c17803e7SJed Brown PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */ 3043b21953SEmil Constantinescu PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/ 3161692a83SJed Brown PetscReal *b; /* Step completion table */ 32fe7e6d57SJed Brown PetscReal *bembed; /* Step completion table for embedded method of order one less */ 3361692a83SJed Brown PetscReal *ASum; /* Row sum of A */ 3461692a83SJed Brown PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */ 3561692a83SJed Brown PetscReal *At; /* Propagation table in transformed variables */ 3661692a83SJed Brown PetscReal *bt; /* Step completion table in transformed variables */ 37fe7e6d57SJed Brown PetscReal *bembedt; /* Step completion table of order one less in transformed variables */ 3861692a83SJed Brown PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */ 398d59e960SJed Brown PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 403ca35412SEmil Constantinescu PetscReal *binterpt; /* Dense output formula */ 41e27a552bSJed Brown }; 4261692a83SJed Brown typedef struct _RosWTableauLink *RosWTableauLink; 4361692a83SJed Brown struct _RosWTableauLink { 4461692a83SJed Brown struct _RosWTableau tab; 4561692a83SJed Brown RosWTableauLink next; 46e27a552bSJed Brown }; 4761692a83SJed Brown static RosWTableauLink RosWTableauList; 48e27a552bSJed Brown 49e27a552bSJed Brown typedef struct { 5061692a83SJed Brown RosWTableau tableau; 5161692a83SJed Brown Vec *Y; /* States computed during the step, used to complete the step */ 52e27a552bSJed Brown Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 5361692a83SJed Brown Vec Ystage; /* Work vector for the state value at each stage */ 5461692a83SJed Brown Vec Zdot; /* Ydot = Zdot + shift*Y */ 5561692a83SJed Brown Vec Zstage; /* Y = Zstage + Y */ 563ca35412SEmil Constantinescu Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation)*/ 571c3436cfSJed Brown PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 58b296d7d5SJed Brown PetscReal scoeff; /* shift = scoeff/dt */ 59e27a552bSJed Brown PetscReal stage_time; 60c17803e7SJed Brown PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */ 6161692a83SJed Brown PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 62108c343cSJed Brown TSStepStatus status; 63e27a552bSJed Brown } TS_RosW; 64e27a552bSJed Brown 65fe7e6d57SJed Brown /*MC 663606a31eSEmil Constantinescu TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method). 673606a31eSEmil Constantinescu 683606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 693606a31eSEmil Constantinescu 703606a31eSEmil Constantinescu Level: intermediate 713606a31eSEmil Constantinescu 723606a31eSEmil Constantinescu .seealso: TSROSW 733606a31eSEmil Constantinescu M*/ 743606a31eSEmil Constantinescu 753606a31eSEmil Constantinescu /*MC 763606a31eSEmil Constantinescu TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method). 773606a31eSEmil Constantinescu 783606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 793606a31eSEmil Constantinescu 803606a31eSEmil Constantinescu Level: intermediate 813606a31eSEmil Constantinescu 823606a31eSEmil Constantinescu .seealso: TSROSW 833606a31eSEmil Constantinescu M*/ 843606a31eSEmil Constantinescu 853606a31eSEmil Constantinescu /*MC 86fe7e6d57SJed Brown TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme. 87fe7e6d57SJed Brown 88fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P. 89fe7e6d57SJed Brown 90fe7e6d57SJed Brown Level: intermediate 91fe7e6d57SJed Brown 92fe7e6d57SJed Brown .seealso: TSROSW 93fe7e6d57SJed Brown M*/ 94fe7e6d57SJed Brown 95fe7e6d57SJed Brown /*MC 96fe7e6d57SJed Brown TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme. 97fe7e6d57SJed Brown 98fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M. 99fe7e6d57SJed Brown 100fe7e6d57SJed Brown Level: intermediate 101fe7e6d57SJed Brown 102fe7e6d57SJed Brown .seealso: TSROSW 103fe7e6d57SJed Brown M*/ 104fe7e6d57SJed Brown 105fe7e6d57SJed Brown /*MC 106fe7e6d57SJed Brown TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1. 107fe7e6d57SJed Brown 108fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 109fe7e6d57SJed Brown 110fe7e6d57SJed Brown This is strongly A-stable with R(infty) = 0.73. The embedded method of order 2 is strongly A-stable with R(infty) = 0.73. 111fe7e6d57SJed Brown 112fe7e6d57SJed Brown References: 113fe7e6d57SJed Brown Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 114fe7e6d57SJed Brown 115fe7e6d57SJed Brown Level: intermediate 116fe7e6d57SJed Brown 117fe7e6d57SJed Brown .seealso: TSROSW 118fe7e6d57SJed Brown M*/ 119fe7e6d57SJed Brown 120fe7e6d57SJed Brown /*MC 121fe7e6d57SJed Brown TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1. 122fe7e6d57SJed Brown 123fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 124fe7e6d57SJed Brown 125fe7e6d57SJed Brown This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48. 126fe7e6d57SJed Brown 127fe7e6d57SJed Brown References: 128fe7e6d57SJed Brown Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 129fe7e6d57SJed Brown 130fe7e6d57SJed Brown Level: intermediate 131fe7e6d57SJed Brown 132fe7e6d57SJed Brown .seealso: TSROSW 133fe7e6d57SJed Brown M*/ 134fe7e6d57SJed Brown 135ef3c5b88SJed Brown /*MC 136ef3c5b88SJed Brown TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme 137ef3c5b88SJed Brown 138ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 139ef3c5b88SJed Brown 140ef3c5b88SJed Brown Both the third order and embedded second order methods are stiffly accurate and L-stable. 141ef3c5b88SJed Brown 142ef3c5b88SJed Brown References: 143ef3c5b88SJed Brown Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 144ef3c5b88SJed Brown 145ef3c5b88SJed Brown Level: intermediate 146ef3c5b88SJed Brown 147ef3c5b88SJed Brown .seealso: TSROSW, TSROSWSANDU3 148ef3c5b88SJed Brown M*/ 149ef3c5b88SJed Brown 150ef3c5b88SJed Brown /*MC 151ef3c5b88SJed Brown TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme 152ef3c5b88SJed Brown 153ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 154ef3c5b88SJed Brown 155ef3c5b88SJed Brown The third order method is L-stable, but not stiffly accurate. 156ef3c5b88SJed Brown The second order embedded method is strongly A-stable with R(infty) = 0.5. 157ef3c5b88SJed Brown The internal stages are L-stable. 158ef3c5b88SJed Brown This method is called ROS3 in the paper. 159ef3c5b88SJed Brown 160ef3c5b88SJed Brown References: 161ef3c5b88SJed Brown Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 162ef3c5b88SJed Brown 163ef3c5b88SJed Brown Level: intermediate 164ef3c5b88SJed Brown 165ef3c5b88SJed Brown .seealso: TSROSW, TSROSWRODAS3 166ef3c5b88SJed Brown M*/ 167ef3c5b88SJed Brown 168961f28d0SJed Brown /*MC 169961f28d0SJed Brown TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 170961f28d0SJed Brown 171961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 172961f28d0SJed Brown 173961f28d0SJed Brown A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 174961f28d0SJed Brown 175961f28d0SJed Brown References: 176961f28d0SJed Brown Emil Constantinescu 177961f28d0SJed Brown 178961f28d0SJed Brown Level: intermediate 179961f28d0SJed Brown 18043b21953SEmil Constantinescu .seealso: TSROSW, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, SSP 181961f28d0SJed Brown M*/ 182961f28d0SJed Brown 183961f28d0SJed Brown /*MC 184998eb97aSJed Brown TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 185961f28d0SJed Brown 186961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 187961f28d0SJed Brown 188961f28d0SJed Brown L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 189961f28d0SJed Brown 190961f28d0SJed Brown References: 191961f28d0SJed Brown Emil Constantinescu 192961f28d0SJed Brown 193961f28d0SJed Brown Level: intermediate 194961f28d0SJed Brown 19543b21953SEmil Constantinescu .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLLSSP3P4S2C, TSSSP 196961f28d0SJed Brown M*/ 197961f28d0SJed Brown 198961f28d0SJed Brown /*MC 199998eb97aSJed Brown TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 200961f28d0SJed Brown 201961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 202961f28d0SJed Brown 203961f28d0SJed Brown L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 204961f28d0SJed Brown 205961f28d0SJed Brown References: 206961f28d0SJed Brown Emil Constantinescu 207961f28d0SJed Brown 208961f28d0SJed Brown Level: intermediate 209961f28d0SJed Brown 210961f28d0SJed Brown .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSSSP 211961f28d0SJed Brown M*/ 212961f28d0SJed Brown 21342faf41dSJed Brown /*MC 21442faf41dSJed Brown TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop 21542faf41dSJed Brown 21642faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 21742faf41dSJed Brown 21842faf41dSJed Brown A(89.3 degrees)-stable, |R(infty)| = 0.454. 21942faf41dSJed Brown 22042faf41dSJed Brown This method does not provide a dense output formula. 22142faf41dSJed Brown 22242faf41dSJed Brown References: 22342faf41dSJed Brown Kaps and Rentrop, Generalized Runge-Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979. 22442faf41dSJed Brown 22542faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 22642faf41dSJed Brown 22742faf41dSJed Brown Hairer's code ros4.f 22842faf41dSJed Brown 22942faf41dSJed Brown Level: intermediate 23042faf41dSJed Brown 23142faf41dSJed Brown .seealso: TSROSW, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 23242faf41dSJed Brown M*/ 23342faf41dSJed Brown 23442faf41dSJed Brown /*MC 23542faf41dSJed Brown TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine 23642faf41dSJed Brown 23742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 23842faf41dSJed Brown 23942faf41dSJed Brown A-stable, |R(infty)| = 1/3. 24042faf41dSJed Brown 24142faf41dSJed Brown This method does not provide a dense output formula. 24242faf41dSJed Brown 24342faf41dSJed Brown References: 24442faf41dSJed Brown Shampine, Implementation of Rosenbrock methods, 1982. 24542faf41dSJed Brown 24642faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 24742faf41dSJed Brown 24842faf41dSJed Brown Hairer's code ros4.f 24942faf41dSJed Brown 25042faf41dSJed Brown Level: intermediate 25142faf41dSJed Brown 25242faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWVELDD4, TSROSW4L 25342faf41dSJed Brown M*/ 25442faf41dSJed Brown 25542faf41dSJed Brown /*MC 25642faf41dSJed Brown TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen 25742faf41dSJed Brown 25842faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 25942faf41dSJed Brown 26042faf41dSJed Brown A(89.5 degrees)-stable, |R(infty)| = 0.24. 26142faf41dSJed Brown 26242faf41dSJed Brown This method does not provide a dense output formula. 26342faf41dSJed Brown 26442faf41dSJed Brown References: 26542faf41dSJed Brown van Veldhuizen, D-stability and Kaps-Rentrop methods, 1984. 26642faf41dSJed Brown 26742faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 26842faf41dSJed Brown 26942faf41dSJed Brown Hairer's code ros4.f 27042faf41dSJed Brown 27142faf41dSJed Brown Level: intermediate 27242faf41dSJed Brown 27342faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 27442faf41dSJed Brown M*/ 27542faf41dSJed Brown 27642faf41dSJed Brown /*MC 27742faf41dSJed Brown TSROSW4L - four stage, fourth order Rosenbrock (not W) method 27842faf41dSJed Brown 27942faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 28042faf41dSJed Brown 28142faf41dSJed Brown A-stable and L-stable 28242faf41dSJed Brown 28342faf41dSJed Brown This method does not provide a dense output formula. 28442faf41dSJed Brown 28542faf41dSJed Brown References: 28642faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 28742faf41dSJed Brown 28842faf41dSJed Brown Hairer's code ros4.f 28942faf41dSJed Brown 29042faf41dSJed Brown Level: intermediate 29142faf41dSJed Brown 29242faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 29342faf41dSJed Brown M*/ 29442faf41dSJed Brown 295e27a552bSJed Brown #undef __FUNCT__ 296e27a552bSJed Brown #define __FUNCT__ "TSRosWRegisterAll" 297e27a552bSJed Brown /*@C 298e27a552bSJed Brown TSRosWRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSRosW 299e27a552bSJed Brown 300e27a552bSJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 301e27a552bSJed Brown 302e27a552bSJed Brown Level: advanced 303e27a552bSJed Brown 304e27a552bSJed Brown .keywords: TS, TSRosW, register, all 305e27a552bSJed Brown 306e27a552bSJed Brown .seealso: TSRosWRegisterDestroy() 307e27a552bSJed Brown @*/ 308e27a552bSJed Brown PetscErrorCode TSRosWRegisterAll(void) 309e27a552bSJed Brown { 310e27a552bSJed Brown PetscErrorCode ierr; 311e27a552bSJed Brown 312e27a552bSJed Brown PetscFunctionBegin; 313e27a552bSJed Brown if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 314e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_TRUE; 315e27a552bSJed Brown 316e27a552bSJed Brown { 3173606a31eSEmil Constantinescu const PetscReal 3183606a31eSEmil Constantinescu A = 0, 3193606a31eSEmil Constantinescu Gamma = 1, 3201f80e275SEmil Constantinescu b = 1, 3211f80e275SEmil Constantinescu binterpt=1; 3221f80e275SEmil Constantinescu 3231f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);CHKERRQ(ierr); 3243606a31eSEmil Constantinescu } 3253606a31eSEmil Constantinescu 3263606a31eSEmil Constantinescu { 3273606a31eSEmil Constantinescu const PetscReal 3283606a31eSEmil Constantinescu A= 0, 3293606a31eSEmil Constantinescu Gamma = 0.5, 3301f80e275SEmil Constantinescu b = 1, 3311f80e275SEmil Constantinescu binterpt=1; 3321f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);CHKERRQ(ierr); 3333606a31eSEmil Constantinescu } 3343606a31eSEmil Constantinescu 3353606a31eSEmil Constantinescu { 336*da80777bSKarl Rupp /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0); Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */ 337e27a552bSJed Brown const PetscReal 33861692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 339*da80777bSKarl Rupp Gamma[2][2] = {{1.707106781186547524401,0}, {-2.*1.707106781186547524401,1.707106781186547524401}}, 3401c3436cfSJed Brown b[2] = {0.5,0.5}, 3411c3436cfSJed Brown b1[2] = {1.0,0.0}; 3421f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 343*da80777bSKarl Rupp binterpt[0][0] = 1.707106781186547524401 - 1.0; 344*da80777bSKarl Rupp binterpt[1][0] = 2.0 - 1.707106781186547524401; 345*da80777bSKarl Rupp binterpt[0][1] = 1.707106781186547524401 - 1.5; 346*da80777bSKarl Rupp binterpt[1][1] = 1.5 - 1.707106781186547524401; 3471f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 348e27a552bSJed Brown } 349e27a552bSJed Brown { 350*da80777bSKarl Rupp /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0); Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */ 351e27a552bSJed Brown const PetscReal 35261692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 353*da80777bSKarl Rupp Gamma[2][2] = {{0.2928932188134524755992,0}, {-2.*0.2928932188134524755992,0.2928932188134524755992}}, 3541c3436cfSJed Brown b[2] = {0.5,0.5}, 3551c3436cfSJed Brown b1[2] = {1.0,0.0}; 3561f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 357*da80777bSKarl Rupp binterpt[0][0] = 0.2928932188134524755992 - 1.0; 358*da80777bSKarl Rupp binterpt[1][0] = 2.0 - 0.2928932188134524755992; 359*da80777bSKarl Rupp binterpt[0][1] = 0.2928932188134524755992 - 1.5; 360*da80777bSKarl Rupp binterpt[1][1] = 1.5 - 0.2928932188134524755992; 3611f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 362fe7e6d57SJed Brown } 363fe7e6d57SJed Brown { 364*da80777bSKarl Rupp /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */ 3651f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 366fe7e6d57SJed Brown const PetscReal 367fe7e6d57SJed Brown A[3][3] = {{0,0,0}, 368fe7e6d57SJed Brown {1.5773502691896257e+00,0,0}, 369fe7e6d57SJed Brown {0.5,0,0}}, 370*da80777bSKarl Rupp Gamma[3][3] = {{7.8867513459481287e-01,0,0}, 371*da80777bSKarl Rupp {-1.5773502691896257e+00,7.8867513459481287e-01,0}, 372*da80777bSKarl Rupp {-6.7075317547305480e-01,-1.7075317547305482e-01,7.8867513459481287e-01}}, 373fe7e6d57SJed Brown b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01}, 374fe7e6d57SJed Brown b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01}; 3751f80e275SEmil Constantinescu 3761f80e275SEmil Constantinescu binterpt[0][0]=-0.8094010767585034; 3771f80e275SEmil Constantinescu binterpt[1][0]=-0.5; 3781f80e275SEmil Constantinescu binterpt[2][0]=2.3094010767585034; 3791f80e275SEmil Constantinescu binterpt[0][1]=0.9641016151377548; 3801f80e275SEmil Constantinescu binterpt[1][1]=0.5; 3811f80e275SEmil Constantinescu binterpt[2][1]=-1.4641016151377548; 3821f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 383fe7e6d57SJed Brown } 384fe7e6d57SJed Brown { 3853ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 386*da80777bSKarl Rupp /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 387fe7e6d57SJed Brown const PetscReal 388fe7e6d57SJed Brown A[4][4] = {{0,0,0,0}, 389fe7e6d57SJed Brown {8.7173304301691801e-01,0,0,0}, 390fe7e6d57SJed Brown {8.4457060015369423e-01,-1.1299064236484185e-01,0,0}, 391fe7e6d57SJed Brown {0,0,1.,0}}, 392*da80777bSKarl Rupp Gamma[4][4] = {{4.3586652150845900e-01,0,0,0}, 393*da80777bSKarl Rupp {-8.7173304301691801e-01,4.3586652150845900e-01,0,0}, 394*da80777bSKarl Rupp {-9.0338057013044082e-01,5.4180672388095326e-02,4.3586652150845900e-01,0}, 395*da80777bSKarl Rupp {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,4.3586652150845900e-01}}, 396fe7e6d57SJed Brown b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01}, 3973ca35412SEmil Constantinescu b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01}; 3983ca35412SEmil Constantinescu 3993ca35412SEmil Constantinescu binterpt[0][0]=1.0564298455794094; 4003ca35412SEmil Constantinescu binterpt[1][0]=2.296429974281067; 4013ca35412SEmil Constantinescu binterpt[2][0]=-1.307599564525376; 4023ca35412SEmil Constantinescu binterpt[3][0]=-1.045260255335102; 4033ca35412SEmil Constantinescu binterpt[0][1]=-1.3864882699759573; 4043ca35412SEmil Constantinescu binterpt[1][1]=-8.262611700275677; 4053ca35412SEmil Constantinescu binterpt[2][1]=7.250979895056055; 4063ca35412SEmil Constantinescu binterpt[3][1]=2.398120075195581; 4073ca35412SEmil Constantinescu binterpt[0][2]=0.5721822314575016; 4083ca35412SEmil Constantinescu binterpt[1][2]=4.742931142090097; 4093ca35412SEmil Constantinescu binterpt[2][2]=-4.398120075195578; 4103ca35412SEmil Constantinescu binterpt[3][2]=-0.9169932983520199; 4113ca35412SEmil Constantinescu 4123ca35412SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 413e27a552bSJed Brown } 414ef3c5b88SJed Brown { 415*da80777bSKarl Rupp /* const PetscReal g = 0.5; Directly written in-place below */ 416ef3c5b88SJed Brown const PetscReal 417ef3c5b88SJed Brown A[4][4] = {{0,0,0,0}, 418ef3c5b88SJed Brown {0,0,0,0}, 419ef3c5b88SJed Brown {1.,0,0,0}, 420ef3c5b88SJed Brown {0.75,-0.25,0.5,0}}, 421*da80777bSKarl Rupp Gamma[4][4] = {{0.5,0,0,0}, 422*da80777bSKarl Rupp {1.,0.5,0,0}, 423*da80777bSKarl Rupp {-0.25,-0.25,0.5,0}, 424*da80777bSKarl Rupp {1./12,1./12,-2./3,0.5}}, 425ef3c5b88SJed Brown b[4] = {5./6,-1./6,-1./6,0.5}, 426ef3c5b88SJed Brown b2[4] = {0.75,-0.25,0.5,0}; 427f4aed992SEmil Constantinescu ierr = TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,PETSC_NULL);CHKERRQ(ierr); 428ef3c5b88SJed Brown } 429ef3c5b88SJed Brown { 430*da80777bSKarl Rupp /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */ 431ef3c5b88SJed Brown const PetscReal 432ef3c5b88SJed Brown A[3][3] = {{0,0,0}, 433*da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}, 434*da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}}, 435*da80777bSKarl Rupp Gamma[3][3] = {{0.43586652150845899941601945119356,0,0}, 436*da80777bSKarl Rupp {-0.19294655696029095575009695436041,0.43586652150845899941601945119356,0}, 437*da80777bSKarl Rupp {0,1.74927148125794685173529749738960,0.43586652150845899941601945119356}}, 438ef3c5b88SJed Brown b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829}, 439ef3c5b88SJed Brown b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619}; 4401f80e275SEmil Constantinescu 4411f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 4421f80e275SEmil Constantinescu binterpt[0][0]=3.793692883777660870425141387941; 4431f80e275SEmil Constantinescu binterpt[1][0]=-2.918692883777660870425141387941; 4441f80e275SEmil Constantinescu binterpt[2][0]=0.125; 4451f80e275SEmil Constantinescu binterpt[0][1]=-0.725741064379812106687651020584; 4461f80e275SEmil Constantinescu binterpt[1][1]=0.559074397713145440020984353917; 4471f80e275SEmil Constantinescu binterpt[2][1]=0.16666666666666666666666666666667; 4481f80e275SEmil Constantinescu 4491f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 450ef3c5b88SJed Brown } 451b1c69cc3SEmil Constantinescu { 452*da80777bSKarl Rupp /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 453*da80777bSKarl Rupp * Direct evaluation: s3 = 1.732050807568877293527; 454*da80777bSKarl Rupp * g = 0.7886751345948128822546; 455*da80777bSKarl Rupp * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */ 456b1c69cc3SEmil Constantinescu const PetscReal 457b1c69cc3SEmil Constantinescu A[3][3] = {{0,0,0}, 458b1c69cc3SEmil Constantinescu {1,0,0}, 459b1c69cc3SEmil Constantinescu {0.25,0.25,0}}, 460b1c69cc3SEmil Constantinescu Gamma[3][3] = {{0,0,0}, 461*da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/6.0,0.7886751345948128822546,0}, 462*da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/24.0,(-3.0-1.732050807568877293527)/8.0,0.7886751345948128822546}}, 463b1c69cc3SEmil Constantinescu b[3] = {1./6.,1./6.,2./3.}, 464b1c69cc3SEmil Constantinescu b2[3] = {1./4.,1./4.,1./2.}; 465c0cb691aSEmil Constantinescu PetscReal binterpt[3][2]; 466*da80777bSKarl Rupp 467c0cb691aSEmil Constantinescu binterpt[0][0]=0.089316397477040902157517886164709; 468c0cb691aSEmil Constantinescu binterpt[1][0]=-0.91068360252295909784248211383529; 469c0cb691aSEmil Constantinescu binterpt[2][0]=1.8213672050459181956849642276706; 470c0cb691aSEmil Constantinescu binterpt[0][1]=0.077350269189625764509148780501957; 471c0cb691aSEmil Constantinescu binterpt[1][1]=1.077350269189625764509148780502; 472c0cb691aSEmil Constantinescu binterpt[2][1]=-1.1547005383792515290182975610039; 473c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 474b1c69cc3SEmil Constantinescu } 475b1c69cc3SEmil Constantinescu 476b1c69cc3SEmil Constantinescu { 477b1c69cc3SEmil Constantinescu const PetscReal 478b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 479b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 480b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 481b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 482b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 483b1c69cc3SEmil Constantinescu {0.0,1./4.,0,0}, 484b1c69cc3SEmil Constantinescu {-2.,-2./3.,2./3.,0}, 485b1c69cc3SEmil Constantinescu {1./2.,5./36.,-2./9,0}}, 486b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 487b1c69cc3SEmil Constantinescu b2[4] = {1./8.,3./4.,1./8.,0}; 488c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 489*da80777bSKarl Rupp 490c0cb691aSEmil Constantinescu binterpt[0][0]=6.25; 491c0cb691aSEmil Constantinescu binterpt[1][0]=-30.25; 492c0cb691aSEmil Constantinescu binterpt[2][0]=1.75; 493c0cb691aSEmil Constantinescu binterpt[3][0]=23.25; 494c0cb691aSEmil Constantinescu binterpt[0][1]=-9.75; 495c0cb691aSEmil Constantinescu binterpt[1][1]=58.75; 496c0cb691aSEmil Constantinescu binterpt[2][1]=-3.25; 497c0cb691aSEmil Constantinescu binterpt[3][1]=-45.75; 498c0cb691aSEmil Constantinescu binterpt[0][2]=3.6666666666666666666666666666667; 499c0cb691aSEmil Constantinescu binterpt[1][2]=-28.333333333333333333333333333333; 500c0cb691aSEmil Constantinescu binterpt[2][2]=1.6666666666666666666666666666667; 501c0cb691aSEmil Constantinescu binterpt[3][2]=23.; 502c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 503b1c69cc3SEmil Constantinescu } 504b1c69cc3SEmil Constantinescu 505b1c69cc3SEmil Constantinescu { 506b1c69cc3SEmil Constantinescu const PetscReal 507b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 508b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 509b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 510b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 511b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 512b1c69cc3SEmil Constantinescu {0.0,3./4.,0,0}, 513b1c69cc3SEmil Constantinescu {-2./3.,-23./9.,2./9.,0}, 514b1c69cc3SEmil Constantinescu {1./18.,65./108.,-2./27,0}}, 515b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 516b1c69cc3SEmil Constantinescu b2[4] = {3./16.,10./16.,3./16.,0}; 517c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 518*da80777bSKarl Rupp 519c0cb691aSEmil Constantinescu binterpt[0][0]=1.6911764705882352941176470588235; 520c0cb691aSEmil Constantinescu binterpt[1][0]=3.6813725490196078431372549019608; 521c0cb691aSEmil Constantinescu binterpt[2][0]=0.23039215686274509803921568627451; 522c0cb691aSEmil Constantinescu binterpt[3][0]=-4.6029411764705882352941176470588; 523c0cb691aSEmil Constantinescu binterpt[0][1]=-0.95588235294117647058823529411765; 524c0cb691aSEmil Constantinescu binterpt[1][1]=-6.2401960784313725490196078431373; 525c0cb691aSEmil Constantinescu binterpt[2][1]=-0.31862745098039215686274509803922; 526c0cb691aSEmil Constantinescu binterpt[3][1]=7.5147058823529411764705882352941; 527c0cb691aSEmil Constantinescu binterpt[0][2]=-0.56862745098039215686274509803922; 528c0cb691aSEmil Constantinescu binterpt[1][2]=2.7254901960784313725490196078431; 529c0cb691aSEmil Constantinescu binterpt[2][2]=0.25490196078431372549019607843137; 530c0cb691aSEmil Constantinescu binterpt[3][2]=-2.4117647058823529411764705882353; 531c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 532b1c69cc3SEmil Constantinescu } 533753f8adbSEmil Constantinescu 534753f8adbSEmil Constantinescu { 535753f8adbSEmil Constantinescu PetscReal A[4][4],Gamma[4][4],b[4],b2[4]; 5363ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 537753f8adbSEmil Constantinescu 538753f8adbSEmil Constantinescu Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816; 53905e8e825SJed Brown Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0; 540753f8adbSEmil Constantinescu Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476; 541753f8adbSEmil Constantinescu Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816; 54205e8e825SJed Brown Gamma[1][2]=0; Gamma[1][3]=0; 543753f8adbSEmil Constantinescu Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903; 544753f8adbSEmil Constantinescu Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131; 545753f8adbSEmil Constantinescu Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816; 54605e8e825SJed Brown Gamma[2][3]=0; 547753f8adbSEmil Constantinescu Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783; 548753f8adbSEmil Constantinescu Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984; 549753f8adbSEmil Constantinescu Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198; 550753f8adbSEmil Constantinescu Gamma[3][3]=0; 551753f8adbSEmil Constantinescu 55205e8e825SJed Brown A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0; 553753f8adbSEmil Constantinescu A[1][0]=0.8717330430169179988320388950590125027645343373957631; 55405e8e825SJed Brown A[1][1]=0; A[1][2]=0; A[1][3]=0; 555753f8adbSEmil Constantinescu A[2][0]=0.5275890119763004115618079766722914408876108660811028; 556753f8adbSEmil Constantinescu A[2][1]=0.07241098802369958843819203208518599088698057726988732; 55705e8e825SJed Brown A[2][2]=0; A[2][3]=0; 558753f8adbSEmil Constantinescu A[3][0]=0.3990960076760701320627260685975778145384666450351314; 559753f8adbSEmil Constantinescu A[3][1]=-0.4375576546135194437228463747348862825846903771419953; 560753f8adbSEmil Constantinescu A[3][2]=1.038461646937449311660120300601880176655352737312713; 56105e8e825SJed Brown A[3][3]=0; 562753f8adbSEmil Constantinescu 563753f8adbSEmil Constantinescu b[0]=0.1876410243467238251612921333138006734899663569186926; 564753f8adbSEmil Constantinescu b[1]=-0.5952974735769549480478230473706443582188442040780541; 565753f8adbSEmil Constantinescu b[2]=0.9717899277217721234705114616271378792182450260943198; 566753f8adbSEmil Constantinescu b[3]=0.4358665215084589994160194475295062513822671686978816; 567753f8adbSEmil Constantinescu 568753f8adbSEmil Constantinescu b2[0]=0.2147402862233891404862383521089097657790734483804460; 569753f8adbSEmil Constantinescu b2[1]=-0.4851622638849390928209050538171743017757490232519684; 570753f8adbSEmil Constantinescu b2[2]=0.8687250025203875511662123688667549217531982787600080; 571753f8adbSEmil Constantinescu b2[3]=0.4016969751411624011684543450940068201770721128357014; 572753f8adbSEmil Constantinescu 5733ca35412SEmil Constantinescu binterpt[0][0]=2.2565812720167954547104627844105; 5743ca35412SEmil Constantinescu binterpt[1][0]=1.349166413351089573796243820819; 5753ca35412SEmil Constantinescu binterpt[2][0]=-2.4695174540533503758652847586647; 5763ca35412SEmil Constantinescu binterpt[3][0]=-0.13623023131453465264142184656474; 5773ca35412SEmil Constantinescu binterpt[0][1]=-3.0826699111559187902922463354557; 5783ca35412SEmil Constantinescu binterpt[1][1]=-2.4689115685996042534544925650515; 5793ca35412SEmil Constantinescu binterpt[2][1]=5.7428279814696677152129332773553; 5803ca35412SEmil Constantinescu binterpt[3][1]=-0.19124650171414467146619437684812; 5813ca35412SEmil Constantinescu binterpt[0][2]=1.0137296634858471607430756831148; 5823ca35412SEmil Constantinescu binterpt[1][2]=0.52444768167155973161042570784064; 5833ca35412SEmil Constantinescu binterpt[2][2]=-2.3015205996945452158771370439586; 5843ca35412SEmil Constantinescu binterpt[3][2]=0.76334325453713832352363565300308; 585f4aed992SEmil Constantinescu 586f73f8d2cSSatish Balay ierr = TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 587753f8adbSEmil Constantinescu } 58842faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWGRK4T,0.231,PETSC_DEFAULT,PETSC_DEFAULT,0,-0.1282612945269037e+01);CHKERRQ(ierr); 58942faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWSHAMP4,0.5,PETSC_DEFAULT,PETSC_DEFAULT,0,125./108.);CHKERRQ(ierr); 59042faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWVELDD4,0.22570811482256823492,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.355958941201148);CHKERRQ(ierr); 59142faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSW4L,0.57282,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.093502252409163);CHKERRQ(ierr); 592e27a552bSJed Brown PetscFunctionReturn(0); 593e27a552bSJed Brown } 594e27a552bSJed Brown 595d5e6173cSPeter Brune 596d5e6173cSPeter Brune 597e27a552bSJed Brown #undef __FUNCT__ 598e27a552bSJed Brown #define __FUNCT__ "TSRosWRegisterDestroy" 599e27a552bSJed Brown /*@C 600e27a552bSJed Brown TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 601e27a552bSJed Brown 602e27a552bSJed Brown Not Collective 603e27a552bSJed Brown 604e27a552bSJed Brown Level: advanced 605e27a552bSJed Brown 606e27a552bSJed Brown .keywords: TSRosW, register, destroy 607e27a552bSJed Brown .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic() 608e27a552bSJed Brown @*/ 609e27a552bSJed Brown PetscErrorCode TSRosWRegisterDestroy(void) 610e27a552bSJed Brown { 611e27a552bSJed Brown PetscErrorCode ierr; 61261692a83SJed Brown RosWTableauLink link; 613e27a552bSJed Brown 614e27a552bSJed Brown PetscFunctionBegin; 61561692a83SJed Brown while ((link = RosWTableauList)) { 61661692a83SJed Brown RosWTableau t = &link->tab; 61761692a83SJed Brown RosWTableauList = link->next; 61861692a83SJed Brown ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 61943b21953SEmil Constantinescu ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 620fe7e6d57SJed Brown ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 621f4aed992SEmil Constantinescu ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 622e27a552bSJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 623e27a552bSJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 624e27a552bSJed Brown } 625e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 626e27a552bSJed Brown PetscFunctionReturn(0); 627e27a552bSJed Brown } 628e27a552bSJed Brown 629e27a552bSJed Brown #undef __FUNCT__ 630e27a552bSJed Brown #define __FUNCT__ "TSRosWInitializePackage" 631e27a552bSJed Brown /*@C 632e27a552bSJed Brown TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 633e27a552bSJed Brown from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW() 634e27a552bSJed Brown when using static libraries. 635e27a552bSJed Brown 636e27a552bSJed Brown Input Parameter: 637e27a552bSJed Brown path - The dynamic library path, or PETSC_NULL 638e27a552bSJed Brown 639e27a552bSJed Brown Level: developer 640e27a552bSJed Brown 641e27a552bSJed Brown .keywords: TS, TSRosW, initialize, package 642e27a552bSJed Brown .seealso: PetscInitialize() 643e27a552bSJed Brown @*/ 644e27a552bSJed Brown PetscErrorCode TSRosWInitializePackage(const char path[]) 645e27a552bSJed Brown { 646e27a552bSJed Brown PetscErrorCode ierr; 647e27a552bSJed Brown 648e27a552bSJed Brown PetscFunctionBegin; 649e27a552bSJed Brown if (TSRosWPackageInitialized) PetscFunctionReturn(0); 650e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 651e27a552bSJed Brown ierr = TSRosWRegisterAll();CHKERRQ(ierr); 652e27a552bSJed Brown ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 653e27a552bSJed Brown PetscFunctionReturn(0); 654e27a552bSJed Brown } 655e27a552bSJed Brown 656e27a552bSJed Brown #undef __FUNCT__ 657e27a552bSJed Brown #define __FUNCT__ "TSRosWFinalizePackage" 658e27a552bSJed Brown /*@C 659e27a552bSJed Brown TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 660e27a552bSJed Brown called from PetscFinalize(). 661e27a552bSJed Brown 662e27a552bSJed Brown Level: developer 663e27a552bSJed Brown 664e27a552bSJed Brown .keywords: Petsc, destroy, package 665e27a552bSJed Brown .seealso: PetscFinalize() 666e27a552bSJed Brown @*/ 667e27a552bSJed Brown PetscErrorCode TSRosWFinalizePackage(void) 668e27a552bSJed Brown { 669e27a552bSJed Brown PetscErrorCode ierr; 670e27a552bSJed Brown 671e27a552bSJed Brown PetscFunctionBegin; 672e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 673e27a552bSJed Brown ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 674e27a552bSJed Brown PetscFunctionReturn(0); 675e27a552bSJed Brown } 676e27a552bSJed Brown 677e27a552bSJed Brown #undef __FUNCT__ 678e27a552bSJed Brown #define __FUNCT__ "TSRosWRegister" 679e27a552bSJed Brown /*@C 68061692a83SJed Brown TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 681e27a552bSJed Brown 682e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 683e27a552bSJed Brown 684e27a552bSJed Brown Input Parameters: 685e27a552bSJed Brown + name - identifier for method 686e27a552bSJed Brown . order - approximation order of method 687e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 68861692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 68961692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 690fe7e6d57SJed Brown . b - Step completion table (dimension s) 69142faf41dSJed Brown . bembed - Step completion table for a scheme of order one less (dimension s, PETSC_NULL if no embedded scheme is available) 692f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 69342faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 694e27a552bSJed Brown 695e27a552bSJed Brown Notes: 69661692a83SJed Brown Several Rosenbrock W methods are provided, this function is only needed to create new methods. 697e27a552bSJed Brown 698e27a552bSJed Brown Level: advanced 699e27a552bSJed Brown 700e27a552bSJed Brown .keywords: TS, register 701e27a552bSJed Brown 702e27a552bSJed Brown .seealso: TSRosW 703e27a552bSJed Brown @*/ 704f9c1d6abSBarry Smith PetscErrorCode TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 705f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 706e27a552bSJed Brown { 707e27a552bSJed Brown PetscErrorCode ierr; 70861692a83SJed Brown RosWTableauLink link; 70961692a83SJed Brown RosWTableau t; 71061692a83SJed Brown PetscInt i,j,k; 71161692a83SJed Brown PetscScalar *GammaInv; 712e27a552bSJed Brown 713e27a552bSJed Brown PetscFunctionBegin; 714fe7e6d57SJed Brown PetscValidCharPointer(name,1); 715fe7e6d57SJed Brown PetscValidPointer(A,4); 716fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 717fe7e6d57SJed Brown PetscValidPointer(b,6); 718fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 719fe7e6d57SJed Brown 720e27a552bSJed Brown ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 721e27a552bSJed Brown ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 722e27a552bSJed Brown t = &link->tab; 723e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 724e27a552bSJed Brown t->order = order; 725e27a552bSJed Brown t->s = s; 72661692a83SJed Brown ierr = PetscMalloc5(s*s,PetscReal,&t->A,s*s,PetscReal,&t->Gamma,s,PetscReal,&t->b,s,PetscReal,&t->ASum,s,PetscReal,&t->GammaSum);CHKERRQ(ierr); 72743b21953SEmil Constantinescu ierr = PetscMalloc5(s*s,PetscReal,&t->At,s,PetscReal,&t->bt,s*s,PetscReal,&t->GammaInv,s,PetscBool,&t->GammaZeroDiag,s*s,PetscReal,&t->GammaExplicitCorr);CHKERRQ(ierr); 728e27a552bSJed Brown ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 72961692a83SJed Brown ierr = PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 73043b21953SEmil Constantinescu ierr = PetscMemcpy(t->GammaExplicitCorr,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 73161692a83SJed Brown ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 732fe7e6d57SJed Brown if (bembed) { 733fe7e6d57SJed Brown ierr = PetscMalloc2(s,PetscReal,&t->bembed,s,PetscReal,&t->bembedt);CHKERRQ(ierr); 734fe7e6d57SJed Brown ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 735fe7e6d57SJed Brown } 73661692a83SJed Brown for (i=0; i<s; i++) { 73761692a83SJed Brown t->ASum[i] = 0; 73861692a83SJed Brown t->GammaSum[i] = 0; 73961692a83SJed Brown for (j=0; j<s; j++) { 74061692a83SJed Brown t->ASum[i] += A[i*s+j]; 741fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 74261692a83SJed Brown } 74361692a83SJed Brown } 74461692a83SJed Brown ierr = PetscMalloc(s*s*sizeof(PetscScalar),&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 74561692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 746fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 747fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 748fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 749c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 750fd96d5b0SEmil Constantinescu } else { 751c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 752fd96d5b0SEmil Constantinescu } 753fd96d5b0SEmil Constantinescu } 754fd96d5b0SEmil Constantinescu 75561692a83SJed Brown switch (s) { 75661692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 75796b95a6bSBarry Smith case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0);CHKERRQ(ierr); break; 75896b95a6bSBarry Smith case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0);CHKERRQ(ierr); break; 75996b95a6bSBarry Smith case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0);CHKERRQ(ierr); break; 76061692a83SJed Brown case 5: { 76161692a83SJed Brown PetscInt ipvt5[5]; 76261692a83SJed Brown MatScalar work5[5*5]; 76396b95a6bSBarry Smith ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0);CHKERRQ(ierr); break; 76461692a83SJed Brown } 76596b95a6bSBarry Smith case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0);CHKERRQ(ierr); break; 76696b95a6bSBarry Smith case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0);CHKERRQ(ierr); break; 76761692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 76861692a83SJed Brown } 76961692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 77061692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 77143b21953SEmil Constantinescu 77243b21953SEmil Constantinescu for (i=0; i<s; i++) { 77343b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 77443b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 77543b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 77643b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 77743b21953SEmil Constantinescu } 77843b21953SEmil Constantinescu } 77943b21953SEmil Constantinescu } 78043b21953SEmil Constantinescu 78161692a83SJed Brown for (i=0; i<s; i++) { 78261692a83SJed Brown for (j=0; j<s; j++) { 78361692a83SJed Brown t->At[i*s+j] = 0; 78461692a83SJed Brown for (k=0; k<s; k++) { 78561692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 78661692a83SJed Brown } 78761692a83SJed Brown } 78861692a83SJed Brown t->bt[i] = 0; 78961692a83SJed Brown for (j=0; j<s; j++) { 79061692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 79161692a83SJed Brown } 792fe7e6d57SJed Brown if (bembed) { 793fe7e6d57SJed Brown t->bembedt[i] = 0; 794fe7e6d57SJed Brown for (j=0; j<s; j++) { 795fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 796fe7e6d57SJed Brown } 797fe7e6d57SJed Brown } 79861692a83SJed Brown } 7998d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 8008d59e960SJed Brown 801f4aed992SEmil Constantinescu t->pinterp = pinterp; 8023ca35412SEmil Constantinescu ierr = PetscMalloc(s*pinterp*sizeof(binterpt[0]),&t->binterpt);CHKERRQ(ierr); 8033ca35412SEmil Constantinescu ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 80461692a83SJed Brown link->next = RosWTableauList; 80561692a83SJed Brown RosWTableauList = link; 806e27a552bSJed Brown PetscFunctionReturn(0); 807e27a552bSJed Brown } 808e27a552bSJed Brown 809e27a552bSJed Brown #undef __FUNCT__ 81042faf41dSJed Brown #define __FUNCT__ "TSRosWRegisterRos4" 81142faf41dSJed Brown /*@C 81242faf41dSJed Brown TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing paramter choices 81342faf41dSJed Brown 81442faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 81542faf41dSJed Brown 81642faf41dSJed Brown Input Parameters: 81742faf41dSJed Brown + name - identifier for method 81842faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 81942faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 82042faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 82142faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 82242faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 82342faf41dSJed Brown . e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 82442faf41dSJed Brown 82542faf41dSJed Brown Notes: 82642faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 82742faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 82842faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 82942faf41dSJed Brown 83042faf41dSJed Brown Level: developer 83142faf41dSJed Brown 83242faf41dSJed Brown .keywords: TS, register 83342faf41dSJed Brown 83442faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 83542faf41dSJed Brown @*/ 83619fd82e9SBarry Smith PetscErrorCode TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 83742faf41dSJed Brown { 83842faf41dSJed Brown PetscErrorCode ierr; 83942faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 84042faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 84142faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 84242faf41dSJed Brown p42 = one/eight - gamma/three, 84342faf41dSJed Brown p43 = one/twelve - gamma/three, 84442faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 84542faf41dSJed Brown p56 = one/twenty - gamma/four; 84642faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 84742faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 84842faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 84942faf41dSJed Brown 85042faf41dSJed Brown PetscFunctionBegin; 85142faf41dSJed Brown /* Step 1: choose Gamma (input) */ 85242faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 85342faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 85442faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 85542faf41dSJed Brown 85642faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 85742faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 85842faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 85942faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 86042faf41dSJed Brown rhs[0] = one - b3; 86142faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 86242faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 86342faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 86442faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 86542faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 86642faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 86742faf41dSJed Brown 86842faf41dSJed Brown /* Step 3 */ 86942faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 87042faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 87142faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 87242faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 87342faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 87442faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 87542faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 87642faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 87742faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 87842faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 87942faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 88042faf41dSJed Brown 88142faf41dSJed Brown /* Step 4: back-substitute */ 88242faf41dSJed Brown beta32 = beta32beta2p / beta2p; 88342faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 88442faf41dSJed Brown 88542faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 88642faf41dSJed Brown a43 = 0; 88742faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 88842faf41dSJed Brown a42 = a32; 88942faf41dSJed Brown 89042faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 89142faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 89242faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 89342faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 89442faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 89542faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 89642faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 89742faf41dSJed Brown Gamma[3][0] = beta4p-beta42-beta43-A[3][0]; Gamma[3][1] = beta42-A[3][1]; Gamma[3][2] = beta43-A[3][2]; Gamma[3][3] = gamma; 89842faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 89942faf41dSJed Brown 90042faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 90142faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 90242faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 90342faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 90442faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 90542faf41dSJed Brown 90642faf41dSJed Brown { 90742faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 90842faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 90942faf41dSJed Brown } 91042faf41dSJed Brown ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,PETSC_NULL);CHKERRQ(ierr); 91142faf41dSJed Brown PetscFunctionReturn(0); 91242faf41dSJed Brown } 91342faf41dSJed Brown 91442faf41dSJed Brown #undef __FUNCT__ 9151c3436cfSJed Brown #define __FUNCT__ "TSEvaluateStep_RosW" 9161c3436cfSJed Brown /* 9171c3436cfSJed Brown The step completion formula is 9181c3436cfSJed Brown 9191c3436cfSJed Brown x1 = x0 + b^T Y 9201c3436cfSJed Brown 9211c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 9221c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 9231c3436cfSJed Brown 9241c3436cfSJed Brown x1e = x0 + be^T Y 9251c3436cfSJed Brown = x1 - b^T Y + be^T Y 9261c3436cfSJed Brown = x1 + (be - b)^T Y 9271c3436cfSJed Brown 9281c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9291c3436cfSJed Brown */ 930f9c1d6abSBarry Smith static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done) 9311c3436cfSJed Brown { 9321c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9331c3436cfSJed Brown RosWTableau tab = ros->tableau; 9341c3436cfSJed Brown PetscScalar *w = ros->work; 9351c3436cfSJed Brown PetscInt i; 9361c3436cfSJed Brown PetscErrorCode ierr; 9371c3436cfSJed Brown 9381c3436cfSJed Brown PetscFunctionBegin; 9391c3436cfSJed Brown if (order == tab->order) { 940108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 941f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 942de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 943f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 944f9c1d6abSBarry Smith } else {ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr);} 9451c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9461c3436cfSJed Brown PetscFunctionReturn(0); 9471c3436cfSJed Brown } else if (order == tab->order-1) { 9481c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 949108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 950f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 951de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 952f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 953108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 954108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 955f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 956f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 9571c3436cfSJed Brown } 9581c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9591c3436cfSJed Brown PetscFunctionReturn(0); 9601c3436cfSJed Brown } 9611c3436cfSJed Brown unavailable: 9621c3436cfSJed Brown if (done) *done = PETSC_FALSE; 9631c3436cfSJed Brown else SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 9641c3436cfSJed Brown PetscFunctionReturn(0); 9651c3436cfSJed Brown } 9661c3436cfSJed Brown 9671c3436cfSJed Brown #undef __FUNCT__ 968e27a552bSJed Brown #define __FUNCT__ "TSStep_RosW" 969e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 970e27a552bSJed Brown { 97161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 97261692a83SJed Brown RosWTableau tab = ros->tableau; 973e27a552bSJed Brown const PetscInt s = tab->s; 9741c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9750feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 976c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 97761692a83SJed Brown PetscScalar *w = ros->work; 9787d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 979e27a552bSJed Brown SNES snes; 9801c3436cfSJed Brown TSAdapt adapt; 9811c3436cfSJed Brown PetscInt i,j,its,lits,reject,next_scheme; 982cdbf8f93SLisandro Dalcin PetscReal next_time_step; 9831c3436cfSJed Brown PetscBool accept; 984e27a552bSJed Brown PetscErrorCode ierr; 9850feba352SEmil Constantinescu MatStructure str; 986e27a552bSJed Brown 987e27a552bSJed Brown PetscFunctionBegin; 988e27a552bSJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 989cdbf8f93SLisandro Dalcin next_time_step = ts->time_step; 9901c3436cfSJed Brown accept = PETSC_TRUE; 991108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 992e27a552bSJed Brown 99397335746SJed Brown for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 9941c3436cfSJed Brown const PetscReal h = ts->time_step; 995b8123daeSJed Brown ierr = TSPreStep(ts);CHKERRQ(ierr); 9963ca35412SEmil Constantinescu ierr = VecCopy(ts->vec_sol,ros->VecSolPrev);CHKERRQ(ierr);/*move this at the end*/ 997e27a552bSJed Brown for (i=0; i<s; i++) { 9981c3436cfSJed Brown ros->stage_time = ts->ptime + h*ASum[i]; 999b8123daeSJed Brown ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 1000c17803e7SJed Brown if (GammaZeroDiag[i]) { 1001c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 1002b296d7d5SJed Brown ros->scoeff = 1.; 1003c17803e7SJed Brown } else { 1004c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 1005b296d7d5SJed Brown ros->scoeff = 1./Gamma[i*s+i]; 1006fd96d5b0SEmil Constantinescu } 100761692a83SJed Brown 100861692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 1009de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 1010de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 101161692a83SJed Brown 101261692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 101361692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 101461692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 101561692a83SJed Brown 1016e27a552bSJed Brown /* Initial guess taken from last stage */ 101761692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 101861692a83SJed Brown 10197d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 102061692a83SJed Brown if (!ros->recompute_jacobian && !i) { 102161692a83SJed Brown ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 102261692a83SJed Brown } 102361692a83SJed Brown ierr = SNESSolve(snes,PETSC_NULL,Y[i]);CHKERRQ(ierr); 1024e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1025e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10265ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 1027ad6bc421SBarry Smith ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 102897335746SJed Brown ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 102997335746SJed Brown if (!accept) goto reject_step; 10307d4bf2deSEmil Constantinescu } else { 10311ce71dffSSatish Balay Mat J,Jp; 10320feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10330feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 103422d28d08SBarry Smith ierr = VecScale(Y[i],-1.0);CHKERRQ(ierr); 10350feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/ 10360feba352SEmil Constantinescu 10370feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10380feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10390feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 10400feba352SEmil Constantinescu /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10410feba352SEmil Constantinescu str = SAME_NONZERO_PATTERN; 1042ccbc64bcSJed Brown ierr = TSGetIJacobian(ts,&J,&Jp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 10430feba352SEmil Constantinescu ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);CHKERRQ(ierr); 104422d28d08SBarry Smith ierr = MatMult(J,Zstage,Zdot);CHKERRQ(ierr); 10450feba352SEmil Constantinescu 10460feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10470feba352SEmil Constantinescu ierr = VecScale(Y[i],h); 10485ef26d82SJed Brown ts->ksp_its += 1; 10497d4bf2deSEmil Constantinescu } 1050e27a552bSJed Brown } 10511c3436cfSJed Brown ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 1052108c343cSJed Brown ros->status = TS_STEP_PENDING; 1053e27a552bSJed Brown 10541c3436cfSJed Brown /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 1055ad6bc421SBarry Smith ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 10561c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10578d59e960SJed Brown ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 10581c3436cfSJed Brown ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 10591c3436cfSJed Brown if (accept) { 10601c3436cfSJed Brown /* ignore next_scheme for now */ 1061e27a552bSJed Brown ts->ptime += ts->time_step; 1062cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 1063e27a552bSJed Brown ts->steps++; 1064108c343cSJed Brown ros->status = TS_STEP_COMPLETE; 10651c3436cfSJed Brown break; 10661c3436cfSJed Brown } else { /* Roll back the current step */ 10671c3436cfSJed Brown for (i=0; i<s; i++) w[i] = -tab->bt[i]; 10681c3436cfSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 10691c3436cfSJed Brown ts->time_step = next_time_step; 1070108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 10711c3436cfSJed Brown } 1072476b6736SJed Brown reject_step: continue; 10731c3436cfSJed Brown } 1074b2ce242eSJed Brown if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 1075e27a552bSJed Brown PetscFunctionReturn(0); 1076e27a552bSJed Brown } 1077e27a552bSJed Brown 1078e27a552bSJed Brown #undef __FUNCT__ 1079e27a552bSJed Brown #define __FUNCT__ "TSInterpolate_RosW" 1080f9c1d6abSBarry Smith static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1081e27a552bSJed Brown { 108261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1083f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1084f4aed992SEmil Constantinescu PetscReal h; 1085f4aed992SEmil Constantinescu PetscReal tt,t; 1086f4aed992SEmil Constantinescu PetscScalar *bt; 1087f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1088f4aed992SEmil Constantinescu PetscErrorCode ierr; 1089f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1090f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1091f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1092e27a552bSJed Brown 1093e27a552bSJed Brown PetscFunctionBegin; 1094f4aed992SEmil Constantinescu if (!Bt) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1095f4aed992SEmil Constantinescu 1096f4aed992SEmil Constantinescu switch (ros->status) { 1097f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1098f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1099f4aed992SEmil Constantinescu h = ts->time_step; 1100f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1101f4aed992SEmil Constantinescu break; 1102f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1103f4aed992SEmil Constantinescu h = ts->time_step_prev; 1104f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1105f4aed992SEmil Constantinescu break; 1106f4aed992SEmil Constantinescu default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1107f4aed992SEmil Constantinescu } 11083ca35412SEmil Constantinescu ierr = PetscMalloc(s*sizeof(bt[0]),&bt);CHKERRQ(ierr); 1109f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1110f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1111f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 11123ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1113f4aed992SEmil Constantinescu } 1114f4aed992SEmil Constantinescu } 1115f4aed992SEmil Constantinescu 1116f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1117f9c1d6abSBarry Smith /*U<-0*/ 1118f9c1d6abSBarry Smith ierr = VecZeroEntries(U);CHKERRQ(ierr); 1119f4aed992SEmil Constantinescu 1120f9c1d6abSBarry Smith /*U<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11213ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j]=0; 11223ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11233ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11243ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1125f4aed992SEmil Constantinescu } 11263ca35412SEmil Constantinescu } 1127f9c1d6abSBarry Smith ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1128f4aed992SEmil Constantinescu 1129f4aed992SEmil Constantinescu /*X<-y(t) + X*/ 1130f9c1d6abSBarry Smith ierr = VecAXPY(U,1.0,ros->VecSolPrev);CHKERRQ(ierr); 1131f4aed992SEmil Constantinescu 1132f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1133f4aed992SEmil Constantinescu 1134e27a552bSJed Brown PetscFunctionReturn(0); 1135e27a552bSJed Brown } 1136e27a552bSJed Brown 1137e27a552bSJed Brown /*------------------------------------------------------------*/ 1138e27a552bSJed Brown #undef __FUNCT__ 1139e27a552bSJed Brown #define __FUNCT__ "TSReset_RosW" 1140e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1141e27a552bSJed Brown { 114261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1143e27a552bSJed Brown PetscInt s; 1144e27a552bSJed Brown PetscErrorCode ierr; 1145e27a552bSJed Brown 1146e27a552bSJed Brown PetscFunctionBegin; 114761692a83SJed Brown if (!ros->tableau) PetscFunctionReturn(0); 114861692a83SJed Brown s = ros->tableau->s; 114961692a83SJed Brown ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 115061692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 115161692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 115261692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 115361692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 11543ca35412SEmil Constantinescu ierr = VecDestroy(&ros->VecSolPrev);CHKERRQ(ierr); 115561692a83SJed Brown ierr = PetscFree(ros->work);CHKERRQ(ierr); 1156e27a552bSJed Brown PetscFunctionReturn(0); 1157e27a552bSJed Brown } 1158e27a552bSJed Brown 1159e27a552bSJed Brown #undef __FUNCT__ 1160e27a552bSJed Brown #define __FUNCT__ "TSDestroy_RosW" 1161e27a552bSJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1162e27a552bSJed Brown { 1163e27a552bSJed Brown PetscErrorCode ierr; 1164e27a552bSJed Brown 1165e27a552bSJed Brown PetscFunctionBegin; 1166e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1167e27a552bSJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 1168e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 1169e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);CHKERRQ(ierr); 117061692a83SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",PETSC_NULL);CHKERRQ(ierr); 1171e27a552bSJed Brown PetscFunctionReturn(0); 1172e27a552bSJed Brown } 1173e27a552bSJed Brown 1174d5e6173cSPeter Brune 1175d5e6173cSPeter Brune #undef __FUNCT__ 1176d5e6173cSPeter Brune #define __FUNCT__ "TSRosWGetVecs" 1177d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1178d5e6173cSPeter Brune { 1179d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1180d5e6173cSPeter Brune PetscErrorCode ierr; 1181d5e6173cSPeter Brune 1182d5e6173cSPeter Brune PetscFunctionBegin; 1183d5e6173cSPeter Brune if (Ydot) { 1184d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1185d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1186d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1187d5e6173cSPeter Brune } 1188d5e6173cSPeter Brune if (Zdot) { 1189d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1190d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1191d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1192d5e6173cSPeter Brune } 1193d5e6173cSPeter Brune if (Ystage) { 1194d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1195d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1196d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1197d5e6173cSPeter Brune } 1198d5e6173cSPeter Brune if (Zstage) { 1199d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1200d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1201d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1202d5e6173cSPeter Brune } 1203d5e6173cSPeter Brune 1204d5e6173cSPeter Brune PetscFunctionReturn(0); 1205d5e6173cSPeter Brune } 1206d5e6173cSPeter Brune 1207d5e6173cSPeter Brune 1208d5e6173cSPeter Brune #undef __FUNCT__ 1209d5e6173cSPeter Brune #define __FUNCT__ "TSRosWRestoreVecs" 1210d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1211d5e6173cSPeter Brune { 1212d5e6173cSPeter Brune PetscErrorCode ierr; 1213d5e6173cSPeter Brune 1214d5e6173cSPeter Brune PetscFunctionBegin; 1215d5e6173cSPeter Brune if (Ydot) { 1216d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1217d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1218d5e6173cSPeter Brune } 1219d5e6173cSPeter Brune } 1220d5e6173cSPeter Brune if (Zdot) { 1221d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1222d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1223d5e6173cSPeter Brune } 1224d5e6173cSPeter Brune } 1225d5e6173cSPeter Brune if (Ystage) { 1226d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1227d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1228d5e6173cSPeter Brune } 1229d5e6173cSPeter Brune } 1230d5e6173cSPeter Brune if (Zstage) { 1231d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1232d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1233d5e6173cSPeter Brune } 1234d5e6173cSPeter Brune } 1235d5e6173cSPeter Brune PetscFunctionReturn(0); 1236d5e6173cSPeter Brune } 1237d5e6173cSPeter Brune 1238d5e6173cSPeter Brune #undef __FUNCT__ 1239d5e6173cSPeter Brune #define __FUNCT__ "DMCoarsenHook_TSRosW" 1240d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1241d5e6173cSPeter Brune { 1242d5e6173cSPeter Brune 1243d5e6173cSPeter Brune PetscFunctionBegin; 1244d5e6173cSPeter Brune PetscFunctionReturn(0); 1245d5e6173cSPeter Brune } 1246d5e6173cSPeter Brune 1247d5e6173cSPeter Brune #undef __FUNCT__ 1248d5e6173cSPeter Brune #define __FUNCT__ "DMRestrictHook_TSRosW" 1249d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1250d5e6173cSPeter Brune { 1251d5e6173cSPeter Brune TS ts = (TS)ctx; 1252d5e6173cSPeter Brune PetscErrorCode ierr; 1253d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1254d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1255d5e6173cSPeter Brune 1256d5e6173cSPeter Brune PetscFunctionBegin; 1257d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1258d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1259d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1260d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1261d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1262d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1263d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1264d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1265d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1266d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1267d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1268d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1269d5e6173cSPeter Brune PetscFunctionReturn(0); 1270d5e6173cSPeter Brune } 1271d5e6173cSPeter Brune 1272258e1594SPeter Brune 1273258e1594SPeter Brune #undef __FUNCT__ 1274258e1594SPeter Brune #define __FUNCT__ "DMSubDomainHook_TSRosW" 1275258e1594SPeter Brune static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx) 1276258e1594SPeter Brune { 1277258e1594SPeter Brune 1278258e1594SPeter Brune PetscFunctionBegin; 1279258e1594SPeter Brune PetscFunctionReturn(0); 1280258e1594SPeter Brune } 1281258e1594SPeter Brune 1282258e1594SPeter Brune #undef __FUNCT__ 1283258e1594SPeter Brune #define __FUNCT__ "DMSubDomainRestrictHook_TSRosW" 1284258e1594SPeter Brune static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1285258e1594SPeter Brune { 1286258e1594SPeter Brune TS ts = (TS)ctx; 1287258e1594SPeter Brune PetscErrorCode ierr; 1288258e1594SPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1289258e1594SPeter Brune Vec Ydots,Zdots,Ystages,Zstages; 1290258e1594SPeter Brune 1291258e1594SPeter Brune PetscFunctionBegin; 1292258e1594SPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1293258e1594SPeter Brune ierr = TSRosWGetVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1294258e1594SPeter Brune 1295258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1296258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1297258e1594SPeter Brune 1298258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1299258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1300258e1594SPeter Brune 1301258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1302258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1303258e1594SPeter Brune 1304258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1305258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1306258e1594SPeter Brune 1307258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1308258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1309258e1594SPeter Brune PetscFunctionReturn(0); 1310258e1594SPeter Brune } 1311258e1594SPeter Brune 1312e27a552bSJed Brown /* 1313e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1314e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1315e27a552bSJed Brown */ 1316e27a552bSJed Brown #undef __FUNCT__ 1317e27a552bSJed Brown #define __FUNCT__ "SNESTSFormFunction_RosW" 1318f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1319e27a552bSJed Brown { 132061692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1321e27a552bSJed Brown PetscErrorCode ierr; 1322d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1323b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1324d5e6173cSPeter Brune DM dm,dmsave; 1325e27a552bSJed Brown 1326e27a552bSJed Brown PetscFunctionBegin; 1327d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1328d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1329b296d7d5SJed Brown ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1330f9c1d6abSBarry Smith ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1331d5e6173cSPeter Brune dmsave = ts->dm; 1332d5e6173cSPeter Brune ts->dm = dm; 1333d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1334d5e6173cSPeter Brune ts->dm = dmsave; 1335d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1336e27a552bSJed Brown PetscFunctionReturn(0); 1337e27a552bSJed Brown } 1338e27a552bSJed Brown 1339e27a552bSJed Brown #undef __FUNCT__ 1340e27a552bSJed Brown #define __FUNCT__ "SNESTSFormJacobian_RosW" 1341f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *str,TS ts) 1342e27a552bSJed Brown { 134361692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1344d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1345b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1346e27a552bSJed Brown PetscErrorCode ierr; 1347d5e6173cSPeter Brune DM dm,dmsave; 1348e27a552bSJed Brown 1349e27a552bSJed Brown PetscFunctionBegin; 135061692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1351d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1352d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1353d5e6173cSPeter Brune dmsave = ts->dm; 1354d5e6173cSPeter Brune ts->dm = dm; 1355b296d7d5SJed Brown ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 1356d5e6173cSPeter Brune ts->dm = dmsave; 1357d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1358e27a552bSJed Brown PetscFunctionReturn(0); 1359e27a552bSJed Brown } 1360e27a552bSJed Brown 1361e27a552bSJed Brown #undef __FUNCT__ 1362e27a552bSJed Brown #define __FUNCT__ "TSSetUp_RosW" 1363e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1364e27a552bSJed Brown { 136561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 136661692a83SJed Brown RosWTableau tab = ros->tableau; 1367e27a552bSJed Brown PetscInt s = tab->s; 1368e27a552bSJed Brown PetscErrorCode ierr; 1369d5e6173cSPeter Brune DM dm; 1370e27a552bSJed Brown 1371e27a552bSJed Brown PetscFunctionBegin; 137261692a83SJed Brown if (!ros->tableau) { 1373e27a552bSJed Brown ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1374e27a552bSJed Brown } 137561692a83SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 137661692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 137761692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 137861692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 137961692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 13803ca35412SEmil Constantinescu ierr = VecDuplicate(ts->vec_sol,&ros->VecSolPrev);CHKERRQ(ierr); 138161692a83SJed Brown ierr = PetscMalloc(s*sizeof(ros->work[0]),&ros->work);CHKERRQ(ierr); 138222d28d08SBarry Smith ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1383d5e6173cSPeter Brune if (dm) { 1384d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1385258e1594SPeter Brune ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1386d5e6173cSPeter Brune } 1387e27a552bSJed Brown PetscFunctionReturn(0); 1388e27a552bSJed Brown } 1389e27a552bSJed Brown /*------------------------------------------------------------*/ 1390e27a552bSJed Brown 1391e27a552bSJed Brown #undef __FUNCT__ 1392e27a552bSJed Brown #define __FUNCT__ "TSSetFromOptions_RosW" 1393e27a552bSJed Brown static PetscErrorCode TSSetFromOptions_RosW(TS ts) 1394e27a552bSJed Brown { 139561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1396e27a552bSJed Brown PetscErrorCode ierr; 139761692a83SJed Brown char rostype[256]; 1398e27a552bSJed Brown 1399e27a552bSJed Brown PetscFunctionBegin; 1400e27a552bSJed Brown ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 1401e27a552bSJed Brown { 140261692a83SJed Brown RosWTableauLink link; 1403e27a552bSJed Brown PetscInt count,choice; 1404e27a552bSJed Brown PetscBool flg; 1405e27a552bSJed Brown const char **namelist; 140661692a83SJed Brown SNES snes; 140761692a83SJed Brown 14088caf3d72SBarry Smith ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof(rostype));CHKERRQ(ierr); 140961692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1410e27a552bSJed Brown ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 141161692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 141261692a83SJed Brown ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 141361692a83SJed Brown ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 1414e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 141561692a83SJed Brown 141661692a83SJed Brown ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,PETSC_NULL);CHKERRQ(ierr); 141761692a83SJed Brown 141861692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 141961692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 142061692a83SJed Brown if (!((PetscObject)snes)->type_name) { 142161692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 142261692a83SJed Brown } 142361692a83SJed Brown ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1424e27a552bSJed Brown } 1425e27a552bSJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 1426e27a552bSJed Brown PetscFunctionReturn(0); 1427e27a552bSJed Brown } 1428e27a552bSJed Brown 1429e27a552bSJed Brown #undef __FUNCT__ 1430e27a552bSJed Brown #define __FUNCT__ "PetscFormatRealArray" 1431e27a552bSJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1432e27a552bSJed Brown { 1433e27a552bSJed Brown PetscErrorCode ierr; 1434e408995aSJed Brown PetscInt i; 1435e408995aSJed Brown size_t left,count; 1436e27a552bSJed Brown char *p; 1437e27a552bSJed Brown 1438e27a552bSJed Brown PetscFunctionBegin; 1439e408995aSJed Brown for (i=0,p=buf,left=len; i<n; i++) { 1440e408995aSJed Brown ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1441e27a552bSJed Brown if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1442e27a552bSJed Brown left -= count; 1443e27a552bSJed Brown p += count; 1444e27a552bSJed Brown *p++ = ' '; 1445e27a552bSJed Brown } 1446e27a552bSJed Brown p[i ? 0 : -1] = 0; 1447e27a552bSJed Brown PetscFunctionReturn(0); 1448e27a552bSJed Brown } 1449e27a552bSJed Brown 1450e27a552bSJed Brown #undef __FUNCT__ 1451e27a552bSJed Brown #define __FUNCT__ "TSView_RosW" 1452e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1453e27a552bSJed Brown { 145461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 145561692a83SJed Brown RosWTableau tab = ros->tableau; 1456e27a552bSJed Brown PetscBool iascii; 1457e27a552bSJed Brown PetscErrorCode ierr; 1458ef20d060SBarry Smith TSAdapt adapt; 1459e27a552bSJed Brown 1460e27a552bSJed Brown PetscFunctionBegin; 1461251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1462e27a552bSJed Brown if (iascii) { 146319fd82e9SBarry Smith TSRosWType rostype; 1464e408995aSJed Brown PetscInt i; 1465e408995aSJed Brown PetscReal abscissa[512]; 1466e27a552bSJed Brown char buf[512]; 146761692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 146861692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 14698caf3d72SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 147061692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1471e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 14728caf3d72SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1473e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1474e27a552bSJed Brown } 1475ad6bc421SBarry Smith ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 1476ef20d060SBarry Smith ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1477e27a552bSJed Brown ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1478e27a552bSJed Brown PetscFunctionReturn(0); 1479e27a552bSJed Brown } 1480e27a552bSJed Brown 1481e27a552bSJed Brown #undef __FUNCT__ 1482e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType" 1483e27a552bSJed Brown /*@C 148461692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1485e27a552bSJed Brown 1486e27a552bSJed Brown Logically collective 1487e27a552bSJed Brown 1488e27a552bSJed Brown Input Parameter: 1489e27a552bSJed Brown + ts - timestepping context 149061692a83SJed Brown - rostype - type of Rosenbrock-W scheme 1491e27a552bSJed Brown 1492020d8f30SJed Brown Level: beginner 1493e27a552bSJed Brown 1494020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1495e27a552bSJed Brown @*/ 149619fd82e9SBarry Smith PetscErrorCode TSRosWSetType(TS ts,TSRosWType rostype) 1497e27a552bSJed Brown { 1498e27a552bSJed Brown PetscErrorCode ierr; 1499e27a552bSJed Brown 1500e27a552bSJed Brown PetscFunctionBegin; 1501e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 150219fd82e9SBarry Smith ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,rostype));CHKERRQ(ierr); 1503e27a552bSJed Brown PetscFunctionReturn(0); 1504e27a552bSJed Brown } 1505e27a552bSJed Brown 1506e27a552bSJed Brown #undef __FUNCT__ 1507e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType" 1508e27a552bSJed Brown /*@C 150961692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1510e27a552bSJed Brown 1511e27a552bSJed Brown Logically collective 1512e27a552bSJed Brown 1513e27a552bSJed Brown Input Parameter: 1514e27a552bSJed Brown . ts - timestepping context 1515e27a552bSJed Brown 1516e27a552bSJed Brown Output Parameter: 151761692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1518e27a552bSJed Brown 1519e27a552bSJed Brown Level: intermediate 1520e27a552bSJed Brown 1521e27a552bSJed Brown .seealso: TSRosWGetType() 1522e27a552bSJed Brown @*/ 152319fd82e9SBarry Smith PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1524e27a552bSJed Brown { 1525e27a552bSJed Brown PetscErrorCode ierr; 1526e27a552bSJed Brown 1527e27a552bSJed Brown PetscFunctionBegin; 1528e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 152919fd82e9SBarry Smith ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1530e27a552bSJed Brown PetscFunctionReturn(0); 1531e27a552bSJed Brown } 1532e27a552bSJed Brown 1533e27a552bSJed Brown #undef __FUNCT__ 153461692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian" 1535e27a552bSJed Brown /*@C 153661692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1537e27a552bSJed Brown 1538e27a552bSJed Brown Logically collective 1539e27a552bSJed Brown 1540e27a552bSJed Brown Input Parameter: 1541e27a552bSJed Brown + ts - timestepping context 154261692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1543e27a552bSJed Brown 1544e27a552bSJed Brown Level: intermediate 1545e27a552bSJed Brown 1546e27a552bSJed Brown .seealso: TSRosWGetType() 1547e27a552bSJed Brown @*/ 154861692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1549e27a552bSJed Brown { 1550e27a552bSJed Brown PetscErrorCode ierr; 1551e27a552bSJed Brown 1552e27a552bSJed Brown PetscFunctionBegin; 1553e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 155461692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1555e27a552bSJed Brown PetscFunctionReturn(0); 1556e27a552bSJed Brown } 1557e27a552bSJed Brown 1558e27a552bSJed Brown EXTERN_C_BEGIN 1559e27a552bSJed Brown #undef __FUNCT__ 1560e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType_RosW" 156119fd82e9SBarry Smith PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1562e27a552bSJed Brown { 156361692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1564e27a552bSJed Brown PetscErrorCode ierr; 1565e27a552bSJed Brown 1566e27a552bSJed Brown PetscFunctionBegin; 156761692a83SJed Brown if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 156861692a83SJed Brown *rostype = ros->tableau->name; 1569e27a552bSJed Brown PetscFunctionReturn(0); 1570e27a552bSJed Brown } 1571ef20d060SBarry Smith 1572e27a552bSJed Brown #undef __FUNCT__ 1573e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType_RosW" 157419fd82e9SBarry Smith PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1575e27a552bSJed Brown { 157661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1577e27a552bSJed Brown PetscErrorCode ierr; 1578e27a552bSJed Brown PetscBool match; 157961692a83SJed Brown RosWTableauLink link; 1580e27a552bSJed Brown 1581e27a552bSJed Brown PetscFunctionBegin; 158261692a83SJed Brown if (ros->tableau) { 158361692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1584e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1585e27a552bSJed Brown } 158661692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 158761692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1588e27a552bSJed Brown if (match) { 1589e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 159061692a83SJed Brown ros->tableau = &link->tab; 1591e27a552bSJed Brown PetscFunctionReturn(0); 1592e27a552bSJed Brown } 1593e27a552bSJed Brown } 159461692a83SJed Brown SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1595e27a552bSJed Brown PetscFunctionReturn(0); 1596e27a552bSJed Brown } 159761692a83SJed Brown 1598e27a552bSJed Brown #undef __FUNCT__ 159961692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 160061692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1601e27a552bSJed Brown { 160261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1603e27a552bSJed Brown 1604e27a552bSJed Brown PetscFunctionBegin; 160561692a83SJed Brown ros->recompute_jacobian = flg; 1606e27a552bSJed Brown PetscFunctionReturn(0); 1607e27a552bSJed Brown } 1608e27a552bSJed Brown EXTERN_C_END 1609e27a552bSJed Brown 1610d5e6173cSPeter Brune 1611e27a552bSJed Brown /* ------------------------------------------------------------ */ 1612e27a552bSJed Brown /*MC 1613020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1614e27a552bSJed Brown 1615e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1616e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1617e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1618e27a552bSJed Brown 1619e27a552bSJed Brown Notes: 162061692a83SJed Brown This method currently only works with autonomous ODE and DAE. 162161692a83SJed Brown 162261692a83SJed Brown Developer notes: 162361692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 162461692a83SJed Brown 1625f9c1d6abSBarry Smith $ udot = f(u) 162661692a83SJed Brown 162761692a83SJed Brown by the stage equations 162861692a83SJed Brown 1629f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 163061692a83SJed Brown 163161692a83SJed Brown and step completion formula 163261692a83SJed Brown 1633f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 163461692a83SJed Brown 1635f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 163661692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 163761692a83SJed Brown we define new variables for the stage equations 163861692a83SJed Brown 163961692a83SJed Brown $ y_i = gamma_ij k_j 164061692a83SJed Brown 164161692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 164261692a83SJed Brown 164361692a83SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 164461692a83SJed Brown 164561692a83SJed Brown to rewrite the method as 164661692a83SJed Brown 1647f9c1d6abSBarry Smith $ [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1648f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j bt_j y_j 164961692a83SJed Brown 165061692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 165161692a83SJed Brown 165261692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 165361692a83SJed Brown 165461692a83SJed Brown or, more compactly in tensor notation 165561692a83SJed Brown 165661692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 165761692a83SJed Brown 165861692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 165961692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 166061692a83SJed Brown equation 166161692a83SJed Brown 1662f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 166361692a83SJed Brown 166461692a83SJed Brown with initial guess y_i = 0. 1665e27a552bSJed Brown 1666e27a552bSJed Brown Level: beginner 1667e27a552bSJed Brown 1668a4386c9eSJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1669a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1670e27a552bSJed Brown M*/ 1671e27a552bSJed Brown EXTERN_C_BEGIN 1672e27a552bSJed Brown #undef __FUNCT__ 1673e27a552bSJed Brown #define __FUNCT__ "TSCreate_RosW" 1674e27a552bSJed Brown PetscErrorCode TSCreate_RosW(TS ts) 1675e27a552bSJed Brown { 167661692a83SJed Brown TS_RosW *ros; 1677e27a552bSJed Brown PetscErrorCode ierr; 1678e27a552bSJed Brown 1679e27a552bSJed Brown PetscFunctionBegin; 1680e27a552bSJed Brown #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1681e27a552bSJed Brown ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1682e27a552bSJed Brown #endif 1683e27a552bSJed Brown 1684e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1685e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1686e27a552bSJed Brown ts->ops->view = TSView_RosW; 1687e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1688e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1689e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16901c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 1691e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1692e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1693e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1694e27a552bSJed Brown 169561692a83SJed Brown ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 169661692a83SJed Brown ts->data = (void*)ros; 1697e27a552bSJed Brown 1698e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 1699e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 170061692a83SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1701e27a552bSJed Brown PetscFunctionReturn(0); 1702e27a552bSJed Brown } 1703e27a552bSJed Brown EXTERN_C_END 1704