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*/ 141e25c274SJed Brown #include <petscdm.h> 15e27a552bSJed Brown 1661692a83SJed Brown #include <../src/mat/blockinvert.h> 1761692a83SJed Brown 1819fd82e9SBarry Smith static TSRosWType TSRosWDefault = TSROSWRA34PW2; 19e27a552bSJed Brown static PetscBool TSRosWRegisterAllCalled; 20e27a552bSJed Brown static PetscBool TSRosWPackageInitialized; 21e27a552bSJed Brown 2261692a83SJed Brown typedef struct _RosWTableau *RosWTableau; 2361692a83SJed Brown struct _RosWTableau { 24e27a552bSJed Brown char *name; 25e27a552bSJed Brown PetscInt order; /* Classical approximation order of the method */ 26e27a552bSJed Brown PetscInt s; /* Number of stages */ 27f4aed992SEmil Constantinescu PetscInt pinterp; /* Interpolation order */ 2861692a83SJed Brown PetscReal *A; /* Propagation table, strictly lower triangular */ 2961692a83SJed Brown PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */ 30c17803e7SJed Brown PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */ 3143b21953SEmil Constantinescu PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/ 3261692a83SJed Brown PetscReal *b; /* Step completion table */ 33fe7e6d57SJed Brown PetscReal *bembed; /* Step completion table for embedded method of order one less */ 3461692a83SJed Brown PetscReal *ASum; /* Row sum of A */ 3561692a83SJed Brown PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */ 3661692a83SJed Brown PetscReal *At; /* Propagation table in transformed variables */ 3761692a83SJed Brown PetscReal *bt; /* Step completion table in transformed variables */ 38fe7e6d57SJed Brown PetscReal *bembedt; /* Step completion table of order one less in transformed variables */ 3961692a83SJed Brown PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */ 408d59e960SJed Brown PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 413ca35412SEmil Constantinescu PetscReal *binterpt; /* Dense output formula */ 42e27a552bSJed Brown }; 4361692a83SJed Brown typedef struct _RosWTableauLink *RosWTableauLink; 4461692a83SJed Brown struct _RosWTableauLink { 4561692a83SJed Brown struct _RosWTableau tab; 4661692a83SJed Brown RosWTableauLink next; 47e27a552bSJed Brown }; 4861692a83SJed Brown static RosWTableauLink RosWTableauList; 49e27a552bSJed Brown 50e27a552bSJed Brown typedef struct { 5161692a83SJed Brown RosWTableau tableau; 5261692a83SJed Brown Vec *Y; /* States computed during the step, used to complete the step */ 53e27a552bSJed Brown Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 5461692a83SJed Brown Vec Ystage; /* Work vector for the state value at each stage */ 5561692a83SJed Brown Vec Zdot; /* Ydot = Zdot + shift*Y */ 5661692a83SJed Brown Vec Zstage; /* Y = Zstage + Y */ 573ca35412SEmil Constantinescu Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation)*/ 581c3436cfSJed Brown PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 59b296d7d5SJed Brown PetscReal scoeff; /* shift = scoeff/dt */ 60e27a552bSJed Brown PetscReal stage_time; 61c17803e7SJed Brown PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */ 6261692a83SJed Brown PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 63108c343cSJed Brown TSStepStatus status; 64e27a552bSJed Brown } TS_RosW; 65e27a552bSJed Brown 66fe7e6d57SJed Brown /*MC 673606a31eSEmil Constantinescu TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method). 683606a31eSEmil Constantinescu 693606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 703606a31eSEmil Constantinescu 713606a31eSEmil Constantinescu Level: intermediate 723606a31eSEmil Constantinescu 733606a31eSEmil Constantinescu .seealso: TSROSW 743606a31eSEmil Constantinescu M*/ 753606a31eSEmil Constantinescu 763606a31eSEmil Constantinescu /*MC 773606a31eSEmil Constantinescu TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method). 783606a31eSEmil Constantinescu 793606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 803606a31eSEmil Constantinescu 813606a31eSEmil Constantinescu Level: intermediate 823606a31eSEmil Constantinescu 833606a31eSEmil Constantinescu .seealso: TSROSW 843606a31eSEmil Constantinescu M*/ 853606a31eSEmil Constantinescu 863606a31eSEmil Constantinescu /*MC 87fe7e6d57SJed Brown TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme. 88fe7e6d57SJed Brown 89fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P. 90fe7e6d57SJed Brown 91fe7e6d57SJed Brown Level: intermediate 92fe7e6d57SJed Brown 93fe7e6d57SJed Brown .seealso: TSROSW 94fe7e6d57SJed Brown M*/ 95fe7e6d57SJed Brown 96fe7e6d57SJed Brown /*MC 97fe7e6d57SJed Brown TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme. 98fe7e6d57SJed Brown 99fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M. 100fe7e6d57SJed Brown 101fe7e6d57SJed Brown Level: intermediate 102fe7e6d57SJed Brown 103fe7e6d57SJed Brown .seealso: TSROSW 104fe7e6d57SJed Brown M*/ 105fe7e6d57SJed Brown 106fe7e6d57SJed Brown /*MC 107fe7e6d57SJed Brown TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1. 108fe7e6d57SJed Brown 109fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 110fe7e6d57SJed Brown 111fe7e6d57SJed 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. 112fe7e6d57SJed Brown 113fe7e6d57SJed Brown References: 114fe7e6d57SJed Brown Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 115fe7e6d57SJed Brown 116fe7e6d57SJed Brown Level: intermediate 117fe7e6d57SJed Brown 118fe7e6d57SJed Brown .seealso: TSROSW 119fe7e6d57SJed Brown M*/ 120fe7e6d57SJed Brown 121fe7e6d57SJed Brown /*MC 122fe7e6d57SJed Brown TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1. 123fe7e6d57SJed Brown 124fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 125fe7e6d57SJed Brown 126fe7e6d57SJed 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. 127fe7e6d57SJed Brown 128fe7e6d57SJed Brown References: 129fe7e6d57SJed Brown Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 130fe7e6d57SJed Brown 131fe7e6d57SJed Brown Level: intermediate 132fe7e6d57SJed Brown 133fe7e6d57SJed Brown .seealso: TSROSW 134fe7e6d57SJed Brown M*/ 135fe7e6d57SJed Brown 136ef3c5b88SJed Brown /*MC 137ef3c5b88SJed Brown TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme 138ef3c5b88SJed Brown 139ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 140ef3c5b88SJed Brown 141ef3c5b88SJed Brown Both the third order and embedded second order methods are stiffly accurate and L-stable. 142ef3c5b88SJed Brown 143ef3c5b88SJed Brown References: 144ef3c5b88SJed Brown Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 145ef3c5b88SJed Brown 146ef3c5b88SJed Brown Level: intermediate 147ef3c5b88SJed Brown 148ef3c5b88SJed Brown .seealso: TSROSW, TSROSWSANDU3 149ef3c5b88SJed Brown M*/ 150ef3c5b88SJed Brown 151ef3c5b88SJed Brown /*MC 152ef3c5b88SJed Brown TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme 153ef3c5b88SJed Brown 154ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 155ef3c5b88SJed Brown 156ef3c5b88SJed Brown The third order method is L-stable, but not stiffly accurate. 157ef3c5b88SJed Brown The second order embedded method is strongly A-stable with R(infty) = 0.5. 158ef3c5b88SJed Brown The internal stages are L-stable. 159ef3c5b88SJed Brown This method is called ROS3 in the paper. 160ef3c5b88SJed Brown 161ef3c5b88SJed Brown References: 162ef3c5b88SJed Brown Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 163ef3c5b88SJed Brown 164ef3c5b88SJed Brown Level: intermediate 165ef3c5b88SJed Brown 166ef3c5b88SJed Brown .seealso: TSROSW, TSROSWRODAS3 167ef3c5b88SJed Brown M*/ 168ef3c5b88SJed Brown 169961f28d0SJed Brown /*MC 170961f28d0SJed Brown TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 171961f28d0SJed Brown 172961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 173961f28d0SJed Brown 174961f28d0SJed Brown A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 175961f28d0SJed Brown 176961f28d0SJed Brown References: 177961f28d0SJed Brown Emil Constantinescu 178961f28d0SJed Brown 179961f28d0SJed Brown Level: intermediate 180961f28d0SJed Brown 18143b21953SEmil Constantinescu .seealso: TSROSW, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, SSP 182961f28d0SJed Brown M*/ 183961f28d0SJed Brown 184961f28d0SJed Brown /*MC 185998eb97aSJed Brown TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 186961f28d0SJed Brown 187961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 188961f28d0SJed Brown 189961f28d0SJed Brown L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 190961f28d0SJed Brown 191961f28d0SJed Brown References: 192961f28d0SJed Brown Emil Constantinescu 193961f28d0SJed Brown 194961f28d0SJed Brown Level: intermediate 195961f28d0SJed Brown 19643b21953SEmil Constantinescu .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLLSSP3P4S2C, TSSSP 197961f28d0SJed Brown M*/ 198961f28d0SJed Brown 199961f28d0SJed Brown /*MC 200998eb97aSJed Brown TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 201961f28d0SJed Brown 202961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 203961f28d0SJed Brown 204961f28d0SJed Brown L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 205961f28d0SJed Brown 206961f28d0SJed Brown References: 207961f28d0SJed Brown Emil Constantinescu 208961f28d0SJed Brown 209961f28d0SJed Brown Level: intermediate 210961f28d0SJed Brown 211961f28d0SJed Brown .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSSSP 212961f28d0SJed Brown M*/ 213961f28d0SJed Brown 21442faf41dSJed Brown /*MC 21542faf41dSJed Brown TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop 21642faf41dSJed Brown 21742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 21842faf41dSJed Brown 21942faf41dSJed Brown A(89.3 degrees)-stable, |R(infty)| = 0.454. 22042faf41dSJed Brown 22142faf41dSJed Brown This method does not provide a dense output formula. 22242faf41dSJed Brown 22342faf41dSJed Brown References: 22442faf41dSJed Brown Kaps and Rentrop, Generalized Runge-Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979. 22542faf41dSJed Brown 22642faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 22742faf41dSJed Brown 22842faf41dSJed Brown Hairer's code ros4.f 22942faf41dSJed Brown 23042faf41dSJed Brown Level: intermediate 23142faf41dSJed Brown 23242faf41dSJed Brown .seealso: TSROSW, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 23342faf41dSJed Brown M*/ 23442faf41dSJed Brown 23542faf41dSJed Brown /*MC 23642faf41dSJed Brown TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine 23742faf41dSJed Brown 23842faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 23942faf41dSJed Brown 24042faf41dSJed Brown A-stable, |R(infty)| = 1/3. 24142faf41dSJed Brown 24242faf41dSJed Brown This method does not provide a dense output formula. 24342faf41dSJed Brown 24442faf41dSJed Brown References: 24542faf41dSJed Brown Shampine, Implementation of Rosenbrock methods, 1982. 24642faf41dSJed Brown 24742faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 24842faf41dSJed Brown 24942faf41dSJed Brown Hairer's code ros4.f 25042faf41dSJed Brown 25142faf41dSJed Brown Level: intermediate 25242faf41dSJed Brown 25342faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWVELDD4, TSROSW4L 25442faf41dSJed Brown M*/ 25542faf41dSJed Brown 25642faf41dSJed Brown /*MC 25742faf41dSJed Brown TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen 25842faf41dSJed Brown 25942faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 26042faf41dSJed Brown 26142faf41dSJed Brown A(89.5 degrees)-stable, |R(infty)| = 0.24. 26242faf41dSJed Brown 26342faf41dSJed Brown This method does not provide a dense output formula. 26442faf41dSJed Brown 26542faf41dSJed Brown References: 26642faf41dSJed Brown van Veldhuizen, D-stability and Kaps-Rentrop methods, 1984. 26742faf41dSJed Brown 26842faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 26942faf41dSJed Brown 27042faf41dSJed Brown Hairer's code ros4.f 27142faf41dSJed Brown 27242faf41dSJed Brown Level: intermediate 27342faf41dSJed Brown 27442faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 27542faf41dSJed Brown M*/ 27642faf41dSJed Brown 27742faf41dSJed Brown /*MC 27842faf41dSJed Brown TSROSW4L - four stage, fourth order Rosenbrock (not W) method 27942faf41dSJed Brown 28042faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 28142faf41dSJed Brown 28242faf41dSJed Brown A-stable and L-stable 28342faf41dSJed Brown 28442faf41dSJed Brown This method does not provide a dense output formula. 28542faf41dSJed Brown 28642faf41dSJed Brown References: 28742faf41dSJed Brown Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 28842faf41dSJed Brown 28942faf41dSJed Brown Hairer's code ros4.f 29042faf41dSJed Brown 29142faf41dSJed Brown Level: intermediate 29242faf41dSJed Brown 29342faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 29442faf41dSJed Brown M*/ 29542faf41dSJed Brown 296e27a552bSJed Brown #undef __FUNCT__ 297e27a552bSJed Brown #define __FUNCT__ "TSRosWRegisterAll" 298e27a552bSJed Brown /*@C 299e27a552bSJed Brown TSRosWRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSRosW 300e27a552bSJed Brown 301e27a552bSJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 302e27a552bSJed Brown 303e27a552bSJed Brown Level: advanced 304e27a552bSJed Brown 305e27a552bSJed Brown .keywords: TS, TSRosW, register, all 306e27a552bSJed Brown 307e27a552bSJed Brown .seealso: TSRosWRegisterDestroy() 308e27a552bSJed Brown @*/ 309e27a552bSJed Brown PetscErrorCode TSRosWRegisterAll(void) 310e27a552bSJed Brown { 311e27a552bSJed Brown PetscErrorCode ierr; 312e27a552bSJed Brown 313e27a552bSJed Brown PetscFunctionBegin; 314e27a552bSJed Brown if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 315e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_TRUE; 316e27a552bSJed Brown 317e27a552bSJed Brown { 318bbd56ea5SKarl Rupp const PetscReal A = 0; 319bbd56ea5SKarl Rupp const PetscReal Gamma = 1; 320bbd56ea5SKarl Rupp const PetscReal b = 1; 321bbd56ea5SKarl Rupp const PetscReal binterpt=1; 3221f80e275SEmil Constantinescu 3230298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3243606a31eSEmil Constantinescu } 3253606a31eSEmil Constantinescu 3263606a31eSEmil Constantinescu { 327bbd56ea5SKarl Rupp const PetscReal A = 0; 328bbd56ea5SKarl Rupp const PetscReal Gamma = 0.5; 329bbd56ea5SKarl Rupp const PetscReal b = 1; 330bbd56ea5SKarl Rupp const PetscReal binterpt=1; 331bbd56ea5SKarl Rupp 3320298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3333606a31eSEmil Constantinescu } 3343606a31eSEmil Constantinescu 3353606a31eSEmil Constantinescu { 336da80777bSKarl Rupp /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0); Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */ 337e27a552bSJed Brown const PetscReal 33861692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 339da80777bSKarl Rupp Gamma[2][2] = {{1.707106781186547524401,0}, {-2.*1.707106781186547524401,1.707106781186547524401}}, 3401c3436cfSJed Brown b[2] = {0.5,0.5}, 3411c3436cfSJed Brown b1[2] = {1.0,0.0}; 3421f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 343da80777bSKarl Rupp binterpt[0][0] = 1.707106781186547524401 - 1.0; 344da80777bSKarl Rupp binterpt[1][0] = 2.0 - 1.707106781186547524401; 345da80777bSKarl Rupp binterpt[0][1] = 1.707106781186547524401 - 1.5; 346da80777bSKarl Rupp binterpt[1][1] = 1.5 - 1.707106781186547524401; 347bbd56ea5SKarl Rupp 3481f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 349e27a552bSJed Brown } 350e27a552bSJed Brown { 351da80777bSKarl Rupp /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0); Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */ 352e27a552bSJed Brown const PetscReal 35361692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 354da80777bSKarl Rupp Gamma[2][2] = {{0.2928932188134524755992,0}, {-2.*0.2928932188134524755992,0.2928932188134524755992}}, 3551c3436cfSJed Brown b[2] = {0.5,0.5}, 3561c3436cfSJed Brown b1[2] = {1.0,0.0}; 3571f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 358da80777bSKarl Rupp binterpt[0][0] = 0.2928932188134524755992 - 1.0; 359da80777bSKarl Rupp binterpt[1][0] = 2.0 - 0.2928932188134524755992; 360da80777bSKarl Rupp binterpt[0][1] = 0.2928932188134524755992 - 1.5; 361da80777bSKarl Rupp binterpt[1][1] = 1.5 - 0.2928932188134524755992; 362bbd56ea5SKarl Rupp 3631f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 364fe7e6d57SJed Brown } 365fe7e6d57SJed Brown { 366da80777bSKarl Rupp /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */ 3671f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 368fe7e6d57SJed Brown const PetscReal 369fe7e6d57SJed Brown A[3][3] = {{0,0,0}, 370fe7e6d57SJed Brown {1.5773502691896257e+00,0,0}, 371fe7e6d57SJed Brown {0.5,0,0}}, 372da80777bSKarl Rupp Gamma[3][3] = {{7.8867513459481287e-01,0,0}, 373da80777bSKarl Rupp {-1.5773502691896257e+00,7.8867513459481287e-01,0}, 374da80777bSKarl Rupp {-6.7075317547305480e-01,-1.7075317547305482e-01,7.8867513459481287e-01}}, 375fe7e6d57SJed Brown b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01}, 376fe7e6d57SJed Brown b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01}; 3771f80e275SEmil Constantinescu 3781f80e275SEmil Constantinescu binterpt[0][0] = -0.8094010767585034; 3791f80e275SEmil Constantinescu binterpt[1][0] = -0.5; 3801f80e275SEmil Constantinescu binterpt[2][0] = 2.3094010767585034; 3811f80e275SEmil Constantinescu binterpt[0][1] = 0.9641016151377548; 3821f80e275SEmil Constantinescu binterpt[1][1] = 0.5; 3831f80e275SEmil Constantinescu binterpt[2][1] = -1.4641016151377548; 384bbd56ea5SKarl Rupp 3851f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 386fe7e6d57SJed Brown } 387fe7e6d57SJed Brown { 3883ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 389da80777bSKarl Rupp /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 390fe7e6d57SJed Brown const PetscReal 391fe7e6d57SJed Brown A[4][4] = {{0,0,0,0}, 392fe7e6d57SJed Brown {8.7173304301691801e-01,0,0,0}, 393fe7e6d57SJed Brown {8.4457060015369423e-01,-1.1299064236484185e-01,0,0}, 394fe7e6d57SJed Brown {0,0,1.,0}}, 395da80777bSKarl Rupp Gamma[4][4] = {{4.3586652150845900e-01,0,0,0}, 396da80777bSKarl Rupp {-8.7173304301691801e-01,4.3586652150845900e-01,0,0}, 397da80777bSKarl Rupp {-9.0338057013044082e-01,5.4180672388095326e-02,4.3586652150845900e-01,0}, 398da80777bSKarl Rupp {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,4.3586652150845900e-01}}, 399fe7e6d57SJed Brown b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01}, 4003ca35412SEmil Constantinescu b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01}; 4013ca35412SEmil Constantinescu 4023ca35412SEmil Constantinescu binterpt[0][0]=1.0564298455794094; 4033ca35412SEmil Constantinescu binterpt[1][0]=2.296429974281067; 4043ca35412SEmil Constantinescu binterpt[2][0]=-1.307599564525376; 4053ca35412SEmil Constantinescu binterpt[3][0]=-1.045260255335102; 4063ca35412SEmil Constantinescu binterpt[0][1]=-1.3864882699759573; 4073ca35412SEmil Constantinescu binterpt[1][1]=-8.262611700275677; 4083ca35412SEmil Constantinescu binterpt[2][1]=7.250979895056055; 4093ca35412SEmil Constantinescu binterpt[3][1]=2.398120075195581; 4103ca35412SEmil Constantinescu binterpt[0][2]=0.5721822314575016; 4113ca35412SEmil Constantinescu binterpt[1][2]=4.742931142090097; 4123ca35412SEmil Constantinescu binterpt[2][2]=-4.398120075195578; 4133ca35412SEmil Constantinescu binterpt[3][2]=-0.9169932983520199; 4143ca35412SEmil Constantinescu 4153ca35412SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 416e27a552bSJed Brown } 417ef3c5b88SJed Brown { 418da80777bSKarl Rupp /* const PetscReal g = 0.5; Directly written in-place below */ 419ef3c5b88SJed Brown const PetscReal 420ef3c5b88SJed Brown A[4][4] = {{0,0,0,0}, 421ef3c5b88SJed Brown {0,0,0,0}, 422ef3c5b88SJed Brown {1.,0,0,0}, 423ef3c5b88SJed Brown {0.75,-0.25,0.5,0}}, 424da80777bSKarl Rupp Gamma[4][4] = {{0.5,0,0,0}, 425da80777bSKarl Rupp {1.,0.5,0,0}, 426da80777bSKarl Rupp {-0.25,-0.25,0.5,0}, 427da80777bSKarl Rupp {1./12,1./12,-2./3,0.5}}, 428ef3c5b88SJed Brown b[4] = {5./6,-1./6,-1./6,0.5}, 429ef3c5b88SJed Brown b2[4] = {0.75,-0.25,0.5,0}; 430bbd56ea5SKarl Rupp 4310298fd71SBarry Smith ierr = TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,NULL);CHKERRQ(ierr); 432ef3c5b88SJed Brown } 433ef3c5b88SJed Brown { 434da80777bSKarl Rupp /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */ 435ef3c5b88SJed Brown const PetscReal 436ef3c5b88SJed Brown A[3][3] = {{0,0,0}, 437da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}, 438da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}}, 439da80777bSKarl Rupp Gamma[3][3] = {{0.43586652150845899941601945119356,0,0}, 440da80777bSKarl Rupp {-0.19294655696029095575009695436041,0.43586652150845899941601945119356,0}, 441da80777bSKarl Rupp {0,1.74927148125794685173529749738960,0.43586652150845899941601945119356}}, 442ef3c5b88SJed Brown b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829}, 443ef3c5b88SJed Brown b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619}; 4441f80e275SEmil Constantinescu 4451f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 4461f80e275SEmil Constantinescu binterpt[0][0] = 3.793692883777660870425141387941; 4471f80e275SEmil Constantinescu binterpt[1][0] = -2.918692883777660870425141387941; 4481f80e275SEmil Constantinescu binterpt[2][0] = 0.125; 4491f80e275SEmil Constantinescu binterpt[0][1] = -0.725741064379812106687651020584; 4501f80e275SEmil Constantinescu binterpt[1][1] = 0.559074397713145440020984353917; 4511f80e275SEmil Constantinescu binterpt[2][1] = 0.16666666666666666666666666666667; 4521f80e275SEmil Constantinescu 4531f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 454ef3c5b88SJed Brown } 455b1c69cc3SEmil Constantinescu { 456da80777bSKarl Rupp /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 457da80777bSKarl Rupp * Direct evaluation: s3 = 1.732050807568877293527; 458da80777bSKarl Rupp * g = 0.7886751345948128822546; 459da80777bSKarl Rupp * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */ 460b1c69cc3SEmil Constantinescu const PetscReal 461b1c69cc3SEmil Constantinescu A[3][3] = {{0,0,0}, 462b1c69cc3SEmil Constantinescu {1,0,0}, 463b1c69cc3SEmil Constantinescu {0.25,0.25,0}}, 464b1c69cc3SEmil Constantinescu Gamma[3][3] = {{0,0,0}, 465da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/6.0,0.7886751345948128822546,0}, 466da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/24.0,(-3.0-1.732050807568877293527)/8.0,0.7886751345948128822546}}, 467b1c69cc3SEmil Constantinescu b[3] = {1./6.,1./6.,2./3.}, 468b1c69cc3SEmil Constantinescu b2[3] = {1./4.,1./4.,1./2.}; 469c0cb691aSEmil Constantinescu PetscReal binterpt[3][2]; 470da80777bSKarl Rupp 471c0cb691aSEmil Constantinescu binterpt[0][0]=0.089316397477040902157517886164709; 472c0cb691aSEmil Constantinescu binterpt[1][0]=-0.91068360252295909784248211383529; 473c0cb691aSEmil Constantinescu binterpt[2][0]=1.8213672050459181956849642276706; 474c0cb691aSEmil Constantinescu binterpt[0][1]=0.077350269189625764509148780501957; 475c0cb691aSEmil Constantinescu binterpt[1][1]=1.077350269189625764509148780502; 476c0cb691aSEmil Constantinescu binterpt[2][1]=-1.1547005383792515290182975610039; 477bbd56ea5SKarl Rupp 478c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 479b1c69cc3SEmil Constantinescu } 480b1c69cc3SEmil Constantinescu 481b1c69cc3SEmil Constantinescu { 482b1c69cc3SEmil Constantinescu const PetscReal 483b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 484b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 485b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 486b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 487b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 488b1c69cc3SEmil Constantinescu {0.0,1./4.,0,0}, 489b1c69cc3SEmil Constantinescu {-2.,-2./3.,2./3.,0}, 490b1c69cc3SEmil Constantinescu {1./2.,5./36.,-2./9,0}}, 491b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 492b1c69cc3SEmil Constantinescu b2[4] = {1./8.,3./4.,1./8.,0}; 493c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 494da80777bSKarl Rupp 495c0cb691aSEmil Constantinescu binterpt[0][0]=6.25; 496c0cb691aSEmil Constantinescu binterpt[1][0]=-30.25; 497c0cb691aSEmil Constantinescu binterpt[2][0]=1.75; 498c0cb691aSEmil Constantinescu binterpt[3][0]=23.25; 499c0cb691aSEmil Constantinescu binterpt[0][1]=-9.75; 500c0cb691aSEmil Constantinescu binterpt[1][1]=58.75; 501c0cb691aSEmil Constantinescu binterpt[2][1]=-3.25; 502c0cb691aSEmil Constantinescu binterpt[3][1]=-45.75; 503c0cb691aSEmil Constantinescu binterpt[0][2]=3.6666666666666666666666666666667; 504c0cb691aSEmil Constantinescu binterpt[1][2]=-28.333333333333333333333333333333; 505c0cb691aSEmil Constantinescu binterpt[2][2]=1.6666666666666666666666666666667; 506c0cb691aSEmil Constantinescu binterpt[3][2]=23.; 507bbd56ea5SKarl Rupp 508c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 509b1c69cc3SEmil Constantinescu } 510b1c69cc3SEmil Constantinescu 511b1c69cc3SEmil Constantinescu { 512b1c69cc3SEmil Constantinescu const PetscReal 513b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 514b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 515b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 516b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 517b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 518b1c69cc3SEmil Constantinescu {0.0,3./4.,0,0}, 519b1c69cc3SEmil Constantinescu {-2./3.,-23./9.,2./9.,0}, 520b1c69cc3SEmil Constantinescu {1./18.,65./108.,-2./27,0}}, 521b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 522b1c69cc3SEmil Constantinescu b2[4] = {3./16.,10./16.,3./16.,0}; 523c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 524da80777bSKarl Rupp 525c0cb691aSEmil Constantinescu binterpt[0][0]=1.6911764705882352941176470588235; 526c0cb691aSEmil Constantinescu binterpt[1][0]=3.6813725490196078431372549019608; 527c0cb691aSEmil Constantinescu binterpt[2][0]=0.23039215686274509803921568627451; 528c0cb691aSEmil Constantinescu binterpt[3][0]=-4.6029411764705882352941176470588; 529c0cb691aSEmil Constantinescu binterpt[0][1]=-0.95588235294117647058823529411765; 530c0cb691aSEmil Constantinescu binterpt[1][1]=-6.2401960784313725490196078431373; 531c0cb691aSEmil Constantinescu binterpt[2][1]=-0.31862745098039215686274509803922; 532c0cb691aSEmil Constantinescu binterpt[3][1]=7.5147058823529411764705882352941; 533c0cb691aSEmil Constantinescu binterpt[0][2]=-0.56862745098039215686274509803922; 534c0cb691aSEmil Constantinescu binterpt[1][2]=2.7254901960784313725490196078431; 535c0cb691aSEmil Constantinescu binterpt[2][2]=0.25490196078431372549019607843137; 536c0cb691aSEmil Constantinescu binterpt[3][2]=-2.4117647058823529411764705882353; 537bbd56ea5SKarl Rupp 538c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 539b1c69cc3SEmil Constantinescu } 540753f8adbSEmil Constantinescu 541753f8adbSEmil Constantinescu { 542753f8adbSEmil Constantinescu PetscReal A[4][4],Gamma[4][4],b[4],b2[4]; 5433ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 544753f8adbSEmil Constantinescu 545753f8adbSEmil Constantinescu Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816; 54605e8e825SJed Brown Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0; 547753f8adbSEmil Constantinescu Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476; 548753f8adbSEmil Constantinescu Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816; 54905e8e825SJed Brown Gamma[1][2]=0; Gamma[1][3]=0; 550753f8adbSEmil Constantinescu Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903; 551753f8adbSEmil Constantinescu Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131; 552753f8adbSEmil Constantinescu Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816; 55305e8e825SJed Brown Gamma[2][3]=0; 554753f8adbSEmil Constantinescu Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783; 555753f8adbSEmil Constantinescu Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984; 556753f8adbSEmil Constantinescu Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198; 557753f8adbSEmil Constantinescu Gamma[3][3]=0; 558753f8adbSEmil Constantinescu 55905e8e825SJed Brown A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0; 560753f8adbSEmil Constantinescu A[1][0]=0.8717330430169179988320388950590125027645343373957631; 56105e8e825SJed Brown A[1][1]=0; A[1][2]=0; A[1][3]=0; 562753f8adbSEmil Constantinescu A[2][0]=0.5275890119763004115618079766722914408876108660811028; 563753f8adbSEmil Constantinescu A[2][1]=0.07241098802369958843819203208518599088698057726988732; 56405e8e825SJed Brown A[2][2]=0; A[2][3]=0; 565753f8adbSEmil Constantinescu A[3][0]=0.3990960076760701320627260685975778145384666450351314; 566753f8adbSEmil Constantinescu A[3][1]=-0.4375576546135194437228463747348862825846903771419953; 567753f8adbSEmil Constantinescu A[3][2]=1.038461646937449311660120300601880176655352737312713; 56805e8e825SJed Brown A[3][3]=0; 569753f8adbSEmil Constantinescu 570753f8adbSEmil Constantinescu b[0]=0.1876410243467238251612921333138006734899663569186926; 571753f8adbSEmil Constantinescu b[1]=-0.5952974735769549480478230473706443582188442040780541; 572753f8adbSEmil Constantinescu b[2]=0.9717899277217721234705114616271378792182450260943198; 573753f8adbSEmil Constantinescu b[3]=0.4358665215084589994160194475295062513822671686978816; 574753f8adbSEmil Constantinescu 575753f8adbSEmil Constantinescu b2[0]=0.2147402862233891404862383521089097657790734483804460; 576753f8adbSEmil Constantinescu b2[1]=-0.4851622638849390928209050538171743017757490232519684; 577753f8adbSEmil Constantinescu b2[2]=0.8687250025203875511662123688667549217531982787600080; 578753f8adbSEmil Constantinescu b2[3]=0.4016969751411624011684543450940068201770721128357014; 579753f8adbSEmil Constantinescu 5803ca35412SEmil Constantinescu binterpt[0][0]=2.2565812720167954547104627844105; 5813ca35412SEmil Constantinescu binterpt[1][0]=1.349166413351089573796243820819; 5823ca35412SEmil Constantinescu binterpt[2][0]=-2.4695174540533503758652847586647; 5833ca35412SEmil Constantinescu binterpt[3][0]=-0.13623023131453465264142184656474; 5843ca35412SEmil Constantinescu binterpt[0][1]=-3.0826699111559187902922463354557; 5853ca35412SEmil Constantinescu binterpt[1][1]=-2.4689115685996042534544925650515; 5863ca35412SEmil Constantinescu binterpt[2][1]=5.7428279814696677152129332773553; 5873ca35412SEmil Constantinescu binterpt[3][1]=-0.19124650171414467146619437684812; 5883ca35412SEmil Constantinescu binterpt[0][2]=1.0137296634858471607430756831148; 5893ca35412SEmil Constantinescu binterpt[1][2]=0.52444768167155973161042570784064; 5903ca35412SEmil Constantinescu binterpt[2][2]=-2.3015205996945452158771370439586; 5913ca35412SEmil Constantinescu binterpt[3][2]=0.76334325453713832352363565300308; 592f4aed992SEmil Constantinescu 593f73f8d2cSSatish Balay ierr = TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 594753f8adbSEmil Constantinescu } 59542faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWGRK4T,0.231,PETSC_DEFAULT,PETSC_DEFAULT,0,-0.1282612945269037e+01);CHKERRQ(ierr); 59642faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWSHAMP4,0.5,PETSC_DEFAULT,PETSC_DEFAULT,0,125./108.);CHKERRQ(ierr); 59742faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWVELDD4,0.22570811482256823492,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.355958941201148);CHKERRQ(ierr); 59842faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSW4L,0.57282,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.093502252409163);CHKERRQ(ierr); 599e27a552bSJed Brown PetscFunctionReturn(0); 600e27a552bSJed Brown } 601e27a552bSJed Brown 602d5e6173cSPeter Brune 603d5e6173cSPeter Brune 604e27a552bSJed Brown #undef __FUNCT__ 605e27a552bSJed Brown #define __FUNCT__ "TSRosWRegisterDestroy" 606e27a552bSJed Brown /*@C 607e27a552bSJed Brown TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 608e27a552bSJed Brown 609e27a552bSJed Brown Not Collective 610e27a552bSJed Brown 611e27a552bSJed Brown Level: advanced 612e27a552bSJed Brown 613e27a552bSJed Brown .keywords: TSRosW, register, destroy 614e27a552bSJed Brown .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic() 615e27a552bSJed Brown @*/ 616e27a552bSJed Brown PetscErrorCode TSRosWRegisterDestroy(void) 617e27a552bSJed Brown { 618e27a552bSJed Brown PetscErrorCode ierr; 61961692a83SJed Brown RosWTableauLink link; 620e27a552bSJed Brown 621e27a552bSJed Brown PetscFunctionBegin; 62261692a83SJed Brown while ((link = RosWTableauList)) { 62361692a83SJed Brown RosWTableau t = &link->tab; 62461692a83SJed Brown RosWTableauList = link->next; 62561692a83SJed Brown ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 62643b21953SEmil Constantinescu ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 627fe7e6d57SJed Brown ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 628f4aed992SEmil Constantinescu ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 629e27a552bSJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 630e27a552bSJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 631e27a552bSJed Brown } 632e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 633e27a552bSJed Brown PetscFunctionReturn(0); 634e27a552bSJed Brown } 635e27a552bSJed Brown 636e27a552bSJed Brown #undef __FUNCT__ 637e27a552bSJed Brown #define __FUNCT__ "TSRosWInitializePackage" 638e27a552bSJed Brown /*@C 639e27a552bSJed Brown TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 640e27a552bSJed Brown from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW() 641e27a552bSJed Brown when using static libraries. 642e27a552bSJed Brown 643e27a552bSJed Brown Input Parameter: 6440298fd71SBarry Smith path - The dynamic library path, or NULL 645e27a552bSJed Brown 646e27a552bSJed Brown Level: developer 647e27a552bSJed Brown 648e27a552bSJed Brown .keywords: TS, TSRosW, initialize, package 649e27a552bSJed Brown .seealso: PetscInitialize() 650e27a552bSJed Brown @*/ 651e27a552bSJed Brown PetscErrorCode TSRosWInitializePackage(const char path[]) 652e27a552bSJed Brown { 653e27a552bSJed Brown PetscErrorCode ierr; 654e27a552bSJed Brown 655e27a552bSJed Brown PetscFunctionBegin; 656e27a552bSJed Brown if (TSRosWPackageInitialized) PetscFunctionReturn(0); 657e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 658e27a552bSJed Brown ierr = TSRosWRegisterAll();CHKERRQ(ierr); 659e27a552bSJed Brown ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 660e27a552bSJed Brown PetscFunctionReturn(0); 661e27a552bSJed Brown } 662e27a552bSJed Brown 663e27a552bSJed Brown #undef __FUNCT__ 664e27a552bSJed Brown #define __FUNCT__ "TSRosWFinalizePackage" 665e27a552bSJed Brown /*@C 666e27a552bSJed Brown TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 667e27a552bSJed Brown called from PetscFinalize(). 668e27a552bSJed Brown 669e27a552bSJed Brown Level: developer 670e27a552bSJed Brown 671e27a552bSJed Brown .keywords: Petsc, destroy, package 672e27a552bSJed Brown .seealso: PetscFinalize() 673e27a552bSJed Brown @*/ 674e27a552bSJed Brown PetscErrorCode TSRosWFinalizePackage(void) 675e27a552bSJed Brown { 676e27a552bSJed Brown PetscErrorCode ierr; 677e27a552bSJed Brown 678e27a552bSJed Brown PetscFunctionBegin; 679e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 680e27a552bSJed Brown ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 681e27a552bSJed Brown PetscFunctionReturn(0); 682e27a552bSJed Brown } 683e27a552bSJed Brown 684e27a552bSJed Brown #undef __FUNCT__ 685e27a552bSJed Brown #define __FUNCT__ "TSRosWRegister" 686e27a552bSJed Brown /*@C 68761692a83SJed Brown TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 688e27a552bSJed Brown 689e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 690e27a552bSJed Brown 691e27a552bSJed Brown Input Parameters: 692e27a552bSJed Brown + name - identifier for method 693e27a552bSJed Brown . order - approximation order of method 694e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 69561692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 69661692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 697fe7e6d57SJed Brown . b - Step completion table (dimension s) 6980298fd71SBarry Smith . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available) 699f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 70042faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 701e27a552bSJed Brown 702e27a552bSJed Brown Notes: 70361692a83SJed Brown Several Rosenbrock W methods are provided, this function is only needed to create new methods. 704e27a552bSJed Brown 705e27a552bSJed Brown Level: advanced 706e27a552bSJed Brown 707e27a552bSJed Brown .keywords: TS, register 708e27a552bSJed Brown 709e27a552bSJed Brown .seealso: TSRosW 710e27a552bSJed Brown @*/ 711f9c1d6abSBarry Smith PetscErrorCode TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 712f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 713e27a552bSJed Brown { 714e27a552bSJed Brown PetscErrorCode ierr; 71561692a83SJed Brown RosWTableauLink link; 71661692a83SJed Brown RosWTableau t; 71761692a83SJed Brown PetscInt i,j,k; 71861692a83SJed Brown PetscScalar *GammaInv; 719e27a552bSJed Brown 720e27a552bSJed Brown PetscFunctionBegin; 721fe7e6d57SJed Brown PetscValidCharPointer(name,1); 722fe7e6d57SJed Brown PetscValidPointer(A,4); 723fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 724fe7e6d57SJed Brown PetscValidPointer(b,6); 725fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 726fe7e6d57SJed Brown 727e27a552bSJed Brown ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 728e27a552bSJed Brown ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 729e27a552bSJed Brown t = &link->tab; 730e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 731e27a552bSJed Brown t->order = order; 732e27a552bSJed Brown t->s = s; 73361692a83SJed 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); 73443b21953SEmil 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); 735e27a552bSJed Brown ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 73661692a83SJed Brown ierr = PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 73743b21953SEmil Constantinescu ierr = PetscMemcpy(t->GammaExplicitCorr,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 73861692a83SJed Brown ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 739fe7e6d57SJed Brown if (bembed) { 740fe7e6d57SJed Brown ierr = PetscMalloc2(s,PetscReal,&t->bembed,s,PetscReal,&t->bembedt);CHKERRQ(ierr); 741fe7e6d57SJed Brown ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 742fe7e6d57SJed Brown } 74361692a83SJed Brown for (i=0; i<s; i++) { 74461692a83SJed Brown t->ASum[i] = 0; 74561692a83SJed Brown t->GammaSum[i] = 0; 74661692a83SJed Brown for (j=0; j<s; j++) { 74761692a83SJed Brown t->ASum[i] += A[i*s+j]; 748fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 74961692a83SJed Brown } 75061692a83SJed Brown } 75161692a83SJed Brown ierr = PetscMalloc(s*s*sizeof(PetscScalar),&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 75261692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 753fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 754fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 755fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 756c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 757fd96d5b0SEmil Constantinescu } else { 758c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 759fd96d5b0SEmil Constantinescu } 760fd96d5b0SEmil Constantinescu } 761fd96d5b0SEmil Constantinescu 76261692a83SJed Brown switch (s) { 76361692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 76496b95a6bSBarry Smith case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0);CHKERRQ(ierr); break; 76596b95a6bSBarry Smith case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0);CHKERRQ(ierr); break; 76696b95a6bSBarry Smith case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0);CHKERRQ(ierr); break; 76761692a83SJed Brown case 5: { 76861692a83SJed Brown PetscInt ipvt5[5]; 76961692a83SJed Brown MatScalar work5[5*5]; 77096b95a6bSBarry Smith ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0);CHKERRQ(ierr); break; 77161692a83SJed Brown } 77296b95a6bSBarry Smith case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0);CHKERRQ(ierr); break; 77396b95a6bSBarry Smith case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0);CHKERRQ(ierr); break; 77461692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 77561692a83SJed Brown } 77661692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 77761692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 77843b21953SEmil Constantinescu 77943b21953SEmil Constantinescu for (i=0; i<s; i++) { 78043b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 78143b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 78243b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 78343b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 78443b21953SEmil Constantinescu } 78543b21953SEmil Constantinescu } 78643b21953SEmil Constantinescu } 78743b21953SEmil Constantinescu 78861692a83SJed Brown for (i=0; i<s; i++) { 78961692a83SJed Brown for (j=0; j<s; j++) { 79061692a83SJed Brown t->At[i*s+j] = 0; 79161692a83SJed Brown for (k=0; k<s; k++) { 79261692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 79361692a83SJed Brown } 79461692a83SJed Brown } 79561692a83SJed Brown t->bt[i] = 0; 79661692a83SJed Brown for (j=0; j<s; j++) { 79761692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 79861692a83SJed Brown } 799fe7e6d57SJed Brown if (bembed) { 800fe7e6d57SJed Brown t->bembedt[i] = 0; 801fe7e6d57SJed Brown for (j=0; j<s; j++) { 802fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 803fe7e6d57SJed Brown } 804fe7e6d57SJed Brown } 80561692a83SJed Brown } 8068d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 8078d59e960SJed Brown 808f4aed992SEmil Constantinescu t->pinterp = pinterp; 8093ca35412SEmil Constantinescu ierr = PetscMalloc(s*pinterp*sizeof(binterpt[0]),&t->binterpt);CHKERRQ(ierr); 8103ca35412SEmil Constantinescu ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 81161692a83SJed Brown link->next = RosWTableauList; 81261692a83SJed Brown RosWTableauList = link; 813e27a552bSJed Brown PetscFunctionReturn(0); 814e27a552bSJed Brown } 815e27a552bSJed Brown 816e27a552bSJed Brown #undef __FUNCT__ 81742faf41dSJed Brown #define __FUNCT__ "TSRosWRegisterRos4" 81842faf41dSJed Brown /*@C 81942faf41dSJed Brown TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing paramter choices 82042faf41dSJed Brown 82142faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 82242faf41dSJed Brown 82342faf41dSJed Brown Input Parameters: 82442faf41dSJed Brown + name - identifier for method 82542faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 82642faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 82742faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 82842faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 82942faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 83042faf41dSJed Brown . e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 83142faf41dSJed Brown 83242faf41dSJed Brown Notes: 83342faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 83442faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 83542faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 83642faf41dSJed Brown 83742faf41dSJed Brown Level: developer 83842faf41dSJed Brown 83942faf41dSJed Brown .keywords: TS, register 84042faf41dSJed Brown 84142faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 84242faf41dSJed Brown @*/ 84319fd82e9SBarry Smith PetscErrorCode TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 84442faf41dSJed Brown { 84542faf41dSJed Brown PetscErrorCode ierr; 84642faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 84742faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 84842faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 84942faf41dSJed Brown p42 = one/eight - gamma/three, 85042faf41dSJed Brown p43 = one/twelve - gamma/three, 85142faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 85242faf41dSJed Brown p56 = one/twenty - gamma/four; 85342faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 85442faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 85542faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 85642faf41dSJed Brown 85742faf41dSJed Brown PetscFunctionBegin; 85842faf41dSJed Brown /* Step 1: choose Gamma (input) */ 85942faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 86042faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 86142faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 86242faf41dSJed Brown 86342faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 86442faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 86542faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 86642faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 86742faf41dSJed Brown rhs[0] = one - b3; 86842faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 86942faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 87042faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 87142faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 87242faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 87342faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 87442faf41dSJed Brown 87542faf41dSJed Brown /* Step 3 */ 87642faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 87742faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 87842faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 87942faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 88042faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 88142faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 88242faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 88342faf41dSJed Brown ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0);CHKERRQ(ierr); 88442faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 88542faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 88642faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 88742faf41dSJed Brown 88842faf41dSJed Brown /* Step 4: back-substitute */ 88942faf41dSJed Brown beta32 = beta32beta2p / beta2p; 89042faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 89142faf41dSJed Brown 89242faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 89342faf41dSJed Brown a43 = 0; 89442faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 89542faf41dSJed Brown a42 = a32; 89642faf41dSJed Brown 89742faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 89842faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 89942faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 90042faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 90142faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 90242faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 90342faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 90442faf41dSJed 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; 90542faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 90642faf41dSJed Brown 90742faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 90842faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 90942faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 91042faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 91142faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 91242faf41dSJed Brown 91342faf41dSJed Brown { 91442faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 91542faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 91642faf41dSJed Brown } 9170298fd71SBarry Smith ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,NULL);CHKERRQ(ierr); 91842faf41dSJed Brown PetscFunctionReturn(0); 91942faf41dSJed Brown } 92042faf41dSJed Brown 92142faf41dSJed Brown #undef __FUNCT__ 9221c3436cfSJed Brown #define __FUNCT__ "TSEvaluateStep_RosW" 9231c3436cfSJed Brown /* 9241c3436cfSJed Brown The step completion formula is 9251c3436cfSJed Brown 9261c3436cfSJed Brown x1 = x0 + b^T Y 9271c3436cfSJed Brown 9281c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 9291c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 9301c3436cfSJed Brown 9311c3436cfSJed Brown x1e = x0 + be^T Y 9321c3436cfSJed Brown = x1 - b^T Y + be^T Y 9331c3436cfSJed Brown = x1 + (be - b)^T Y 9341c3436cfSJed Brown 9351c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9361c3436cfSJed Brown */ 937f9c1d6abSBarry Smith static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done) 9381c3436cfSJed Brown { 9391c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9401c3436cfSJed Brown RosWTableau tab = ros->tableau; 9411c3436cfSJed Brown PetscScalar *w = ros->work; 9421c3436cfSJed Brown PetscInt i; 9431c3436cfSJed Brown PetscErrorCode ierr; 9441c3436cfSJed Brown 9451c3436cfSJed Brown PetscFunctionBegin; 9461c3436cfSJed Brown if (order == tab->order) { 947108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 948f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 949de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 950f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 951f9c1d6abSBarry Smith } else {ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr);} 9521c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9531c3436cfSJed Brown PetscFunctionReturn(0); 9541c3436cfSJed Brown } else if (order == tab->order-1) { 9551c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 956108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 957f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 958de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 959f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 960108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 961108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 962f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 963f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 9641c3436cfSJed Brown } 9651c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9661c3436cfSJed Brown PetscFunctionReturn(0); 9671c3436cfSJed Brown } 9681c3436cfSJed Brown unavailable: 9691c3436cfSJed Brown if (done) *done = PETSC_FALSE; 970ce94432eSBarry Smith else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 9711c3436cfSJed Brown PetscFunctionReturn(0); 9721c3436cfSJed Brown } 9731c3436cfSJed Brown 9741c3436cfSJed Brown #undef __FUNCT__ 975e27a552bSJed Brown #define __FUNCT__ "TSStep_RosW" 976e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 977e27a552bSJed Brown { 97861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 97961692a83SJed Brown RosWTableau tab = ros->tableau; 980e27a552bSJed Brown const PetscInt s = tab->s; 9811c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9820feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 983c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 98461692a83SJed Brown PetscScalar *w = ros->work; 9857d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 986e27a552bSJed Brown SNES snes; 9871c3436cfSJed Brown TSAdapt adapt; 9881c3436cfSJed Brown PetscInt i,j,its,lits,reject,next_scheme; 989cdbf8f93SLisandro Dalcin PetscReal next_time_step; 9901c3436cfSJed Brown PetscBool accept; 991e27a552bSJed Brown PetscErrorCode ierr; 9920feba352SEmil Constantinescu MatStructure str; 993e27a552bSJed Brown 994e27a552bSJed Brown PetscFunctionBegin; 995e27a552bSJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 996cdbf8f93SLisandro Dalcin next_time_step = ts->time_step; 9971c3436cfSJed Brown accept = PETSC_TRUE; 998108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 999e27a552bSJed Brown 100097335746SJed Brown for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 10011c3436cfSJed Brown const PetscReal h = ts->time_step; 1002b8123daeSJed Brown ierr = TSPreStep(ts);CHKERRQ(ierr); 10033ca35412SEmil Constantinescu ierr = VecCopy(ts->vec_sol,ros->VecSolPrev);CHKERRQ(ierr); /*move this at the end*/ 1004e27a552bSJed Brown for (i=0; i<s; i++) { 10051c3436cfSJed Brown ros->stage_time = ts->ptime + h*ASum[i]; 1006b8123daeSJed Brown ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 1007c17803e7SJed Brown if (GammaZeroDiag[i]) { 1008c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 1009b296d7d5SJed Brown ros->scoeff = 1.; 1010c17803e7SJed Brown } else { 1011c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 1012b296d7d5SJed Brown ros->scoeff = 1./Gamma[i*s+i]; 1013fd96d5b0SEmil Constantinescu } 101461692a83SJed Brown 101561692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 1016de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 1017de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 101861692a83SJed Brown 101961692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 102061692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 102161692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 102261692a83SJed Brown 1023e27a552bSJed Brown /* Initial guess taken from last stage */ 102461692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 102561692a83SJed Brown 10267d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 102761692a83SJed Brown if (!ros->recompute_jacobian && !i) { 102861692a83SJed Brown ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 102961692a83SJed Brown } 10300298fd71SBarry Smith ierr = SNESSolve(snes,NULL,Y[i]);CHKERRQ(ierr); 1031e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1032e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10335ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 1034ad6bc421SBarry Smith ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 103597335746SJed Brown ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 103697335746SJed Brown if (!accept) goto reject_step; 10377d4bf2deSEmil Constantinescu } else { 10381ce71dffSSatish Balay Mat J,Jp; 10390feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10400feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 104122d28d08SBarry Smith ierr = VecScale(Y[i],-1.0);CHKERRQ(ierr); 10420feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i]=F(Zstage)-Zdot[=GammaInv*Y]*/ 10430feba352SEmil Constantinescu 10440feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10450feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10460feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 10470feba352SEmil Constantinescu /*Y[i] += Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10480feba352SEmil Constantinescu str = SAME_NONZERO_PATTERN; 10490298fd71SBarry Smith ierr = TSGetIJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 10500feba352SEmil Constantinescu ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,&J,&Jp,&str,PETSC_FALSE);CHKERRQ(ierr); 105122d28d08SBarry Smith ierr = MatMult(J,Zstage,Zdot);CHKERRQ(ierr); 10520feba352SEmil Constantinescu 10530feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10540feba352SEmil Constantinescu ierr = VecScale(Y[i],h); 10555ef26d82SJed Brown ts->ksp_its += 1; 10567d4bf2deSEmil Constantinescu } 1057e27a552bSJed Brown } 10580298fd71SBarry Smith ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 1059108c343cSJed Brown ros->status = TS_STEP_PENDING; 1060e27a552bSJed Brown 10611c3436cfSJed Brown /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 1062ad6bc421SBarry Smith ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 10631c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10648d59e960SJed Brown ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 10651c3436cfSJed Brown ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 10661c3436cfSJed Brown if (accept) { 10671c3436cfSJed Brown /* ignore next_scheme for now */ 1068e27a552bSJed Brown ts->ptime += ts->time_step; 1069cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 1070e27a552bSJed Brown ts->steps++; 1071108c343cSJed Brown ros->status = TS_STEP_COMPLETE; 10721c3436cfSJed Brown break; 10731c3436cfSJed Brown } else { /* Roll back the current step */ 10741c3436cfSJed Brown for (i=0; i<s; i++) w[i] = -tab->bt[i]; 10751c3436cfSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 10761c3436cfSJed Brown ts->time_step = next_time_step; 1077108c343cSJed Brown ros->status = TS_STEP_INCOMPLETE; 10781c3436cfSJed Brown } 1079476b6736SJed Brown reject_step: continue; 10801c3436cfSJed Brown } 1081b2ce242eSJed Brown if (ros->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 1082e27a552bSJed Brown PetscFunctionReturn(0); 1083e27a552bSJed Brown } 1084e27a552bSJed Brown 1085e27a552bSJed Brown #undef __FUNCT__ 1086e27a552bSJed Brown #define __FUNCT__ "TSInterpolate_RosW" 1087f9c1d6abSBarry Smith static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1088e27a552bSJed Brown { 108961692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1090f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1091f4aed992SEmil Constantinescu PetscReal h; 1092f4aed992SEmil Constantinescu PetscReal tt,t; 1093f4aed992SEmil Constantinescu PetscScalar *bt; 1094f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1095f4aed992SEmil Constantinescu PetscErrorCode ierr; 1096f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1097f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1098f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1099e27a552bSJed Brown 1100e27a552bSJed Brown PetscFunctionBegin; 1101ce94432eSBarry Smith if (!Bt) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1102f4aed992SEmil Constantinescu 1103f4aed992SEmil Constantinescu switch (ros->status) { 1104f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1105f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1106f4aed992SEmil Constantinescu h = ts->time_step; 1107f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1108f4aed992SEmil Constantinescu break; 1109f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1110f4aed992SEmil Constantinescu h = ts->time_step_prev; 1111f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1112f4aed992SEmil Constantinescu break; 1113ce94432eSBarry Smith default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1114f4aed992SEmil Constantinescu } 11153ca35412SEmil Constantinescu ierr = PetscMalloc(s*sizeof(bt[0]),&bt);CHKERRQ(ierr); 1116f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1117f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1118f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 11193ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1120f4aed992SEmil Constantinescu } 1121f4aed992SEmil Constantinescu } 1122f4aed992SEmil Constantinescu 1123f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1124f9c1d6abSBarry Smith /*U<-0*/ 1125f9c1d6abSBarry Smith ierr = VecZeroEntries(U);CHKERRQ(ierr); 1126f4aed992SEmil Constantinescu 1127f9c1d6abSBarry Smith /*U<- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11283ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j]=0; 11293ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11303ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11313ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1132f4aed992SEmil Constantinescu } 11333ca35412SEmil Constantinescu } 1134f9c1d6abSBarry Smith ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1135f4aed992SEmil Constantinescu 1136f4aed992SEmil Constantinescu /*X<-y(t) + X*/ 1137f9c1d6abSBarry Smith ierr = VecAXPY(U,1.0,ros->VecSolPrev);CHKERRQ(ierr); 1138f4aed992SEmil Constantinescu 1139f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1140e27a552bSJed Brown PetscFunctionReturn(0); 1141e27a552bSJed Brown } 1142e27a552bSJed Brown 1143e27a552bSJed Brown /*------------------------------------------------------------*/ 1144e27a552bSJed Brown #undef __FUNCT__ 1145e27a552bSJed Brown #define __FUNCT__ "TSReset_RosW" 1146e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1147e27a552bSJed Brown { 114861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1149e27a552bSJed Brown PetscInt s; 1150e27a552bSJed Brown PetscErrorCode ierr; 1151e27a552bSJed Brown 1152e27a552bSJed Brown PetscFunctionBegin; 115361692a83SJed Brown if (!ros->tableau) PetscFunctionReturn(0); 115461692a83SJed Brown s = ros->tableau->s; 115561692a83SJed Brown ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 115661692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 115761692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 115861692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 115961692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 11603ca35412SEmil Constantinescu ierr = VecDestroy(&ros->VecSolPrev);CHKERRQ(ierr); 116161692a83SJed Brown ierr = PetscFree(ros->work);CHKERRQ(ierr); 1162e27a552bSJed Brown PetscFunctionReturn(0); 1163e27a552bSJed Brown } 1164e27a552bSJed Brown 1165e27a552bSJed Brown #undef __FUNCT__ 1166e27a552bSJed Brown #define __FUNCT__ "TSDestroy_RosW" 1167e27a552bSJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1168e27a552bSJed Brown { 1169e27a552bSJed Brown PetscErrorCode ierr; 1170e27a552bSJed Brown 1171e27a552bSJed Brown PetscFunctionBegin; 1172e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1173e27a552bSJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 117400de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C","",NULL);CHKERRQ(ierr); 117500de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C","",NULL);CHKERRQ(ierr); 117600de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",NULL);CHKERRQ(ierr); 1177e27a552bSJed Brown PetscFunctionReturn(0); 1178e27a552bSJed Brown } 1179e27a552bSJed Brown 1180d5e6173cSPeter Brune 1181d5e6173cSPeter Brune #undef __FUNCT__ 1182d5e6173cSPeter Brune #define __FUNCT__ "TSRosWGetVecs" 1183d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1184d5e6173cSPeter Brune { 1185d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1186d5e6173cSPeter Brune PetscErrorCode ierr; 1187d5e6173cSPeter Brune 1188d5e6173cSPeter Brune PetscFunctionBegin; 1189d5e6173cSPeter Brune if (Ydot) { 1190d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1191d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1192d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1193d5e6173cSPeter Brune } 1194d5e6173cSPeter Brune if (Zdot) { 1195d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1196d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1197d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1198d5e6173cSPeter Brune } 1199d5e6173cSPeter Brune if (Ystage) { 1200d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1201d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1202d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1203d5e6173cSPeter Brune } 1204d5e6173cSPeter Brune if (Zstage) { 1205d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1206d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1207d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1208d5e6173cSPeter Brune } 1209d5e6173cSPeter Brune PetscFunctionReturn(0); 1210d5e6173cSPeter Brune } 1211d5e6173cSPeter Brune 1212d5e6173cSPeter Brune 1213d5e6173cSPeter Brune #undef __FUNCT__ 1214d5e6173cSPeter Brune #define __FUNCT__ "TSRosWRestoreVecs" 1215d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1216d5e6173cSPeter Brune { 1217d5e6173cSPeter Brune PetscErrorCode ierr; 1218d5e6173cSPeter Brune 1219d5e6173cSPeter Brune PetscFunctionBegin; 1220d5e6173cSPeter Brune if (Ydot) { 1221d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1222d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1223d5e6173cSPeter Brune } 1224d5e6173cSPeter Brune } 1225d5e6173cSPeter Brune if (Zdot) { 1226d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1227d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1228d5e6173cSPeter Brune } 1229d5e6173cSPeter Brune } 1230d5e6173cSPeter Brune if (Ystage) { 1231d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1232d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1233d5e6173cSPeter Brune } 1234d5e6173cSPeter Brune } 1235d5e6173cSPeter Brune if (Zstage) { 1236d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1237d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1238d5e6173cSPeter Brune } 1239d5e6173cSPeter Brune } 1240d5e6173cSPeter Brune PetscFunctionReturn(0); 1241d5e6173cSPeter Brune } 1242d5e6173cSPeter Brune 1243d5e6173cSPeter Brune #undef __FUNCT__ 1244d5e6173cSPeter Brune #define __FUNCT__ "DMCoarsenHook_TSRosW" 1245d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1246d5e6173cSPeter Brune { 1247d5e6173cSPeter Brune PetscFunctionBegin; 1248d5e6173cSPeter Brune PetscFunctionReturn(0); 1249d5e6173cSPeter Brune } 1250d5e6173cSPeter Brune 1251d5e6173cSPeter Brune #undef __FUNCT__ 1252d5e6173cSPeter Brune #define __FUNCT__ "DMRestrictHook_TSRosW" 1253d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1254d5e6173cSPeter Brune { 1255d5e6173cSPeter Brune TS ts = (TS)ctx; 1256d5e6173cSPeter Brune PetscErrorCode ierr; 1257d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1258d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1259d5e6173cSPeter Brune 1260d5e6173cSPeter Brune PetscFunctionBegin; 1261d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1262d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1263d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1264d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1265d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1266d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1267d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1268d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1269d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1270d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1271d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1272d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1273d5e6173cSPeter Brune PetscFunctionReturn(0); 1274d5e6173cSPeter Brune } 1275d5e6173cSPeter Brune 1276258e1594SPeter Brune 1277258e1594SPeter Brune #undef __FUNCT__ 1278258e1594SPeter Brune #define __FUNCT__ "DMSubDomainHook_TSRosW" 1279258e1594SPeter Brune static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx) 1280258e1594SPeter Brune { 1281258e1594SPeter Brune PetscFunctionBegin; 1282258e1594SPeter Brune PetscFunctionReturn(0); 1283258e1594SPeter Brune } 1284258e1594SPeter Brune 1285258e1594SPeter Brune #undef __FUNCT__ 1286258e1594SPeter Brune #define __FUNCT__ "DMSubDomainRestrictHook_TSRosW" 1287258e1594SPeter Brune static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1288258e1594SPeter Brune { 1289258e1594SPeter Brune TS ts = (TS)ctx; 1290258e1594SPeter Brune PetscErrorCode ierr; 1291258e1594SPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1292258e1594SPeter Brune Vec Ydots,Zdots,Ystages,Zstages; 1293258e1594SPeter Brune 1294258e1594SPeter Brune PetscFunctionBegin; 1295258e1594SPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1296258e1594SPeter Brune ierr = TSRosWGetVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1297258e1594SPeter Brune 1298258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1299258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1300258e1594SPeter Brune 1301258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1302258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1303258e1594SPeter Brune 1304258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1305258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1306258e1594SPeter Brune 1307258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1308258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1309258e1594SPeter Brune 1310258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1311258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1312258e1594SPeter Brune PetscFunctionReturn(0); 1313258e1594SPeter Brune } 1314258e1594SPeter Brune 1315e27a552bSJed Brown /* 1316e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1317e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1318e27a552bSJed Brown */ 1319e27a552bSJed Brown #undef __FUNCT__ 1320e27a552bSJed Brown #define __FUNCT__ "SNESTSFormFunction_RosW" 1321f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1322e27a552bSJed Brown { 132361692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1324e27a552bSJed Brown PetscErrorCode ierr; 1325d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1326b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1327d5e6173cSPeter Brune DM dm,dmsave; 1328e27a552bSJed Brown 1329e27a552bSJed Brown PetscFunctionBegin; 1330d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1331d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1332b296d7d5SJed Brown ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1333f9c1d6abSBarry Smith ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1334d5e6173cSPeter Brune dmsave = ts->dm; 1335d5e6173cSPeter Brune ts->dm = dm; 1336d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1337d5e6173cSPeter Brune ts->dm = dmsave; 1338d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1339e27a552bSJed Brown PetscFunctionReturn(0); 1340e27a552bSJed Brown } 1341e27a552bSJed Brown 1342e27a552bSJed Brown #undef __FUNCT__ 1343e27a552bSJed Brown #define __FUNCT__ "SNESTSFormJacobian_RosW" 1344f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *str,TS ts) 1345e27a552bSJed Brown { 134661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1347d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1348b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1349e27a552bSJed Brown PetscErrorCode ierr; 1350d5e6173cSPeter Brune DM dm,dmsave; 1351e27a552bSJed Brown 1352e27a552bSJed Brown PetscFunctionBegin; 135361692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1354d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1355d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1356d5e6173cSPeter Brune dmsave = ts->dm; 1357d5e6173cSPeter Brune ts->dm = dm; 1358b296d7d5SJed Brown ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 1359d5e6173cSPeter Brune ts->dm = dmsave; 1360d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1361e27a552bSJed Brown PetscFunctionReturn(0); 1362e27a552bSJed Brown } 1363e27a552bSJed Brown 1364e27a552bSJed Brown #undef __FUNCT__ 1365e27a552bSJed Brown #define __FUNCT__ "TSSetUp_RosW" 1366e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1367e27a552bSJed Brown { 136861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 136961692a83SJed Brown RosWTableau tab = ros->tableau; 1370e27a552bSJed Brown PetscInt s = tab->s; 1371e27a552bSJed Brown PetscErrorCode ierr; 1372d5e6173cSPeter Brune DM dm; 1373e27a552bSJed Brown 1374e27a552bSJed Brown PetscFunctionBegin; 137561692a83SJed Brown if (!ros->tableau) { 1376e27a552bSJed Brown ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1377e27a552bSJed Brown } 137861692a83SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 137961692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 138061692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 138161692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 138261692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 13833ca35412SEmil Constantinescu ierr = VecDuplicate(ts->vec_sol,&ros->VecSolPrev);CHKERRQ(ierr); 138461692a83SJed Brown ierr = PetscMalloc(s*sizeof(ros->work[0]),&ros->work);CHKERRQ(ierr); 138522d28d08SBarry Smith ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1386d5e6173cSPeter Brune if (dm) { 1387d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1388258e1594SPeter Brune ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1389d5e6173cSPeter Brune } 1390e27a552bSJed Brown PetscFunctionReturn(0); 1391e27a552bSJed Brown } 1392e27a552bSJed Brown /*------------------------------------------------------------*/ 1393e27a552bSJed Brown 1394e27a552bSJed Brown #undef __FUNCT__ 1395e27a552bSJed Brown #define __FUNCT__ "TSSetFromOptions_RosW" 1396e27a552bSJed Brown static PetscErrorCode TSSetFromOptions_RosW(TS ts) 1397e27a552bSJed Brown { 139861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1399e27a552bSJed Brown PetscErrorCode ierr; 140061692a83SJed Brown char rostype[256]; 1401e27a552bSJed Brown 1402e27a552bSJed Brown PetscFunctionBegin; 1403e27a552bSJed Brown ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 1404e27a552bSJed Brown { 140561692a83SJed Brown RosWTableauLink link; 1406e27a552bSJed Brown PetscInt count,choice; 1407e27a552bSJed Brown PetscBool flg; 1408e27a552bSJed Brown const char **namelist; 140961692a83SJed Brown SNES snes; 141061692a83SJed Brown 14118caf3d72SBarry Smith ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof(rostype));CHKERRQ(ierr); 141261692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1413e27a552bSJed Brown ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 141461692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 141561692a83SJed Brown ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 141661692a83SJed Brown ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 1417e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 141861692a83SJed Brown 14190298fd71SBarry Smith ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,NULL);CHKERRQ(ierr); 142061692a83SJed Brown 142161692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 142261692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 142361692a83SJed Brown if (!((PetscObject)snes)->type_name) { 142461692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 142561692a83SJed Brown } 142661692a83SJed Brown ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 1427e27a552bSJed Brown } 1428e27a552bSJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 1429e27a552bSJed Brown PetscFunctionReturn(0); 1430e27a552bSJed Brown } 1431e27a552bSJed Brown 1432e27a552bSJed Brown #undef __FUNCT__ 1433e27a552bSJed Brown #define __FUNCT__ "PetscFormatRealArray" 1434e27a552bSJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1435e27a552bSJed Brown { 1436e27a552bSJed Brown PetscErrorCode ierr; 1437e408995aSJed Brown PetscInt i; 1438e408995aSJed Brown size_t left,count; 1439e27a552bSJed Brown char *p; 1440e27a552bSJed Brown 1441e27a552bSJed Brown PetscFunctionBegin; 1442e408995aSJed Brown for (i=0,p=buf,left=len; i<n; i++) { 1443e408995aSJed Brown ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 1444e27a552bSJed Brown if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1445e27a552bSJed Brown left -= count; 1446e27a552bSJed Brown p += count; 1447e27a552bSJed Brown *p++ = ' '; 1448e27a552bSJed Brown } 1449e27a552bSJed Brown p[i ? 0 : -1] = 0; 1450e27a552bSJed Brown PetscFunctionReturn(0); 1451e27a552bSJed Brown } 1452e27a552bSJed Brown 1453e27a552bSJed Brown #undef __FUNCT__ 1454e27a552bSJed Brown #define __FUNCT__ "TSView_RosW" 1455e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1456e27a552bSJed Brown { 145761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 145861692a83SJed Brown RosWTableau tab = ros->tableau; 1459e27a552bSJed Brown PetscBool iascii; 1460e27a552bSJed Brown PetscErrorCode ierr; 1461ef20d060SBarry Smith TSAdapt adapt; 1462e27a552bSJed Brown 1463e27a552bSJed Brown PetscFunctionBegin; 1464251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1465e27a552bSJed Brown if (iascii) { 146619fd82e9SBarry Smith TSRosWType rostype; 1467e408995aSJed Brown PetscInt i; 1468e408995aSJed Brown PetscReal abscissa[512]; 1469e27a552bSJed Brown char buf[512]; 147061692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 147161692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 14728caf3d72SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 147361692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1474e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 14758caf3d72SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1476e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1477e27a552bSJed Brown } 1478ad6bc421SBarry Smith ierr = TSGetTSAdapt(ts,&adapt);CHKERRQ(ierr); 1479ef20d060SBarry Smith ierr = TSAdaptView(adapt,viewer);CHKERRQ(ierr); 1480e27a552bSJed Brown ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 1481e27a552bSJed Brown PetscFunctionReturn(0); 1482e27a552bSJed Brown } 1483e27a552bSJed Brown 1484e27a552bSJed Brown #undef __FUNCT__ 1485e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType" 1486e27a552bSJed Brown /*@C 148761692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1488e27a552bSJed Brown 1489e27a552bSJed Brown Logically collective 1490e27a552bSJed Brown 1491e27a552bSJed Brown Input Parameter: 1492e27a552bSJed Brown + ts - timestepping context 149361692a83SJed Brown - rostype - type of Rosenbrock-W scheme 1494e27a552bSJed Brown 1495020d8f30SJed Brown Level: beginner 1496e27a552bSJed Brown 1497020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1498e27a552bSJed Brown @*/ 149919fd82e9SBarry Smith PetscErrorCode TSRosWSetType(TS ts,TSRosWType rostype) 1500e27a552bSJed Brown { 1501e27a552bSJed Brown PetscErrorCode ierr; 1502e27a552bSJed Brown 1503e27a552bSJed Brown PetscFunctionBegin; 1504e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 150519fd82e9SBarry Smith ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,rostype));CHKERRQ(ierr); 1506e27a552bSJed Brown PetscFunctionReturn(0); 1507e27a552bSJed Brown } 1508e27a552bSJed Brown 1509e27a552bSJed Brown #undef __FUNCT__ 1510e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType" 1511e27a552bSJed Brown /*@C 151261692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1513e27a552bSJed Brown 1514e27a552bSJed Brown Logically collective 1515e27a552bSJed Brown 1516e27a552bSJed Brown Input Parameter: 1517e27a552bSJed Brown . ts - timestepping context 1518e27a552bSJed Brown 1519e27a552bSJed Brown Output Parameter: 152061692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1521e27a552bSJed Brown 1522e27a552bSJed Brown Level: intermediate 1523e27a552bSJed Brown 1524e27a552bSJed Brown .seealso: TSRosWGetType() 1525e27a552bSJed Brown @*/ 152619fd82e9SBarry Smith PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1527e27a552bSJed Brown { 1528e27a552bSJed Brown PetscErrorCode ierr; 1529e27a552bSJed Brown 1530e27a552bSJed Brown PetscFunctionBegin; 1531e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 153219fd82e9SBarry Smith ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1533e27a552bSJed Brown PetscFunctionReturn(0); 1534e27a552bSJed Brown } 1535e27a552bSJed Brown 1536e27a552bSJed Brown #undef __FUNCT__ 153761692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian" 1538e27a552bSJed Brown /*@C 153961692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1540e27a552bSJed Brown 1541e27a552bSJed Brown Logically collective 1542e27a552bSJed Brown 1543e27a552bSJed Brown Input Parameter: 1544e27a552bSJed Brown + ts - timestepping context 154561692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1546e27a552bSJed Brown 1547e27a552bSJed Brown Level: intermediate 1548e27a552bSJed Brown 1549e27a552bSJed Brown .seealso: TSRosWGetType() 1550e27a552bSJed Brown @*/ 155161692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1552e27a552bSJed Brown { 1553e27a552bSJed Brown PetscErrorCode ierr; 1554e27a552bSJed Brown 1555e27a552bSJed Brown PetscFunctionBegin; 1556e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 155761692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1558e27a552bSJed Brown PetscFunctionReturn(0); 1559e27a552bSJed Brown } 1560e27a552bSJed Brown 1561e27a552bSJed Brown #undef __FUNCT__ 1562e27a552bSJed Brown #define __FUNCT__ "TSRosWGetType_RosW" 156319fd82e9SBarry Smith PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1564e27a552bSJed Brown { 156561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1566e27a552bSJed Brown PetscErrorCode ierr; 1567e27a552bSJed Brown 1568e27a552bSJed Brown PetscFunctionBegin; 156961692a83SJed Brown if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 157061692a83SJed Brown *rostype = ros->tableau->name; 1571e27a552bSJed Brown PetscFunctionReturn(0); 1572e27a552bSJed Brown } 1573ef20d060SBarry Smith 1574e27a552bSJed Brown #undef __FUNCT__ 1575e27a552bSJed Brown #define __FUNCT__ "TSRosWSetType_RosW" 157619fd82e9SBarry Smith PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1577e27a552bSJed Brown { 157861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1579e27a552bSJed Brown PetscErrorCode ierr; 1580e27a552bSJed Brown PetscBool match; 158161692a83SJed Brown RosWTableauLink link; 1582e27a552bSJed Brown 1583e27a552bSJed Brown PetscFunctionBegin; 158461692a83SJed Brown if (ros->tableau) { 158561692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1586e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1587e27a552bSJed Brown } 158861692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 158961692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1590e27a552bSJed Brown if (match) { 1591e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 159261692a83SJed Brown ros->tableau = &link->tab; 1593e27a552bSJed Brown PetscFunctionReturn(0); 1594e27a552bSJed Brown } 1595e27a552bSJed Brown } 1596ce94432eSBarry Smith SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1597e27a552bSJed Brown PetscFunctionReturn(0); 1598e27a552bSJed Brown } 159961692a83SJed Brown 1600e27a552bSJed Brown #undef __FUNCT__ 160161692a83SJed Brown #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 160261692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1603e27a552bSJed Brown { 160461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1605e27a552bSJed Brown 1606e27a552bSJed Brown PetscFunctionBegin; 160761692a83SJed Brown ros->recompute_jacobian = flg; 1608e27a552bSJed Brown PetscFunctionReturn(0); 1609e27a552bSJed Brown } 1610e27a552bSJed Brown 1611d5e6173cSPeter Brune 1612e27a552bSJed Brown /* ------------------------------------------------------------ */ 1613e27a552bSJed Brown /*MC 1614020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1615e27a552bSJed Brown 1616e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1617e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1618e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1619e27a552bSJed Brown 1620e27a552bSJed Brown Notes: 162161692a83SJed Brown This method currently only works with autonomous ODE and DAE. 162261692a83SJed Brown 162361692a83SJed Brown Developer notes: 162461692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 162561692a83SJed Brown 1626f9c1d6abSBarry Smith $ udot = f(u) 162761692a83SJed Brown 162861692a83SJed Brown by the stage equations 162961692a83SJed Brown 1630f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 163161692a83SJed Brown 163261692a83SJed Brown and step completion formula 163361692a83SJed Brown 1634f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 163561692a83SJed Brown 1636f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 163761692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 163861692a83SJed Brown we define new variables for the stage equations 163961692a83SJed Brown 164061692a83SJed Brown $ y_i = gamma_ij k_j 164161692a83SJed Brown 164261692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 164361692a83SJed Brown 164461692a83SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 164561692a83SJed Brown 164661692a83SJed Brown to rewrite the method as 164761692a83SJed Brown 1648f9c1d6abSBarry 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 1649f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j bt_j y_j 165061692a83SJed Brown 165161692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 165261692a83SJed Brown 165361692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 165461692a83SJed Brown 165561692a83SJed Brown or, more compactly in tensor notation 165661692a83SJed Brown 165761692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 165861692a83SJed Brown 165961692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 166061692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 166161692a83SJed Brown equation 166261692a83SJed Brown 1663f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 166461692a83SJed Brown 166561692a83SJed Brown with initial guess y_i = 0. 1666e27a552bSJed Brown 1667e27a552bSJed Brown Level: beginner 1668e27a552bSJed Brown 1669a4386c9eSJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1670a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1671e27a552bSJed Brown M*/ 1672e27a552bSJed Brown #undef __FUNCT__ 1673e27a552bSJed Brown #define __FUNCT__ "TSCreate_RosW" 1674*8cc058d9SJed Brown PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1675e27a552bSJed Brown { 167661692a83SJed Brown TS_RosW *ros; 1677e27a552bSJed Brown PetscErrorCode ierr; 1678e27a552bSJed Brown 1679e27a552bSJed Brown PetscFunctionBegin; 1680e27a552bSJed Brown #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 16810298fd71SBarry Smith ierr = TSRosWInitializePackage(NULL);CHKERRQ(ierr); 1682e27a552bSJed Brown #endif 1683e27a552bSJed Brown 1684e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1685e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1686e27a552bSJed Brown ts->ops->view = TSView_RosW; 1687e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1688e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1689e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16901c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 1691e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1692e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1693e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1694e27a552bSJed Brown 169561692a83SJed Brown ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 169661692a83SJed Brown ts->data = (void*)ros; 1697e27a552bSJed Brown 169800de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 169900de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 170000de8ff0SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1701e27a552bSJed Brown PetscFunctionReturn(0); 1702e27a552bSJed Brown } 1703