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 .keywords: TS, TSRosW, register, all 301e27a552bSJed Brown 302e27a552bSJed Brown .seealso: TSRosWRegisterDestroy() 303e27a552bSJed Brown @*/ 304e27a552bSJed Brown PetscErrorCode TSRosWRegisterAll(void) 305e27a552bSJed Brown { 306e27a552bSJed Brown PetscErrorCode ierr; 307e27a552bSJed Brown 308e27a552bSJed Brown PetscFunctionBegin; 309e27a552bSJed Brown if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 310e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_TRUE; 311e27a552bSJed Brown 312e27a552bSJed Brown { 313bbd56ea5SKarl Rupp const PetscReal A = 0; 314bbd56ea5SKarl Rupp const PetscReal Gamma = 1; 315bbd56ea5SKarl Rupp const PetscReal b = 1; 316bbd56ea5SKarl Rupp const PetscReal binterpt=1; 3171f80e275SEmil Constantinescu 3180298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3193606a31eSEmil Constantinescu } 3203606a31eSEmil Constantinescu 3213606a31eSEmil Constantinescu { 322bbd56ea5SKarl Rupp const PetscReal A = 0; 323bbd56ea5SKarl Rupp const PetscReal Gamma = 0.5; 324bbd56ea5SKarl Rupp const PetscReal b = 1; 325bbd56ea5SKarl Rupp const PetscReal binterpt=1; 326bbd56ea5SKarl Rupp 3270298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3283606a31eSEmil Constantinescu } 3293606a31eSEmil Constantinescu 3303606a31eSEmil Constantinescu { 331da80777bSKarl 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. */ 332e27a552bSJed Brown const PetscReal 33361692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 334da80777bSKarl Rupp Gamma[2][2] = {{1.707106781186547524401,0}, {-2.*1.707106781186547524401,1.707106781186547524401}}, 3351c3436cfSJed Brown b[2] = {0.5,0.5}, 3361c3436cfSJed Brown b1[2] = {1.0,0.0}; 3371f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 338da80777bSKarl Rupp binterpt[0][0] = 1.707106781186547524401 - 1.0; 339da80777bSKarl Rupp binterpt[1][0] = 2.0 - 1.707106781186547524401; 340da80777bSKarl Rupp binterpt[0][1] = 1.707106781186547524401 - 1.5; 341da80777bSKarl Rupp binterpt[1][1] = 1.5 - 1.707106781186547524401; 342bbd56ea5SKarl Rupp 3431f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 344e27a552bSJed Brown } 345e27a552bSJed Brown { 346da80777bSKarl 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. */ 347e27a552bSJed Brown const PetscReal 34861692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 349da80777bSKarl Rupp Gamma[2][2] = {{0.2928932188134524755992,0}, {-2.*0.2928932188134524755992,0.2928932188134524755992}}, 3501c3436cfSJed Brown b[2] = {0.5,0.5}, 3511c3436cfSJed Brown b1[2] = {1.0,0.0}; 3521f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 353da80777bSKarl Rupp binterpt[0][0] = 0.2928932188134524755992 - 1.0; 354da80777bSKarl Rupp binterpt[1][0] = 2.0 - 0.2928932188134524755992; 355da80777bSKarl Rupp binterpt[0][1] = 0.2928932188134524755992 - 1.5; 356da80777bSKarl Rupp binterpt[1][1] = 1.5 - 0.2928932188134524755992; 357bbd56ea5SKarl Rupp 3581f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 359fe7e6d57SJed Brown } 360fe7e6d57SJed Brown { 361da80777bSKarl Rupp /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */ 3621f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 363fe7e6d57SJed Brown const PetscReal 364fe7e6d57SJed Brown A[3][3] = {{0,0,0}, 365fe7e6d57SJed Brown {1.5773502691896257e+00,0,0}, 366fe7e6d57SJed Brown {0.5,0,0}}, 367da80777bSKarl Rupp Gamma[3][3] = {{7.8867513459481287e-01,0,0}, 368da80777bSKarl Rupp {-1.5773502691896257e+00,7.8867513459481287e-01,0}, 369da80777bSKarl Rupp {-6.7075317547305480e-01,-1.7075317547305482e-01,7.8867513459481287e-01}}, 370fe7e6d57SJed Brown b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01}, 371fe7e6d57SJed Brown b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01}; 3721f80e275SEmil Constantinescu 3731f80e275SEmil Constantinescu binterpt[0][0] = -0.8094010767585034; 3741f80e275SEmil Constantinescu binterpt[1][0] = -0.5; 3751f80e275SEmil Constantinescu binterpt[2][0] = 2.3094010767585034; 3761f80e275SEmil Constantinescu binterpt[0][1] = 0.9641016151377548; 3771f80e275SEmil Constantinescu binterpt[1][1] = 0.5; 3781f80e275SEmil Constantinescu binterpt[2][1] = -1.4641016151377548; 379bbd56ea5SKarl Rupp 3801f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 381fe7e6d57SJed Brown } 382fe7e6d57SJed Brown { 3833ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 384da80777bSKarl Rupp /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 385fe7e6d57SJed Brown const PetscReal 386fe7e6d57SJed Brown A[4][4] = {{0,0,0,0}, 387fe7e6d57SJed Brown {8.7173304301691801e-01,0,0,0}, 388fe7e6d57SJed Brown {8.4457060015369423e-01,-1.1299064236484185e-01,0,0}, 389fe7e6d57SJed Brown {0,0,1.,0}}, 390da80777bSKarl Rupp Gamma[4][4] = {{4.3586652150845900e-01,0,0,0}, 391da80777bSKarl Rupp {-8.7173304301691801e-01,4.3586652150845900e-01,0,0}, 392da80777bSKarl Rupp {-9.0338057013044082e-01,5.4180672388095326e-02,4.3586652150845900e-01,0}, 393da80777bSKarl Rupp {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,4.3586652150845900e-01}}, 394fe7e6d57SJed Brown b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01}, 3953ca35412SEmil Constantinescu b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01}; 3963ca35412SEmil Constantinescu 3973ca35412SEmil Constantinescu binterpt[0][0]=1.0564298455794094; 3983ca35412SEmil Constantinescu binterpt[1][0]=2.296429974281067; 3993ca35412SEmil Constantinescu binterpt[2][0]=-1.307599564525376; 4003ca35412SEmil Constantinescu binterpt[3][0]=-1.045260255335102; 4013ca35412SEmil Constantinescu binterpt[0][1]=-1.3864882699759573; 4023ca35412SEmil Constantinescu binterpt[1][1]=-8.262611700275677; 4033ca35412SEmil Constantinescu binterpt[2][1]=7.250979895056055; 4043ca35412SEmil Constantinescu binterpt[3][1]=2.398120075195581; 4053ca35412SEmil Constantinescu binterpt[0][2]=0.5721822314575016; 4063ca35412SEmil Constantinescu binterpt[1][2]=4.742931142090097; 4073ca35412SEmil Constantinescu binterpt[2][2]=-4.398120075195578; 4083ca35412SEmil Constantinescu binterpt[3][2]=-0.9169932983520199; 4093ca35412SEmil Constantinescu 4103ca35412SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 411e27a552bSJed Brown } 412ef3c5b88SJed Brown { 413da80777bSKarl Rupp /* const PetscReal g = 0.5; Directly written in-place below */ 414ef3c5b88SJed Brown const PetscReal 415ef3c5b88SJed Brown A[4][4] = {{0,0,0,0}, 416ef3c5b88SJed Brown {0,0,0,0}, 417ef3c5b88SJed Brown {1.,0,0,0}, 418ef3c5b88SJed Brown {0.75,-0.25,0.5,0}}, 419da80777bSKarl Rupp Gamma[4][4] = {{0.5,0,0,0}, 420da80777bSKarl Rupp {1.,0.5,0,0}, 421da80777bSKarl Rupp {-0.25,-0.25,0.5,0}, 422da80777bSKarl Rupp {1./12,1./12,-2./3,0.5}}, 423ef3c5b88SJed Brown b[4] = {5./6,-1./6,-1./6,0.5}, 424ef3c5b88SJed Brown b2[4] = {0.75,-0.25,0.5,0}; 425bbd56ea5SKarl Rupp 4260298fd71SBarry Smith ierr = TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,NULL);CHKERRQ(ierr); 427ef3c5b88SJed Brown } 428ef3c5b88SJed Brown { 429da80777bSKarl Rupp /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */ 430ef3c5b88SJed Brown const PetscReal 431ef3c5b88SJed Brown A[3][3] = {{0,0,0}, 432da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}, 433da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}}, 434da80777bSKarl Rupp Gamma[3][3] = {{0.43586652150845899941601945119356,0,0}, 435da80777bSKarl Rupp {-0.19294655696029095575009695436041,0.43586652150845899941601945119356,0}, 436da80777bSKarl Rupp {0,1.74927148125794685173529749738960,0.43586652150845899941601945119356}}, 437ef3c5b88SJed Brown b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829}, 438ef3c5b88SJed Brown b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619}; 4391f80e275SEmil Constantinescu 4401f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 4411f80e275SEmil Constantinescu binterpt[0][0] = 3.793692883777660870425141387941; 4421f80e275SEmil Constantinescu binterpt[1][0] = -2.918692883777660870425141387941; 4431f80e275SEmil Constantinescu binterpt[2][0] = 0.125; 4441f80e275SEmil Constantinescu binterpt[0][1] = -0.725741064379812106687651020584; 4451f80e275SEmil Constantinescu binterpt[1][1] = 0.559074397713145440020984353917; 4461f80e275SEmil Constantinescu binterpt[2][1] = 0.16666666666666666666666666666667; 4471f80e275SEmil Constantinescu 4481f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 449ef3c5b88SJed Brown } 450b1c69cc3SEmil Constantinescu { 451da80777bSKarl Rupp /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 452da80777bSKarl Rupp * Direct evaluation: s3 = 1.732050807568877293527; 453da80777bSKarl Rupp * g = 0.7886751345948128822546; 454da80777bSKarl Rupp * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */ 455b1c69cc3SEmil Constantinescu const PetscReal 456b1c69cc3SEmil Constantinescu A[3][3] = {{0,0,0}, 457b1c69cc3SEmil Constantinescu {1,0,0}, 458b1c69cc3SEmil Constantinescu {0.25,0.25,0}}, 459b1c69cc3SEmil Constantinescu Gamma[3][3] = {{0,0,0}, 460da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/6.0,0.7886751345948128822546,0}, 461da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/24.0,(-3.0-1.732050807568877293527)/8.0,0.7886751345948128822546}}, 462b1c69cc3SEmil Constantinescu b[3] = {1./6.,1./6.,2./3.}, 463b1c69cc3SEmil Constantinescu b2[3] = {1./4.,1./4.,1./2.}; 464c0cb691aSEmil Constantinescu PetscReal binterpt[3][2]; 465da80777bSKarl Rupp 466c0cb691aSEmil Constantinescu binterpt[0][0]=0.089316397477040902157517886164709; 467c0cb691aSEmil Constantinescu binterpt[1][0]=-0.91068360252295909784248211383529; 468c0cb691aSEmil Constantinescu binterpt[2][0]=1.8213672050459181956849642276706; 469c0cb691aSEmil Constantinescu binterpt[0][1]=0.077350269189625764509148780501957; 470c0cb691aSEmil Constantinescu binterpt[1][1]=1.077350269189625764509148780502; 471c0cb691aSEmil Constantinescu binterpt[2][1]=-1.1547005383792515290182975610039; 472bbd56ea5SKarl Rupp 473c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 474b1c69cc3SEmil Constantinescu } 475b1c69cc3SEmil Constantinescu 476b1c69cc3SEmil Constantinescu { 477b1c69cc3SEmil Constantinescu const PetscReal 478b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 479b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 480b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 481b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 482b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 483b1c69cc3SEmil Constantinescu {0.0,1./4.,0,0}, 484b1c69cc3SEmil Constantinescu {-2.,-2./3.,2./3.,0}, 485b1c69cc3SEmil Constantinescu {1./2.,5./36.,-2./9,0}}, 486b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 487b1c69cc3SEmil Constantinescu b2[4] = {1./8.,3./4.,1./8.,0}; 488c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 489da80777bSKarl Rupp 490c0cb691aSEmil Constantinescu binterpt[0][0]=6.25; 491c0cb691aSEmil Constantinescu binterpt[1][0]=-30.25; 492c0cb691aSEmil Constantinescu binterpt[2][0]=1.75; 493c0cb691aSEmil Constantinescu binterpt[3][0]=23.25; 494c0cb691aSEmil Constantinescu binterpt[0][1]=-9.75; 495c0cb691aSEmil Constantinescu binterpt[1][1]=58.75; 496c0cb691aSEmil Constantinescu binterpt[2][1]=-3.25; 497c0cb691aSEmil Constantinescu binterpt[3][1]=-45.75; 498c0cb691aSEmil Constantinescu binterpt[0][2]=3.6666666666666666666666666666667; 499c0cb691aSEmil Constantinescu binterpt[1][2]=-28.333333333333333333333333333333; 500c0cb691aSEmil Constantinescu binterpt[2][2]=1.6666666666666666666666666666667; 501c0cb691aSEmil Constantinescu binterpt[3][2]=23.; 502bbd56ea5SKarl Rupp 503c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 504b1c69cc3SEmil Constantinescu } 505b1c69cc3SEmil Constantinescu 506b1c69cc3SEmil Constantinescu { 507b1c69cc3SEmil Constantinescu const PetscReal 508b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 509b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 510b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 511b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 512b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 513b1c69cc3SEmil Constantinescu {0.0,3./4.,0,0}, 514b1c69cc3SEmil Constantinescu {-2./3.,-23./9.,2./9.,0}, 515b1c69cc3SEmil Constantinescu {1./18.,65./108.,-2./27,0}}, 516b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 517b1c69cc3SEmil Constantinescu b2[4] = {3./16.,10./16.,3./16.,0}; 518c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 519da80777bSKarl Rupp 520c0cb691aSEmil Constantinescu binterpt[0][0]=1.6911764705882352941176470588235; 521c0cb691aSEmil Constantinescu binterpt[1][0]=3.6813725490196078431372549019608; 522c0cb691aSEmil Constantinescu binterpt[2][0]=0.23039215686274509803921568627451; 523c0cb691aSEmil Constantinescu binterpt[3][0]=-4.6029411764705882352941176470588; 524c0cb691aSEmil Constantinescu binterpt[0][1]=-0.95588235294117647058823529411765; 525c0cb691aSEmil Constantinescu binterpt[1][1]=-6.2401960784313725490196078431373; 526c0cb691aSEmil Constantinescu binterpt[2][1]=-0.31862745098039215686274509803922; 527c0cb691aSEmil Constantinescu binterpt[3][1]=7.5147058823529411764705882352941; 528c0cb691aSEmil Constantinescu binterpt[0][2]=-0.56862745098039215686274509803922; 529c0cb691aSEmil Constantinescu binterpt[1][2]=2.7254901960784313725490196078431; 530c0cb691aSEmil Constantinescu binterpt[2][2]=0.25490196078431372549019607843137; 531c0cb691aSEmil Constantinescu binterpt[3][2]=-2.4117647058823529411764705882353; 532bbd56ea5SKarl Rupp 533c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 534b1c69cc3SEmil Constantinescu } 535753f8adbSEmil Constantinescu 536753f8adbSEmil Constantinescu { 537753f8adbSEmil Constantinescu PetscReal A[4][4],Gamma[4][4],b[4],b2[4]; 5383ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 539753f8adbSEmil Constantinescu 540753f8adbSEmil Constantinescu Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816; 54105e8e825SJed Brown Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0; 542753f8adbSEmil Constantinescu Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476; 543753f8adbSEmil Constantinescu Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816; 54405e8e825SJed Brown Gamma[1][2]=0; Gamma[1][3]=0; 545753f8adbSEmil Constantinescu Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903; 546753f8adbSEmil Constantinescu Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131; 547753f8adbSEmil Constantinescu Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816; 54805e8e825SJed Brown Gamma[2][3]=0; 549753f8adbSEmil Constantinescu Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783; 550753f8adbSEmil Constantinescu Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984; 551753f8adbSEmil Constantinescu Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198; 552753f8adbSEmil Constantinescu Gamma[3][3]=0; 553753f8adbSEmil Constantinescu 55405e8e825SJed Brown A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0; 555753f8adbSEmil Constantinescu A[1][0]=0.8717330430169179988320388950590125027645343373957631; 55605e8e825SJed Brown A[1][1]=0; A[1][2]=0; A[1][3]=0; 557753f8adbSEmil Constantinescu A[2][0]=0.5275890119763004115618079766722914408876108660811028; 558753f8adbSEmil Constantinescu A[2][1]=0.07241098802369958843819203208518599088698057726988732; 55905e8e825SJed Brown A[2][2]=0; A[2][3]=0; 560753f8adbSEmil Constantinescu A[3][0]=0.3990960076760701320627260685975778145384666450351314; 561753f8adbSEmil Constantinescu A[3][1]=-0.4375576546135194437228463747348862825846903771419953; 562753f8adbSEmil Constantinescu A[3][2]=1.038461646937449311660120300601880176655352737312713; 56305e8e825SJed Brown A[3][3]=0; 564753f8adbSEmil Constantinescu 565753f8adbSEmil Constantinescu b[0]=0.1876410243467238251612921333138006734899663569186926; 566753f8adbSEmil Constantinescu b[1]=-0.5952974735769549480478230473706443582188442040780541; 567753f8adbSEmil Constantinescu b[2]=0.9717899277217721234705114616271378792182450260943198; 568753f8adbSEmil Constantinescu b[3]=0.4358665215084589994160194475295062513822671686978816; 569753f8adbSEmil Constantinescu 570753f8adbSEmil Constantinescu b2[0]=0.2147402862233891404862383521089097657790734483804460; 571753f8adbSEmil Constantinescu b2[1]=-0.4851622638849390928209050538171743017757490232519684; 572753f8adbSEmil Constantinescu b2[2]=0.8687250025203875511662123688667549217531982787600080; 573753f8adbSEmil Constantinescu b2[3]=0.4016969751411624011684543450940068201770721128357014; 574753f8adbSEmil Constantinescu 5753ca35412SEmil Constantinescu binterpt[0][0]=2.2565812720167954547104627844105; 5763ca35412SEmil Constantinescu binterpt[1][0]=1.349166413351089573796243820819; 5773ca35412SEmil Constantinescu binterpt[2][0]=-2.4695174540533503758652847586647; 5783ca35412SEmil Constantinescu binterpt[3][0]=-0.13623023131453465264142184656474; 5793ca35412SEmil Constantinescu binterpt[0][1]=-3.0826699111559187902922463354557; 5803ca35412SEmil Constantinescu binterpt[1][1]=-2.4689115685996042534544925650515; 5813ca35412SEmil Constantinescu binterpt[2][1]=5.7428279814696677152129332773553; 5823ca35412SEmil Constantinescu binterpt[3][1]=-0.19124650171414467146619437684812; 5833ca35412SEmil Constantinescu binterpt[0][2]=1.0137296634858471607430756831148; 5843ca35412SEmil Constantinescu binterpt[1][2]=0.52444768167155973161042570784064; 5853ca35412SEmil Constantinescu binterpt[2][2]=-2.3015205996945452158771370439586; 5863ca35412SEmil Constantinescu binterpt[3][2]=0.76334325453713832352363565300308; 587f4aed992SEmil Constantinescu 588f73f8d2cSSatish Balay ierr = TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 589753f8adbSEmil Constantinescu } 59042faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWGRK4T,0.231,PETSC_DEFAULT,PETSC_DEFAULT,0,-0.1282612945269037e+01);CHKERRQ(ierr); 59142faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWSHAMP4,0.5,PETSC_DEFAULT,PETSC_DEFAULT,0,125./108.);CHKERRQ(ierr); 59242faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWVELDD4,0.22570811482256823492,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.355958941201148);CHKERRQ(ierr); 59342faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSW4L,0.57282,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.093502252409163);CHKERRQ(ierr); 594e27a552bSJed Brown PetscFunctionReturn(0); 595e27a552bSJed Brown } 596e27a552bSJed Brown 597d5e6173cSPeter Brune 598d5e6173cSPeter Brune 599e27a552bSJed Brown /*@C 600e27a552bSJed Brown TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 601e27a552bSJed Brown 602e27a552bSJed Brown Not Collective 603e27a552bSJed Brown 604e27a552bSJed Brown Level: advanced 605e27a552bSJed Brown 606e27a552bSJed Brown .keywords: TSRosW, register, destroy 607607a6623SBarry Smith .seealso: TSRosWRegister(), TSRosWRegisterAll() 608e27a552bSJed Brown @*/ 609e27a552bSJed Brown PetscErrorCode TSRosWRegisterDestroy(void) 610e27a552bSJed Brown { 611e27a552bSJed Brown PetscErrorCode ierr; 61261692a83SJed Brown RosWTableauLink link; 613e27a552bSJed Brown 614e27a552bSJed Brown PetscFunctionBegin; 61561692a83SJed Brown while ((link = RosWTableauList)) { 61661692a83SJed Brown RosWTableau t = &link->tab; 61761692a83SJed Brown RosWTableauList = link->next; 61861692a83SJed Brown ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 61943b21953SEmil Constantinescu ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 620fe7e6d57SJed Brown ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 621f4aed992SEmil Constantinescu ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 622e27a552bSJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 623e27a552bSJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 624e27a552bSJed Brown } 625e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 626e27a552bSJed Brown PetscFunctionReturn(0); 627e27a552bSJed Brown } 628e27a552bSJed Brown 629e27a552bSJed Brown /*@C 630e27a552bSJed Brown TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 631e27a552bSJed Brown from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW() 632e27a552bSJed Brown when using static libraries. 633e27a552bSJed Brown 634e27a552bSJed Brown Level: developer 635e27a552bSJed Brown 636e27a552bSJed Brown .keywords: TS, TSRosW, initialize, package 637e27a552bSJed Brown .seealso: PetscInitialize() 638e27a552bSJed Brown @*/ 639607a6623SBarry Smith PetscErrorCode TSRosWInitializePackage(void) 640e27a552bSJed Brown { 641e27a552bSJed Brown PetscErrorCode ierr; 642e27a552bSJed Brown 643e27a552bSJed Brown PetscFunctionBegin; 644e27a552bSJed Brown if (TSRosWPackageInitialized) PetscFunctionReturn(0); 645e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 646e27a552bSJed Brown ierr = TSRosWRegisterAll();CHKERRQ(ierr); 647e27a552bSJed Brown ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 648e27a552bSJed Brown PetscFunctionReturn(0); 649e27a552bSJed Brown } 650e27a552bSJed Brown 651e27a552bSJed Brown /*@C 652e27a552bSJed Brown TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 653e27a552bSJed Brown called from PetscFinalize(). 654e27a552bSJed Brown 655e27a552bSJed Brown Level: developer 656e27a552bSJed Brown 657e27a552bSJed Brown .keywords: Petsc, destroy, package 658e27a552bSJed Brown .seealso: PetscFinalize() 659e27a552bSJed Brown @*/ 660e27a552bSJed Brown PetscErrorCode TSRosWFinalizePackage(void) 661e27a552bSJed Brown { 662e27a552bSJed Brown PetscErrorCode ierr; 663e27a552bSJed Brown 664e27a552bSJed Brown PetscFunctionBegin; 665e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 666e27a552bSJed Brown ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 667e27a552bSJed Brown PetscFunctionReturn(0); 668e27a552bSJed Brown } 669e27a552bSJed Brown 670e27a552bSJed Brown /*@C 67161692a83SJed Brown TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 672e27a552bSJed Brown 673e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 674e27a552bSJed Brown 675e27a552bSJed Brown Input Parameters: 676e27a552bSJed Brown + name - identifier for method 677e27a552bSJed Brown . order - approximation order of method 678e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 67961692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 68061692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 681fe7e6d57SJed Brown . b - Step completion table (dimension s) 6820298fd71SBarry Smith . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available) 683f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 68442faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 685e27a552bSJed Brown 686e27a552bSJed Brown Notes: 68761692a83SJed Brown Several Rosenbrock W methods are provided, this function is only needed to create new methods. 688e27a552bSJed Brown 689e27a552bSJed Brown Level: advanced 690e27a552bSJed Brown 691e27a552bSJed Brown .keywords: TS, register 692e27a552bSJed Brown 693e27a552bSJed Brown .seealso: TSRosW 694e27a552bSJed Brown @*/ 695f9c1d6abSBarry Smith PetscErrorCode TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 696f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 697e27a552bSJed Brown { 698e27a552bSJed Brown PetscErrorCode ierr; 69961692a83SJed Brown RosWTableauLink link; 70061692a83SJed Brown RosWTableau t; 70161692a83SJed Brown PetscInt i,j,k; 70261692a83SJed Brown PetscScalar *GammaInv; 703e27a552bSJed Brown 704e27a552bSJed Brown PetscFunctionBegin; 705fe7e6d57SJed Brown PetscValidCharPointer(name,1); 706fe7e6d57SJed Brown PetscValidPointer(A,4); 707fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 708fe7e6d57SJed Brown PetscValidPointer(b,6); 709fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 710fe7e6d57SJed Brown 7111795a4d1SJed Brown ierr = PetscCalloc1(1,&link);CHKERRQ(ierr); 712e27a552bSJed Brown t = &link->tab; 713e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 714e27a552bSJed Brown t->order = order; 715e27a552bSJed Brown t->s = s; 716dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->A,s*s,&t->Gamma,s,&t->b,s,&t->ASum,s,&t->GammaSum);CHKERRQ(ierr); 717dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->At,s,&t->bt,s*s,&t->GammaInv,s,&t->GammaZeroDiag,s*s,&t->GammaExplicitCorr);CHKERRQ(ierr); 718e27a552bSJed Brown ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 71961692a83SJed Brown ierr = PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 72043b21953SEmil Constantinescu ierr = PetscMemcpy(t->GammaExplicitCorr,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 72161692a83SJed Brown ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 722fe7e6d57SJed Brown if (bembed) { 723dcca6d9dSJed Brown ierr = PetscMalloc2(s,&t->bembed,s,&t->bembedt);CHKERRQ(ierr); 724fe7e6d57SJed Brown ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 725fe7e6d57SJed Brown } 72661692a83SJed Brown for (i=0; i<s; i++) { 72761692a83SJed Brown t->ASum[i] = 0; 72861692a83SJed Brown t->GammaSum[i] = 0; 72961692a83SJed Brown for (j=0; j<s; j++) { 73061692a83SJed Brown t->ASum[i] += A[i*s+j]; 731fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 73261692a83SJed Brown } 73361692a83SJed Brown } 734785e854fSJed Brown ierr = PetscMalloc1(s*s,&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 73561692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 736fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 737fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 738fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 739c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 740fd96d5b0SEmil Constantinescu } else { 741c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 742fd96d5b0SEmil Constantinescu } 743fd96d5b0SEmil Constantinescu } 744fd96d5b0SEmil Constantinescu 74561692a83SJed Brown switch (s) { 74661692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 7472e92ee13SHong Zhang case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7486baedc03SHong Zhang case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7492e92ee13SHong Zhang case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 75061692a83SJed Brown case 5: { 75161692a83SJed Brown PetscInt ipvt5[5]; 75261692a83SJed Brown MatScalar work5[5*5]; 7532e92ee13SHong Zhang ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 75461692a83SJed Brown } 7552e92ee13SHong Zhang case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7562e92ee13SHong Zhang case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 75761692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 75861692a83SJed Brown } 75961692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 76061692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 76143b21953SEmil Constantinescu 76243b21953SEmil Constantinescu for (i=0; i<s; i++) { 76343b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 76443b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 76543b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 76643b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 76743b21953SEmil Constantinescu } 76843b21953SEmil Constantinescu } 76943b21953SEmil Constantinescu } 77043b21953SEmil Constantinescu 77161692a83SJed Brown for (i=0; i<s; i++) { 77261692a83SJed Brown for (j=0; j<s; j++) { 77361692a83SJed Brown t->At[i*s+j] = 0; 77461692a83SJed Brown for (k=0; k<s; k++) { 77561692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 77661692a83SJed Brown } 77761692a83SJed Brown } 77861692a83SJed Brown t->bt[i] = 0; 77961692a83SJed Brown for (j=0; j<s; j++) { 78061692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 78161692a83SJed Brown } 782fe7e6d57SJed Brown if (bembed) { 783fe7e6d57SJed Brown t->bembedt[i] = 0; 784fe7e6d57SJed Brown for (j=0; j<s; j++) { 785fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 786fe7e6d57SJed Brown } 787fe7e6d57SJed Brown } 78861692a83SJed Brown } 7898d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 7908d59e960SJed Brown 791f4aed992SEmil Constantinescu t->pinterp = pinterp; 792785e854fSJed Brown ierr = PetscMalloc1(s*pinterp,&t->binterpt);CHKERRQ(ierr); 7933ca35412SEmil Constantinescu ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 79461692a83SJed Brown link->next = RosWTableauList; 79561692a83SJed Brown RosWTableauList = link; 796e27a552bSJed Brown PetscFunctionReturn(0); 797e27a552bSJed Brown } 798e27a552bSJed Brown 79942faf41dSJed Brown /*@C 80042faf41dSJed Brown TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing paramter choices 80142faf41dSJed Brown 80242faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 80342faf41dSJed Brown 80442faf41dSJed Brown Input Parameters: 80542faf41dSJed Brown + name - identifier for method 80642faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 80742faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 80842faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 80942faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 81042faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 81142faf41dSJed Brown . e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 81242faf41dSJed Brown 81342faf41dSJed Brown Notes: 81442faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 81542faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 81642faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 81742faf41dSJed Brown 81842faf41dSJed Brown Level: developer 81942faf41dSJed Brown 82042faf41dSJed Brown .keywords: TS, register 82142faf41dSJed Brown 82242faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 82342faf41dSJed Brown @*/ 82419fd82e9SBarry Smith PetscErrorCode TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 82542faf41dSJed Brown { 82642faf41dSJed Brown PetscErrorCode ierr; 82742faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 82842faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 82942faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 83042faf41dSJed Brown p42 = one/eight - gamma/three, 83142faf41dSJed Brown p43 = one/twelve - gamma/three, 83242faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 83342faf41dSJed Brown p56 = one/twenty - gamma/four; 83442faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 83542faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 83642faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 83742faf41dSJed Brown 83842faf41dSJed Brown PetscFunctionBegin; 83942faf41dSJed Brown /* Step 1: choose Gamma (input) */ 84042faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 84142faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 84242faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 84342faf41dSJed Brown 84442faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 84542faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 84642faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 84742faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 84842faf41dSJed Brown rhs[0] = one - b3; 84942faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 85042faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 8516baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 85242faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 85342faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 85442faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 85542faf41dSJed Brown 85642faf41dSJed Brown /* Step 3 */ 85742faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 85842faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 85942faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 86042faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 86142faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 86242faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 86342faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 8646baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 86542faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 86642faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 86742faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 86842faf41dSJed Brown 86942faf41dSJed Brown /* Step 4: back-substitute */ 87042faf41dSJed Brown beta32 = beta32beta2p / beta2p; 87142faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 87242faf41dSJed Brown 87342faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 87442faf41dSJed Brown a43 = 0; 87542faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 87642faf41dSJed Brown a42 = a32; 87742faf41dSJed Brown 87842faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 87942faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 88042faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 88142faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 88242faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 88342faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 88442faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 88542faf41dSJed 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; 88642faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 88742faf41dSJed Brown 88842faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 88942faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 89042faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 89142faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 89242faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 89342faf41dSJed Brown 89442faf41dSJed Brown { 89542faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 89642faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 89742faf41dSJed Brown } 8980298fd71SBarry Smith ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,NULL);CHKERRQ(ierr); 89942faf41dSJed Brown PetscFunctionReturn(0); 90042faf41dSJed Brown } 90142faf41dSJed Brown 9021c3436cfSJed Brown /* 9031c3436cfSJed Brown The step completion formula is 9041c3436cfSJed Brown 9051c3436cfSJed Brown x1 = x0 + b^T Y 9061c3436cfSJed Brown 9071c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 9081c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 9091c3436cfSJed Brown 9101c3436cfSJed Brown x1e = x0 + be^T Y 9111c3436cfSJed Brown = x1 - b^T Y + be^T Y 9121c3436cfSJed Brown = x1 + (be - b)^T Y 9131c3436cfSJed Brown 9141c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9151c3436cfSJed Brown */ 916f9c1d6abSBarry Smith static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done) 9171c3436cfSJed Brown { 9181c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9191c3436cfSJed Brown RosWTableau tab = ros->tableau; 9201c3436cfSJed Brown PetscScalar *w = ros->work; 9211c3436cfSJed Brown PetscInt i; 9221c3436cfSJed Brown PetscErrorCode ierr; 9231c3436cfSJed Brown 9241c3436cfSJed Brown PetscFunctionBegin; 9251c3436cfSJed Brown if (order == tab->order) { 926108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 927f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 928de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 929f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 930f9c1d6abSBarry Smith } else {ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr);} 9311c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9321c3436cfSJed Brown PetscFunctionReturn(0); 9331c3436cfSJed Brown } else if (order == tab->order-1) { 9341c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 935108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 936f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 937de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 938f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 939108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 940108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 941f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 942f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 9431c3436cfSJed Brown } 9441c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9451c3436cfSJed Brown PetscFunctionReturn(0); 9461c3436cfSJed Brown } 9471c3436cfSJed Brown unavailable: 9481c3436cfSJed Brown if (done) *done = PETSC_FALSE; 949a7fac7c2SEmil 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); 9501c3436cfSJed Brown PetscFunctionReturn(0); 9511c3436cfSJed Brown } 9521c3436cfSJed Brown 953560360afSLisandro Dalcin static PetscErrorCode TSRollBack_RosW(TS ts) 95424655328SShri { 95524655328SShri TS_RosW *ros = (TS_RosW*)ts->data; 95624655328SShri PetscErrorCode ierr; 95724655328SShri 95824655328SShri PetscFunctionBegin; 959be5899b3SLisandro Dalcin ierr = VecCopy(ros->vec_sol_prev,ts->vec_sol);CHKERRQ(ierr); 96024655328SShri PetscFunctionReturn(0); 96124655328SShri } 96224655328SShri 963e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 964e27a552bSJed Brown { 96561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 96661692a83SJed Brown RosWTableau tab = ros->tableau; 967e27a552bSJed Brown const PetscInt s = tab->s; 9681c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9690feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 970c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 97161692a83SJed Brown PetscScalar *w = ros->work; 9727d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 973e27a552bSJed Brown SNES snes; 9741c3436cfSJed Brown TSAdapt adapt; 975fecfb714SLisandro Dalcin PetscInt i,j,its,lits; 976be5899b3SLisandro Dalcin PetscInt rejections = 0; 977b39943a6SLisandro Dalcin PetscBool stageok,accept = PETSC_TRUE; 978b39943a6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 979e27a552bSJed Brown PetscErrorCode ierr; 980e27a552bSJed Brown 981e27a552bSJed Brown PetscFunctionBegin; 982b39943a6SLisandro Dalcin if (!ts->steprollback) { 983be5899b3SLisandro Dalcin ierr = VecCopy(ts->vec_sol,ros->vec_sol_prev);CHKERRQ(ierr); 984b39943a6SLisandro Dalcin } 985e27a552bSJed Brown 986b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 987b39943a6SLisandro Dalcin while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 9881c3436cfSJed Brown const PetscReal h = ts->time_step; 989e27a552bSJed Brown for (i=0; i<s; i++) { 9901c3436cfSJed Brown ros->stage_time = ts->ptime + h*ASum[i]; 991b8123daeSJed Brown ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 992c17803e7SJed Brown if (GammaZeroDiag[i]) { 993c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 994b296d7d5SJed Brown ros->scoeff = 1.; 995c17803e7SJed Brown } else { 996c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 997b296d7d5SJed Brown ros->scoeff = 1./Gamma[i*s+i]; 998fd96d5b0SEmil Constantinescu } 99961692a83SJed Brown 100061692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 1001de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 1002de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 100361692a83SJed Brown 100461692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 100561692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 100661692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 100761692a83SJed Brown 1008e27a552bSJed Brown /* Initial guess taken from last stage */ 100961692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 101061692a83SJed Brown 10117d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 1012b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 101361692a83SJed Brown if (!ros->recompute_jacobian && !i) { 101461692a83SJed Brown ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 101561692a83SJed Brown } 10160298fd71SBarry Smith ierr = SNESSolve(snes,NULL,Y[i]);CHKERRQ(ierr); 1017e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1018e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10195ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 10207d4bf2deSEmil Constantinescu } else { 10211ce71dffSSatish Balay Mat J,Jp; 10220feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10230feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 102422d28d08SBarry Smith ierr = VecScale(Y[i],-1.0);CHKERRQ(ierr); 10250feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 10260feba352SEmil Constantinescu 10270feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10280feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10290feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 1030fecfb714SLisandro Dalcin 1031fecfb714SLisandro Dalcin /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10320298fd71SBarry Smith ierr = TSGetIJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 1033d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,J,Jp,PETSC_FALSE);CHKERRQ(ierr); 103422d28d08SBarry Smith ierr = MatMult(J,Zstage,Zdot);CHKERRQ(ierr); 10350feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10365ef26d82SJed Brown ts->ksp_its += 1; 1037fecfb714SLisandro Dalcin 1038fecfb714SLisandro Dalcin ierr = VecScale(Y[i],h);CHKERRQ(ierr); 10397d4bf2deSEmil Constantinescu } 10409be3e283SDebojyoti Ghosh ierr = TSPostStage(ts,ros->stage_time,i,Y);CHKERRQ(ierr); 1041fecfb714SLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1042fecfb714SLisandro Dalcin ierr = TSAdaptCheckStage(adapt,ts,ros->stage_time,Y[i],&stageok);CHKERRQ(ierr); 1043fecfb714SLisandro Dalcin if (!stageok) goto reject_step; 1044e27a552bSJed Brown } 1045e27a552bSJed Brown 1046b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 1047b39943a6SLisandro Dalcin ierr = TSEvaluateStep_RosW(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 1048b39943a6SLisandro Dalcin ros->status = TS_STEP_PENDING; 1049552698daSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 10501c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10511917a363SLisandro Dalcin ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 1052fecfb714SLisandro Dalcin ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 1053b39943a6SLisandro Dalcin ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1054b39943a6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 1055b39943a6SLisandro Dalcin ierr = TSRollBack_RosW(ts);CHKERRQ(ierr); 1056be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1057be5899b3SLisandro Dalcin goto reject_step; 1058b39943a6SLisandro Dalcin } 1059b39943a6SLisandro Dalcin 1060e27a552bSJed Brown ts->ptime += ts->time_step; 1061cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 10621c3436cfSJed Brown break; 1063b39943a6SLisandro Dalcin 1064b39943a6SLisandro Dalcin reject_step: 1065fecfb714SLisandro Dalcin ts->reject++; accept = PETSC_FALSE; 1066be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1067b39943a6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 1068be5899b3SLisandro Dalcin ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 10691c3436cfSJed Brown } 10701c3436cfSJed Brown } 1071e27a552bSJed Brown PetscFunctionReturn(0); 1072e27a552bSJed Brown } 1073e27a552bSJed Brown 1074f9c1d6abSBarry Smith static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1075e27a552bSJed Brown { 107661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1077f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1078f4aed992SEmil Constantinescu PetscReal h; 1079f4aed992SEmil Constantinescu PetscReal tt,t; 1080f4aed992SEmil Constantinescu PetscScalar *bt; 1081f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1082f4aed992SEmil Constantinescu PetscErrorCode ierr; 1083f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1084f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1085f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1086e27a552bSJed Brown 1087e27a552bSJed Brown PetscFunctionBegin; 1088ce94432eSBarry Smith if (!Bt) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1089f4aed992SEmil Constantinescu 1090f4aed992SEmil Constantinescu switch (ros->status) { 1091f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1092f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1093f4aed992SEmil Constantinescu h = ts->time_step; 1094f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1095f4aed992SEmil Constantinescu break; 1096f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1097be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1098f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1099f4aed992SEmil Constantinescu break; 1100ce94432eSBarry Smith default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1101f4aed992SEmil Constantinescu } 1102785e854fSJed Brown ierr = PetscMalloc1(s,&bt);CHKERRQ(ierr); 1103f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1104f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1105f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 11063ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1107f4aed992SEmil Constantinescu } 1108f4aed992SEmil Constantinescu } 1109f4aed992SEmil Constantinescu 1110f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1111f9c1d6abSBarry Smith /* U <- 0*/ 1112f9c1d6abSBarry Smith ierr = VecZeroEntries(U);CHKERRQ(ierr); 1113f9c1d6abSBarry Smith /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11143ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j] = 0; 11153ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11163ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11173ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1118f4aed992SEmil Constantinescu } 11193ca35412SEmil Constantinescu } 1120f9c1d6abSBarry Smith ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1121be5899b3SLisandro Dalcin /* U <- y(t) + U */ 1122be5899b3SLisandro Dalcin ierr = VecAXPY(U,1,ros->vec_sol_prev);CHKERRQ(ierr); 1123f4aed992SEmil Constantinescu 1124f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1125e27a552bSJed Brown PetscFunctionReturn(0); 1126e27a552bSJed Brown } 1127e27a552bSJed Brown 1128e27a552bSJed Brown /*------------------------------------------------------------*/ 1129b39943a6SLisandro Dalcin 1130b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauReset(TS ts) 1131b39943a6SLisandro Dalcin { 1132b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1133b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1134b39943a6SLisandro Dalcin PetscErrorCode ierr; 1135b39943a6SLisandro Dalcin 1136b39943a6SLisandro Dalcin PetscFunctionBegin; 1137b39943a6SLisandro Dalcin if (!tab) PetscFunctionReturn(0); 1138b39943a6SLisandro Dalcin ierr = VecDestroyVecs(tab->s,&ros->Y);CHKERRQ(ierr); 1139b39943a6SLisandro Dalcin ierr = PetscFree(ros->work);CHKERRQ(ierr); 1140b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1141b39943a6SLisandro Dalcin } 1142b39943a6SLisandro Dalcin 1143e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1144e27a552bSJed Brown { 114561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1146e27a552bSJed Brown PetscErrorCode ierr; 1147e27a552bSJed Brown 1148e27a552bSJed Brown PetscFunctionBegin; 1149b39943a6SLisandro Dalcin ierr = TSRosWTableauReset(ts);CHKERRQ(ierr); 115061692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 115161692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 115261692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 115361692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 1154be5899b3SLisandro Dalcin ierr = VecDestroy(&ros->vec_sol_prev);CHKERRQ(ierr); 1155e27a552bSJed Brown PetscFunctionReturn(0); 1156e27a552bSJed Brown } 1157e27a552bSJed Brown 1158e27a552bSJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1159e27a552bSJed Brown { 1160e27a552bSJed Brown PetscErrorCode ierr; 1161e27a552bSJed Brown 1162e27a552bSJed Brown PetscFunctionBegin; 1163e27a552bSJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1164e27a552bSJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 1165bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",NULL);CHKERRQ(ierr); 1166bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",NULL);CHKERRQ(ierr); 1167bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",NULL);CHKERRQ(ierr); 1168e27a552bSJed Brown PetscFunctionReturn(0); 1169e27a552bSJed Brown } 1170e27a552bSJed Brown 1171d5e6173cSPeter Brune 1172d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1173d5e6173cSPeter Brune { 1174d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1175d5e6173cSPeter Brune PetscErrorCode ierr; 1176d5e6173cSPeter Brune 1177d5e6173cSPeter Brune PetscFunctionBegin; 1178d5e6173cSPeter Brune if (Ydot) { 1179d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1180d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1181d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1182d5e6173cSPeter Brune } 1183d5e6173cSPeter Brune if (Zdot) { 1184d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1185d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1186d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1187d5e6173cSPeter Brune } 1188d5e6173cSPeter Brune if (Ystage) { 1189d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1190d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1191d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1192d5e6173cSPeter Brune } 1193d5e6173cSPeter Brune if (Zstage) { 1194d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1195d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1196d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1197d5e6173cSPeter Brune } 1198d5e6173cSPeter Brune PetscFunctionReturn(0); 1199d5e6173cSPeter Brune } 1200d5e6173cSPeter Brune 1201d5e6173cSPeter Brune 1202d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1203d5e6173cSPeter Brune { 1204d5e6173cSPeter Brune PetscErrorCode ierr; 1205d5e6173cSPeter Brune 1206d5e6173cSPeter Brune PetscFunctionBegin; 1207d5e6173cSPeter Brune if (Ydot) { 1208d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1209d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1210d5e6173cSPeter Brune } 1211d5e6173cSPeter Brune } 1212d5e6173cSPeter Brune if (Zdot) { 1213d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1214d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1215d5e6173cSPeter Brune } 1216d5e6173cSPeter Brune } 1217d5e6173cSPeter Brune if (Ystage) { 1218d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1219d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1220d5e6173cSPeter Brune } 1221d5e6173cSPeter Brune } 1222d5e6173cSPeter Brune if (Zstage) { 1223d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1224d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1225d5e6173cSPeter Brune } 1226d5e6173cSPeter Brune } 1227d5e6173cSPeter Brune PetscFunctionReturn(0); 1228d5e6173cSPeter Brune } 1229d5e6173cSPeter Brune 1230d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1231d5e6173cSPeter Brune { 1232d5e6173cSPeter Brune PetscFunctionBegin; 1233d5e6173cSPeter Brune PetscFunctionReturn(0); 1234d5e6173cSPeter Brune } 1235d5e6173cSPeter Brune 1236d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1237d5e6173cSPeter Brune { 1238d5e6173cSPeter Brune TS ts = (TS)ctx; 1239d5e6173cSPeter Brune PetscErrorCode ierr; 1240d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1241d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1242d5e6173cSPeter Brune 1243d5e6173cSPeter Brune PetscFunctionBegin; 1244d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1245d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1246d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1247d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1248d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1249d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1250d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1251d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1252d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1253d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1254d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1255d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1256d5e6173cSPeter Brune PetscFunctionReturn(0); 1257d5e6173cSPeter Brune } 1258d5e6173cSPeter Brune 1259258e1594SPeter Brune 1260258e1594SPeter Brune static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx) 1261258e1594SPeter Brune { 1262258e1594SPeter Brune PetscFunctionBegin; 1263258e1594SPeter Brune PetscFunctionReturn(0); 1264258e1594SPeter Brune } 1265258e1594SPeter Brune 1266258e1594SPeter Brune static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1267258e1594SPeter Brune { 1268258e1594SPeter Brune TS ts = (TS)ctx; 1269258e1594SPeter Brune PetscErrorCode ierr; 1270258e1594SPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1271258e1594SPeter Brune Vec Ydots,Zdots,Ystages,Zstages; 1272258e1594SPeter Brune 1273258e1594SPeter Brune PetscFunctionBegin; 1274258e1594SPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1275258e1594SPeter Brune ierr = TSRosWGetVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1276258e1594SPeter Brune 1277258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1278258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1279258e1594SPeter Brune 1280258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1281258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1282258e1594SPeter Brune 1283258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1284258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1285258e1594SPeter Brune 1286258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1287258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1288258e1594SPeter Brune 1289258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1290258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1291258e1594SPeter Brune PetscFunctionReturn(0); 1292258e1594SPeter Brune } 1293258e1594SPeter Brune 1294e27a552bSJed Brown /* 1295e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1296e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1297e27a552bSJed Brown */ 1298f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1299e27a552bSJed Brown { 130061692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1301e27a552bSJed Brown PetscErrorCode ierr; 1302d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1303b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1304d5e6173cSPeter Brune DM dm,dmsave; 1305e27a552bSJed Brown 1306e27a552bSJed Brown PetscFunctionBegin; 1307d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1308d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1309b296d7d5SJed Brown ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1310f9c1d6abSBarry Smith ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1311d5e6173cSPeter Brune dmsave = ts->dm; 1312d5e6173cSPeter Brune ts->dm = dm; 1313d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1314d5e6173cSPeter Brune ts->dm = dmsave; 1315d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1316e27a552bSJed Brown PetscFunctionReturn(0); 1317e27a552bSJed Brown } 1318e27a552bSJed Brown 1319d1e9a80fSBarry Smith static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat A,Mat B,TS ts) 1320e27a552bSJed Brown { 132161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1322d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1323b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1324e27a552bSJed Brown PetscErrorCode ierr; 1325d5e6173cSPeter Brune DM dm,dmsave; 1326e27a552bSJed Brown 1327e27a552bSJed Brown PetscFunctionBegin; 132861692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1329d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1330d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1331d5e6173cSPeter Brune dmsave = ts->dm; 1332d5e6173cSPeter Brune ts->dm = dm; 1333d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,PETSC_TRUE);CHKERRQ(ierr); 1334d5e6173cSPeter Brune ts->dm = dmsave; 1335d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1336e27a552bSJed Brown PetscFunctionReturn(0); 1337e27a552bSJed Brown } 1338e27a552bSJed Brown 1339b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauSetUp(TS ts) 1340b39943a6SLisandro Dalcin { 1341b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1342b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1343b39943a6SLisandro Dalcin PetscErrorCode ierr; 1344b39943a6SLisandro Dalcin 1345b39943a6SLisandro Dalcin PetscFunctionBegin; 1346b39943a6SLisandro Dalcin ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ros->Y);CHKERRQ(ierr); 1347b39943a6SLisandro Dalcin ierr = PetscMalloc1(tab->s,&ros->work);CHKERRQ(ierr); 1348b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1349b39943a6SLisandro Dalcin } 1350b39943a6SLisandro Dalcin 1351e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1352e27a552bSJed Brown { 135361692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1354e27a552bSJed Brown PetscErrorCode ierr; 1355d5e6173cSPeter Brune DM dm; 1356b39943a6SLisandro Dalcin SNES snes; 1357e27a552bSJed Brown 1358e27a552bSJed Brown PetscFunctionBegin; 1359b39943a6SLisandro Dalcin ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr); 136061692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 136161692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 136261692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 136361692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 1364be5899b3SLisandro Dalcin ierr = VecDuplicate(ts->vec_sol,&ros->vec_sol_prev);CHKERRQ(ierr); 136522d28d08SBarry Smith ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1366d5e6173cSPeter Brune if (dm) { 1367d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1368258e1594SPeter Brune ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1369d5e6173cSPeter Brune } 1370b39943a6SLisandro Dalcin /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1371b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1372b39943a6SLisandro Dalcin if (!((PetscObject)snes)->type_name) { 1373b39943a6SLisandro Dalcin ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1374b39943a6SLisandro Dalcin } 1375e27a552bSJed Brown PetscFunctionReturn(0); 1376e27a552bSJed Brown } 1377e27a552bSJed Brown /*------------------------------------------------------------*/ 1378e27a552bSJed Brown 13794416b707SBarry Smith static PetscErrorCode TSSetFromOptions_RosW(PetscOptionItems *PetscOptionsObject,TS ts) 1380e27a552bSJed Brown { 138161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1382e27a552bSJed Brown PetscErrorCode ierr; 1383b39943a6SLisandro Dalcin SNES snes; 1384e27a552bSJed Brown 1385e27a552bSJed Brown PetscFunctionBegin; 1386e55864a3SBarry Smith ierr = PetscOptionsHead(PetscOptionsObject,"RosW ODE solver options");CHKERRQ(ierr); 1387e27a552bSJed Brown { 138861692a83SJed Brown RosWTableauLink link; 1389e27a552bSJed Brown PetscInt count,choice; 1390e27a552bSJed Brown PetscBool flg; 1391e27a552bSJed Brown const char **namelist; 139261692a83SJed Brown 139361692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1394785e854fSJed Brown ierr = PetscMalloc1(count,&namelist);CHKERRQ(ierr); 139561692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1396b39943a6SLisandro Dalcin ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,ros->tableau->name,&choice,&flg);CHKERRQ(ierr); 1397b39943a6SLisandro Dalcin if (flg) {ierr = TSRosWSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1398e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 139961692a83SJed Brown 14000298fd71SBarry Smith ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,NULL);CHKERRQ(ierr); 1401b39943a6SLisandro Dalcin } 1402b39943a6SLisandro Dalcin ierr = PetscOptionsTail();CHKERRQ(ierr); 140361692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 140461692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 140561692a83SJed Brown if (!((PetscObject)snes)->type_name) { 140661692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 140761692a83SJed Brown } 1408e27a552bSJed Brown PetscFunctionReturn(0); 1409e27a552bSJed Brown } 1410e27a552bSJed Brown 1411e27a552bSJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 1412e27a552bSJed Brown { 1413e27a552bSJed Brown PetscErrorCode ierr; 1414e408995aSJed Brown PetscInt i; 1415e408995aSJed Brown size_t left,count; 1416e27a552bSJed Brown char *p; 1417e27a552bSJed Brown 1418e27a552bSJed Brown PetscFunctionBegin; 1419e408995aSJed Brown for (i=0,p=buf,left=len; i<n; i++) { 1420fc6138e5SBarry Smith ierr = PetscSNPrintfCount(p,left,fmt,&count,(double)x[i]);CHKERRQ(ierr); 1421e27a552bSJed Brown if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 1422e27a552bSJed Brown left -= count; 1423e27a552bSJed Brown p += count; 1424e27a552bSJed Brown *p++ = ' '; 1425e27a552bSJed Brown } 1426e27a552bSJed Brown p[i ? 0 : -1] = 0; 1427e27a552bSJed Brown PetscFunctionReturn(0); 1428e27a552bSJed Brown } 1429e27a552bSJed Brown 1430e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1431e27a552bSJed Brown { 143261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1433e27a552bSJed Brown PetscBool iascii; 1434e27a552bSJed Brown PetscErrorCode ierr; 1435e27a552bSJed Brown 1436e27a552bSJed Brown PetscFunctionBegin; 1437251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1438e27a552bSJed Brown if (iascii) { 14399c334d8fSLisandro Dalcin RosWTableau tab = ros->tableau; 144019fd82e9SBarry Smith TSRosWType rostype; 14419c334d8fSLisandro Dalcin char buf[512]; 1442e408995aSJed Brown PetscInt i; 1443e408995aSJed Brown PetscReal abscissa[512]; 144461692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 144561692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1446de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 144761692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1448e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1449de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1450e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1451e27a552bSJed Brown } 1452e27a552bSJed Brown PetscFunctionReturn(0); 1453e27a552bSJed Brown } 1454e27a552bSJed Brown 14559200755eSBarry Smith static PetscErrorCode TSLoad_RosW(TS ts,PetscViewer viewer) 14569200755eSBarry Smith { 14579200755eSBarry Smith PetscErrorCode ierr; 14589200755eSBarry Smith SNES snes; 14599c334d8fSLisandro Dalcin TSAdapt adapt; 14609200755eSBarry Smith 14619200755eSBarry Smith PetscFunctionBegin; 14629c334d8fSLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 14639c334d8fSLisandro Dalcin ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 14649200755eSBarry Smith ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 14659200755eSBarry Smith ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 14669200755eSBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 14679200755eSBarry Smith ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 14689200755eSBarry Smith ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 14699200755eSBarry Smith PetscFunctionReturn(0); 14709200755eSBarry Smith } 14719200755eSBarry Smith 1472e27a552bSJed Brown /*@C 147361692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1474e27a552bSJed Brown 1475e27a552bSJed Brown Logically collective 1476e27a552bSJed Brown 1477e27a552bSJed Brown Input Parameter: 1478e27a552bSJed Brown + ts - timestepping context 1479b92453a8SLisandro Dalcin - roswtype - type of Rosenbrock-W scheme 1480e27a552bSJed Brown 1481020d8f30SJed Brown Level: beginner 1482e27a552bSJed Brown 1483020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1484e27a552bSJed Brown @*/ 1485b92453a8SLisandro Dalcin PetscErrorCode TSRosWSetType(TS ts,TSRosWType roswtype) 1486e27a552bSJed Brown { 1487e27a552bSJed Brown PetscErrorCode ierr; 1488e27a552bSJed Brown 1489e27a552bSJed Brown PetscFunctionBegin; 1490e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1491b92453a8SLisandro Dalcin PetscValidCharPointer(roswtype,2); 1492b92453a8SLisandro Dalcin ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,roswtype));CHKERRQ(ierr); 1493e27a552bSJed Brown PetscFunctionReturn(0); 1494e27a552bSJed Brown } 1495e27a552bSJed Brown 1496e27a552bSJed Brown /*@C 149761692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1498e27a552bSJed Brown 1499e27a552bSJed Brown Logically collective 1500e27a552bSJed Brown 1501e27a552bSJed Brown Input Parameter: 1502e27a552bSJed Brown . ts - timestepping context 1503e27a552bSJed Brown 1504e27a552bSJed Brown Output Parameter: 150561692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1506e27a552bSJed Brown 1507e27a552bSJed Brown Level: intermediate 1508e27a552bSJed Brown 1509e27a552bSJed Brown .seealso: TSRosWGetType() 1510e27a552bSJed Brown @*/ 151119fd82e9SBarry Smith PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1512e27a552bSJed Brown { 1513e27a552bSJed Brown PetscErrorCode ierr; 1514e27a552bSJed Brown 1515e27a552bSJed Brown PetscFunctionBegin; 1516e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 151719fd82e9SBarry Smith ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1518e27a552bSJed Brown PetscFunctionReturn(0); 1519e27a552bSJed Brown } 1520e27a552bSJed Brown 1521e27a552bSJed Brown /*@C 152261692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1523e27a552bSJed Brown 1524e27a552bSJed Brown Logically collective 1525e27a552bSJed Brown 1526e27a552bSJed Brown Input Parameter: 1527e27a552bSJed Brown + ts - timestepping context 152861692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1529e27a552bSJed Brown 1530e27a552bSJed Brown Level: intermediate 1531e27a552bSJed Brown 1532e27a552bSJed Brown .seealso: TSRosWGetType() 1533e27a552bSJed Brown @*/ 153461692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1535e27a552bSJed Brown { 1536e27a552bSJed Brown PetscErrorCode ierr; 1537e27a552bSJed Brown 1538e27a552bSJed Brown PetscFunctionBegin; 1539e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 154061692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1541e27a552bSJed Brown PetscFunctionReturn(0); 1542e27a552bSJed Brown } 1543e27a552bSJed Brown 1544560360afSLisandro Dalcin static PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1545e27a552bSJed Brown { 154661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1547e27a552bSJed Brown 1548e27a552bSJed Brown PetscFunctionBegin; 154961692a83SJed Brown *rostype = ros->tableau->name; 1550e27a552bSJed Brown PetscFunctionReturn(0); 1551e27a552bSJed Brown } 1552ef20d060SBarry Smith 1553560360afSLisandro Dalcin static PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1554e27a552bSJed Brown { 155561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1556e27a552bSJed Brown PetscErrorCode ierr; 1557e27a552bSJed Brown PetscBool match; 155861692a83SJed Brown RosWTableauLink link; 1559e27a552bSJed Brown 1560e27a552bSJed Brown PetscFunctionBegin; 156161692a83SJed Brown if (ros->tableau) { 156261692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1563e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1564e27a552bSJed Brown } 156561692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 156661692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1567e27a552bSJed Brown if (match) { 1568b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauReset(ts);CHKERRQ(ierr);} 156961692a83SJed Brown ros->tableau = &link->tab; 1570b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr);} 15712ffb9264SLisandro Dalcin ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1572e27a552bSJed Brown PetscFunctionReturn(0); 1573e27a552bSJed Brown } 1574e27a552bSJed Brown } 1575ce94432eSBarry Smith SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1576e27a552bSJed Brown PetscFunctionReturn(0); 1577e27a552bSJed Brown } 157861692a83SJed Brown 1579560360afSLisandro Dalcin static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1580e27a552bSJed Brown { 158161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1582e27a552bSJed Brown 1583e27a552bSJed Brown PetscFunctionBegin; 158461692a83SJed Brown ros->recompute_jacobian = flg; 1585e27a552bSJed Brown PetscFunctionReturn(0); 1586e27a552bSJed Brown } 1587e27a552bSJed Brown 1588d5e6173cSPeter Brune 1589e27a552bSJed Brown /* ------------------------------------------------------------ */ 1590e27a552bSJed Brown /*MC 1591020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1592e27a552bSJed Brown 1593e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1594e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1595e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1596e27a552bSJed Brown 1597e27a552bSJed Brown Notes: 159861692a83SJed Brown This method currently only works with autonomous ODE and DAE. 159961692a83SJed Brown 1600d0685a90SJed Brown Consider trying TSARKIMEX if the stiff part is strongly nonlinear. 1601d0685a90SJed Brown 160261692a83SJed Brown Developer notes: 160361692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 160461692a83SJed Brown 1605f9c1d6abSBarry Smith $ udot = f(u) 160661692a83SJed Brown 160761692a83SJed Brown by the stage equations 160861692a83SJed Brown 1609f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 161061692a83SJed Brown 161161692a83SJed Brown and step completion formula 161261692a83SJed Brown 1613f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 161461692a83SJed Brown 1615f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 161661692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 161761692a83SJed Brown we define new variables for the stage equations 161861692a83SJed Brown 161961692a83SJed Brown $ y_i = gamma_ij k_j 162061692a83SJed Brown 162161692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 162261692a83SJed Brown 1623b70472e2SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1} 162461692a83SJed Brown 162561692a83SJed Brown to rewrite the method as 162661692a83SJed Brown 1627f9c1d6abSBarry 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 1628f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j bt_j y_j 162961692a83SJed Brown 163061692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 163161692a83SJed Brown 163261692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 163361692a83SJed Brown 163461692a83SJed Brown or, more compactly in tensor notation 163561692a83SJed Brown 163661692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 163761692a83SJed Brown 163861692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 163961692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 164061692a83SJed Brown equation 164161692a83SJed Brown 1642f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 164361692a83SJed Brown 164461692a83SJed Brown with initial guess y_i = 0. 1645e27a552bSJed Brown 1646e27a552bSJed Brown Level: beginner 1647e27a552bSJed Brown 1648d0685a90SJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSWTHETA1, TSROSWTHETA2, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1649a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1650e27a552bSJed Brown M*/ 16518cc058d9SJed Brown PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1652e27a552bSJed Brown { 165361692a83SJed Brown TS_RosW *ros; 1654e27a552bSJed Brown PetscErrorCode ierr; 1655e27a552bSJed Brown 1656e27a552bSJed Brown PetscFunctionBegin; 1657607a6623SBarry Smith ierr = TSRosWInitializePackage();CHKERRQ(ierr); 1658e27a552bSJed Brown 1659e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1660e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1661e27a552bSJed Brown ts->ops->view = TSView_RosW; 16629200755eSBarry Smith ts->ops->load = TSLoad_RosW; 1663e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1664e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1665e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16661c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 166724655328SShri ts->ops->rollback = TSRollBack_RosW; 1668e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1669e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1670e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1671e27a552bSJed Brown 1672*efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 1673*efd4aadfSBarry Smith 1674b00a9115SJed Brown ierr = PetscNewLog(ts,&ros);CHKERRQ(ierr); 167561692a83SJed Brown ts->data = (void*)ros; 1676e27a552bSJed Brown 1677bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",TSRosWGetType_RosW);CHKERRQ(ierr); 1678bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",TSRosWSetType_RosW);CHKERRQ(ierr); 1679bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1680b39943a6SLisandro Dalcin 1681b39943a6SLisandro Dalcin ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1682e27a552bSJed Brown PetscFunctionReturn(0); 1683e27a552bSJed Brown } 1684