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 761692a83SJed Brown G(t,X,Xdot) = F(t,X) 8e27a552bSJed Brown 961692a83SJed Brown where G represents the stiff part of the physics and F represents the non-stiff part. 1061692a83SJed Brown This method is designed to be linearly implicit on G 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 17ae6f9490SJed Brown static const 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() */ 58e27a552bSJed Brown PetscReal shift; 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 { 33661692a83SJed Brown const PetscReal g = 1. + 1./PetscSqrtReal(2.0); 337e27a552bSJed Brown const PetscReal 33861692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 33961692a83SJed Brown Gamma[2][2] = {{g,0}, {-2.*g,g}}, 3401c3436cfSJed Brown b[2] = {0.5,0.5}, 3411c3436cfSJed Brown b1[2] = {1.0,0.0}; 3421f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 3431f80e275SEmil Constantinescu binterpt[0][0]=g-1.0; 3441f80e275SEmil Constantinescu binterpt[1][0]=2.0-g; 3451f80e275SEmil Constantinescu binterpt[0][1]=g-1.5; 3461f80e275SEmil Constantinescu binterpt[1][1]=1.5-g; 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 { 35061692a83SJed Brown const PetscReal g = 1. - 1./PetscSqrtReal(2.0); 351e27a552bSJed Brown const PetscReal 35261692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 35361692a83SJed Brown Gamma[2][2] = {{g,0}, {-2.*g,g}}, 3541c3436cfSJed Brown b[2] = {0.5,0.5}, 3551c3436cfSJed Brown b1[2] = {1.0,0.0}; 3561f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 3571f80e275SEmil Constantinescu binterpt[0][0]=g-1.0; 3581f80e275SEmil Constantinescu binterpt[1][0]=2.0-g; 3591f80e275SEmil Constantinescu binterpt[0][1]=g-1.5; 3601f80e275SEmil Constantinescu binterpt[1][1]=1.5-g; 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 { 364fe7e6d57SJed Brown const PetscReal g = 7.8867513459481287e-01; 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}}, 370fe7e6d57SJed Brown Gamma[3][3] = {{g,0,0}, 371fe7e6d57SJed Brown {-1.5773502691896257e+00,g,0}, 37225833a80SEmil Constantinescu {-6.7075317547305480e-01,-1.7075317547305482e-01,g}}, 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]; 386fe7e6d57SJed Brown const PetscReal g = 4.3586652150845900e-01; 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}}, 392fe7e6d57SJed Brown Gamma[4][4] = {{g,0,0,0}, 393fe7e6d57SJed Brown {-8.7173304301691801e-01,g,0,0}, 394fe7e6d57SJed Brown {-9.0338057013044082e-01,5.4180672388095326e-02,g,0}, 395fe7e6d57SJed Brown {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,g}}, 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 { 415ef3c5b88SJed Brown const PetscReal g = 0.5; 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}}, 421ef3c5b88SJed Brown Gamma[4][4] = {{g,0,0,0}, 422ef3c5b88SJed Brown {1.,g,0,0}, 423ef3c5b88SJed Brown {-0.25,-0.25,g,0}, 424ef3c5b88SJed Brown {1./12,1./12,-2./3,g}}, 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 { 430ef3c5b88SJed Brown const PetscReal g = 0.43586652150845899941601945119356; 431ef3c5b88SJed Brown const PetscReal 432ef3c5b88SJed Brown A[3][3] = {{0,0,0}, 433ef3c5b88SJed Brown {g,0,0}, 434ef3c5b88SJed Brown {g,0,0}}, 435ef3c5b88SJed Brown Gamma[3][3] = {{g,0,0}, 436ef3c5b88SJed Brown {-0.19294655696029095575009695436041,g,0}, 437ef3c5b88SJed Brown {0,1.74927148125794685173529749738960,g}}, 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 { 452aaf9cf16SJed Brown const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 453b1c69cc3SEmil Constantinescu const PetscReal 454b1c69cc3SEmil Constantinescu A[3][3] = {{0,0,0}, 455b1c69cc3SEmil Constantinescu {1,0,0}, 456b1c69cc3SEmil Constantinescu {0.25,0.25,0}}, 457b1c69cc3SEmil Constantinescu Gamma[3][3] = {{0,0,0}, 458aaf9cf16SJed Brown {(-3.0-s3)/6.0,g,0}, 459aaf9cf16SJed Brown {(-3.0-s3)/24.0,(-3.0-s3)/8.0,g}}, 460b1c69cc3SEmil Constantinescu b[3] = {1./6.,1./6.,2./3.}, 461b1c69cc3SEmil Constantinescu b2[3] = {1./4.,1./4.,1./2.}; 462c0cb691aSEmil Constantinescu 463c0cb691aSEmil Constantinescu PetscReal binterpt[3][2]; 464c0cb691aSEmil Constantinescu binterpt[0][0]=0.089316397477040902157517886164709; 465c0cb691aSEmil Constantinescu binterpt[1][0]=-0.91068360252295909784248211383529; 466c0cb691aSEmil Constantinescu binterpt[2][0]=1.8213672050459181956849642276706; 467c0cb691aSEmil Constantinescu binterpt[0][1]=0.077350269189625764509148780501957; 468c0cb691aSEmil Constantinescu binterpt[1][1]=1.077350269189625764509148780502; 469c0cb691aSEmil Constantinescu binterpt[2][1]=-1.1547005383792515290182975610039; 470c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 471b1c69cc3SEmil Constantinescu } 472b1c69cc3SEmil Constantinescu 473b1c69cc3SEmil Constantinescu { 474b1c69cc3SEmil Constantinescu const PetscReal 475b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 476b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 477b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 478b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 479b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 480b1c69cc3SEmil Constantinescu {0.0,1./4.,0,0}, 481b1c69cc3SEmil Constantinescu {-2.,-2./3.,2./3.,0}, 482b1c69cc3SEmil Constantinescu {1./2.,5./36.,-2./9,0}}, 483b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 484b1c69cc3SEmil Constantinescu b2[4] = {1./8.,3./4.,1./8.,0}; 485c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 486c0cb691aSEmil Constantinescu binterpt[0][0]=6.25; 487c0cb691aSEmil Constantinescu binterpt[1][0]=-30.25; 488c0cb691aSEmil Constantinescu binterpt[2][0]=1.75; 489c0cb691aSEmil Constantinescu binterpt[3][0]=23.25; 490c0cb691aSEmil Constantinescu binterpt[0][1]=-9.75; 491c0cb691aSEmil Constantinescu binterpt[1][1]=58.75; 492c0cb691aSEmil Constantinescu binterpt[2][1]=-3.25; 493c0cb691aSEmil Constantinescu binterpt[3][1]=-45.75; 494c0cb691aSEmil Constantinescu binterpt[0][2]=3.6666666666666666666666666666667; 495c0cb691aSEmil Constantinescu binterpt[1][2]=-28.333333333333333333333333333333; 496c0cb691aSEmil Constantinescu binterpt[2][2]=1.6666666666666666666666666666667; 497c0cb691aSEmil Constantinescu binterpt[3][2]=23.; 498c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 499b1c69cc3SEmil Constantinescu } 500b1c69cc3SEmil Constantinescu 501b1c69cc3SEmil Constantinescu { 502b1c69cc3SEmil Constantinescu const PetscReal 503b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 504b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 505b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 506b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 507b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 508b1c69cc3SEmil Constantinescu {0.0,3./4.,0,0}, 509b1c69cc3SEmil Constantinescu {-2./3.,-23./9.,2./9.,0}, 510b1c69cc3SEmil Constantinescu {1./18.,65./108.,-2./27,0}}, 511b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 512b1c69cc3SEmil Constantinescu b2[4] = {3./16.,10./16.,3./16.,0}; 513c0cb691aSEmil Constantinescu 514c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 515c0cb691aSEmil Constantinescu binterpt[0][0]=1.6911764705882352941176470588235; 516c0cb691aSEmil Constantinescu binterpt[1][0]=3.6813725490196078431372549019608; 517c0cb691aSEmil Constantinescu binterpt[2][0]=0.23039215686274509803921568627451; 518c0cb691aSEmil Constantinescu binterpt[3][0]=-4.6029411764705882352941176470588; 519c0cb691aSEmil Constantinescu binterpt[0][1]=-0.95588235294117647058823529411765; 520c0cb691aSEmil Constantinescu binterpt[1][1]=-6.2401960784313725490196078431373; 521c0cb691aSEmil Constantinescu binterpt[2][1]=-0.31862745098039215686274509803922; 522c0cb691aSEmil Constantinescu binterpt[3][1]=7.5147058823529411764705882352941; 523c0cb691aSEmil Constantinescu binterpt[0][2]=-0.56862745098039215686274509803922; 524c0cb691aSEmil Constantinescu binterpt[1][2]=2.7254901960784313725490196078431; 525c0cb691aSEmil Constantinescu binterpt[2][2]=0.25490196078431372549019607843137; 526c0cb691aSEmil Constantinescu binterpt[3][2]=-2.4117647058823529411764705882353; 527c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 528b1c69cc3SEmil Constantinescu } 529753f8adbSEmil Constantinescu 530753f8adbSEmil Constantinescu { 531753f8adbSEmil Constantinescu PetscReal A[4][4],Gamma[4][4],b[4],b2[4]; 5323ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 533753f8adbSEmil Constantinescu 534753f8adbSEmil Constantinescu Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816; 53505e8e825SJed Brown Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0; 536753f8adbSEmil Constantinescu Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476; 537753f8adbSEmil Constantinescu Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816; 53805e8e825SJed Brown Gamma[1][2]=0; Gamma[1][3]=0; 539753f8adbSEmil Constantinescu Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903; 540753f8adbSEmil Constantinescu Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131; 541753f8adbSEmil Constantinescu Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816; 54205e8e825SJed Brown Gamma[2][3]=0; 543753f8adbSEmil Constantinescu Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783; 544753f8adbSEmil Constantinescu Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984; 545753f8adbSEmil Constantinescu Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198; 546753f8adbSEmil Constantinescu Gamma[3][3]=0; 547753f8adbSEmil Constantinescu 54805e8e825SJed Brown A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0; 549753f8adbSEmil Constantinescu A[1][0]=0.8717330430169179988320388950590125027645343373957631; 55005e8e825SJed Brown A[1][1]=0; A[1][2]=0; A[1][3]=0; 551753f8adbSEmil Constantinescu A[2][0]=0.5275890119763004115618079766722914408876108660811028; 552753f8adbSEmil Constantinescu A[2][1]=0.07241098802369958843819203208518599088698057726988732; 55305e8e825SJed Brown A[2][2]=0; A[2][3]=0; 554753f8adbSEmil Constantinescu A[3][0]=0.3990960076760701320627260685975778145384666450351314; 555753f8adbSEmil Constantinescu A[3][1]=-0.4375576546135194437228463747348862825846903771419953; 556753f8adbSEmil Constantinescu A[3][2]=1.038461646937449311660120300601880176655352737312713; 55705e8e825SJed Brown A[3][3]=0; 558753f8adbSEmil Constantinescu 559753f8adbSEmil Constantinescu b[0]=0.1876410243467238251612921333138006734899663569186926; 560753f8adbSEmil Constantinescu b[1]=-0.5952974735769549480478230473706443582188442040780541; 561753f8adbSEmil Constantinescu b[2]=0.9717899277217721234705114616271378792182450260943198; 562753f8adbSEmil Constantinescu b[3]=0.4358665215084589994160194475295062513822671686978816; 563753f8adbSEmil Constantinescu 564753f8adbSEmil Constantinescu b2[0]=0.2147402862233891404862383521089097657790734483804460; 565753f8adbSEmil Constantinescu b2[1]=-0.4851622638849390928209050538171743017757490232519684; 566753f8adbSEmil Constantinescu b2[2]=0.8687250025203875511662123688667549217531982787600080; 567753f8adbSEmil Constantinescu b2[3]=0.4016969751411624011684543450940068201770721128357014; 568753f8adbSEmil Constantinescu 5693ca35412SEmil Constantinescu binterpt[0][0]=2.2565812720167954547104627844105; 5703ca35412SEmil Constantinescu binterpt[1][0]=1.349166413351089573796243820819; 5713ca35412SEmil Constantinescu binterpt[2][0]=-2.4695174540533503758652847586647; 5723ca35412SEmil Constantinescu binterpt[3][0]=-0.13623023131453465264142184656474; 5733ca35412SEmil Constantinescu binterpt[0][1]=-3.0826699111559187902922463354557; 5743ca35412SEmil Constantinescu binterpt[1][1]=-2.4689115685996042534544925650515; 5753ca35412SEmil Constantinescu binterpt[2][1]=5.7428279814696677152129332773553; 5763ca35412SEmil Constantinescu binterpt[3][1]=-0.19124650171414467146619437684812; 5773ca35412SEmil Constantinescu binterpt[0][2]=1.0137296634858471607430756831148; 5783ca35412SEmil Constantinescu binterpt[1][2]=0.52444768167155973161042570784064; 5793ca35412SEmil Constantinescu binterpt[2][2]=-2.3015205996945452158771370439586; 5803ca35412SEmil Constantinescu binterpt[3][2]=0.76334325453713832352363565300308; 581f4aed992SEmil Constantinescu 582f73f8d2cSSatish Balay ierr = TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 583753f8adbSEmil Constantinescu } 58442faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWGRK4T,0.231,PETSC_DEFAULT,PETSC_DEFAULT,0,-0.1282612945269037e+01);CHKERRQ(ierr); 58542faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWSHAMP4,0.5,PETSC_DEFAULT,PETSC_DEFAULT,0,125./108.);CHKERRQ(ierr); 58642faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWVELDD4,0.22570811482256823492,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.355958941201148);CHKERRQ(ierr); 58742faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSW4L,0.57282,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.093502252409163);CHKERRQ(ierr); 588e27a552bSJed Brown PetscFunctionReturn(0); 589e27a552bSJed Brown } 590e27a552bSJed Brown 591d5e6173cSPeter Brune 592d5e6173cSPeter Brune 593e27a552bSJed Brown #undef __FUNCT__ 594e27a552bSJed Brown #define __FUNCT__ "TSRosWRegisterDestroy" 595e27a552bSJed Brown /*@C 596e27a552bSJed Brown TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 597e27a552bSJed Brown 598e27a552bSJed Brown Not Collective 599e27a552bSJed Brown 600e27a552bSJed Brown Level: advanced 601e27a552bSJed Brown 602e27a552bSJed Brown .keywords: TSRosW, register, destroy 603e27a552bSJed Brown .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic() 604e27a552bSJed Brown @*/ 605e27a552bSJed Brown PetscErrorCode TSRosWRegisterDestroy(void) 606e27a552bSJed Brown { 607e27a552bSJed Brown PetscErrorCode ierr; 60861692a83SJed Brown RosWTableauLink link; 609e27a552bSJed Brown 610e27a552bSJed Brown PetscFunctionBegin; 61161692a83SJed Brown while ((link = RosWTableauList)) { 61261692a83SJed Brown RosWTableau t = &link->tab; 61361692a83SJed Brown RosWTableauList = link->next; 61461692a83SJed Brown ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 61543b21953SEmil Constantinescu ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 616fe7e6d57SJed Brown ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 617f4aed992SEmil Constantinescu ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 618e27a552bSJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 619e27a552bSJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 620e27a552bSJed Brown } 621e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 622e27a552bSJed Brown PetscFunctionReturn(0); 623e27a552bSJed Brown } 624e27a552bSJed Brown 625e27a552bSJed Brown #undef __FUNCT__ 626e27a552bSJed Brown #define __FUNCT__ "TSRosWInitializePackage" 627e27a552bSJed Brown /*@C 628e27a552bSJed Brown TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 629e27a552bSJed Brown from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW() 630e27a552bSJed Brown when using static libraries. 631e27a552bSJed Brown 632e27a552bSJed Brown Input Parameter: 633e27a552bSJed Brown path - The dynamic library path, or PETSC_NULL 634e27a552bSJed Brown 635e27a552bSJed Brown Level: developer 636e27a552bSJed Brown 637e27a552bSJed Brown .keywords: TS, TSRosW, initialize, package 638e27a552bSJed Brown .seealso: PetscInitialize() 639e27a552bSJed Brown @*/ 640e27a552bSJed Brown PetscErrorCode TSRosWInitializePackage(const char path[]) 641e27a552bSJed Brown { 642e27a552bSJed Brown PetscErrorCode ierr; 643e27a552bSJed Brown 644e27a552bSJed Brown PetscFunctionBegin; 645e27a552bSJed Brown if (TSRosWPackageInitialized) PetscFunctionReturn(0); 646e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 647e27a552bSJed Brown ierr = TSRosWRegisterAll();CHKERRQ(ierr); 648e27a552bSJed Brown ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 649e27a552bSJed Brown PetscFunctionReturn(0); 650e27a552bSJed Brown } 651e27a552bSJed Brown 652e27a552bSJed Brown #undef __FUNCT__ 653e27a552bSJed Brown #define __FUNCT__ "TSRosWFinalizePackage" 654e27a552bSJed Brown /*@C 655e27a552bSJed Brown TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 656e27a552bSJed Brown called from PetscFinalize(). 657e27a552bSJed Brown 658e27a552bSJed Brown Level: developer 659e27a552bSJed Brown 660e27a552bSJed Brown .keywords: Petsc, destroy, package 661e27a552bSJed Brown .seealso: PetscFinalize() 662e27a552bSJed Brown @*/ 663e27a552bSJed Brown PetscErrorCode TSRosWFinalizePackage(void) 664e27a552bSJed Brown { 665e27a552bSJed Brown PetscErrorCode ierr; 666e27a552bSJed Brown 667e27a552bSJed Brown PetscFunctionBegin; 668e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 669e27a552bSJed Brown ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 670e27a552bSJed Brown PetscFunctionReturn(0); 671e27a552bSJed Brown } 672e27a552bSJed Brown 673e27a552bSJed Brown #undef __FUNCT__ 674e27a552bSJed Brown #define __FUNCT__ "TSRosWRegister" 675e27a552bSJed Brown /*@C 67661692a83SJed Brown TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 677e27a552bSJed Brown 678e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 679e27a552bSJed Brown 680e27a552bSJed Brown Input Parameters: 681e27a552bSJed Brown + name - identifier for method 682e27a552bSJed Brown . order - approximation order of method 683e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 68461692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 68561692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 686fe7e6d57SJed Brown . b - Step completion table (dimension s) 68742faf41dSJed Brown . bembed - Step completion table for a scheme of order one less (dimension s, PETSC_NULL if no embedded scheme is available) 688f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 68942faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 690e27a552bSJed Brown 691e27a552bSJed Brown Notes: 69261692a83SJed Brown Several Rosenbrock W methods are provided, this function is only needed to create new methods. 693e27a552bSJed Brown 694e27a552bSJed Brown Level: advanced 695e27a552bSJed Brown 696e27a552bSJed Brown .keywords: TS, register 697e27a552bSJed Brown 698e27a552bSJed Brown .seealso: TSRosW 699e27a552bSJed Brown @*/ 700e27a552bSJed Brown PetscErrorCode TSRosWRegister(const TSRosWType name,PetscInt order,PetscInt s, 701f4aed992SEmil Constantinescu const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 702f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 703e27a552bSJed Brown { 704e27a552bSJed Brown PetscErrorCode ierr; 70561692a83SJed Brown RosWTableauLink link; 70661692a83SJed Brown RosWTableau t; 70761692a83SJed Brown PetscInt i,j,k; 70861692a83SJed Brown PetscScalar *GammaInv; 709e27a552bSJed Brown 710e27a552bSJed Brown PetscFunctionBegin; 711fe7e6d57SJed Brown PetscValidCharPointer(name,1); 712fe7e6d57SJed Brown PetscValidPointer(A,4); 713fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 714fe7e6d57SJed Brown PetscValidPointer(b,6); 715fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 716fe7e6d57SJed Brown 717e27a552bSJed Brown ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 718e27a552bSJed Brown ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 719e27a552bSJed Brown t = &link->tab; 720e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 721e27a552bSJed Brown t->order = order; 722e27a552bSJed Brown t->s = s; 72361692a83SJed 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); 72443b21953SEmil 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); 725e27a552bSJed Brown ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 72661692a83SJed Brown ierr = PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 72743b21953SEmil Constantinescu ierr = PetscMemcpy(t->GammaExplicitCorr,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 72861692a83SJed Brown ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 729fe7e6d57SJed Brown if (bembed) { 730fe7e6d57SJed Brown ierr = PetscMalloc2(s,PetscReal,&t->bembed,s,PetscReal,&t->bembedt);CHKERRQ(ierr); 731fe7e6d57SJed Brown ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 732fe7e6d57SJed Brown } 73361692a83SJed Brown for (i=0; i<s; i++) { 73461692a83SJed Brown t->ASum[i] = 0; 73561692a83SJed Brown t->GammaSum[i] = 0; 73661692a83SJed Brown for (j=0; j<s; j++) { 73761692a83SJed Brown t->ASum[i] += A[i*s+j]; 738fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 73961692a83SJed Brown } 74061692a83SJed Brown } 74161692a83SJed Brown ierr = PetscMalloc(s*s*sizeof(PetscScalar),&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 74261692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 743fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 744fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 745fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 746c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 747fd96d5b0SEmil Constantinescu } else { 748c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 749fd96d5b0SEmil Constantinescu } 750fd96d5b0SEmil Constantinescu } 751fd96d5b0SEmil Constantinescu 75261692a83SJed Brown switch (s) { 75361692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 75496b95a6bSBarry Smith case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0);CHKERRQ(ierr); break; 75596b95a6bSBarry Smith case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0);CHKERRQ(ierr); break; 75696b95a6bSBarry Smith case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0);CHKERRQ(ierr); break; 75761692a83SJed Brown case 5: { 75861692a83SJed Brown PetscInt ipvt5[5]; 75961692a83SJed Brown MatScalar work5[5*5]; 76096b95a6bSBarry Smith ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0);CHKERRQ(ierr); break; 76161692a83SJed Brown } 76296b95a6bSBarry Smith case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0);CHKERRQ(ierr); break; 76396b95a6bSBarry Smith case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0);CHKERRQ(ierr); break; 76461692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 76561692a83SJed Brown } 76661692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 76761692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 76843b21953SEmil Constantinescu 76943b21953SEmil Constantinescu for (i=0; i<s; i++) { 77043b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 77143b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 77243b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 77343b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 77443b21953SEmil Constantinescu } 77543b21953SEmil Constantinescu } 77643b21953SEmil Constantinescu } 77743b21953SEmil Constantinescu 77861692a83SJed Brown for (i=0; i<s; i++) { 77961692a83SJed Brown for (j=0; j<s; j++) { 78061692a83SJed Brown t->At[i*s+j] = 0; 78161692a83SJed Brown for (k=0; k<s; k++) { 78261692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 78361692a83SJed Brown } 78461692a83SJed Brown } 78561692a83SJed Brown t->bt[i] = 0; 78661692a83SJed Brown for (j=0; j<s; j++) { 78761692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 78861692a83SJed Brown } 789fe7e6d57SJed Brown if (bembed) { 790fe7e6d57SJed Brown t->bembedt[i] = 0; 791fe7e6d57SJed Brown for (j=0; j<s; j++) { 792fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 793fe7e6d57SJed Brown } 794fe7e6d57SJed Brown } 79561692a83SJed Brown } 7968d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 7978d59e960SJed Brown 798f4aed992SEmil Constantinescu t->pinterp = pinterp; 7993ca35412SEmil Constantinescu ierr = PetscMalloc(s*pinterp*sizeof(binterpt[0]),&t->binterpt);CHKERRQ(ierr); 8003ca35412SEmil Constantinescu ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 80161692a83SJed Brown link->next = RosWTableauList; 80261692a83SJed Brown RosWTableauList = link; 803e27a552bSJed Brown PetscFunctionReturn(0); 804e27a552bSJed Brown } 805e27a552bSJed Brown 806e27a552bSJed Brown #undef __FUNCT__ 80742faf41dSJed Brown #define __FUNCT__ "TSRosWRegisterRos4" 80842faf41dSJed Brown /*@C 80942faf41dSJed Brown TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing paramter choices 81042faf41dSJed Brown 81142faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 81242faf41dSJed Brown 81342faf41dSJed Brown Input Parameters: 81442faf41dSJed Brown + name - identifier for method 81542faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 81642faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 81742faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 81842faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 81942faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 82042faf41dSJed Brown . e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 82142faf41dSJed Brown 82242faf41dSJed Brown Notes: 82342faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 82442faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 82542faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 82642faf41dSJed Brown 82742faf41dSJed Brown Level: developer 82842faf41dSJed Brown 82942faf41dSJed Brown .keywords: TS, register 83042faf41dSJed Brown 83142faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 83242faf41dSJed Brown @*/ 83342faf41dSJed Brown PetscErrorCode TSRosWRegisterRos4(const TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 83442faf41dSJed Brown { 83542faf41dSJed Brown PetscErrorCode ierr; 83642faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 83742faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 83842faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 83942faf41dSJed Brown p42 = one/eight - gamma/three, 84042faf41dSJed Brown p43 = one/twelve - gamma/three, 84142faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 84242faf41dSJed Brown p56 = one/twenty - gamma/four; 84342faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 84442faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 84542faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 84642faf41dSJed Brown 84742faf41dSJed Brown PetscFunctionBegin; 84842faf41dSJed Brown /* Step 1: choose Gamma (input) */ 84942faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 85042faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 85142faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 85242faf41dSJed Brown 85342faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 85442faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 85542faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 85642faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 85742faf41dSJed Brown rhs[0] = one - b3; 85842faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 85942faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 86042faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 86142faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 86242faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 86342faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 86442faf41dSJed Brown 86542faf41dSJed Brown /* Step 3 */ 86642faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 86742faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 86842faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 86942faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 87042faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 87142faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 87242faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 87342faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 87442faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 87542faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 87642faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 87742faf41dSJed Brown 87842faf41dSJed Brown /* Step 4: back-substitute */ 87942faf41dSJed Brown beta32 = beta32beta2p / beta2p; 88042faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 88142faf41dSJed Brown 88242faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 88342faf41dSJed Brown a43 = 0; 88442faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 88542faf41dSJed Brown a42 = a32; 88642faf41dSJed Brown 88742faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 88842faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 88942faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 89042faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 89142faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 89242faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 89342faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 89442faf41dSJed 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; 89542faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 89642faf41dSJed Brown 89742faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 89842faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 89942faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 90042faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 90142faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 90242faf41dSJed Brown 90342faf41dSJed Brown { 90442faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 90542faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 90642faf41dSJed Brown } 90742faf41dSJed Brown ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,PETSC_NULL);CHKERRQ(ierr); 90842faf41dSJed Brown PetscFunctionReturn(0); 90942faf41dSJed Brown } 91042faf41dSJed Brown 91142faf41dSJed Brown #undef __FUNCT__ 9121c3436cfSJed Brown #define __FUNCT__ "TSEvaluateStep_RosW" 9131c3436cfSJed Brown /* 9141c3436cfSJed Brown The step completion formula is 9151c3436cfSJed Brown 9161c3436cfSJed Brown x1 = x0 + b^T Y 9171c3436cfSJed Brown 9181c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 9191c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 9201c3436cfSJed Brown 9211c3436cfSJed Brown x1e = x0 + be^T Y 9221c3436cfSJed Brown = x1 - b^T Y + be^T Y 9231c3436cfSJed Brown = x1 + (be - b)^T Y 9241c3436cfSJed Brown 9251c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9261c3436cfSJed Brown */ 9271c3436cfSJed Brown static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec X,PetscBool *done) 9281c3436cfSJed Brown { 9291c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9301c3436cfSJed Brown RosWTableau tab = ros->tableau; 9311c3436cfSJed Brown PetscScalar *w = ros->work; 9321c3436cfSJed Brown PetscInt i; 9331c3436cfSJed Brown PetscErrorCode ierr; 9341c3436cfSJed Brown 9351c3436cfSJed Brown PetscFunctionBegin; 9361c3436cfSJed Brown if (order == tab->order) { 937108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 9381c3436cfSJed Brown ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 939de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 940de19f811SJed Brown ierr = VecMAXPY(X,tab->s,w,ros->Y);CHKERRQ(ierr); 941108c343cSJed Brown } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 9421c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9431c3436cfSJed Brown PetscFunctionReturn(0); 9441c3436cfSJed Brown } else if (order == tab->order-1) { 9451c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 946108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 9471c3436cfSJed Brown ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 948de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 949de19f811SJed Brown ierr = VecMAXPY(X,tab->s,w,ros->Y);CHKERRQ(ierr); 950108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 951108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 952108c343cSJed Brown ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 953108c343cSJed Brown ierr = VecMAXPY(X,tab->s,w,ros->Y);CHKERRQ(ierr); 9541c3436cfSJed Brown } 9551c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9561c3436cfSJed Brown PetscFunctionReturn(0); 9571c3436cfSJed Brown } 9581c3436cfSJed Brown unavailable: 9591c3436cfSJed Brown if (done) *done = PETSC_FALSE; 9601c3436cfSJed 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); 9611c3436cfSJed Brown PetscFunctionReturn(0); 9621c3436cfSJed Brown } 9631c3436cfSJed Brown 9641c3436cfSJed Brown #undef __FUNCT__ 965e27a552bSJed Brown #define __FUNCT__ "TSStep_RosW" 966e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 967e27a552bSJed Brown { 96861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 96961692a83SJed Brown RosWTableau tab = ros->tableau; 970e27a552bSJed Brown const PetscInt s = tab->s; 9711c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9720feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 973c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 97461692a83SJed Brown PetscScalar *w = ros->work; 9757d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 976e27a552bSJed Brown SNES snes; 9771c3436cfSJed Brown TSAdapt adapt; 9781c3436cfSJed Brown PetscInt i,j,its,lits,reject,next_scheme; 979cdbf8f93SLisandro Dalcin PetscReal next_time_step; 9801c3436cfSJed Brown PetscBool accept; 981e27a552bSJed Brown PetscErrorCode ierr; 9820feba352SEmil Constantinescu MatStructure str; 983e27a552bSJed Brown 984e27a552bSJed Brown PetscFunctionBegin; 985e27a552bSJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 986cdbf8f93SLisandro Dalcin next_time_step = ts->time_step; 9871c3436cfSJed Brown accept = PETSC_TRUE; 988108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 989e27a552bSJed Brown 99097335746SJed Brown for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 9911c3436cfSJed Brown const PetscReal h = ts->time_step; 9923ca35412SEmil Constantinescu ierr = VecCopy(ts->vec_sol,ros->VecSolPrev);CHKERRQ(ierr);/*move this at the end*/ 993e27a552bSJed Brown for (i=0; i<s; i++) { 9941c3436cfSJed Brown ros->stage_time = ts->ptime + h*ASum[i]; 995c17803e7SJed Brown if (GammaZeroDiag[i]) { 996c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 997fd96d5b0SEmil Constantinescu ros->shift = 1./h; 998c17803e7SJed Brown } else { 999c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 100061692a83SJed Brown ros->shift = 1./(h*Gamma[i*s+i]); 1001fd96d5b0SEmil Constantinescu } 100261692a83SJed Brown 100361692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 1004de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 1005de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 100661692a83SJed Brown 100761692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 100861692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 100961692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 101061692a83SJed Brown 1011e27a552bSJed Brown /* Initial guess taken from last stage */ 101261692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 101361692a83SJed Brown 10147d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 101561692a83SJed Brown if (!ros->recompute_jacobian && !i) { 101661692a83SJed Brown ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 101761692a83SJed Brown } 101861692a83SJed Brown ierr = SNESSolve(snes,PETSC_NULL,Y[i]);CHKERRQ(ierr); 1019e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1020e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10215ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 102297335746SJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 102397335746SJed Brown ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 102497335746SJed Brown if (!accept) goto reject_step; 10257d4bf2deSEmil Constantinescu } else { 10261ce71dffSSatish Balay Mat J,Jp; 10270feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10280feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 10290feba352SEmil Constantinescu ierr = VecScale(Y[i],-1.0); 10300feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/ 10310feba352SEmil Constantinescu 10320feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10330feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10340feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 10350feba352SEmil Constantinescu /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10360feba352SEmil Constantinescu str = SAME_NONZERO_PATTERN; 1037ccbc64bcSJed Brown ierr = TSGetIJacobian(ts,&J,&Jp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 10380feba352SEmil Constantinescu ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);CHKERRQ(ierr); 10390feba352SEmil Constantinescu ierr = MatMult(J,Zstage,Zdot); 10400feba352SEmil Constantinescu 10410feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10420feba352SEmil Constantinescu ierr = VecScale(Y[i],h); 10435ef26d82SJed Brown ts->ksp_its += 1; 10447d4bf2deSEmil Constantinescu } 1045e27a552bSJed Brown } 10461c3436cfSJed Brown ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 1047108c343cSJed Brown ros->status = TS_STEP_PENDING; 1048e27a552bSJed Brown 10491c3436cfSJed Brown /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 10501c3436cfSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 10511c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10528d59e960SJed Brown ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 10531c3436cfSJed Brown ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 10541c3436cfSJed Brown if (accept) { 10551c3436cfSJed Brown /* ignore next_scheme for now */ 1056e27a552bSJed Brown ts->ptime += ts->time_step; 1057cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 1058e27a552bSJed Brown ts->steps++; 1059108c343cSJed Brown ros->status = TS_STEP_COMPLETE; 10601c3436cfSJed Brown break; 10611c3436cfSJed Brown } else { /* Roll back the current step */ 10621c3436cfSJed Brown for (i=0; i<s; i++) w[i] = -tab->bt[i]; 10631c3436cfSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 10641c3436cfSJed Brown ts->time_step = next_time_step; 1065108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 10661c3436cfSJed Brown } 1067476b6736SJed Brown reject_step: continue; 10681c3436cfSJed Brown } 1069b2ce242eSJed Brown if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 1070e27a552bSJed Brown PetscFunctionReturn(0); 1071e27a552bSJed Brown } 1072e27a552bSJed Brown 1073e27a552bSJed Brown #undef __FUNCT__ 1074e27a552bSJed Brown #define __FUNCT__ "TSInterpolate_RosW" 1075e27a552bSJed Brown static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec X) 1076e27a552bSJed Brown { 107761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1078f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1079f4aed992SEmil Constantinescu PetscReal h; 1080f4aed992SEmil Constantinescu PetscReal tt,t; 1081f4aed992SEmil Constantinescu PetscScalar *bt; 1082f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1083f4aed992SEmil Constantinescu PetscErrorCode ierr; 1084f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1085f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1086f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1087e27a552bSJed Brown 1088e27a552bSJed Brown PetscFunctionBegin; 1089f4aed992SEmil Constantinescu if (!Bt) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1090f4aed992SEmil Constantinescu 1091f4aed992SEmil Constantinescu switch (ros->status) { 1092f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1093f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1094f4aed992SEmil Constantinescu h = ts->time_step; 1095f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1096f4aed992SEmil Constantinescu break; 1097f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1098f4aed992SEmil Constantinescu h = ts->time_step_prev; 1099f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1100f4aed992SEmil Constantinescu break; 1101f4aed992SEmil Constantinescu default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1102f4aed992SEmil Constantinescu } 11033ca35412SEmil Constantinescu ierr = PetscMalloc(s*sizeof(bt[0]),&bt);CHKERRQ(ierr); 1104f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1105f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1106f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 11073ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1108f4aed992SEmil Constantinescu } 1109f4aed992SEmil Constantinescu } 1110f4aed992SEmil Constantinescu 1111f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1112f4aed992SEmil Constantinescu /*X<-0*/ 1113f4aed992SEmil Constantinescu ierr = VecZeroEntries(X);CHKERRQ(ierr); 1114f4aed992SEmil Constantinescu 1115f4aed992SEmil Constantinescu /*X<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11163ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j]=0; 11173ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11183ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11193ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1120f4aed992SEmil Constantinescu } 11213ca35412SEmil Constantinescu } 11223ca35412SEmil Constantinescu ierr = VecMAXPY(X,i,w,Y);CHKERRQ(ierr); 1123f4aed992SEmil Constantinescu 1124f4aed992SEmil Constantinescu /*X<-y(t) + X*/ 11253ca35412SEmil Constantinescu ierr = VecAXPY(X,1.0,ros->VecSolPrev);CHKERRQ(ierr); 1126f4aed992SEmil Constantinescu 1127f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1128f4aed992SEmil Constantinescu 1129e27a552bSJed Brown PetscFunctionReturn(0); 1130e27a552bSJed Brown } 1131e27a552bSJed Brown 1132e27a552bSJed Brown /*------------------------------------------------------------*/ 1133e27a552bSJed Brown #undef __FUNCT__ 1134e27a552bSJed Brown #define __FUNCT__ "TSReset_RosW" 1135e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1136e27a552bSJed Brown { 113761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1138e27a552bSJed Brown PetscInt s; 1139e27a552bSJed Brown PetscErrorCode ierr; 1140e27a552bSJed Brown 1141e27a552bSJed Brown PetscFunctionBegin; 114261692a83SJed Brown if (!ros->tableau) PetscFunctionReturn(0); 114361692a83SJed Brown s = ros->tableau->s; 114461692a83SJed Brown ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 114561692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 114661692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 114761692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 114861692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 11493ca35412SEmil Constantinescu ierr = VecDestroy(&ros->VecSolPrev);CHKERRQ(ierr); 115061692a83SJed Brown ierr = PetscFree(ros->work);CHKERRQ(ierr); 1151e27a552bSJed Brown PetscFunctionReturn(0); 1152e27a552bSJed Brown } 1153e27a552bSJed Brown 1154e27a552bSJed Brown #undef __FUNCT__ 1155e27a552bSJed Brown #define __FUNCT__ "TSDestroy_RosW" 1156e27a552bSJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1157e27a552bSJed Brown { 1158e27a552bSJed Brown PetscErrorCode ierr; 1159e27a552bSJed Brown 1160e27a552bSJed Brown PetscFunctionBegin; 1161e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1162e27a552bSJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 1163e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 1164e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);CHKERRQ(ierr); 116561692a83SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",PETSC_NULL);CHKERRQ(ierr); 1166e27a552bSJed Brown PetscFunctionReturn(0); 1167e27a552bSJed Brown } 1168e27a552bSJed Brown 1169d5e6173cSPeter Brune 1170d5e6173cSPeter Brune #undef __FUNCT__ 1171d5e6173cSPeter Brune #define __FUNCT__ "TSRosWGetVecs" 1172d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1173d5e6173cSPeter Brune { 1174d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1175d5e6173cSPeter Brune PetscErrorCode ierr; 1176d5e6173cSPeter Brune 1177d5e6173cSPeter Brune PetscFunctionBegin; 1178d5e6173cSPeter Brune if (Ydot) { 1179d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1180d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1181d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1182d5e6173cSPeter Brune } 1183d5e6173cSPeter Brune if (Zdot) { 1184d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1185d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1186d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1187d5e6173cSPeter Brune } 1188d5e6173cSPeter Brune if (Ystage) { 1189d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1190d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1191d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1192d5e6173cSPeter Brune } 1193d5e6173cSPeter Brune if (Zstage) { 1194d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1195d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1196d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1197d5e6173cSPeter Brune } 1198d5e6173cSPeter Brune 1199d5e6173cSPeter Brune PetscFunctionReturn(0); 1200d5e6173cSPeter Brune } 1201d5e6173cSPeter Brune 1202d5e6173cSPeter Brune 1203d5e6173cSPeter Brune #undef __FUNCT__ 1204d5e6173cSPeter Brune #define __FUNCT__ "TSRosWRestoreVecs" 1205d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1206d5e6173cSPeter Brune { 1207d5e6173cSPeter Brune PetscErrorCode ierr; 1208d5e6173cSPeter Brune 1209d5e6173cSPeter Brune PetscFunctionBegin; 1210d5e6173cSPeter Brune if (Ydot) { 1211d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1212d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1213d5e6173cSPeter Brune } 1214d5e6173cSPeter Brune } 1215d5e6173cSPeter Brune if (Zdot) { 1216d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1217d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1218d5e6173cSPeter Brune } 1219d5e6173cSPeter Brune } 1220d5e6173cSPeter Brune if (Ystage) { 1221d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1222d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1223d5e6173cSPeter Brune } 1224d5e6173cSPeter Brune } 1225d5e6173cSPeter Brune if (Zstage) { 1226d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1227d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1228d5e6173cSPeter Brune } 1229d5e6173cSPeter Brune } 1230d5e6173cSPeter Brune PetscFunctionReturn(0); 1231d5e6173cSPeter Brune } 1232d5e6173cSPeter Brune 1233d5e6173cSPeter Brune #undef __FUNCT__ 1234d5e6173cSPeter Brune #define __FUNCT__ "DMCoarsenHook_TSRosW" 1235d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1236d5e6173cSPeter Brune { 1237d5e6173cSPeter Brune 1238d5e6173cSPeter Brune PetscFunctionBegin; 1239d5e6173cSPeter Brune PetscFunctionReturn(0); 1240d5e6173cSPeter Brune } 1241d5e6173cSPeter Brune 1242d5e6173cSPeter Brune #undef __FUNCT__ 1243d5e6173cSPeter Brune #define __FUNCT__ "DMRestrictHook_TSRosW" 1244d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1245d5e6173cSPeter Brune { 1246d5e6173cSPeter Brune TS ts = (TS)ctx; 1247d5e6173cSPeter Brune PetscErrorCode ierr; 1248d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1249d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1250d5e6173cSPeter Brune 1251d5e6173cSPeter Brune PetscFunctionBegin; 1252d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1253d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1254d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1255d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1256d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1257d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1258d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1259d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1260d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1261d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1262d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1263d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1264d5e6173cSPeter Brune PetscFunctionReturn(0); 1265d5e6173cSPeter Brune } 1266d5e6173cSPeter Brune 1267e27a552bSJed Brown /* 1268e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1269e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1270e27a552bSJed Brown */ 1271e27a552bSJed Brown #undef __FUNCT__ 1272e27a552bSJed Brown #define __FUNCT__ "SNESTSFormFunction_RosW" 1273e27a552bSJed Brown static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec X,Vec F,TS ts) 1274e27a552bSJed Brown { 127561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1276e27a552bSJed Brown PetscErrorCode ierr; 1277d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1278d5e6173cSPeter Brune DM dm,dmsave; 1279e27a552bSJed Brown 1280e27a552bSJed Brown PetscFunctionBegin; 1281d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1282d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1283d5e6173cSPeter Brune ierr = VecWAXPY(Ydot,ros->shift,X,Zdot);CHKERRQ(ierr); /* Ydot = shift*X + Zdot */ 1284d5e6173cSPeter Brune ierr = VecWAXPY(Ystage,1.0,X,Zstage);CHKERRQ(ierr); /* Ystage = X + Zstage */ 1285d5e6173cSPeter Brune dmsave = ts->dm; 1286d5e6173cSPeter Brune ts->dm = dm; 1287d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1288d5e6173cSPeter Brune ts->dm = dmsave; 1289d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1290e27a552bSJed Brown PetscFunctionReturn(0); 1291e27a552bSJed Brown } 1292e27a552bSJed Brown 1293e27a552bSJed Brown #undef __FUNCT__ 1294e27a552bSJed Brown #define __FUNCT__ "SNESTSFormJacobian_RosW" 1295e27a552bSJed Brown static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 1296e27a552bSJed Brown { 129761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1298d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1299e27a552bSJed Brown PetscErrorCode ierr; 1300d5e6173cSPeter Brune DM dm,dmsave; 1301e27a552bSJed Brown 1302e27a552bSJed Brown PetscFunctionBegin; 130361692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1304d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1305d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1306d5e6173cSPeter Brune dmsave = ts->dm; 1307d5e6173cSPeter Brune ts->dm = dm; 1308d5e6173cSPeter Brune ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,ros->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 1309d5e6173cSPeter Brune ts->dm = dmsave; 1310d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1311e27a552bSJed Brown PetscFunctionReturn(0); 1312e27a552bSJed Brown } 1313e27a552bSJed Brown 1314e27a552bSJed Brown #undef __FUNCT__ 1315e27a552bSJed Brown #define __FUNCT__ "TSSetUp_RosW" 1316e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1317e27a552bSJed Brown { 131861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 131961692a83SJed Brown RosWTableau tab = ros->tableau; 1320e27a552bSJed Brown PetscInt s = tab->s; 1321e27a552bSJed Brown PetscErrorCode ierr; 1322d5e6173cSPeter Brune DM dm; 1323e27a552bSJed Brown 1324e27a552bSJed Brown PetscFunctionBegin; 132561692a83SJed Brown if (!ros->tableau) { 1326e27a552bSJed Brown ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1327e27a552bSJed Brown } 132861692a83SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 132961692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 133061692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 133161692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 133261692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 13333ca35412SEmil Constantinescu ierr = VecDuplicate(ts->vec_sol,&ros->VecSolPrev);CHKERRQ(ierr); 133461692a83SJed Brown ierr = PetscMalloc(s*sizeof(ros->work[0]),&ros->work);CHKERRQ(ierr); 1335d5e6173cSPeter Brune ierr = TSGetDM(ts,&dm); 1336d5e6173cSPeter Brune if (dm) { 1337d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1338d5e6173cSPeter Brune } 1339e27a552bSJed Brown PetscFunctionReturn(0); 1340e27a552bSJed Brown } 1341e27a552bSJed Brown /*------------------------------------------------------------*/ 1342e27a552bSJed Brown 1343e27a552bSJed Brown #undef __FUNCT__ 1344e27a552bSJed Brown #define __FUNCT__ "TSSetFromOptions_RosW" 1345e27a552bSJed Brown static PetscErrorCode TSSetFromOptions_RosW(TS ts) 1346e27a552bSJed Brown { 134761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1348e27a552bSJed Brown PetscErrorCode ierr; 134961692a83SJed Brown char rostype[256]; 1350e27a552bSJed Brown 1351e27a552bSJed Brown PetscFunctionBegin; 1352e27a552bSJed Brown ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 1353e27a552bSJed Brown { 135461692a83SJed Brown RosWTableauLink link; 1355e27a552bSJed Brown PetscInt count,choice; 1356e27a552bSJed Brown PetscBool flg; 1357e27a552bSJed Brown const char **namelist; 135861692a83SJed Brown SNES snes; 135961692a83SJed Brown 136061692a83SJed Brown ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof rostype);CHKERRQ(ierr); 136161692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1362e27a552bSJed Brown ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 136361692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 136461692a83SJed Brown ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 136561692a83SJed Brown ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 1366e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 136761692a83SJed Brown 136861692a83SJed Brown ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,PETSC_NULL);CHKERRQ(ierr); 136961692a83SJed Brown 137061692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 137161692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 137261692a83SJed Brown if (!((PetscObject)snes)->type_name) { 137361692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 137461692a83SJed Brown } 137561692a83SJed Brown ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1376e27a552bSJed Brown } 1377e27a552bSJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 1378e27a552bSJed Brown PetscFunctionReturn(0); 1379e27a552bSJed Brown } 1380e27a552bSJed Brown 1381e27a552bSJed Brown #undef __FUNCT__ 1382e27a552bSJed Brown #define __FUNCT__ "PetscFormatRealArray" 1383e27a552bSJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1384e27a552bSJed Brown { 1385e27a552bSJed Brown PetscErrorCode ierr; 1386e408995aSJed Brown PetscInt i; 1387e408995aSJed Brown size_t left,count; 1388e27a552bSJed Brown char *p; 1389e27a552bSJed Brown 1390e27a552bSJed Brown PetscFunctionBegin; 1391e408995aSJed Brown for (i=0,p=buf,left=len; i<n; i++) { 1392e408995aSJed Brown ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1393e27a552bSJed Brown if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1394e27a552bSJed Brown left -= count; 1395e27a552bSJed Brown p += count; 1396e27a552bSJed Brown *p++ = ' '; 1397e27a552bSJed Brown } 1398e27a552bSJed Brown p[i ? 0 : -1] = 0; 1399e27a552bSJed Brown PetscFunctionReturn(0); 1400e27a552bSJed Brown } 1401e27a552bSJed Brown 1402e27a552bSJed Brown #undef __FUNCT__ 1403e27a552bSJed Brown #define __FUNCT__ "TSView_RosW" 1404e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1405e27a552bSJed Brown { 140661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 140761692a83SJed Brown RosWTableau tab = ros->tableau; 1408e27a552bSJed Brown PetscBool iascii; 1409e27a552bSJed Brown PetscErrorCode ierr; 1410e27a552bSJed Brown 1411e27a552bSJed Brown PetscFunctionBegin; 1412251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1413e27a552bSJed Brown if (iascii) { 141461692a83SJed Brown const TSRosWType rostype; 1415e408995aSJed Brown PetscInt i; 1416e408995aSJed Brown PetscReal abscissa[512]; 1417e27a552bSJed Brown char buf[512]; 141861692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 141961692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1420e408995aSJed Brown ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 142161692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1422e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1423e408995aSJed Brown ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1424e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1425e27a552bSJed Brown } 1426e27a552bSJed Brown ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1427e27a552bSJed Brown PetscFunctionReturn(0); 1428e27a552bSJed Brown } 1429e27a552bSJed Brown 1430e27a552bSJed Brown #undef __FUNCT__ 1431e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType" 1432e27a552bSJed Brown /*@C 143361692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1434e27a552bSJed Brown 1435e27a552bSJed Brown Logically collective 1436e27a552bSJed Brown 1437e27a552bSJed Brown Input Parameter: 1438e27a552bSJed Brown + ts - timestepping context 143961692a83SJed Brown - rostype - type of Rosenbrock-W scheme 1440e27a552bSJed Brown 1441020d8f30SJed Brown Level: beginner 1442e27a552bSJed Brown 1443020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1444e27a552bSJed Brown @*/ 144561692a83SJed Brown PetscErrorCode TSRosWSetType(TS ts,const TSRosWType rostype) 1446e27a552bSJed Brown { 1447e27a552bSJed Brown PetscErrorCode ierr; 1448e27a552bSJed Brown 1449e27a552bSJed Brown PetscFunctionBegin; 1450e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 145161692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,const TSRosWType),(ts,rostype));CHKERRQ(ierr); 1452e27a552bSJed Brown PetscFunctionReturn(0); 1453e27a552bSJed Brown } 1454e27a552bSJed Brown 1455e27a552bSJed Brown #undef __FUNCT__ 1456e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType" 1457e27a552bSJed Brown /*@C 145861692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1459e27a552bSJed Brown 1460e27a552bSJed Brown Logically collective 1461e27a552bSJed Brown 1462e27a552bSJed Brown Input Parameter: 1463e27a552bSJed Brown . ts - timestepping context 1464e27a552bSJed Brown 1465e27a552bSJed Brown Output Parameter: 146661692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1467e27a552bSJed Brown 1468e27a552bSJed Brown Level: intermediate 1469e27a552bSJed Brown 1470e27a552bSJed Brown .seealso: TSRosWGetType() 1471e27a552bSJed Brown @*/ 147261692a83SJed Brown PetscErrorCode TSRosWGetType(TS ts,const TSRosWType *rostype) 1473e27a552bSJed Brown { 1474e27a552bSJed Brown PetscErrorCode ierr; 1475e27a552bSJed Brown 1476e27a552bSJed Brown PetscFunctionBegin; 1477e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 147861692a83SJed Brown ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,const TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1479e27a552bSJed Brown PetscFunctionReturn(0); 1480e27a552bSJed Brown } 1481e27a552bSJed Brown 1482e27a552bSJed Brown #undef __FUNCT__ 148361692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian" 1484e27a552bSJed Brown /*@C 148561692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1486e27a552bSJed Brown 1487e27a552bSJed Brown Logically collective 1488e27a552bSJed Brown 1489e27a552bSJed Brown Input Parameter: 1490e27a552bSJed Brown + ts - timestepping context 149161692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1492e27a552bSJed Brown 1493e27a552bSJed Brown Level: intermediate 1494e27a552bSJed Brown 1495e27a552bSJed Brown .seealso: TSRosWGetType() 1496e27a552bSJed Brown @*/ 149761692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1498e27a552bSJed Brown { 1499e27a552bSJed Brown PetscErrorCode ierr; 1500e27a552bSJed Brown 1501e27a552bSJed Brown PetscFunctionBegin; 1502e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 150361692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1504e27a552bSJed Brown PetscFunctionReturn(0); 1505e27a552bSJed Brown } 1506e27a552bSJed Brown 1507e27a552bSJed Brown EXTERN_C_BEGIN 1508e27a552bSJed Brown #undef __FUNCT__ 1509e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType_RosW" 151061692a83SJed Brown PetscErrorCode TSRosWGetType_RosW(TS ts,const TSRosWType *rostype) 1511e27a552bSJed Brown { 151261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1513e27a552bSJed Brown PetscErrorCode ierr; 1514e27a552bSJed Brown 1515e27a552bSJed Brown PetscFunctionBegin; 151661692a83SJed Brown if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 151761692a83SJed Brown *rostype = ros->tableau->name; 1518e27a552bSJed Brown PetscFunctionReturn(0); 1519e27a552bSJed Brown } 1520e27a552bSJed Brown #undef __FUNCT__ 1521e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType_RosW" 152261692a83SJed Brown PetscErrorCode TSRosWSetType_RosW(TS ts,const TSRosWType rostype) 1523e27a552bSJed Brown { 152461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1525e27a552bSJed Brown PetscErrorCode ierr; 1526e27a552bSJed Brown PetscBool match; 152761692a83SJed Brown RosWTableauLink link; 1528e27a552bSJed Brown 1529e27a552bSJed Brown PetscFunctionBegin; 153061692a83SJed Brown if (ros->tableau) { 153161692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1532e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1533e27a552bSJed Brown } 153461692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 153561692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1536e27a552bSJed Brown if (match) { 1537e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 153861692a83SJed Brown ros->tableau = &link->tab; 1539e27a552bSJed Brown PetscFunctionReturn(0); 1540e27a552bSJed Brown } 1541e27a552bSJed Brown } 154261692a83SJed Brown SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1543e27a552bSJed Brown PetscFunctionReturn(0); 1544e27a552bSJed Brown } 154561692a83SJed Brown 1546e27a552bSJed Brown #undef __FUNCT__ 154761692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 154861692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1549e27a552bSJed Brown { 155061692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1551e27a552bSJed Brown 1552e27a552bSJed Brown PetscFunctionBegin; 155361692a83SJed Brown ros->recompute_jacobian = flg; 1554e27a552bSJed Brown PetscFunctionReturn(0); 1555e27a552bSJed Brown } 1556e27a552bSJed Brown EXTERN_C_END 1557e27a552bSJed Brown 1558d5e6173cSPeter Brune 1559e27a552bSJed Brown /* ------------------------------------------------------------ */ 1560e27a552bSJed Brown /*MC 1561020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1562e27a552bSJed Brown 1563e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1564e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1565e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1566e27a552bSJed Brown 1567e27a552bSJed Brown Notes: 156861692a83SJed Brown This method currently only works with autonomous ODE and DAE. 156961692a83SJed Brown 157061692a83SJed Brown Developer notes: 157161692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 157261692a83SJed Brown 157361692a83SJed Brown $ xdot = f(x) 157461692a83SJed Brown 157561692a83SJed Brown by the stage equations 157661692a83SJed Brown 157761692a83SJed Brown $ k_i = h f(x_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 157861692a83SJed Brown 157961692a83SJed Brown and step completion formula 158061692a83SJed Brown 158161692a83SJed Brown $ x_1 = x_0 + sum_j b_j k_j 158261692a83SJed Brown 158361692a83SJed Brown with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(x) 158461692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 158561692a83SJed Brown we define new variables for the stage equations 158661692a83SJed Brown 158761692a83SJed Brown $ y_i = gamma_ij k_j 158861692a83SJed Brown 158961692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 159061692a83SJed Brown 159161692a83SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 159261692a83SJed Brown 159361692a83SJed Brown to rewrite the method as 159461692a83SJed Brown 159561692a83SJed Brown $ [M/(h gamma_ii) - J] y_i = f(x_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 159661692a83SJed Brown $ x_1 = x_0 + sum_j bt_j y_j 159761692a83SJed Brown 159861692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 159961692a83SJed Brown 160061692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 160161692a83SJed Brown 160261692a83SJed Brown or, more compactly in tensor notation 160361692a83SJed Brown 160461692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 160561692a83SJed Brown 160661692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 160761692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 160861692a83SJed Brown equation 160961692a83SJed Brown 161061692a83SJed Brown $ g(x_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 161161692a83SJed Brown 161261692a83SJed Brown with initial guess y_i = 0. 1613e27a552bSJed Brown 1614e27a552bSJed Brown Level: beginner 1615e27a552bSJed Brown 1616*a4386c9eSJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1617*a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1618e27a552bSJed Brown M*/ 1619e27a552bSJed Brown EXTERN_C_BEGIN 1620e27a552bSJed Brown #undef __FUNCT__ 1621e27a552bSJed Brown #define __FUNCT__ "TSCreate_RosW" 1622e27a552bSJed Brown PetscErrorCode TSCreate_RosW(TS ts) 1623e27a552bSJed Brown { 162461692a83SJed Brown TS_RosW *ros; 1625e27a552bSJed Brown PetscErrorCode ierr; 1626e27a552bSJed Brown 1627e27a552bSJed Brown PetscFunctionBegin; 1628e27a552bSJed Brown #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1629e27a552bSJed Brown ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1630e27a552bSJed Brown #endif 1631e27a552bSJed Brown 1632e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1633e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1634e27a552bSJed Brown ts->ops->view = TSView_RosW; 1635e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1636e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1637e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16381c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 1639e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1640e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1641e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1642e27a552bSJed Brown 164361692a83SJed Brown ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 164461692a83SJed Brown ts->data = (void*)ros; 1645e27a552bSJed Brown 1646e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 1647e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 164861692a83SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1649e27a552bSJed Brown PetscFunctionReturn(0); 1650e27a552bSJed Brown } 1651e27a552bSJed Brown EXTERN_C_END 1652