1e27a552bSJed Brown /* 261692a83SJed Brown Code for timestepping with Rosenbrock W methods 3e27a552bSJed Brown 4e27a552bSJed Brown Notes: 5e27a552bSJed Brown The general system is written as 6e27a552bSJed Brown 7f9c1d6abSBarry Smith F(t,U,Udot) = G(t,U) 8e27a552bSJed Brown 9f9c1d6abSBarry Smith where F represents the stiff part of the physics and G represents the non-stiff part. 10f9c1d6abSBarry Smith This method is designed to be linearly implicit on F and can use an approximate and lagged Jacobian. 11e27a552bSJed Brown 12e27a552bSJed Brown */ 13af0996ceSBarry Smith #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 141e25c274SJed Brown #include <petscdm.h> 15e27a552bSJed Brown 16af0996ceSBarry Smith #include <petsc/private/kernels/blockinvert.h> 1761692a83SJed Brown 1819fd82e9SBarry Smith static TSRosWType TSRosWDefault = TSROSWRA34PW2; 19e27a552bSJed Brown static PetscBool TSRosWRegisterAllCalled; 20e27a552bSJed Brown static PetscBool TSRosWPackageInitialized; 21e27a552bSJed Brown 2261692a83SJed Brown typedef struct _RosWTableau *RosWTableau; 2361692a83SJed Brown struct _RosWTableau { 24e27a552bSJed Brown char *name; 25e27a552bSJed Brown PetscInt order; /* Classical approximation order of the method */ 26e27a552bSJed Brown PetscInt s; /* Number of stages */ 27f4aed992SEmil Constantinescu PetscInt pinterp; /* Interpolation order */ 2861692a83SJed Brown PetscReal *A; /* Propagation table, strictly lower triangular */ 2961692a83SJed Brown PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */ 30c17803e7SJed Brown PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */ 3143b21953SEmil Constantinescu PetscReal *GammaExplicitCorr; /* Coefficients for correction terms needed for explicit stages in transformed variables*/ 3261692a83SJed Brown PetscReal *b; /* Step completion table */ 33fe7e6d57SJed Brown PetscReal *bembed; /* Step completion table for embedded method of order one less */ 3461692a83SJed Brown PetscReal *ASum; /* Row sum of A */ 3561692a83SJed Brown PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */ 3661692a83SJed Brown PetscReal *At; /* Propagation table in transformed variables */ 3761692a83SJed Brown PetscReal *bt; /* Step completion table in transformed variables */ 38fe7e6d57SJed Brown PetscReal *bembedt; /* Step completion table of order one less in transformed variables */ 3961692a83SJed Brown PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */ 408d59e960SJed Brown PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 413ca35412SEmil Constantinescu PetscReal *binterpt; /* Dense output formula */ 42e27a552bSJed Brown }; 4361692a83SJed Brown typedef struct _RosWTableauLink *RosWTableauLink; 4461692a83SJed Brown struct _RosWTableauLink { 4561692a83SJed Brown struct _RosWTableau tab; 4661692a83SJed Brown RosWTableauLink next; 47e27a552bSJed Brown }; 4861692a83SJed Brown static RosWTableauLink RosWTableauList; 49e27a552bSJed Brown 50e27a552bSJed Brown typedef struct { 5161692a83SJed Brown RosWTableau tableau; 5261692a83SJed Brown Vec *Y; /* States computed during the step, used to complete the step */ 53e27a552bSJed Brown Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 5461692a83SJed Brown Vec Ystage; /* Work vector for the state value at each stage */ 5561692a83SJed Brown Vec Zdot; /* Ydot = Zdot + shift*Y */ 5661692a83SJed Brown Vec Zstage; /* Y = Zstage + Y */ 57be5899b3SLisandro Dalcin Vec vec_sol_prev; /* Solution from the previous step (used for interpolation and rollback)*/ 581c3436cfSJed Brown PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 59b296d7d5SJed Brown PetscReal scoeff; /* shift = scoeff/dt */ 60e27a552bSJed Brown PetscReal stage_time; 61c17803e7SJed Brown PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */ 6261692a83SJed Brown PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 63108c343cSJed Brown TSStepStatus status; 64e27a552bSJed Brown } TS_RosW; 65e27a552bSJed Brown 66fe7e6d57SJed Brown /*MC 673606a31eSEmil Constantinescu TSROSWTHETA1 - One stage first order L-stable Rosenbrock-W scheme (aka theta method). 683606a31eSEmil Constantinescu 693606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 703606a31eSEmil Constantinescu 713606a31eSEmil Constantinescu Level: intermediate 723606a31eSEmil Constantinescu 733606a31eSEmil Constantinescu .seealso: TSROSW 743606a31eSEmil Constantinescu M*/ 753606a31eSEmil Constantinescu 763606a31eSEmil Constantinescu /*MC 773606a31eSEmil Constantinescu TSROSWTHETA2 - One stage second order A-stable Rosenbrock-W scheme (aka theta method). 783606a31eSEmil Constantinescu 793606a31eSEmil Constantinescu Only an approximate Jacobian is needed. 803606a31eSEmil Constantinescu 813606a31eSEmil Constantinescu Level: intermediate 823606a31eSEmil Constantinescu 833606a31eSEmil Constantinescu .seealso: TSROSW 843606a31eSEmil Constantinescu M*/ 853606a31eSEmil Constantinescu 863606a31eSEmil Constantinescu /*MC 87fe7e6d57SJed Brown TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme. 88fe7e6d57SJed Brown 89fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P. 90fe7e6d57SJed Brown 91fe7e6d57SJed Brown Level: intermediate 92fe7e6d57SJed Brown 93fe7e6d57SJed Brown .seealso: TSROSW 94fe7e6d57SJed Brown M*/ 95fe7e6d57SJed Brown 96fe7e6d57SJed Brown /*MC 97fe7e6d57SJed Brown TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme. 98fe7e6d57SJed Brown 99fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M. 100fe7e6d57SJed Brown 101fe7e6d57SJed Brown Level: intermediate 102fe7e6d57SJed Brown 103fe7e6d57SJed Brown .seealso: TSROSW 104fe7e6d57SJed Brown M*/ 105fe7e6d57SJed Brown 106fe7e6d57SJed Brown /*MC 107fe7e6d57SJed Brown TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1. 108fe7e6d57SJed Brown 109fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 110fe7e6d57SJed Brown 111fe7e6d57SJed Brown This is strongly A-stable with R(infty) = 0.73. The embedded method of order 2 is strongly A-stable with R(infty) = 0.73. 112fe7e6d57SJed Brown 113fe7e6d57SJed Brown References: 11496a0c994SBarry Smith . 1. - Rang and Angermann, New Rosenbrock W methods of order 3 for partial differential algebraic equations of index 1, 2005. 115fe7e6d57SJed Brown 116fe7e6d57SJed Brown Level: intermediate 117fe7e6d57SJed Brown 118fe7e6d57SJed Brown .seealso: TSROSW 119fe7e6d57SJed Brown M*/ 120fe7e6d57SJed Brown 121fe7e6d57SJed Brown /*MC 122fe7e6d57SJed Brown TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1. 123fe7e6d57SJed Brown 124fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 125fe7e6d57SJed Brown 126fe7e6d57SJed Brown This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48. 127fe7e6d57SJed Brown 128fe7e6d57SJed Brown References: 12996a0c994SBarry Smith . 1. - Rang and Angermann, New Rosenbrock W methods of order 3 for partial differential algebraic equations of index 1, 2005. 130fe7e6d57SJed Brown 131fe7e6d57SJed Brown Level: intermediate 132fe7e6d57SJed Brown 133fe7e6d57SJed Brown .seealso: TSROSW 134fe7e6d57SJed Brown M*/ 135fe7e6d57SJed Brown 136ef3c5b88SJed Brown /*MC 137ef3c5b88SJed Brown TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme 138ef3c5b88SJed Brown 139ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 140ef3c5b88SJed Brown 141ef3c5b88SJed Brown Both the third order and embedded second order methods are stiffly accurate and L-stable. 142ef3c5b88SJed Brown 143ef3c5b88SJed Brown References: 14496a0c994SBarry Smith . 1. - Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 145ef3c5b88SJed Brown 146ef3c5b88SJed Brown Level: intermediate 147ef3c5b88SJed Brown 148ef3c5b88SJed Brown .seealso: TSROSW, TSROSWSANDU3 149ef3c5b88SJed Brown M*/ 150ef3c5b88SJed Brown 151ef3c5b88SJed Brown /*MC 152ef3c5b88SJed Brown TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme 153ef3c5b88SJed Brown 154ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 155ef3c5b88SJed Brown 156ef3c5b88SJed Brown The third order method is L-stable, but not stiffly accurate. 157ef3c5b88SJed Brown The second order embedded method is strongly A-stable with R(infty) = 0.5. 158ef3c5b88SJed Brown The internal stages are L-stable. 159ef3c5b88SJed Brown This method is called ROS3 in the paper. 160ef3c5b88SJed Brown 161ef3c5b88SJed Brown References: 16296a0c994SBarry Smith . 1. - Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 163ef3c5b88SJed Brown 164ef3c5b88SJed Brown Level: intermediate 165ef3c5b88SJed Brown 166ef3c5b88SJed Brown .seealso: TSROSW, TSROSWRODAS3 167ef3c5b88SJed Brown M*/ 168ef3c5b88SJed Brown 169961f28d0SJed Brown /*MC 170961f28d0SJed Brown TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 171961f28d0SJed Brown 172961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 173961f28d0SJed Brown 174961f28d0SJed Brown A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 175961f28d0SJed Brown 176961f28d0SJed Brown References: 17796a0c994SBarry Smith . Emil Constantinescu 178961f28d0SJed Brown 179961f28d0SJed Brown Level: intermediate 180961f28d0SJed Brown 18143b21953SEmil Constantinescu .seealso: TSROSW, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, SSP 182961f28d0SJed Brown M*/ 183961f28d0SJed Brown 184961f28d0SJed Brown /*MC 185998eb97aSJed Brown TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 186961f28d0SJed Brown 187961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 188961f28d0SJed Brown 189961f28d0SJed Brown L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 190961f28d0SJed Brown 191961f28d0SJed Brown References: 19296a0c994SBarry Smith . Emil Constantinescu 193961f28d0SJed Brown 194961f28d0SJed Brown Level: intermediate 195961f28d0SJed Brown 19643b21953SEmil Constantinescu .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLLSSP3P4S2C, TSSSP 197961f28d0SJed Brown M*/ 198961f28d0SJed Brown 199961f28d0SJed Brown /*MC 200998eb97aSJed Brown TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 201961f28d0SJed Brown 202961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 203961f28d0SJed Brown 204961f28d0SJed Brown L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 205961f28d0SJed Brown 206961f28d0SJed Brown References: 20796a0c994SBarry Smith . Emil Constantinescu 208961f28d0SJed Brown 209961f28d0SJed Brown Level: intermediate 210961f28d0SJed Brown 211961f28d0SJed Brown .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSSSP 212961f28d0SJed Brown M*/ 213961f28d0SJed Brown 21442faf41dSJed Brown /*MC 21542faf41dSJed Brown TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop 21642faf41dSJed Brown 21742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 21842faf41dSJed Brown 21942faf41dSJed Brown A(89.3 degrees)-stable, |R(infty)| = 0.454. 22042faf41dSJed Brown 22142faf41dSJed Brown This method does not provide a dense output formula. 22242faf41dSJed Brown 22342faf41dSJed Brown References: 22496a0c994SBarry Smith + 1. - Kaps and Rentrop, Generalized Runge Kutta methods of order four with stepsize control for stiff ordinary differential equations, 1979. 22596a0c994SBarry Smith - 2. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 22642faf41dSJed Brown 22742faf41dSJed Brown Hairer's code ros4.f 22842faf41dSJed Brown 22942faf41dSJed Brown Level: intermediate 23042faf41dSJed Brown 23142faf41dSJed Brown .seealso: TSROSW, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 23242faf41dSJed Brown M*/ 23342faf41dSJed Brown 23442faf41dSJed Brown /*MC 23542faf41dSJed Brown TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine 23642faf41dSJed Brown 23742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 23842faf41dSJed Brown 23942faf41dSJed Brown A-stable, |R(infty)| = 1/3. 24042faf41dSJed Brown 24142faf41dSJed Brown This method does not provide a dense output formula. 24242faf41dSJed Brown 24342faf41dSJed Brown References: 24496a0c994SBarry Smith + 1. - Shampine, Implementation of Rosenbrock methods, 1982. 24596a0c994SBarry Smith - 2. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 24642faf41dSJed Brown 24742faf41dSJed Brown Hairer's code ros4.f 24842faf41dSJed Brown 24942faf41dSJed Brown Level: intermediate 25042faf41dSJed Brown 25142faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWVELDD4, TSROSW4L 25242faf41dSJed Brown M*/ 25342faf41dSJed Brown 25442faf41dSJed Brown /*MC 25542faf41dSJed Brown TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen 25642faf41dSJed Brown 25742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 25842faf41dSJed Brown 25942faf41dSJed Brown A(89.5 degrees)-stable, |R(infty)| = 0.24. 26042faf41dSJed Brown 26142faf41dSJed Brown This method does not provide a dense output formula. 26242faf41dSJed Brown 26342faf41dSJed Brown References: 26496a0c994SBarry Smith + 1. - van Veldhuizen, D stability and Kaps Rentrop methods, 1984. 26596a0c994SBarry Smith - 2. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 26642faf41dSJed Brown 26742faf41dSJed Brown Hairer's code ros4.f 26842faf41dSJed Brown 26942faf41dSJed Brown Level: intermediate 27042faf41dSJed Brown 27142faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 27242faf41dSJed Brown M*/ 27342faf41dSJed Brown 27442faf41dSJed Brown /*MC 27542faf41dSJed Brown TSROSW4L - four stage, fourth order Rosenbrock (not W) method 27642faf41dSJed Brown 27742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 27842faf41dSJed Brown 27942faf41dSJed Brown A-stable and L-stable 28042faf41dSJed Brown 28142faf41dSJed Brown This method does not provide a dense output formula. 28242faf41dSJed Brown 28342faf41dSJed Brown References: 28496a0c994SBarry Smith . 1. - Hairer and Wanner, Solving Ordinary Differential Equations II, Section 4 Table 7.2. 28542faf41dSJed Brown 28642faf41dSJed Brown Hairer's code ros4.f 28742faf41dSJed Brown 28842faf41dSJed Brown Level: intermediate 28942faf41dSJed Brown 29042faf41dSJed Brown .seealso: TSROSW, TSROSWGRK4T, TSROSWSHAMP4, TSROSW4L 29142faf41dSJed Brown M*/ 29242faf41dSJed Brown 293e27a552bSJed Brown /*@C 294be5899b3SLisandro Dalcin TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in TSRosW 295e27a552bSJed Brown 296e27a552bSJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 297e27a552bSJed Brown 298e27a552bSJed Brown Level: advanced 299e27a552bSJed Brown 300e27a552bSJed Brown .seealso: TSRosWRegisterDestroy() 301e27a552bSJed Brown @*/ 302e27a552bSJed Brown PetscErrorCode TSRosWRegisterAll(void) 303e27a552bSJed Brown { 304e27a552bSJed Brown PetscErrorCode ierr; 305e27a552bSJed Brown 306e27a552bSJed Brown PetscFunctionBegin; 307e27a552bSJed Brown if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 308e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_TRUE; 309e27a552bSJed Brown 310e27a552bSJed Brown { 311bbd56ea5SKarl Rupp const PetscReal A = 0; 312bbd56ea5SKarl Rupp const PetscReal Gamma = 1; 313bbd56ea5SKarl Rupp const PetscReal b = 1; 314bbd56ea5SKarl Rupp const PetscReal binterpt=1; 3151f80e275SEmil Constantinescu 3160298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA1,1,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3173606a31eSEmil Constantinescu } 3183606a31eSEmil Constantinescu 3193606a31eSEmil Constantinescu { 320bbd56ea5SKarl Rupp const PetscReal A = 0; 321bbd56ea5SKarl Rupp const PetscReal Gamma = 0.5; 322bbd56ea5SKarl Rupp const PetscReal b = 1; 323bbd56ea5SKarl Rupp const PetscReal binterpt=1; 324bbd56ea5SKarl Rupp 3250298fd71SBarry Smith ierr = TSRosWRegister(TSROSWTHETA2,2,1,&A,&Gamma,&b,NULL,1,&binterpt);CHKERRQ(ierr); 3263606a31eSEmil Constantinescu } 3273606a31eSEmil Constantinescu 3283606a31eSEmil Constantinescu { 329da80777bSKarl Rupp /*const PetscReal g = 1. + 1./PetscSqrtReal(2.0); Direct evaluation: 1.707106781186547524401. Used for setting up arrays of values known at compile time below. */ 330e27a552bSJed Brown const PetscReal 33161692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 332da80777bSKarl Rupp Gamma[2][2] = {{1.707106781186547524401,0}, {-2.*1.707106781186547524401,1.707106781186547524401}}, 3331c3436cfSJed Brown b[2] = {0.5,0.5}, 3341c3436cfSJed Brown b1[2] = {1.0,0.0}; 3351f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 336da80777bSKarl Rupp binterpt[0][0] = 1.707106781186547524401 - 1.0; 337da80777bSKarl Rupp binterpt[1][0] = 2.0 - 1.707106781186547524401; 338da80777bSKarl Rupp binterpt[0][1] = 1.707106781186547524401 - 1.5; 339da80777bSKarl Rupp binterpt[1][1] = 1.5 - 1.707106781186547524401; 340bbd56ea5SKarl Rupp 3411f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 342e27a552bSJed Brown } 343e27a552bSJed Brown { 344da80777bSKarl Rupp /*const PetscReal g = 1. - 1./PetscSqrtReal(2.0); Direct evaluation: 0.2928932188134524755992. Used for setting up arrays of values known at compile time below. */ 345e27a552bSJed Brown const PetscReal 34661692a83SJed Brown A[2][2] = {{0,0}, {1.,0}}, 347da80777bSKarl Rupp Gamma[2][2] = {{0.2928932188134524755992,0}, {-2.*0.2928932188134524755992,0.2928932188134524755992}}, 3481c3436cfSJed Brown b[2] = {0.5,0.5}, 3491c3436cfSJed Brown b1[2] = {1.0,0.0}; 3501f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 351da80777bSKarl Rupp binterpt[0][0] = 0.2928932188134524755992 - 1.0; 352da80777bSKarl Rupp binterpt[1][0] = 2.0 - 0.2928932188134524755992; 353da80777bSKarl Rupp binterpt[0][1] = 0.2928932188134524755992 - 1.5; 354da80777bSKarl Rupp binterpt[1][1] = 1.5 - 0.2928932188134524755992; 355bbd56ea5SKarl Rupp 3561f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1,2,&binterpt[0][0]);CHKERRQ(ierr); 357fe7e6d57SJed Brown } 358fe7e6d57SJed Brown { 359da80777bSKarl Rupp /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */ 3601f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 361fe7e6d57SJed Brown const PetscReal 362fe7e6d57SJed Brown A[3][3] = {{0,0,0}, 363fe7e6d57SJed Brown {1.5773502691896257e+00,0,0}, 364fe7e6d57SJed Brown {0.5,0,0}}, 365da80777bSKarl Rupp Gamma[3][3] = {{7.8867513459481287e-01,0,0}, 366da80777bSKarl Rupp {-1.5773502691896257e+00,7.8867513459481287e-01,0}, 367da80777bSKarl Rupp {-6.7075317547305480e-01,-1.7075317547305482e-01,7.8867513459481287e-01}}, 368fe7e6d57SJed Brown b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01}, 369fe7e6d57SJed Brown b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01}; 3701f80e275SEmil Constantinescu 3711f80e275SEmil Constantinescu binterpt[0][0] = -0.8094010767585034; 3721f80e275SEmil Constantinescu binterpt[1][0] = -0.5; 3731f80e275SEmil Constantinescu binterpt[2][0] = 2.3094010767585034; 3741f80e275SEmil Constantinescu binterpt[0][1] = 0.9641016151377548; 3751f80e275SEmil Constantinescu binterpt[1][1] = 0.5; 3761f80e275SEmil Constantinescu binterpt[2][1] = -1.4641016151377548; 377bbd56ea5SKarl Rupp 3781f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 379fe7e6d57SJed Brown } 380fe7e6d57SJed Brown { 3813ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 382da80777bSKarl Rupp /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 383fe7e6d57SJed Brown const PetscReal 384fe7e6d57SJed Brown A[4][4] = {{0,0,0,0}, 385fe7e6d57SJed Brown {8.7173304301691801e-01,0,0,0}, 386fe7e6d57SJed Brown {8.4457060015369423e-01,-1.1299064236484185e-01,0,0}, 387fe7e6d57SJed Brown {0,0,1.,0}}, 388da80777bSKarl Rupp Gamma[4][4] = {{4.3586652150845900e-01,0,0,0}, 389da80777bSKarl Rupp {-8.7173304301691801e-01,4.3586652150845900e-01,0,0}, 390da80777bSKarl Rupp {-9.0338057013044082e-01,5.4180672388095326e-02,4.3586652150845900e-01,0}, 391da80777bSKarl Rupp {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,4.3586652150845900e-01}}, 392fe7e6d57SJed Brown b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01}, 3933ca35412SEmil Constantinescu b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01}; 3943ca35412SEmil Constantinescu 3953ca35412SEmil Constantinescu binterpt[0][0]=1.0564298455794094; 3963ca35412SEmil Constantinescu binterpt[1][0]=2.296429974281067; 3973ca35412SEmil Constantinescu binterpt[2][0]=-1.307599564525376; 3983ca35412SEmil Constantinescu binterpt[3][0]=-1.045260255335102; 3993ca35412SEmil Constantinescu binterpt[0][1]=-1.3864882699759573; 4003ca35412SEmil Constantinescu binterpt[1][1]=-8.262611700275677; 4013ca35412SEmil Constantinescu binterpt[2][1]=7.250979895056055; 4023ca35412SEmil Constantinescu binterpt[3][1]=2.398120075195581; 4033ca35412SEmil Constantinescu binterpt[0][2]=0.5721822314575016; 4043ca35412SEmil Constantinescu binterpt[1][2]=4.742931142090097; 4053ca35412SEmil Constantinescu binterpt[2][2]=-4.398120075195578; 4063ca35412SEmil Constantinescu binterpt[3][2]=-0.9169932983520199; 4073ca35412SEmil Constantinescu 4083ca35412SEmil Constantinescu ierr = TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 409e27a552bSJed Brown } 410ef3c5b88SJed Brown { 411da80777bSKarl Rupp /* const PetscReal g = 0.5; Directly written in-place below */ 412ef3c5b88SJed Brown const PetscReal 413ef3c5b88SJed Brown A[4][4] = {{0,0,0,0}, 414ef3c5b88SJed Brown {0,0,0,0}, 415ef3c5b88SJed Brown {1.,0,0,0}, 416ef3c5b88SJed Brown {0.75,-0.25,0.5,0}}, 417da80777bSKarl Rupp Gamma[4][4] = {{0.5,0,0,0}, 418da80777bSKarl Rupp {1.,0.5,0,0}, 419da80777bSKarl Rupp {-0.25,-0.25,0.5,0}, 420da80777bSKarl Rupp {1./12,1./12,-2./3,0.5}}, 421ef3c5b88SJed Brown b[4] = {5./6,-1./6,-1./6,0.5}, 422ef3c5b88SJed Brown b2[4] = {0.75,-0.25,0.5,0}; 423bbd56ea5SKarl Rupp 4240298fd71SBarry Smith ierr = TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2,0,NULL);CHKERRQ(ierr); 425ef3c5b88SJed Brown } 426ef3c5b88SJed Brown { 427da80777bSKarl Rupp /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */ 428ef3c5b88SJed Brown const PetscReal 429ef3c5b88SJed Brown A[3][3] = {{0,0,0}, 430da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}, 431da80777bSKarl Rupp {0.43586652150845899941601945119356,0,0}}, 432da80777bSKarl Rupp Gamma[3][3] = {{0.43586652150845899941601945119356,0,0}, 433da80777bSKarl Rupp {-0.19294655696029095575009695436041,0.43586652150845899941601945119356,0}, 434da80777bSKarl Rupp {0,1.74927148125794685173529749738960,0.43586652150845899941601945119356}}, 435ef3c5b88SJed Brown b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829}, 436ef3c5b88SJed Brown b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619}; 4371f80e275SEmil Constantinescu 4381f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 4391f80e275SEmil Constantinescu binterpt[0][0] = 3.793692883777660870425141387941; 4401f80e275SEmil Constantinescu binterpt[1][0] = -2.918692883777660870425141387941; 4411f80e275SEmil Constantinescu binterpt[2][0] = 0.125; 4421f80e275SEmil Constantinescu binterpt[0][1] = -0.725741064379812106687651020584; 4431f80e275SEmil Constantinescu binterpt[1][1] = 0.559074397713145440020984353917; 4441f80e275SEmil Constantinescu binterpt[2][1] = 0.16666666666666666666666666666667; 4451f80e275SEmil Constantinescu 4461f80e275SEmil Constantinescu ierr = TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 447ef3c5b88SJed Brown } 448b1c69cc3SEmil Constantinescu { 449da80777bSKarl Rupp /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 450da80777bSKarl Rupp * Direct evaluation: s3 = 1.732050807568877293527; 451da80777bSKarl Rupp * g = 0.7886751345948128822546; 452da80777bSKarl Rupp * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */ 453b1c69cc3SEmil Constantinescu const PetscReal 454b1c69cc3SEmil Constantinescu A[3][3] = {{0,0,0}, 455b1c69cc3SEmil Constantinescu {1,0,0}, 456b1c69cc3SEmil Constantinescu {0.25,0.25,0}}, 457b1c69cc3SEmil Constantinescu Gamma[3][3] = {{0,0,0}, 458da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/6.0,0.7886751345948128822546,0}, 459da80777bSKarl Rupp {(-3.0-1.732050807568877293527)/24.0,(-3.0-1.732050807568877293527)/8.0,0.7886751345948128822546}}, 460b1c69cc3SEmil Constantinescu b[3] = {1./6.,1./6.,2./3.}, 461b1c69cc3SEmil Constantinescu b2[3] = {1./4.,1./4.,1./2.}; 462c0cb691aSEmil Constantinescu PetscReal binterpt[3][2]; 463da80777bSKarl Rupp 464c0cb691aSEmil Constantinescu binterpt[0][0]=0.089316397477040902157517886164709; 465c0cb691aSEmil Constantinescu binterpt[1][0]=-0.91068360252295909784248211383529; 466c0cb691aSEmil Constantinescu binterpt[2][0]=1.8213672050459181956849642276706; 467c0cb691aSEmil Constantinescu binterpt[0][1]=0.077350269189625764509148780501957; 468c0cb691aSEmil Constantinescu binterpt[1][1]=1.077350269189625764509148780502; 469c0cb691aSEmil Constantinescu binterpt[2][1]=-1.1547005383792515290182975610039; 470bbd56ea5SKarl Rupp 471c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2,2,&binterpt[0][0]);CHKERRQ(ierr); 472b1c69cc3SEmil Constantinescu } 473b1c69cc3SEmil Constantinescu 474b1c69cc3SEmil Constantinescu { 475b1c69cc3SEmil Constantinescu const PetscReal 476b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 477b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 478b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 479b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 480b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 481b1c69cc3SEmil Constantinescu {0.0,1./4.,0,0}, 482b1c69cc3SEmil Constantinescu {-2.,-2./3.,2./3.,0}, 483b1c69cc3SEmil Constantinescu {1./2.,5./36.,-2./9,0}}, 484b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 485b1c69cc3SEmil Constantinescu b2[4] = {1./8.,3./4.,1./8.,0}; 486c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 487da80777bSKarl Rupp 488c0cb691aSEmil Constantinescu binterpt[0][0]=6.25; 489c0cb691aSEmil Constantinescu binterpt[1][0]=-30.25; 490c0cb691aSEmil Constantinescu binterpt[2][0]=1.75; 491c0cb691aSEmil Constantinescu binterpt[3][0]=23.25; 492c0cb691aSEmil Constantinescu binterpt[0][1]=-9.75; 493c0cb691aSEmil Constantinescu binterpt[1][1]=58.75; 494c0cb691aSEmil Constantinescu binterpt[2][1]=-3.25; 495c0cb691aSEmil Constantinescu binterpt[3][1]=-45.75; 496c0cb691aSEmil Constantinescu binterpt[0][2]=3.6666666666666666666666666666667; 497c0cb691aSEmil Constantinescu binterpt[1][2]=-28.333333333333333333333333333333; 498c0cb691aSEmil Constantinescu binterpt[2][2]=1.6666666666666666666666666666667; 499c0cb691aSEmil Constantinescu binterpt[3][2]=23.; 500bbd56ea5SKarl Rupp 501c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 502b1c69cc3SEmil Constantinescu } 503b1c69cc3SEmil Constantinescu 504b1c69cc3SEmil Constantinescu { 505b1c69cc3SEmil Constantinescu const PetscReal 506b1c69cc3SEmil Constantinescu A[4][4] = {{0,0,0,0}, 507b1c69cc3SEmil Constantinescu {1./2.,0,0,0}, 508b1c69cc3SEmil Constantinescu {1./2.,1./2.,0,0}, 509b1c69cc3SEmil Constantinescu {1./6.,1./6.,1./6.,0}}, 510b1c69cc3SEmil Constantinescu Gamma[4][4] = {{1./2.,0,0,0}, 511b1c69cc3SEmil Constantinescu {0.0,3./4.,0,0}, 512b1c69cc3SEmil Constantinescu {-2./3.,-23./9.,2./9.,0}, 513b1c69cc3SEmil Constantinescu {1./18.,65./108.,-2./27,0}}, 514b1c69cc3SEmil Constantinescu b[4] = {1./6.,1./6.,1./6.,1./2.}, 515b1c69cc3SEmil Constantinescu b2[4] = {3./16.,10./16.,3./16.,0}; 516c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 517da80777bSKarl Rupp 518c0cb691aSEmil Constantinescu binterpt[0][0]=1.6911764705882352941176470588235; 519c0cb691aSEmil Constantinescu binterpt[1][0]=3.6813725490196078431372549019608; 520c0cb691aSEmil Constantinescu binterpt[2][0]=0.23039215686274509803921568627451; 521c0cb691aSEmil Constantinescu binterpt[3][0]=-4.6029411764705882352941176470588; 522c0cb691aSEmil Constantinescu binterpt[0][1]=-0.95588235294117647058823529411765; 523c0cb691aSEmil Constantinescu binterpt[1][1]=-6.2401960784313725490196078431373; 524c0cb691aSEmil Constantinescu binterpt[2][1]=-0.31862745098039215686274509803922; 525c0cb691aSEmil Constantinescu binterpt[3][1]=7.5147058823529411764705882352941; 526c0cb691aSEmil Constantinescu binterpt[0][2]=-0.56862745098039215686274509803922; 527c0cb691aSEmil Constantinescu binterpt[1][2]=2.7254901960784313725490196078431; 528c0cb691aSEmil Constantinescu binterpt[2][2]=0.25490196078431372549019607843137; 529c0cb691aSEmil Constantinescu binterpt[3][2]=-2.4117647058823529411764705882353; 530bbd56ea5SKarl Rupp 531c0cb691aSEmil Constantinescu ierr = TSRosWRegister(TSROSWLLSSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 532b1c69cc3SEmil Constantinescu } 533753f8adbSEmil Constantinescu 534753f8adbSEmil Constantinescu { 535753f8adbSEmil Constantinescu PetscReal A[4][4],Gamma[4][4],b[4],b2[4]; 5363ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 537753f8adbSEmil Constantinescu 538753f8adbSEmil Constantinescu Gamma[0][0]=0.4358665215084589994160194475295062513822671686978816; 53905e8e825SJed Brown Gamma[0][1]=0; Gamma[0][2]=0; Gamma[0][3]=0; 540753f8adbSEmil Constantinescu Gamma[1][0]=-1.997527830934941248426324674704153457289527280554476; 541753f8adbSEmil Constantinescu Gamma[1][1]=0.4358665215084589994160194475295062513822671686978816; 54205e8e825SJed Brown Gamma[1][2]=0; Gamma[1][3]=0; 543753f8adbSEmil Constantinescu Gamma[2][0]=-1.007948511795029620852002345345404191008352770119903; 544753f8adbSEmil Constantinescu Gamma[2][1]=-0.004648958462629345562774289390054679806993396798458131; 545753f8adbSEmil Constantinescu Gamma[2][2]=0.4358665215084589994160194475295062513822671686978816; 54605e8e825SJed Brown Gamma[2][3]=0; 547753f8adbSEmil Constantinescu Gamma[3][0]=-0.6685429734233467180451604600279552604364311322650783; 548753f8adbSEmil Constantinescu Gamma[3][1]=0.6056625986449338476089525334450053439525178740492984; 549753f8adbSEmil Constantinescu Gamma[3][2]=-0.9717899277217721234705114616271378792182450260943198; 550753f8adbSEmil Constantinescu Gamma[3][3]=0; 551753f8adbSEmil Constantinescu 55205e8e825SJed Brown A[0][0]=0; A[0][1]=0; A[0][2]=0; A[0][3]=0; 553753f8adbSEmil Constantinescu A[1][0]=0.8717330430169179988320388950590125027645343373957631; 55405e8e825SJed Brown A[1][1]=0; A[1][2]=0; A[1][3]=0; 555753f8adbSEmil Constantinescu A[2][0]=0.5275890119763004115618079766722914408876108660811028; 556753f8adbSEmil Constantinescu A[2][1]=0.07241098802369958843819203208518599088698057726988732; 55705e8e825SJed Brown A[2][2]=0; A[2][3]=0; 558753f8adbSEmil Constantinescu A[3][0]=0.3990960076760701320627260685975778145384666450351314; 559753f8adbSEmil Constantinescu A[3][1]=-0.4375576546135194437228463747348862825846903771419953; 560753f8adbSEmil Constantinescu A[3][2]=1.038461646937449311660120300601880176655352737312713; 56105e8e825SJed Brown A[3][3]=0; 562753f8adbSEmil Constantinescu 563753f8adbSEmil Constantinescu b[0]=0.1876410243467238251612921333138006734899663569186926; 564753f8adbSEmil Constantinescu b[1]=-0.5952974735769549480478230473706443582188442040780541; 565753f8adbSEmil Constantinescu b[2]=0.9717899277217721234705114616271378792182450260943198; 566753f8adbSEmil Constantinescu b[3]=0.4358665215084589994160194475295062513822671686978816; 567753f8adbSEmil Constantinescu 568753f8adbSEmil Constantinescu b2[0]=0.2147402862233891404862383521089097657790734483804460; 569753f8adbSEmil Constantinescu b2[1]=-0.4851622638849390928209050538171743017757490232519684; 570753f8adbSEmil Constantinescu b2[2]=0.8687250025203875511662123688667549217531982787600080; 571753f8adbSEmil Constantinescu b2[3]=0.4016969751411624011684543450940068201770721128357014; 572753f8adbSEmil Constantinescu 5733ca35412SEmil Constantinescu binterpt[0][0]=2.2565812720167954547104627844105; 5743ca35412SEmil Constantinescu binterpt[1][0]=1.349166413351089573796243820819; 5753ca35412SEmil Constantinescu binterpt[2][0]=-2.4695174540533503758652847586647; 5763ca35412SEmil Constantinescu binterpt[3][0]=-0.13623023131453465264142184656474; 5773ca35412SEmil Constantinescu binterpt[0][1]=-3.0826699111559187902922463354557; 5783ca35412SEmil Constantinescu binterpt[1][1]=-2.4689115685996042534544925650515; 5793ca35412SEmil Constantinescu binterpt[2][1]=5.7428279814696677152129332773553; 5803ca35412SEmil Constantinescu binterpt[3][1]=-0.19124650171414467146619437684812; 5813ca35412SEmil Constantinescu binterpt[0][2]=1.0137296634858471607430756831148; 5823ca35412SEmil Constantinescu binterpt[1][2]=0.52444768167155973161042570784064; 5833ca35412SEmil Constantinescu binterpt[2][2]=-2.3015205996945452158771370439586; 5843ca35412SEmil Constantinescu binterpt[3][2]=0.76334325453713832352363565300308; 585f4aed992SEmil Constantinescu 586f73f8d2cSSatish Balay ierr = TSRosWRegister(TSROSWARK3,3,4,&A[0][0],&Gamma[0][0],b,b2,3,&binterpt[0][0]);CHKERRQ(ierr); 587753f8adbSEmil Constantinescu } 58842faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWGRK4T,0.231,PETSC_DEFAULT,PETSC_DEFAULT,0,-0.1282612945269037e+01);CHKERRQ(ierr); 58942faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWSHAMP4,0.5,PETSC_DEFAULT,PETSC_DEFAULT,0,125./108.);CHKERRQ(ierr); 59042faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSWVELDD4,0.22570811482256823492,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.355958941201148);CHKERRQ(ierr); 59142faf41dSJed Brown ierr = TSRosWRegisterRos4(TSROSW4L,0.57282,PETSC_DEFAULT,PETSC_DEFAULT,0,-1.093502252409163);CHKERRQ(ierr); 592e27a552bSJed Brown PetscFunctionReturn(0); 593e27a552bSJed Brown } 594e27a552bSJed Brown 595d5e6173cSPeter Brune 596d5e6173cSPeter Brune 597e27a552bSJed Brown /*@C 598e27a552bSJed Brown TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 599e27a552bSJed Brown 600e27a552bSJed Brown Not Collective 601e27a552bSJed Brown 602e27a552bSJed Brown Level: advanced 603e27a552bSJed Brown 604607a6623SBarry Smith .seealso: TSRosWRegister(), TSRosWRegisterAll() 605e27a552bSJed Brown @*/ 606e27a552bSJed Brown PetscErrorCode TSRosWRegisterDestroy(void) 607e27a552bSJed Brown { 608e27a552bSJed Brown PetscErrorCode ierr; 60961692a83SJed Brown RosWTableauLink link; 610e27a552bSJed Brown 611e27a552bSJed Brown PetscFunctionBegin; 61261692a83SJed Brown while ((link = RosWTableauList)) { 61361692a83SJed Brown RosWTableau t = &link->tab; 61461692a83SJed Brown RosWTableauList = link->next; 61561692a83SJed Brown ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 61643b21953SEmil Constantinescu ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 617fe7e6d57SJed Brown ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 618f4aed992SEmil Constantinescu ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 619e27a552bSJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 620e27a552bSJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 621e27a552bSJed Brown } 622e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 623e27a552bSJed Brown PetscFunctionReturn(0); 624e27a552bSJed Brown } 625e27a552bSJed Brown 626e27a552bSJed Brown /*@C 627e27a552bSJed Brown TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 6288a690491SBarry Smith from TSInitializePackage(). 629e27a552bSJed Brown 630e27a552bSJed Brown Level: developer 631e27a552bSJed Brown 632e27a552bSJed Brown .seealso: PetscInitialize() 633e27a552bSJed Brown @*/ 634607a6623SBarry Smith PetscErrorCode TSRosWInitializePackage(void) 635e27a552bSJed Brown { 636e27a552bSJed Brown PetscErrorCode ierr; 637e27a552bSJed Brown 638e27a552bSJed Brown PetscFunctionBegin; 639e27a552bSJed Brown if (TSRosWPackageInitialized) PetscFunctionReturn(0); 640e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 641e27a552bSJed Brown ierr = TSRosWRegisterAll();CHKERRQ(ierr); 642e27a552bSJed Brown ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 643e27a552bSJed Brown PetscFunctionReturn(0); 644e27a552bSJed Brown } 645e27a552bSJed Brown 646e27a552bSJed Brown /*@C 647e27a552bSJed Brown TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 648e27a552bSJed Brown called from PetscFinalize(). 649e27a552bSJed Brown 650e27a552bSJed Brown Level: developer 651e27a552bSJed Brown 652e27a552bSJed Brown .seealso: PetscFinalize() 653e27a552bSJed Brown @*/ 654e27a552bSJed Brown PetscErrorCode TSRosWFinalizePackage(void) 655e27a552bSJed Brown { 656e27a552bSJed Brown PetscErrorCode ierr; 657e27a552bSJed Brown 658e27a552bSJed Brown PetscFunctionBegin; 659e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 660e27a552bSJed Brown ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 661e27a552bSJed Brown PetscFunctionReturn(0); 662e27a552bSJed Brown } 663e27a552bSJed Brown 664e27a552bSJed Brown /*@C 66561692a83SJed Brown TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 666e27a552bSJed Brown 667e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 668e27a552bSJed Brown 669e27a552bSJed Brown Input Parameters: 670e27a552bSJed Brown + name - identifier for method 671e27a552bSJed Brown . order - approximation order of method 672e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 67361692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 67461692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 675fe7e6d57SJed Brown . b - Step completion table (dimension s) 6760298fd71SBarry Smith . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available) 677f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 67842faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 679e27a552bSJed Brown 680e27a552bSJed Brown Notes: 68161692a83SJed Brown Several Rosenbrock W methods are provided, this function is only needed to create new methods. 682e27a552bSJed Brown 683e27a552bSJed Brown Level: advanced 684e27a552bSJed Brown 685e27a552bSJed Brown .seealso: TSRosW 686e27a552bSJed Brown @*/ 687f9c1d6abSBarry Smith PetscErrorCode TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 688f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 689e27a552bSJed Brown { 690e27a552bSJed Brown PetscErrorCode ierr; 69161692a83SJed Brown RosWTableauLink link; 69261692a83SJed Brown RosWTableau t; 69361692a83SJed Brown PetscInt i,j,k; 69461692a83SJed Brown PetscScalar *GammaInv; 695e27a552bSJed Brown 696e27a552bSJed Brown PetscFunctionBegin; 697fe7e6d57SJed Brown PetscValidCharPointer(name,1); 698fe7e6d57SJed Brown PetscValidPointer(A,4); 699fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 700fe7e6d57SJed Brown PetscValidPointer(b,6); 701fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 702fe7e6d57SJed Brown 7031d36bdfdSBarry Smith ierr = TSRosWInitializePackage();CHKERRQ(ierr); 704071fcb05SBarry Smith ierr = PetscNew(&link);CHKERRQ(ierr); 705e27a552bSJed Brown t = &link->tab; 706e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 707e27a552bSJed Brown t->order = order; 708e27a552bSJed Brown t->s = s; 709dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->A,s*s,&t->Gamma,s,&t->b,s,&t->ASum,s,&t->GammaSum);CHKERRQ(ierr); 710dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->At,s,&t->bt,s*s,&t->GammaInv,s,&t->GammaZeroDiag,s*s,&t->GammaExplicitCorr);CHKERRQ(ierr); 711580bdb30SBarry Smith ierr = PetscArraycpy(t->A,A,s*s);CHKERRQ(ierr); 712580bdb30SBarry Smith ierr = PetscArraycpy(t->Gamma,Gamma,s*s);CHKERRQ(ierr); 713580bdb30SBarry Smith ierr = PetscArraycpy(t->GammaExplicitCorr,Gamma,s*s);CHKERRQ(ierr); 714580bdb30SBarry Smith ierr = PetscArraycpy(t->b,b,s);CHKERRQ(ierr); 715fe7e6d57SJed Brown if (bembed) { 716dcca6d9dSJed Brown ierr = PetscMalloc2(s,&t->bembed,s,&t->bembedt);CHKERRQ(ierr); 717580bdb30SBarry Smith ierr = PetscArraycpy(t->bembed,bembed,s);CHKERRQ(ierr); 718fe7e6d57SJed Brown } 71961692a83SJed Brown for (i=0; i<s; i++) { 72061692a83SJed Brown t->ASum[i] = 0; 72161692a83SJed Brown t->GammaSum[i] = 0; 72261692a83SJed Brown for (j=0; j<s; j++) { 72361692a83SJed Brown t->ASum[i] += A[i*s+j]; 724fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 72561692a83SJed Brown } 72661692a83SJed Brown } 727785e854fSJed Brown ierr = PetscMalloc1(s*s,&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 72861692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 729fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 730fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 731fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 732c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 733fd96d5b0SEmil Constantinescu } else { 734c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 735fd96d5b0SEmil Constantinescu } 736fd96d5b0SEmil Constantinescu } 737fd96d5b0SEmil Constantinescu 73861692a83SJed Brown switch (s) { 73961692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 7402e92ee13SHong Zhang case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7416baedc03SHong Zhang case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7422e92ee13SHong Zhang case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 74361692a83SJed Brown case 5: { 74461692a83SJed Brown PetscInt ipvt5[5]; 74561692a83SJed Brown MatScalar work5[5*5]; 7462e92ee13SHong Zhang ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 74761692a83SJed Brown } 7482e92ee13SHong Zhang case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7492e92ee13SHong Zhang case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 75061692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 75161692a83SJed Brown } 75261692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 75361692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 75443b21953SEmil Constantinescu 75543b21953SEmil Constantinescu for (i=0; i<s; i++) { 75643b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 75743b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 75843b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 75943b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 76043b21953SEmil Constantinescu } 76143b21953SEmil Constantinescu } 76243b21953SEmil Constantinescu } 76343b21953SEmil Constantinescu 76461692a83SJed Brown for (i=0; i<s; i++) { 76561692a83SJed Brown for (j=0; j<s; j++) { 76661692a83SJed Brown t->At[i*s+j] = 0; 76761692a83SJed Brown for (k=0; k<s; k++) { 76861692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 76961692a83SJed Brown } 77061692a83SJed Brown } 77161692a83SJed Brown t->bt[i] = 0; 77261692a83SJed Brown for (j=0; j<s; j++) { 77361692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 77461692a83SJed Brown } 775fe7e6d57SJed Brown if (bembed) { 776fe7e6d57SJed Brown t->bembedt[i] = 0; 777fe7e6d57SJed Brown for (j=0; j<s; j++) { 778fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 779fe7e6d57SJed Brown } 780fe7e6d57SJed Brown } 78161692a83SJed Brown } 7828d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 7838d59e960SJed Brown 784f4aed992SEmil Constantinescu t->pinterp = pinterp; 785785e854fSJed Brown ierr = PetscMalloc1(s*pinterp,&t->binterpt);CHKERRQ(ierr); 786580bdb30SBarry Smith ierr = PetscArraycpy(t->binterpt,binterpt,s*pinterp);CHKERRQ(ierr); 78761692a83SJed Brown link->next = RosWTableauList; 78861692a83SJed Brown RosWTableauList = link; 789e27a552bSJed Brown PetscFunctionReturn(0); 790e27a552bSJed Brown } 791e27a552bSJed Brown 79242faf41dSJed Brown /*@C 793fd292e60Sprj- TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices 79442faf41dSJed Brown 79542faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 79642faf41dSJed Brown 79742faf41dSJed Brown Input Parameters: 79842faf41dSJed Brown + name - identifier for method 79942faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 80042faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 80142faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 80242faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 80342faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 804a2b725a8SWilliam Gropp - e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 80542faf41dSJed Brown 80642faf41dSJed Brown Notes: 80742faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 80842faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 80942faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 81042faf41dSJed Brown 81142faf41dSJed Brown Level: developer 81242faf41dSJed Brown 81342faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 81442faf41dSJed Brown @*/ 81519fd82e9SBarry Smith PetscErrorCode TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 81642faf41dSJed Brown { 81742faf41dSJed Brown PetscErrorCode ierr; 81842faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 81942faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 82042faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 82142faf41dSJed Brown p42 = one/eight - gamma/three, 82242faf41dSJed Brown p43 = one/twelve - gamma/three, 82342faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 82442faf41dSJed Brown p56 = one/twenty - gamma/four; 82542faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 82642faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 82742faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 82842faf41dSJed Brown 82942faf41dSJed Brown PetscFunctionBegin; 83042faf41dSJed Brown /* Step 1: choose Gamma (input) */ 83142faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 83242faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 83342faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 83442faf41dSJed Brown 83542faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 83642faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 83742faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 83842faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 83942faf41dSJed Brown rhs[0] = one - b3; 84042faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 84142faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 8426baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 84342faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 84442faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 84542faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 84642faf41dSJed Brown 84742faf41dSJed Brown /* Step 3 */ 84842faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 84942faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 85042faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 85142faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 85242faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 85342faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 85442faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 8556baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 85642faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 85742faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 85842faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 85942faf41dSJed Brown 86042faf41dSJed Brown /* Step 4: back-substitute */ 86142faf41dSJed Brown beta32 = beta32beta2p / beta2p; 86242faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 86342faf41dSJed Brown 86442faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 86542faf41dSJed Brown a43 = 0; 86642faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 86742faf41dSJed Brown a42 = a32; 86842faf41dSJed Brown 86942faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 87042faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 87142faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 87242faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 87342faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 87442faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 87542faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 87642faf41dSJed Brown Gamma[3][0] = beta4p-beta42-beta43-A[3][0]; Gamma[3][1] = beta42-A[3][1]; Gamma[3][2] = beta43-A[3][2]; Gamma[3][3] = gamma; 87742faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 87842faf41dSJed Brown 87942faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 88042faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 88142faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 88242faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 88342faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 88442faf41dSJed Brown 88542faf41dSJed Brown { 88642faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 88742faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 88842faf41dSJed Brown } 8890298fd71SBarry Smith ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,NULL);CHKERRQ(ierr); 89042faf41dSJed Brown PetscFunctionReturn(0); 89142faf41dSJed Brown } 89242faf41dSJed Brown 8931c3436cfSJed Brown /* 8941c3436cfSJed Brown The step completion formula is 8951c3436cfSJed Brown 8961c3436cfSJed Brown x1 = x0 + b^T Y 8971c3436cfSJed Brown 8981c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 8991c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 9001c3436cfSJed Brown 9011c3436cfSJed Brown x1e = x0 + be^T Y 9021c3436cfSJed Brown = x1 - b^T Y + be^T Y 9031c3436cfSJed Brown = x1 + (be - b)^T Y 9041c3436cfSJed Brown 9051c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9061c3436cfSJed Brown */ 907f9c1d6abSBarry Smith static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done) 9081c3436cfSJed Brown { 9091c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9101c3436cfSJed Brown RosWTableau tab = ros->tableau; 9111c3436cfSJed Brown PetscScalar *w = ros->work; 9121c3436cfSJed Brown PetscInt i; 9131c3436cfSJed Brown PetscErrorCode ierr; 9141c3436cfSJed Brown 9151c3436cfSJed Brown PetscFunctionBegin; 9161c3436cfSJed Brown if (order == tab->order) { 917108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 918f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 919de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 920f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 921f9c1d6abSBarry Smith } else {ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr);} 9221c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9231c3436cfSJed Brown PetscFunctionReturn(0); 9241c3436cfSJed Brown } else if (order == tab->order-1) { 9251c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 926108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 927f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 928de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 929f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 930108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 931108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 932f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 933f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 9341c3436cfSJed Brown } 9351c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9361c3436cfSJed Brown PetscFunctionReturn(0); 9371c3436cfSJed Brown } 9381c3436cfSJed Brown unavailable: 9391c3436cfSJed Brown if (done) *done = PETSC_FALSE; 940a7fac7c2SEmil Constantinescu else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order); 9411c3436cfSJed Brown PetscFunctionReturn(0); 9421c3436cfSJed Brown } 9431c3436cfSJed Brown 944560360afSLisandro Dalcin static PetscErrorCode TSRollBack_RosW(TS ts) 94524655328SShri { 94624655328SShri TS_RosW *ros = (TS_RosW*)ts->data; 94724655328SShri PetscErrorCode ierr; 94824655328SShri 94924655328SShri PetscFunctionBegin; 950be5899b3SLisandro Dalcin ierr = VecCopy(ros->vec_sol_prev,ts->vec_sol);CHKERRQ(ierr); 95124655328SShri PetscFunctionReturn(0); 95224655328SShri } 95324655328SShri 954e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 955e27a552bSJed Brown { 95661692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 95761692a83SJed Brown RosWTableau tab = ros->tableau; 958e27a552bSJed Brown const PetscInt s = tab->s; 9591c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9600feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 961c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 96261692a83SJed Brown PetscScalar *w = ros->work; 9637d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 964e27a552bSJed Brown SNES snes; 9651c3436cfSJed Brown TSAdapt adapt; 966fecfb714SLisandro Dalcin PetscInt i,j,its,lits; 967be5899b3SLisandro Dalcin PetscInt rejections = 0; 968b39943a6SLisandro Dalcin PetscBool stageok,accept = PETSC_TRUE; 969b39943a6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 970e27a552bSJed Brown PetscErrorCode ierr; 971f7f07198SBarry Smith PetscInt lag; 972e27a552bSJed Brown 973e27a552bSJed Brown PetscFunctionBegin; 974b39943a6SLisandro Dalcin if (!ts->steprollback) { 975be5899b3SLisandro Dalcin ierr = VecCopy(ts->vec_sol,ros->vec_sol_prev);CHKERRQ(ierr); 976b39943a6SLisandro Dalcin } 977e27a552bSJed Brown 978b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 979b39943a6SLisandro Dalcin while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 9801c3436cfSJed Brown const PetscReal h = ts->time_step; 981e27a552bSJed Brown for (i=0; i<s; i++) { 9821c3436cfSJed Brown ros->stage_time = ts->ptime + h*ASum[i]; 983b8123daeSJed Brown ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 984c17803e7SJed Brown if (GammaZeroDiag[i]) { 985c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 986b296d7d5SJed Brown ros->scoeff = 1.; 987c17803e7SJed Brown } else { 988c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 989b296d7d5SJed Brown ros->scoeff = 1./Gamma[i*s+i]; 990fd96d5b0SEmil Constantinescu } 99161692a83SJed Brown 99261692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 993de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 994de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 99561692a83SJed Brown 99661692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 99761692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 99861692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 99961692a83SJed Brown 1000e27a552bSJed Brown /* Initial guess taken from last stage */ 100161692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 100261692a83SJed Brown 10037d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 1004b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 100561692a83SJed Brown if (!ros->recompute_jacobian && !i) { 1006f7f07198SBarry Smith ierr = SNESGetLagJacobian(snes,&lag);CHKERRQ(ierr); 1007f7f07198SBarry Smith if (lag == 1) { /* use did not set a nontrival lag, so lag over all stages */ 1008f7f07198SBarry Smith ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */ 1009f7f07198SBarry Smith } 101061692a83SJed Brown } 10110298fd71SBarry Smith ierr = SNESSolve(snes,NULL,Y[i]);CHKERRQ(ierr); 1012f7f07198SBarry Smith if (!ros->recompute_jacobian && i == s-1 && lag == 1) { 1013f7f07198SBarry Smith ierr = SNESSetLagJacobian(snes,lag);CHKERRQ(ierr); /* Set lag back to 1 so we know user did not set it */ 1014f7f07198SBarry Smith } 1015e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1016e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10175ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 10187d4bf2deSEmil Constantinescu } else { 10191ce71dffSSatish Balay Mat J,Jp; 10200feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10210feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 102222d28d08SBarry Smith ierr = VecScale(Y[i],-1.0);CHKERRQ(ierr); 10230feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 10240feba352SEmil Constantinescu 10250feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10260feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10270feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 1028fecfb714SLisandro Dalcin 1029fecfb714SLisandro Dalcin /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10300298fd71SBarry Smith ierr = TSGetIJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 1031d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,J,Jp,PETSC_FALSE);CHKERRQ(ierr); 103222d28d08SBarry Smith ierr = MatMult(J,Zstage,Zdot);CHKERRQ(ierr); 10330feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10345ef26d82SJed Brown ts->ksp_its += 1; 1035fecfb714SLisandro Dalcin 1036fecfb714SLisandro Dalcin ierr = VecScale(Y[i],h);CHKERRQ(ierr); 10377d4bf2deSEmil Constantinescu } 10389be3e283SDebojyoti Ghosh ierr = TSPostStage(ts,ros->stage_time,i,Y);CHKERRQ(ierr); 1039fecfb714SLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1040fecfb714SLisandro Dalcin ierr = TSAdaptCheckStage(adapt,ts,ros->stage_time,Y[i],&stageok);CHKERRQ(ierr); 1041fecfb714SLisandro Dalcin if (!stageok) goto reject_step; 1042e27a552bSJed Brown } 1043e27a552bSJed Brown 1044b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 1045b39943a6SLisandro Dalcin ierr = TSEvaluateStep_RosW(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 1046b39943a6SLisandro Dalcin ros->status = TS_STEP_PENDING; 1047552698daSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 10481c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10491917a363SLisandro Dalcin ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 1050fecfb714SLisandro Dalcin ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 1051b39943a6SLisandro Dalcin ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1052b39943a6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 1053b39943a6SLisandro Dalcin ierr = TSRollBack_RosW(ts);CHKERRQ(ierr); 1054be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1055be5899b3SLisandro Dalcin goto reject_step; 1056b39943a6SLisandro Dalcin } 1057b39943a6SLisandro Dalcin 1058e27a552bSJed Brown ts->ptime += ts->time_step; 1059cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 10601c3436cfSJed Brown break; 1061b39943a6SLisandro Dalcin 1062b39943a6SLisandro Dalcin reject_step: 1063fecfb714SLisandro Dalcin ts->reject++; accept = PETSC_FALSE; 1064be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1065b39943a6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 1066be5899b3SLisandro Dalcin ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 10671c3436cfSJed Brown } 10681c3436cfSJed Brown } 1069e27a552bSJed Brown PetscFunctionReturn(0); 1070e27a552bSJed Brown } 1071e27a552bSJed Brown 1072f9c1d6abSBarry Smith static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1073e27a552bSJed Brown { 107461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1075f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1076f4aed992SEmil Constantinescu PetscReal h; 1077f4aed992SEmil Constantinescu PetscReal tt,t; 1078f4aed992SEmil Constantinescu PetscScalar *bt; 1079f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1080f4aed992SEmil Constantinescu PetscErrorCode ierr; 1081f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1082f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1083f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1084e27a552bSJed Brown 1085e27a552bSJed Brown PetscFunctionBegin; 1086ce94432eSBarry Smith if (!Bt) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1087f4aed992SEmil Constantinescu 1088f4aed992SEmil Constantinescu switch (ros->status) { 1089f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1090f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1091f4aed992SEmil Constantinescu h = ts->time_step; 1092f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1093f4aed992SEmil Constantinescu break; 1094f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1095be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1096f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1097f4aed992SEmil Constantinescu break; 1098ce94432eSBarry Smith default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1099f4aed992SEmil Constantinescu } 1100785e854fSJed Brown ierr = PetscMalloc1(s,&bt);CHKERRQ(ierr); 1101f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1102f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1103f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 11043ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1105f4aed992SEmil Constantinescu } 1106f4aed992SEmil Constantinescu } 1107f4aed992SEmil Constantinescu 1108f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1109f9c1d6abSBarry Smith /* U <- 0*/ 1110f9c1d6abSBarry Smith ierr = VecZeroEntries(U);CHKERRQ(ierr); 1111f9c1d6abSBarry Smith /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11123ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j] = 0; 11133ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11143ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11153ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1116f4aed992SEmil Constantinescu } 11173ca35412SEmil Constantinescu } 1118f9c1d6abSBarry Smith ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1119be5899b3SLisandro Dalcin /* U <- y(t) + U */ 1120be5899b3SLisandro Dalcin ierr = VecAXPY(U,1,ros->vec_sol_prev);CHKERRQ(ierr); 1121f4aed992SEmil Constantinescu 1122f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1123e27a552bSJed Brown PetscFunctionReturn(0); 1124e27a552bSJed Brown } 1125e27a552bSJed Brown 1126e27a552bSJed Brown /*------------------------------------------------------------*/ 1127b39943a6SLisandro Dalcin 1128b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauReset(TS ts) 1129b39943a6SLisandro Dalcin { 1130b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1131b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1132b39943a6SLisandro Dalcin PetscErrorCode ierr; 1133b39943a6SLisandro Dalcin 1134b39943a6SLisandro Dalcin PetscFunctionBegin; 1135b39943a6SLisandro Dalcin if (!tab) PetscFunctionReturn(0); 1136b39943a6SLisandro Dalcin ierr = VecDestroyVecs(tab->s,&ros->Y);CHKERRQ(ierr); 1137b39943a6SLisandro Dalcin ierr = PetscFree(ros->work);CHKERRQ(ierr); 1138b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1139b39943a6SLisandro Dalcin } 1140b39943a6SLisandro Dalcin 1141e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1142e27a552bSJed Brown { 114361692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1144e27a552bSJed Brown PetscErrorCode ierr; 1145e27a552bSJed Brown 1146e27a552bSJed Brown PetscFunctionBegin; 1147b39943a6SLisandro Dalcin ierr = TSRosWTableauReset(ts);CHKERRQ(ierr); 114861692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 114961692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 115061692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 115161692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 1152be5899b3SLisandro Dalcin ierr = VecDestroy(&ros->vec_sol_prev);CHKERRQ(ierr); 1153e27a552bSJed Brown PetscFunctionReturn(0); 1154e27a552bSJed Brown } 1155e27a552bSJed Brown 1156d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1157d5e6173cSPeter Brune { 1158d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1159d5e6173cSPeter Brune PetscErrorCode ierr; 1160d5e6173cSPeter Brune 1161d5e6173cSPeter Brune PetscFunctionBegin; 1162d5e6173cSPeter Brune if (Ydot) { 1163d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1164d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1165d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1166d5e6173cSPeter Brune } 1167d5e6173cSPeter Brune if (Zdot) { 1168d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1169d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1170d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1171d5e6173cSPeter Brune } 1172d5e6173cSPeter Brune if (Ystage) { 1173d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1174d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1175d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1176d5e6173cSPeter Brune } 1177d5e6173cSPeter Brune if (Zstage) { 1178d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1179d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1180d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1181d5e6173cSPeter Brune } 1182d5e6173cSPeter Brune PetscFunctionReturn(0); 1183d5e6173cSPeter Brune } 1184d5e6173cSPeter Brune 1185d5e6173cSPeter Brune 1186d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1187d5e6173cSPeter Brune { 1188d5e6173cSPeter Brune PetscErrorCode ierr; 1189d5e6173cSPeter Brune 1190d5e6173cSPeter Brune PetscFunctionBegin; 1191d5e6173cSPeter Brune if (Ydot) { 1192d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1193d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1194d5e6173cSPeter Brune } 1195d5e6173cSPeter Brune } 1196d5e6173cSPeter Brune if (Zdot) { 1197d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1198d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1199d5e6173cSPeter Brune } 1200d5e6173cSPeter Brune } 1201d5e6173cSPeter Brune if (Ystage) { 1202d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1203d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1204d5e6173cSPeter Brune } 1205d5e6173cSPeter Brune } 1206d5e6173cSPeter Brune if (Zstage) { 1207d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1208d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1209d5e6173cSPeter Brune } 1210d5e6173cSPeter Brune } 1211d5e6173cSPeter Brune PetscFunctionReturn(0); 1212d5e6173cSPeter Brune } 1213d5e6173cSPeter Brune 1214d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1215d5e6173cSPeter Brune { 1216d5e6173cSPeter Brune PetscFunctionBegin; 1217d5e6173cSPeter Brune PetscFunctionReturn(0); 1218d5e6173cSPeter Brune } 1219d5e6173cSPeter Brune 1220d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1221d5e6173cSPeter Brune { 1222d5e6173cSPeter Brune TS ts = (TS)ctx; 1223d5e6173cSPeter Brune PetscErrorCode ierr; 1224d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1225d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1226d5e6173cSPeter Brune 1227d5e6173cSPeter Brune PetscFunctionBegin; 1228d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1229d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1230d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1231d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1232d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1233d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1234d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1235d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1236d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1237d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1238d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1239d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1240d5e6173cSPeter Brune PetscFunctionReturn(0); 1241d5e6173cSPeter Brune } 1242d5e6173cSPeter Brune 1243258e1594SPeter Brune 1244258e1594SPeter Brune static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx) 1245258e1594SPeter Brune { 1246258e1594SPeter Brune PetscFunctionBegin; 1247258e1594SPeter Brune PetscFunctionReturn(0); 1248258e1594SPeter Brune } 1249258e1594SPeter Brune 1250258e1594SPeter Brune static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1251258e1594SPeter Brune { 1252258e1594SPeter Brune TS ts = (TS)ctx; 1253258e1594SPeter Brune PetscErrorCode ierr; 1254258e1594SPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1255258e1594SPeter Brune Vec Ydots,Zdots,Ystages,Zstages; 1256258e1594SPeter Brune 1257258e1594SPeter Brune PetscFunctionBegin; 1258258e1594SPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1259258e1594SPeter Brune ierr = TSRosWGetVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1260258e1594SPeter Brune 1261258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1262258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1263258e1594SPeter Brune 1264258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1265258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1266258e1594SPeter Brune 1267258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1268258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1269258e1594SPeter Brune 1270258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1271258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1272258e1594SPeter Brune 1273258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1274258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1275258e1594SPeter Brune PetscFunctionReturn(0); 1276258e1594SPeter Brune } 1277258e1594SPeter Brune 1278e27a552bSJed Brown /* 1279e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1280e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1281e27a552bSJed Brown */ 1282f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1283e27a552bSJed Brown { 128461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1285e27a552bSJed Brown PetscErrorCode ierr; 1286d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1287b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1288d5e6173cSPeter Brune DM dm,dmsave; 1289e27a552bSJed Brown 1290e27a552bSJed Brown PetscFunctionBegin; 1291d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1292d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1293b296d7d5SJed Brown ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1294f9c1d6abSBarry Smith ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1295d5e6173cSPeter Brune dmsave = ts->dm; 1296d5e6173cSPeter Brune ts->dm = dm; 1297d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1298d5e6173cSPeter Brune ts->dm = dmsave; 1299d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1300e27a552bSJed Brown PetscFunctionReturn(0); 1301e27a552bSJed Brown } 1302e27a552bSJed Brown 1303d1e9a80fSBarry Smith static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat A,Mat B,TS ts) 1304e27a552bSJed Brown { 130561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1306d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1307b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1308e27a552bSJed Brown PetscErrorCode ierr; 1309d5e6173cSPeter Brune DM dm,dmsave; 1310e27a552bSJed Brown 1311e27a552bSJed Brown PetscFunctionBegin; 131261692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1313d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1314d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1315d5e6173cSPeter Brune dmsave = ts->dm; 1316d5e6173cSPeter Brune ts->dm = dm; 1317d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,PETSC_TRUE);CHKERRQ(ierr); 1318d5e6173cSPeter Brune ts->dm = dmsave; 1319d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1320e27a552bSJed Brown PetscFunctionReturn(0); 1321e27a552bSJed Brown } 1322e27a552bSJed Brown 1323b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauSetUp(TS ts) 1324b39943a6SLisandro Dalcin { 1325b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1326b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1327b39943a6SLisandro Dalcin PetscErrorCode ierr; 1328b39943a6SLisandro Dalcin 1329b39943a6SLisandro Dalcin PetscFunctionBegin; 1330b39943a6SLisandro Dalcin ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ros->Y);CHKERRQ(ierr); 1331b39943a6SLisandro Dalcin ierr = PetscMalloc1(tab->s,&ros->work);CHKERRQ(ierr); 1332b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1333b39943a6SLisandro Dalcin } 1334b39943a6SLisandro Dalcin 1335e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1336e27a552bSJed Brown { 133761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1338e27a552bSJed Brown PetscErrorCode ierr; 1339d5e6173cSPeter Brune DM dm; 1340b39943a6SLisandro Dalcin SNES snes; 1341e27a552bSJed Brown 1342e27a552bSJed Brown PetscFunctionBegin; 1343b39943a6SLisandro Dalcin ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr); 134461692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 134561692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 134661692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 134761692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 1348be5899b3SLisandro Dalcin ierr = VecDuplicate(ts->vec_sol,&ros->vec_sol_prev);CHKERRQ(ierr); 134922d28d08SBarry Smith ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1350d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1351258e1594SPeter Brune ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1352b39943a6SLisandro Dalcin /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1353b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1354b39943a6SLisandro Dalcin if (!((PetscObject)snes)->type_name) { 1355b39943a6SLisandro Dalcin ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1356b39943a6SLisandro Dalcin } 1357e27a552bSJed Brown PetscFunctionReturn(0); 1358e27a552bSJed Brown } 1359e27a552bSJed Brown /*------------------------------------------------------------*/ 1360e27a552bSJed Brown 13614416b707SBarry Smith static PetscErrorCode TSSetFromOptions_RosW(PetscOptionItems *PetscOptionsObject,TS ts) 1362e27a552bSJed Brown { 136361692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1364e27a552bSJed Brown PetscErrorCode ierr; 1365b39943a6SLisandro Dalcin SNES snes; 1366e27a552bSJed Brown 1367e27a552bSJed Brown PetscFunctionBegin; 1368e55864a3SBarry Smith ierr = PetscOptionsHead(PetscOptionsObject,"RosW ODE solver options");CHKERRQ(ierr); 1369e27a552bSJed Brown { 137061692a83SJed Brown RosWTableauLink link; 1371e27a552bSJed Brown PetscInt count,choice; 1372e27a552bSJed Brown PetscBool flg; 1373e27a552bSJed Brown const char **namelist; 137461692a83SJed Brown 137561692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1376f489ac74SBarry Smith ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 137761692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1378b39943a6SLisandro Dalcin ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,ros->tableau->name,&choice,&flg);CHKERRQ(ierr); 1379b39943a6SLisandro Dalcin if (flg) {ierr = TSRosWSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1380e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 138161692a83SJed Brown 13820298fd71SBarry Smith ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,NULL);CHKERRQ(ierr); 1383b39943a6SLisandro Dalcin } 1384b39943a6SLisandro Dalcin ierr = PetscOptionsTail();CHKERRQ(ierr); 138561692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 138661692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 138761692a83SJed Brown if (!((PetscObject)snes)->type_name) { 138861692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 138961692a83SJed Brown } 1390e27a552bSJed Brown PetscFunctionReturn(0); 1391e27a552bSJed Brown } 1392e27a552bSJed Brown 1393e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1394e27a552bSJed Brown { 139561692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1396e27a552bSJed Brown PetscBool iascii; 1397e27a552bSJed Brown PetscErrorCode ierr; 1398e27a552bSJed Brown 1399e27a552bSJed Brown PetscFunctionBegin; 1400251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1401e27a552bSJed Brown if (iascii) { 14029c334d8fSLisandro Dalcin RosWTableau tab = ros->tableau; 140319fd82e9SBarry Smith TSRosWType rostype; 14049c334d8fSLisandro Dalcin char buf[512]; 1405e408995aSJed Brown PetscInt i; 1406e408995aSJed Brown PetscReal abscissa[512]; 140761692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 140861692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1409de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 141061692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1411e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1412de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1413e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1414e27a552bSJed Brown } 1415e27a552bSJed Brown PetscFunctionReturn(0); 1416e27a552bSJed Brown } 1417e27a552bSJed Brown 14189200755eSBarry Smith static PetscErrorCode TSLoad_RosW(TS ts,PetscViewer viewer) 14199200755eSBarry Smith { 14209200755eSBarry Smith PetscErrorCode ierr; 14219200755eSBarry Smith SNES snes; 14229c334d8fSLisandro Dalcin TSAdapt adapt; 14239200755eSBarry Smith 14249200755eSBarry Smith PetscFunctionBegin; 14259c334d8fSLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 14269c334d8fSLisandro Dalcin ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 14279200755eSBarry Smith ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 14289200755eSBarry Smith ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 14299200755eSBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 14309200755eSBarry Smith ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 14319200755eSBarry Smith ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 14329200755eSBarry Smith PetscFunctionReturn(0); 14339200755eSBarry Smith } 14349200755eSBarry Smith 1435e27a552bSJed Brown /*@C 143661692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1437e27a552bSJed Brown 1438e27a552bSJed Brown Logically collective 1439e27a552bSJed Brown 1440e27a552bSJed Brown Input Parameter: 1441e27a552bSJed Brown + ts - timestepping context 1442b92453a8SLisandro Dalcin - roswtype - type of Rosenbrock-W scheme 1443e27a552bSJed Brown 1444020d8f30SJed Brown Level: beginner 1445e27a552bSJed Brown 1446020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1447e27a552bSJed Brown @*/ 1448b92453a8SLisandro Dalcin PetscErrorCode TSRosWSetType(TS ts,TSRosWType roswtype) 1449e27a552bSJed Brown { 1450e27a552bSJed Brown PetscErrorCode ierr; 1451e27a552bSJed Brown 1452e27a552bSJed Brown PetscFunctionBegin; 1453e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1454b92453a8SLisandro Dalcin PetscValidCharPointer(roswtype,2); 1455b92453a8SLisandro Dalcin ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,roswtype));CHKERRQ(ierr); 1456e27a552bSJed Brown PetscFunctionReturn(0); 1457e27a552bSJed Brown } 1458e27a552bSJed Brown 1459e27a552bSJed Brown /*@C 146061692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1461e27a552bSJed Brown 1462e27a552bSJed Brown Logically collective 1463e27a552bSJed Brown 1464e27a552bSJed Brown Input Parameter: 1465e27a552bSJed Brown . ts - timestepping context 1466e27a552bSJed Brown 1467e27a552bSJed Brown Output Parameter: 146861692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1469e27a552bSJed Brown 1470e27a552bSJed Brown Level: intermediate 1471e27a552bSJed Brown 1472e27a552bSJed Brown .seealso: TSRosWGetType() 1473e27a552bSJed Brown @*/ 147419fd82e9SBarry Smith PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1475e27a552bSJed Brown { 1476e27a552bSJed Brown PetscErrorCode ierr; 1477e27a552bSJed Brown 1478e27a552bSJed Brown PetscFunctionBegin; 1479e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 148019fd82e9SBarry Smith ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1481e27a552bSJed Brown PetscFunctionReturn(0); 1482e27a552bSJed Brown } 1483e27a552bSJed Brown 1484e27a552bSJed Brown /*@C 148561692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1486e27a552bSJed Brown 1487e27a552bSJed Brown Logically collective 1488e27a552bSJed Brown 1489e27a552bSJed Brown Input Parameter: 1490e27a552bSJed Brown + ts - timestepping context 149161692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1492e27a552bSJed Brown 1493e27a552bSJed Brown Level: intermediate 1494e27a552bSJed Brown 1495e27a552bSJed Brown .seealso: TSRosWGetType() 1496e27a552bSJed Brown @*/ 149761692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1498e27a552bSJed Brown { 1499e27a552bSJed Brown PetscErrorCode ierr; 1500e27a552bSJed Brown 1501e27a552bSJed Brown PetscFunctionBegin; 1502e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 150361692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1504e27a552bSJed Brown PetscFunctionReturn(0); 1505e27a552bSJed Brown } 1506e27a552bSJed Brown 1507560360afSLisandro Dalcin static PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1508e27a552bSJed Brown { 150961692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1510e27a552bSJed Brown 1511e27a552bSJed Brown PetscFunctionBegin; 151261692a83SJed Brown *rostype = ros->tableau->name; 1513e27a552bSJed Brown PetscFunctionReturn(0); 1514e27a552bSJed Brown } 1515ef20d060SBarry Smith 1516560360afSLisandro Dalcin static PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1517e27a552bSJed Brown { 151861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1519e27a552bSJed Brown PetscErrorCode ierr; 1520e27a552bSJed Brown PetscBool match; 152161692a83SJed Brown RosWTableauLink link; 1522e27a552bSJed Brown 1523e27a552bSJed Brown PetscFunctionBegin; 152461692a83SJed Brown if (ros->tableau) { 152561692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1526e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1527e27a552bSJed Brown } 152861692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 152961692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1530e27a552bSJed Brown if (match) { 1531b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauReset(ts);CHKERRQ(ierr);} 153261692a83SJed Brown ros->tableau = &link->tab; 1533b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr);} 15342ffb9264SLisandro Dalcin ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1535e27a552bSJed Brown PetscFunctionReturn(0); 1536e27a552bSJed Brown } 1537e27a552bSJed Brown } 1538ce94432eSBarry Smith SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1539e27a552bSJed Brown PetscFunctionReturn(0); 1540e27a552bSJed Brown } 154161692a83SJed Brown 1542560360afSLisandro Dalcin static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1543e27a552bSJed Brown { 154461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1545e27a552bSJed Brown 1546e27a552bSJed Brown PetscFunctionBegin; 154761692a83SJed Brown ros->recompute_jacobian = flg; 1548e27a552bSJed Brown PetscFunctionReturn(0); 1549e27a552bSJed Brown } 1550e27a552bSJed Brown 1551b3a6b972SJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1552b3a6b972SJed Brown { 1553b3a6b972SJed Brown PetscErrorCode ierr; 1554b3a6b972SJed Brown 1555b3a6b972SJed Brown PetscFunctionBegin; 1556b3a6b972SJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1557b3a6b972SJed Brown if (ts->dm) { 1558b3a6b972SJed Brown ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1559b3a6b972SJed Brown ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1560b3a6b972SJed Brown } 1561b3a6b972SJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 1562b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",NULL);CHKERRQ(ierr); 1563b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",NULL);CHKERRQ(ierr); 1564b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",NULL);CHKERRQ(ierr); 1565b3a6b972SJed Brown PetscFunctionReturn(0); 1566b3a6b972SJed Brown } 1567d5e6173cSPeter Brune 1568e27a552bSJed Brown /* ------------------------------------------------------------ */ 1569e27a552bSJed Brown /*MC 1570020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1571e27a552bSJed Brown 1572e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1573e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1574e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1575e27a552bSJed Brown 1576e27a552bSJed Brown Notes: 157761692a83SJed Brown This method currently only works with autonomous ODE and DAE. 157861692a83SJed Brown 1579d0685a90SJed Brown Consider trying TSARKIMEX if the stiff part is strongly nonlinear. 1580d0685a90SJed Brown 1581*3d5a8a6aSBarry Smith Since this uses a single linear solve per time-step if you wish to lag the jacobian or preconditioner computation you must use also -snes_lag_jacobian_persists true or -snes_lag_jacobian_preconditioner true 1582*3d5a8a6aSBarry Smith 1583e94cfbe0SPatrick Sanan Developer Notes: 158461692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 158561692a83SJed Brown 1586f9c1d6abSBarry Smith $ udot = f(u) 158761692a83SJed Brown 158861692a83SJed Brown by the stage equations 158961692a83SJed Brown 1590f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 159161692a83SJed Brown 159261692a83SJed Brown and step completion formula 159361692a83SJed Brown 1594f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 159561692a83SJed Brown 1596f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 159761692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 159861692a83SJed Brown we define new variables for the stage equations 159961692a83SJed Brown 160061692a83SJed Brown $ y_i = gamma_ij k_j 160161692a83SJed Brown 160261692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 160361692a83SJed Brown 1604b70472e2SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1} 160561692a83SJed Brown 160661692a83SJed Brown to rewrite the method as 160761692a83SJed Brown 1608f9c1d6abSBarry 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 1609f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j bt_j y_j 161061692a83SJed Brown 161161692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 161261692a83SJed Brown 161361692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 161461692a83SJed Brown 161561692a83SJed Brown or, more compactly in tensor notation 161661692a83SJed Brown 161761692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 161861692a83SJed Brown 161961692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 162061692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 162161692a83SJed Brown equation 162261692a83SJed Brown 1623f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 162461692a83SJed Brown 162561692a83SJed Brown with initial guess y_i = 0. 1626e27a552bSJed Brown 1627e27a552bSJed Brown Level: beginner 1628e27a552bSJed Brown 1629d0685a90SJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSWTHETA1, TSROSWTHETA2, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1630a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1631e27a552bSJed Brown M*/ 16328cc058d9SJed Brown PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1633e27a552bSJed Brown { 163461692a83SJed Brown TS_RosW *ros; 1635e27a552bSJed Brown PetscErrorCode ierr; 1636e27a552bSJed Brown 1637e27a552bSJed Brown PetscFunctionBegin; 1638607a6623SBarry Smith ierr = TSRosWInitializePackage();CHKERRQ(ierr); 1639e27a552bSJed Brown 1640e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1641e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1642e27a552bSJed Brown ts->ops->view = TSView_RosW; 16439200755eSBarry Smith ts->ops->load = TSLoad_RosW; 1644e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1645e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1646e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16471c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 164824655328SShri ts->ops->rollback = TSRollBack_RosW; 1649e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1650e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1651e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1652e27a552bSJed Brown 1653efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 1654efd4aadfSBarry Smith 1655b00a9115SJed Brown ierr = PetscNewLog(ts,&ros);CHKERRQ(ierr); 165661692a83SJed Brown ts->data = (void*)ros; 1657e27a552bSJed Brown 1658bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",TSRosWGetType_RosW);CHKERRQ(ierr); 1659bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",TSRosWSetType_RosW);CHKERRQ(ierr); 1660bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1661b39943a6SLisandro Dalcin 1662b39943a6SLisandro Dalcin ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1663e27a552bSJed Brown PetscFunctionReturn(0); 1664e27a552bSJed Brown } 1665