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 959*bfe80ac4SPierre 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: 11191c3436cfSJed Brown if (done) *done = PETSC_FALSE; 11209371c9d4SSatish Balay else 11219371c9d4SSatish Balay SETERRQ(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.", tab->name, 11229371c9d4SSatish Balay tab->order, order); 11233ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 11241c3436cfSJed Brown } 11251c3436cfSJed Brown 1126d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSStep_RosW(TS ts) 1127d71ae5a4SJacob Faibussowitsch { 112861692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 112961692a83SJed Brown RosWTableau tab = ros->tableau; 1130e27a552bSJed Brown const PetscInt s = tab->s; 11311c3436cfSJed Brown const PetscReal *At = tab->At, *Gamma = tab->Gamma, *ASum = tab->ASum, *GammaInv = tab->GammaInv; 11320feba352SEmil Constantinescu const PetscReal *GammaExplicitCorr = tab->GammaExplicitCorr; 1133c17803e7SJed Brown const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 113461692a83SJed Brown PetscScalar *w = ros->work; 11357d4bf2deSEmil Constantinescu Vec *Y = ros->Y, Ydot = ros->Ydot, Zdot = ros->Zdot, Zstage = ros->Zstage; 1136e27a552bSJed Brown SNES snes; 11371c3436cfSJed Brown TSAdapt adapt; 1138fecfb714SLisandro Dalcin PetscInt i, j, its, lits; 1139be5899b3SLisandro Dalcin PetscInt rejections = 0; 1140b39943a6SLisandro Dalcin PetscBool stageok, accept = PETSC_TRUE; 1141b39943a6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 1142f7f07198SBarry Smith PetscInt lag; 1143e27a552bSJed Brown 1144e27a552bSJed Brown PetscFunctionBegin; 114548a46eb9SPierre Jolivet if (!ts->steprollback) PetscCall(VecCopy(ts->vec_sol, ros->vec_sol_prev)); 1146e27a552bSJed Brown 1147b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 1148b39943a6SLisandro Dalcin while (!ts->reason && ros->status != TS_STEP_COMPLETE) { 11491c3436cfSJed Brown const PetscReal h = ts->time_step; 1150e27a552bSJed Brown for (i = 0; i < s; i++) { 11511c3436cfSJed Brown ros->stage_time = ts->ptime + h * ASum[i]; 11529566063dSJacob Faibussowitsch PetscCall(TSPreStage(ts, ros->stage_time)); 1153c17803e7SJed Brown if (GammaZeroDiag[i]) { 1154c17803e7SJed Brown ros->stage_explicit = PETSC_TRUE; 1155b296d7d5SJed Brown ros->scoeff = 1.; 1156c17803e7SJed Brown } else { 1157c17803e7SJed Brown ros->stage_explicit = PETSC_FALSE; 1158b296d7d5SJed Brown ros->scoeff = 1. / Gamma[i * s + i]; 1159fd96d5b0SEmil Constantinescu } 116061692a83SJed Brown 11619566063dSJacob Faibussowitsch PetscCall(VecCopy(ts->vec_sol, Zstage)); 1162de19f811SJed Brown for (j = 0; j < i; j++) w[j] = At[i * s + j]; 11639566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Zstage, i, w, Y)); 116461692a83SJed Brown 116561692a83SJed Brown for (j = 0; j < i; j++) w[j] = 1. / h * GammaInv[i * s + j]; 11669566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Zdot)); 11679566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Zdot, i, w, Y)); 116861692a83SJed Brown 1169e27a552bSJed Brown /* Initial guess taken from last stage */ 11709566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Y[i])); 117161692a83SJed Brown 11727d4bf2deSEmil Constantinescu if (!ros->stage_explicit) { 11739566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 117461692a83SJed Brown if (!ros->recompute_jacobian && !i) { 11759566063dSJacob Faibussowitsch PetscCall(SNESGetLagJacobian(snes, &lag)); 11766aad120cSJose E. Roman if (lag == 1) { /* use did not set a nontrivial lag, so lag over all stages */ 11779566063dSJacob Faibussowitsch PetscCall(SNESSetLagJacobian(snes, -2)); /* Recompute the Jacobian on this solve, but not again for the rest of the stages */ 1178f7f07198SBarry Smith } 117961692a83SJed Brown } 11809566063dSJacob Faibussowitsch PetscCall(SNESSolve(snes, NULL, Y[i])); 11819371c9d4SSatish Balay 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 */ } 11829566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(snes, &its)); 11839566063dSJacob Faibussowitsch PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 11849371c9d4SSatish Balay ts->snes_its += its; 11859371c9d4SSatish Balay ts->ksp_its += lits; 11867d4bf2deSEmil Constantinescu } else { 11871ce71dffSSatish Balay Mat J, Jp; 11889566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Ydot)); /* Evaluate Y[i]=G(t,Ydot=0,Zstage) */ 11899566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ros->stage_time, Zstage, Ydot, Y[i], PETSC_FALSE)); 11909566063dSJacob Faibussowitsch PetscCall(VecScale(Y[i], -1.0)); 11919566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y[i], -1.0, Zdot)); /*Y[i] = F(Zstage)-Zdot[=GammaInv*Y]*/ 11920feba352SEmil Constantinescu 11939566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(Zstage)); /* Zstage = GammaExplicitCorr[i,j] * Y[j] */ 11940feba352SEmil Constantinescu for (j = 0; j < i; j++) w[j] = GammaExplicitCorr[i * s + j]; 11959566063dSJacob Faibussowitsch PetscCall(VecMAXPY(Zstage, i, w, Y)); 1196fecfb714SLisandro Dalcin 1197fecfb714SLisandro Dalcin /* Y[i] = Y[i] + Jac*Zstage[=Jac*GammaExplicitCorr[i,j] * Y[j]] */ 11989566063dSJacob Faibussowitsch PetscCall(TSGetIJacobian(ts, &J, &Jp, NULL, NULL)); 11999566063dSJacob Faibussowitsch PetscCall(TSComputeIJacobian(ts, ros->stage_time, ts->vec_sol, Ydot, 0, J, Jp, PETSC_FALSE)); 12009566063dSJacob Faibussowitsch PetscCall(MatMult(J, Zstage, Zdot)); 12019566063dSJacob Faibussowitsch PetscCall(VecAXPY(Y[i], -1.0, Zdot)); 12025ef26d82SJed Brown ts->ksp_its += 1; 1203fecfb714SLisandro Dalcin 12049566063dSJacob Faibussowitsch PetscCall(VecScale(Y[i], h)); 12057d4bf2deSEmil Constantinescu } 12069566063dSJacob Faibussowitsch PetscCall(TSPostStage(ts, ros->stage_time, i, Y)); 12079566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 12089566063dSJacob Faibussowitsch PetscCall(TSAdaptCheckStage(adapt, ts, ros->stage_time, Y[i], &stageok)); 1209fecfb714SLisandro Dalcin if (!stageok) goto reject_step; 1210e27a552bSJed Brown } 1211e27a552bSJed Brown 1212b39943a6SLisandro Dalcin ros->status = TS_STEP_INCOMPLETE; 12139566063dSJacob Faibussowitsch PetscCall(TSEvaluateStep_RosW(ts, tab->order, ts->vec_sol, NULL)); 1214b39943a6SLisandro Dalcin ros->status = TS_STEP_PENDING; 12159566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 12169566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidatesClear(adapt)); 12179566063dSJacob Faibussowitsch PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 12189566063dSJacob Faibussowitsch PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 1219b39943a6SLisandro Dalcin ros->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 1220b39943a6SLisandro Dalcin if (!accept) { /* Roll back the current step */ 1221c61711c8SStefano Zampini PetscCall(VecCopy(ts->vec_sol0, ts->vec_sol)); 1222be5899b3SLisandro Dalcin ts->time_step = next_time_step; 1223be5899b3SLisandro Dalcin goto reject_step; 1224b39943a6SLisandro Dalcin } 1225b39943a6SLisandro Dalcin 1226e27a552bSJed Brown ts->ptime += ts->time_step; 1227cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 12281c3436cfSJed Brown break; 1229b39943a6SLisandro Dalcin 1230b39943a6SLisandro Dalcin reject_step: 12319371c9d4SSatish Balay ts->reject++; 12329371c9d4SSatish Balay accept = PETSC_FALSE; 1233be5899b3SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 1234b39943a6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 123563a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 12361c3436cfSJed Brown } 12371c3436cfSJed Brown } 12383ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1239e27a552bSJed Brown } 1240e27a552bSJed Brown 1241d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSInterpolate_RosW(TS ts, PetscReal itime, Vec U) 1242d71ae5a4SJacob Faibussowitsch { 124361692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1244f4aed992SEmil Constantinescu PetscInt s = ros->tableau->s, pinterp = ros->tableau->pinterp, i, j; 1245f4aed992SEmil Constantinescu PetscReal h; 1246f4aed992SEmil Constantinescu PetscReal tt, t; 1247f4aed992SEmil Constantinescu PetscScalar *bt; 1248f4aed992SEmil Constantinescu const PetscReal *Bt = ros->tableau->binterpt; 1249f4aed992SEmil Constantinescu const PetscReal *GammaInv = ros->tableau->GammaInv; 1250f4aed992SEmil Constantinescu PetscScalar *w = ros->work; 1251f4aed992SEmil Constantinescu Vec *Y = ros->Y; 1252e27a552bSJed Brown 1253e27a552bSJed Brown PetscFunctionBegin; 12543c633725SBarry Smith PetscCheck(Bt, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRosW %s does not have an interpolation formula", ros->tableau->name); 1255f4aed992SEmil Constantinescu 1256f4aed992SEmil Constantinescu switch (ros->status) { 1257f4aed992SEmil Constantinescu case TS_STEP_INCOMPLETE: 1258f4aed992SEmil Constantinescu case TS_STEP_PENDING: 1259f4aed992SEmil Constantinescu h = ts->time_step; 1260f4aed992SEmil Constantinescu t = (itime - ts->ptime) / h; 1261f4aed992SEmil Constantinescu break; 1262f4aed992SEmil Constantinescu case TS_STEP_COMPLETE: 1263be5899b3SLisandro Dalcin h = ts->ptime - ts->ptime_prev; 1264f4aed992SEmil Constantinescu t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1265f4aed992SEmil Constantinescu break; 1266d71ae5a4SJacob Faibussowitsch default: 1267d71ae5a4SJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1268f4aed992SEmil Constantinescu } 12699566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(s, &bt)); 1270f4aed992SEmil Constantinescu for (i = 0; i < s; i++) bt[i] = 0; 1271f4aed992SEmil Constantinescu for (j = 0, tt = t; j < pinterp; j++, tt *= t) { 1272ad540459SPierre Jolivet for (i = 0; i < s; i++) bt[i] += Bt[i * pinterp + j] * tt; 1273f4aed992SEmil Constantinescu } 1274f4aed992SEmil Constantinescu 1275f4aed992SEmil Constantinescu /* y(t+tt*h) = y(t) + Sum bt(tt) * GammaInv * Ydot */ 1276f9c1d6abSBarry Smith /* U <- 0*/ 12779566063dSJacob Faibussowitsch PetscCall(VecZeroEntries(U)); 1278f9c1d6abSBarry Smith /* U <- Sum bt_i * GammaInv(i,1:i) * Y(1:i) */ 12793ca35412SEmil Constantinescu for (j = 0; j < s; j++) w[j] = 0; 12803ca35412SEmil Constantinescu for (j = 0; j < s; j++) { 1281ad540459SPierre Jolivet for (i = j; i < s; i++) w[j] += bt[i] * GammaInv[i * s + j]; 12823ca35412SEmil Constantinescu } 12839566063dSJacob Faibussowitsch PetscCall(VecMAXPY(U, i, w, Y)); 1284be5899b3SLisandro Dalcin /* U <- y(t) + U */ 12859566063dSJacob Faibussowitsch PetscCall(VecAXPY(U, 1, ros->vec_sol_prev)); 1286f4aed992SEmil Constantinescu 12879566063dSJacob Faibussowitsch PetscCall(PetscFree(bt)); 12883ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1289e27a552bSJed Brown } 1290e27a552bSJed Brown 1291e27a552bSJed Brown /*------------------------------------------------------------*/ 1292b39943a6SLisandro Dalcin 1293d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWTableauReset(TS ts) 1294d71ae5a4SJacob Faibussowitsch { 1295b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW *)ts->data; 1296b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1297b39943a6SLisandro Dalcin 1298b39943a6SLisandro Dalcin PetscFunctionBegin; 12993ba16761SJacob Faibussowitsch if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 13009566063dSJacob Faibussowitsch PetscCall(VecDestroyVecs(tab->s, &ros->Y)); 13019566063dSJacob Faibussowitsch PetscCall(PetscFree(ros->work)); 13023ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1303b39943a6SLisandro Dalcin } 1304b39943a6SLisandro Dalcin 1305d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSReset_RosW(TS ts) 1306d71ae5a4SJacob Faibussowitsch { 130761692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1308e27a552bSJed Brown 1309e27a552bSJed Brown PetscFunctionBegin; 13109566063dSJacob Faibussowitsch PetscCall(TSRosWTableauReset(ts)); 13119566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Ydot)); 13129566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Ystage)); 13139566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Zdot)); 13149566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->Zstage)); 13159566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ros->vec_sol_prev)); 13163ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1317e27a552bSJed Brown } 1318e27a552bSJed Brown 1319d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWGetVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1320d71ae5a4SJacob Faibussowitsch { 1321d5e6173cSPeter Brune TS_RosW *rw = (TS_RosW *)ts->data; 1322d5e6173cSPeter Brune 1323d5e6173cSPeter Brune PetscFunctionBegin; 1324d5e6173cSPeter Brune if (Ydot) { 1325d5e6173cSPeter Brune if (dm && dm != ts->dm) { 13269566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1327d5e6173cSPeter Brune } else *Ydot = rw->Ydot; 1328d5e6173cSPeter Brune } 1329d5e6173cSPeter Brune if (Zdot) { 1330d5e6173cSPeter Brune if (dm && dm != ts->dm) { 13319566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1332d5e6173cSPeter Brune } else *Zdot = rw->Zdot; 1333d5e6173cSPeter Brune } 1334d5e6173cSPeter Brune if (Ystage) { 1335d5e6173cSPeter Brune if (dm && dm != ts->dm) { 13369566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1337d5e6173cSPeter Brune } else *Ystage = rw->Ystage; 1338d5e6173cSPeter Brune } 1339d5e6173cSPeter Brune if (Zstage) { 1340d5e6173cSPeter Brune if (dm && dm != ts->dm) { 13419566063dSJacob Faibussowitsch PetscCall(DMGetNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1342d5e6173cSPeter Brune } else *Zstage = rw->Zstage; 1343d5e6173cSPeter Brune } 13443ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1345d5e6173cSPeter Brune } 1346d5e6173cSPeter Brune 1347d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWRestoreVecs(TS ts, DM dm, Vec *Ydot, Vec *Zdot, Vec *Ystage, Vec *Zstage) 1348d71ae5a4SJacob Faibussowitsch { 1349d5e6173cSPeter Brune PetscFunctionBegin; 1350d5e6173cSPeter Brune if (Ydot) { 135148a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ydot", Ydot)); 1352d5e6173cSPeter Brune } 1353d5e6173cSPeter Brune if (Zdot) { 135448a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zdot", Zdot)); 1355d5e6173cSPeter Brune } 1356d5e6173cSPeter Brune if (Ystage) { 135748a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Ystage", Ystage)); 1358d5e6173cSPeter Brune } 1359d5e6173cSPeter Brune if (Zstage) { 136048a46eb9SPierre Jolivet if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSRosW_Zstage", Zstage)); 1361d5e6173cSPeter Brune } 13623ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1363d5e6173cSPeter Brune } 1364d5e6173cSPeter Brune 1365d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMCoarsenHook_TSRosW(DM fine, DM coarse, void *ctx) 1366d71ae5a4SJacob Faibussowitsch { 1367d5e6173cSPeter Brune PetscFunctionBegin; 13683ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1369d5e6173cSPeter Brune } 1370d5e6173cSPeter Brune 1371d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMRestrictHook_TSRosW(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1372d71ae5a4SJacob Faibussowitsch { 1373d5e6173cSPeter Brune TS ts = (TS)ctx; 1374d5e6173cSPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1375d5e6173cSPeter Brune Vec Ydotc, Zdotc, Ystagec, Zstagec; 1376d5e6173cSPeter Brune 1377d5e6173cSPeter Brune PetscFunctionBegin; 13789566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 13799566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 13809566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Ydot, Ydotc)); 13819566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Ydotc, rscale, Ydotc)); 13829566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Ystage, Ystagec)); 13839566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Ystagec, rscale, Ystagec)); 13849566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Zdot, Zdotc)); 13859566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Zdotc, rscale, Zdotc)); 13869566063dSJacob Faibussowitsch PetscCall(MatRestrict(restrct, Zstage, Zstagec)); 13879566063dSJacob Faibussowitsch PetscCall(VecPointwiseMult(Zstagec, rscale, Zstagec)); 13889566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, fine, &Ydot, &Ystage, &Zdot, &Zstage)); 13899566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, coarse, &Ydotc, &Ystagec, &Zdotc, &Zstagec)); 13903ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1391d5e6173cSPeter Brune } 1392d5e6173cSPeter Brune 1393d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainHook_TSRosW(DM fine, DM coarse, void *ctx) 1394d71ae5a4SJacob Faibussowitsch { 1395258e1594SPeter Brune PetscFunctionBegin; 13963ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1397258e1594SPeter Brune } 1398258e1594SPeter Brune 1399d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMSubDomainRestrictHook_TSRosW(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1400d71ae5a4SJacob Faibussowitsch { 1401258e1594SPeter Brune TS ts = (TS)ctx; 1402258e1594SPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1403258e1594SPeter Brune Vec Ydots, Zdots, Ystages, Zstages; 1404258e1594SPeter Brune 1405258e1594SPeter Brune PetscFunctionBegin; 14069566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 14079566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 1408258e1594SPeter Brune 14099566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 14109566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Ydot, Ydots, INSERT_VALUES, SCATTER_FORWARD)); 1411258e1594SPeter Brune 14129566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 14139566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Ystage, Ystages, INSERT_VALUES, SCATTER_FORWARD)); 1414258e1594SPeter Brune 14159566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 14169566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Zdot, Zdots, INSERT_VALUES, SCATTER_FORWARD)); 1417258e1594SPeter Brune 14189566063dSJacob Faibussowitsch PetscCall(VecScatterBegin(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 14199566063dSJacob Faibussowitsch PetscCall(VecScatterEnd(gscat, Zstage, Zstages, INSERT_VALUES, SCATTER_FORWARD)); 1420258e1594SPeter Brune 14219566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Ystage, &Zdot, &Zstage)); 14229566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, subdm, &Ydots, &Ystages, &Zdots, &Zstages)); 14233ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1424258e1594SPeter Brune } 1425258e1594SPeter Brune 1426d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormFunction_RosW(SNES snes, Vec U, Vec F, TS ts) 1427d71ae5a4SJacob Faibussowitsch { 142861692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1429d5e6173cSPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1430b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1431d5e6173cSPeter Brune DM dm, dmsave; 1432e27a552bSJed Brown 1433e27a552bSJed Brown PetscFunctionBegin; 14349566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 14359566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 14369566063dSJacob Faibussowitsch PetscCall(VecWAXPY(Ydot, shift, U, Zdot)); /* Ydot = shift*U + Zdot */ 14379566063dSJacob Faibussowitsch PetscCall(VecWAXPY(Ystage, 1.0, U, Zstage)); /* Ystage = U + Zstage */ 1438d5e6173cSPeter Brune dmsave = ts->dm; 1439d5e6173cSPeter Brune ts->dm = dm; 14409566063dSJacob Faibussowitsch PetscCall(TSComputeIFunction(ts, ros->stage_time, Ystage, Ydot, F, PETSC_FALSE)); 1441d5e6173cSPeter Brune ts->dm = dmsave; 14429566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 14433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1444e27a552bSJed Brown } 1445e27a552bSJed Brown 1446d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes, Vec U, Mat A, Mat B, TS ts) 1447d71ae5a4SJacob Faibussowitsch { 144861692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1449d5e6173cSPeter Brune Vec Ydot, Zdot, Ystage, Zstage; 1450b296d7d5SJed Brown PetscReal shift = ros->scoeff / ts->time_step; 1451d5e6173cSPeter Brune DM dm, dmsave; 1452e27a552bSJed Brown 1453e27a552bSJed Brown PetscFunctionBegin; 145461692a83SJed Brown /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 14559566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &dm)); 14569566063dSJacob Faibussowitsch PetscCall(TSRosWGetVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 1457d5e6173cSPeter Brune dmsave = ts->dm; 1458d5e6173cSPeter Brune ts->dm = dm; 14599566063dSJacob Faibussowitsch PetscCall(TSComputeIJacobian(ts, ros->stage_time, Ystage, Ydot, shift, A, B, PETSC_TRUE)); 1460d5e6173cSPeter Brune ts->dm = dmsave; 14619566063dSJacob Faibussowitsch PetscCall(TSRosWRestoreVecs(ts, dm, &Ydot, &Zdot, &Ystage, &Zstage)); 14623ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1463e27a552bSJed Brown } 1464e27a552bSJed Brown 1465d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWTableauSetUp(TS ts) 1466d71ae5a4SJacob Faibussowitsch { 1467b39943a6SLisandro Dalcin TS_RosW *ros = (TS_RosW *)ts->data; 1468b39943a6SLisandro Dalcin RosWTableau tab = ros->tableau; 1469b39943a6SLisandro Dalcin 1470b39943a6SLisandro Dalcin PetscFunctionBegin; 14719566063dSJacob Faibussowitsch PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &ros->Y)); 14729566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(tab->s, &ros->work)); 14733ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1474b39943a6SLisandro Dalcin } 1475b39943a6SLisandro Dalcin 1476d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSSetUp_RosW(TS ts) 1477d71ae5a4SJacob Faibussowitsch { 147861692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1479d5e6173cSPeter Brune DM dm; 1480b39943a6SLisandro Dalcin SNES snes; 14818434afd1SBarry Smith TSRHSJacobianFn *rhsjacobian; 1482e27a552bSJed Brown 1483e27a552bSJed Brown PetscFunctionBegin; 14849566063dSJacob Faibussowitsch PetscCall(TSRosWTableauSetUp(ts)); 14859566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Ydot)); 14869566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Ystage)); 14879566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Zdot)); 14889566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->Zstage)); 14899566063dSJacob Faibussowitsch PetscCall(VecDuplicate(ts->vec_sol, &ros->vec_sol_prev)); 14909566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &dm)); 14919566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 14929566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1493b39943a6SLisandro Dalcin /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 14949566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 149548a46eb9SPierre Jolivet if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 14969566063dSJacob Faibussowitsch PetscCall(DMTSGetRHSJacobian(dm, &rhsjacobian, NULL)); 1497a3ab5968SHong Zhang if (rhsjacobian == TSComputeRHSJacobianConstant) { 1498a3ab5968SHong Zhang Mat Amat, Pmat; 1499a3ab5968SHong Zhang 1500a3ab5968SHong Zhang /* Set the SNES matrix to be different from the RHS matrix because there is no way to reconstruct shift*M-J */ 15019566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &Amat, &Pmat, NULL, NULL)); 1502a3ab5968SHong Zhang if (Amat && Amat == ts->Arhs) { 1503a3ab5968SHong Zhang if (Amat == Pmat) { 15049566063dSJacob Faibussowitsch PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 15059566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, Amat, Amat, NULL, NULL)); 1506a3ab5968SHong Zhang } else { 15079566063dSJacob Faibussowitsch PetscCall(MatDuplicate(ts->Arhs, MAT_COPY_VALUES, &Amat)); 15089566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, Amat, NULL, NULL, NULL)); 1509a3ab5968SHong Zhang if (Pmat && Pmat == ts->Brhs) { 15109566063dSJacob Faibussowitsch PetscCall(MatDuplicate(ts->Brhs, MAT_COPY_VALUES, &Pmat)); 15119566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, Pmat, NULL, NULL)); 15129566063dSJacob Faibussowitsch PetscCall(MatDestroy(&Pmat)); 1513a3ab5968SHong Zhang } 1514a3ab5968SHong Zhang } 15159566063dSJacob Faibussowitsch PetscCall(MatDestroy(&Amat)); 1516a3ab5968SHong Zhang } 1517a3ab5968SHong Zhang } 15183ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1519e27a552bSJed Brown } 1520e27a552bSJed Brown /*------------------------------------------------------------*/ 1521e27a552bSJed Brown 1522ce78bad3SBarry Smith static PetscErrorCode TSSetFromOptions_RosW(TS ts, PetscOptionItems PetscOptionsObject) 1523d71ae5a4SJacob Faibussowitsch { 152461692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1525b39943a6SLisandro Dalcin SNES snes; 1526e27a552bSJed Brown 1527e27a552bSJed Brown PetscFunctionBegin; 1528d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "RosW ODE solver options"); 1529e27a552bSJed Brown { 153061692a83SJed Brown RosWTableauLink link; 1531e27a552bSJed Brown PetscInt count, choice; 1532e27a552bSJed Brown PetscBool flg; 1533e27a552bSJed Brown const char **namelist; 153461692a83SJed Brown 1535fbccb6d4SPierre Jolivet for (link = RosWTableauList, count = 0; link; link = link->next, count++); 15369566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(count, (char ***)&namelist)); 153761692a83SJed Brown for (link = RosWTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 15389566063dSJacob Faibussowitsch PetscCall(PetscOptionsEList("-ts_rosw_type", "Family of Rosenbrock-W method", "TSRosWSetType", (const char *const *)namelist, count, ros->tableau->name, &choice, &flg)); 15399566063dSJacob Faibussowitsch if (flg) PetscCall(TSRosWSetType(ts, namelist[choice])); 15409566063dSJacob Faibussowitsch PetscCall(PetscFree(namelist)); 154161692a83SJed Brown 15429566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-ts_rosw_recompute_jacobian", "Recompute the Jacobian at each stage", "TSRosWSetRecomputeJacobian", ros->recompute_jacobian, &ros->recompute_jacobian, NULL)); 1543b39943a6SLisandro Dalcin } 1544d0609cedSBarry Smith PetscOptionsHeadEnd(); 154561692a83SJed Brown /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 15469566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 154748a46eb9SPierre Jolivet if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESKSPONLY)); 15483ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1549e27a552bSJed Brown } 1550e27a552bSJed Brown 1551d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSView_RosW(TS ts, PetscViewer viewer) 1552d71ae5a4SJacob Faibussowitsch { 155361692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1554e27a552bSJed Brown PetscBool iascii; 1555e27a552bSJed Brown 1556e27a552bSJed Brown PetscFunctionBegin; 15579566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1558e27a552bSJed Brown if (iascii) { 15599c334d8fSLisandro Dalcin RosWTableau tab = ros->tableau; 156019fd82e9SBarry Smith TSRosWType rostype; 15619c334d8fSLisandro Dalcin char buf[512]; 1562e408995aSJed Brown PetscInt i; 1563e408995aSJed Brown PetscReal abscissa[512]; 15649566063dSJacob Faibussowitsch PetscCall(TSRosWGetType(ts, &rostype)); 15659566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Rosenbrock-W %s\n", rostype)); 15669566063dSJacob Faibussowitsch PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, tab->ASum)); 15679566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A = %s\n", buf)); 1568e408995aSJed Brown for (i = 0; i < tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 15699566063dSJacob Faibussowitsch PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", tab->s, abscissa)); 15709566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa of A+Gamma = %s\n", buf)); 1571e27a552bSJed Brown } 15723ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1573e27a552bSJed Brown } 1574e27a552bSJed Brown 1575d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSLoad_RosW(TS ts, PetscViewer viewer) 1576d71ae5a4SJacob Faibussowitsch { 15779200755eSBarry Smith SNES snes; 15789c334d8fSLisandro Dalcin TSAdapt adapt; 15799200755eSBarry Smith 15809200755eSBarry Smith PetscFunctionBegin; 15819566063dSJacob Faibussowitsch PetscCall(TSGetAdapt(ts, &adapt)); 15829566063dSJacob Faibussowitsch PetscCall(TSAdaptLoad(adapt, viewer)); 15839566063dSJacob Faibussowitsch PetscCall(TSGetSNES(ts, &snes)); 15849566063dSJacob Faibussowitsch PetscCall(SNESLoad(snes, viewer)); 15859200755eSBarry Smith /* function and Jacobian context for SNES when used with TS is always ts object */ 15869566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, NULL, NULL, ts)); 15879566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, ts)); 15883ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 15899200755eSBarry Smith } 15909200755eSBarry Smith 1591cc4c1da9SBarry Smith /*@ 1592bcf0153eSBarry Smith TSRosWSetType - Set the type of Rosenbrock-W, `TSROSW`, scheme 1593e27a552bSJed Brown 159420f4b53cSBarry Smith Logically Collective 1595e27a552bSJed Brown 1596d8d19677SJose E. Roman Input Parameters: 1597e27a552bSJed Brown + ts - timestepping context 1598b92453a8SLisandro Dalcin - roswtype - type of Rosenbrock-W scheme 1599e27a552bSJed Brown 1600020d8f30SJed Brown Level: beginner 1601e27a552bSJed Brown 16021cc06b55SBarry Smith .seealso: [](ch_ts), `TSRosWGetType()`, `TSROSW`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWARK3` 1603e27a552bSJed Brown @*/ 1604d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWSetType(TS ts, TSRosWType roswtype) 1605d71ae5a4SJacob Faibussowitsch { 1606e27a552bSJed Brown PetscFunctionBegin; 1607e27a552bSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 16084f572ea9SToby Isaac PetscAssertPointer(roswtype, 2); 1609cac4c232SBarry Smith PetscTryMethod(ts, "TSRosWSetType_C", (TS, TSRosWType), (ts, roswtype)); 16103ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1611e27a552bSJed Brown } 1612e27a552bSJed Brown 1613cc4c1da9SBarry Smith /*@ 161461692a83SJed Brown TSRosWGetType - Get the type of Rosenbrock-W scheme 1615e27a552bSJed Brown 161620f4b53cSBarry Smith Logically Collective 1617e27a552bSJed Brown 1618e27a552bSJed Brown Input Parameter: 1619e27a552bSJed Brown . ts - timestepping context 1620e27a552bSJed Brown 1621e27a552bSJed Brown Output Parameter: 162261692a83SJed Brown . rostype - type of Rosenbrock-W scheme 1623e27a552bSJed Brown 1624e27a552bSJed Brown Level: intermediate 1625e27a552bSJed Brown 16261cc06b55SBarry Smith .seealso: [](ch_ts), `TSRosWType`, `TSRosWSetType()` 1627e27a552bSJed Brown @*/ 1628d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWGetType(TS ts, TSRosWType *rostype) 1629d71ae5a4SJacob Faibussowitsch { 1630e27a552bSJed Brown PetscFunctionBegin; 1631e27a552bSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1632cac4c232SBarry Smith PetscUseMethod(ts, "TSRosWGetType_C", (TS, TSRosWType *), (ts, rostype)); 16333ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1634e27a552bSJed Brown } 1635e27a552bSJed Brown 1636cc4c1da9SBarry Smith /*@ 163761692a83SJed Brown TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 1638e27a552bSJed Brown 163920f4b53cSBarry Smith Logically Collective 1640e27a552bSJed Brown 1641d8d19677SJose E. Roman Input Parameters: 1642e27a552bSJed Brown + ts - timestepping context 1643bcf0153eSBarry Smith - flg - `PETSC_TRUE` to recompute the Jacobian at each stage 1644e27a552bSJed Brown 1645e27a552bSJed Brown Level: intermediate 1646e27a552bSJed Brown 16471cc06b55SBarry Smith .seealso: [](ch_ts), `TSRosWType`, `TSRosWGetType()` 1648e27a552bSJed Brown @*/ 1649d71ae5a4SJacob Faibussowitsch PetscErrorCode TSRosWSetRecomputeJacobian(TS ts, PetscBool flg) 1650d71ae5a4SJacob Faibussowitsch { 1651e27a552bSJed Brown PetscFunctionBegin; 1652e27a552bSJed Brown PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1653cac4c232SBarry Smith PetscTryMethod(ts, "TSRosWSetRecomputeJacobian_C", (TS, PetscBool), (ts, flg)); 16543ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1655e27a552bSJed Brown } 1656e27a552bSJed Brown 1657d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWGetType_RosW(TS ts, TSRosWType *rostype) 1658d71ae5a4SJacob Faibussowitsch { 165961692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1660e27a552bSJed Brown 1661e27a552bSJed Brown PetscFunctionBegin; 166261692a83SJed Brown *rostype = ros->tableau->name; 16633ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1664e27a552bSJed Brown } 1665ef20d060SBarry Smith 1666d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWSetType_RosW(TS ts, TSRosWType rostype) 1667d71ae5a4SJacob Faibussowitsch { 166861692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1669e27a552bSJed Brown PetscBool match; 167061692a83SJed Brown RosWTableauLink link; 1671e27a552bSJed Brown 1672e27a552bSJed Brown PetscFunctionBegin; 167361692a83SJed Brown if (ros->tableau) { 16749566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(ros->tableau->name, rostype, &match)); 16753ba16761SJacob Faibussowitsch if (match) PetscFunctionReturn(PETSC_SUCCESS); 1676e27a552bSJed Brown } 167761692a83SJed Brown for (link = RosWTableauList; link; link = link->next) { 16789566063dSJacob Faibussowitsch PetscCall(PetscStrcmp(link->tab.name, rostype, &match)); 1679e27a552bSJed Brown if (match) { 16809566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSRosWTableauReset(ts)); 168161692a83SJed Brown ros->tableau = &link->tab; 16829566063dSJacob Faibussowitsch if (ts->setupcalled) PetscCall(TSRosWTableauSetUp(ts)); 16832ffb9264SLisandro Dalcin ts->default_adapt_type = ros->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 16843ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1685e27a552bSJed Brown } 1686e27a552bSJed Brown } 168798921bdaSJacob Faibussowitsch SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rostype); 1688e27a552bSJed Brown } 168961692a83SJed Brown 1690d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts, PetscBool flg) 1691d71ae5a4SJacob Faibussowitsch { 169261692a83SJed Brown TS_RosW *ros = (TS_RosW *)ts->data; 1693e27a552bSJed Brown 1694e27a552bSJed Brown PetscFunctionBegin; 169561692a83SJed Brown ros->recompute_jacobian = flg; 16963ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1697e27a552bSJed Brown } 1698e27a552bSJed Brown 1699d71ae5a4SJacob Faibussowitsch static PetscErrorCode TSDestroy_RosW(TS ts) 1700d71ae5a4SJacob Faibussowitsch { 1701b3a6b972SJed Brown PetscFunctionBegin; 17029566063dSJacob Faibussowitsch PetscCall(TSReset_RosW(ts)); 1703b3a6b972SJed Brown if (ts->dm) { 17049566063dSJacob Faibussowitsch PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRosW, DMRestrictHook_TSRosW, ts)); 17059566063dSJacob Faibussowitsch PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRosW, DMSubDomainRestrictHook_TSRosW, ts)); 1706b3a6b972SJed Brown } 17079566063dSJacob Faibussowitsch PetscCall(PetscFree(ts->data)); 17089566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", NULL)); 17099566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", NULL)); 17109566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", NULL)); 17113ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1712b3a6b972SJed Brown } 1713d5e6173cSPeter Brune 1714e27a552bSJed Brown /* ------------------------------------------------------------ */ 1715e27a552bSJed Brown /*MC 1716020d8f30SJed Brown TSROSW - ODE solver using Rosenbrock-W schemes 1717e27a552bSJed Brown 1718e27a552bSJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 1719e27a552bSJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 1720bcf0153eSBarry Smith of the equation using `TSSetIFunction()` and the non-stiff part with `TSSetRHSFunction()`. 1721bcf0153eSBarry Smith 1722bcf0153eSBarry Smith Level: beginner 1723e27a552bSJed Brown 1724e27a552bSJed Brown Notes: 172561692a83SJed Brown This method currently only works with autonomous ODE and DAE. 172661692a83SJed Brown 1727bcf0153eSBarry Smith Consider trying `TSARKIMEX` if the stiff part is strongly nonlinear. 1728d0685a90SJed Brown 17293d5a8a6aSBarry Smith Since this uses a single linear solve per time-step if you wish to lag the jacobian or preconditioner computation you must use also -snes_lag_jacobian_persists true or -snes_lag_jacobian_preconditioner true 17303d5a8a6aSBarry Smith 1731e94cfbe0SPatrick Sanan Developer Notes: 173261692a83SJed Brown Rosenbrock-W methods are typically specified for autonomous ODE 173361692a83SJed Brown 1734f9c1d6abSBarry Smith $ udot = f(u) 173561692a83SJed Brown 173661692a83SJed Brown by the stage equations 173761692a83SJed Brown 1738f9c1d6abSBarry Smith $ k_i = h f(u_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 173961692a83SJed Brown 174061692a83SJed Brown and step completion formula 174161692a83SJed Brown 1742f9c1d6abSBarry Smith $ u_1 = u_0 + sum_j b_j k_j 174361692a83SJed Brown 1744f9c1d6abSBarry Smith with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(u) 174561692a83SJed Brown and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 174661692a83SJed Brown we define new variables for the stage equations 174761692a83SJed Brown 174861692a83SJed Brown $ y_i = gamma_ij k_j 174961692a83SJed Brown 175061692a83SJed Brown The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 175161692a83SJed Brown 1752b70472e2SJed Brown $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-1} 175361692a83SJed Brown 175461692a83SJed Brown to rewrite the method as 175561692a83SJed Brown 175637fdd005SBarry Smith .vb 175737fdd005SBarry 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 175837fdd005SBarry Smith u_1 = u_0 + sum_j bt_j y_j 175937fdd005SBarry Smith .ve 176061692a83SJed Brown 176161692a83SJed Brown where we have introduced the mass matrix M. Continue by defining 176261692a83SJed Brown 176361692a83SJed Brown $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 176461692a83SJed Brown 176561692a83SJed Brown or, more compactly in tensor notation 176661692a83SJed Brown 176761692a83SJed Brown $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 176861692a83SJed Brown 176961692a83SJed Brown Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 177061692a83SJed Brown stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 177161692a83SJed Brown equation 177261692a83SJed Brown 1773f9c1d6abSBarry Smith $ g(u_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 177461692a83SJed Brown 177561692a83SJed Brown with initial guess y_i = 0. 1776e27a552bSJed Brown 17771cc06b55SBarry Smith .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSRosWSetType()`, `TSRosWRegister()`, `TSROSWTHETA1`, `TSROSWTHETA2`, `TSROSW2M`, `TSROSW2P`, `TSROSWRA3PW`, `TSROSWRA34PW2`, `TSROSWRODAS3`, 1778bcf0153eSBarry Smith `TSROSWSANDU3`, `TSROSWASSP3P3S1C`, `TSROSWLASSP3P4S2C`, `TSROSWLLSSP3P4S2C`, `TSROSWGRK4T`, `TSROSWSHAMP4`, `TSROSWVELDD4`, `TSROSW4L`, `TSType` 1779e27a552bSJed Brown M*/ 1780d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS ts) 1781d71ae5a4SJacob Faibussowitsch { 178261692a83SJed Brown TS_RosW *ros; 1783e27a552bSJed Brown 1784e27a552bSJed Brown PetscFunctionBegin; 17859566063dSJacob Faibussowitsch PetscCall(TSRosWInitializePackage()); 1786e27a552bSJed Brown 1787e27a552bSJed Brown ts->ops->reset = TSReset_RosW; 1788e27a552bSJed Brown ts->ops->destroy = TSDestroy_RosW; 1789e27a552bSJed Brown ts->ops->view = TSView_RosW; 17909200755eSBarry Smith ts->ops->load = TSLoad_RosW; 1791e27a552bSJed Brown ts->ops->setup = TSSetUp_RosW; 1792e27a552bSJed Brown ts->ops->step = TSStep_RosW; 1793e27a552bSJed Brown ts->ops->interpolate = TSInterpolate_RosW; 17941c3436cfSJed Brown ts->ops->evaluatestep = TSEvaluateStep_RosW; 1795e27a552bSJed Brown ts->ops->setfromoptions = TSSetFromOptions_RosW; 1796e27a552bSJed Brown ts->ops->snesfunction = SNESTSFormFunction_RosW; 1797e27a552bSJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1798e27a552bSJed Brown 1799efd4aadfSBarry Smith ts->usessnes = PETSC_TRUE; 1800efd4aadfSBarry Smith 18014dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&ros)); 180261692a83SJed Brown ts->data = (void *)ros; 1803e27a552bSJed Brown 18049566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWGetType_C", TSRosWGetType_RosW)); 18059566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetType_C", TSRosWSetType_RosW)); 18069566063dSJacob Faibussowitsch PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRosWSetRecomputeJacobian_C", TSRosWSetRecomputeJacobian_RosW)); 1807b39943a6SLisandro Dalcin 18089566063dSJacob Faibussowitsch PetscCall(TSRosWSetType(ts, TSRosWDefault)); 18093ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1810e27a552bSJed Brown } 1811