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 */ 13af0996ceSBarry Smith #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 141e25c274SJed Brown #include <petscdm.h> 15e27a552bSJed Brown 16af0996ceSBarry Smith #include <petsc/private/kernels/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 */ 57be5899b3SLisandro Dalcin Vec vec_sol_prev; /* Solution from the previous step (used for interpolation and rollback)*/ 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: 11496a0c994SBarry Smith . 1. - 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: 12996a0c994SBarry Smith . 1. - 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: 14496a0c994SBarry Smith . 1. - 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: 16296a0c994SBarry Smith . 1. - 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: 17796a0c994SBarry Smith . 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: 19296a0c994SBarry Smith . 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: 20796a0c994SBarry Smith . 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: 22496a0c994SBarry Smith + 1. - Kaps and Rentrop, Generalized Runge Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979. 22596a0c994SBarry Smith - 2. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 22642faf41dSJed Brown 22742faf41dSJed Brown Hairer's code ros4.f 22842faf41dSJed Brown 22942faf41dSJed Brown Level: intermediate 23042faf41dSJed Brown 23142faf41dSJed Brown .seealso: TSROSW, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 23242faf41dSJed Brown M*/ 23342faf41dSJed Brown 23442faf41dSJed Brown /*MC 23542faf41dSJed Brown TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine 23642faf41dSJed Brown 23742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 23842faf41dSJed Brown 23942faf41dSJed Brown A-stable, |R(infty)| = 1/3. 24042faf41dSJed Brown 24142faf41dSJed Brown This method does not provide a dense output formula. 24242faf41dSJed Brown 24342faf41dSJed Brown References: 24496a0c994SBarry Smith + 1. - Shampine, Implementation of Rosenbrock methods, 1982. 24596a0c994SBarry Smith - 2. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 24642faf41dSJed Brown 24742faf41dSJed Brown Hairer's code ros4.f 24842faf41dSJed Brown 24942faf41dSJed Brown Level: intermediate 25042faf41dSJed Brown 25142faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWVELDD4, TSROSW4L 25242faf41dSJed Brown M*/ 25342faf41dSJed Brown 25442faf41dSJed Brown /*MC 25542faf41dSJed Brown TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen 25642faf41dSJed Brown 25742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 25842faf41dSJed Brown 25942faf41dSJed Brown A(89.5 degrees)-stable, |R(infty)| = 0.24. 26042faf41dSJed Brown 26142faf41dSJed Brown This method does not provide a dense output formula. 26242faf41dSJed Brown 26342faf41dSJed Brown References: 26496a0c994SBarry Smith + 1. - van Veldhuizen, D stability and Kaps Rentrop methods, 1984. 26596a0c994SBarry Smith - 2. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 26642faf41dSJed Brown 26742faf41dSJed Brown Hairer's code ros4.f 26842faf41dSJed Brown 26942faf41dSJed Brown Level: intermediate 27042faf41dSJed Brown 27142faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 27242faf41dSJed Brown M*/ 27342faf41dSJed Brown 27442faf41dSJed Brown /*MC 27542faf41dSJed Brown TSROSW4L - four stage, fourth order Rosenbrock (not W) method 27642faf41dSJed Brown 27742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 27842faf41dSJed Brown 27942faf41dSJed Brown A-stable and L-stable 28042faf41dSJed Brown 28142faf41dSJed Brown This method does not provide a dense output formula. 28242faf41dSJed Brown 28342faf41dSJed Brown References: 28496a0c994SBarry Smith . 1. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 28542faf41dSJed Brown 28642faf41dSJed Brown Hairer's code ros4.f 28742faf41dSJed Brown 28842faf41dSJed Brown Level: intermediate 28942faf41dSJed Brown 29042faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 29142faf41dSJed Brown M*/ 29242faf41dSJed Brown 293e27a552bSJed Brown /*@C 294be5899b3SLisandro Dalcin TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in TSRosW 295e27a552bSJed Brown 296e27a552bSJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 297e27a552bSJed Brown 298e27a552bSJed Brown Level: advanced 299e27a552bSJed Brown 300e27a552bSJed Brown .seealso: TSRosWRegisterDestroy() 301e27a552bSJed Brown @*/ 302e27a552bSJed Brown PetscErrorCode TSRosWRegisterAll(void) 303e27a552bSJed Brown { 304e27a552bSJed Brown PetscErrorCode ierr; 305e27a552bSJed Brown 306e27a552bSJed Brown PetscFunctionBegin; 307e27a552bSJed Brown if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 308e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_TRUE; 309e27a552bSJed Brown 310e27a552bSJed Brown { 311bbd56ea5SKarl Rupp const PetscReal A = 0; 312bbd56ea5SKarl Rupp const PetscReal Gamma = 1; 313bbd56ea5SKarl Rupp const PetscReal b = 1; 314bbd56ea5SKarl Rupp const PetscReal binterpt=1; 3151f80e275SEmil Constantinescu 3160298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3173606a31eSEmil Constantinescu } 3183606a31eSEmil Constantinescu 3193606a31eSEmil Constantinescu { 320bbd56ea5SKarl Rupp const PetscReal A = 0; 321bbd56ea5SKarl Rupp const PetscReal Gamma = 0.5; 322bbd56ea5SKarl Rupp const PetscReal b = 1; 323bbd56ea5SKarl Rupp const PetscReal binterpt=1; 324bbd56ea5SKarl Rupp 3250298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3263606a31eSEmil Constantinescu } 3273606a31eSEmil Constantinescu 3283606a31eSEmil Constantinescu { 329da80777bSKarl 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. */ 330e27a552bSJed Brown const PetscReal 33161692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 332da80777bSKarl Rupp Gamma[2][2] = {{1.707106781186547524401,0}, {-2.*1.707106781186547524401,1.707106781186547524401}}, 3331c3436cfSJed Brown b[2] = {0.5,0.5}, 3341c3436cfSJed Brown b1[2] = {1.0,0.0}; 3351f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 336da80777bSKarl Rupp binterpt[0][0] = 1.707106781186547524401 - 1.0; 337da80777bSKarl Rupp binterpt[1][0] = 2.0 - 1.707106781186547524401; 338da80777bSKarl Rupp binterpt[0][1] = 1.707106781186547524401 - 1.5; 339da80777bSKarl Rupp binterpt[1][1] = 1.5 - 1.707106781186547524401; 340bbd56ea5SKarl Rupp 3411f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 342e27a552bSJed Brown } 343e27a552bSJed Brown { 344da80777bSKarl 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. */ 345e27a552bSJed Brown const PetscReal 34661692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 347da80777bSKarl Rupp Gamma[2][2] = {{0.2928932188134524755992,0}, {-2.*0.2928932188134524755992,0.2928932188134524755992}}, 3481c3436cfSJed Brown b[2] = {0.5,0.5}, 3491c3436cfSJed Brown b1[2] = {1.0,0.0}; 3501f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 351da80777bSKarl Rupp binterpt[0][0] = 0.2928932188134524755992 - 1.0; 352da80777bSKarl Rupp binterpt[1][0] = 2.0 - 0.2928932188134524755992; 353da80777bSKarl Rupp binterpt[0][1] = 0.2928932188134524755992 - 1.5; 354da80777bSKarl Rupp binterpt[1][1] = 1.5 - 0.2928932188134524755992; 355bbd56ea5SKarl Rupp 3561f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 357fe7e6d57SJed Brown } 358fe7e6d57SJed Brown { 359da80777bSKarl Rupp /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */ 3601f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 361fe7e6d57SJed Brown const PetscReal 362fe7e6d57SJed Brown A[3][3] = {{0,0,0}, 363fe7e6d57SJed Brown {1.5773502691896257e+00,0,0}, 364fe7e6d57SJed Brown {0.5,0,0}}, 365da80777bSKarl Rupp Gamma[3][3] = {{7.8867513459481287e-01,0,0}, 366da80777bSKarl Rupp {-1.5773502691896257e+00,7.8867513459481287e-01,0}, 367da80777bSKarl Rupp {-6.7075317547305480e-01,-1.7075317547305482e-01,7.8867513459481287e-01}}, 368fe7e6d57SJed Brown b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01}, 369fe7e6d57SJed Brown b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01}; 3701f80e275SEmil Constantinescu 3711f80e275SEmil Constantinescu binterpt[0][0] = -0.8094010767585034; 3721f80e275SEmil Constantinescu binterpt[1][0] = -0.5; 3731f80e275SEmil Constantinescu binterpt[2][0] = 2.3094010767585034; 3741f80e275SEmil Constantinescu binterpt[0][1] = 0.9641016151377548; 3751f80e275SEmil Constantinescu binterpt[1][1] = 0.5; 3761f80e275SEmil Constantinescu binterpt[2][1] = -1.4641016151377548; 377bbd56ea5SKarl Rupp 3781f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 379fe7e6d57SJed Brown } 380fe7e6d57SJed Brown { 3813ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 382da80777bSKarl Rupp /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 383fe7e6d57SJed Brown const PetscReal 384fe7e6d57SJed Brown A[4][4] = {{0,0,0,0}, 385fe7e6d57SJed Brown {8.7173304301691801e-01,0,0,0}, 386fe7e6d57SJed Brown {8.4457060015369423e-01,-1.1299064236484185e-01,0,0}, 387fe7e6d57SJed Brown {0,0,1.,0}}, 388da80777bSKarl Rupp Gamma[4][4] = {{4.3586652150845900e-01,0,0,0}, 389da80777bSKarl Rupp {-8.7173304301691801e-01,4.3586652150845900e-01,0,0}, 390da80777bSKarl Rupp {-9.0338057013044082e-01,5.4180672388095326e-02,4.3586652150845900e-01,0}, 391da80777bSKarl Rupp {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,4.3586652150845900e-01}}, 392fe7e6d57SJed Brown b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01}, 3933ca35412SEmil Constantinescu b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01}; 3943ca35412SEmil Constantinescu 3953ca35412SEmil Constantinescu binterpt[0][0]=1.0564298455794094; 3963ca35412SEmil Constantinescu binterpt[1][0]=2.296429974281067; 3973ca35412SEmil Constantinescu binterpt[2][0]=-1.307599564525376; 3983ca35412SEmil Constantinescu binterpt[3][0]=-1.045260255335102; 3993ca35412SEmil Constantinescu binterpt[0][1]=-1.3864882699759573; 4003ca35412SEmil Constantinescu binterpt[1][1]=-8.262611700275677; 4013ca35412SEmil Constantinescu binterpt[2][1]=7.250979895056055; 4023ca35412SEmil Constantinescu binterpt[3][1]=2.398120075195581; 4033ca35412SEmil Constantinescu binterpt[0][2]=0.5721822314575016; 4043ca35412SEmil Constantinescu binterpt[1][2]=4.742931142090097; 4053ca35412SEmil Constantinescu binterpt[2][2]=-4.398120075195578; 4063ca35412SEmil Constantinescu binterpt[3][2]=-0.9169932983520199; 4073ca35412SEmil Constantinescu 4083ca35412SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 409e27a552bSJed Brown } 410ef3c5b88SJed Brown { 411da80777bSKarl Rupp /* const PetscReal g = 0.5; Directly written in-place below */ 412ef3c5b88SJed Brown const PetscReal 413ef3c5b88SJed Brown A[4][4] = {{0,0,0,0}, 414ef3c5b88SJed Brown {0,0,0,0}, 415ef3c5b88SJed Brown {1.,0,0,0}, 416ef3c5b88SJed Brown {0.75,-0.25,0.5,0}}, 417da80777bSKarl Rupp Gamma[4][4] = {{0.5,0,0,0}, 418da80777bSKarl Rupp {1.,0.5,0,0}, 419da80777bSKarl Rupp {-0.25,-0.25,0.5,0}, 420da80777bSKarl Rupp {1./12,1./12,-2./3,0.5}}, 421ef3c5b88SJed Brown b[4] = {5./6,-1./6,-1./6,0.5}, 422ef3c5b88SJed Brown b2[4] = {0.75,-0.25,0.5,0}; 423bbd56ea5SKarl Rupp 4240298fd71SBarry Smith ierr = TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,NULL);CHKERRQ(ierr); 425ef3c5b88SJed Brown } 426ef3c5b88SJed Brown { 427da80777bSKarl Rupp /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */ 428ef3c5b88SJed Brown const PetscReal 429ef3c5b88SJed Brown A[3][3] = {{0,0,0}, 430da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}, 431da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}}, 432da80777bSKarl Rupp Gamma[3][3] = {{0.43586652150845899941601945119356,0,0}, 433da80777bSKarl Rupp {-0.19294655696029095575009695436041,0.43586652150845899941601945119356,0}, 434da80777bSKarl Rupp {0,1.74927148125794685173529749738960,0.43586652150845899941601945119356}}, 435ef3c5b88SJed Brown b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829}, 436ef3c5b88SJed Brown b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619}; 4371f80e275SEmil Constantinescu 4381f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 4391f80e275SEmil Constantinescu binterpt[0][0] = 3.793692883777660870425141387941; 4401f80e275SEmil Constantinescu binterpt[1][0] = -2.918692883777660870425141387941; 4411f80e275SEmil Constantinescu binterpt[2][0] = 0.125; 4421f80e275SEmil Constantinescu binterpt[0][1] = -0.725741064379812106687651020584; 4431f80e275SEmil Constantinescu binterpt[1][1] = 0.559074397713145440020984353917; 4441f80e275SEmil Constantinescu binterpt[2][1] = 0.16666666666666666666666666666667; 4451f80e275SEmil Constantinescu 4461f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 447ef3c5b88SJed Brown } 448b1c69cc3SEmil Constantinescu { 449da80777bSKarl Rupp /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 450da80777bSKarl Rupp * Direct evaluation: s3 = 1.732050807568877293527; 451da80777bSKarl Rupp * g = 0.7886751345948128822546; 452da80777bSKarl Rupp * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */ 453b1c69cc3SEmil Constantinescu const PetscReal 454b1c69cc3SEmil Constantinescu A[3][3] = {{0,0,0}, 455b1c69cc3SEmil Constantinescu {1,0,0}, 456b1c69cc3SEmil Constantinescu {0.25,0.25,0}}, 457b1c69cc3SEmil Constantinescu Gamma[3][3] = {{0,0,0}, 458da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/6.0,0.7886751345948128822546,0}, 459da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/24.0,(-3.0-1.732050807568877293527)/8.0,0.7886751345948128822546}}, 460b1c69cc3SEmil Constantinescu b[3] = {1./6.,1./6.,2./3.}, 461b1c69cc3SEmil Constantinescu b2[3] = {1./4.,1./4.,1./2.}; 462c0cb691aSEmil Constantinescu PetscReal binterpt[3][2]; 463da80777bSKarl Rupp 464c0cb691aSEmil Constantinescu binterpt[0][0]=0.089316397477040902157517886164709; 465c0cb691aSEmil Constantinescu binterpt[1][0]=-0.91068360252295909784248211383529; 466c0cb691aSEmil Constantinescu binterpt[2][0]=1.8213672050459181956849642276706; 467c0cb691aSEmil Constantinescu binterpt[0][1]=0.077350269189625764509148780501957; 468c0cb691aSEmil Constantinescu binterpt[1][1]=1.077350269189625764509148780502; 469c0cb691aSEmil Constantinescu binterpt[2][1]=-1.1547005383792515290182975610039; 470bbd56ea5SKarl Rupp 471c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 472b1c69cc3SEmil Constantinescu } 473b1c69cc3SEmil Constantinescu 474b1c69cc3SEmil Constantinescu { 475b1c69cc3SEmil Constantinescu const PetscReal 476b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 477b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 478b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 479b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 480b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 481b1c69cc3SEmil Constantinescu {0.0,1./4.,0,0}, 482b1c69cc3SEmil Constantinescu {-2.,-2./3.,2./3.,0}, 483b1c69cc3SEmil Constantinescu {1./2.,5./36.,-2./9,0}}, 484b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 485b1c69cc3SEmil Constantinescu b2[4] = {1./8.,3./4.,1./8.,0}; 486c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 487da80777bSKarl Rupp 488c0cb691aSEmil Constantinescu binterpt[0][0]=6.25; 489c0cb691aSEmil Constantinescu binterpt[1][0]=-30.25; 490c0cb691aSEmil Constantinescu binterpt[2][0]=1.75; 491c0cb691aSEmil Constantinescu binterpt[3][0]=23.25; 492c0cb691aSEmil Constantinescu binterpt[0][1]=-9.75; 493c0cb691aSEmil Constantinescu binterpt[1][1]=58.75; 494c0cb691aSEmil Constantinescu binterpt[2][1]=-3.25; 495c0cb691aSEmil Constantinescu binterpt[3][1]=-45.75; 496c0cb691aSEmil Constantinescu binterpt[0][2]=3.6666666666666666666666666666667; 497c0cb691aSEmil Constantinescu binterpt[1][2]=-28.333333333333333333333333333333; 498c0cb691aSEmil Constantinescu binterpt[2][2]=1.6666666666666666666666666666667; 499c0cb691aSEmil Constantinescu binterpt[3][2]=23.; 500bbd56ea5SKarl Rupp 501c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 502b1c69cc3SEmil Constantinescu } 503b1c69cc3SEmil Constantinescu 504b1c69cc3SEmil Constantinescu { 505b1c69cc3SEmil Constantinescu const PetscReal 506b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 507b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 508b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 509b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 510b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 511b1c69cc3SEmil Constantinescu {0.0,3./4.,0,0}, 512b1c69cc3SEmil Constantinescu {-2./3.,-23./9.,2./9.,0}, 513b1c69cc3SEmil Constantinescu {1./18.,65./108.,-2./27,0}}, 514b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 515b1c69cc3SEmil Constantinescu b2[4] = {3./16.,10./16.,3./16.,0}; 516c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 517da80777bSKarl Rupp 518c0cb691aSEmil Constantinescu binterpt[0][0]=1.6911764705882352941176470588235; 519c0cb691aSEmil Constantinescu binterpt[1][0]=3.6813725490196078431372549019608; 520c0cb691aSEmil Constantinescu binterpt[2][0]=0.23039215686274509803921568627451; 521c0cb691aSEmil Constantinescu binterpt[3][0]=-4.6029411764705882352941176470588; 522c0cb691aSEmil Constantinescu binterpt[0][1]=-0.95588235294117647058823529411765; 523c0cb691aSEmil Constantinescu binterpt[1][1]=-6.2401960784313725490196078431373; 524c0cb691aSEmil Constantinescu binterpt[2][1]=-0.31862745098039215686274509803922; 525c0cb691aSEmil Constantinescu binterpt[3][1]=7.5147058823529411764705882352941; 526c0cb691aSEmil Constantinescu binterpt[0][2]=-0.56862745098039215686274509803922; 527c0cb691aSEmil Constantinescu binterpt[1][2]=2.7254901960784313725490196078431; 528c0cb691aSEmil Constantinescu binterpt[2][2]=0.25490196078431372549019607843137; 529c0cb691aSEmil Constantinescu binterpt[3][2]=-2.4117647058823529411764705882353; 530bbd56ea5SKarl Rupp 531c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 532b1c69cc3SEmil Constantinescu } 533753f8adbSEmil Constantinescu 534753f8adbSEmil Constantinescu { 535753f8adbSEmil Constantinescu PetscReal A[4][4],Gamma[4][4],b[4],b2[4]; 5363ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 537753f8adbSEmil Constantinescu 538753f8adbSEmil Constantinescu Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816; 53905e8e825SJed Brown Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0; 540753f8adbSEmil Constantinescu Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476; 541753f8adbSEmil Constantinescu Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816; 54205e8e825SJed Brown Gamma[1][2]=0; Gamma[1][3]=0; 543753f8adbSEmil Constantinescu Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903; 544753f8adbSEmil Constantinescu Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131; 545753f8adbSEmil Constantinescu Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816; 54605e8e825SJed Brown Gamma[2][3]=0; 547753f8adbSEmil Constantinescu Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783; 548753f8adbSEmil Constantinescu Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984; 549753f8adbSEmil Constantinescu Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198; 550753f8adbSEmil Constantinescu Gamma[3][3]=0; 551753f8adbSEmil Constantinescu 55205e8e825SJed Brown A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0; 553753f8adbSEmil Constantinescu A[1][0]=0.8717330430169179988320388950590125027645343373957631; 55405e8e825SJed Brown A[1][1]=0; A[1][2]=0; A[1][3]=0; 555753f8adbSEmil Constantinescu A[2][0]=0.5275890119763004115618079766722914408876108660811028; 556753f8adbSEmil Constantinescu A[2][1]=0.07241098802369958843819203208518599088698057726988732; 55705e8e825SJed Brown A[2][2]=0; A[2][3]=0; 558753f8adbSEmil Constantinescu A[3][0]=0.3990960076760701320627260685975778145384666450351314; 559753f8adbSEmil Constantinescu A[3][1]=-0.4375576546135194437228463747348862825846903771419953; 560753f8adbSEmil Constantinescu A[3][2]=1.038461646937449311660120300601880176655352737312713; 56105e8e825SJed Brown A[3][3]=0; 562753f8adbSEmil Constantinescu 563753f8adbSEmil Constantinescu b[0]=0.1876410243467238251612921333138006734899663569186926; 564753f8adbSEmil Constantinescu b[1]=-0.5952974735769549480478230473706443582188442040780541; 565753f8adbSEmil Constantinescu b[2]=0.9717899277217721234705114616271378792182450260943198; 566753f8adbSEmil Constantinescu b[3]=0.4358665215084589994160194475295062513822671686978816; 567753f8adbSEmil Constantinescu 568753f8adbSEmil Constantinescu b2[0]=0.2147402862233891404862383521089097657790734483804460; 569753f8adbSEmil Constantinescu b2[1]=-0.4851622638849390928209050538171743017757490232519684; 570753f8adbSEmil Constantinescu b2[2]=0.8687250025203875511662123688667549217531982787600080; 571753f8adbSEmil Constantinescu b2[3]=0.4016969751411624011684543450940068201770721128357014; 572753f8adbSEmil Constantinescu 5733ca35412SEmil Constantinescu binterpt[0][0]=2.2565812720167954547104627844105; 5743ca35412SEmil Constantinescu binterpt[1][0]=1.349166413351089573796243820819; 5753ca35412SEmil Constantinescu binterpt[2][0]=-2.4695174540533503758652847586647; 5763ca35412SEmil Constantinescu binterpt[3][0]=-0.13623023131453465264142184656474; 5773ca35412SEmil Constantinescu binterpt[0][1]=-3.0826699111559187902922463354557; 5783ca35412SEmil Constantinescu binterpt[1][1]=-2.4689115685996042534544925650515; 5793ca35412SEmil Constantinescu binterpt[2][1]=5.7428279814696677152129332773553; 5803ca35412SEmil Constantinescu binterpt[3][1]=-0.19124650171414467146619437684812; 5813ca35412SEmil Constantinescu binterpt[0][2]=1.0137296634858471607430756831148; 5823ca35412SEmil Constantinescu binterpt[1][2]=0.52444768167155973161042570784064; 5833ca35412SEmil Constantinescu binterpt[2][2]=-2.3015205996945452158771370439586; 5843ca35412SEmil Constantinescu binterpt[3][2]=0.76334325453713832352363565300308; 585f4aed992SEmil Constantinescu 586f73f8d2cSSatish Balay ierr = TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 587753f8adbSEmil Constantinescu } 58842faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWGRK4T,0.231,PETSC_DEFAULT,PETSC_DEFAULT,0,-0.1282612945269037e+01);CHKERRQ(ierr); 58942faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWSHAMP4,0.5,PETSC_DEFAULT,PETSC_DEFAULT,0,125./108.);CHKERRQ(ierr); 59042faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWVELDD4,0.22570811482256823492,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.355958941201148);CHKERRQ(ierr); 59142faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSW4L,0.57282,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.093502252409163);CHKERRQ(ierr); 592e27a552bSJed Brown PetscFunctionReturn(0); 593e27a552bSJed Brown } 594e27a552bSJed Brown 595d5e6173cSPeter Brune 596d5e6173cSPeter Brune 597e27a552bSJed Brown /*@C 598e27a552bSJed Brown TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 599e27a552bSJed Brown 600e27a552bSJed Brown Not Collective 601e27a552bSJed Brown 602e27a552bSJed Brown Level: advanced 603e27a552bSJed Brown 604607a6623SBarry Smith .seealso: TSRosWRegister(), TSRosWRegisterAll() 605e27a552bSJed Brown @*/ 606e27a552bSJed Brown PetscErrorCode TSRosWRegisterDestroy(void) 607e27a552bSJed Brown { 608e27a552bSJed Brown PetscErrorCode ierr; 60961692a83SJed Brown RosWTableauLink link; 610e27a552bSJed Brown 611e27a552bSJed Brown PetscFunctionBegin; 61261692a83SJed Brown while ((link = RosWTableauList)) { 61361692a83SJed Brown RosWTableau t = &link->tab; 61461692a83SJed Brown RosWTableauList = link->next; 61561692a83SJed Brown ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 61643b21953SEmil Constantinescu ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 617fe7e6d57SJed Brown ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 618f4aed992SEmil Constantinescu ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 619e27a552bSJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 620e27a552bSJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 621e27a552bSJed Brown } 622e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 623e27a552bSJed Brown PetscFunctionReturn(0); 624e27a552bSJed Brown } 625e27a552bSJed Brown 626e27a552bSJed Brown /*@C 627e27a552bSJed Brown TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 6288a690491SBarry Smith from TSInitializePackage(). 629e27a552bSJed Brown 630e27a552bSJed Brown Level: developer 631e27a552bSJed Brown 632e27a552bSJed Brown .seealso: PetscInitialize() 633e27a552bSJed Brown @*/ 634607a6623SBarry Smith PetscErrorCode TSRosWInitializePackage(void) 635e27a552bSJed Brown { 636e27a552bSJed Brown PetscErrorCode ierr; 637e27a552bSJed Brown 638e27a552bSJed Brown PetscFunctionBegin; 639e27a552bSJed Brown if (TSRosWPackageInitialized) PetscFunctionReturn(0); 640e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 641e27a552bSJed Brown ierr = TSRosWRegisterAll();CHKERRQ(ierr); 642e27a552bSJed Brown ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 643e27a552bSJed Brown PetscFunctionReturn(0); 644e27a552bSJed Brown } 645e27a552bSJed Brown 646e27a552bSJed Brown /*@C 647e27a552bSJed Brown TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 648e27a552bSJed Brown called from PetscFinalize(). 649e27a552bSJed Brown 650e27a552bSJed Brown Level: developer 651e27a552bSJed Brown 652e27a552bSJed Brown .seealso: PetscFinalize() 653e27a552bSJed Brown @*/ 654e27a552bSJed Brown PetscErrorCode TSRosWFinalizePackage(void) 655e27a552bSJed Brown { 656e27a552bSJed Brown PetscErrorCode ierr; 657e27a552bSJed Brown 658e27a552bSJed Brown PetscFunctionBegin; 659e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 660e27a552bSJed Brown ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 661e27a552bSJed Brown PetscFunctionReturn(0); 662e27a552bSJed Brown } 663e27a552bSJed Brown 664e27a552bSJed Brown /*@C 66561692a83SJed Brown TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 666e27a552bSJed Brown 667e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 668e27a552bSJed Brown 669e27a552bSJed Brown Input Parameters: 670e27a552bSJed Brown + name - identifier for method 671e27a552bSJed Brown . order - approximation order of method 672e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 67361692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 67461692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 675fe7e6d57SJed Brown . b - Step completion table (dimension s) 6760298fd71SBarry Smith . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available) 677f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 67842faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 679e27a552bSJed Brown 680e27a552bSJed Brown Notes: 68161692a83SJed Brown Several Rosenbrock W methods are provided, this function is only needed to create new methods. 682e27a552bSJed Brown 683e27a552bSJed Brown Level: advanced 684e27a552bSJed Brown 685e27a552bSJed Brown .seealso: TSRosW 686e27a552bSJed Brown @*/ 687f9c1d6abSBarry Smith PetscErrorCode TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 688f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 689e27a552bSJed Brown { 690e27a552bSJed Brown PetscErrorCode ierr; 69161692a83SJed Brown RosWTableauLink link; 69261692a83SJed Brown RosWTableau t; 69361692a83SJed Brown PetscInt i,j,k; 69461692a83SJed Brown PetscScalar *GammaInv; 695e27a552bSJed Brown 696e27a552bSJed Brown PetscFunctionBegin; 697fe7e6d57SJed Brown PetscValidCharPointer(name,1); 698fe7e6d57SJed Brown PetscValidPointer(A,4); 699fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 700fe7e6d57SJed Brown PetscValidPointer(b,6); 701fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 702fe7e6d57SJed Brown 7031d36bdfdSBarry Smith ierr = TSRosWInitializePackage();CHKERRQ(ierr); 704*071fcb05SBarry Smith ierr = PetscNew(&link);CHKERRQ(ierr); 705e27a552bSJed Brown t = &link->tab; 706e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 707e27a552bSJed Brown t->order = order; 708e27a552bSJed Brown t->s = s; 709dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->A,s*s,&t->Gamma,s,&t->b,s,&t->ASum,s,&t->GammaSum);CHKERRQ(ierr); 710dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->At,s,&t->bt,s*s,&t->GammaInv,s,&t->GammaZeroDiag,s*s,&t->GammaExplicitCorr);CHKERRQ(ierr); 711580bdb30SBarry Smith ierr = PetscArraycpy(t->A,A,s*s);CHKERRQ(ierr); 712580bdb30SBarry Smith ierr = PetscArraycpy(t->Gamma,Gamma,s*s);CHKERRQ(ierr); 713580bdb30SBarry Smith ierr = PetscArraycpy(t->GammaExplicitCorr,Gamma,s*s);CHKERRQ(ierr); 714580bdb30SBarry Smith ierr = PetscArraycpy(t->b,b,s);CHKERRQ(ierr); 715fe7e6d57SJed Brown if (bembed) { 716dcca6d9dSJed Brown ierr = PetscMalloc2(s,&t->bembed,s,&t->bembedt);CHKERRQ(ierr); 717580bdb30SBarry Smith ierr = PetscArraycpy(t->bembed,bembed,s);CHKERRQ(ierr); 718fe7e6d57SJed Brown } 71961692a83SJed Brown for (i=0; i<s; i++) { 72061692a83SJed Brown t->ASum[i] = 0; 72161692a83SJed Brown t->GammaSum[i] = 0; 72261692a83SJed Brown for (j=0; j<s; j++) { 72361692a83SJed Brown t->ASum[i] += A[i*s+j]; 724fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 72561692a83SJed Brown } 72661692a83SJed Brown } 727785e854fSJed Brown ierr = PetscMalloc1(s*s,&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 72861692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 729fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 730fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 731fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 732c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 733fd96d5b0SEmil Constantinescu } else { 734c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 735fd96d5b0SEmil Constantinescu } 736fd96d5b0SEmil Constantinescu } 737fd96d5b0SEmil Constantinescu 73861692a83SJed Brown switch (s) { 73961692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 7402e92ee13SHong Zhang case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7416baedc03SHong Zhang case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7422e92ee13SHong Zhang case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 74361692a83SJed Brown case 5: { 74461692a83SJed Brown PetscInt ipvt5[5]; 74561692a83SJed Brown MatScalar work5[5*5]; 7462e92ee13SHong Zhang ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 74761692a83SJed Brown } 7482e92ee13SHong Zhang case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7492e92ee13SHong Zhang case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 75061692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 75161692a83SJed Brown } 75261692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 75361692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 75443b21953SEmil Constantinescu 75543b21953SEmil Constantinescu for (i=0; i<s; i++) { 75643b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 75743b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 75843b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 75943b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 76043b21953SEmil Constantinescu } 76143b21953SEmil Constantinescu } 76243b21953SEmil Constantinescu } 76343b21953SEmil Constantinescu 76461692a83SJed Brown for (i=0; i<s; i++) { 76561692a83SJed Brown for (j=0; j<s; j++) { 76661692a83SJed Brown t->At[i*s+j] = 0; 76761692a83SJed Brown for (k=0; k<s; k++) { 76861692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 76961692a83SJed Brown } 77061692a83SJed Brown } 77161692a83SJed Brown t->bt[i] = 0; 77261692a83SJed Brown for (j=0; j<s; j++) { 77361692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 77461692a83SJed Brown } 775fe7e6d57SJed Brown if (bembed) { 776fe7e6d57SJed Brown t->bembedt[i] = 0; 777fe7e6d57SJed Brown for (j=0; j<s; j++) { 778fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 779fe7e6d57SJed Brown } 780fe7e6d57SJed Brown } 78161692a83SJed Brown } 7828d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 7838d59e960SJed Brown 784f4aed992SEmil Constantinescu t->pinterp = pinterp; 785785e854fSJed Brown ierr = PetscMalloc1(s*pinterp,&t->binterpt);CHKERRQ(ierr); 786580bdb30SBarry Smith ierr = PetscArraycpy(t->binterpt,binterpt,s*pinterp);CHKERRQ(ierr); 78761692a83SJed Brown link->next = RosWTableauList; 78861692a83SJed Brown RosWTableauList = link; 789e27a552bSJed Brown PetscFunctionReturn(0); 790e27a552bSJed Brown } 791e27a552bSJed Brown 79242faf41dSJed Brown /*@C 79342faf41dSJed Brown TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing paramter choices 79442faf41dSJed Brown 79542faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 79642faf41dSJed Brown 79742faf41dSJed Brown Input Parameters: 79842faf41dSJed Brown + name - identifier for method 79942faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 80042faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 80142faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 80242faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 80342faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 80442faf41dSJed Brown . e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 80542faf41dSJed Brown 80642faf41dSJed Brown Notes: 80742faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 80842faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 80942faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 81042faf41dSJed Brown 81142faf41dSJed Brown Level: developer 81242faf41dSJed Brown 81342faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 81442faf41dSJed Brown @*/ 81519fd82e9SBarry Smith PetscErrorCode TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 81642faf41dSJed Brown { 81742faf41dSJed Brown PetscErrorCode ierr; 81842faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 81942faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 82042faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 82142faf41dSJed Brown p42 = one/eight - gamma/three, 82242faf41dSJed Brown p43 = one/twelve - gamma/three, 82342faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 82442faf41dSJed Brown p56 = one/twenty - gamma/four; 82542faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 82642faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 82742faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 82842faf41dSJed Brown 82942faf41dSJed Brown PetscFunctionBegin; 83042faf41dSJed Brown /* Step 1: choose Gamma (input) */ 83142faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 83242faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 83342faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 83442faf41dSJed Brown 83542faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 83642faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 83742faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 83842faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 83942faf41dSJed Brown rhs[0] = one - b3; 84042faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 84142faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 8426baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 84342faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 84442faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 84542faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 84642faf41dSJed Brown 84742faf41dSJed Brown /* Step 3 */ 84842faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 84942faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 85042faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 85142faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 85242faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 85342faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 85442faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 8556baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 85642faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 85742faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 85842faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 85942faf41dSJed Brown 86042faf41dSJed Brown /* Step 4: back-substitute */ 86142faf41dSJed Brown beta32 = beta32beta2p / beta2p; 86242faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 86342faf41dSJed Brown 86442faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 86542faf41dSJed Brown a43 = 0; 86642faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 86742faf41dSJed Brown a42 = a32; 86842faf41dSJed Brown 86942faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 87042faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 87142faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 87242faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 87342faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 87442faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 87542faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 87642faf41dSJed 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; 87742faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 87842faf41dSJed Brown 87942faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 88042faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 88142faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 88242faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 88342faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 88442faf41dSJed Brown 88542faf41dSJed Brown { 88642faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 88742faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 88842faf41dSJed Brown } 8890298fd71SBarry Smith ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,NULL);CHKERRQ(ierr); 89042faf41dSJed Brown PetscFunctionReturn(0); 89142faf41dSJed Brown } 89242faf41dSJed Brown 8931c3436cfSJed Brown /* 8941c3436cfSJed Brown The step completion formula is 8951c3436cfSJed Brown 8961c3436cfSJed Brown x1 = x0 + b^T Y 8971c3436cfSJed Brown 8981c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 8991c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 9001c3436cfSJed Brown 9011c3436cfSJed Brown x1e = x0 + be^T Y 9021c3436cfSJed Brown = x1 - b^T Y + be^T Y 9031c3436cfSJed Brown = x1 + (be - b)^T Y 9041c3436cfSJed Brown 9051c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9061c3436cfSJed Brown */ 907f9c1d6abSBarry Smith static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done) 9081c3436cfSJed Brown { 9091c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9101c3436cfSJed Brown RosWTableau tab = ros->tableau; 9111c3436cfSJed Brown PetscScalar *w = ros->work; 9121c3436cfSJed Brown PetscInt i; 9131c3436cfSJed Brown PetscErrorCode ierr; 9141c3436cfSJed Brown 9151c3436cfSJed Brown PetscFunctionBegin; 9161c3436cfSJed Brown if (order == tab->order) { 917108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 918f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 919de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 920f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 921f9c1d6abSBarry Smith } else {ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr);} 9221c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9231c3436cfSJed Brown PetscFunctionReturn(0); 9241c3436cfSJed Brown } else if (order == tab->order-1) { 9251c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 926108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 927f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 928de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 929f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 930108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 931108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 932f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 933f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 9341c3436cfSJed Brown } 9351c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9361c3436cfSJed Brown PetscFunctionReturn(0); 9371c3436cfSJed Brown } 9381c3436cfSJed Brown unavailable: 9391c3436cfSJed Brown if (done) *done = PETSC_FALSE; 940a7fac7c2SEmil Constantinescu else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order); 9411c3436cfSJed Brown PetscFunctionReturn(0); 9421c3436cfSJed Brown } 9431c3436cfSJed Brown 944560360afSLisandro Dalcin static PetscErrorCode TSRollBack_RosW(TS ts) 94524655328SShri { 94624655328SShri TS_RosW *ros = (TS_RosW*)ts->data; 94724655328SShri PetscErrorCode ierr; 94824655328SShri 94924655328SShri PetscFunctionBegin; 950be5899b3SLisandro Dalcin ierr = VecCopy(ros->vec_sol_prev,ts->vec_sol);CHKERRQ(ierr); 95124655328SShri PetscFunctionReturn(0); 95224655328SShri } 95324655328SShri 954e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 955e27a552bSJed Brown { 95661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 95761692a83SJed Brown RosWTableau tab = ros->tableau; 958e27a552bSJed Brown const PetscInt s = tab->s; 9591c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9600feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 961c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 96261692a83SJed Brown PetscScalar *w = ros->work; 9637d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 964e27a552bSJed Brown SNES snes; 9651c3436cfSJed Brown TSAdapt adapt; 966fecfb714SLisandro Dalcin PetscInt i,j,its,lits; 967be5899b3SLisandro Dalcin PetscInt rejections = 0; 968b39943a6SLisandro Dalcin PetscBool stageok,accept = PETSC_TRUE; 969b39943a6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 970e27a552bSJed Brown PetscErrorCode ierr; 971e27a552bSJed Brown 972e27a552bSJed Brown PetscFunctionBegin; 973b39943a6SLisandro Dalcin if (!ts->steprollback) { 974be5899b3SLisandro Dalcin ierr = VecCopy(ts->vec_sol,ros->vec_sol_prev);CHKERRQ(ierr); 975b39943a6SLisandro Dalcin } 976e27a552bSJed Brown 977b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 978b39943a6SLisandro Dalcin while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 9791c3436cfSJed Brown const PetscReal h = ts->time_step; 980e27a552bSJed Brown for (i=0; i<s; i++) { 9811c3436cfSJed Brown ros->stage_time = ts->ptime + h*ASum[i]; 982b8123daeSJed Brown ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 983c17803e7SJed Brown if (GammaZeroDiag[i]) { 984c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 985b296d7d5SJed Brown ros->scoeff = 1.; 986c17803e7SJed Brown } else { 987c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 988b296d7d5SJed Brown ros->scoeff = 1./Gamma[i*s+i]; 989fd96d5b0SEmil Constantinescu } 99061692a83SJed Brown 99161692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 992de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 993de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 99461692a83SJed Brown 99561692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 99661692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 99761692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 99861692a83SJed Brown 999e27a552bSJed Brown /* Initial guess taken from last stage */ 100061692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 100161692a83SJed Brown 10027d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 1003b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 100461692a83SJed Brown if (!ros->recompute_jacobian && !i) { 100561692a83SJed Brown ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 100661692a83SJed Brown } 10070298fd71SBarry Smith ierr = SNESSolve(snes,NULL,Y[i]);CHKERRQ(ierr); 1008e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1009e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10105ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 10117d4bf2deSEmil Constantinescu } else { 10121ce71dffSSatish Balay Mat J,Jp; 10130feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10140feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 101522d28d08SBarry Smith ierr = VecScale(Y[i],-1.0);CHKERRQ(ierr); 10160feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 10170feba352SEmil Constantinescu 10180feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10190feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10200feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 1021fecfb714SLisandro Dalcin 1022fecfb714SLisandro Dalcin /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10230298fd71SBarry Smith ierr = TSGetIJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 1024d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,J,Jp,PETSC_FALSE);CHKERRQ(ierr); 102522d28d08SBarry Smith ierr = MatMult(J,Zstage,Zdot);CHKERRQ(ierr); 10260feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10275ef26d82SJed Brown ts->ksp_its += 1; 1028fecfb714SLisandro Dalcin 1029fecfb714SLisandro Dalcin ierr = VecScale(Y[i],h);CHKERRQ(ierr); 10307d4bf2deSEmil Constantinescu } 10319be3e283SDebojyoti Ghosh ierr = TSPostStage(ts,ros->stage_time,i,Y);CHKERRQ(ierr); 1032fecfb714SLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1033fecfb714SLisandro Dalcin ierr = TSAdaptCheckStage(adapt,ts,ros->stage_time,Y[i],&stageok);CHKERRQ(ierr); 1034fecfb714SLisandro Dalcin if (!stageok) goto reject_step; 1035e27a552bSJed Brown } 1036e27a552bSJed Brown 1037b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 1038b39943a6SLisandro Dalcin ierr = TSEvaluateStep_RosW(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 1039b39943a6SLisandro Dalcin ros->status = TS_STEP_PENDING; 1040552698daSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 10411c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10421917a363SLisandro Dalcin ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 1043fecfb714SLisandro Dalcin ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 1044b39943a6SLisandro Dalcin ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1045b39943a6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 1046b39943a6SLisandro Dalcin ierr = TSRollBack_RosW(ts);CHKERRQ(ierr); 1047be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1048be5899b3SLisandro Dalcin goto reject_step; 1049b39943a6SLisandro Dalcin } 1050b39943a6SLisandro Dalcin 1051e27a552bSJed Brown ts->ptime += ts->time_step; 1052cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 10531c3436cfSJed Brown break; 1054b39943a6SLisandro Dalcin 1055b39943a6SLisandro Dalcin reject_step: 1056fecfb714SLisandro Dalcin ts->reject++; accept = PETSC_FALSE; 1057be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1058b39943a6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 1059be5899b3SLisandro Dalcin ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 10601c3436cfSJed Brown } 10611c3436cfSJed Brown } 1062e27a552bSJed Brown PetscFunctionReturn(0); 1063e27a552bSJed Brown } 1064e27a552bSJed Brown 1065f9c1d6abSBarry Smith static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1066e27a552bSJed Brown { 106761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1068f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1069f4aed992SEmil Constantinescu PetscReal h; 1070f4aed992SEmil Constantinescu PetscReal tt,t; 1071f4aed992SEmil Constantinescu PetscScalar *bt; 1072f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1073f4aed992SEmil Constantinescu PetscErrorCode ierr; 1074f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1075f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1076f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1077e27a552bSJed Brown 1078e27a552bSJed Brown PetscFunctionBegin; 1079ce94432eSBarry Smith if (!Bt) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1080f4aed992SEmil Constantinescu 1081f4aed992SEmil Constantinescu switch (ros->status) { 1082f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1083f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1084f4aed992SEmil Constantinescu h = ts->time_step; 1085f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1086f4aed992SEmil Constantinescu break; 1087f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1088be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1089f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1090f4aed992SEmil Constantinescu break; 1091ce94432eSBarry Smith default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1092f4aed992SEmil Constantinescu } 1093785e854fSJed Brown ierr = PetscMalloc1(s,&bt);CHKERRQ(ierr); 1094f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1095f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1096f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 10973ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1098f4aed992SEmil Constantinescu } 1099f4aed992SEmil Constantinescu } 1100f4aed992SEmil Constantinescu 1101f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1102f9c1d6abSBarry Smith /* U <- 0*/ 1103f9c1d6abSBarry Smith ierr = VecZeroEntries(U);CHKERRQ(ierr); 1104f9c1d6abSBarry Smith /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11053ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j] = 0; 11063ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11073ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11083ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1109f4aed992SEmil Constantinescu } 11103ca35412SEmil Constantinescu } 1111f9c1d6abSBarry Smith ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1112be5899b3SLisandro Dalcin /* U <- y(t) + U */ 1113be5899b3SLisandro Dalcin ierr = VecAXPY(U,1,ros->vec_sol_prev);CHKERRQ(ierr); 1114f4aed992SEmil Constantinescu 1115f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1116e27a552bSJed Brown PetscFunctionReturn(0); 1117e27a552bSJed Brown } 1118e27a552bSJed Brown 1119e27a552bSJed Brown /*------------------------------------------------------------*/ 1120b39943a6SLisandro Dalcin 1121b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauReset(TS ts) 1122b39943a6SLisandro Dalcin { 1123b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1124b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1125b39943a6SLisandro Dalcin PetscErrorCode ierr; 1126b39943a6SLisandro Dalcin 1127b39943a6SLisandro Dalcin PetscFunctionBegin; 1128b39943a6SLisandro Dalcin if (!tab) PetscFunctionReturn(0); 1129b39943a6SLisandro Dalcin ierr = VecDestroyVecs(tab->s,&ros->Y);CHKERRQ(ierr); 1130b39943a6SLisandro Dalcin ierr = PetscFree(ros->work);CHKERRQ(ierr); 1131b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1132b39943a6SLisandro Dalcin } 1133b39943a6SLisandro Dalcin 1134e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1135e27a552bSJed Brown { 113661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1137e27a552bSJed Brown PetscErrorCode ierr; 1138e27a552bSJed Brown 1139e27a552bSJed Brown PetscFunctionBegin; 1140b39943a6SLisandro Dalcin ierr = TSRosWTableauReset(ts);CHKERRQ(ierr); 114161692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 114261692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 114361692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 114461692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 1145be5899b3SLisandro Dalcin ierr = VecDestroy(&ros->vec_sol_prev);CHKERRQ(ierr); 1146e27a552bSJed Brown PetscFunctionReturn(0); 1147e27a552bSJed Brown } 1148e27a552bSJed Brown 1149d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1150d5e6173cSPeter Brune { 1151d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1152d5e6173cSPeter Brune PetscErrorCode ierr; 1153d5e6173cSPeter Brune 1154d5e6173cSPeter Brune PetscFunctionBegin; 1155d5e6173cSPeter Brune if (Ydot) { 1156d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1157d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1158d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1159d5e6173cSPeter Brune } 1160d5e6173cSPeter Brune if (Zdot) { 1161d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1162d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1163d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1164d5e6173cSPeter Brune } 1165d5e6173cSPeter Brune if (Ystage) { 1166d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1167d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1168d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1169d5e6173cSPeter Brune } 1170d5e6173cSPeter Brune if (Zstage) { 1171d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1172d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1173d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1174d5e6173cSPeter Brune } 1175d5e6173cSPeter Brune PetscFunctionReturn(0); 1176d5e6173cSPeter Brune } 1177d5e6173cSPeter Brune 1178d5e6173cSPeter Brune 1179d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1180d5e6173cSPeter Brune { 1181d5e6173cSPeter Brune PetscErrorCode ierr; 1182d5e6173cSPeter Brune 1183d5e6173cSPeter Brune PetscFunctionBegin; 1184d5e6173cSPeter Brune if (Ydot) { 1185d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1186d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1187d5e6173cSPeter Brune } 1188d5e6173cSPeter Brune } 1189d5e6173cSPeter Brune if (Zdot) { 1190d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1191d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1192d5e6173cSPeter Brune } 1193d5e6173cSPeter Brune } 1194d5e6173cSPeter Brune if (Ystage) { 1195d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1196d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1197d5e6173cSPeter Brune } 1198d5e6173cSPeter Brune } 1199d5e6173cSPeter Brune if (Zstage) { 1200d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1201d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1202d5e6173cSPeter Brune } 1203d5e6173cSPeter Brune } 1204d5e6173cSPeter Brune PetscFunctionReturn(0); 1205d5e6173cSPeter Brune } 1206d5e6173cSPeter Brune 1207d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1208d5e6173cSPeter Brune { 1209d5e6173cSPeter Brune PetscFunctionBegin; 1210d5e6173cSPeter Brune PetscFunctionReturn(0); 1211d5e6173cSPeter Brune } 1212d5e6173cSPeter Brune 1213d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1214d5e6173cSPeter Brune { 1215d5e6173cSPeter Brune TS ts = (TS)ctx; 1216d5e6173cSPeter Brune PetscErrorCode ierr; 1217d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1218d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1219d5e6173cSPeter Brune 1220d5e6173cSPeter Brune PetscFunctionBegin; 1221d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1222d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1223d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1224d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1225d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1226d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1227d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1228d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1229d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1230d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1231d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1232d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1233d5e6173cSPeter Brune PetscFunctionReturn(0); 1234d5e6173cSPeter Brune } 1235d5e6173cSPeter Brune 1236258e1594SPeter Brune 1237258e1594SPeter Brune static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx) 1238258e1594SPeter Brune { 1239258e1594SPeter Brune PetscFunctionBegin; 1240258e1594SPeter Brune PetscFunctionReturn(0); 1241258e1594SPeter Brune } 1242258e1594SPeter Brune 1243258e1594SPeter Brune static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1244258e1594SPeter Brune { 1245258e1594SPeter Brune TS ts = (TS)ctx; 1246258e1594SPeter Brune PetscErrorCode ierr; 1247258e1594SPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1248258e1594SPeter Brune Vec Ydots,Zdots,Ystages,Zstages; 1249258e1594SPeter Brune 1250258e1594SPeter Brune PetscFunctionBegin; 1251258e1594SPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1252258e1594SPeter Brune ierr = TSRosWGetVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1253258e1594SPeter Brune 1254258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1255258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1256258e1594SPeter Brune 1257258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1258258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1259258e1594SPeter Brune 1260258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1261258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1262258e1594SPeter Brune 1263258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1264258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1265258e1594SPeter Brune 1266258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1267258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1268258e1594SPeter Brune PetscFunctionReturn(0); 1269258e1594SPeter Brune } 1270258e1594SPeter Brune 1271e27a552bSJed Brown /* 1272e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1273e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1274e27a552bSJed Brown */ 1275f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1276e27a552bSJed Brown { 127761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1278e27a552bSJed Brown PetscErrorCode ierr; 1279d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1280b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1281d5e6173cSPeter Brune DM dm,dmsave; 1282e27a552bSJed Brown 1283e27a552bSJed Brown PetscFunctionBegin; 1284d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1285d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1286b296d7d5SJed Brown ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1287f9c1d6abSBarry Smith ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1288d5e6173cSPeter Brune dmsave = ts->dm; 1289d5e6173cSPeter Brune ts->dm = dm; 1290d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1291d5e6173cSPeter Brune ts->dm = dmsave; 1292d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1293e27a552bSJed Brown PetscFunctionReturn(0); 1294e27a552bSJed Brown } 1295e27a552bSJed Brown 1296d1e9a80fSBarry Smith static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat A,Mat B,TS ts) 1297e27a552bSJed Brown { 129861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1299d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1300b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1301e27a552bSJed Brown PetscErrorCode ierr; 1302d5e6173cSPeter Brune DM dm,dmsave; 1303e27a552bSJed Brown 1304e27a552bSJed Brown PetscFunctionBegin; 130561692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1306d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1307d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1308d5e6173cSPeter Brune dmsave = ts->dm; 1309d5e6173cSPeter Brune ts->dm = dm; 1310d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,PETSC_TRUE);CHKERRQ(ierr); 1311d5e6173cSPeter Brune ts->dm = dmsave; 1312d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1313e27a552bSJed Brown PetscFunctionReturn(0); 1314e27a552bSJed Brown } 1315e27a552bSJed Brown 1316b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauSetUp(TS ts) 1317b39943a6SLisandro Dalcin { 1318b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1319b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1320b39943a6SLisandro Dalcin PetscErrorCode ierr; 1321b39943a6SLisandro Dalcin 1322b39943a6SLisandro Dalcin PetscFunctionBegin; 1323b39943a6SLisandro Dalcin ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ros->Y);CHKERRQ(ierr); 1324b39943a6SLisandro Dalcin ierr = PetscMalloc1(tab->s,&ros->work);CHKERRQ(ierr); 1325b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1326b39943a6SLisandro Dalcin } 1327b39943a6SLisandro Dalcin 1328e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1329e27a552bSJed Brown { 133061692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1331e27a552bSJed Brown PetscErrorCode ierr; 1332d5e6173cSPeter Brune DM dm; 1333b39943a6SLisandro Dalcin SNES snes; 1334e27a552bSJed Brown 1335e27a552bSJed Brown PetscFunctionBegin; 1336b39943a6SLisandro Dalcin ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr); 133761692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 133861692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 133961692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 134061692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 1341be5899b3SLisandro Dalcin ierr = VecDuplicate(ts->vec_sol,&ros->vec_sol_prev);CHKERRQ(ierr); 134222d28d08SBarry Smith ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1343d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1344258e1594SPeter Brune ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1345b39943a6SLisandro Dalcin /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1346b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1347b39943a6SLisandro Dalcin if (!((PetscObject)snes)->type_name) { 1348b39943a6SLisandro Dalcin ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1349b39943a6SLisandro Dalcin } 1350e27a552bSJed Brown PetscFunctionReturn(0); 1351e27a552bSJed Brown } 1352e27a552bSJed Brown /*------------------------------------------------------------*/ 1353e27a552bSJed Brown 13544416b707SBarry Smith static PetscErrorCode TSSetFromOptions_RosW(PetscOptionItems *PetscOptionsObject,TS ts) 1355e27a552bSJed Brown { 135661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1357e27a552bSJed Brown PetscErrorCode ierr; 1358b39943a6SLisandro Dalcin SNES snes; 1359e27a552bSJed Brown 1360e27a552bSJed Brown PetscFunctionBegin; 1361e55864a3SBarry Smith ierr = PetscOptionsHead(PetscOptionsObject,"RosW ODE solver options");CHKERRQ(ierr); 1362e27a552bSJed Brown { 136361692a83SJed Brown RosWTableauLink link; 1364e27a552bSJed Brown PetscInt count,choice; 1365e27a552bSJed Brown PetscBool flg; 1366e27a552bSJed Brown const char **namelist; 136761692a83SJed Brown 136861692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1369f489ac74SBarry Smith ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 137061692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1371b39943a6SLisandro Dalcin ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,ros->tableau->name,&choice,&flg);CHKERRQ(ierr); 1372b39943a6SLisandro Dalcin if (flg) {ierr = TSRosWSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1373e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 137461692a83SJed Brown 13750298fd71SBarry Smith ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,NULL);CHKERRQ(ierr); 1376b39943a6SLisandro Dalcin } 1377b39943a6SLisandro Dalcin ierr = PetscOptionsTail();CHKERRQ(ierr); 137861692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 137961692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 138061692a83SJed Brown if (!((PetscObject)snes)->type_name) { 138161692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 138261692a83SJed Brown } 1383e27a552bSJed Brown PetscFunctionReturn(0); 1384e27a552bSJed Brown } 1385e27a552bSJed Brown 1386e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1387e27a552bSJed Brown { 138861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1389e27a552bSJed Brown PetscBool iascii; 1390e27a552bSJed Brown PetscErrorCode ierr; 1391e27a552bSJed Brown 1392e27a552bSJed Brown PetscFunctionBegin; 1393251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1394e27a552bSJed Brown if (iascii) { 13959c334d8fSLisandro Dalcin RosWTableau tab = ros->tableau; 139619fd82e9SBarry Smith TSRosWType rostype; 13979c334d8fSLisandro Dalcin char buf[512]; 1398e408995aSJed Brown PetscInt i; 1399e408995aSJed Brown PetscReal abscissa[512]; 140061692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 140161692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1402de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 140361692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1404e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1405de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1406e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1407e27a552bSJed Brown } 1408e27a552bSJed Brown PetscFunctionReturn(0); 1409e27a552bSJed Brown } 1410e27a552bSJed Brown 14119200755eSBarry Smith static PetscErrorCode TSLoad_RosW(TS ts,PetscViewer viewer) 14129200755eSBarry Smith { 14139200755eSBarry Smith PetscErrorCode ierr; 14149200755eSBarry Smith SNES snes; 14159c334d8fSLisandro Dalcin TSAdapt adapt; 14169200755eSBarry Smith 14179200755eSBarry Smith PetscFunctionBegin; 14189c334d8fSLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 14199c334d8fSLisandro Dalcin ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 14209200755eSBarry Smith ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 14219200755eSBarry Smith ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 14229200755eSBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 14239200755eSBarry Smith ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 14249200755eSBarry Smith ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 14259200755eSBarry Smith PetscFunctionReturn(0); 14269200755eSBarry Smith } 14279200755eSBarry Smith 1428e27a552bSJed Brown /*@C 142961692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1430e27a552bSJed Brown 1431e27a552bSJed Brown Logically collective 1432e27a552bSJed Brown 1433e27a552bSJed Brown Input Parameter: 1434e27a552bSJed Brown + ts - timestepping context 1435b92453a8SLisandro Dalcin - roswtype - type of Rosenbrock-W scheme 1436e27a552bSJed Brown 1437020d8f30SJed Brown Level: beginner 1438e27a552bSJed Brown 1439020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1440e27a552bSJed Brown @*/ 1441b92453a8SLisandro Dalcin PetscErrorCode TSRosWSetType(TS ts,TSRosWType roswtype) 1442e27a552bSJed Brown { 1443e27a552bSJed Brown PetscErrorCode ierr; 1444e27a552bSJed Brown 1445e27a552bSJed Brown PetscFunctionBegin; 1446e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1447b92453a8SLisandro Dalcin PetscValidCharPointer(roswtype,2); 1448b92453a8SLisandro Dalcin ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,roswtype));CHKERRQ(ierr); 1449e27a552bSJed Brown PetscFunctionReturn(0); 1450e27a552bSJed Brown } 1451e27a552bSJed Brown 1452e27a552bSJed Brown /*@C 145361692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1454e27a552bSJed Brown 1455e27a552bSJed Brown Logically collective 1456e27a552bSJed Brown 1457e27a552bSJed Brown Input Parameter: 1458e27a552bSJed Brown . ts - timestepping context 1459e27a552bSJed Brown 1460e27a552bSJed Brown Output Parameter: 146161692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1462e27a552bSJed Brown 1463e27a552bSJed Brown Level: intermediate 1464e27a552bSJed Brown 1465e27a552bSJed Brown .seealso: TSRosWGetType() 1466e27a552bSJed Brown @*/ 146719fd82e9SBarry Smith PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1468e27a552bSJed Brown { 1469e27a552bSJed Brown PetscErrorCode ierr; 1470e27a552bSJed Brown 1471e27a552bSJed Brown PetscFunctionBegin; 1472e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 147319fd82e9SBarry Smith ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1474e27a552bSJed Brown PetscFunctionReturn(0); 1475e27a552bSJed Brown } 1476e27a552bSJed Brown 1477e27a552bSJed Brown /*@C 147861692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1479e27a552bSJed Brown 1480e27a552bSJed Brown Logically collective 1481e27a552bSJed Brown 1482e27a552bSJed Brown Input Parameter: 1483e27a552bSJed Brown + ts - timestepping context 148461692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1485e27a552bSJed Brown 1486e27a552bSJed Brown Level: intermediate 1487e27a552bSJed Brown 1488e27a552bSJed Brown .seealso: TSRosWGetType() 1489e27a552bSJed Brown @*/ 149061692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1491e27a552bSJed Brown { 1492e27a552bSJed Brown PetscErrorCode ierr; 1493e27a552bSJed Brown 1494e27a552bSJed Brown PetscFunctionBegin; 1495e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 149661692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1497e27a552bSJed Brown PetscFunctionReturn(0); 1498e27a552bSJed Brown } 1499e27a552bSJed Brown 1500560360afSLisandro Dalcin static PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1501e27a552bSJed Brown { 150261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1503e27a552bSJed Brown 1504e27a552bSJed Brown PetscFunctionBegin; 150561692a83SJed Brown *rostype = ros->tableau->name; 1506e27a552bSJed Brown PetscFunctionReturn(0); 1507e27a552bSJed Brown } 1508ef20d060SBarry Smith 1509560360afSLisandro Dalcin static PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1510e27a552bSJed Brown { 151161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1512e27a552bSJed Brown PetscErrorCode ierr; 1513e27a552bSJed Brown PetscBool match; 151461692a83SJed Brown RosWTableauLink link; 1515e27a552bSJed Brown 1516e27a552bSJed Brown PetscFunctionBegin; 151761692a83SJed Brown if (ros->tableau) { 151861692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1519e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1520e27a552bSJed Brown } 152161692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 152261692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1523e27a552bSJed Brown if (match) { 1524b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauReset(ts);CHKERRQ(ierr);} 152561692a83SJed Brown ros->tableau = &link->tab; 1526b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr);} 15272ffb9264SLisandro Dalcin ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1528e27a552bSJed Brown PetscFunctionReturn(0); 1529e27a552bSJed Brown } 1530e27a552bSJed Brown } 1531ce94432eSBarry Smith SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1532e27a552bSJed Brown PetscFunctionReturn(0); 1533e27a552bSJed Brown } 153461692a83SJed Brown 1535560360afSLisandro Dalcin static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1536e27a552bSJed Brown { 153761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1538e27a552bSJed Brown 1539e27a552bSJed Brown PetscFunctionBegin; 154061692a83SJed Brown ros->recompute_jacobian = flg; 1541e27a552bSJed Brown PetscFunctionReturn(0); 1542e27a552bSJed Brown } 1543e27a552bSJed Brown 1544b3a6b972SJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1545b3a6b972SJed Brown { 1546b3a6b972SJed Brown PetscErrorCode ierr; 1547b3a6b972SJed Brown 1548b3a6b972SJed Brown PetscFunctionBegin; 1549b3a6b972SJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1550b3a6b972SJed Brown if (ts->dm) { 1551b3a6b972SJed Brown ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1552b3a6b972SJed Brown ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1553b3a6b972SJed Brown } 1554b3a6b972SJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 1555b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",NULL);CHKERRQ(ierr); 1556b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",NULL);CHKERRQ(ierr); 1557b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",NULL);CHKERRQ(ierr); 1558b3a6b972SJed Brown PetscFunctionReturn(0); 1559b3a6b972SJed Brown } 1560d5e6173cSPeter Brune 1561e27a552bSJed Brown /* ------------------------------------------------------------ */ 1562e27a552bSJed Brown /*MC 1563020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1564e27a552bSJed Brown 1565e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1566e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1567e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1568e27a552bSJed Brown 1569e27a552bSJed Brown Notes: 157061692a83SJed Brown This method currently only works with autonomous ODE and DAE. 157161692a83SJed Brown 1572d0685a90SJed Brown Consider trying TSARKIMEX if the stiff part is strongly nonlinear. 1573d0685a90SJed Brown 1574e94cfbe0SPatrick Sanan Developer Notes: 157561692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 157661692a83SJed Brown 1577f9c1d6abSBarry Smith $ udot = f(u) 157861692a83SJed Brown 157961692a83SJed Brown by the stage equations 158061692a83SJed Brown 1581f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 158261692a83SJed Brown 158361692a83SJed Brown and step completion formula 158461692a83SJed Brown 1585f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 158661692a83SJed Brown 1587f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 158861692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 158961692a83SJed Brown we define new variables for the stage equations 159061692a83SJed Brown 159161692a83SJed Brown $ y_i = gamma_ij k_j 159261692a83SJed Brown 159361692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 159461692a83SJed Brown 1595b70472e2SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1} 159661692a83SJed Brown 159761692a83SJed Brown to rewrite the method as 159861692a83SJed Brown 1599f9c1d6abSBarry 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 1600f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j bt_j y_j 160161692a83SJed Brown 160261692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 160361692a83SJed Brown 160461692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 160561692a83SJed Brown 160661692a83SJed Brown or, more compactly in tensor notation 160761692a83SJed Brown 160861692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 160961692a83SJed Brown 161061692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 161161692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 161261692a83SJed Brown equation 161361692a83SJed Brown 1614f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 161561692a83SJed Brown 161661692a83SJed Brown with initial guess y_i = 0. 1617e27a552bSJed Brown 1618e27a552bSJed Brown Level: beginner 1619e27a552bSJed Brown 1620d0685a90SJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSWTHETA1, TSROSWTHETA2, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1621a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1622e27a552bSJed Brown M*/ 16238cc058d9SJed Brown PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1624e27a552bSJed Brown { 162561692a83SJed Brown TS_RosW *ros; 1626e27a552bSJed Brown PetscErrorCode ierr; 1627e27a552bSJed Brown 1628e27a552bSJed Brown PetscFunctionBegin; 1629607a6623SBarry Smith ierr = TSRosWInitializePackage();CHKERRQ(ierr); 1630e27a552bSJed Brown 1631e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1632e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1633e27a552bSJed Brown ts->ops->view = TSView_RosW; 16349200755eSBarry Smith ts->ops->load = TSLoad_RosW; 1635e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1636e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1637e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16381c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 163924655328SShri ts->ops->rollback = TSRollBack_RosW; 1640e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1641e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1642e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1643e27a552bSJed Brown 1644efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 1645efd4aadfSBarry Smith 1646b00a9115SJed Brown ierr = PetscNewLog(ts,&ros);CHKERRQ(ierr); 164761692a83SJed Brown ts->data = (void*)ros; 1648e27a552bSJed Brown 1649bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",TSRosWGetType_RosW);CHKERRQ(ierr); 1650bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",TSRosWSetType_RosW);CHKERRQ(ierr); 1651bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1652b39943a6SLisandro Dalcin 1653b39943a6SLisandro Dalcin ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1654e27a552bSJed Brown PetscFunctionReturn(0); 1655e27a552bSJed Brown } 1656