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 595e27a552bSJed Brown /*@C 596e27a552bSJed Brown TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 597e27a552bSJed Brown 598e27a552bSJed Brown Not Collective 599e27a552bSJed Brown 600e27a552bSJed Brown Level: advanced 601e27a552bSJed Brown 602607a6623SBarry Smith .seealso: TSRosWRegister(), TSRosWRegisterAll() 603e27a552bSJed Brown @*/ 604e27a552bSJed Brown PetscErrorCode TSRosWRegisterDestroy(void) 605e27a552bSJed Brown { 606e27a552bSJed Brown PetscErrorCode ierr; 60761692a83SJed Brown RosWTableauLink link; 608e27a552bSJed Brown 609e27a552bSJed Brown PetscFunctionBegin; 61061692a83SJed Brown while ((link = RosWTableauList)) { 61161692a83SJed Brown RosWTableau t = &link->tab; 61261692a83SJed Brown RosWTableauList = link->next; 61361692a83SJed Brown ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 61443b21953SEmil Constantinescu ierr = PetscFree5(t->At,t->bt,t->GammaInv,t->GammaZeroDiag,t->GammaExplicitCorr);CHKERRQ(ierr); 615fe7e6d57SJed Brown ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 616f4aed992SEmil Constantinescu ierr = PetscFree(t->binterpt);CHKERRQ(ierr); 617e27a552bSJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 618e27a552bSJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 619e27a552bSJed Brown } 620e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 621e27a552bSJed Brown PetscFunctionReturn(0); 622e27a552bSJed Brown } 623e27a552bSJed Brown 624e27a552bSJed Brown /*@C 625e27a552bSJed Brown TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 6268a690491SBarry Smith from TSInitializePackage(). 627e27a552bSJed Brown 628e27a552bSJed Brown Level: developer 629e27a552bSJed Brown 630e27a552bSJed Brown .seealso: PetscInitialize() 631e27a552bSJed Brown @*/ 632607a6623SBarry Smith PetscErrorCode TSRosWInitializePackage(void) 633e27a552bSJed Brown { 634e27a552bSJed Brown PetscErrorCode ierr; 635e27a552bSJed Brown 636e27a552bSJed Brown PetscFunctionBegin; 637e27a552bSJed Brown if (TSRosWPackageInitialized) PetscFunctionReturn(0); 638e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 639e27a552bSJed Brown ierr = TSRosWRegisterAll();CHKERRQ(ierr); 640e27a552bSJed Brown ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 641e27a552bSJed Brown PetscFunctionReturn(0); 642e27a552bSJed Brown } 643e27a552bSJed Brown 644e27a552bSJed Brown /*@C 645e27a552bSJed Brown TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 646e27a552bSJed Brown called from PetscFinalize(). 647e27a552bSJed Brown 648e27a552bSJed Brown Level: developer 649e27a552bSJed Brown 650e27a552bSJed Brown .seealso: PetscFinalize() 651e27a552bSJed Brown @*/ 652e27a552bSJed Brown PetscErrorCode TSRosWFinalizePackage(void) 653e27a552bSJed Brown { 654e27a552bSJed Brown PetscErrorCode ierr; 655e27a552bSJed Brown 656e27a552bSJed Brown PetscFunctionBegin; 657e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 658e27a552bSJed Brown ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 659e27a552bSJed Brown PetscFunctionReturn(0); 660e27a552bSJed Brown } 661e27a552bSJed Brown 662e27a552bSJed Brown /*@C 66361692a83SJed Brown TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 664e27a552bSJed Brown 665e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 666e27a552bSJed Brown 667e27a552bSJed Brown Input Parameters: 668e27a552bSJed Brown + name - identifier for method 669e27a552bSJed Brown . order - approximation order of method 670e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 67161692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 67261692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 673fe7e6d57SJed Brown . b - Step completion table (dimension s) 6740298fd71SBarry Smith . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available) 675f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 67642faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 677e27a552bSJed Brown 678e27a552bSJed Brown Notes: 67961692a83SJed Brown Several Rosenbrock W methods are provided, this function is only needed to create new methods. 680e27a552bSJed Brown 681e27a552bSJed Brown Level: advanced 682e27a552bSJed Brown 683e27a552bSJed Brown .seealso: TSRosW 684e27a552bSJed Brown @*/ 685f9c1d6abSBarry Smith PetscErrorCode TSRosWRegister(TSRosWType name,PetscInt order,PetscInt s,const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[], 686f4aed992SEmil Constantinescu PetscInt pinterp,const PetscReal binterpt[]) 687e27a552bSJed Brown { 688e27a552bSJed Brown PetscErrorCode ierr; 68961692a83SJed Brown RosWTableauLink link; 69061692a83SJed Brown RosWTableau t; 69161692a83SJed Brown PetscInt i,j,k; 69261692a83SJed Brown PetscScalar *GammaInv; 693e27a552bSJed Brown 694e27a552bSJed Brown PetscFunctionBegin; 695fe7e6d57SJed Brown PetscValidCharPointer(name,1); 696fe7e6d57SJed Brown PetscValidPointer(A,4); 697fe7e6d57SJed Brown PetscValidPointer(Gamma,5); 698fe7e6d57SJed Brown PetscValidPointer(b,6); 699fe7e6d57SJed Brown if (bembed) PetscValidPointer(bembed,7); 700fe7e6d57SJed Brown 7011d36bdfdSBarry Smith ierr = TSRosWInitializePackage();CHKERRQ(ierr); 702071fcb05SBarry Smith ierr = PetscNew(&link);CHKERRQ(ierr); 703e27a552bSJed Brown t = &link->tab; 704e27a552bSJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 705e27a552bSJed Brown t->order = order; 706e27a552bSJed Brown t->s = s; 707dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->A,s*s,&t->Gamma,s,&t->b,s,&t->ASum,s,&t->GammaSum);CHKERRQ(ierr); 708dcca6d9dSJed Brown ierr = PetscMalloc5(s*s,&t->At,s,&t->bt,s*s,&t->GammaInv,s,&t->GammaZeroDiag,s*s,&t->GammaExplicitCorr);CHKERRQ(ierr); 709580bdb30SBarry Smith ierr = PetscArraycpy(t->A,A,s*s);CHKERRQ(ierr); 710580bdb30SBarry Smith ierr = PetscArraycpy(t->Gamma,Gamma,s*s);CHKERRQ(ierr); 711580bdb30SBarry Smith ierr = PetscArraycpy(t->GammaExplicitCorr,Gamma,s*s);CHKERRQ(ierr); 712580bdb30SBarry Smith ierr = PetscArraycpy(t->b,b,s);CHKERRQ(ierr); 713fe7e6d57SJed Brown if (bembed) { 714dcca6d9dSJed Brown ierr = PetscMalloc2(s,&t->bembed,s,&t->bembedt);CHKERRQ(ierr); 715580bdb30SBarry Smith ierr = PetscArraycpy(t->bembed,bembed,s);CHKERRQ(ierr); 716fe7e6d57SJed Brown } 71761692a83SJed Brown for (i=0; i<s; i++) { 71861692a83SJed Brown t->ASum[i] = 0; 71961692a83SJed Brown t->GammaSum[i] = 0; 72061692a83SJed Brown for (j=0; j<s; j++) { 72161692a83SJed Brown t->ASum[i] += A[i*s+j]; 722fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i*s+j]; 72361692a83SJed Brown } 72461692a83SJed Brown } 725785e854fSJed Brown ierr = PetscMalloc1(s*s,&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 72661692a83SJed Brown for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 727fd96d5b0SEmil Constantinescu for (i=0; i<s; i++) { 728fd96d5b0SEmil Constantinescu if (Gamma[i*s+i] == 0.0) { 729fd96d5b0SEmil Constantinescu GammaInv[i*s+i] = 1.0; 730c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 731fd96d5b0SEmil Constantinescu } else { 732c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 733fd96d5b0SEmil Constantinescu } 734fd96d5b0SEmil Constantinescu } 735fd96d5b0SEmil Constantinescu 73661692a83SJed Brown switch (s) { 73761692a83SJed Brown case 1: GammaInv[0] = 1./GammaInv[0]; break; 7382e92ee13SHong Zhang case 2: ierr = PetscKernel_A_gets_inverse_A_2(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7396baedc03SHong Zhang case 3: ierr = PetscKernel_A_gets_inverse_A_3(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7402e92ee13SHong Zhang case 4: ierr = PetscKernel_A_gets_inverse_A_4(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 74161692a83SJed Brown case 5: { 74261692a83SJed Brown PetscInt ipvt5[5]; 74361692a83SJed Brown MatScalar work5[5*5]; 7442e92ee13SHong Zhang ierr = PetscKernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 74561692a83SJed Brown } 7462e92ee13SHong Zhang case 6: ierr = PetscKernel_A_gets_inverse_A_6(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 7472e92ee13SHong Zhang case 7: ierr = PetscKernel_A_gets_inverse_A_7(GammaInv,0,PETSC_FALSE,NULL);CHKERRQ(ierr); break; 74861692a83SJed Brown default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 74961692a83SJed Brown } 75061692a83SJed Brown for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 75161692a83SJed Brown ierr = PetscFree(GammaInv);CHKERRQ(ierr); 75243b21953SEmil Constantinescu 75343b21953SEmil Constantinescu for (i=0; i<s; i++) { 75443b21953SEmil Constantinescu for (k=0; k<i+1; k++) { 75543b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]=(t->GammaExplicitCorr[i*s+k])*(t->GammaInv[k*s+k]); 75643b21953SEmil Constantinescu for (j=k+1; j<i+1; j++) { 75743b21953SEmil Constantinescu t->GammaExplicitCorr[i*s+k]+=(t->GammaExplicitCorr[i*s+j])*(t->GammaInv[j*s+k]); 75843b21953SEmil Constantinescu } 75943b21953SEmil Constantinescu } 76043b21953SEmil Constantinescu } 76143b21953SEmil Constantinescu 76261692a83SJed Brown for (i=0; i<s; i++) { 76361692a83SJed Brown for (j=0; j<s; j++) { 76461692a83SJed Brown t->At[i*s+j] = 0; 76561692a83SJed Brown for (k=0; k<s; k++) { 76661692a83SJed Brown t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 76761692a83SJed Brown } 76861692a83SJed Brown } 76961692a83SJed Brown t->bt[i] = 0; 77061692a83SJed Brown for (j=0; j<s; j++) { 77161692a83SJed Brown t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 77261692a83SJed Brown } 773fe7e6d57SJed Brown if (bembed) { 774fe7e6d57SJed Brown t->bembedt[i] = 0; 775fe7e6d57SJed Brown for (j=0; j<s; j++) { 776fe7e6d57SJed Brown t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 777fe7e6d57SJed Brown } 778fe7e6d57SJed Brown } 77961692a83SJed Brown } 7808d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 7818d59e960SJed Brown 782f4aed992SEmil Constantinescu t->pinterp = pinterp; 783785e854fSJed Brown ierr = PetscMalloc1(s*pinterp,&t->binterpt);CHKERRQ(ierr); 784580bdb30SBarry Smith ierr = PetscArraycpy(t->binterpt,binterpt,s*pinterp);CHKERRQ(ierr); 78561692a83SJed Brown link->next = RosWTableauList; 78661692a83SJed Brown RosWTableauList = link; 787e27a552bSJed Brown PetscFunctionReturn(0); 788e27a552bSJed Brown } 789e27a552bSJed Brown 79042faf41dSJed Brown /*@C 791fd292e60Sprj- TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices 79242faf41dSJed Brown 79342faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 79442faf41dSJed Brown 79542faf41dSJed Brown Input Parameters: 79642faf41dSJed Brown + name - identifier for method 79742faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 79842faf41dSJed Brown . a2 - design parameter, see Table 7.2 of Hairer&Wanner 79942faf41dSJed Brown . a3 - design parameter or PETSC_DEFAULT to satisfy one of the order five conditions (Eq 7.22) 80042faf41dSJed Brown . b3 - design parameter, see Table 7.2 of Hairer&Wanner 80142faf41dSJed Brown . beta43 - design parameter or PETSC_DEFAULT to use Equation 7.21 of Hairer&Wanner 802a2b725a8SWilliam Gropp - e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 80342faf41dSJed Brown 80442faf41dSJed Brown Notes: 80542faf41dSJed Brown This routine encodes the design of fourth order Rosenbrock methods as described in Hairer and Wanner volume 2. 80642faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 80742faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 80842faf41dSJed Brown 80942faf41dSJed Brown Level: developer 81042faf41dSJed Brown 81142faf41dSJed Brown .seealso: TSRosW, TSRosWRegister() 81242faf41dSJed Brown @*/ 81319fd82e9SBarry Smith PetscErrorCode TSRosWRegisterRos4(TSRosWType name,PetscReal gamma,PetscReal a2,PetscReal a3,PetscReal b3,PetscReal e4) 81442faf41dSJed Brown { 81542faf41dSJed Brown PetscErrorCode ierr; 81642faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 81742faf41dSJed Brown const PetscReal one = 1,two = 2,three = 3,four = 4,five = 5,six = 6,eight = 8,twelve = 12,twenty = 20,twentyfour = 24, 81842faf41dSJed Brown p32 = one/six - gamma + gamma*gamma, 81942faf41dSJed Brown p42 = one/eight - gamma/three, 82042faf41dSJed Brown p43 = one/twelve - gamma/three, 82142faf41dSJed Brown p44 = one/twentyfour - gamma/two + three/two*gamma*gamma - gamma*gamma*gamma, 82242faf41dSJed Brown p56 = one/twenty - gamma/four; 82342faf41dSJed Brown PetscReal a4,a32,a42,a43,b1,b2,b4,beta2p,beta3p,beta4p,beta32,beta42,beta43,beta32beta2p,beta4jbetajp; 82442faf41dSJed Brown PetscReal A[4][4],Gamma[4][4],b[4],bm[4]; 82542faf41dSJed Brown PetscScalar M[3][3],rhs[3]; 82642faf41dSJed Brown 82742faf41dSJed Brown PetscFunctionBegin; 82842faf41dSJed Brown /* Step 1: choose Gamma (input) */ 82942faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 83042faf41dSJed Brown if (a3 == PETSC_DEFAULT) a3 = (one/five - a2/four)/(one/four - a2/three); /* Eq 7.22 */ 83142faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 83242faf41dSJed Brown 83342faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 83442faf41dSJed Brown M[0][0] = one; M[0][1] = one; M[0][2] = one; /* 7.15a */ 83542faf41dSJed Brown M[1][0] = 0.0; M[1][1] = a2*a2; M[1][2] = a4*a4; /* 7.15c */ 83642faf41dSJed Brown M[2][0] = 0.0; M[2][1] = a2*a2*a2; M[2][2] = a4*a4*a4; /* 7.15e */ 83742faf41dSJed Brown rhs[0] = one - b3; 83842faf41dSJed Brown rhs[1] = one/three - a3*a3*b3; 83942faf41dSJed Brown rhs[2] = one/four - a3*a3*a3*b3; 8406baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 84142faf41dSJed Brown b1 = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 84242faf41dSJed Brown b2 = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 84342faf41dSJed Brown b4 = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 84442faf41dSJed Brown 84542faf41dSJed Brown /* Step 3 */ 84642faf41dSJed Brown beta43 = (p56 - a2*p43) / (b4*a3*a3*(a3 - a2)); /* 7.21 */ 84742faf41dSJed Brown beta32beta2p = p44 / (b4*beta43); /* 7.15h */ 84842faf41dSJed Brown beta4jbetajp = (p32 - b3*beta32beta2p) / b4; 84942faf41dSJed Brown M[0][0] = b2; M[0][1] = b3; M[0][2] = b4; 85042faf41dSJed Brown M[1][0] = a4*a4*beta32beta2p-a3*a3*beta4jbetajp; M[1][1] = a2*a2*beta4jbetajp; M[1][2] = -a2*a2*beta32beta2p; 85142faf41dSJed Brown M[2][0] = b4*beta43*a3*a3-p43; M[2][1] = -b4*beta43*a2*a2; M[2][2] = 0; 85242faf41dSJed Brown rhs[0] = one/two - gamma; rhs[1] = 0; rhs[2] = -a2*a2*p32; 8536baedc03SHong Zhang ierr = PetscKernel_A_gets_inverse_A_3(&M[0][0],0,PETSC_FALSE,NULL);CHKERRQ(ierr); 85442faf41dSJed Brown beta2p = PetscRealPart(M[0][0]*rhs[0] + M[0][1]*rhs[1] + M[0][2]*rhs[2]); 85542faf41dSJed Brown beta3p = PetscRealPart(M[1][0]*rhs[0] + M[1][1]*rhs[1] + M[1][2]*rhs[2]); 85642faf41dSJed Brown beta4p = PetscRealPart(M[2][0]*rhs[0] + M[2][1]*rhs[1] + M[2][2]*rhs[2]); 85742faf41dSJed Brown 85842faf41dSJed Brown /* Step 4: back-substitute */ 85942faf41dSJed Brown beta32 = beta32beta2p / beta2p; 86042faf41dSJed Brown beta42 = (beta4jbetajp - beta43*beta3p) / beta2p; 86142faf41dSJed Brown 86242faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 86342faf41dSJed Brown a43 = 0; 86442faf41dSJed Brown a32 = p42 / (b3*a3*beta2p + b4*a4*beta2p); 86542faf41dSJed Brown a42 = a32; 86642faf41dSJed Brown 86742faf41dSJed Brown A[0][0] = 0; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0; 86842faf41dSJed Brown A[1][0] = a2; A[1][1] = 0; A[1][2] = 0; A[1][3] = 0; 86942faf41dSJed Brown A[2][0] = a3-a32; A[2][1] = a32; A[2][2] = 0; A[2][3] = 0; 87042faf41dSJed Brown A[3][0] = a4-a43-a42; A[3][1] = a42; A[3][2] = a43; A[3][3] = 0; 87142faf41dSJed Brown Gamma[0][0] = gamma; Gamma[0][1] = 0; Gamma[0][2] = 0; Gamma[0][3] = 0; 87242faf41dSJed Brown Gamma[1][0] = beta2p-A[1][0]; Gamma[1][1] = gamma; Gamma[1][2] = 0; Gamma[1][3] = 0; 87342faf41dSJed Brown Gamma[2][0] = beta3p-beta32-A[2][0]; Gamma[2][1] = beta32-A[2][1]; Gamma[2][2] = gamma; Gamma[2][3] = 0; 87442faf41dSJed 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; 87542faf41dSJed Brown b[0] = b1; b[1] = b2; b[2] = b3; b[3] = b4; 87642faf41dSJed Brown 87742faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 87842faf41dSJed Brown bm[3] = b[3] - e4*gamma; /* using definition of E4 */ 87942faf41dSJed Brown bm[2] = (p32 - beta4jbetajp*bm[3]) / (beta32*beta2p); /* fourth row of 7.18 */ 88042faf41dSJed Brown bm[1] = (one/two - gamma - beta3p*bm[2] - beta4p*bm[3]) / beta2p; /* second row */ 88142faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 88242faf41dSJed Brown 88342faf41dSJed Brown { 88442faf41dSJed Brown const PetscReal misfit = a2*a2*bm[1] + a3*a3*bm[2] + a4*a4*bm[3] - one/three; 88542faf41dSJed Brown if (PetscAbs(misfit) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Assumptions violated, could not construct a third order embedded method"); 88642faf41dSJed Brown } 8870298fd71SBarry Smith ierr = TSRosWRegister(name,4,4,&A[0][0],&Gamma[0][0],b,bm,0,NULL);CHKERRQ(ierr); 88842faf41dSJed Brown PetscFunctionReturn(0); 88942faf41dSJed Brown } 89042faf41dSJed Brown 8911c3436cfSJed Brown /* 8921c3436cfSJed Brown The step completion formula is 8931c3436cfSJed Brown 8941c3436cfSJed Brown x1 = x0 + b^T Y 8951c3436cfSJed Brown 8961c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 8971c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 8981c3436cfSJed Brown 8991c3436cfSJed Brown x1e = x0 + be^T Y 9001c3436cfSJed Brown = x1 - b^T Y + be^T Y 9011c3436cfSJed Brown = x1 + (be - b)^T Y 9021c3436cfSJed Brown 9031c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 9041c3436cfSJed Brown */ 905f9c1d6abSBarry Smith static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec U,PetscBool *done) 9061c3436cfSJed Brown { 9071c3436cfSJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 9081c3436cfSJed Brown RosWTableau tab = ros->tableau; 9091c3436cfSJed Brown PetscScalar *w = ros->work; 9101c3436cfSJed Brown PetscInt i; 9111c3436cfSJed Brown PetscErrorCode ierr; 9121c3436cfSJed Brown 9131c3436cfSJed Brown PetscFunctionBegin; 9141c3436cfSJed Brown if (order == tab->order) { 915108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 916f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 917de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bt[i]; 918f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 919f9c1d6abSBarry Smith } else {ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr);} 9201c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9211c3436cfSJed Brown PetscFunctionReturn(0); 9221c3436cfSJed Brown } else if (order == tab->order-1) { 9231c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 924108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 925f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 926de19f811SJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i]; 927f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 928108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 929108c343cSJed Brown for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 930f9c1d6abSBarry Smith ierr = VecCopy(ts->vec_sol,U);CHKERRQ(ierr); 931f9c1d6abSBarry Smith ierr = VecMAXPY(U,tab->s,w,ros->Y);CHKERRQ(ierr); 9321c3436cfSJed Brown } 9331c3436cfSJed Brown if (done) *done = PETSC_TRUE; 9341c3436cfSJed Brown PetscFunctionReturn(0); 9351c3436cfSJed Brown } 9361c3436cfSJed Brown unavailable: 9371c3436cfSJed Brown if (done) *done = PETSC_FALSE; 938a7fac7c2SEmil 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); 9391c3436cfSJed Brown PetscFunctionReturn(0); 9401c3436cfSJed Brown } 9411c3436cfSJed Brown 942560360afSLisandro Dalcin static PetscErrorCode TSRollBack_RosW(TS ts) 94324655328SShri { 94424655328SShri TS_RosW *ros = (TS_RosW*)ts->data; 94524655328SShri PetscErrorCode ierr; 94624655328SShri 94724655328SShri PetscFunctionBegin; 948be5899b3SLisandro Dalcin ierr = VecCopy(ros->vec_sol_prev,ts->vec_sol);CHKERRQ(ierr); 94924655328SShri PetscFunctionReturn(0); 95024655328SShri } 95124655328SShri 952e27a552bSJed Brown static PetscErrorCode TSStep_RosW(TS ts) 953e27a552bSJed Brown { 95461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 95561692a83SJed Brown RosWTableau tab = ros->tableau; 956e27a552bSJed Brown const PetscInt s = tab->s; 9571c3436cfSJed Brown const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 9580feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 959c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 96061692a83SJed Brown PetscScalar *w = ros->work; 9617d4bf2deSEmil Constantinescu Vec *Y = ros->Y,Ydot = ros->Ydot,Zdot = ros->Zdot,Zstage = ros->Zstage; 962e27a552bSJed Brown SNES snes; 9631c3436cfSJed Brown TSAdapt adapt; 964fecfb714SLisandro Dalcin PetscInt i,j,its,lits; 965be5899b3SLisandro Dalcin PetscInt rejections = 0; 966b39943a6SLisandro Dalcin PetscBool stageok,accept = PETSC_TRUE; 967b39943a6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 968e27a552bSJed Brown PetscErrorCode ierr; 969f7f07198SBarry Smith PetscInt lag; 970e27a552bSJed Brown 971e27a552bSJed Brown PetscFunctionBegin; 972b39943a6SLisandro Dalcin if (!ts->steprollback) { 973be5899b3SLisandro Dalcin ierr = VecCopy(ts->vec_sol,ros->vec_sol_prev);CHKERRQ(ierr); 974b39943a6SLisandro Dalcin } 975e27a552bSJed Brown 976b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 977b39943a6SLisandro Dalcin while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 9781c3436cfSJed Brown const PetscReal h = ts->time_step; 979e27a552bSJed Brown for (i=0; i<s; i++) { 9801c3436cfSJed Brown ros->stage_time = ts->ptime + h*ASum[i]; 981b8123daeSJed Brown ierr = TSPreStage(ts,ros->stage_time);CHKERRQ(ierr); 982c17803e7SJed Brown if (GammaZeroDiag[i]) { 983c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 984b296d7d5SJed Brown ros->scoeff = 1.; 985c17803e7SJed Brown } else { 986c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 987b296d7d5SJed Brown ros->scoeff = 1./Gamma[i*s+i]; 988fd96d5b0SEmil Constantinescu } 98961692a83SJed Brown 99061692a83SJed Brown ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 991de19f811SJed Brown for (j=0; j<i; j++) w[j] = At[i*s+j]; 992de19f811SJed Brown ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 99361692a83SJed Brown 99461692a83SJed Brown for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 99561692a83SJed Brown ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 99661692a83SJed Brown ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 99761692a83SJed Brown 998e27a552bSJed Brown /* Initial guess taken from last stage */ 99961692a83SJed Brown ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 100061692a83SJed Brown 10017d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 1002b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 100361692a83SJed Brown if (!ros->recompute_jacobian && !i) { 1004f7f07198SBarry Smith ierr = SNESGetLagJacobian(snes,&lag);CHKERRQ(ierr); 1005f7f07198SBarry Smith if (lag == 1) { /* use did not set a nontrival lag, so lag over all stages */ 1006f7f07198SBarry Smith ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */ 1007f7f07198SBarry Smith } 100861692a83SJed Brown } 10090298fd71SBarry Smith ierr = SNESSolve(snes,NULL,Y[i]);CHKERRQ(ierr); 1010f7f07198SBarry Smith if (!ros->recompute_jacobian && i == s-1 && lag == 1) { 1011f7f07198SBarry Smith ierr = SNESSetLagJacobian(snes,lag);CHKERRQ(ierr); /* Set lag back to 1 so we know user did not set it */ 1012f7f07198SBarry Smith } 1013e27a552bSJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 1014e27a552bSJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 10155ef26d82SJed Brown ts->snes_its += its; ts->ksp_its += lits; 10167d4bf2deSEmil Constantinescu } else { 10171ce71dffSSatish Balay Mat J,Jp; 10180feba352SEmil Constantinescu ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 10190feba352SEmil Constantinescu ierr = TSComputeIFunction(ts,ros->stage_time,Zstage,Ydot,Y[i],PETSC_FALSE);CHKERRQ(ierr); 102022d28d08SBarry Smith ierr = VecScale(Y[i],-1.0);CHKERRQ(ierr); 10210feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 10220feba352SEmil Constantinescu 10230feba352SEmil Constantinescu ierr = VecZeroEntries(Zstage);CHKERRQ(ierr); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 10240feba352SEmil Constantinescu for (j=0; j<i; j++) w[j] = GammaExplicitCorr[i*s+j]; 10250feba352SEmil Constantinescu ierr = VecMAXPY(Zstage,i,w,Y);CHKERRQ(ierr); 1026fecfb714SLisandro Dalcin 1027fecfb714SLisandro Dalcin /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 10280298fd71SBarry Smith ierr = TSGetIJacobian(ts,&J,&Jp,NULL,NULL);CHKERRQ(ierr); 1029d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,ts->vec_sol,Ydot,0,J,Jp,PETSC_FALSE);CHKERRQ(ierr); 103022d28d08SBarry Smith ierr = MatMult(J,Zstage,Zdot);CHKERRQ(ierr); 10310feba352SEmil Constantinescu ierr = VecAXPY(Y[i],-1.0,Zdot);CHKERRQ(ierr); 10325ef26d82SJed Brown ts->ksp_its += 1; 1033fecfb714SLisandro Dalcin 1034fecfb714SLisandro Dalcin ierr = VecScale(Y[i],h);CHKERRQ(ierr); 10357d4bf2deSEmil Constantinescu } 10369be3e283SDebojyoti Ghosh ierr = TSPostStage(ts,ros->stage_time,i,Y);CHKERRQ(ierr); 1037fecfb714SLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1038fecfb714SLisandro Dalcin ierr = TSAdaptCheckStage(adapt,ts,ros->stage_time,Y[i],&stageok);CHKERRQ(ierr); 1039fecfb714SLisandro Dalcin if (!stageok) goto reject_step; 1040e27a552bSJed Brown } 1041e27a552bSJed Brown 1042b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 1043b39943a6SLisandro Dalcin ierr = TSEvaluateStep_RosW(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr); 1044b39943a6SLisandro Dalcin ros->status = TS_STEP_PENDING; 1045552698daSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 10461c3436cfSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 10471917a363SLisandro Dalcin ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr); 1048fecfb714SLisandro Dalcin ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 1049b39943a6SLisandro Dalcin ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1050b39943a6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 1051b39943a6SLisandro Dalcin ierr = TSRollBack_RosW(ts);CHKERRQ(ierr); 1052be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1053be5899b3SLisandro Dalcin goto reject_step; 1054b39943a6SLisandro Dalcin } 1055b39943a6SLisandro Dalcin 1056e27a552bSJed Brown ts->ptime += ts->time_step; 1057cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 10581c3436cfSJed Brown break; 1059b39943a6SLisandro Dalcin 1060b39943a6SLisandro Dalcin reject_step: 1061fecfb714SLisandro Dalcin ts->reject++; accept = PETSC_FALSE; 1062be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1063b39943a6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 1064be5899b3SLisandro Dalcin ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 10651c3436cfSJed Brown } 10661c3436cfSJed Brown } 1067e27a552bSJed Brown PetscFunctionReturn(0); 1068e27a552bSJed Brown } 1069e27a552bSJed Brown 1070f9c1d6abSBarry Smith static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec U) 1071e27a552bSJed Brown { 107261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1073f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s,pinterp = ros->tableau->pinterp,i,j; 1074f4aed992SEmil Constantinescu PetscReal h; 1075f4aed992SEmil Constantinescu PetscReal tt,t; 1076f4aed992SEmil Constantinescu PetscScalar *bt; 1077f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1078f4aed992SEmil Constantinescu PetscErrorCode ierr; 1079f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1080f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1081f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1082e27a552bSJed Brown 1083e27a552bSJed Brown PetscFunctionBegin; 1084ce94432eSBarry Smith if (!Bt) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 1085f4aed992SEmil Constantinescu 1086f4aed992SEmil Constantinescu switch (ros->status) { 1087f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1088f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1089f4aed992SEmil Constantinescu h = ts->time_step; 1090f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h; 1091f4aed992SEmil Constantinescu break; 1092f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1093be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1094f4aed992SEmil Constantinescu t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1095f4aed992SEmil Constantinescu break; 1096ce94432eSBarry Smith default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1097f4aed992SEmil Constantinescu } 1098785e854fSJed Brown ierr = PetscMalloc1(s,&bt);CHKERRQ(ierr); 1099f4aed992SEmil Constantinescu for (i=0; i<s; i++) bt[i] = 0; 1100f4aed992SEmil Constantinescu for (j=0,tt=t; j<pinterp; j++,tt*=t) { 1101f4aed992SEmil Constantinescu for (i=0; i<s; i++) { 11023ca35412SEmil Constantinescu bt[i] += Bt[i*pinterp+j] * tt; 1103f4aed992SEmil Constantinescu } 1104f4aed992SEmil Constantinescu } 1105f4aed992SEmil Constantinescu 1106f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1107f9c1d6abSBarry Smith /* U <- 0*/ 1108f9c1d6abSBarry Smith ierr = VecZeroEntries(U);CHKERRQ(ierr); 1109f9c1d6abSBarry Smith /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 11103ca35412SEmil Constantinescu for (j=0; j<s; j++) w[j] = 0; 11113ca35412SEmil Constantinescu for (j=0; j<s; j++) { 11123ca35412SEmil Constantinescu for (i=j; i<s; i++) { 11133ca35412SEmil Constantinescu w[j] += bt[i]*GammaInv[i*s+j]; 1114f4aed992SEmil Constantinescu } 11153ca35412SEmil Constantinescu } 1116f9c1d6abSBarry Smith ierr = VecMAXPY(U,i,w,Y);CHKERRQ(ierr); 1117be5899b3SLisandro Dalcin /* U <- y(t) + U */ 1118be5899b3SLisandro Dalcin ierr = VecAXPY(U,1,ros->vec_sol_prev);CHKERRQ(ierr); 1119f4aed992SEmil Constantinescu 1120f4aed992SEmil Constantinescu ierr = PetscFree(bt);CHKERRQ(ierr); 1121e27a552bSJed Brown PetscFunctionReturn(0); 1122e27a552bSJed Brown } 1123e27a552bSJed Brown 1124e27a552bSJed Brown /*------------------------------------------------------------*/ 1125b39943a6SLisandro Dalcin 1126b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauReset(TS ts) 1127b39943a6SLisandro Dalcin { 1128b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1129b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1130b39943a6SLisandro Dalcin PetscErrorCode ierr; 1131b39943a6SLisandro Dalcin 1132b39943a6SLisandro Dalcin PetscFunctionBegin; 1133b39943a6SLisandro Dalcin if (!tab) PetscFunctionReturn(0); 1134b39943a6SLisandro Dalcin ierr = VecDestroyVecs(tab->s,&ros->Y);CHKERRQ(ierr); 1135b39943a6SLisandro Dalcin ierr = PetscFree(ros->work);CHKERRQ(ierr); 1136b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1137b39943a6SLisandro Dalcin } 1138b39943a6SLisandro Dalcin 1139e27a552bSJed Brown static PetscErrorCode TSReset_RosW(TS ts) 1140e27a552bSJed Brown { 114161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1142e27a552bSJed Brown PetscErrorCode ierr; 1143e27a552bSJed Brown 1144e27a552bSJed Brown PetscFunctionBegin; 1145b39943a6SLisandro Dalcin ierr = TSRosWTableauReset(ts);CHKERRQ(ierr); 114661692a83SJed Brown ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 114761692a83SJed Brown ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 114861692a83SJed Brown ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 114961692a83SJed Brown ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 1150be5899b3SLisandro Dalcin ierr = VecDestroy(&ros->vec_sol_prev);CHKERRQ(ierr); 1151e27a552bSJed Brown PetscFunctionReturn(0); 1152e27a552bSJed Brown } 1153e27a552bSJed Brown 1154d5e6173cSPeter Brune static PetscErrorCode TSRosWGetVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot,Vec *Ystage,Vec *Zstage) 1155d5e6173cSPeter Brune { 1156d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW*)ts->data; 1157d5e6173cSPeter Brune PetscErrorCode ierr; 1158d5e6173cSPeter Brune 1159d5e6173cSPeter Brune PetscFunctionBegin; 1160d5e6173cSPeter Brune if (Ydot) { 1161d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1162d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1163d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1164d5e6173cSPeter Brune } 1165d5e6173cSPeter Brune if (Zdot) { 1166d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1167d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1168d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1169d5e6173cSPeter Brune } 1170d5e6173cSPeter Brune if (Ystage) { 1171d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1172d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1173d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1174d5e6173cSPeter Brune } 1175d5e6173cSPeter Brune if (Zstage) { 1176d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1177d5e6173cSPeter Brune ierr = DMGetNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1178d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1179d5e6173cSPeter Brune } 1180d5e6173cSPeter Brune PetscFunctionReturn(0); 1181d5e6173cSPeter Brune } 1182d5e6173cSPeter Brune 1183d5e6173cSPeter Brune static PetscErrorCode TSRosWRestoreVecs(TS ts,DM dm,Vec *Ydot,Vec *Zdot, Vec *Ystage, Vec *Zstage) 1184d5e6173cSPeter Brune { 1185d5e6173cSPeter Brune PetscErrorCode ierr; 1186d5e6173cSPeter Brune 1187d5e6173cSPeter Brune PetscFunctionBegin; 1188d5e6173cSPeter Brune if (Ydot) { 1189d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1190d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ydot",Ydot);CHKERRQ(ierr); 1191d5e6173cSPeter Brune } 1192d5e6173cSPeter Brune } 1193d5e6173cSPeter Brune if (Zdot) { 1194d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1195d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zdot",Zdot);CHKERRQ(ierr); 1196d5e6173cSPeter Brune } 1197d5e6173cSPeter Brune } 1198d5e6173cSPeter Brune if (Ystage) { 1199d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1200d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Ystage",Ystage);CHKERRQ(ierr); 1201d5e6173cSPeter Brune } 1202d5e6173cSPeter Brune } 1203d5e6173cSPeter Brune if (Zstage) { 1204d5e6173cSPeter Brune if (dm && dm != ts->dm) { 1205d5e6173cSPeter Brune ierr = DMRestoreNamedGlobalVector(dm,"TSRosW_Zstage",Zstage);CHKERRQ(ierr); 1206d5e6173cSPeter Brune } 1207d5e6173cSPeter Brune } 1208d5e6173cSPeter Brune PetscFunctionReturn(0); 1209d5e6173cSPeter Brune } 1210d5e6173cSPeter Brune 1211d5e6173cSPeter Brune static PetscErrorCode DMCoarsenHook_TSRosW(DM fine,DM coarse,void *ctx) 1212d5e6173cSPeter Brune { 1213d5e6173cSPeter Brune PetscFunctionBegin; 1214d5e6173cSPeter Brune PetscFunctionReturn(0); 1215d5e6173cSPeter Brune } 1216d5e6173cSPeter Brune 1217d5e6173cSPeter Brune static PetscErrorCode DMRestrictHook_TSRosW(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1218d5e6173cSPeter Brune { 1219d5e6173cSPeter Brune TS ts = (TS)ctx; 1220d5e6173cSPeter Brune PetscErrorCode ierr; 1221d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1222d5e6173cSPeter Brune Vec Ydotc,Zdotc,Ystagec,Zstagec; 1223d5e6173cSPeter Brune 1224d5e6173cSPeter Brune PetscFunctionBegin; 1225d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1226d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1227d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ydot,Ydotc);CHKERRQ(ierr); 1228d5e6173cSPeter Brune ierr = VecPointwiseMult(Ydotc,rscale,Ydotc);CHKERRQ(ierr); 1229d5e6173cSPeter Brune ierr = MatRestrict(restrct,Ystage,Ystagec);CHKERRQ(ierr); 1230d5e6173cSPeter Brune ierr = VecPointwiseMult(Ystagec,rscale,Ystagec);CHKERRQ(ierr); 1231d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zdot,Zdotc);CHKERRQ(ierr); 1232d5e6173cSPeter Brune ierr = VecPointwiseMult(Zdotc,rscale,Zdotc);CHKERRQ(ierr); 1233d5e6173cSPeter Brune ierr = MatRestrict(restrct,Zstage,Zstagec);CHKERRQ(ierr); 1234d5e6173cSPeter Brune ierr = VecPointwiseMult(Zstagec,rscale,Zstagec);CHKERRQ(ierr); 1235d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,fine,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1236d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,coarse,&Ydotc,&Ystagec,&Zdotc,&Zstagec);CHKERRQ(ierr); 1237d5e6173cSPeter Brune PetscFunctionReturn(0); 1238d5e6173cSPeter Brune } 1239d5e6173cSPeter Brune 1240258e1594SPeter Brune static PetscErrorCode DMSubDomainHook_TSRosW(DM fine,DM coarse,void *ctx) 1241258e1594SPeter Brune { 1242258e1594SPeter Brune PetscFunctionBegin; 1243258e1594SPeter Brune PetscFunctionReturn(0); 1244258e1594SPeter Brune } 1245258e1594SPeter Brune 1246258e1594SPeter Brune static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1247258e1594SPeter Brune { 1248258e1594SPeter Brune TS ts = (TS)ctx; 1249258e1594SPeter Brune PetscErrorCode ierr; 1250258e1594SPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1251258e1594SPeter Brune Vec Ydots,Zdots,Ystages,Zstages; 1252258e1594SPeter Brune 1253258e1594SPeter Brune PetscFunctionBegin; 1254258e1594SPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1255258e1594SPeter Brune ierr = TSRosWGetVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1256258e1594SPeter Brune 1257258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1258258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ydot,Ydots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1259258e1594SPeter Brune 1260258e1594SPeter Brune ierr = VecScatterBegin(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1261258e1594SPeter Brune ierr = VecScatterEnd(gscat,Ystage,Ystages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1262258e1594SPeter Brune 1263258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1264258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zdot,Zdots,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1265258e1594SPeter Brune 1266258e1594SPeter Brune ierr = VecScatterBegin(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1267258e1594SPeter Brune ierr = VecScatterEnd(gscat,Zstage,Zstages,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1268258e1594SPeter Brune 1269258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Ystage,&Zdot,&Zstage);CHKERRQ(ierr); 1270258e1594SPeter Brune ierr = TSRosWRestoreVecs(ts,subdm,&Ydots,&Ystages,&Zdots,&Zstages);CHKERRQ(ierr); 1271258e1594SPeter Brune PetscFunctionReturn(0); 1272258e1594SPeter Brune } 1273258e1594SPeter Brune 1274e27a552bSJed Brown /* 1275e27a552bSJed Brown This defines the nonlinear equation that is to be solved with SNES 1276e27a552bSJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 1277e27a552bSJed Brown */ 1278f9c1d6abSBarry Smith static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec U,Vec F,TS ts) 1279e27a552bSJed Brown { 128061692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1281e27a552bSJed Brown PetscErrorCode ierr; 1282d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1283b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1284d5e6173cSPeter Brune DM dm,dmsave; 1285e27a552bSJed Brown 1286e27a552bSJed Brown PetscFunctionBegin; 1287d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1288d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1289b296d7d5SJed Brown ierr = VecWAXPY(Ydot,shift,U,Zdot);CHKERRQ(ierr); /* Ydot = shift*U + Zdot */ 1290f9c1d6abSBarry Smith ierr = VecWAXPY(Ystage,1.0,U,Zstage);CHKERRQ(ierr); /* Ystage = U + Zstage */ 1291d5e6173cSPeter Brune dmsave = ts->dm; 1292d5e6173cSPeter Brune ts->dm = dm; 1293d5e6173cSPeter Brune ierr = TSComputeIFunction(ts,ros->stage_time,Ystage,Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 1294d5e6173cSPeter Brune ts->dm = dmsave; 1295d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1296e27a552bSJed Brown PetscFunctionReturn(0); 1297e27a552bSJed Brown } 1298e27a552bSJed Brown 1299d1e9a80fSBarry Smith static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec U,Mat A,Mat B,TS ts) 1300e27a552bSJed Brown { 130161692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1302d5e6173cSPeter Brune Vec Ydot,Zdot,Ystage,Zstage; 1303b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1304e27a552bSJed Brown PetscErrorCode ierr; 1305d5e6173cSPeter Brune DM dm,dmsave; 1306e27a552bSJed Brown 1307e27a552bSJed Brown PetscFunctionBegin; 130861692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 1309d5e6173cSPeter Brune ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1310d5e6173cSPeter Brune ierr = TSRosWGetVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1311d5e6173cSPeter Brune dmsave = ts->dm; 1312d5e6173cSPeter Brune ts->dm = dm; 1313d1e9a80fSBarry Smith ierr = TSComputeIJacobian(ts,ros->stage_time,Ystage,Ydot,shift,A,B,PETSC_TRUE);CHKERRQ(ierr); 1314d5e6173cSPeter Brune ts->dm = dmsave; 1315d5e6173cSPeter Brune ierr = TSRosWRestoreVecs(ts,dm,&Ydot,&Zdot,&Ystage,&Zstage);CHKERRQ(ierr); 1316e27a552bSJed Brown PetscFunctionReturn(0); 1317e27a552bSJed Brown } 1318e27a552bSJed Brown 1319b39943a6SLisandro Dalcin static PetscErrorCode TSRosWTableauSetUp(TS ts) 1320b39943a6SLisandro Dalcin { 1321b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW*)ts->data; 1322b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1323b39943a6SLisandro Dalcin PetscErrorCode ierr; 1324b39943a6SLisandro Dalcin 1325b39943a6SLisandro Dalcin PetscFunctionBegin; 1326b39943a6SLisandro Dalcin ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ros->Y);CHKERRQ(ierr); 1327b39943a6SLisandro Dalcin ierr = PetscMalloc1(tab->s,&ros->work);CHKERRQ(ierr); 1328b39943a6SLisandro Dalcin PetscFunctionReturn(0); 1329b39943a6SLisandro Dalcin } 1330b39943a6SLisandro Dalcin 1331e27a552bSJed Brown static PetscErrorCode TSSetUp_RosW(TS ts) 1332e27a552bSJed Brown { 133361692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1334e27a552bSJed Brown PetscErrorCode ierr; 1335d5e6173cSPeter Brune DM dm; 1336b39943a6SLisandro Dalcin SNES snes; 1337a3ab5968SHong Zhang TSRHSJacobian rhsjacobian; 1338e27a552bSJed Brown 1339e27a552bSJed Brown PetscFunctionBegin; 1340b39943a6SLisandro Dalcin ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr); 134161692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 134261692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 134361692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 134461692a83SJed Brown ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 1345be5899b3SLisandro Dalcin ierr = VecDuplicate(ts->vec_sol,&ros->vec_sol_prev);CHKERRQ(ierr); 134622d28d08SBarry Smith ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1347d5e6173cSPeter Brune ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1348258e1594SPeter Brune ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1349b39943a6SLisandro Dalcin /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 1350b39943a6SLisandro Dalcin ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1351b39943a6SLisandro Dalcin if (!((PetscObject)snes)->type_name) { 1352b39943a6SLisandro Dalcin ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1353b39943a6SLisandro Dalcin } 1354a3ab5968SHong Zhang ierr = DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);CHKERRQ(ierr); 1355a3ab5968SHong Zhang if (rhsjacobian == TSComputeRHSJacobianConstant) { 1356a3ab5968SHong Zhang Mat Amat,Pmat; 1357a3ab5968SHong Zhang 1358a3ab5968SHong Zhang /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */ 1359a3ab5968SHong Zhang ierr = SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);CHKERRQ(ierr); 1360a3ab5968SHong Zhang if (Amat && Amat == ts->Arhs) { 1361a3ab5968SHong Zhang if (Amat == Pmat) { 1362a3ab5968SHong Zhang ierr = MatDuplicate(ts->Arhs,MAT_COPY_VALUES,&Amat);CHKERRQ(ierr); 1363a3ab5968SHong Zhang ierr = SNESSetJacobian(snes,Amat,Amat,NULL,NULL);CHKERRQ(ierr); 1364a3ab5968SHong Zhang } else { 1365a3ab5968SHong Zhang ierr = MatDuplicate(ts->Arhs,MAT_COPY_VALUES,&Amat);CHKERRQ(ierr); 1366a3ab5968SHong Zhang ierr = SNESSetJacobian(snes,Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1367a3ab5968SHong Zhang if (Pmat && Pmat == ts->Brhs) { 1368a3ab5968SHong Zhang ierr = MatDuplicate(ts->Brhs,MAT_COPY_VALUES,&Pmat);CHKERRQ(ierr); 1369a3ab5968SHong Zhang ierr = SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);CHKERRQ(ierr); 1370a3ab5968SHong Zhang ierr = MatDestroy(&Pmat);CHKERRQ(ierr); 1371a3ab5968SHong Zhang } 1372a3ab5968SHong Zhang } 1373a3ab5968SHong Zhang ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1374a3ab5968SHong Zhang } 1375a3ab5968SHong Zhang } 1376e27a552bSJed Brown PetscFunctionReturn(0); 1377e27a552bSJed Brown } 1378e27a552bSJed Brown /*------------------------------------------------------------*/ 1379e27a552bSJed Brown 13804416b707SBarry Smith static PetscErrorCode TSSetFromOptions_RosW(PetscOptionItems *PetscOptionsObject,TS ts) 1381e27a552bSJed Brown { 138261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1383e27a552bSJed Brown PetscErrorCode ierr; 1384b39943a6SLisandro Dalcin SNES snes; 1385e27a552bSJed Brown 1386e27a552bSJed Brown PetscFunctionBegin; 1387e55864a3SBarry Smith ierr = PetscOptionsHead(PetscOptionsObject,"RosW ODE solver options");CHKERRQ(ierr); 1388e27a552bSJed Brown { 138961692a83SJed Brown RosWTableauLink link; 1390e27a552bSJed Brown PetscInt count,choice; 1391e27a552bSJed Brown PetscBool flg; 1392e27a552bSJed Brown const char **namelist; 139361692a83SJed Brown 139461692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 1395f489ac74SBarry Smith ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 139661692a83SJed Brown for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1397b39943a6SLisandro Dalcin ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,ros->tableau->name,&choice,&flg);CHKERRQ(ierr); 1398b39943a6SLisandro Dalcin if (flg) {ierr = TSRosWSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1399e27a552bSJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 140061692a83SJed Brown 14010298fd71SBarry Smith ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,NULL);CHKERRQ(ierr); 1402b39943a6SLisandro Dalcin } 1403b39943a6SLisandro Dalcin ierr = PetscOptionsTail();CHKERRQ(ierr); 140461692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 140561692a83SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 140661692a83SJed Brown if (!((PetscObject)snes)->type_name) { 140761692a83SJed Brown ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 140861692a83SJed Brown } 1409e27a552bSJed Brown PetscFunctionReturn(0); 1410e27a552bSJed Brown } 1411e27a552bSJed Brown 1412e27a552bSJed Brown static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 1413e27a552bSJed Brown { 141461692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1415e27a552bSJed Brown PetscBool iascii; 1416e27a552bSJed Brown PetscErrorCode ierr; 1417e27a552bSJed Brown 1418e27a552bSJed Brown PetscFunctionBegin; 1419251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1420e27a552bSJed Brown if (iascii) { 14219c334d8fSLisandro Dalcin RosWTableau tab = ros->tableau; 142219fd82e9SBarry Smith TSRosWType rostype; 14239c334d8fSLisandro Dalcin char buf[512]; 1424e408995aSJed Brown PetscInt i; 1425e408995aSJed Brown PetscReal abscissa[512]; 142661692a83SJed Brown ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 142761692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 1428de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 142961692a83SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 1430e408995aSJed Brown for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 1431de043e34SBarry Smith ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 1432e408995aSJed Brown ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 1433e27a552bSJed Brown } 1434e27a552bSJed Brown PetscFunctionReturn(0); 1435e27a552bSJed Brown } 1436e27a552bSJed Brown 14379200755eSBarry Smith static PetscErrorCode TSLoad_RosW(TS ts,PetscViewer viewer) 14389200755eSBarry Smith { 14399200755eSBarry Smith PetscErrorCode ierr; 14409200755eSBarry Smith SNES snes; 14419c334d8fSLisandro Dalcin TSAdapt adapt; 14429200755eSBarry Smith 14439200755eSBarry Smith PetscFunctionBegin; 14449c334d8fSLisandro Dalcin ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 14459c334d8fSLisandro Dalcin ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 14469200755eSBarry Smith ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 14479200755eSBarry Smith ierr = SNESLoad(snes,viewer);CHKERRQ(ierr); 14489200755eSBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 14499200755eSBarry Smith ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr); 14509200755eSBarry Smith ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr); 14519200755eSBarry Smith PetscFunctionReturn(0); 14529200755eSBarry Smith } 14539200755eSBarry Smith 1454e27a552bSJed Brown /*@C 145561692a83SJed Brown TSRosWSetType - Set the type of Rosenbrock-W scheme 1456e27a552bSJed Brown 1457e27a552bSJed Brown Logically collective 1458e27a552bSJed Brown 1459*d8d19677SJose E. Roman Input Parameters: 1460e27a552bSJed Brown + ts - timestepping context 1461b92453a8SLisandro Dalcin - roswtype - type of Rosenbrock-W scheme 1462e27a552bSJed Brown 1463020d8f30SJed Brown Level: beginner 1464e27a552bSJed Brown 1465020d8f30SJed Brown .seealso: TSRosWGetType(), TSROSW, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWARK3 1466e27a552bSJed Brown @*/ 1467b92453a8SLisandro Dalcin PetscErrorCode TSRosWSetType(TS ts,TSRosWType roswtype) 1468e27a552bSJed Brown { 1469e27a552bSJed Brown PetscErrorCode ierr; 1470e27a552bSJed Brown 1471e27a552bSJed Brown PetscFunctionBegin; 1472e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1473b92453a8SLisandro Dalcin PetscValidCharPointer(roswtype,2); 1474b92453a8SLisandro Dalcin ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,TSRosWType),(ts,roswtype));CHKERRQ(ierr); 1475e27a552bSJed Brown PetscFunctionReturn(0); 1476e27a552bSJed Brown } 1477e27a552bSJed Brown 1478e27a552bSJed Brown /*@C 147961692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1480e27a552bSJed Brown 1481e27a552bSJed Brown Logically collective 1482e27a552bSJed Brown 1483e27a552bSJed Brown Input Parameter: 1484e27a552bSJed Brown . ts - timestepping context 1485e27a552bSJed Brown 1486e27a552bSJed Brown Output Parameter: 148761692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1488e27a552bSJed Brown 1489e27a552bSJed Brown Level: intermediate 1490e27a552bSJed Brown 1491e27a552bSJed Brown .seealso: TSRosWGetType() 1492e27a552bSJed Brown @*/ 149319fd82e9SBarry Smith PetscErrorCode TSRosWGetType(TS ts,TSRosWType *rostype) 1494e27a552bSJed Brown { 1495e27a552bSJed Brown PetscErrorCode ierr; 1496e27a552bSJed Brown 1497e27a552bSJed Brown PetscFunctionBegin; 1498e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 149919fd82e9SBarry Smith ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,TSRosWType*),(ts,rostype));CHKERRQ(ierr); 1500e27a552bSJed Brown PetscFunctionReturn(0); 1501e27a552bSJed Brown } 1502e27a552bSJed Brown 1503e27a552bSJed Brown /*@C 150461692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1505e27a552bSJed Brown 1506e27a552bSJed Brown Logically collective 1507e27a552bSJed Brown 1508*d8d19677SJose E. Roman Input Parameters: 1509e27a552bSJed Brown + ts - timestepping context 151061692a83SJed Brown - flg - PETSC_TRUE to recompute the Jacobian at each stage 1511e27a552bSJed Brown 1512e27a552bSJed Brown Level: intermediate 1513e27a552bSJed Brown 1514e27a552bSJed Brown .seealso: TSRosWGetType() 1515e27a552bSJed Brown @*/ 151661692a83SJed Brown PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 1517e27a552bSJed Brown { 1518e27a552bSJed Brown PetscErrorCode ierr; 1519e27a552bSJed Brown 1520e27a552bSJed Brown PetscFunctionBegin; 1521e27a552bSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 152261692a83SJed Brown ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 1523e27a552bSJed Brown PetscFunctionReturn(0); 1524e27a552bSJed Brown } 1525e27a552bSJed Brown 1526560360afSLisandro Dalcin static PetscErrorCode TSRosWGetType_RosW(TS ts,TSRosWType *rostype) 1527e27a552bSJed Brown { 152861692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1529e27a552bSJed Brown 1530e27a552bSJed Brown PetscFunctionBegin; 153161692a83SJed Brown *rostype = ros->tableau->name; 1532e27a552bSJed Brown PetscFunctionReturn(0); 1533e27a552bSJed Brown } 1534ef20d060SBarry Smith 1535560360afSLisandro Dalcin static PetscErrorCode TSRosWSetType_RosW(TS ts,TSRosWType rostype) 1536e27a552bSJed Brown { 153761692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1538e27a552bSJed Brown PetscErrorCode ierr; 1539e27a552bSJed Brown PetscBool match; 154061692a83SJed Brown RosWTableauLink link; 1541e27a552bSJed Brown 1542e27a552bSJed Brown PetscFunctionBegin; 154361692a83SJed Brown if (ros->tableau) { 154461692a83SJed Brown ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 1545e27a552bSJed Brown if (match) PetscFunctionReturn(0); 1546e27a552bSJed Brown } 154761692a83SJed Brown for (link = RosWTableauList; link; link=link->next) { 154861692a83SJed Brown ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 1549e27a552bSJed Brown if (match) { 1550b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauReset(ts);CHKERRQ(ierr);} 155161692a83SJed Brown ros->tableau = &link->tab; 1552b39943a6SLisandro Dalcin if (ts->setupcalled) {ierr = TSRosWTableauSetUp(ts);CHKERRQ(ierr);} 15532ffb9264SLisandro Dalcin ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1554e27a552bSJed Brown PetscFunctionReturn(0); 1555e27a552bSJed Brown } 1556e27a552bSJed Brown } 1557ce94432eSBarry Smith SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 1558e27a552bSJed Brown } 155961692a83SJed Brown 1560560360afSLisandro Dalcin static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 1561e27a552bSJed Brown { 156261692a83SJed Brown TS_RosW *ros = (TS_RosW*)ts->data; 1563e27a552bSJed Brown 1564e27a552bSJed Brown PetscFunctionBegin; 156561692a83SJed Brown ros->recompute_jacobian = flg; 1566e27a552bSJed Brown PetscFunctionReturn(0); 1567e27a552bSJed Brown } 1568e27a552bSJed Brown 1569b3a6b972SJed Brown static PetscErrorCode TSDestroy_RosW(TS ts) 1570b3a6b972SJed Brown { 1571b3a6b972SJed Brown PetscErrorCode ierr; 1572b3a6b972SJed Brown 1573b3a6b972SJed Brown PetscFunctionBegin; 1574b3a6b972SJed Brown ierr = TSReset_RosW(ts);CHKERRQ(ierr); 1575b3a6b972SJed Brown if (ts->dm) { 1576b3a6b972SJed Brown ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRosW,DMRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1577b3a6b972SJed Brown ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRosW,DMSubDomainRestrictHook_TSRosW,ts);CHKERRQ(ierr); 1578b3a6b972SJed Brown } 1579b3a6b972SJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 1580b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",NULL);CHKERRQ(ierr); 1581b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",NULL);CHKERRQ(ierr); 1582b3a6b972SJed Brown ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",NULL);CHKERRQ(ierr); 1583b3a6b972SJed Brown PetscFunctionReturn(0); 1584b3a6b972SJed Brown } 1585d5e6173cSPeter Brune 1586e27a552bSJed Brown /* ------------------------------------------------------------ */ 1587e27a552bSJed Brown /*MC 1588020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1589e27a552bSJed Brown 1590e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1591e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1592e27a552bSJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 1593e27a552bSJed Brown 1594e27a552bSJed Brown Notes: 159561692a83SJed Brown This method currently only works with autonomous ODE and DAE. 159661692a83SJed Brown 1597d0685a90SJed Brown Consider trying TSARKIMEX if the stiff part is strongly nonlinear. 1598d0685a90SJed Brown 15993d5a8a6aSBarry 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 16003d5a8a6aSBarry Smith 1601e94cfbe0SPatrick Sanan Developer Notes: 160261692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 160361692a83SJed Brown 1604f9c1d6abSBarry Smith $ udot = f(u) 160561692a83SJed Brown 160661692a83SJed Brown by the stage equations 160761692a83SJed Brown 1608f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 160961692a83SJed Brown 161061692a83SJed Brown and step completion formula 161161692a83SJed Brown 1612f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 161361692a83SJed Brown 1614f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 161561692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 161661692a83SJed Brown we define new variables for the stage equations 161761692a83SJed Brown 161861692a83SJed Brown $ y_i = gamma_ij k_j 161961692a83SJed Brown 162061692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 162161692a83SJed Brown 1622b70472e2SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1} 162361692a83SJed Brown 162461692a83SJed Brown to rewrite the method as 162561692a83SJed Brown 1626f9c1d6abSBarry 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 1627f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j bt_j y_j 162861692a83SJed Brown 162961692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 163061692a83SJed Brown 163161692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 163261692a83SJed Brown 163361692a83SJed Brown or, more compactly in tensor notation 163461692a83SJed Brown 163561692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 163661692a83SJed Brown 163761692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 163861692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 163961692a83SJed Brown equation 164061692a83SJed Brown 1641f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 164261692a83SJed Brown 164361692a83SJed Brown with initial guess y_i = 0. 1644e27a552bSJed Brown 1645e27a552bSJed Brown Level: beginner 1646e27a552bSJed Brown 1647d0685a90SJed Brown .seealso: TSCreate(), TS, TSSetType(), TSRosWSetType(), TSRosWRegister(), TSROSWTHETA1, TSROSWTHETA2, TSROSW2M, TSROSW2P, TSROSWRA3PW, TSROSWRA34PW2, TSROSWRODAS3, 1648a4386c9eSJed Brown TSROSWSANDU3, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSROSWLLSSP3P4S2C, TSROSWGRK4T, TSROSWSHAMP4, TSROSWVELDD4, TSROSW4L 1649e27a552bSJed Brown M*/ 16508cc058d9SJed Brown PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1651e27a552bSJed Brown { 165261692a83SJed Brown TS_RosW *ros; 1653e27a552bSJed Brown PetscErrorCode ierr; 1654e27a552bSJed Brown 1655e27a552bSJed Brown PetscFunctionBegin; 1656607a6623SBarry Smith ierr = TSRosWInitializePackage();CHKERRQ(ierr); 1657e27a552bSJed Brown 1658e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1659e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1660e27a552bSJed Brown ts->ops->view = TSView_RosW; 16619200755eSBarry Smith ts->ops->load = TSLoad_RosW; 1662e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1663e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1664e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 16651c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 166624655328SShri ts->ops->rollback = TSRollBack_RosW; 1667e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1668e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1669e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1670e27a552bSJed Brown 1671efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 1672efd4aadfSBarry Smith 1673b00a9115SJed Brown ierr = PetscNewLog(ts,&ros);CHKERRQ(ierr); 167461692a83SJed Brown ts->data = (void*)ros; 1675e27a552bSJed Brown 1676bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWGetType_C",TSRosWGetType_RosW);CHKERRQ(ierr); 1677bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetType_C",TSRosWSetType_RosW);CHKERRQ(ierr); 1678bdf89e91SBarry Smith ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRosWSetRecomputeJacobian_C",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1679b39943a6SLisandro Dalcin 1680b39943a6SLisandro Dalcin ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 1681e27a552bSJed Brown PetscFunctionReturn(0); 1682e27a552bSJed Brown } 1683