1e27a552bSJed Brown /* 261692a83SJed Brown Code for timestepping with Rosenbrock W methods 3e27a552bSJed Brown 4e27a552bSJed Brown Notes: 5e27a552bSJed Brown The general system is written as 6e27a552bSJed Brown 7f9c1d6abSBarry Smith F(t,U,Udot) = G(t,U) 8e27a552bSJed Brown 9f9c1d6abSBarry Smith where F represents the stiff part of the physics and G represents the non-stiff part. 10f9c1d6abSBarry Smith This method is designed to be linearly implicit on F and can use an approximate and lagged Jacobian. 11e27a552bSJed Brown 12e27a552bSJed Brown */ 13b45d2f2cSJed Brown #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 14e27a552bSJed Brown 1561692a83SJed Brown #include <../src/mat/blockinvert.h> 1661692a83SJed Brown 1719fd82e9SBarry Smith static TSRosWType TSRosWDefault = TSROSWRA34PW2; 18e27a552bSJed Brown static PetscBool TSRosWRegisterAllCalled; 19e27a552bSJed Brown static PetscBool TSRosWPackageInitialized; 20e27a552bSJed Brown 2161692a83SJed Brown typedef struct _RosWTableau *RosWTableau; 2261692a83SJed Brown struct _RosWTableau { 23e27a552bSJed Brown char *name; 24e27a552bSJed Brown PetscInt order; /* Classical approximation order of the method */ 25e27a552bSJed Brown PetscInt s; /* Number of stages */ 26f4aed992SEmil Constantinescu PetscInt pinterp; /* Interpolation order */ 2761692a83SJed Brown PetscReal *A; /* Propagation table, strictly lower triangular */ 2861692a83SJed Brown PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */ 29c17803e7SJed Brown PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */ 3043b21953SEmil Constantinescu PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/ 3161692a83SJed Brown PetscReal *b; /* Step completion table */ 32fe7e6d57SJed Brown PetscReal *bembed; /* Step completion table for embedded method of order one less */ 3361692a83SJed Brown PetscReal *ASum; /* Row sum of A */ 3461692a83SJed Brown PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */ 3561692a83SJed Brown PetscReal *At; /* Propagation table in transformed variables */ 3661692a83SJed Brown PetscReal *bt; /* Step completion table in transformed variables */ 37fe7e6d57SJed Brown PetscReal *bembedt; /* Step completion table of order one less in transformed variables */ 3861692a83SJed Brown PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */ 398d59e960SJed Brown PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 403ca35412SEmil Constantinescu PetscReal *binterpt; /* Dense output formula */ 41e27a552bSJed Brown }; 4261692a83SJed Brown typedef struct _RosWTableauLink *RosWTableauLink; 4361692a83SJed Brown struct _RosWTableauLink { 4461692a83SJed Brown struct _RosWTableau tab; 4561692a83SJed Brown RosWTableauLink next; 46e27a552bSJed Brown }; 4761692a83SJed Brown static RosWTableauLink RosWTableauList; 48e27a552bSJed Brown 49e27a552bSJed Brown typedef struct { 5061692a83SJed Brown RosWTableau tableau; 5161692a83SJed Brown Vec *Y; /* States computed during the step, used to complete the step */ 52e27a552bSJed Brown Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 5361692a83SJed Brown Vec Ystage; /* Work vector for the state value at each stage */ 5461692a83SJed Brown Vec Zdot; /* Ydot = Zdot + shift*Y */ 5561692a83SJed Brown Vec Zstage; /* Y = Zstage + Y */ 563ca35412SEmil Constantinescu Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation)*/ 571c3436cfSJed Brown PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 58*b296d7d5SJed Brown PetscReal scoeff; /* shift = scoeff/dt */ 59e27a552bSJed Brown PetscReal stage_time; 60c17803e7SJed Brown PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */ 6161692a83SJed Brown PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 62108c343cSJed Brown TSStepStatus status; 63e27a552bSJed Brown } TS_RosW; 64e27a552bSJed Brown 65fe7e6d57SJed Brown /*MC 663606a31eSEmil Constantinescu TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method). 673606a31eSEmil Constantinescu 683606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 693606a31eSEmil Constantinescu 703606a31eSEmil Constantinescu Level: intermediate 713606a31eSEmil Constantinescu 723606a31eSEmil Constantinescu .seealso: TSROSW 733606a31eSEmil Constantinescu M*/ 743606a31eSEmil Constantinescu 753606a31eSEmil Constantinescu /*MC 763606a31eSEmil Constantinescu TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method). 773606a31eSEmil Constantinescu 783606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 793606a31eSEmil Constantinescu 803606a31eSEmil Constantinescu Level: intermediate 813606a31eSEmil Constantinescu 823606a31eSEmil Constantinescu .seealso: TSROSW 833606a31eSEmil Constantinescu M*/ 843606a31eSEmil Constantinescu 853606a31eSEmil Constantinescu /*MC 86fe7e6d57SJed Brown TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme. 87fe7e6d57SJed Brown 88fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P. 89fe7e6d57SJed Brown 90fe7e6d57SJed Brown Level: intermediate 91fe7e6d57SJed Brown 92fe7e6d57SJed Brown .seealso: TSROSW 93fe7e6d57SJed Brown M*/ 94fe7e6d57SJed Brown 95fe7e6d57SJed Brown /*MC 96fe7e6d57SJed Brown TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme. 97fe7e6d57SJed Brown 98fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M. 99fe7e6d57SJed Brown 100fe7e6d57SJed Brown Level: intermediate 101fe7e6d57SJed Brown 102fe7e6d57SJed Brown .seealso: TSROSW 103fe7e6d57SJed Brown M*/ 104fe7e6d57SJed Brown 105fe7e6d57SJed Brown /*MC 106fe7e6d57SJed Brown TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1. 107fe7e6d57SJed Brown 108fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 109fe7e6d57SJed Brown 110fe7e6d57SJed Brown This is strongly A-stable with R(infty) = 0.73. The embedded method of order 2 is strongly A-stable with R(infty) = 0.73. 111fe7e6d57SJed Brown 112fe7e6d57SJed Brown References: 113fe7e6d57SJed Brown Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 114fe7e6d57SJed Brown 115fe7e6d57SJed Brown Level: intermediate 116fe7e6d57SJed Brown 117fe7e6d57SJed Brown .seealso: TSROSW 118fe7e6d57SJed Brown M*/ 119fe7e6d57SJed Brown 120fe7e6d57SJed Brown /*MC 121fe7e6d57SJed Brown TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1. 122fe7e6d57SJed Brown 123fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 124fe7e6d57SJed Brown 125fe7e6d57SJed Brown This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48. 126fe7e6d57SJed Brown 127fe7e6d57SJed Brown References: 128fe7e6d57SJed Brown Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 129fe7e6d57SJed Brown 130fe7e6d57SJed Brown Level: intermediate 131fe7e6d57SJed Brown 132fe7e6d57SJed Brown .seealso: TSROSW 133fe7e6d57SJed Brown M*/ 134fe7e6d57SJed Brown 135ef3c5b88SJed Brown /*MC 136ef3c5b88SJed Brown TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme 137ef3c5b88SJed Brown 138ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 139ef3c5b88SJed Brown 140ef3c5b88SJed Brown Both the third order and embedded second order methods are stiffly accurate and L-stable. 141ef3c5b88SJed Brown 142ef3c5b88SJed Brown References: 143ef3c5b88SJed Brown Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 144ef3c5b88SJed Brown 145ef3c5b88SJed Brown Level: intermediate 146ef3c5b88SJed Brown 147ef3c5b88SJed Brown .seealso: TSROSW, TSROSWSANDU3 148ef3c5b88SJed Brown M*/ 149ef3c5b88SJed Brown 150ef3c5b88SJed Brown /*MC 151ef3c5b88SJed Brown TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme 152ef3c5b88SJed Brown 153ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 154ef3c5b88SJed Brown 155ef3c5b88SJed Brown The third order method is L-stable, but not stiffly accurate. 156ef3c5b88SJed Brown The second order embedded method is strongly A-stable with R(infty) = 0.5. 157ef3c5b88SJed Brown The internal stages are L-stable. 158ef3c5b88SJed Brown This method is called ROS3 in the paper. 159ef3c5b88SJed Brown 160ef3c5b88SJed Brown References: 161ef3c5b88SJed Brown Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 162ef3c5b88SJed Brown 163ef3c5b88SJed Brown Level: intermediate 164ef3c5b88SJed Brown 165ef3c5b88SJed Brown .seealso: TSROSW, TSROSWRODAS3 166ef3c5b88SJed Brown M*/ 167ef3c5b88SJed Brown 168961f28d0SJed Brown /*MC 169961f28d0SJed Brown TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 170961f28d0SJed Brown 171961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 172961f28d0SJed Brown 173961f28d0SJed Brown A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 174961f28d0SJed Brown 175961f28d0SJed Brown References: 176961f28d0SJed Brown Emil Constantinescu 177961f28d0SJed Brown 178961f28d0SJed Brown Level: intermediate 179961f28d0SJed Brown 18043b21953SEmil Constantinescu .seealso: TSROSW, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, SSP 181961f28d0SJed Brown M*/ 182961f28d0SJed Brown 183961f28d0SJed Brown /*MC 184998eb97aSJed Brown TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 185961f28d0SJed Brown 186961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 187961f28d0SJed Brown 188961f28d0SJed Brown L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 189961f28d0SJed Brown 190961f28d0SJed Brown References: 191961f28d0SJed Brown Emil Constantinescu 192961f28d0SJed Brown 193961f28d0SJed Brown Level: intermediate 194961f28d0SJed Brown 19543b21953SEmil Constantinescu .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLLSSP3P4S2C, TSSSP 196961f28d0SJed Brown M*/ 197961f28d0SJed Brown 198961f28d0SJed Brown /*MC 199998eb97aSJed Brown TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 200961f28d0SJed Brown 201961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 202961f28d0SJed Brown 203961f28d0SJed Brown L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 204961f28d0SJed Brown 205961f28d0SJed Brown References: 206961f28d0SJed Brown Emil Constantinescu 207961f28d0SJed Brown 208961f28d0SJed Brown Level: intermediate 209961f28d0SJed Brown 210961f28d0SJed Brown .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSSSP 211961f28d0SJed Brown M*/ 212961f28d0SJed Brown 21342faf41dSJed Brown /*MC 21442faf41dSJed Brown TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop 21542faf41dSJed Brown 21642faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 21742faf41dSJed Brown 21842faf41dSJed Brown A(89.3 degrees)-stable, |R(infty)| = 0.454. 21942faf41dSJed Brown 22042faf41dSJed Brown This method does not provide a dense output formula. 22142faf41dSJed Brown 22242faf41dSJed Brown References: 22342faf41dSJed Brown Kaps and Rentrop, Generalized Runge-Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979. 22442faf41dSJed Brown 22542faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 22642faf41dSJed Brown 22742faf41dSJed Brown Hairer's code ros4.f 22842faf41dSJed Brown 22942faf41dSJed Brown Level: intermediate 23042faf41dSJed Brown 23142faf41dSJed Brown .seealso: TSROSW, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 23242faf41dSJed Brown M*/ 23342faf41dSJed Brown 23442faf41dSJed Brown /*MC 23542faf41dSJed Brown TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine 23642faf41dSJed Brown 23742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 23842faf41dSJed Brown 23942faf41dSJed Brown A-stable, |R(infty)| = 1/3. 24042faf41dSJed Brown 24142faf41dSJed Brown This method does not provide a dense output formula. 24242faf41dSJed Brown 24342faf41dSJed Brown References: 24442faf41dSJed Brown Shampine, Implementation of Rosenbrock methods, 1982. 24542faf41dSJed Brown 24642faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 24742faf41dSJed Brown 24842faf41dSJed Brown Hairer's code ros4.f 24942faf41dSJed Brown 25042faf41dSJed Brown Level: intermediate 25142faf41dSJed Brown 25242faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWVELDD4, TSROSW4L 25342faf41dSJed Brown M*/ 25442faf41dSJed Brown 25542faf41dSJed Brown /*MC 25642faf41dSJed Brown TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen 25742faf41dSJed Brown 25842faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 25942faf41dSJed Brown 26042faf41dSJed Brown A(89.5 degrees)-stable, |R(infty)| = 0.24. 26142faf41dSJed Brown 26242faf41dSJed Brown This method does not provide a dense output formula. 26342faf41dSJed Brown 26442faf41dSJed Brown References: 26542faf41dSJed Brown van Veldhuizen, D-stability and Kaps-Rentrop methods, 1984. 26642faf41dSJed Brown 26742faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 26842faf41dSJed Brown 26942faf41dSJed Brown Hairer's code ros4.f 27042faf41dSJed Brown 27142faf41dSJed Brown Level: intermediate 27242faf41dSJed Brown 27342faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 27442faf41dSJed Brown M*/ 27542faf41dSJed Brown 27642faf41dSJed Brown /*MC 27742faf41dSJed Brown TSROSW4L - four stage, fourth order Rosenbrock (not W) method 27842faf41dSJed Brown 27942faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 28042faf41dSJed Brown 28142faf41dSJed Brown A-stable and L-stable 28242faf41dSJed Brown 28342faf41dSJed Brown This method does not provide a dense output formula. 28442faf41dSJed Brown 28542faf41dSJed Brown References: 28642faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 28742faf41dSJed Brown 28842faf41dSJed Brown Hairer's code ros4.f 28942faf41dSJed Brown 29042faf41dSJed Brown Level: intermediate 29142faf41dSJed Brown 29242faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 29342faf41dSJed Brown M*/ 29442faf41dSJed Brown 295e27a552bSJed Brown #undef __FUNCT__ 296e27a552bSJed Brown #define __FUNCT__ "TSRosWRegisterAll" 297e27a552bSJed Brown /*@C 298e27a552bSJed Brown TSRosWRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSRosW 299e27a552bSJed Brown 300e27a552bSJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 301e27a552bSJed Brown 302e27a552bSJed Brown Level: advanced 303e27a552bSJed Brown 304e27a552bSJed Brown .keywords: TS, TSRosW, register, all 305e27a552bSJed Brown 306e27a552bSJed Brown .seealso: TSRosWRegisterDestroy() 307e27a552bSJed Brown @*/ 308e27a552bSJed Brown PetscErrorCode TSRosWRegisterAll(void) 309e27a552bSJed Brown { 310e27a552bSJed Brown PetscErrorCode ierr; 311e27a552bSJed Brown 312e27a552bSJed Brown PetscFunctionBegin; 313e27a552bSJed Brown if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 314e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_TRUE; 315e27a552bSJed Brown 316e27a552bSJed Brown { 3173606a31eSEmil Constantinescu const PetscReal 3183606a31eSEmil Constantinescu A = 0, 3193606a31eSEmil Constantinescu Gamma = 1, 3201f80e275SEmil Constantinescu b = 1, 3211f80e275SEmil Constantinescu binterpt=1; 3221f80e275SEmil Constantinescu 3231f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);CHKERRQ(ierr); 3243606a31eSEmil Constantinescu } 3253606a31eSEmil Constantinescu 3263606a31eSEmil Constantinescu { 3273606a31eSEmil Constantinescu const PetscReal 3283606a31eSEmil Constantinescu A= 0, 3293606a31eSEmil Constantinescu Gamma = 0.5, 3301f80e275SEmil Constantinescu b = 1, 3311f80e275SEmil Constantinescu binterpt=1; 3321f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,PETSC_NULL,1,&binterpt);CHKERRQ(ierr); 3333606a31eSEmil Constantinescu } 3343606a31eSEmil Constantinescu 3353606a31eSEmil Constantinescu { 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 @*/ 700f9c1d6abSBarry Smith PetscErrorCode TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 701f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 702e27a552bSJed Brown { 703e27a552bSJed Brown PetscErrorCode ierr; 70461692a83SJed Brown RosWTableauLink link; 70561692a83SJed Brown RosWTableau t; 70661692a83SJed Brown PetscInt i,j,k; 70761692a83SJed Brown PetscScalar *GammaInv; 708e27a552bSJed Brown 709e27a552bSJed Brown PetscFunctionBegin; 710fe7e6d57SJed Brown PetscValidCharPointer(name,1); 711fe7e6d57SJed Brown PetscValidPointer(A,4); 712fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 713fe7e6d57SJed Brown PetscValidPointer(b,6); 714fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 715fe7e6d57SJed Brown 716e27a552bSJed Brown ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 717e27a552bSJed Brown ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 718e27a552bSJed Brown t = &link->tab; 719e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 720e27a552bSJed Brown t->order = order; 721e27a552bSJed Brown t->s = s; 72261692a83SJed 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); 72343b21953SEmil 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); 724e27a552bSJed Brown ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 72561692a83SJed Brown ierr = PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 72643b21953SEmil Constantinescu ierr = PetscMemcpy(t->GammaExplicitCorr,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 72761692a83SJed Brown ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 728fe7e6d57SJed Brown if (bembed) { 729fe7e6d57SJed Brown ierr = PetscMalloc2(s,PetscReal,&t->bembed,s,PetscReal,&t->bembedt);CHKERRQ(ierr); 730fe7e6d57SJed Brown ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 731fe7e6d57SJed Brown } 73261692a83SJed Brown for (i=0; i<s; i++) { 73361692a83SJed Brown t->ASum[i] = 0; 73461692a83SJed Brown t->GammaSum[i] = 0; 73561692a83SJed Brown for (j=0; j<s; j++) { 73661692a83SJed Brown t->ASum[i] += A[i*s+j]; 737fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 73861692a83SJed Brown } 73961692a83SJed Brown } 74061692a83SJed Brown ierr = PetscMalloc(s*s*sizeof(PetscScalar),&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 74161692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 742fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 743fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 744fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 745c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 746fd96d5b0SEmil Constantinescu } else { 747c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 748fd96d5b0SEmil Constantinescu } 749fd96d5b0SEmil Constantinescu } 750fd96d5b0SEmil Constantinescu 75161692a83SJed Brown switch (s) { 75261692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 75396b95a6bSBarry Smith case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0);CHKERRQ(ierr); break; 75496b95a6bSBarry Smith case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0);CHKERRQ(ierr); break; 75596b95a6bSBarry Smith case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0);CHKERRQ(ierr); break; 75661692a83SJed Brown case 5: { 75761692a83SJed Brown PetscInt ipvt5[5]; 75861692a83SJed Brown MatScalar work5[5*5]; 75996b95a6bSBarry Smith ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0);CHKERRQ(ierr); break; 76061692a83SJed Brown } 76196b95a6bSBarry Smith case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0);CHKERRQ(ierr); break; 76296b95a6bSBarry Smith case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0);CHKERRQ(ierr); break; 76361692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 76461692a83SJed Brown } 76561692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 76661692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 76743b21953SEmil Constantinescu 76843b21953SEmil Constantinescu for (i=0; i<s; i++) { 76943b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 77043b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 77143b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 77243b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 77343b21953SEmil Constantinescu } 77443b21953SEmil Constantinescu } 77543b21953SEmil Constantinescu } 77643b21953SEmil Constantinescu 77761692a83SJed Brown for (i=0; i<s; i++) { 77861692a83SJed Brown for (j=0; j<s; j++) { 77961692a83SJed Brown t->At[i*s+j] = 0; 78061692a83SJed Brown for (k=0; k<s; k++) { 78161692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 78261692a83SJed Brown } 78361692a83SJed Brown } 78461692a83SJed Brown t->bt[i] = 0; 78561692a83SJed Brown for (j=0; j<s; j++) { 78661692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 78761692a83SJed Brown } 788fe7e6d57SJed Brown if (bembed) { 789fe7e6d57SJed Brown t->bembedt[i] = 0; 790fe7e6d57SJed Brown for (j=0; j<s; j++) { 791fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 792fe7e6d57SJed Brown } 793fe7e6d57SJed Brown } 79461692a83SJed Brown } 7958d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 7968d59e960SJed Brown 797f4aed992SEmil Constantinescu t->pinterp = pinterp; 7983ca35412SEmil Constantinescu ierr = PetscMalloc(s*pinterp*sizeof(binterpt[0]),&t->binterpt);CHKERRQ(ierr); 7993ca35412SEmil Constantinescu ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 80061692a83SJed Brown link->next = RosWTableauList; 80161692a83SJed Brown RosWTableauList = link; 802e27a552bSJed Brown PetscFunctionReturn(0); 803e27a552bSJed Brown } 804e27a552bSJed Brown 805e27a552bSJed Brown #undef __FUNCT__ 80642faf41dSJed Brown #define __FUNCT__ "TSRosWRegisterRos4" 80742faf41dSJed Brown /*@C 80842faf41dSJed Brown TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing paramter choices 80942faf41dSJed Brown 81042faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 81142faf41dSJed Brown 81242faf41dSJed Brown Input Parameters: 81342faf41dSJed Brown + name - identifier for method 81442faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 81542faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 81642faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 81742faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 81842faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 81942faf41dSJed Brown . e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 82042faf41dSJed Brown 82142faf41dSJed Brown Notes: 82242faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 82342faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 82442faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 82542faf41dSJed Brown 82642faf41dSJed Brown Level: developer 82742faf41dSJed Brown 82842faf41dSJed Brown .keywords: TS, register 82942faf41dSJed Brown 83042faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 83142faf41dSJed Brown @*/ 83219fd82e9SBarry Smith PetscErrorCode TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 83342faf41dSJed Brown { 83442faf41dSJed Brown PetscErrorCode ierr; 83542faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 83642faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 83742faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 83842faf41dSJed Brown p42 = one/eight - gamma/three, 83942faf41dSJed Brown p43 = one/twelve - gamma/three, 84042faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 84142faf41dSJed Brown p56 = one/twenty - gamma/four; 84242faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 84342faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 84442faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 84542faf41dSJed Brown 84642faf41dSJed Brown PetscFunctionBegin; 84742faf41dSJed Brown /* Step 1: choose Gamma (input) */ 84842faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 84942faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 85042faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 85142faf41dSJed Brown 85242faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 85342faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 85442faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 85542faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 85642faf41dSJed Brown rhs[0] = one - b3; 85742faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 85842faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 85942faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 86042faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 86142faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 86242faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 86342faf41dSJed Brown 86442faf41dSJed Brown /* Step 3 */ 86542faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 86642faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 86742faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 86842faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 86942faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 87042faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 87142faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 87242faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 87342faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 87442faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 87542faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 87642faf41dSJed Brown 87742faf41dSJed Brown /* Step 4: back-substitute */ 87842faf41dSJed Brown beta32 = beta32beta2p / beta2p; 87942faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 88042faf41dSJed Brown 88142faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 88242faf41dSJed Brown a43 = 0; 88342faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 88442faf41dSJed Brown a42 = a32; 88542faf41dSJed Brown 88642faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 88742faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 88842faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 88942faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 89042faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 89142faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 89242faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 89342faf41dSJed 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; 89442faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 89542faf41dSJed Brown 89642faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 89742faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 89842faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 89942faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 90042faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 90142faf41dSJed Brown 90242faf41dSJed Brown { 90342faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 90442faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 90542faf41dSJed Brown } 90642faf41dSJed Brown ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,PETSC_NULL);CHKERRQ(ierr); 90742faf41dSJed Brown PetscFunctionReturn(0); 90842faf41dSJed Brown } 90942faf41dSJed Brown 91042faf41dSJed Brown #undef __FUNCT__ 9111c3436cfSJed Brown #define __FUNCT__ "TSEvaluateStep_RosW" 9121c3436cfSJed Brown /* 9131c3436cfSJed Brown The step completion formula is 9141c3436cfSJed Brown 9151c3436cfSJed Brown x1 = x0 + b^T Y 9161c3436cfSJed Brown 9171c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 9181c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 9191c3436cfSJed Brown 9201c3436cfSJed Brown x1e = x0 + be^T Y 9211c3436cfSJed Brown = x1 - b^T Y + be^T Y 9221c3436cfSJed Brown = x1 + (be - b)^T Y 9231c3436cfSJed Brown 9241c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9251c3436cfSJed Brown */ 926f9c1d6abSBarry Smith static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done) 9271c3436cfSJed Brown { 9281c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9291c3436cfSJed Brown RosWTableau tab = ros->tableau; 9301c3436cfSJed Brown PetscScalar *w = ros->work; 9311c3436cfSJed Brown PetscInt i; 9321c3436cfSJed Brown PetscErrorCode ierr; 9331c3436cfSJed Brown 9341c3436cfSJed Brown PetscFunctionBegin; 9351c3436cfSJed Brown if (order == tab->order) { 936108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 937f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 938de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 939f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 940f9c1d6abSBarry Smith } else {ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr);} 9411c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9421c3436cfSJed Brown PetscFunctionReturn(0); 9431c3436cfSJed Brown } else if (order == tab->order-1) { 9441c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 945108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 946f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 947de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 948f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 949108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 950108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 951f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 952f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 9531c3436cfSJed Brown } 9541c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9551c3436cfSJed Brown PetscFunctionReturn(0); 9561c3436cfSJed Brown } 9571c3436cfSJed Brown unavailable: 9581c3436cfSJed Brown if (done) *done = PETSC_FALSE; 9591c3436cfSJed 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); 9601c3436cfSJed Brown PetscFunctionReturn(0); 9611c3436cfSJed Brown } 9621c3436cfSJed Brown 9631c3436cfSJed Brown #undef __FUNCT__ 964e27a552bSJed Brown #define __FUNCT__ "TSStep_RosW" 965e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 966e27a552bSJed Brown { 96761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 96861692a83SJed Brown RosWTableau tab = ros->tableau; 969e27a552bSJed Brown const PetscInt s = tab->s; 9701c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9710feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 972c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 97361692a83SJed Brown PetscScalar *w = ros->work; 9747d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 975e27a552bSJed Brown SNES snes; 9761c3436cfSJed Brown TSAdapt adapt; 9771c3436cfSJed Brown PetscInt i,j,its,lits,reject,next_scheme; 978cdbf8f93SLisandro Dalcin PetscReal next_time_step; 9791c3436cfSJed Brown PetscBool accept; 980e27a552bSJed Brown PetscErrorCode ierr; 9810feba352SEmil Constantinescu MatStructure str; 982e27a552bSJed Brown 983e27a552bSJed Brown PetscFunctionBegin; 984e27a552bSJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 985cdbf8f93SLisandro Dalcin next_time_step = ts->time_step; 9861c3436cfSJed Brown accept = PETSC_TRUE; 987108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 988e27a552bSJed Brown 98997335746SJed Brown for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 9901c3436cfSJed Brown const PetscReal h = ts->time_step; 991b8123daeSJed Brown ierr = TSPreStep(ts);CHKERRQ(ierr); 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]; 995b8123daeSJed Brown ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 996c17803e7SJed Brown if (GammaZeroDiag[i]) { 997c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 998*b296d7d5SJed Brown ros->scoeff = 1.; 999c17803e7SJed Brown } else { 1000c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 1001*b296d7d5SJed Brown ros->scoeff = 1./Gamma[i*s+i]; 1002fd96d5b0SEmil Constantinescu } 100361692a83SJed Brown 100461692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 1005de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 1006de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 100761692a83SJed Brown 100861692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 100961692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 101061692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 101161692a83SJed Brown 1012e27a552bSJed Brown /* Initial guess taken from last stage */ 101361692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 101461692a83SJed Brown 10157d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 101661692a83SJed Brown if (!ros->recompute_jacobian && !i) { 101761692a83SJed Brown ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 101861692a83SJed Brown } 101961692a83SJed Brown ierr = SNESSolve(snes,PETSC_NULL,Y[i]);CHKERRQ(ierr); 1020e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1021e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10225ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 102397335746SJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 102497335746SJed Brown ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 102597335746SJed Brown if (!accept) goto reject_step; 10267d4bf2deSEmil Constantinescu } else { 10271ce71dffSSatish Balay Mat J,Jp; 10280feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10290feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 10300feba352SEmil Constantinescu ierr = VecScale(Y[i],-1.0); 10310feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/ 10320feba352SEmil Constantinescu 10330feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10340feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10350feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 10360feba352SEmil Constantinescu /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10370feba352SEmil Constantinescu str = SAME_NONZERO_PATTERN; 1038ccbc64bcSJed Brown ierr = TSGetIJacobian(ts,&J,&Jp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 10390feba352SEmil Constantinescu ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);CHKERRQ(ierr); 10400feba352SEmil Constantinescu ierr = MatMult(J,Zstage,Zdot); 10410feba352SEmil Constantinescu 10420feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10430feba352SEmil Constantinescu ierr = VecScale(Y[i],h); 10445ef26d82SJed Brown ts->ksp_its += 1; 10457d4bf2deSEmil Constantinescu } 1046e27a552bSJed Brown } 10471c3436cfSJed Brown ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 1048108c343cSJed Brown ros->status = TS_STEP_PENDING; 1049e27a552bSJed Brown 10501c3436cfSJed Brown /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 10511c3436cfSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 10521c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10538d59e960SJed Brown ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 10541c3436cfSJed Brown ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 10551c3436cfSJed Brown if (accept) { 10561c3436cfSJed Brown /* ignore next_scheme for now */ 1057e27a552bSJed Brown ts->ptime += ts->time_step; 1058cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 1059e27a552bSJed Brown ts->steps++; 1060108c343cSJed Brown ros->status = TS_STEP_COMPLETE; 10611c3436cfSJed Brown break; 10621c3436cfSJed Brown } else { /* Roll back the current step */ 10631c3436cfSJed Brown for (i=0; i<s; i++) w[i] = -tab->bt[i]; 10641c3436cfSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 10651c3436cfSJed Brown ts->time_step = next_time_step; 1066108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 10671c3436cfSJed Brown } 1068476b6736SJed Brown reject_step: continue; 10691c3436cfSJed Brown } 1070b2ce242eSJed Brown if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 1071e27a552bSJed Brown PetscFunctionReturn(0); 1072e27a552bSJed Brown } 1073e27a552bSJed Brown 1074e27a552bSJed Brown #undef __FUNCT__ 1075e27a552bSJed Brown #define __FUNCT__ "TSInterpolate_RosW" 1076f9c1d6abSBarry Smith static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1077e27a552bSJed Brown { 107861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1079f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1080f4aed992SEmil Constantinescu PetscReal h; 1081f4aed992SEmil Constantinescu PetscReal tt,t; 1082f4aed992SEmil Constantinescu PetscScalar *bt; 1083f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1084f4aed992SEmil Constantinescu PetscErrorCode ierr; 1085f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1086f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1087f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1088e27a552bSJed Brown 1089e27a552bSJed Brown PetscFunctionBegin; 1090f4aed992SEmil Constantinescu if (!Bt) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1091f4aed992SEmil Constantinescu 1092f4aed992SEmil Constantinescu switch (ros->status) { 1093f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1094f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1095f4aed992SEmil Constantinescu h = ts->time_step; 1096f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1097f4aed992SEmil Constantinescu break; 1098f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1099f4aed992SEmil Constantinescu h = ts->time_step_prev; 1100f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1101f4aed992SEmil Constantinescu break; 1102f4aed992SEmil Constantinescu default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1103f4aed992SEmil Constantinescu } 11043ca35412SEmil Constantinescu ierr = PetscMalloc(s*sizeof(bt[0]),&bt);CHKERRQ(ierr); 1105f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1106f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1107f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 11083ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1109f4aed992SEmil Constantinescu } 1110f4aed992SEmil Constantinescu } 1111f4aed992SEmil Constantinescu 1112f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1113f9c1d6abSBarry Smith /*U<-0*/ 1114f9c1d6abSBarry Smith ierr = VecZeroEntries(U);CHKERRQ(ierr); 1115f4aed992SEmil Constantinescu 1116f9c1d6abSBarry Smith /*U<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11173ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j]=0; 11183ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11193ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11203ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1121f4aed992SEmil Constantinescu } 11223ca35412SEmil Constantinescu } 1123f9c1d6abSBarry Smith ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1124f4aed992SEmil Constantinescu 1125f4aed992SEmil Constantinescu /*X<-y(t) + X*/ 1126f9c1d6abSBarry Smith ierr = VecAXPY(U,1.0,ros->VecSolPrev);CHKERRQ(ierr); 1127f4aed992SEmil Constantinescu 1128f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1129f4aed992SEmil Constantinescu 1130e27a552bSJed Brown PetscFunctionReturn(0); 1131e27a552bSJed Brown } 1132e27a552bSJed Brown 1133e27a552bSJed Brown /*------------------------------------------------------------*/ 1134e27a552bSJed Brown #undef __FUNCT__ 1135e27a552bSJed Brown #define __FUNCT__ "TSReset_RosW" 1136e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1137e27a552bSJed Brown { 113861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1139e27a552bSJed Brown PetscInt s; 1140e27a552bSJed Brown PetscErrorCode ierr; 1141e27a552bSJed Brown 1142e27a552bSJed Brown PetscFunctionBegin; 114361692a83SJed Brown if (!ros->tableau) PetscFunctionReturn(0); 114461692a83SJed Brown s = ros->tableau->s; 114561692a83SJed Brown ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 114661692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 114761692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 114861692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 114961692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 11503ca35412SEmil Constantinescu ierr = VecDestroy(&ros->VecSolPrev);CHKERRQ(ierr); 115161692a83SJed Brown ierr = PetscFree(ros->work);CHKERRQ(ierr); 1152e27a552bSJed Brown PetscFunctionReturn(0); 1153e27a552bSJed Brown } 1154e27a552bSJed Brown 1155e27a552bSJed Brown #undef __FUNCT__ 1156e27a552bSJed Brown #define __FUNCT__ "TSDestroy_RosW" 1157e27a552bSJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1158e27a552bSJed Brown { 1159e27a552bSJed Brown PetscErrorCode ierr; 1160e27a552bSJed Brown 1161e27a552bSJed Brown PetscFunctionBegin; 1162e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1163e27a552bSJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 1164e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 1165e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);CHKERRQ(ierr); 116661692a83SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",PETSC_NULL);CHKERRQ(ierr); 1167e27a552bSJed Brown PetscFunctionReturn(0); 1168e27a552bSJed Brown } 1169e27a552bSJed Brown 1170d5e6173cSPeter Brune 1171d5e6173cSPeter Brune #undef __FUNCT__ 1172d5e6173cSPeter Brune #define __FUNCT__ "TSRosWGetVecs" 1173d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1174d5e6173cSPeter Brune { 1175d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1176d5e6173cSPeter Brune PetscErrorCode ierr; 1177d5e6173cSPeter Brune 1178d5e6173cSPeter Brune PetscFunctionBegin; 1179d5e6173cSPeter Brune if (Ydot) { 1180d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1181d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1182d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1183d5e6173cSPeter Brune } 1184d5e6173cSPeter Brune if (Zdot) { 1185d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1186d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1187d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1188d5e6173cSPeter Brune } 1189d5e6173cSPeter Brune if (Ystage) { 1190d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1191d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1192d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1193d5e6173cSPeter Brune } 1194d5e6173cSPeter Brune if (Zstage) { 1195d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1196d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1197d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1198d5e6173cSPeter Brune } 1199d5e6173cSPeter Brune 1200d5e6173cSPeter Brune PetscFunctionReturn(0); 1201d5e6173cSPeter Brune } 1202d5e6173cSPeter Brune 1203d5e6173cSPeter Brune 1204d5e6173cSPeter Brune #undef __FUNCT__ 1205d5e6173cSPeter Brune #define __FUNCT__ "TSRosWRestoreVecs" 1206d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1207d5e6173cSPeter Brune { 1208d5e6173cSPeter Brune PetscErrorCode ierr; 1209d5e6173cSPeter Brune 1210d5e6173cSPeter Brune PetscFunctionBegin; 1211d5e6173cSPeter Brune if (Ydot) { 1212d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1213d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1214d5e6173cSPeter Brune } 1215d5e6173cSPeter Brune } 1216d5e6173cSPeter Brune if (Zdot) { 1217d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1218d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1219d5e6173cSPeter Brune } 1220d5e6173cSPeter Brune } 1221d5e6173cSPeter Brune if (Ystage) { 1222d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1223d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1224d5e6173cSPeter Brune } 1225d5e6173cSPeter Brune } 1226d5e6173cSPeter Brune if (Zstage) { 1227d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1228d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1229d5e6173cSPeter Brune } 1230d5e6173cSPeter Brune } 1231d5e6173cSPeter Brune PetscFunctionReturn(0); 1232d5e6173cSPeter Brune } 1233d5e6173cSPeter Brune 1234d5e6173cSPeter Brune #undef __FUNCT__ 1235d5e6173cSPeter Brune #define __FUNCT__ "DMCoarsenHook_TSRosW" 1236d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1237d5e6173cSPeter Brune { 1238d5e6173cSPeter Brune 1239d5e6173cSPeter Brune PetscFunctionBegin; 1240d5e6173cSPeter Brune PetscFunctionReturn(0); 1241d5e6173cSPeter Brune } 1242d5e6173cSPeter Brune 1243d5e6173cSPeter Brune #undef __FUNCT__ 1244d5e6173cSPeter Brune #define __FUNCT__ "DMRestrictHook_TSRosW" 1245d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1246d5e6173cSPeter Brune { 1247d5e6173cSPeter Brune TS ts = (TS)ctx; 1248d5e6173cSPeter Brune PetscErrorCode ierr; 1249d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1250d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1251d5e6173cSPeter Brune 1252d5e6173cSPeter Brune PetscFunctionBegin; 1253d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1254d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1255d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1256d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1257d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1258d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1259d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1260d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1261d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1262d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1263d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1264d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1265d5e6173cSPeter Brune PetscFunctionReturn(0); 1266d5e6173cSPeter Brune } 1267d5e6173cSPeter Brune 1268e27a552bSJed Brown /* 1269e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1270e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1271e27a552bSJed Brown */ 1272e27a552bSJed Brown #undef __FUNCT__ 1273e27a552bSJed Brown #define __FUNCT__ "SNESTSFormFunction_RosW" 1274f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1275e27a552bSJed Brown { 127661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1277e27a552bSJed Brown PetscErrorCode ierr; 1278d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1279*b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1280d5e6173cSPeter Brune DM dm,dmsave; 1281e27a552bSJed Brown 1282e27a552bSJed Brown PetscFunctionBegin; 1283d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1284d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1285*b296d7d5SJed Brown ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1286f9c1d6abSBarry Smith ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1287d5e6173cSPeter Brune dmsave = ts->dm; 1288d5e6173cSPeter Brune ts->dm = dm; 1289d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1290d5e6173cSPeter Brune ts->dm = dmsave; 1291d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1292e27a552bSJed Brown PetscFunctionReturn(0); 1293e27a552bSJed Brown } 1294e27a552bSJed Brown 1295e27a552bSJed Brown #undef __FUNCT__ 1296e27a552bSJed Brown #define __FUNCT__ "SNESTSFormJacobian_RosW" 1297f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *str,TS ts) 1298e27a552bSJed Brown { 129961692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1300d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1301*b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1302e27a552bSJed Brown PetscErrorCode ierr; 1303d5e6173cSPeter Brune DM dm,dmsave; 1304e27a552bSJed Brown 1305e27a552bSJed Brown PetscFunctionBegin; 130661692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1307d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1308d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1309d5e6173cSPeter Brune dmsave = ts->dm; 1310d5e6173cSPeter Brune ts->dm = dm; 1311*b296d7d5SJed Brown ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 1312d5e6173cSPeter Brune ts->dm = dmsave; 1313d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1314e27a552bSJed Brown PetscFunctionReturn(0); 1315e27a552bSJed Brown } 1316e27a552bSJed Brown 1317e27a552bSJed Brown #undef __FUNCT__ 1318e27a552bSJed Brown #define __FUNCT__ "TSSetUp_RosW" 1319e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1320e27a552bSJed Brown { 132161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 132261692a83SJed Brown RosWTableau tab = ros->tableau; 1323e27a552bSJed Brown PetscInt s = tab->s; 1324e27a552bSJed Brown PetscErrorCode ierr; 1325d5e6173cSPeter Brune DM dm; 1326e27a552bSJed Brown 1327e27a552bSJed Brown PetscFunctionBegin; 132861692a83SJed Brown if (!ros->tableau) { 1329e27a552bSJed Brown ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1330e27a552bSJed Brown } 133161692a83SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 133261692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 133361692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 133461692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 133561692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 13363ca35412SEmil Constantinescu ierr = VecDuplicate(ts->vec_sol,&ros->VecSolPrev);CHKERRQ(ierr); 133761692a83SJed Brown ierr = PetscMalloc(s*sizeof(ros->work[0]),&ros->work);CHKERRQ(ierr); 1338d5e6173cSPeter Brune ierr = TSGetDM(ts,&dm); 1339d5e6173cSPeter Brune if (dm) { 1340d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1341d5e6173cSPeter Brune } 1342e27a552bSJed Brown PetscFunctionReturn(0); 1343e27a552bSJed Brown } 1344e27a552bSJed Brown /*------------------------------------------------------------*/ 1345e27a552bSJed Brown 1346e27a552bSJed Brown #undef __FUNCT__ 1347e27a552bSJed Brown #define __FUNCT__ "TSSetFromOptions_RosW" 1348e27a552bSJed Brown static PetscErrorCode TSSetFromOptions_RosW(TS ts) 1349e27a552bSJed Brown { 135061692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1351e27a552bSJed Brown PetscErrorCode ierr; 135261692a83SJed Brown char rostype[256]; 1353e27a552bSJed Brown 1354e27a552bSJed Brown PetscFunctionBegin; 1355e27a552bSJed Brown ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 1356e27a552bSJed Brown { 135761692a83SJed Brown RosWTableauLink link; 1358e27a552bSJed Brown PetscInt count,choice; 1359e27a552bSJed Brown PetscBool flg; 1360e27a552bSJed Brown const char **namelist; 136161692a83SJed Brown SNES snes; 136261692a83SJed Brown 13638caf3d72SBarry Smith ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof(rostype));CHKERRQ(ierr); 136461692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1365e27a552bSJed Brown ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 136661692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 136761692a83SJed Brown ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 136861692a83SJed Brown ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 1369e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 137061692a83SJed Brown 137161692a83SJed Brown ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,PETSC_NULL);CHKERRQ(ierr); 137261692a83SJed Brown 137361692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 137461692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 137561692a83SJed Brown if (!((PetscObject)snes)->type_name) { 137661692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 137761692a83SJed Brown } 137861692a83SJed Brown ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1379e27a552bSJed Brown } 1380e27a552bSJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 1381e27a552bSJed Brown PetscFunctionReturn(0); 1382e27a552bSJed Brown } 1383e27a552bSJed Brown 1384e27a552bSJed Brown #undef __FUNCT__ 1385e27a552bSJed Brown #define __FUNCT__ "PetscFormatRealArray" 1386e27a552bSJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1387e27a552bSJed Brown { 1388e27a552bSJed Brown PetscErrorCode ierr; 1389e408995aSJed Brown PetscInt i; 1390e408995aSJed Brown size_t left,count; 1391e27a552bSJed Brown char *p; 1392e27a552bSJed Brown 1393e27a552bSJed Brown PetscFunctionBegin; 1394e408995aSJed Brown for (i=0,p=buf,left=len; i<n; i++) { 1395e408995aSJed Brown ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1396e27a552bSJed Brown if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1397e27a552bSJed Brown left -= count; 1398e27a552bSJed Brown p += count; 1399e27a552bSJed Brown *p++ = ' '; 1400e27a552bSJed Brown } 1401e27a552bSJed Brown p[i ? 0 : -1] = 0; 1402e27a552bSJed Brown PetscFunctionReturn(0); 1403e27a552bSJed Brown } 1404e27a552bSJed Brown 1405e27a552bSJed Brown #undef __FUNCT__ 1406e27a552bSJed Brown #define __FUNCT__ "TSView_RosW" 1407e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1408e27a552bSJed Brown { 140961692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 141061692a83SJed Brown RosWTableau tab = ros->tableau; 1411e27a552bSJed Brown PetscBool iascii; 1412e27a552bSJed Brown PetscErrorCode ierr; 1413ef20d060SBarry Smith TSAdapt adapt; 1414e27a552bSJed Brown 1415e27a552bSJed Brown PetscFunctionBegin; 1416251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1417e27a552bSJed Brown if (iascii) { 141819fd82e9SBarry Smith TSRosWType rostype; 1419e408995aSJed Brown PetscInt i; 1420e408995aSJed Brown PetscReal abscissa[512]; 1421e27a552bSJed Brown char buf[512]; 142261692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 142361692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 14248caf3d72SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 142561692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1426e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 14278caf3d72SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1428e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1429e27a552bSJed Brown } 1430ef20d060SBarry Smith ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1431ef20d060SBarry Smith ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1432e27a552bSJed Brown ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1433e27a552bSJed Brown PetscFunctionReturn(0); 1434e27a552bSJed Brown } 1435e27a552bSJed Brown 1436e27a552bSJed Brown #undef __FUNCT__ 1437e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType" 1438e27a552bSJed Brown /*@C 143961692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1440e27a552bSJed Brown 1441e27a552bSJed Brown Logically collective 1442e27a552bSJed Brown 1443e27a552bSJed Brown Input Parameter: 1444e27a552bSJed Brown + ts - timestepping context 144561692a83SJed Brown - rostype - type of Rosenbrock-W scheme 1446e27a552bSJed Brown 1447020d8f30SJed Brown Level: beginner 1448e27a552bSJed Brown 1449020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1450e27a552bSJed Brown @*/ 145119fd82e9SBarry Smith PetscErrorCode TSRosWSetType(TS ts,TSRosWType rostype) 1452e27a552bSJed Brown { 1453e27a552bSJed Brown PetscErrorCode ierr; 1454e27a552bSJed Brown 1455e27a552bSJed Brown PetscFunctionBegin; 1456e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 145719fd82e9SBarry Smith ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,rostype));CHKERRQ(ierr); 1458e27a552bSJed Brown PetscFunctionReturn(0); 1459e27a552bSJed Brown } 1460e27a552bSJed Brown 1461e27a552bSJed Brown #undef __FUNCT__ 1462e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType" 1463e27a552bSJed Brown /*@C 146461692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1465e27a552bSJed Brown 1466e27a552bSJed Brown Logically collective 1467e27a552bSJed Brown 1468e27a552bSJed Brown Input Parameter: 1469e27a552bSJed Brown . ts - timestepping context 1470e27a552bSJed Brown 1471e27a552bSJed Brown Output Parameter: 147261692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1473e27a552bSJed Brown 1474e27a552bSJed Brown Level: intermediate 1475e27a552bSJed Brown 1476e27a552bSJed Brown .seealso: TSRosWGetType() 1477e27a552bSJed Brown @*/ 147819fd82e9SBarry Smith PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1479e27a552bSJed Brown { 1480e27a552bSJed Brown PetscErrorCode ierr; 1481e27a552bSJed Brown 1482e27a552bSJed Brown PetscFunctionBegin; 1483e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 148419fd82e9SBarry Smith ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1485e27a552bSJed Brown PetscFunctionReturn(0); 1486e27a552bSJed Brown } 1487e27a552bSJed Brown 1488e27a552bSJed Brown #undef __FUNCT__ 148961692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian" 1490e27a552bSJed Brown /*@C 149161692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1492e27a552bSJed Brown 1493e27a552bSJed Brown Logically collective 1494e27a552bSJed Brown 1495e27a552bSJed Brown Input Parameter: 1496e27a552bSJed Brown + ts - timestepping context 149761692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1498e27a552bSJed Brown 1499e27a552bSJed Brown Level: intermediate 1500e27a552bSJed Brown 1501e27a552bSJed Brown .seealso: TSRosWGetType() 1502e27a552bSJed Brown @*/ 150361692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1504e27a552bSJed Brown { 1505e27a552bSJed Brown PetscErrorCode ierr; 1506e27a552bSJed Brown 1507e27a552bSJed Brown PetscFunctionBegin; 1508e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 150961692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1510e27a552bSJed Brown PetscFunctionReturn(0); 1511e27a552bSJed Brown } 1512e27a552bSJed Brown 1513e27a552bSJed Brown EXTERN_C_BEGIN 1514e27a552bSJed Brown #undef __FUNCT__ 1515e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType_RosW" 151619fd82e9SBarry Smith PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1517e27a552bSJed Brown { 151861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1519e27a552bSJed Brown PetscErrorCode ierr; 1520e27a552bSJed Brown 1521e27a552bSJed Brown PetscFunctionBegin; 152261692a83SJed Brown if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 152361692a83SJed Brown *rostype = ros->tableau->name; 1524e27a552bSJed Brown PetscFunctionReturn(0); 1525e27a552bSJed Brown } 1526ef20d060SBarry Smith 1527e27a552bSJed Brown #undef __FUNCT__ 1528e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType_RosW" 152919fd82e9SBarry Smith PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1530e27a552bSJed Brown { 153161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1532e27a552bSJed Brown PetscErrorCode ierr; 1533e27a552bSJed Brown PetscBool match; 153461692a83SJed Brown RosWTableauLink link; 1535e27a552bSJed Brown 1536e27a552bSJed Brown PetscFunctionBegin; 153761692a83SJed Brown if (ros->tableau) { 153861692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1539e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1540e27a552bSJed Brown } 154161692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 154261692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1543e27a552bSJed Brown if (match) { 1544e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 154561692a83SJed Brown ros->tableau = &link->tab; 1546e27a552bSJed Brown PetscFunctionReturn(0); 1547e27a552bSJed Brown } 1548e27a552bSJed Brown } 154961692a83SJed Brown SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1550e27a552bSJed Brown PetscFunctionReturn(0); 1551e27a552bSJed Brown } 155261692a83SJed Brown 1553e27a552bSJed Brown #undef __FUNCT__ 155461692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 155561692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1556e27a552bSJed Brown { 155761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1558e27a552bSJed Brown 1559e27a552bSJed Brown PetscFunctionBegin; 156061692a83SJed Brown ros->recompute_jacobian = flg; 1561e27a552bSJed Brown PetscFunctionReturn(0); 1562e27a552bSJed Brown } 1563e27a552bSJed Brown EXTERN_C_END 1564e27a552bSJed Brown 1565d5e6173cSPeter Brune 1566e27a552bSJed Brown /* ------------------------------------------------------------ */ 1567e27a552bSJed Brown /*MC 1568020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1569e27a552bSJed Brown 1570e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1571e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1572e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1573e27a552bSJed Brown 1574e27a552bSJed Brown Notes: 157561692a83SJed Brown This method currently only works with autonomous ODE and DAE. 157661692a83SJed Brown 157761692a83SJed Brown Developer notes: 157861692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 157961692a83SJed Brown 1580f9c1d6abSBarry Smith $ udot = f(u) 158161692a83SJed Brown 158261692a83SJed Brown by the stage equations 158361692a83SJed Brown 1584f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 158561692a83SJed Brown 158661692a83SJed Brown and step completion formula 158761692a83SJed Brown 1588f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 158961692a83SJed Brown 1590f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 159161692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 159261692a83SJed Brown we define new variables for the stage equations 159361692a83SJed Brown 159461692a83SJed Brown $ y_i = gamma_ij k_j 159561692a83SJed Brown 159661692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 159761692a83SJed Brown 159861692a83SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 159961692a83SJed Brown 160061692a83SJed Brown to rewrite the method as 160161692a83SJed Brown 1602f9c1d6abSBarry Smith $ [M/(h gamma_ii) - J] y_i = f(u_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1603f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j bt_j y_j 160461692a83SJed Brown 160561692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 160661692a83SJed Brown 160761692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 160861692a83SJed Brown 160961692a83SJed Brown or, more compactly in tensor notation 161061692a83SJed Brown 161161692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 161261692a83SJed Brown 161361692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 161461692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 161561692a83SJed Brown equation 161661692a83SJed Brown 1617f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 161861692a83SJed Brown 161961692a83SJed Brown with initial guess y_i = 0. 1620e27a552bSJed Brown 1621e27a552bSJed Brown Level: beginner 1622e27a552bSJed Brown 1623a4386c9eSJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1624a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1625e27a552bSJed Brown M*/ 1626e27a552bSJed Brown EXTERN_C_BEGIN 1627e27a552bSJed Brown #undef __FUNCT__ 1628e27a552bSJed Brown #define __FUNCT__ "TSCreate_RosW" 1629e27a552bSJed Brown PetscErrorCode TSCreate_RosW(TS ts) 1630e27a552bSJed Brown { 163161692a83SJed Brown TS_RosW *ros; 1632e27a552bSJed Brown PetscErrorCode ierr; 1633e27a552bSJed Brown 1634e27a552bSJed Brown PetscFunctionBegin; 1635e27a552bSJed Brown #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1636e27a552bSJed Brown ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1637e27a552bSJed Brown #endif 1638e27a552bSJed Brown 1639e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1640e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1641e27a552bSJed Brown ts->ops->view = TSView_RosW; 1642e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1643e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1644e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16451c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 1646e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1647e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1648e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1649e27a552bSJed Brown 165061692a83SJed Brown ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 165161692a83SJed Brown ts->data = (void*)ros; 1652e27a552bSJed Brown 1653e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 1654e27a552bSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 165561692a83SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1656e27a552bSJed Brown PetscFunctionReturn(0); 1657e27a552bSJed Brown } 1658e27a552bSJed Brown EXTERN_C_END 1659