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 731cc06b55SBarry Smith .seealso: [](ch_ts), `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 831cc06b55SBarry Smith .seealso: [](ch_ts), `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 931cc06b55SBarry Smith .seealso: [](ch_ts), `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 1031cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW` 104fe7e6d57SJed Brown M*/ 105fe7e6d57SJed Brown 106fe7e6d57SJed Brown /*MC 1071d27aa22SBarry Smith TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005` 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 113bcf0153eSBarry Smith Level: intermediate 114bcf0153eSBarry Smith 1151cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW` 116fe7e6d57SJed Brown M*/ 117fe7e6d57SJed Brown 118fe7e6d57SJed Brown /*MC 1191d27aa22SBarry Smith TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang_2005`. 120fe7e6d57SJed Brown 121fe7e6d57SJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 122fe7e6d57SJed Brown 123fe7e6d57SJed 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. 124fe7e6d57SJed Brown 125bcf0153eSBarry Smith Level: intermediate 126bcf0153eSBarry Smith 1271cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW` 128fe7e6d57SJed Brown M*/ 129fe7e6d57SJed Brown 130ef3c5b88SJed Brown /*MC 131337ef93bSJed Brown TSROSWR34PRW - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`. 132337ef93bSJed Brown 133337ef93bSJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 134337ef93bSJed Brown 135337ef93bSJed Brown This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.25. 136337ef93bSJed Brown This method is B_{PR} consistent of order 3. 137337ef93bSJed Brown This method is spelled "ROS34PRw" in the paper, an improvement to an earlier "ROS34PRW" method from the same author with B_{PR} order 2. 138337ef93bSJed Brown 139337ef93bSJed Brown Level: intermediate 140337ef93bSJed Brown 141337ef93bSJed Brown .seealso: [](ch_ts), `TSROSW` 142337ef93bSJed Brown M*/ 143337ef93bSJed Brown 144337ef93bSJed Brown /*MC 145337ef93bSJed Brown TSROSWR3PRL2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1 {cite}`rang2015improved`. 146337ef93bSJed Brown 147337ef93bSJed Brown Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 148337ef93bSJed Brown 149337ef93bSJed Brown This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.25. 150337ef93bSJed Brown This method is B_{PR} consistent of order 3. 151337ef93bSJed Brown 152337ef93bSJed Brown Level: intermediate 153337ef93bSJed Brown 154337ef93bSJed Brown .seealso: [](ch_ts), `TSROSW` 155337ef93bSJed Brown M*/ 156337ef93bSJed Brown 157337ef93bSJed Brown /*MC 1581d27aa22SBarry Smith TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme {cite}`sandu_1997` 159ef3c5b88SJed Brown 160ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 161ef3c5b88SJed Brown 162ef3c5b88SJed Brown Both the third order and embedded second order methods are stiffly accurate and L-stable. 163ef3c5b88SJed Brown 164bcf0153eSBarry Smith Level: intermediate 165bcf0153eSBarry Smith 1661cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWSANDU3` 167ef3c5b88SJed Brown M*/ 168ef3c5b88SJed Brown 169ef3c5b88SJed Brown /*MC 17048665b53SJed Brown TSROSWRODASPR - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved` 17148665b53SJed Brown 17248665b53SJed Brown By default, the Jacobian is only recomputed once per step. 17348665b53SJed Brown 17448665b53SJed Brown Both the fourth order and embedded third order methods are stiffly accurate and L-stable. 17548665b53SJed Brown The method is B_{PR} consistent of order 3, which ensures convergence order for non-stiff, medium stiff, and stiff problems. 17648665b53SJed Brown 17748665b53SJed Brown Level: intermediate 17848665b53SJed Brown 179337ef93bSJed Brown .seealso: [](ch_ts), `TSROSW`, `TSROSWR34PRW`, `TSROSWR3PRL2` 180337ef93bSJed Brown M*/ 181337ef93bSJed Brown 182337ef93bSJed Brown /*MC 183337ef93bSJed Brown TSROSWRODASPR2 - Six stage fourth order L-stable Rosenbrock scheme {cite}`rang2015improved` 184337ef93bSJed Brown 185337ef93bSJed Brown By default, the Jacobian is only recomputed once per step. 186337ef93bSJed Brown 187337ef93bSJed Brown Both the fourth order and embedded third order methods are stiffly accurate and L-stable. 188337ef93bSJed Brown The method is B_{PR} consistent of order 3, which ensures convergence order for non-stiff, medium stiff, and stiff problems. 189337ef93bSJed Brown This method is similar to `TSROSWRODASPR`, but satisfies one extra B_{PR} order condition. 190337ef93bSJed Brown 191337ef93bSJed Brown Developer Note: 192337ef93bSJed Brown In numerical experiments with ts/tutorials/ex22.c, I (Jed) find this to produce surprisingly poor results. 193337ef93bSJed Brown Although the coefficients pass basic smoke tests, I'm not confident it was tabulated correctly in the paper. 194337ef93bSJed Brown It would be informative if someone could reproduce tests from the paper and/or reach out to the author to understand why it fails on this test problem. 195337ef93bSJed Brown If the method is implemented correctly, doing so might shed light on an additional analysis lens (or further conditions) for robustness on such problems. 196337ef93bSJed Brown 197337ef93bSJed Brown Level: intermediate 198337ef93bSJed Brown 199337ef93bSJed Brown .seealso: [](ch_ts), `TSROSW`, `TSROSWRODASPR` 20048665b53SJed Brown M*/ 20148665b53SJed Brown 20248665b53SJed Brown /*MC 2031d27aa22SBarry Smith TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme {cite}`sandu_1997` 204ef3c5b88SJed Brown 205ef3c5b88SJed Brown By default, the Jacobian is only recomputed once per step. 206ef3c5b88SJed Brown 207ef3c5b88SJed Brown The third order method is L-stable, but not stiffly accurate. 208ef3c5b88SJed Brown The second order embedded method is strongly A-stable with R(infty) = 0.5. 209ef3c5b88SJed Brown The internal stages are L-stable. 2101d27aa22SBarry Smith This method is called ROS3 in {cite}`sandu_1997`. 211ef3c5b88SJed Brown 212bcf0153eSBarry Smith Level: intermediate 213bcf0153eSBarry Smith 2141cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWRODAS3` 215ef3c5b88SJed Brown M*/ 216ef3c5b88SJed Brown 217961f28d0SJed Brown /*MC 218961f28d0SJed Brown TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 219961f28d0SJed Brown 220961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 221961f28d0SJed Brown 222961f28d0SJed Brown A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 223961f28d0SJed Brown 224bcf0153eSBarry Smith Level: intermediate 225bcf0153eSBarry Smith 2261cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `SSP` 227961f28d0SJed Brown M*/ 228961f28d0SJed Brown 229961f28d0SJed Brown /*MC 230998eb97aSJed Brown TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 231961f28d0SJed Brown 232961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 233961f28d0SJed Brown 234961f28d0SJed Brown L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 235961f28d0SJed Brown 236bcf0153eSBarry Smith Level: intermediate 237bcf0153eSBarry Smith 2381cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLLSSP3P4S2C`, `TSSSP` 239961f28d0SJed Brown M*/ 240961f28d0SJed Brown 241961f28d0SJed Brown /*MC 242998eb97aSJed Brown TSROSWLLSSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, four stages 243961f28d0SJed Brown 244961f28d0SJed Brown By default, the Jacobian is only recomputed once per step. 245961f28d0SJed Brown 246961f28d0SJed Brown L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 247961f28d0SJed Brown 248bcf0153eSBarry Smith Level: intermediate 249bcf0153eSBarry Smith 2501cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSSSP` 251961f28d0SJed Brown M*/ 252961f28d0SJed Brown 25342faf41dSJed Brown /*MC 2541d27aa22SBarry Smith TSROSWGRK4T - four stage, fourth order Rosenbrock (not W) method from Kaps and Rentrop {cite}`kaps1979generalized` 25542faf41dSJed Brown 25642faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 25742faf41dSJed Brown 25842faf41dSJed Brown A(89.3 degrees)-stable, |R(infty)| = 0.454. 25942faf41dSJed Brown 26042faf41dSJed Brown This method does not provide a dense output formula. 26142faf41dSJed Brown 262bcf0153eSBarry Smith Level: intermediate 263bcf0153eSBarry Smith 2641d27aa22SBarry Smith Note: 2651d27aa22SBarry Smith See Section 4 Table 7.2 in {cite}`wanner1996solving` 26642faf41dSJed Brown 2671cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L` 26842faf41dSJed Brown M*/ 26942faf41dSJed Brown 27042faf41dSJed Brown /*MC 2711d27aa22SBarry Smith TSROSWSHAMP4 - four stage, fourth order Rosenbrock (not W) method from Shampine {cite}`shampine1982implementation` 27242faf41dSJed Brown 27342faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 27442faf41dSJed Brown 27542faf41dSJed Brown A-stable, |R(infty)| = 1/3. 27642faf41dSJed Brown 27742faf41dSJed Brown This method does not provide a dense output formula. 27842faf41dSJed Brown 279bcf0153eSBarry Smith Level: intermediate 280bcf0153eSBarry Smith 2811d27aa22SBarry Smith Note: 28215229ffcSPierre Jolivet See Section 4 Table 7.2 in {cite}`wanner1996solving` 28342faf41dSJed Brown 2841cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWVELDD4`, `TSROSW4L` 28542faf41dSJed Brown M*/ 28642faf41dSJed Brown 28742faf41dSJed Brown /*MC 2881d27aa22SBarry Smith TSROSWVELDD4 - four stage, fourth order Rosenbrock (not W) method from van Veldhuizen {cite}`veldhuizen1984d` 28942faf41dSJed Brown 29042faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 29142faf41dSJed Brown 29242faf41dSJed Brown A(89.5 degrees)-stable, |R(infty)| = 0.24. 29342faf41dSJed Brown 29442faf41dSJed Brown This method does not provide a dense output formula. 29542faf41dSJed Brown 296bcf0153eSBarry Smith Level: intermediate 297bcf0153eSBarry Smith 2981d27aa22SBarry Smith Note: 2991d27aa22SBarry Smith See Section 4 Table 7.2 in {cite}`wanner1996solving` 30042faf41dSJed Brown 3011cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L` 30242faf41dSJed Brown M*/ 30342faf41dSJed Brown 30442faf41dSJed Brown /*MC 30542faf41dSJed Brown TSROSW4L - four stage, fourth order Rosenbrock (not W) method 30642faf41dSJed Brown 30742faf41dSJed Brown By default, the Jacobian is only recomputed once per step. 30842faf41dSJed Brown 30942faf41dSJed Brown A-stable and L-stable 31042faf41dSJed Brown 31142faf41dSJed Brown This method does not provide a dense output formula. 31242faf41dSJed Brown 313bcf0153eSBarry Smith Level: intermediate 314bcf0153eSBarry Smith 3151d27aa22SBarry Smith Note: 31615229ffcSPierre Jolivet See Section 4 Table 7.2 in {cite}`wanner1996solving` 31742faf41dSJed Brown 3181cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSW4L` 31942faf41dSJed Brown M*/ 32042faf41dSJed Brown 321e27a552bSJed Brown /*@C 322bcf0153eSBarry Smith TSRosWRegisterAll - Registers all of the Rosenbrock-W methods in `TSROSW` 323e27a552bSJed Brown 3241d27aa22SBarry Smith Not Collective, but should be called by all MPI processes which will need the schemes to be registered 325e27a552bSJed Brown 326e27a552bSJed Brown Level: advanced 327e27a552bSJed Brown 3281cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `TSRosWRegisterDestroy()` 329e27a552bSJed Brown @*/ 330d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWRegisterAll(void) 331d71ae5a4SJacob Faibussowitsch { 332e27a552bSJed Brown PetscFunctionBegin; 3333ba16761SJacob Faibussowitsch if (TSRosWRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS); 334e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_TRUE; 335e27a552bSJed Brown 336e27a552bSJed Brown { 337bbd56ea5SKarl Rupp const PetscReal A = 0; 338bbd56ea5SKarl Rupp const PetscReal Gamma = 1; 339bbd56ea5SKarl Rupp const PetscReal b = 1; 340bbd56ea5SKarl Rupp const PetscReal binterpt = 1; 3411f80e275SEmil Constantinescu 3429566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWTHETA1, 1, 1, &A, &Gamma, &b, NULL, 1, &binterpt)); 3433606a31eSEmil Constantinescu } 3443606a31eSEmil Constantinescu 3453606a31eSEmil Constantinescu { 346bbd56ea5SKarl Rupp const PetscReal A = 0; 347bbd56ea5SKarl Rupp const PetscReal Gamma = 0.5; 348bbd56ea5SKarl Rupp const PetscReal b = 1; 349bbd56ea5SKarl Rupp const PetscReal binterpt = 1; 350bbd56ea5SKarl Rupp 3519566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWTHETA2, 2, 1, &A, &Gamma, &b, NULL, 1, &binterpt)); 3523606a31eSEmil Constantinescu } 3533606a31eSEmil Constantinescu 3543606a31eSEmil Constantinescu { 355da80777bSKarl 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. */ 356b16ce868SStefano Zampini const PetscReal A[2][2] = { 3579371c9d4SSatish Balay {0, 0}, 3589371c9d4SSatish Balay {1., 0} 359b16ce868SStefano Zampini }; 360b16ce868SStefano Zampini const PetscReal Gamma[2][2] = { 361b16ce868SStefano Zampini {1.707106781186547524401, 0 }, 362b16ce868SStefano Zampini {-2. * 1.707106781186547524401, 1.707106781186547524401} 363b16ce868SStefano Zampini }; 364b16ce868SStefano Zampini const PetscReal b[2] = {0.5, 0.5}; 365b16ce868SStefano Zampini const PetscReal b1[2] = {1.0, 0.0}; 3661f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 367da80777bSKarl Rupp binterpt[0][0] = 1.707106781186547524401 - 1.0; 368da80777bSKarl Rupp binterpt[1][0] = 2.0 - 1.707106781186547524401; 369da80777bSKarl Rupp binterpt[0][1] = 1.707106781186547524401 - 1.5; 370da80777bSKarl Rupp binterpt[1][1] = 1.5 - 1.707106781186547524401; 371bbd56ea5SKarl Rupp 3729566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSW2P, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0])); 373e27a552bSJed Brown } 374e27a552bSJed Brown { 375da80777bSKarl 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. */ 376b16ce868SStefano Zampini const PetscReal A[2][2] = { 3779371c9d4SSatish Balay {0, 0}, 3789371c9d4SSatish Balay {1., 0} 379b16ce868SStefano Zampini }; 380b16ce868SStefano Zampini const PetscReal Gamma[2][2] = { 381b16ce868SStefano Zampini {0.2928932188134524755992, 0 }, 382b16ce868SStefano Zampini {-2. * 0.2928932188134524755992, 0.2928932188134524755992} 383b16ce868SStefano Zampini }; 384b16ce868SStefano Zampini const PetscReal b[2] = {0.5, 0.5}; 385b16ce868SStefano Zampini const PetscReal b1[2] = {1.0, 0.0}; 3861f80e275SEmil Constantinescu PetscReal binterpt[2][2]; 387da80777bSKarl Rupp binterpt[0][0] = 0.2928932188134524755992 - 1.0; 388da80777bSKarl Rupp binterpt[1][0] = 2.0 - 0.2928932188134524755992; 389da80777bSKarl Rupp binterpt[0][1] = 0.2928932188134524755992 - 1.5; 390da80777bSKarl Rupp binterpt[1][1] = 1.5 - 0.2928932188134524755992; 391bbd56ea5SKarl Rupp 3929566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSW2M, 2, 2, &A[0][0], &Gamma[0][0], b, b1, 2, &binterpt[0][0])); 393fe7e6d57SJed Brown } 394fe7e6d57SJed Brown { 395da80777bSKarl Rupp /*const PetscReal g = 7.8867513459481287e-01; Directly written in-place below */ 3961f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 397b16ce868SStefano Zampini const PetscReal A[3][3] = { 3989371c9d4SSatish Balay {0, 0, 0}, 399fe7e6d57SJed Brown {1.5773502691896257e+00, 0, 0}, 4009371c9d4SSatish Balay {0.5, 0, 0} 401b16ce868SStefano Zampini }; 402b16ce868SStefano Zampini const PetscReal Gamma[3][3] = { 403b16ce868SStefano Zampini {7.8867513459481287e-01, 0, 0 }, 404b16ce868SStefano Zampini {-1.5773502691896257e+00, 7.8867513459481287e-01, 0 }, 405b16ce868SStefano Zampini {-6.7075317547305480e-01, -1.7075317547305482e-01, 7.8867513459481287e-01} 406b16ce868SStefano Zampini }; 407b16ce868SStefano Zampini const PetscReal b[3] = {1.0566243270259355e-01, 4.9038105676657971e-02, 8.4529946162074843e-01}; 408b16ce868SStefano Zampini const PetscReal b2[3] = {-1.7863279495408180e-01, 1. / 3., 8.4529946162074843e-01}; 4091f80e275SEmil Constantinescu 4101f80e275SEmil Constantinescu binterpt[0][0] = -0.8094010767585034; 4111f80e275SEmil Constantinescu binterpt[1][0] = -0.5; 4121f80e275SEmil Constantinescu binterpt[2][0] = 2.3094010767585034; 4131f80e275SEmil Constantinescu binterpt[0][1] = 0.9641016151377548; 4141f80e275SEmil Constantinescu binterpt[1][1] = 0.5; 4151f80e275SEmil Constantinescu binterpt[2][1] = -1.4641016151377548; 416bbd56ea5SKarl Rupp 4179566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWRA3PW, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0])); 418fe7e6d57SJed Brown } 419fe7e6d57SJed Brown { 4203ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 421da80777bSKarl Rupp /*const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 422b16ce868SStefano Zampini const PetscReal A[4][4] = { 4239371c9d4SSatish Balay {0, 0, 0, 0}, 424fe7e6d57SJed Brown {8.7173304301691801e-01, 0, 0, 0}, 425fe7e6d57SJed Brown {8.4457060015369423e-01, -1.1299064236484185e-01, 0, 0}, 4269371c9d4SSatish Balay {0, 0, 1., 0} 427b16ce868SStefano Zampini }; 428b16ce868SStefano Zampini const PetscReal Gamma[4][4] = { 429b16ce868SStefano Zampini {4.3586652150845900e-01, 0, 0, 0 }, 430b16ce868SStefano Zampini {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0 }, 431b16ce868SStefano Zampini {-9.0338057013044082e-01, 5.4180672388095326e-02, 4.3586652150845900e-01, 0 }, 432b16ce868SStefano Zampini {2.4212380706095346e-01, -1.2232505839045147e+00, 5.4526025533510214e-01, 4.3586652150845900e-01} 433b16ce868SStefano Zampini }; 434b16ce868SStefano Zampini const PetscReal b[4] = {2.4212380706095346e-01, -1.2232505839045147e+00, 1.5452602553351020e+00, 4.3586652150845900e-01}; 435b16ce868SStefano Zampini const PetscReal b2[4] = {3.7810903145819369e-01, -9.6042292212423178e-02, 5.0000000000000000e-01, 2.1793326075422950e-01}; 4363ca35412SEmil Constantinescu 4373ca35412SEmil Constantinescu binterpt[0][0] = 1.0564298455794094; 4383ca35412SEmil Constantinescu binterpt[1][0] = 2.296429974281067; 4393ca35412SEmil Constantinescu binterpt[2][0] = -1.307599564525376; 4403ca35412SEmil Constantinescu binterpt[3][0] = -1.045260255335102; 4413ca35412SEmil Constantinescu binterpt[0][1] = -1.3864882699759573; 4423ca35412SEmil Constantinescu binterpt[1][1] = -8.262611700275677; 4433ca35412SEmil Constantinescu binterpt[2][1] = 7.250979895056055; 4443ca35412SEmil Constantinescu binterpt[3][1] = 2.398120075195581; 4453ca35412SEmil Constantinescu binterpt[0][2] = 0.5721822314575016; 4463ca35412SEmil Constantinescu binterpt[1][2] = 4.742931142090097; 4473ca35412SEmil Constantinescu binterpt[2][2] = -4.398120075195578; 4483ca35412SEmil Constantinescu binterpt[3][2] = -0.9169932983520199; 4493ca35412SEmil Constantinescu 4509566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWRA34PW2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 451e27a552bSJed Brown } 452ef3c5b88SJed Brown { 453337ef93bSJed Brown /* const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 454337ef93bSJed Brown const PetscReal A[4][4] = { 455337ef93bSJed Brown {0, 0, 0, 0}, 456337ef93bSJed Brown {8.7173304301691801e-01, 0, 0, 0}, 457337ef93bSJed Brown {1.4722022879435914e+00, -3.1840250568090289e-01, 0, 0}, 458337ef93bSJed Brown {8.1505192016694938e-01, 5.0000000000000000e-01, -3.1505192016694938e-01, 0} 459337ef93bSJed Brown }; 460337ef93bSJed Brown const PetscReal Gamma[4][4] = { 461337ef93bSJed Brown {4.3586652150845900e-01, 0, 0, 0 }, 462337ef93bSJed Brown {-8.7173304301691801e-01, 4.3586652150845900e-01, 0, 0 }, 463337ef93bSJed Brown {-1.2855347382089872e+00, 5.0507005541550687e-01, 4.3586652150845900e-01, 0 }, 464337ef93bSJed Brown {-4.8201449182864348e-01, 2.1793326075422950e-01, -1.7178529043404503e-01, 4.3586652150845900e-01} 465337ef93bSJed Brown }; 466337ef93bSJed Brown const PetscReal b[4] = {3.3303742833830591e-01, 7.1793326075422947e-01, -4.8683721060099439e-01, 4.3586652150845900e-01}; 467337ef93bSJed Brown const PetscReal b2[4] = {2.5000000000000000e-01, 7.4276119608319180e-01, -3.1472922970066219e-01, 3.2196803361747034e-01}; 468337ef93bSJed Brown 469337ef93bSJed Brown PetscCall(TSRosWRegister(TSROSWR34PRW, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 470337ef93bSJed Brown } 471337ef93bSJed Brown { 472337ef93bSJed Brown /* const PetscReal g = 4.3586652150845900e-01; Directly written in-place below */ 473337ef93bSJed Brown const PetscReal A[4][4] = { 474337ef93bSJed Brown {0, 0, 0, 0}, 475337ef93bSJed Brown {1.3075995645253771e+00, 0, 0, 0}, 476337ef93bSJed Brown {0.5, 0.5, 0, 0}, 477337ef93bSJed Brown {0.5, 0.5, 0, 0} 478337ef93bSJed Brown }; 479337ef93bSJed Brown const PetscReal Gamma[4][4] = { 480337ef93bSJed Brown {4.3586652150845900e-01, 0, 0, 0 }, 481337ef93bSJed Brown {-1.3075995645253771e+00, 4.3586652150845900e-01, 0, 0 }, 482337ef93bSJed Brown {-7.0988575860972170e-01, -5.5996735960277766e-01, 4.3586652150845900e-01, 0 }, 483337ef93bSJed Brown {-1.5550856807552085e-01, -9.5388516575112225e-01, 6.7352721231818413e-01, 4.3586652150845900e-01} 484337ef93bSJed Brown }; 485337ef93bSJed Brown const PetscReal b[4] = {3.4449143192447917e-01, -4.5388516575112231e-01, 6.7352721231818413e-01, 4.3586652150845900e-01}; 486337ef93bSJed Brown const PetscReal b2[4] = {5.0000000000000000e-01, -2.5738812086522078e-01, 4.3542008724775044e-01, 3.2196803361747034e-01}; 487337ef93bSJed Brown 488337ef93bSJed Brown PetscCall(TSRosWRegister(TSROSWR3PRL2, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 489337ef93bSJed Brown } 490337ef93bSJed Brown { 491da80777bSKarl Rupp /* const PetscReal g = 0.5; Directly written in-place below */ 492b16ce868SStefano Zampini const PetscReal A[4][4] = { 4939371c9d4SSatish Balay {0, 0, 0, 0}, 494ef3c5b88SJed Brown {0, 0, 0, 0}, 495ef3c5b88SJed Brown {1., 0, 0, 0}, 4969371c9d4SSatish Balay {0.75, -0.25, 0.5, 0} 497b16ce868SStefano Zampini }; 498b16ce868SStefano Zampini const PetscReal Gamma[4][4] = { 499b16ce868SStefano Zampini {0.5, 0, 0, 0 }, 500b16ce868SStefano Zampini {1., 0.5, 0, 0 }, 501b16ce868SStefano Zampini {-0.25, -0.25, 0.5, 0 }, 502b16ce868SStefano Zampini {1. / 12, 1. / 12, -2. / 3, 0.5} 503b16ce868SStefano Zampini }; 504b16ce868SStefano Zampini const PetscReal b[4] = {5. / 6, -1. / 6, -1. / 6, 0.5}; 505b16ce868SStefano Zampini const PetscReal b2[4] = {0.75, -0.25, 0.5, 0}; 506bbd56ea5SKarl Rupp 5079566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWRODAS3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 508ef3c5b88SJed Brown } 509ef3c5b88SJed Brown { 51048665b53SJed Brown /* const PetscReal g = 0.25; Directly written in-place below */ 51148665b53SJed Brown const PetscReal A[6][6] = { 51248665b53SJed Brown {0, 0, 0, 0, 0, 0}, 51348665b53SJed Brown {0.75, 0, 0, 0, 0, 0}, 51448665b53SJed Brown {7.5162877593868457e-02, 2.4837122406131545e-02, 0, 0, 0, 0}, 51548665b53SJed Brown {1.6532708886396510e+00, 2.1545706385445562e-01, -1.3157488872766792e+00, 0, 0, 0}, 51648665b53SJed Brown {1.9385003738039885e+01, 1.2007117225835324e+00, -1.9337924059522791e+01, -2.4779140110062559e-01, 0, 0}, 51748665b53SJed Brown {-7.3844531665375115e+00, -3.0593419030174646e-01, 7.8622074209377981e+00, 5.7817993590145966e-01, 0.25, 0} 51848665b53SJed Brown }; 51948665b53SJed Brown const PetscReal Gamma[6][6] = { 52048665b53SJed Brown {0.25, 0, 0, 0, 0, 0 }, 52148665b53SJed Brown {-7.5000000000000000e-01, 0.25, 0, 0, 0, 0 }, 52248665b53SJed Brown {-8.8644359075349941e-02, -2.8688974257983398e-02, 0.25, 0, 0, 0 }, 52348665b53SJed Brown {-4.8470034585330284e+00, -3.1583244269672095e-01, 4.9536568360123221e+00, 0.25, 0, 0 }, 52448665b53SJed Brown {-2.6769456904577400e+01, -1.5066459128852787e+00, 2.7200131480460591e+01, 8.2597133700208525e-01, 0.25, 0 }, 52548665b53SJed Brown {6.5876206496361416e+00, 3.6807059172993878e-01, -6.7423520694658121e+00, -1.0619631475741095e-01, -3.5714285714285715e-01, 0.25} 52648665b53SJed Brown }; 52748665b53SJed Brown const PetscReal b[6] = {-7.9683251690137014e-01, 6.2136401428192344e-02, 1.1198553514719862e+00, 4.7198362114404874e-01, -1.0714285714285714e-01, 0.25}; 52848665b53SJed Brown const PetscReal b2[6] = {-7.3844531665375115e+00, -3.0593419030174646e-01, 7.8622074209377981e+00, 5.7817993590145966e-01, 0.25, 0.0}; 52948665b53SJed Brown 53048665b53SJed Brown PetscCall(TSRosWRegister(TSROSWRODASPR, 4, 6, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 53148665b53SJed Brown } 53248665b53SJed Brown { 533337ef93bSJed Brown /* const PetscReal g = 0.3125; Directly written in-place below */ 534337ef93bSJed Brown const PetscReal A[6][6] = { 535337ef93bSJed Brown {0, 0, 0, 0, 0, 0}, 536337ef93bSJed Brown {9.3750000000000000e-01, 0, 0, 0, 0, 0}, 537337ef93bSJed Brown {-4.7145892646261345e-02, 5.4531286650471122e-01, 0, 0, 0, 0}, 538337ef93bSJed Brown {4.6915543899742240e-01, 4.4490537602383673e-01, -2.2498239334061121e-01, 0, 0, 0}, 539337ef93bSJed Brown {1.0950372887345903e+00, 6.3223023457294381e-01, -8.9232966090485821e-01, 1.6506213759732410e-01, 0, 0}, 540337ef93bSJed Brown {-1.7746585073632790e-01, -5.8241418952602364e-01, 6.8180612588238165e-01, 7.6557391437996980e-01, 3.1250000000000000e-01, 0} 541337ef93bSJed Brown }; 542337ef93bSJed Brown const PetscReal Gamma[6][6] = { 543337ef93bSJed Brown {0.3125, 0, 0, 0, 0, 0 }, 544337ef93bSJed Brown {-9.3750000000000000e-01, 0.3125, 0, 0, 0, 0 }, 545337ef93bSJed Brown {-9.7580572085994507e-02, -5.8666328499964138e-01, 0.3125, 0, 0, 0 }, 546337ef93bSJed Brown {-4.9407065013256957e-01, -5.6819726428975503e-01, 5.0318949274167679e-01, 0.3125, 0, 0 }, 547337ef93bSJed Brown {-1.2725031394709183e+00, -1.2146444240989676e+00, 1.5741357867872399e+00, 6.0051177678264578e-01, 0.3125, 0 }, 548337ef93bSJed Brown {6.9690744901421153e-01, 6.2237005730756434e-01, -1.1553701989197045e+00, 1.8350029013386296e-01, -6.5990759753593431e-01, 0.3125} 549337ef93bSJed Brown }; 550337ef93bSJed Brown const PetscReal b[6] = {5.1944159827788361e-01, 3.9955867781540699e-02, -4.7356407303732290e-01, 9.4907420451383284e-01, -3.4740759753593431e-01, 0.3125}; 551337ef93bSJed Brown const PetscReal b2[6] = {-1.7746585073632790e-01, -5.8241418952602364e-01, 6.8180612588238165e-01, 7.6557391437996980e-01, 0.3125, 0}; 552337ef93bSJed Brown 553337ef93bSJed Brown PetscCall(TSRosWRegister(TSROSWRODASPR2, 4, 6, &A[0][0], &Gamma[0][0], b, b2, 0, NULL)); 554337ef93bSJed Brown } 555337ef93bSJed Brown { 556da80777bSKarl Rupp /*const PetscReal g = 0.43586652150845899941601945119356; Directly written in-place below */ 557b16ce868SStefano Zampini const PetscReal A[3][3] = { 5589371c9d4SSatish Balay {0, 0, 0}, 559da80777bSKarl Rupp {0.43586652150845899941601945119356, 0, 0}, 5609371c9d4SSatish Balay {0.43586652150845899941601945119356, 0, 0} 561b16ce868SStefano Zampini }; 562b16ce868SStefano Zampini const PetscReal Gamma[3][3] = { 563b16ce868SStefano Zampini {0.43586652150845899941601945119356, 0, 0 }, 564b16ce868SStefano Zampini {-0.19294655696029095575009695436041, 0.43586652150845899941601945119356, 0 }, 565b16ce868SStefano Zampini {0, 1.74927148125794685173529749738960, 0.43586652150845899941601945119356} 566b16ce868SStefano Zampini }; 567b16ce868SStefano Zampini const PetscReal b[3] = {-0.75457412385404315829818998646589, 1.94100407061964420292840123379419, -0.18642994676560104463021124732829}; 568b16ce868SStefano Zampini const PetscReal b2[3] = {-1.53358745784149585370766523913002, 2.81745131148625772213931745457622, -0.28386385364476186843165221544619}; 5691f80e275SEmil Constantinescu 5701f80e275SEmil Constantinescu PetscReal binterpt[3][2]; 5711f80e275SEmil Constantinescu binterpt[0][0] = 3.793692883777660870425141387941; 5721f80e275SEmil Constantinescu binterpt[1][0] = -2.918692883777660870425141387941; 5731f80e275SEmil Constantinescu binterpt[2][0] = 0.125; 5741f80e275SEmil Constantinescu binterpt[0][1] = -0.725741064379812106687651020584; 5751f80e275SEmil Constantinescu binterpt[1][1] = 0.559074397713145440020984353917; 5761f80e275SEmil Constantinescu binterpt[2][1] = 0.16666666666666666666666666666667; 5771f80e275SEmil Constantinescu 5789566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWSANDU3, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0])); 579ef3c5b88SJed Brown } 580b1c69cc3SEmil Constantinescu { 581da80777bSKarl Rupp /*const PetscReal s3 = PetscSqrtReal(3.),g = (3.0+s3)/6.0; 582da80777bSKarl Rupp * Direct evaluation: s3 = 1.732050807568877293527; 583da80777bSKarl Rupp * g = 0.7886751345948128822546; 584da80777bSKarl Rupp * Values are directly inserted below to ensure availability at compile time (compiler warnings otherwise...) */ 585b16ce868SStefano Zampini const PetscReal A[3][3] = { 5869371c9d4SSatish Balay {0, 0, 0}, 587b1c69cc3SEmil Constantinescu {1, 0, 0}, 5889371c9d4SSatish Balay {0.25, 0.25, 0} 589b16ce868SStefano Zampini }; 590b16ce868SStefano Zampini const PetscReal Gamma[3][3] = { 591b16ce868SStefano Zampini {0, 0, 0 }, 592b16ce868SStefano Zampini {(-3.0 - 1.732050807568877293527) / 6.0, 0.7886751345948128822546, 0 }, 593b16ce868SStefano Zampini {(-3.0 - 1.732050807568877293527) / 24.0, (-3.0 - 1.732050807568877293527) / 8.0, 0.7886751345948128822546} 594b16ce868SStefano Zampini }; 595b16ce868SStefano Zampini const PetscReal b[3] = {1. / 6., 1. / 6., 2. / 3.}; 596b16ce868SStefano Zampini const PetscReal b2[3] = {1. / 4., 1. / 4., 1. / 2.}; 597c0cb691aSEmil Constantinescu PetscReal binterpt[3][2]; 598da80777bSKarl Rupp 599c0cb691aSEmil Constantinescu binterpt[0][0] = 0.089316397477040902157517886164709; 600c0cb691aSEmil Constantinescu binterpt[1][0] = -0.91068360252295909784248211383529; 601c0cb691aSEmil Constantinescu binterpt[2][0] = 1.8213672050459181956849642276706; 602c0cb691aSEmil Constantinescu binterpt[0][1] = 0.077350269189625764509148780501957; 603c0cb691aSEmil Constantinescu binterpt[1][1] = 1.077350269189625764509148780502; 604c0cb691aSEmil Constantinescu binterpt[2][1] = -1.1547005383792515290182975610039; 605bbd56ea5SKarl Rupp 6069566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWASSP3P3S1C, 3, 3, &A[0][0], &Gamma[0][0], b, b2, 2, &binterpt[0][0])); 607b1c69cc3SEmil Constantinescu } 608b1c69cc3SEmil Constantinescu 609b1c69cc3SEmil Constantinescu { 610b16ce868SStefano Zampini const PetscReal A[4][4] = { 6119371c9d4SSatish Balay {0, 0, 0, 0}, 612b1c69cc3SEmil Constantinescu {1. / 2., 0, 0, 0}, 613b1c69cc3SEmil Constantinescu {1. / 2., 1. / 2., 0, 0}, 6149371c9d4SSatish Balay {1. / 6., 1. / 6., 1. / 6., 0} 615b16ce868SStefano Zampini }; 616b16ce868SStefano Zampini const PetscReal Gamma[4][4] = { 617b16ce868SStefano Zampini {1. / 2., 0, 0, 0}, 618b16ce868SStefano Zampini {0.0, 1. / 4., 0, 0}, 619b16ce868SStefano Zampini {-2., -2. / 3., 2. / 3., 0}, 620b16ce868SStefano Zampini {1. / 2., 5. / 36., -2. / 9, 0} 621b16ce868SStefano Zampini }; 622b16ce868SStefano Zampini const PetscReal b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}; 623b16ce868SStefano Zampini const PetscReal b2[4] = {1. / 8., 3. / 4., 1. / 8., 0}; 624c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 625da80777bSKarl Rupp 626c0cb691aSEmil Constantinescu binterpt[0][0] = 6.25; 627c0cb691aSEmil Constantinescu binterpt[1][0] = -30.25; 628c0cb691aSEmil Constantinescu binterpt[2][0] = 1.75; 629c0cb691aSEmil Constantinescu binterpt[3][0] = 23.25; 630c0cb691aSEmil Constantinescu binterpt[0][1] = -9.75; 631c0cb691aSEmil Constantinescu binterpt[1][1] = 58.75; 632c0cb691aSEmil Constantinescu binterpt[2][1] = -3.25; 633c0cb691aSEmil Constantinescu binterpt[3][1] = -45.75; 634c0cb691aSEmil Constantinescu binterpt[0][2] = 3.6666666666666666666666666666667; 635c0cb691aSEmil Constantinescu binterpt[1][2] = -28.333333333333333333333333333333; 636c0cb691aSEmil Constantinescu binterpt[2][2] = 1.6666666666666666666666666666667; 637c0cb691aSEmil Constantinescu binterpt[3][2] = 23.; 638bbd56ea5SKarl Rupp 6399566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWLASSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 640b1c69cc3SEmil Constantinescu } 641b1c69cc3SEmil Constantinescu 642b1c69cc3SEmil Constantinescu { 643b16ce868SStefano Zampini const PetscReal A[4][4] = { 6449371c9d4SSatish Balay {0, 0, 0, 0}, 645b1c69cc3SEmil Constantinescu {1. / 2., 0, 0, 0}, 646b1c69cc3SEmil Constantinescu {1. / 2., 1. / 2., 0, 0}, 6479371c9d4SSatish Balay {1. / 6., 1. / 6., 1. / 6., 0} 648b16ce868SStefano Zampini }; 649b16ce868SStefano Zampini const PetscReal Gamma[4][4] = { 650b16ce868SStefano Zampini {1. / 2., 0, 0, 0}, 651b16ce868SStefano Zampini {0.0, 3. / 4., 0, 0}, 652b16ce868SStefano Zampini {-2. / 3., -23. / 9., 2. / 9., 0}, 653b16ce868SStefano Zampini {1. / 18., 65. / 108., -2. / 27, 0} 654b16ce868SStefano Zampini }; 655b16ce868SStefano Zampini const PetscReal b[4] = {1. / 6., 1. / 6., 1. / 6., 1. / 2.}; 656b16ce868SStefano Zampini const PetscReal b2[4] = {3. / 16., 10. / 16., 3. / 16., 0}; 657c0cb691aSEmil Constantinescu PetscReal binterpt[4][3]; 658da80777bSKarl Rupp 659c0cb691aSEmil Constantinescu binterpt[0][0] = 1.6911764705882352941176470588235; 660c0cb691aSEmil Constantinescu binterpt[1][0] = 3.6813725490196078431372549019608; 661c0cb691aSEmil Constantinescu binterpt[2][0] = 0.23039215686274509803921568627451; 662c0cb691aSEmil Constantinescu binterpt[3][0] = -4.6029411764705882352941176470588; 663c0cb691aSEmil Constantinescu binterpt[0][1] = -0.95588235294117647058823529411765; 664c0cb691aSEmil Constantinescu binterpt[1][1] = -6.2401960784313725490196078431373; 665c0cb691aSEmil Constantinescu binterpt[2][1] = -0.31862745098039215686274509803922; 666c0cb691aSEmil Constantinescu binterpt[3][1] = 7.5147058823529411764705882352941; 667c0cb691aSEmil Constantinescu binterpt[0][2] = -0.56862745098039215686274509803922; 668c0cb691aSEmil Constantinescu binterpt[1][2] = 2.7254901960784313725490196078431; 669c0cb691aSEmil Constantinescu binterpt[2][2] = 0.25490196078431372549019607843137; 670c0cb691aSEmil Constantinescu binterpt[3][2] = -2.4117647058823529411764705882353; 671bbd56ea5SKarl Rupp 6729566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWLLSSP3P4S2C, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 673b1c69cc3SEmil Constantinescu } 674753f8adbSEmil Constantinescu 675753f8adbSEmil Constantinescu { 676753f8adbSEmil Constantinescu PetscReal A[4][4], Gamma[4][4], b[4], b2[4]; 6773ca35412SEmil Constantinescu PetscReal binterpt[4][3]; 678753f8adbSEmil Constantinescu 679753f8adbSEmil Constantinescu Gamma[0][0] = 0.4358665215084589994160194475295062513822671686978816; 6809371c9d4SSatish Balay Gamma[0][1] = 0; 6819371c9d4SSatish Balay Gamma[0][2] = 0; 6829371c9d4SSatish Balay Gamma[0][3] = 0; 683753f8adbSEmil Constantinescu Gamma[1][0] = -1.997527830934941248426324674704153457289527280554476; 684753f8adbSEmil Constantinescu Gamma[1][1] = 0.4358665215084589994160194475295062513822671686978816; 6859371c9d4SSatish Balay Gamma[1][2] = 0; 6869371c9d4SSatish Balay Gamma[1][3] = 0; 687753f8adbSEmil Constantinescu Gamma[2][0] = -1.007948511795029620852002345345404191008352770119903; 688753f8adbSEmil Constantinescu Gamma[2][1] = -0.004648958462629345562774289390054679806993396798458131; 689753f8adbSEmil Constantinescu Gamma[2][2] = 0.4358665215084589994160194475295062513822671686978816; 69005e8e825SJed Brown Gamma[2][3] = 0; 691753f8adbSEmil Constantinescu Gamma[3][0] = -0.6685429734233467180451604600279552604364311322650783; 692753f8adbSEmil Constantinescu Gamma[3][1] = 0.6056625986449338476089525334450053439525178740492984; 693753f8adbSEmil Constantinescu Gamma[3][2] = -0.9717899277217721234705114616271378792182450260943198; 694753f8adbSEmil Constantinescu Gamma[3][3] = 0; 695753f8adbSEmil Constantinescu 6969371c9d4SSatish Balay A[0][0] = 0; 6979371c9d4SSatish Balay A[0][1] = 0; 6989371c9d4SSatish Balay A[0][2] = 0; 6999371c9d4SSatish Balay A[0][3] = 0; 700753f8adbSEmil Constantinescu A[1][0] = 0.8717330430169179988320388950590125027645343373957631; 7019371c9d4SSatish Balay A[1][1] = 0; 7029371c9d4SSatish Balay A[1][2] = 0; 7039371c9d4SSatish Balay A[1][3] = 0; 704753f8adbSEmil Constantinescu A[2][0] = 0.5275890119763004115618079766722914408876108660811028; 705753f8adbSEmil Constantinescu A[2][1] = 0.07241098802369958843819203208518599088698057726988732; 7069371c9d4SSatish Balay A[2][2] = 0; 7079371c9d4SSatish Balay A[2][3] = 0; 708753f8adbSEmil Constantinescu A[3][0] = 0.3990960076760701320627260685975778145384666450351314; 709753f8adbSEmil Constantinescu A[3][1] = -0.4375576546135194437228463747348862825846903771419953; 710753f8adbSEmil Constantinescu A[3][2] = 1.038461646937449311660120300601880176655352737312713; 71105e8e825SJed Brown A[3][3] = 0; 712753f8adbSEmil Constantinescu 713753f8adbSEmil Constantinescu b[0] = 0.1876410243467238251612921333138006734899663569186926; 714753f8adbSEmil Constantinescu b[1] = -0.5952974735769549480478230473706443582188442040780541; 715753f8adbSEmil Constantinescu b[2] = 0.9717899277217721234705114616271378792182450260943198; 716753f8adbSEmil Constantinescu b[3] = 0.4358665215084589994160194475295062513822671686978816; 717753f8adbSEmil Constantinescu 718753f8adbSEmil Constantinescu b2[0] = 0.2147402862233891404862383521089097657790734483804460; 719753f8adbSEmil Constantinescu b2[1] = -0.4851622638849390928209050538171743017757490232519684; 720753f8adbSEmil Constantinescu b2[2] = 0.8687250025203875511662123688667549217531982787600080; 721753f8adbSEmil Constantinescu b2[3] = 0.4016969751411624011684543450940068201770721128357014; 722753f8adbSEmil Constantinescu 7233ca35412SEmil Constantinescu binterpt[0][0] = 2.2565812720167954547104627844105; 7243ca35412SEmil Constantinescu binterpt[1][0] = 1.349166413351089573796243820819; 7253ca35412SEmil Constantinescu binterpt[2][0] = -2.4695174540533503758652847586647; 7263ca35412SEmil Constantinescu binterpt[3][0] = -0.13623023131453465264142184656474; 7273ca35412SEmil Constantinescu binterpt[0][1] = -3.0826699111559187902922463354557; 7283ca35412SEmil Constantinescu binterpt[1][1] = -2.4689115685996042534544925650515; 7293ca35412SEmil Constantinescu binterpt[2][1] = 5.7428279814696677152129332773553; 7303ca35412SEmil Constantinescu binterpt[3][1] = -0.19124650171414467146619437684812; 7313ca35412SEmil Constantinescu binterpt[0][2] = 1.0137296634858471607430756831148; 7323ca35412SEmil Constantinescu binterpt[1][2] = 0.52444768167155973161042570784064; 7333ca35412SEmil Constantinescu binterpt[2][2] = -2.3015205996945452158771370439586; 7343ca35412SEmil Constantinescu binterpt[3][2] = 0.76334325453713832352363565300308; 735f4aed992SEmil Constantinescu 7369566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(TSROSWARK3, 3, 4, &A[0][0], &Gamma[0][0], b, b2, 3, &binterpt[0][0])); 737753f8adbSEmil Constantinescu } 73809cb0f53SBarry Smith PetscCall(TSRosWRegisterRos4(TSROSWGRK4T, 0.231, PETSC_CURRENT, PETSC_CURRENT, 0, -0.1282612945269037e+01)); 73909cb0f53SBarry Smith PetscCall(TSRosWRegisterRos4(TSROSWSHAMP4, 0.5, PETSC_CURRENT, PETSC_CURRENT, 0, 125. / 108.)); 74009cb0f53SBarry Smith PetscCall(TSRosWRegisterRos4(TSROSWVELDD4, 0.22570811482256823492, PETSC_CURRENT, PETSC_CURRENT, 0, -1.355958941201148)); 74109cb0f53SBarry Smith PetscCall(TSRosWRegisterRos4(TSROSW4L, 0.57282, PETSC_CURRENT, PETSC_CURRENT, 0, -1.093502252409163)); 7423ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 743e27a552bSJed Brown } 744e27a552bSJed Brown 745e27a552bSJed Brown /*@C 746bcf0153eSBarry Smith TSRosWRegisterDestroy - Frees the list of schemes that were registered by `TSRosWRegister()`. 747e27a552bSJed Brown 748e27a552bSJed Brown Not Collective 749e27a552bSJed Brown 750e27a552bSJed Brown Level: advanced 751e27a552bSJed Brown 7521cc06b55SBarry Smith .seealso: [](ch_ts), `TSRosWRegister()`, `TSRosWRegisterAll()` 753e27a552bSJed Brown @*/ 754d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWRegisterDestroy(void) 755d71ae5a4SJacob Faibussowitsch { 75661692a83SJed Brown RosWTableauLink link; 757e27a552bSJed Brown 758e27a552bSJed Brown PetscFunctionBegin; 75961692a83SJed Brown while ((link = RosWTableauList)) { 76061692a83SJed Brown RosWTableau t = &link->tab; 76161692a83SJed Brown RosWTableauList = link->next; 7629566063dSJacob Faibussowitsch PetscCall(PetscFree5(t->A, t->Gamma, t->b, t->ASum, t->GammaSum)); 7639566063dSJacob Faibussowitsch PetscCall(PetscFree5(t->At, t->bt, t->GammaInv, t->GammaZeroDiag, t->GammaExplicitCorr)); 7649566063dSJacob Faibussowitsch PetscCall(PetscFree2(t->bembed, t->bembedt)); 7659566063dSJacob Faibussowitsch PetscCall(PetscFree(t->binterpt)); 7669566063dSJacob Faibussowitsch PetscCall(PetscFree(t->name)); 7679566063dSJacob Faibussowitsch PetscCall(PetscFree(link)); 768e27a552bSJed Brown } 769e27a552bSJed Brown TSRosWRegisterAllCalled = PETSC_FALSE; 7703ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 771e27a552bSJed Brown } 772e27a552bSJed Brown 773e27a552bSJed Brown /*@C 774bcf0153eSBarry Smith TSRosWInitializePackage - This function initializes everything in the `TSROSW` package. It is called 775bcf0153eSBarry Smith from `TSInitializePackage()`. 776e27a552bSJed Brown 777e27a552bSJed Brown Level: developer 778e27a552bSJed Brown 7791cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `PetscInitialize()`, `TSRosWFinalizePackage()` 780e27a552bSJed Brown @*/ 781d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWInitializePackage(void) 782d71ae5a4SJacob Faibussowitsch { 783e27a552bSJed Brown PetscFunctionBegin; 7843ba16761SJacob Faibussowitsch if (TSRosWPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 785e27a552bSJed Brown TSRosWPackageInitialized = PETSC_TRUE; 7869566063dSJacob Faibussowitsch PetscCall(TSRosWRegisterAll()); 7879566063dSJacob Faibussowitsch PetscCall(PetscRegisterFinalize(TSRosWFinalizePackage)); 7883ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 789e27a552bSJed Brown } 790e27a552bSJed Brown 791e27a552bSJed Brown /*@C 792bcf0153eSBarry Smith TSRosWFinalizePackage - This function destroys everything in the `TSROSW` package. It is 793bcf0153eSBarry Smith called from `PetscFinalize()`. 794e27a552bSJed Brown 795e27a552bSJed Brown Level: developer 796e27a552bSJed Brown 7971cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW`, `PetscFinalize()`, `TSRosWInitializePackage()` 798e27a552bSJed Brown @*/ 799d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWFinalizePackage(void) 800d71ae5a4SJacob Faibussowitsch { 801e27a552bSJed Brown PetscFunctionBegin; 802e27a552bSJed Brown TSRosWPackageInitialized = PETSC_FALSE; 8039566063dSJacob Faibussowitsch PetscCall(TSRosWRegisterDestroy()); 8043ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 805e27a552bSJed Brown } 806e27a552bSJed Brown 807e27a552bSJed Brown /*@C 808bcf0153eSBarry Smith TSRosWRegister - register a `TSROSW`, Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 809e27a552bSJed Brown 810e27a552bSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 811e27a552bSJed Brown 812e27a552bSJed Brown Input Parameters: 813e27a552bSJed Brown + name - identifier for method 814e27a552bSJed Brown . order - approximation order of method 815e27a552bSJed Brown . s - number of stages, this is the dimension of the matrices below 81661692a83SJed Brown . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 81761692a83SJed Brown . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 818fe7e6d57SJed Brown . b - Step completion table (dimension s) 8190298fd71SBarry Smith . bembed - Step completion table for a scheme of order one less (dimension s, NULL if no embedded scheme is available) 820f4aed992SEmil Constantinescu . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt 82142faf41dSJed Brown - binterpt - Coefficients of the interpolation formula (dimension s*pinterp) 822e27a552bSJed Brown 823e27a552bSJed Brown Level: advanced 824e27a552bSJed Brown 825bcf0153eSBarry Smith Note: 826bcf0153eSBarry Smith Several Rosenbrock W methods are provided, this function is only needed to create new methods. 827bcf0153eSBarry Smith 8281cc06b55SBarry Smith .seealso: [](ch_ts), `TSROSW` 829e27a552bSJed Brown @*/ 830d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWRegister(TSRosWType name, PetscInt order, PetscInt s, const PetscReal A[], const PetscReal Gamma[], const PetscReal b[], const PetscReal bembed[], PetscInt pinterp, const PetscReal binterpt[]) 831d71ae5a4SJacob Faibussowitsch { 83261692a83SJed Brown RosWTableauLink link; 83361692a83SJed Brown RosWTableau t; 83461692a83SJed Brown PetscInt i, j, k; 83561692a83SJed Brown PetscScalar *GammaInv; 836e27a552bSJed Brown 837e27a552bSJed Brown PetscFunctionBegin; 8384f572ea9SToby Isaac PetscAssertPointer(name, 1); 8394f572ea9SToby Isaac PetscAssertPointer(A, 4); 8404f572ea9SToby Isaac PetscAssertPointer(Gamma, 5); 8414f572ea9SToby Isaac PetscAssertPointer(b, 6); 8424f572ea9SToby Isaac if (bembed) PetscAssertPointer(bembed, 7); 843fe7e6d57SJed Brown 8449566063dSJacob Faibussowitsch PetscCall(TSRosWInitializePackage()); 8459566063dSJacob Faibussowitsch PetscCall(PetscNew(&link)); 846e27a552bSJed Brown t = &link->tab; 8479566063dSJacob Faibussowitsch PetscCall(PetscStrallocpy(name, &t->name)); 848e27a552bSJed Brown t->order = order; 849e27a552bSJed Brown t->s = s; 8509566063dSJacob Faibussowitsch PetscCall(PetscMalloc5(s * s, &t->A, s * s, &t->Gamma, s, &t->b, s, &t->ASum, s, &t->GammaSum)); 8519566063dSJacob Faibussowitsch PetscCall(PetscMalloc5(s * s, &t->At, s, &t->bt, s * s, &t->GammaInv, s, &t->GammaZeroDiag, s * s, &t->GammaExplicitCorr)); 8529566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->A, A, s * s)); 8539566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->Gamma, Gamma, s * s)); 8549566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->GammaExplicitCorr, Gamma, s * s)); 8559566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->b, b, s)); 856fe7e6d57SJed Brown if (bembed) { 8579566063dSJacob Faibussowitsch PetscCall(PetscMalloc2(s, &t->bembed, s, &t->bembedt)); 8589566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->bembed, bembed, s)); 859fe7e6d57SJed Brown } 86061692a83SJed Brown for (i = 0; i < s; i++) { 86161692a83SJed Brown t->ASum[i] = 0; 86261692a83SJed Brown t->GammaSum[i] = 0; 86361692a83SJed Brown for (j = 0; j < s; j++) { 86461692a83SJed Brown t->ASum[i] += A[i * s + j]; 865fe7e6d57SJed Brown t->GammaSum[i] += Gamma[i * s + j]; 86661692a83SJed Brown } 86761692a83SJed Brown } 8689566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(s * s, &GammaInv)); /* Need to use Scalar for inverse, then convert back to Real */ 86961692a83SJed Brown for (i = 0; i < s * s; i++) GammaInv[i] = Gamma[i]; 870fd96d5b0SEmil Constantinescu for (i = 0; i < s; i++) { 871fd96d5b0SEmil Constantinescu if (Gamma[i * s + i] == 0.0) { 872fd96d5b0SEmil Constantinescu GammaInv[i * s + i] = 1.0; 873c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_TRUE; 874fd96d5b0SEmil Constantinescu } else { 875c17803e7SJed Brown t->GammaZeroDiag[i] = PETSC_FALSE; 876fd96d5b0SEmil Constantinescu } 877fd96d5b0SEmil Constantinescu } 878fd96d5b0SEmil Constantinescu 87961692a83SJed Brown switch (s) { 880d71ae5a4SJacob Faibussowitsch case 1: 881d71ae5a4SJacob Faibussowitsch GammaInv[0] = 1. / GammaInv[0]; 882d71ae5a4SJacob Faibussowitsch break; 883d71ae5a4SJacob Faibussowitsch case 2: 884d71ae5a4SJacob Faibussowitsch PetscCall(PetscKernel_A_gets_inverse_A_2(GammaInv, 0, PETSC_FALSE, NULL)); 885d71ae5a4SJacob Faibussowitsch break; 886d71ae5a4SJacob Faibussowitsch case 3: 887d71ae5a4SJacob Faibussowitsch PetscCall(PetscKernel_A_gets_inverse_A_3(GammaInv, 0, PETSC_FALSE, NULL)); 888d71ae5a4SJacob Faibussowitsch break; 889d71ae5a4SJacob Faibussowitsch case 4: 890d71ae5a4SJacob Faibussowitsch PetscCall(PetscKernel_A_gets_inverse_A_4(GammaInv, 0, PETSC_FALSE, NULL)); 891d71ae5a4SJacob Faibussowitsch break; 89261692a83SJed Brown case 5: { 89361692a83SJed Brown PetscInt ipvt5[5]; 89461692a83SJed Brown MatScalar work5[5 * 5]; 8959371c9d4SSatish Balay PetscCall(PetscKernel_A_gets_inverse_A_5(GammaInv, ipvt5, work5, 0, PETSC_FALSE, NULL)); 8969371c9d4SSatish Balay break; 89761692a83SJed Brown } 898d71ae5a4SJacob Faibussowitsch case 6: 899d71ae5a4SJacob Faibussowitsch PetscCall(PetscKernel_A_gets_inverse_A_6(GammaInv, 0, PETSC_FALSE, NULL)); 900d71ae5a4SJacob Faibussowitsch break; 901d71ae5a4SJacob Faibussowitsch case 7: 902d71ae5a4SJacob Faibussowitsch PetscCall(PetscKernel_A_gets_inverse_A_7(GammaInv, 0, PETSC_FALSE, NULL)); 903d71ae5a4SJacob Faibussowitsch break; 904d71ae5a4SJacob Faibussowitsch default: 905d71ae5a4SJacob Faibussowitsch SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " stages", s); 90661692a83SJed Brown } 90761692a83SJed Brown for (i = 0; i < s * s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 9089566063dSJacob Faibussowitsch PetscCall(PetscFree(GammaInv)); 90943b21953SEmil Constantinescu 91043b21953SEmil Constantinescu for (i = 0; i < s; i++) { 91143b21953SEmil Constantinescu for (k = 0; k < i + 1; k++) { 91243b21953SEmil Constantinescu t->GammaExplicitCorr[i * s + k] = (t->GammaExplicitCorr[i * s + k]) * (t->GammaInv[k * s + k]); 913ad540459SPierre Jolivet for (j = k + 1; j < i + 1; j++) t->GammaExplicitCorr[i * s + k] += (t->GammaExplicitCorr[i * s + j]) * (t->GammaInv[j * s + k]); 91443b21953SEmil Constantinescu } 91543b21953SEmil Constantinescu } 91643b21953SEmil Constantinescu 91761692a83SJed Brown for (i = 0; i < s; i++) { 91861692a83SJed Brown for (j = 0; j < s; j++) { 91961692a83SJed Brown t->At[i * s + j] = 0; 920ad540459SPierre Jolivet for (k = 0; k < s; k++) t->At[i * s + j] += t->A[i * s + k] * t->GammaInv[k * s + j]; 92161692a83SJed Brown } 92261692a83SJed Brown t->bt[i] = 0; 923ad540459SPierre Jolivet for (j = 0; j < s; j++) t->bt[i] += t->b[j] * t->GammaInv[j * s + i]; 924fe7e6d57SJed Brown if (bembed) { 925fe7e6d57SJed Brown t->bembedt[i] = 0; 926ad540459SPierre Jolivet for (j = 0; j < s; j++) t->bembedt[i] += t->bembed[j] * t->GammaInv[j * s + i]; 927fe7e6d57SJed Brown } 92861692a83SJed Brown } 9298d59e960SJed Brown t->ccfl = 1.0; /* Fix this */ 9308d59e960SJed Brown 931f4aed992SEmil Constantinescu t->pinterp = pinterp; 9329566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(s * pinterp, &t->binterpt)); 9339566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(t->binterpt, binterpt, s * pinterp)); 93461692a83SJed Brown link->next = RosWTableauList; 93561692a83SJed Brown RosWTableauList = link; 9363ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 937e27a552bSJed Brown } 938e27a552bSJed Brown 93942faf41dSJed Brown /*@C 940fd292e60Sprj- TSRosWRegisterRos4 - register a fourth order Rosenbrock scheme by providing parameter choices 94142faf41dSJed Brown 94242faf41dSJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 94342faf41dSJed Brown 94442faf41dSJed Brown Input Parameters: 94542faf41dSJed Brown + name - identifier for method 94642faf41dSJed Brown . gamma - leading coefficient (diagonal entry) 94709cb0f53SBarry Smith . a2 - design parameter, see Table 7.2 of {cite}`wanner1996solving` 94809cb0f53SBarry Smith . a3 - design parameter or `PETSC_DETERMINE` to satisfy one of the order five conditions (Eq 7.22) 94909cb0f53SBarry Smith . b3 - design parameter, see Table 7.2 of {cite}`wanner1996solving` 950a2b725a8SWilliam Gropp - e4 - design parameter for embedded method, see coefficient E4 in ros4.f code from Hairer 95142faf41dSJed Brown 952bcf0153eSBarry Smith Level: developer 953bcf0153eSBarry Smith 95442faf41dSJed Brown Notes: 95509cb0f53SBarry Smith This routine encodes the design of fourth order Rosenbrock methods as described in {cite}`wanner1996solving` 95642faf41dSJed Brown It is used here to implement several methods from the book and can be used to experiment with new methods. 95742faf41dSJed Brown It was written this way instead of by copying coefficients in order to provide better than double precision satisfaction of the order conditions. 95842faf41dSJed Brown 959bfe80ac4SPierre Jolivet .seealso: [](ch_ts), `TSROSW`, `TSRosWRegister()` 96042faf41dSJed Brown @*/ 961d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWRegisterRos4(TSRosWType name, PetscReal gamma, PetscReal a2, PetscReal a3, PetscReal b3, PetscReal e4) 962d71ae5a4SJacob Faibussowitsch { 96342faf41dSJed Brown /* Declare numeric constants so they can be quad precision without being truncated at double */ 9649371c9d4SSatish Balay const PetscReal one = 1, two = 2, three = 3, four = 4, five = 5, six = 6, eight = 8, twelve = 12, twenty = 20, twentyfour = 24, p32 = one / six - gamma + gamma * gamma, p42 = one / eight - gamma / three, p43 = one / twelve - gamma / three, p44 = one / twentyfour - gamma / two + three / two * gamma * gamma - gamma * gamma * gamma, p56 = one / twenty - gamma / four; 96542faf41dSJed Brown PetscReal a4, a32, a42, a43, b1, b2, b4, beta2p, beta3p, beta4p, beta32, beta42, beta43, beta32beta2p, beta4jbetajp; 96642faf41dSJed Brown PetscReal A[4][4], Gamma[4][4], b[4], bm[4]; 96742faf41dSJed Brown PetscScalar M[3][3], rhs[3]; 96842faf41dSJed Brown 96942faf41dSJed Brown PetscFunctionBegin; 97042faf41dSJed Brown /* Step 1: choose Gamma (input) */ 97142faf41dSJed Brown /* Step 2: choose a2,a3,a4; b1,b2,b3,b4 to satisfy order conditions */ 97209cb0f53SBarry Smith if (a3 == (PetscReal)PETSC_DEFAULT || a3 == (PetscReal)PETSC_DETERMINE) a3 = (one / five - a2 / four) / (one / four - a2 / three); /* Eq 7.22 */ 97342faf41dSJed Brown a4 = a3; /* consequence of 7.20 */ 97442faf41dSJed Brown 97542faf41dSJed Brown /* Solve order conditions 7.15a, 7.15c, 7.15e */ 9769371c9d4SSatish Balay M[0][0] = one; 9779371c9d4SSatish Balay M[0][1] = one; 9789371c9d4SSatish Balay M[0][2] = one; /* 7.15a */ 9799371c9d4SSatish Balay M[1][0] = 0.0; 9809371c9d4SSatish Balay M[1][1] = a2 * a2; 9819371c9d4SSatish Balay M[1][2] = a4 * a4; /* 7.15c */ 9829371c9d4SSatish Balay M[2][0] = 0.0; 9839371c9d4SSatish Balay M[2][1] = a2 * a2 * a2; 9849371c9d4SSatish Balay M[2][2] = a4 * a4 * a4; /* 7.15e */ 98542faf41dSJed Brown rhs[0] = one - b3; 98642faf41dSJed Brown rhs[1] = one / three - a3 * a3 * b3; 98742faf41dSJed Brown rhs[2] = one / four - a3 * a3 * a3 * b3; 9889566063dSJacob Faibussowitsch PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL)); 98942faf41dSJed Brown b1 = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]); 99042faf41dSJed Brown b2 = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]); 99142faf41dSJed Brown b4 = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]); 99242faf41dSJed Brown 99342faf41dSJed Brown /* Step 3 */ 99442faf41dSJed Brown beta43 = (p56 - a2 * p43) / (b4 * a3 * a3 * (a3 - a2)); /* 7.21 */ 99542faf41dSJed Brown beta32beta2p = p44 / (b4 * beta43); /* 7.15h */ 99642faf41dSJed Brown beta4jbetajp = (p32 - b3 * beta32beta2p) / b4; 9979371c9d4SSatish Balay M[0][0] = b2; 9989371c9d4SSatish Balay M[0][1] = b3; 9999371c9d4SSatish Balay M[0][2] = b4; 10009371c9d4SSatish Balay M[1][0] = a4 * a4 * beta32beta2p - a3 * a3 * beta4jbetajp; 10019371c9d4SSatish Balay M[1][1] = a2 * a2 * beta4jbetajp; 10029371c9d4SSatish Balay M[1][2] = -a2 * a2 * beta32beta2p; 10039371c9d4SSatish Balay M[2][0] = b4 * beta43 * a3 * a3 - p43; 10049371c9d4SSatish Balay M[2][1] = -b4 * beta43 * a2 * a2; 10059371c9d4SSatish Balay M[2][2] = 0; 10069371c9d4SSatish Balay rhs[0] = one / two - gamma; 10079371c9d4SSatish Balay rhs[1] = 0; 10089371c9d4SSatish Balay rhs[2] = -a2 * a2 * p32; 10099566063dSJacob Faibussowitsch PetscCall(PetscKernel_A_gets_inverse_A_3(&M[0][0], 0, PETSC_FALSE, NULL)); 101042faf41dSJed Brown beta2p = PetscRealPart(M[0][0] * rhs[0] + M[0][1] * rhs[1] + M[0][2] * rhs[2]); 101142faf41dSJed Brown beta3p = PetscRealPart(M[1][0] * rhs[0] + M[1][1] * rhs[1] + M[1][2] * rhs[2]); 101242faf41dSJed Brown beta4p = PetscRealPart(M[2][0] * rhs[0] + M[2][1] * rhs[1] + M[2][2] * rhs[2]); 101342faf41dSJed Brown 101442faf41dSJed Brown /* Step 4: back-substitute */ 101542faf41dSJed Brown beta32 = beta32beta2p / beta2p; 101642faf41dSJed Brown beta42 = (beta4jbetajp - beta43 * beta3p) / beta2p; 101742faf41dSJed Brown 101842faf41dSJed Brown /* Step 5: 7.15f and 7.20, then 7.16 */ 101942faf41dSJed Brown a43 = 0; 102042faf41dSJed Brown a32 = p42 / (b3 * a3 * beta2p + b4 * a4 * beta2p); 102142faf41dSJed Brown a42 = a32; 102242faf41dSJed Brown 10239371c9d4SSatish Balay A[0][0] = 0; 10249371c9d4SSatish Balay A[0][1] = 0; 10259371c9d4SSatish Balay A[0][2] = 0; 10269371c9d4SSatish Balay A[0][3] = 0; 10279371c9d4SSatish Balay A[1][0] = a2; 10289371c9d4SSatish Balay A[1][1] = 0; 10299371c9d4SSatish Balay A[1][2] = 0; 10309371c9d4SSatish Balay A[1][3] = 0; 10319371c9d4SSatish Balay A[2][0] = a3 - a32; 10329371c9d4SSatish Balay A[2][1] = a32; 10339371c9d4SSatish Balay A[2][2] = 0; 10349371c9d4SSatish Balay A[2][3] = 0; 10359371c9d4SSatish Balay A[3][0] = a4 - a43 - a42; 10369371c9d4SSatish Balay A[3][1] = a42; 10379371c9d4SSatish Balay A[3][2] = a43; 10389371c9d4SSatish Balay A[3][3] = 0; 10399371c9d4SSatish Balay Gamma[0][0] = gamma; 10409371c9d4SSatish Balay Gamma[0][1] = 0; 10419371c9d4SSatish Balay Gamma[0][2] = 0; 10429371c9d4SSatish Balay Gamma[0][3] = 0; 10439371c9d4SSatish Balay Gamma[1][0] = beta2p - A[1][0]; 10449371c9d4SSatish Balay Gamma[1][1] = gamma; 10459371c9d4SSatish Balay Gamma[1][2] = 0; 10469371c9d4SSatish Balay Gamma[1][3] = 0; 10479371c9d4SSatish Balay Gamma[2][0] = beta3p - beta32 - A[2][0]; 10489371c9d4SSatish Balay Gamma[2][1] = beta32 - A[2][1]; 10499371c9d4SSatish Balay Gamma[2][2] = gamma; 10509371c9d4SSatish Balay Gamma[2][3] = 0; 10519371c9d4SSatish Balay Gamma[3][0] = beta4p - beta42 - beta43 - A[3][0]; 10529371c9d4SSatish Balay Gamma[3][1] = beta42 - A[3][1]; 10539371c9d4SSatish Balay Gamma[3][2] = beta43 - A[3][2]; 10549371c9d4SSatish Balay Gamma[3][3] = gamma; 10559371c9d4SSatish Balay b[0] = b1; 10569371c9d4SSatish Balay b[1] = b2; 10579371c9d4SSatish Balay b[2] = b3; 10589371c9d4SSatish Balay b[3] = b4; 105942faf41dSJed Brown 106042faf41dSJed Brown /* Construct embedded formula using given e4. We are solving Equation 7.18. */ 106142faf41dSJed Brown bm[3] = b[3] - e4 * gamma; /* using definition of E4 */ 106242faf41dSJed Brown bm[2] = (p32 - beta4jbetajp * bm[3]) / (beta32 * beta2p); /* fourth row of 7.18 */ 106342faf41dSJed Brown bm[1] = (one / two - gamma - beta3p * bm[2] - beta4p * bm[3]) / beta2p; /* second row */ 106442faf41dSJed Brown bm[0] = one - bm[1] - bm[2] - bm[3]; /* first row */ 106542faf41dSJed Brown 106642faf41dSJed Brown { 106742faf41dSJed Brown const PetscReal misfit = a2 * a2 * bm[1] + a3 * a3 * bm[2] + a4 * a4 * bm[3] - one / three; 10683c633725SBarry Smith PetscCheck(PetscAbs(misfit) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "Assumptions violated, could not construct a third order embedded method"); 106942faf41dSJed Brown } 10709566063dSJacob Faibussowitsch PetscCall(TSRosWRegister(name, 4, 4, &A[0][0], &Gamma[0][0], b, bm, 0, NULL)); 10713ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 107242faf41dSJed Brown } 107342faf41dSJed Brown 10741c3436cfSJed Brown /* 10751c3436cfSJed Brown The step completion formula is 10761c3436cfSJed Brown 10771c3436cfSJed Brown x1 = x0 + b^T Y 10781c3436cfSJed Brown 10791c3436cfSJed Brown where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 10801c3436cfSJed Brown updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 10811c3436cfSJed Brown 10821c3436cfSJed Brown x1e = x0 + be^T Y 10831c3436cfSJed Brown = x1 - b^T Y + be^T Y 10841c3436cfSJed Brown = x1 + (be - b)^T Y 10851c3436cfSJed Brown 10861c3436cfSJed Brown so we can evaluate the method of different order even after the step has been optimistically completed. 10871c3436cfSJed Brown */ 1088d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSEvaluateStep_RosW(TS ts, PetscInt order, Vec U, PetscBool *done) 1089d71ae5a4SJacob Faibussowitsch { 10901c3436cfSJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 10911c3436cfSJed Brown RosWTableau tab = ros->tableau; 10921c3436cfSJed Brown PetscScalar *w = ros->work; 10931c3436cfSJed Brown PetscInt i; 10941c3436cfSJed Brown 10951c3436cfSJed Brown PetscFunctionBegin; 10961c3436cfSJed Brown if (order == tab->order) { 1097108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use standard completion formula */ 10989566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, U)); 1099de19f811SJed Brown for (i = 0; i < tab->s; i++) w[i] = tab->bt[i]; 11009566063dSJacob Faibussowitsch PetscCall(VecMAXPY(U, tab->s, w, ros->Y)); 11019566063dSJacob Faibussowitsch } else PetscCall(VecCopy(ts->vec_sol, U)); 11021c3436cfSJed Brown if (done) *done = PETSC_TRUE; 11033ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 11041c3436cfSJed Brown } else if (order == tab->order - 1) { 11051c3436cfSJed Brown if (!tab->bembedt) goto unavailable; 1106108c343cSJed Brown if (ros->status == TS_STEP_INCOMPLETE) { /* Use embedded completion formula */ 11079566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, U)); 1108de19f811SJed Brown for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i]; 11099566063dSJacob Faibussowitsch PetscCall(VecMAXPY(U, tab->s, w, ros->Y)); 1110108c343cSJed Brown } else { /* Use rollback-and-recomplete formula (bembedt - bt) */ 1111108c343cSJed Brown for (i = 0; i < tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 11129566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, U)); 11139566063dSJacob Faibussowitsch PetscCall(VecMAXPY(U, tab->s, w, ros->Y)); 11141c3436cfSJed Brown } 11151c3436cfSJed Brown if (done) *done = PETSC_TRUE; 11163ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 11171c3436cfSJed Brown } 11181c3436cfSJed Brown unavailable: 1119966bd95aSPierre Jolivet PetscCheck(done, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Rosenbrock-W '%s' of order %" PetscInt_FMT " cannot evaluate step at order %" PetscInt_FMT ". Consider using -ts_adapt_type none or a different method that has an embedded estimate.", 1120966bd95aSPierre Jolivet tab->name, tab->order, order); 1121966bd95aSPierre Jolivet *done = PETSC_FALSE; 11223ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 11231c3436cfSJed Brown } 11241c3436cfSJed Brown 1125d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSStep_RosW(TS ts) 1126d71ae5a4SJacob Faibussowitsch { 112761692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 112861692a83SJed Brown RosWTableau tab = ros->tableau; 1129e27a552bSJed Brown const PetscInt s = tab->s; 11301c3436cfSJed Brown const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv; 11310feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 1132c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 113361692a83SJed Brown PetscScalar *w = ros->work; 11347d4bf2deSEmil Constantinescu Vec *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage; 1135e27a552bSJed Brown SNES snes; 11361c3436cfSJed Brown TSAdapt adapt; 1137fecfb714SLisandro Dalcin PetscInt i, j, its, lits; 1138be5899b3SLisandro Dalcin PetscInt rejections = 0; 1139b39943a6SLisandro Dalcin PetscBool stageok, accept = PETSC_TRUE; 1140b39943a6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 1141f7f07198SBarry Smith PetscInt lag; 1142e27a552bSJed Brown 1143e27a552bSJed Brown PetscFunctionBegin; 114448a46eb9SPierre Jolivet if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, ros->vec_sol_prev)); 1145e27a552bSJed Brown 1146b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 1147b39943a6SLisandro Dalcin while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 11481c3436cfSJed Brown const PetscReal h = ts->time_step; 1149e27a552bSJed Brown for (i = 0; i < s; i++) { 11501c3436cfSJed Brown ros->stage_time = ts->ptime + h * ASum[i]; 11519566063dSJacob Faibussowitsch PetscCall(TSPreStage(ts, ros->stage_time)); 1152c17803e7SJed Brown if (GammaZeroDiag[i]) { 1153c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 1154b296d7d5SJed Brown ros->scoeff = 1.; 1155c17803e7SJed Brown } else { 1156c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 1157b296d7d5SJed Brown ros->scoeff = 1. / Gamma[i * s + i]; 1158fd96d5b0SEmil Constantinescu } 115961692a83SJed Brown 11609566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, Zstage)); 1161de19f811SJed Brown for (j = 0; j < i; j++) w[j] = At[i * s + j]; 11629566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Zstage, i, w, Y)); 116361692a83SJed Brown 116461692a83SJed Brown for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j]; 11659566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Zdot)); 11669566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Zdot, i, w, Y)); 116761692a83SJed Brown 1168e27a552bSJed Brown /* Initial guess taken from last stage */ 11699566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Y[i])); 117061692a83SJed Brown 11717d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 11729566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 117361692a83SJed Brown if (!ros->recompute_jacobian && !i) { 11749566063dSJacob Faibussowitsch PetscCall(SNESGetLagJacobian(snes, &lag)); 11756aad120cSJose E. Roman if (lag == 1) { /* use did not set a nontrivial lag, so lag over all stages */ 11769566063dSJacob Faibussowitsch PetscCall(SNESSetLagJacobian(snes, -2)); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */ 1177f7f07198SBarry Smith } 117861692a83SJed Brown } 11799566063dSJacob Faibussowitsch PetscCall(SNESSolve(snes, NULL, Y[i])); 1180ac530a7eSPierre Jolivet if (!ros->recompute_jacobian && i == s - 1 && lag == 1) PetscCall(SNESSetLagJacobian(snes, lag)); /* Set lag back to 1 so we know user did not set it */ 11819566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(snes, &its)); 11829566063dSJacob Faibussowitsch PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 11839371c9d4SSatish Balay ts->snes_its += its; 11849371c9d4SSatish Balay ts->ksp_its += lits; 11857d4bf2deSEmil Constantinescu } else { 11861ce71dffSSatish Balay Mat J, Jp; 11879566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Ydot)); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 11889566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE)); 11899566063dSJacob Faibussowitsch PetscCall(VecScale(Y[i], -1.0)); 11909566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y[i], -1.0, Zdot)); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 11910feba352SEmil Constantinescu 11929566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Zstage)); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 11930feba352SEmil Constantinescu for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j]; 11949566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Zstage, i, w, Y)); 1195fecfb714SLisandro Dalcin 1196fecfb714SLisandro Dalcin /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 11979566063dSJacob Faibussowitsch PetscCall(TSGetIJacobian(ts, &J, &Jp, NULL, NULL)); 11989566063dSJacob Faibussowitsch PetscCall(TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE)); 11999566063dSJacob Faibussowitsch PetscCall(MatMult(J, Zstage, Zdot)); 12009566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y[i], -1.0, Zdot)); 12015ef26d82SJed Brown ts->ksp_its += 1; 1202fecfb714SLisandro Dalcin 12039566063dSJacob Faibussowitsch PetscCall(VecScale(Y[i], h)); 12047d4bf2deSEmil Constantinescu } 12059566063dSJacob Faibussowitsch PetscCall(TSPostStage(ts, ros->stage_time, i, Y)); 12069566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 12079566063dSJacob Faibussowitsch PetscCall(TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok)); 1208fecfb714SLisandro Dalcin if (!stageok) goto reject_step; 1209e27a552bSJed Brown } 1210e27a552bSJed Brown 1211b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 12129566063dSJacob Faibussowitsch PetscCall(TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL)); 1213b39943a6SLisandro Dalcin ros->status = TS_STEP_PENDING; 12149566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 12159566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidatesClear(adapt)); 12169566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 12179566063dSJacob Faibussowitsch PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1218b39943a6SLisandro Dalcin ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1219b39943a6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 1220c61711c8SStefano Zampini PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol)); 1221be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1222be5899b3SLisandro Dalcin goto reject_step; 1223b39943a6SLisandro Dalcin } 1224b39943a6SLisandro Dalcin 1225e27a552bSJed Brown ts->ptime += ts->time_step; 1226cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 12271c3436cfSJed Brown break; 1228b39943a6SLisandro Dalcin 1229b39943a6SLisandro Dalcin reject_step: 12309371c9d4SSatish Balay ts->reject++; 12319371c9d4SSatish Balay accept = PETSC_FALSE; 1232be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1233b39943a6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 123463a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 12351c3436cfSJed Brown } 12361c3436cfSJed Brown } 12373ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1238e27a552bSJed Brown } 1239e27a552bSJed Brown 1240d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U) 1241d71ae5a4SJacob Faibussowitsch { 124261692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1243f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j; 1244f4aed992SEmil Constantinescu PetscReal h; 1245f4aed992SEmil Constantinescu PetscReal tt, t; 1246f4aed992SEmil Constantinescu PetscScalar *bt; 1247f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1248f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1249f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1250f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1251e27a552bSJed Brown 1252e27a552bSJed Brown PetscFunctionBegin; 12533c633725SBarry Smith PetscCheck(Bt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRosW %s does not have an interpolation formula", ros->tableau->name); 1254f4aed992SEmil Constantinescu 1255f4aed992SEmil Constantinescu switch (ros->status) { 1256f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1257f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1258f4aed992SEmil Constantinescu h = ts->time_step; 1259f4aed992SEmil Constantinescu t = (itime - ts->ptime) / h; 1260f4aed992SEmil Constantinescu break; 1261f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1262be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1263f4aed992SEmil Constantinescu t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1264f4aed992SEmil Constantinescu break; 1265d71ae5a4SJacob Faibussowitsch default: 1266d71ae5a4SJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1267f4aed992SEmil Constantinescu } 12689566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(s, &bt)); 1269f4aed992SEmil Constantinescu for (i = 0; i < s; i++) bt[i] = 0; 1270f4aed992SEmil Constantinescu for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1271ad540459SPierre Jolivet for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt; 1272f4aed992SEmil Constantinescu } 1273f4aed992SEmil Constantinescu 1274f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1275f9c1d6abSBarry Smith /* U <- 0*/ 12769566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(U)); 1277f9c1d6abSBarry Smith /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 12783ca35412SEmil Constantinescu for (j = 0; j < s; j++) w[j] = 0; 12793ca35412SEmil Constantinescu for (j = 0; j < s; j++) { 1280ad540459SPierre Jolivet for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j]; 12813ca35412SEmil Constantinescu } 12829566063dSJacob Faibussowitsch PetscCall(VecMAXPY(U, i, w, Y)); 1283be5899b3SLisandro Dalcin /* U <- y(t) + U */ 12849566063dSJacob Faibussowitsch PetscCall(VecAXPY(U, 1, ros->vec_sol_prev)); 1285f4aed992SEmil Constantinescu 12869566063dSJacob Faibussowitsch PetscCall(PetscFree(bt)); 12873ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1288e27a552bSJed Brown } 1289e27a552bSJed Brown 1290e27a552bSJed Brown /*------------------------------------------------------------*/ 1291b39943a6SLisandro Dalcin 1292d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWTableauReset(TS ts) 1293d71ae5a4SJacob Faibussowitsch { 1294b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW *)ts->data; 1295b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1296b39943a6SLisandro Dalcin 1297b39943a6SLisandro Dalcin PetscFunctionBegin; 12983ba16761SJacob Faibussowitsch if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 12999566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ros->Y)); 13009566063dSJacob Faibussowitsch PetscCall(PetscFree(ros->work)); 13013ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1302b39943a6SLisandro Dalcin } 1303b39943a6SLisandro Dalcin 1304d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSReset_RosW(TS ts) 1305d71ae5a4SJacob Faibussowitsch { 130661692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1307e27a552bSJed Brown 1308e27a552bSJed Brown PetscFunctionBegin; 13099566063dSJacob Faibussowitsch PetscCall(TSRosWTableauReset(ts)); 13109566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Ydot)); 13119566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Ystage)); 13129566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Zdot)); 13139566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Zstage)); 13149566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->vec_sol_prev)); 13153ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1316e27a552bSJed Brown } 1317e27a552bSJed Brown 1318d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1319d71ae5a4SJacob Faibussowitsch { 1320d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW *)ts->data; 1321d5e6173cSPeter Brune 1322d5e6173cSPeter Brune PetscFunctionBegin; 1323d5e6173cSPeter Brune if (Ydot) { 1324ac530a7eSPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1325ac530a7eSPierre Jolivet else *Ydot = rw->Ydot; 1326d5e6173cSPeter Brune } 1327d5e6173cSPeter Brune if (Zdot) { 1328ac530a7eSPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1329ac530a7eSPierre Jolivet else *Zdot = rw->Zdot; 1330d5e6173cSPeter Brune } 1331d5e6173cSPeter Brune if (Ystage) { 1332ac530a7eSPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1333ac530a7eSPierre Jolivet else *Ystage = rw->Ystage; 1334d5e6173cSPeter Brune } 1335d5e6173cSPeter Brune if (Zstage) { 1336ac530a7eSPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1337ac530a7eSPierre Jolivet else *Zstage = rw->Zstage; 1338d5e6173cSPeter Brune } 13393ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1340d5e6173cSPeter Brune } 1341d5e6173cSPeter Brune 1342d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1343d71ae5a4SJacob Faibussowitsch { 1344d5e6173cSPeter Brune PetscFunctionBegin; 1345d5e6173cSPeter Brune if (Ydot) { 134648a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1347d5e6173cSPeter Brune } 1348d5e6173cSPeter Brune if (Zdot) { 134948a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1350d5e6173cSPeter Brune } 1351d5e6173cSPeter Brune if (Ystage) { 135248a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1353d5e6173cSPeter Brune } 1354d5e6173cSPeter Brune if (Zstage) { 135548a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1356d5e6173cSPeter Brune } 13573ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1358d5e6173cSPeter Brune } 1359d5e6173cSPeter Brune 1360*2a8381b2SBarry Smith static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, PetscCtx ctx) 1361d71ae5a4SJacob Faibussowitsch { 1362d5e6173cSPeter Brune PetscFunctionBegin; 13633ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1364d5e6173cSPeter Brune } 1365d5e6173cSPeter Brune 1366*2a8381b2SBarry Smith static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx) 1367d71ae5a4SJacob Faibussowitsch { 1368d5e6173cSPeter Brune TS ts = (TS)ctx; 1369d5e6173cSPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1370d5e6173cSPeter Brune Vec Ydotc, Zdotc, Ystagec, Zstagec; 1371d5e6173cSPeter Brune 1372d5e6173cSPeter Brune PetscFunctionBegin; 13739566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 13749566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 13759566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Ydot, Ydotc)); 13769566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Ydotc, rscale, Ydotc)); 13779566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Ystage, Ystagec)); 13789566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Ystagec, rscale, Ystagec)); 13799566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Zdot, Zdotc)); 13809566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Zdotc, rscale, Zdotc)); 13819566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Zstage, Zstagec)); 13829566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Zstagec, rscale, Zstagec)); 13839566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 13849566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 13853ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1386d5e6173cSPeter Brune } 1387d5e6173cSPeter Brune 1388*2a8381b2SBarry Smith static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, PetscCtx ctx) 1389d71ae5a4SJacob Faibussowitsch { 1390258e1594SPeter Brune PetscFunctionBegin; 13913ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1392258e1594SPeter Brune } 1393258e1594SPeter Brune 1394*2a8381b2SBarry Smith static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx) 1395d71ae5a4SJacob Faibussowitsch { 1396258e1594SPeter Brune TS ts = (TS)ctx; 1397258e1594SPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1398258e1594SPeter Brune Vec Ydots, Zdots, Ystages, Zstages; 1399258e1594SPeter Brune 1400258e1594SPeter Brune PetscFunctionBegin; 14019566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 14029566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 1403258e1594SPeter Brune 14049566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 14059566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 1406258e1594SPeter Brune 14079566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 14089566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 1409258e1594SPeter Brune 14109566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 14119566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 1412258e1594SPeter Brune 14139566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 14149566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 1415258e1594SPeter Brune 14169566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 14179566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 14183ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1419258e1594SPeter Brune } 1420258e1594SPeter Brune 1421d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts) 1422d71ae5a4SJacob Faibussowitsch { 142361692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1424d5e6173cSPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1425b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1426d5e6173cSPeter Brune DM dm, dmsave; 1427e27a552bSJed Brown 1428e27a552bSJed Brown PetscFunctionBegin; 14299566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 14309566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 14319566063dSJacob Faibussowitsch PetscCall(VecWAXPY(Ydot, shift, U, Zdot)); /* Ydot = shift*U + Zdot */ 14329566063dSJacob Faibussowitsch PetscCall(VecWAXPY(Ystage, 1.0, U, Zstage)); /* Ystage = U + Zstage */ 1433d5e6173cSPeter Brune dmsave = ts->dm; 1434d5e6173cSPeter Brune ts->dm = dm; 14359566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE)); 1436d5e6173cSPeter Brune ts->dm = dmsave; 14379566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 14383ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1439e27a552bSJed Brown } 1440e27a552bSJed Brown 1441d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts) 1442d71ae5a4SJacob Faibussowitsch { 144361692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1444d5e6173cSPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1445b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1446d5e6173cSPeter Brune DM dm, dmsave; 1447e27a552bSJed Brown 1448e27a552bSJed Brown PetscFunctionBegin; 144961692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 14509566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 14519566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1452d5e6173cSPeter Brune dmsave = ts->dm; 1453d5e6173cSPeter Brune ts->dm = dm; 14549566063dSJacob Faibussowitsch PetscCall(TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE)); 1455d5e6173cSPeter Brune ts->dm = dmsave; 14569566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 14573ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1458e27a552bSJed Brown } 1459e27a552bSJed Brown 1460d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWTableauSetUp(TS ts) 1461d71ae5a4SJacob Faibussowitsch { 1462b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW *)ts->data; 1463b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1464b39943a6SLisandro Dalcin 1465b39943a6SLisandro Dalcin PetscFunctionBegin; 14669566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y)); 14679566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(tab->s, &ros->work)); 14683ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1469b39943a6SLisandro Dalcin } 1470b39943a6SLisandro Dalcin 1471d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSSetUp_RosW(TS ts) 1472d71ae5a4SJacob Faibussowitsch { 147361692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1474d5e6173cSPeter Brune DM dm; 1475b39943a6SLisandro Dalcin SNES snes; 14768434afd1SBarry Smith TSRHSJacobianFn *rhsjacobian; 1477e27a552bSJed Brown 1478e27a552bSJed Brown PetscFunctionBegin; 14799566063dSJacob Faibussowitsch PetscCall(TSRosWTableauSetUp(ts)); 14809566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Ydot)); 14819566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Ystage)); 14829566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Zdot)); 14839566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Zstage)); 14849566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->vec_sol_prev)); 14859566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &dm)); 14869566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 14879566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1488b39943a6SLisandro Dalcin /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 14899566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 149048a46eb9SPierre Jolivet if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 14919566063dSJacob Faibussowitsch PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL)); 1492a3ab5968SHong Zhang if (rhsjacobian == TSComputeRHSJacobianConstant) { 1493a3ab5968SHong Zhang Mat Amat, Pmat; 1494a3ab5968SHong Zhang 1495a3ab5968SHong Zhang /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */ 14969566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL)); 1497a3ab5968SHong Zhang if (Amat && Amat == ts->Arhs) { 1498a3ab5968SHong Zhang if (Amat == Pmat) { 14999566063dSJacob Faibussowitsch PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 15009566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, Amat, Amat, NULL, NULL)); 1501a3ab5968SHong Zhang } else { 15029566063dSJacob Faibussowitsch PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 15039566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL)); 1504a3ab5968SHong Zhang if (Pmat && Pmat == ts->Brhs) { 15059566063dSJacob Faibussowitsch PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat)); 15069566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL)); 15079566063dSJacob Faibussowitsch PetscCall(MatDestroy(&Pmat)); 1508a3ab5968SHong Zhang } 1509a3ab5968SHong Zhang } 15109566063dSJacob Faibussowitsch PetscCall(MatDestroy(&Amat)); 1511a3ab5968SHong Zhang } 1512a3ab5968SHong Zhang } 15133ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1514e27a552bSJed Brown } 1515e27a552bSJed Brown /*------------------------------------------------------------*/ 1516e27a552bSJed Brown 1517ce78bad3SBarry Smith static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems PetscOptionsObject) 1518d71ae5a4SJacob Faibussowitsch { 151961692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1520b39943a6SLisandro Dalcin SNES snes; 1521e27a552bSJed Brown 1522e27a552bSJed Brown PetscFunctionBegin; 1523d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options"); 1524e27a552bSJed Brown { 152561692a83SJed Brown RosWTableauLink link; 1526e27a552bSJed Brown PetscInt count, choice; 1527e27a552bSJed Brown PetscBool flg; 1528e27a552bSJed Brown const char **namelist; 152961692a83SJed Brown 1530fbccb6d4SPierre Jolivet for (link = RosWTableauList, count = 0; link; link = link->next, count++); 15319566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(count, (char ***)&namelist)); 153261692a83SJed Brown for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 15339566063dSJacob Faibussowitsch PetscCall(PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg)); 15349566063dSJacob Faibussowitsch if (flg) PetscCall(TSRosWSetType(ts, namelist[choice])); 15359566063dSJacob Faibussowitsch PetscCall(PetscFree(namelist)); 153661692a83SJed Brown 15379566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL)); 1538b39943a6SLisandro Dalcin } 1539d0609cedSBarry Smith PetscOptionsHeadEnd(); 154061692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 15419566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 154248a46eb9SPierre Jolivet if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 15433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1544e27a552bSJed Brown } 1545e27a552bSJed Brown 1546d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer) 1547d71ae5a4SJacob Faibussowitsch { 154861692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 15499f196a02SMartin Diehl PetscBool isascii; 1550e27a552bSJed Brown 1551e27a552bSJed Brown PetscFunctionBegin; 15529f196a02SMartin Diehl PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 15539f196a02SMartin Diehl if (isascii) { 15549c334d8fSLisandro Dalcin RosWTableau tab = ros->tableau; 155519fd82e9SBarry Smith TSRosWType rostype; 15569c334d8fSLisandro Dalcin char buf[512]; 1557e408995aSJed Brown PetscInt i; 1558e408995aSJed Brown PetscReal abscissa[512]; 15599566063dSJacob Faibussowitsch PetscCall(TSRosWGetType(ts, &rostype)); 15609566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Rosenbrock-W %s\n", rostype)); 15619566063dSJacob Faibussowitsch PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum)); 15629566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A = %s\n", buf)); 15637edfa88dSJames Wright for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->GammaSum[i]; 15649566063dSJacob Faibussowitsch PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa)); 15659566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A+Gamma = %s\n", buf)); 1566e27a552bSJed Brown } 15673ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1568e27a552bSJed Brown } 1569e27a552bSJed Brown 1570d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer) 1571d71ae5a4SJacob Faibussowitsch { 15729200755eSBarry Smith SNES snes; 15739c334d8fSLisandro Dalcin TSAdapt adapt; 15749200755eSBarry Smith 15759200755eSBarry Smith PetscFunctionBegin; 15769566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 15779566063dSJacob Faibussowitsch PetscCall(TSAdaptLoad(adapt, viewer)); 15789566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 15799566063dSJacob Faibussowitsch PetscCall(SNESLoad(snes, viewer)); 15809200755eSBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 15819566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 15829566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 15833ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 15849200755eSBarry Smith } 15859200755eSBarry Smith 1586cc4c1da9SBarry Smith /*@ 1587bcf0153eSBarry Smith TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme 1588e27a552bSJed Brown 158920f4b53cSBarry Smith Logically Collective 1590e27a552bSJed Brown 1591d8d19677SJose E. Roman Input Parameters: 1592e27a552bSJed Brown + ts - timestepping context 1593b92453a8SLisandro Dalcin - roswtype - type of Rosenbrock-W scheme 1594e27a552bSJed Brown 1595020d8f30SJed Brown Level: beginner 1596e27a552bSJed Brown 15971cc06b55SBarry Smith .seealso: [](ch_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3` 1598e27a552bSJed Brown @*/ 1599d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype) 1600d71ae5a4SJacob Faibussowitsch { 1601e27a552bSJed Brown PetscFunctionBegin; 1602e27a552bSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 16034f572ea9SToby Isaac PetscAssertPointer(roswtype, 2); 1604cac4c232SBarry Smith PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype)); 16053ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1606e27a552bSJed Brown } 1607e27a552bSJed Brown 1608cc4c1da9SBarry Smith /*@ 160961692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1610e27a552bSJed Brown 161120f4b53cSBarry Smith Logically Collective 1612e27a552bSJed Brown 1613e27a552bSJed Brown Input Parameter: 1614e27a552bSJed Brown . ts - timestepping context 1615e27a552bSJed Brown 1616e27a552bSJed Brown Output Parameter: 161761692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1618e27a552bSJed Brown 1619e27a552bSJed Brown Level: intermediate 1620e27a552bSJed Brown 16211cc06b55SBarry Smith .seealso: [](ch_ts), `TSRosWType`, `TSRosWSetType()` 1622e27a552bSJed Brown @*/ 1623d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype) 1624d71ae5a4SJacob Faibussowitsch { 1625e27a552bSJed Brown PetscFunctionBegin; 1626e27a552bSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1627cac4c232SBarry Smith PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype)); 16283ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1629e27a552bSJed Brown } 1630e27a552bSJed Brown 1631cc4c1da9SBarry Smith /*@ 163261692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1633e27a552bSJed Brown 163420f4b53cSBarry Smith Logically Collective 1635e27a552bSJed Brown 1636d8d19677SJose E. Roman Input Parameters: 1637e27a552bSJed Brown + ts - timestepping context 1638bcf0153eSBarry Smith - flg - `PETSC_TRUE` to recompute the Jacobian at each stage 1639e27a552bSJed Brown 1640e27a552bSJed Brown Level: intermediate 1641e27a552bSJed Brown 16421cc06b55SBarry Smith .seealso: [](ch_ts), `TSRosWType`, `TSRosWGetType()` 1643e27a552bSJed Brown @*/ 1644d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg) 1645d71ae5a4SJacob Faibussowitsch { 1646e27a552bSJed Brown PetscFunctionBegin; 1647e27a552bSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1648cac4c232SBarry Smith PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg)); 16493ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1650e27a552bSJed Brown } 1651e27a552bSJed Brown 1652d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype) 1653d71ae5a4SJacob Faibussowitsch { 165461692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1655e27a552bSJed Brown 1656e27a552bSJed Brown PetscFunctionBegin; 165761692a83SJed Brown *rostype = ros->tableau->name; 16583ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1659e27a552bSJed Brown } 1660ef20d060SBarry Smith 1661d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype) 1662d71ae5a4SJacob Faibussowitsch { 166361692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1664e27a552bSJed Brown PetscBool match; 166561692a83SJed Brown RosWTableauLink link; 1666e27a552bSJed Brown 1667e27a552bSJed Brown PetscFunctionBegin; 166861692a83SJed Brown if (ros->tableau) { 16699566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(ros->tableau->name, rostype, &match)); 16703ba16761SJacob Faibussowitsch if (match) PetscFunctionReturn(PETSC_SUCCESS); 1671e27a552bSJed Brown } 167261692a83SJed Brown for (link = RosWTableauList; link; link = link->next) { 16739566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(link->tab.name, rostype, &match)); 1674e27a552bSJed Brown if (match) { 16759566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSRosWTableauReset(ts)); 167661692a83SJed Brown ros->tableau = &link->tab; 16779566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSRosWTableauSetUp(ts)); 16782ffb9264SLisandro Dalcin ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 16793ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1680e27a552bSJed Brown } 1681e27a552bSJed Brown } 168298921bdaSJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype); 1683e27a552bSJed Brown } 168461692a83SJed Brown 1685d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg) 1686d71ae5a4SJacob Faibussowitsch { 168761692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1688e27a552bSJed Brown 1689e27a552bSJed Brown PetscFunctionBegin; 169061692a83SJed Brown ros->recompute_jacobian = flg; 16913ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1692e27a552bSJed Brown } 1693e27a552bSJed Brown 1694d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSDestroy_RosW(TS ts) 1695d71ae5a4SJacob Faibussowitsch { 1696b3a6b972SJed Brown PetscFunctionBegin; 16979566063dSJacob Faibussowitsch PetscCall(TSReset_RosW(ts)); 1698b3a6b972SJed Brown if (ts->dm) { 16999566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 17009566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1701b3a6b972SJed Brown } 17029566063dSJacob Faibussowitsch PetscCall(PetscFree(ts->data)); 17039566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL)); 17049566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL)); 17059566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL)); 17063ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1707b3a6b972SJed Brown } 1708d5e6173cSPeter Brune 1709e27a552bSJed Brown /*MC 1710020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1711e27a552bSJed Brown 1712e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1713e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1714bcf0153eSBarry Smith of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 1715bcf0153eSBarry Smith 1716bcf0153eSBarry Smith Level: beginner 1717e27a552bSJed Brown 1718e27a552bSJed Brown Notes: 1719efa39862SBarry Smith This is an IMEX method. 1720efa39862SBarry Smith 172161692a83SJed Brown This method currently only works with autonomous ODE and DAE. 172261692a83SJed Brown 1723bcf0153eSBarry Smith Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear. 1724d0685a90SJed Brown 1725efa39862SBarry 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`. 17263d5a8a6aSBarry Smith 1727e94cfbe0SPatrick Sanan Developer Notes: 172861692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 1729efa39862SBarry Smith $$ 1730efa39862SBarry Smith \dot{u} = f(u) 1731efa39862SBarry Smith $$ 173261692a83SJed Brown by the stage equations 1733efa39862SBarry Smith $$ 1734efa39862SBarry Smith k_i = h f(u_0 + \sum_j \alpha_{ij} k_j) + h J \sum_j \gamma_{ij} k_j 1735efa39862SBarry Smith $$ 173661692a83SJed Brown and step completion formula 1737efa39862SBarry Smith $$ 1738efa39862SBarry Smith u_1 = u_0 + \sum_j b_j k_j 1739efa39862SBarry Smith $$ 1740efa39862SBarry Smith with step size $h$ and coefficients $\alpha_{ij}$, $\gamma_{ij}$, and $b_i$. Implementing the method in this form would require $f(u)$ 1741efa39862SBarry Smith and the Jacobian $J$ to be available, in addition to the shifted matrix $I - h \gamma_{ii} J$. Following Hairer and Wanner, 174261692a83SJed Brown we define new variables for the stage equations 1743efa39862SBarry Smith $$ 1744efa39862SBarry Smith y_i = \gamma_{ij} k_j 1745efa39862SBarry Smith $$ 1746efa39862SBarry Smith The $ k_j $ can be recovered because $\Gamma$ is invertible. Let $C$ be the lower triangular part of $\Gamma^{-1}$ and define 1747efa39862SBarry Smith $$ 1748efa39862SBarry Smith A = \Alpha \Gamma^{-1}, bt^T = b^T \Gamma^{-1} 1749efa39862SBarry Smith $$ 175061692a83SJed Brown to rewrite the method as 1751efa39862SBarry Smith $$ 1752efa39862SBarry 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 \\ 1753efa39862SBarry Smith u_1 = u_0 + \sum_j bt_j y_j 1754efa39862SBarry Smith $$ 175561692a83SJed Brown 1756efa39862SBarry Smith where we have introduced the mass matrix $M$. Continue by defining 1757efa39862SBarry Smith $$ 1758efa39862SBarry Smith \dot{y}_i = 1/(h \gamma_ii) y_i - \sum_j (c_{ij}/h) y_j 1759efa39862SBarry Smith $$ 176061692a83SJed Brown or, more compactly in tensor notation 1761efa39862SBarry Smith $$ 1762efa39862SBarry Smith \dot{Y} = 1/h (Gamma^{-1} \otimes I) Y . 1763efa39862SBarry Smith $$ 1764efa39862SBarry Smith Note that $\Gamma^{-1}$ is lower triangular. With this definition of $\dot{Y} in terms of known quantities and the current 1765efa39862SBarry Smith stage $y_i$, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 176661692a83SJed Brown equation 1767efa39862SBarry Smith $$ 1768efa39862SBarry Smith g(u_0 + \sum_j a_{ij} y_j + y_i, \dot{y}_i) = 0 1769efa39862SBarry Smith $$ 1770efa39862SBarry Smith with initial guess $y_i = 0$. 1771e27a552bSJed Brown 1772efa39862SBarry Smith .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, 1773efa39862SBarry Smith `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, 1774efa39862SBarry Smith `TSROSW4L`, `TSType` 1775e27a552bSJed Brown M*/ 1776d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1777d71ae5a4SJacob Faibussowitsch { 177861692a83SJed Brown TS_RosW *ros; 1779e27a552bSJed Brown 1780e27a552bSJed Brown PetscFunctionBegin; 17819566063dSJacob Faibussowitsch PetscCall(TSRosWInitializePackage()); 1782e27a552bSJed Brown 1783e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1784e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1785e27a552bSJed Brown ts->ops->view = TSView_RosW; 17869200755eSBarry Smith ts->ops->load = TSLoad_RosW; 1787e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1788e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1789e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 17901c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 1791e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1792e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1793e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1794e27a552bSJed Brown 1795efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 1796efd4aadfSBarry Smith 17974dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&ros)); 179861692a83SJed Brown ts->data = (void *)ros; 1799e27a552bSJed Brown 18009566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW)); 18019566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW)); 18029566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW)); 1803b39943a6SLisandro Dalcin 18049566063dSJacob Faibussowitsch PetscCall(TSRosWSetType(ts, TSRosWDefault)); 18053ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1806e27a552bSJed Brown } 1807