1 /* 2 Code for time stepping with the Runge-Kutta method 3 4 Notes: 5 The general system is written as 6 7 Udot = F(t,U) 8 9 */ 10 11 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 12 #include <petscdm.h> 13 #include <../src/ts/impls/explicit/rk/rk.h> 14 #include <../src/ts/impls/explicit/rk/mrk.h> 15 16 static TSRKType TSRKDefault = TSRK3BS; 17 static PetscBool TSRKRegisterAllCalled; 18 static PetscBool TSRKPackageInitialized; 19 20 static RKTableauLink RKTableauList; 21 22 /*MC 23 TSRK1FE - First order forward Euler scheme. 24 25 This method has one stage. 26 27 Options database: 28 . -ts_rk_type 1fe - use type 1fe 29 30 Level: advanced 31 32 .seealso: TSRK, TSRKType, TSRKSetType() 33 M*/ 34 /*MC 35 TSRK2A - Second order RK scheme (Heun's method). 36 37 This method has two stages. 38 39 Options database: 40 . -ts_rk_type 2a - use type 2a 41 42 Level: advanced 43 44 .seealso: TSRK, TSRKType, TSRKSetType() 45 M*/ 46 /*MC 47 TSRK2B - Second order RK scheme (the midpoint method). 48 49 This method has two stages. 50 51 Options database: 52 . -ts_rk_type 2b - use type 2b 53 54 Level: advanced 55 56 .seealso: TSRK, TSRKType, TSRKSetType() 57 M*/ 58 /*MC 59 TSRK3 - Third order RK scheme. 60 61 This method has three stages. 62 63 Options database: 64 . -ts_rk_type 3 - use type 3 65 66 Level: advanced 67 68 .seealso: TSRK, TSRKType, TSRKSetType() 69 M*/ 70 /*MC 71 TSRK3BS - Third order RK scheme of Bogacki-Shampine with 2nd order embedded method. 72 73 This method has four stages with the First Same As Last (FSAL) property. 74 75 Options database: 76 . -ts_rk_type 3bs - use type 3bs 77 78 Level: advanced 79 80 References: 81 . * - https://doi.org/10.1016/0893-9659(89)90079-7 82 83 .seealso: TSRK, TSRKType, TSRKSetType() 84 M*/ 85 /*MC 86 TSRK4 - Fourth order RK scheme. 87 88 This is the classical Runge-Kutta method with four stages. 89 90 Options database: 91 . -ts_rk_type 4 - use type 4 92 93 Level: advanced 94 95 .seealso: TSRK, TSRKType, TSRKSetType() 96 M*/ 97 /*MC 98 TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method. 99 100 This method has six stages. 101 102 Options database: 103 . -ts_rk_type 5f - use type 5f 104 105 Level: advanced 106 107 .seealso: TSRK, TSRKType, TSRKSetType() 108 M*/ 109 /*MC 110 TSRK5DP - Fifth order Dormand-Prince RK scheme with the 4th order embedded method. 111 112 This method has seven stages with the First Same As Last (FSAL) property. 113 114 Options database: 115 . -ts_rk_type 5dp - use type 5dp 116 117 Level: advanced 118 119 References: 120 . * - https://doi.org/10.1016/0771-050X(80)90013-3 121 122 .seealso: TSRK, TSRKType, TSRKSetType() 123 M*/ 124 /*MC 125 TSRK5BS - Fifth order Bogacki-Shampine RK scheme with 4th order embedded method. 126 127 This method has eight stages with the First Same As Last (FSAL) property. 128 129 Options database: 130 . -ts_rk_type 5bs - use type 5bs 131 132 Level: advanced 133 134 References: 135 . * - https://doi.org/10.1016/0898-1221(96)00141-1 136 137 .seealso: TSRK, TSRKType, TSRKSetType() 138 M*/ 139 /*MC 140 TSRK6VR - Sixth order robust Verner RK scheme with fifth order embedded method. 141 142 This method has nine stages with the First Same As Last (FSAL) property. 143 144 Options database: 145 . -ts_rk_type 6vr - use type 6vr 146 147 Level: advanced 148 149 References: 150 . * - http://people.math.sfu.ca/~jverner/RKV65.IIIXb.Robust.00010102836.081204.CoeffsOnlyRAT 151 152 .seealso: TSRK, TSRKType, TSRKSetType() 153 M*/ 154 /*MC 155 TSRK7VR - Seventh order robust Verner RK scheme with sixth order embedded method. 156 157 This method has ten stages. 158 159 Options database: 160 . -ts_rk_type 7vr - use type 7vr 161 162 Level: advanced 163 164 References: 165 . * - http://people.math.sfu.ca/~jverner/RKV76.IIa.Robust.000027015646.081206.CoeffsOnlyRAT 166 167 .seealso: TSRK, TSRKType, TSRKSetType() 168 M*/ 169 /*MC 170 TSRK8VR - Eigth order robust Verner RK scheme with seventh order embedded method. 171 172 This method has thirteen stages. 173 174 Options database: 175 . -ts_rk_type 8vr - use type 8vr 176 177 Level: advanced 178 179 References: 180 . * - http://people.math.sfu.ca/~jverner/RKV87.IIa.Robust.00000754677.081208.CoeffsOnlyRATandFLOAT 181 182 .seealso: TSRK, TSRKType, TSRKSetType() 183 M*/ 184 185 /*@C 186 TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in TSRK 187 188 Not Collective, but should be called by all processes which will need the schemes to be registered 189 190 Level: advanced 191 192 .seealso: TSRKRegisterDestroy() 193 @*/ 194 PetscErrorCode TSRKRegisterAll(void) 195 { 196 PetscFunctionBegin; 197 if (TSRKRegisterAllCalled) PetscFunctionReturn(0); 198 TSRKRegisterAllCalled = PETSC_TRUE; 199 200 #define RC PetscRealConstant 201 { 202 const PetscReal 203 A[1][1] = {{0}}, 204 b[1] = {RC(1.0)}; 205 PetscCall(TSRKRegister(TSRK1FE,1,1,&A[0][0],b,NULL,NULL,0,NULL)); 206 } 207 { 208 const PetscReal 209 A[2][2] = {{0,0}, 210 {RC(1.0),0}}, 211 b[2] = {RC(0.5),RC(0.5)}, 212 bembed[2] = {RC(1.0),0}; 213 PetscCall(TSRKRegister(TSRK2A,2,2,&A[0][0],b,NULL,bembed,0,NULL)); 214 } 215 { 216 const PetscReal 217 A[2][2] = {{0,0}, 218 {RC(0.5),0}}, 219 b[2] = {0,RC(1.0)}; 220 PetscCall(TSRKRegister(TSRK2B,2,2,&A[0][0],b,NULL,NULL,0,NULL)); 221 } 222 { 223 const PetscReal 224 A[3][3] = {{0,0,0}, 225 {RC(2.0)/RC(3.0),0,0}, 226 {RC(-1.0)/RC(3.0),RC(1.0),0}}, 227 b[3] = {RC(0.25),RC(0.5),RC(0.25)}; 228 PetscCall(TSRKRegister(TSRK3,3,3,&A[0][0],b,NULL,NULL,0,NULL)); 229 } 230 { 231 const PetscReal 232 A[4][4] = {{0,0,0,0}, 233 {RC(1.0)/RC(2.0),0,0,0}, 234 {0,RC(3.0)/RC(4.0),0,0}, 235 {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}}, 236 b[4] = {RC(2.0)/RC(9.0),RC(1.0)/RC(3.0),RC(4.0)/RC(9.0),0}, 237 bembed[4] = {RC(7.0)/RC(24.0),RC(1.0)/RC(4.0),RC(1.0)/RC(3.0),RC(1.0)/RC(8.0)}; 238 PetscCall(TSRKRegister(TSRK3BS,3,4,&A[0][0],b,NULL,bembed,0,NULL)); 239 } 240 { 241 const PetscReal 242 A[4][4] = {{0,0,0,0}, 243 {RC(0.5),0,0,0}, 244 {0,RC(0.5),0,0}, 245 {0,0,RC(1.0),0}}, 246 b[4] = {RC(1.0)/RC(6.0),RC(1.0)/RC(3.0),RC(1.0)/RC(3.0),RC(1.0)/RC(6.0)}; 247 PetscCall(TSRKRegister(TSRK4,4,4,&A[0][0],b,NULL,NULL,0,NULL)); 248 } 249 { 250 const PetscReal 251 A[6][6] = {{0,0,0,0,0,0}, 252 {RC(0.25),0,0,0,0,0}, 253 {RC(3.0)/RC(32.0),RC(9.0)/RC(32.0),0,0,0,0}, 254 {RC(1932.0)/RC(2197.0),RC(-7200.0)/RC(2197.0),RC(7296.0)/RC(2197.0),0,0,0}, 255 {RC(439.0)/RC(216.0),RC(-8.0),RC(3680.0)/RC(513.0),RC(-845.0)/RC(4104.0),0,0}, 256 {RC(-8.0)/RC(27.0),RC(2.0),RC(-3544.0)/RC(2565.0),RC(1859.0)/RC(4104.0),RC(-11.0)/RC(40.0),0}}, 257 b[6] = {RC(16.0)/RC(135.0),0,RC(6656.0)/RC(12825.0),RC(28561.0)/RC(56430.0),RC(-9.0)/RC(50.0),RC(2.0)/RC(55.0)}, 258 bembed[6] = {RC(25.0)/RC(216.0),0,RC(1408.0)/RC(2565.0),RC(2197.0)/RC(4104.0),RC(-1.0)/RC(5.0),0}; 259 PetscCall(TSRKRegister(TSRK5F,5,6,&A[0][0],b,NULL,bembed,0,NULL)); 260 } 261 { 262 const PetscReal 263 A[7][7] = {{0,0,0,0,0,0,0}, 264 {RC(1.0)/RC(5.0),0,0,0,0,0,0}, 265 {RC(3.0)/RC(40.0),RC(9.0)/RC(40.0),0,0,0,0,0}, 266 {RC(44.0)/RC(45.0),RC(-56.0)/RC(15.0),RC(32.0)/RC(9.0),0,0,0,0}, 267 {RC(19372.0)/RC(6561.0),RC(-25360.0)/RC(2187.0),RC(64448.0)/RC(6561.0),RC(-212.0)/RC(729.0),0,0,0}, 268 {RC(9017.0)/RC(3168.0),RC(-355.0)/RC(33.0),RC(46732.0)/RC(5247.0),RC(49.0)/RC(176.0),RC(-5103.0)/RC(18656.0),0,0}, 269 {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}}, 270 b[7] = {RC(35.0)/RC(384.0),0,RC(500.0)/RC(1113.0),RC(125.0)/RC(192.0),RC(-2187.0)/RC(6784.0),RC(11.0)/RC(84.0),0}, 271 bembed[7] = {RC(5179.0)/RC(57600.0),0,RC(7571.0)/RC(16695.0),RC(393.0)/RC(640.0),RC(-92097.0)/RC(339200.0),RC(187.0)/RC(2100.0),RC(1.0)/RC(40.0)}, 272 binterp[7][5] = {{RC(1.0),RC(-4034104133.0)/RC(1410260304.0),RC(105330401.0)/RC(33982176.0),RC(-13107642775.0)/RC(11282082432.0),RC(6542295.0)/RC(470086768.0)}, 273 {0,0,0,0,0}, 274 {0,RC(132343189600.0)/RC(32700410799.0),RC(-833316000.0)/RC(131326951.0),RC(91412856700.0)/RC(32700410799.0),RC(-523383600.0)/RC(10900136933.0)}, 275 {0,RC(-115792950.0)/RC(29380423.0),RC(185270875.0)/RC(16991088.0),RC(-12653452475.0)/RC(1880347072.0),RC(98134425.0)/RC(235043384.0)}, 276 {0,RC(70805911779.0)/RC(24914598704.0),RC(-4531260609.0)/RC(600351776.0),RC(988140236175.0)/RC(199316789632.0),RC(-14307999165.0)/RC(24914598704.0)}, 277 {0,RC(-331320693.0)/RC(205662961.0),RC(31361737.0)/RC(7433601.0),RC(-2426908385.0)/RC(822651844.0),RC(97305120.0)/RC(205662961.0)}, 278 {0,RC(44764047.0)/RC(29380423.0),RC(-1532549.0)/RC(353981.0),RC(90730570.0)/RC(29380423.0),RC(-8293050.0)/RC(29380423.0)}}; 279 PetscCall(TSRKRegister(TSRK5DP,5,7,&A[0][0],b,NULL,bembed,5,binterp[0])); 280 } 281 { 282 const PetscReal 283 A[8][8] = {{0,0,0,0,0,0,0,0}, 284 {RC(1.0)/RC(6.0),0,0,0,0,0,0,0}, 285 {RC(2.0)/RC(27.0),RC(4.0)/RC(27.0),0,0,0,0,0,0}, 286 {RC(183.0)/RC(1372.0),RC(-162.0)/RC(343.0),RC(1053.0)/RC(1372.0),0,0,0,0,0}, 287 {RC(68.0)/RC(297.0),RC(-4.0)/RC(11.0),RC(42.0)/RC(143.0),RC(1960.0)/RC(3861.0),0,0,0,0}, 288 {RC(597.0)/RC(22528.0),RC(81.0)/RC(352.0),RC(63099.0)/RC(585728.0),RC(58653.0)/RC(366080.0),RC(4617.0)/RC(20480.0),0,0,0}, 289 {RC(174197.0)/RC(959244.0),RC(-30942.0)/RC(79937.0),RC(8152137.0)/RC(19744439.0),RC(666106.0)/RC(1039181.0),RC(-29421.0)/RC(29068.0),RC(482048.0)/RC(414219.0),0,0}, 290 {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}}, 291 b[8] = {RC(587.0)/RC(8064.0),0,RC(4440339.0)/RC(15491840.0),RC(24353.0)/RC(124800.0),RC(387.0)/RC(44800.0),RC(2152.0)/RC(5985.0),RC(7267.0)/RC(94080.0),0}, 292 bembed[8] = {RC(2479.0)/RC(34992.0),0,RC(123.0)/RC(416.0),RC(612941.0)/RC(3411720.0),RC(43.0)/RC(1440.0),RC(2272.0)/RC(6561.0),RC(79937.0)/RC(1113912.0),RC(3293.0)/RC(556956.0)}; 293 PetscCall(TSRKRegister(TSRK5BS,5,8,&A[0][0],b,NULL,bembed,0,NULL)); 294 } 295 { 296 const PetscReal 297 A[9][9] = {{0,0,0,0,0,0,0,0,0}, 298 {RC(1.8000000000000000000000000000000000000000e-01),0,0,0,0,0,0,0,0}, 299 {RC(8.9506172839506172839506172839506172839506e-02),RC(7.7160493827160493827160493827160493827160e-02),0,0,0,0,0,0,0}, 300 {RC(6.2500000000000000000000000000000000000000e-02),0,RC(1.8750000000000000000000000000000000000000e-01),0,0,0,0,0,0}, 301 {RC(3.1651600000000000000000000000000000000000e-01),0,RC(-1.0449480000000000000000000000000000000000e+00),RC(1.2584320000000000000000000000000000000000e+00),0,0,0,0,0}, 302 {RC(2.7232612736485626257225065566674305502508e-01),0,RC(-8.2513360323886639676113360323886639676113e-01),RC(1.0480917678812415654520917678812415654521e+00),RC(1.0471570799276856873679117969088177628396e-01),0,0,0,0}, 303 {RC(-1.6699418599716514314329607278961797333198e-01),0,RC(6.3170850202429149797570850202429149797571e-01),RC(1.7461044552773876082146758838488161796432e-01),RC(-1.0665356459086066122525194734018680677781e+00),RC(1.2272108843537414965986394557823129251701e+00),0,0,0}, 304 {RC(3.6423751686909581646423751686909581646424e-01),0,RC(-2.0404858299595141700404858299595141700405e-01),RC(-3.4883737816068643136312309244640071707741e-01),RC(3.2619323032856867443333608747142581729048e+00),RC(-2.7551020408163265306122448979591836734694e+00),RC(6.8181818181818181818181818181818181818182e-01),0,0}, 305 {RC(7.6388888888888888888888888888888888888889e-02),0,0,RC(3.6940836940836940836940836940836940836941e-01),0,RC(2.4801587301587301587301587301587301587302e-01),RC(2.3674242424242424242424242424242424242424e-01),RC(6.9444444444444444444444444444444444444444e-02),0}}, 306 b[9] = {RC(7.6388888888888888888888888888888888888889e-02),0,0,RC(3.6940836940836940836940836940836940836941e-01),0,RC(2.4801587301587301587301587301587301587302e-01),RC(2.3674242424242424242424242424242424242424e-01),RC(6.9444444444444444444444444444444444444444e-02),0}, 307 bembed[9] = {RC(5.8700209643605870020964360587002096436059e-02),0,0,RC(4.8072562358276643990929705215419501133787e-01),RC(-8.5341242076919085578832094861228313083563e-01),RC(1.2046485260770975056689342403628117913832e+00),0,RC(-5.9242373072160306202859394348756050883710e-02),RC(1.6858043453788134639198468985703028256220e-01)}; 308 PetscCall(TSRKRegister(TSRK6VR,6,9,&A[0][0],b,NULL,bembed,0,NULL)); 309 } 310 { 311 const PetscReal 312 A[10][10] = {{0,0,0,0,0,0,0,0,0,0}, 313 {RC(5.0000000000000000000000000000000000000000e-03),0,0,0,0,0,0,0,0,0}, 314 {RC(-1.0767901234567901234567901234567901234568e+00),RC(1.1856790123456790123456790123456790123457e+00),0,0,0,0,0,0,0,0}, 315 {RC(4.0833333333333333333333333333333333333333e-02),0,RC(1.2250000000000000000000000000000000000000e-01),0,0,0,0,0,0,0}, 316 {RC(6.3607142857142857142857142857142857142857e-01),0,RC(-2.4444642857142857142857142857142857142857e+00),RC(2.2633928571428571428571428571428571428571e+00),0,0,0,0,0,0}, 317 {RC(-2.5351211079349245229256383554660215487207e+00),0,RC(1.0299374654449267920438514460756024913612e+01),RC(-7.9513032885990579949493217458266876536482e+00),RC(7.9301148923100592201226014271115261823800e-01),0,0,0,0,0}, 318 {RC(1.0018765812524632961969196583094999808207e+00),0,RC(-4.1665712824423798331313938005470971453189e+00),RC(3.8343432929128642412552665218251378665197e+00),RC(-5.0233333560710847547464330228611765612403e-01),RC(6.6768474388416077115385092269857695410259e-01),0,0,0,0}, 319 {RC(2.7255018354630767130333963819175005717348e+01),0,RC(-4.2004617278410638355318645443909295369611e+01),RC(-1.0535713126619489917921081600546526103722e+01),RC(8.0495536711411937147983652158926826634202e+01),RC(-6.7343882271790513468549075963212975640927e+01),RC(1.3048657610777937463471187029566964762710e+01),0,0,0}, 320 {RC(-3.0397378057114965146943658658755763226883e+00),0,RC(1.0138161410329801111857946190709700150441e+01),RC(-6.4293056748647215721462825629555298064437e+00),RC(-1.5864371483408276587115312853798610579467e+00),RC(1.8921781841968424410864308909131353365021e+00),RC(1.9699335407608869061292360163336442838006e-02),RC(5.4416989827933235465102724247952572977903e-03),0,0}, 321 {RC(-1.4449518916777735137351003179355712360517e+00),0,RC(8.0318913859955919224117033223019560435041e+00),RC(-7.5831741663401346820798883023671588604984e+00),RC(3.5816169353190074211247685442452878696855e+00),RC(-2.4369722632199529411183809065693752383733e+00),RC(8.5158999992326179339689766032486142173390e-01),0,0,0}}, 322 b[10] = {RC(4.7425837833706756083569172717574534698932e-02),0,0,RC(2.5622361659370562659961727458274623448160e-01),RC(2.6951376833074206619473817258075952886764e-01),RC(1.2686622409092782845989138364739173247882e-01),RC(2.4887225942060071622046449427647492767292e-01),RC(3.0744837408200631335304388479099184768645e-03),RC(4.8023809989496943308189063347143123323209e-02),0}, 323 bembed[10] = {RC(4.7485247699299631037531273805727961552268e-02),0,0,RC(2.5599412588690633297154918245905393870497e-01),RC(2.7058478081067688722530891099268135732387e-01),RC(1.2505618684425992913638822323746917920448e-01),RC(2.5204468723743860507184043820197442562182e-01),0,0,RC(4.8834971521418614557381971303093137592592e-02)}; 324 PetscCall(TSRKRegister(TSRK7VR,7,10,&A[0][0],b,NULL,bembed,0,NULL)); 325 } 326 { 327 const PetscReal 328 A[13][13] = {{0,0,0,0,0,0,0,0,0,0,0,0,0}, 329 {RC(2.5000000000000000000000000000000000000000e-01),0,0,0,0,0,0,0,0,0,0,0,0}, 330 {RC(8.7400846504915232052686327594877411977046e-02),RC(2.5487604938654321753087950620345685135815e-02),0,0,0,0,0,0,0,0,0,0,0}, 331 {RC(4.2333169291338582677165354330708661417323e-02),0,RC(1.2699950787401574803149606299212598425197e-01),0,0,0,0,0,0,0,0,0,0}, 332 {RC(4.2609505888742261494881445237572274090942e-01),0,RC(-1.5987952846591523265427733230657181117089e+00),RC(1.5967002257717297115939588706899953707994e+00),0,0,0,0,0,0,0,0,0}, 333 {RC(5.0719337296713929515090618138513639239329e-02),0,0,RC(2.5433377264600407582754714408877778031369e-01),RC(2.0394689005728199465736223777270858044698e-01),0,0,0,0,0,0,0,0}, 334 {RC(-2.9000374717523110970388379285425896124091e-01),0,0,RC(1.3441873910260789889438681109414337003184e+00),RC(-2.8647779433614427309611103827036562829470e+00),RC(2.6775942995105948517211260646164815438695e+00),0,0,0,0,0,0,0}, 335 {RC(9.8535011337993546469740402980727014284756e-02),0,0,0,RC(2.2192680630751384842024036498197387903583e-01),RC(-1.8140622911806994312690338288073952457474e-01),RC(1.0944411472562548236922614918038631254153e-02),0,0,0,0,0,0}, 336 {RC(3.8711052545731144679444618165166373405645e-01),0,0,RC(-1.4424454974855277571256745553077927767173e+00),RC(2.9053981890699509317691346449233848441744e+00),RC(-1.8537710696301059290843332675811978025183e+00),RC(1.4003648098728154269497325109771241479223e-01),RC(5.7273940811495816575746774624447706488753e-01),0,0,0,0,0}, 337 {RC(-1.6124403444439308100630016197913480595436e-01),0,0,RC(-1.7339602957358984083578404473962567894901e-01),RC(-1.3012892814065147406016812745172492529744e+00),RC(1.1379503751738617308558792131431003472124e+00),RC(-3.1747649663966880106923521138043024698980e-02),RC(9.3351293824933666439811064486056884856590e-01),RC(-8.3786318334733852703300855629616433201504e-02),0,0,0,0}, 338 {RC(-1.9199444881589533281510804651483576073142e-02),0,0,RC(2.7330857265264284907942326254016124275617e-01),RC(-6.7534973206944372919691611210942380856240e-01),RC(3.4151849813846016071738489974728382711981e-01),RC(-6.7950064803375772478920516198524629391910e-02),RC(9.6591752247623878884265586491216376509746e-02),RC(1.3253082511182101180721038466545389951226e-01),RC(3.6854959360386113446906329951531666812946e-01),0,0,0}, 339 {RC(6.0918774036452898676888412111588817784584e-01),0,0,RC(-2.2725690858980016768999800931413088399719e+00),RC(4.7578983426940290068155255881914785497547e+00),RC(-5.5161067066927584824294689667844248244842e+00),RC(2.9005963696801192709095818565946174378180e-01),RC(5.6914239633590368229109858454801849145630e-01),RC(7.9267957603321670271339916205893327579951e-01),RC(1.5473720453288822894126190771849898232047e-01),RC(1.6149708956621816247083215106334544434974e+00),0,0}, 340 {RC(8.8735762208534719663211694051981022704884e-01),0,0,RC(-2.9754597821085367558513632804709301581977e+00),RC(5.6007170094881630597990392548350098923829e+00),RC(-5.9156074505366744680014930189941657351840e+00),RC(2.2029689156134927016879142540807638331238e-01),RC(1.0155097824462216666143271340902996997549e-01),RC(1.1514345647386055909780397752125850553556e+00),RC(1.9297101665271239396134361900805843653065e+00),0,0,0}}, 341 b[13] = {RC(4.4729564666695714203015840429049382466467e-02),0,0,0,0,RC(1.5691033527708199813368698010726645409175e-01),RC(1.8460973408151637740702451873526277892035e-01),RC(2.2516380602086991042479419400350721970920e-01),RC(1.4794615651970234687005179885449141753736e-01),RC(7.6055542444955825269798361910336491012732e-02),RC(1.2277290235018619610824346315921437388535e-01),RC(4.1811958638991631583384842800871882376786e-02),0}, 342 bembed[13] = {RC(4.5847111400495925878664730122010282095875e-02),0,0,0,0,RC(2.6231891404152387437443356584845803392392e-01),RC(1.9169372337852611904485738635688429008025e-01),RC(2.1709172327902618330978407422906448568196e-01),RC(1.2738189624833706796803169450656737867900e-01),RC(1.1510530385365326258240515750043192148894e-01),0,0,RC(4.0561327798437566841823391436583608050053e-02)}; 343 PetscCall(TSRKRegister(TSRK8VR,8,13,&A[0][0],b,NULL,bembed,0,NULL)); 344 } 345 #undef RC 346 PetscFunctionReturn(0); 347 } 348 349 /*@C 350 TSRKRegisterDestroy - Frees the list of schemes that were registered by TSRKRegister(). 351 352 Not Collective 353 354 Level: advanced 355 356 .seealso: TSRKRegister(), TSRKRegisterAll() 357 @*/ 358 PetscErrorCode TSRKRegisterDestroy(void) 359 { 360 RKTableauLink link; 361 362 PetscFunctionBegin; 363 while ((link = RKTableauList)) { 364 RKTableau t = &link->tab; 365 RKTableauList = link->next; 366 PetscCall(PetscFree3(t->A,t->b,t->c)); 367 PetscCall(PetscFree(t->bembed)); 368 PetscCall(PetscFree(t->binterp)); 369 PetscCall(PetscFree(t->name)); 370 PetscCall(PetscFree(link)); 371 } 372 TSRKRegisterAllCalled = PETSC_FALSE; 373 PetscFunctionReturn(0); 374 } 375 376 /*@C 377 TSRKInitializePackage - This function initializes everything in the TSRK package. It is called 378 from TSInitializePackage(). 379 380 Level: developer 381 382 .seealso: PetscInitialize() 383 @*/ 384 PetscErrorCode TSRKInitializePackage(void) 385 { 386 PetscFunctionBegin; 387 if (TSRKPackageInitialized) PetscFunctionReturn(0); 388 TSRKPackageInitialized = PETSC_TRUE; 389 PetscCall(TSRKRegisterAll()); 390 PetscCall(PetscRegisterFinalize(TSRKFinalizePackage)); 391 PetscFunctionReturn(0); 392 } 393 394 /*@C 395 TSRKFinalizePackage - This function destroys everything in the TSRK package. It is 396 called from PetscFinalize(). 397 398 Level: developer 399 400 .seealso: PetscFinalize() 401 @*/ 402 PetscErrorCode TSRKFinalizePackage(void) 403 { 404 PetscFunctionBegin; 405 TSRKPackageInitialized = PETSC_FALSE; 406 PetscCall(TSRKRegisterDestroy()); 407 PetscFunctionReturn(0); 408 } 409 410 /*@C 411 TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 412 413 Not Collective, but the same schemes should be registered on all processes on which they will be used 414 415 Input Parameters: 416 + name - identifier for method 417 . order - approximation order of method 418 . s - number of stages, this is the dimension of the matrices below 419 . A - stage coefficients (dimension s*s, row-major) 420 . b - step completion table (dimension s; NULL to use last row of A) 421 . c - abscissa (dimension s; NULL to use row sums of A) 422 . bembed - completion table for embedded method (dimension s; NULL if not available) 423 . p - Order of the interpolation scheme, equal to the number of columns of binterp 424 - binterp - Coefficients of the interpolation formula (dimension s*p; NULL to reuse b with p=1) 425 426 Notes: 427 Several RK methods are provided, this function is only needed to create new methods. 428 429 Level: advanced 430 431 .seealso: TSRK 432 @*/ 433 PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s, 434 const PetscReal A[],const PetscReal b[],const PetscReal c[], 435 const PetscReal bembed[],PetscInt p,const PetscReal binterp[]) 436 { 437 RKTableauLink link; 438 RKTableau t; 439 PetscInt i,j; 440 441 PetscFunctionBegin; 442 PetscValidCharPointer(name,1); 443 PetscValidRealPointer(A,4); 444 if (b) PetscValidRealPointer(b,5); 445 if (c) PetscValidRealPointer(c,6); 446 if (bembed) PetscValidRealPointer(bembed,7); 447 if (binterp || p > 1) PetscValidRealPointer(binterp,9); 448 449 PetscCall(TSRKInitializePackage()); 450 PetscCall(PetscNew(&link)); 451 t = &link->tab; 452 453 PetscCall(PetscStrallocpy(name,&t->name)); 454 t->order = order; 455 t->s = s; 456 PetscCall(PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c)); 457 PetscCall(PetscArraycpy(t->A,A,s*s)); 458 if (b) PetscCall(PetscArraycpy(t->b,b,s)); 459 else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i]; 460 if (c) PetscCall(PetscArraycpy(t->c,c,s)); 461 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 462 t->FSAL = PETSC_TRUE; 463 for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE; 464 465 if (bembed) { 466 PetscCall(PetscMalloc1(s,&t->bembed)); 467 PetscCall(PetscArraycpy(t->bembed,bembed,s)); 468 } 469 470 if (!binterp) { p = 1; binterp = t->b; } 471 t->p = p; 472 PetscCall(PetscMalloc1(s*p,&t->binterp)); 473 PetscCall(PetscArraycpy(t->binterp,binterp,s*p)); 474 475 link->next = RKTableauList; 476 RKTableauList = link; 477 PetscFunctionReturn(0); 478 } 479 480 PetscErrorCode TSRKGetTableau_RK(TS ts, PetscInt *s, const PetscReal **A, const PetscReal **b, const PetscReal **c, const PetscReal **bembed, 481 PetscInt *p, const PetscReal **binterp, PetscBool *FSAL) 482 { 483 TS_RK *rk = (TS_RK*)ts->data; 484 RKTableau tab = rk->tableau; 485 486 PetscFunctionBegin; 487 if (s) *s = tab->s; 488 if (A) *A = tab->A; 489 if (b) *b = tab->b; 490 if (c) *c = tab->c; 491 if (bembed) *bembed = tab->bembed; 492 if (p) *p = tab->p; 493 if (binterp) *binterp = tab->binterp; 494 if (FSAL) *FSAL = tab->FSAL; 495 PetscFunctionReturn(0); 496 } 497 498 /*@C 499 TSRKGetTableau - Get info on the RK tableau 500 501 Not Collective 502 503 Input Parameter: 504 . ts - timestepping context 505 506 Output Parameters: 507 + s - number of stages, this is the dimension of the matrices below 508 . A - stage coefficients (dimension s*s, row-major) 509 . b - step completion table (dimension s) 510 . c - abscissa (dimension s) 511 . bembed - completion table for embedded method (dimension s; NULL if not available) 512 . p - Order of the interpolation scheme, equal to the number of columns of binterp 513 . binterp - Coefficients of the interpolation formula (dimension s*p) 514 - FSAL - wheather or not the scheme has the First Same As Last property 515 516 Level: developer 517 518 .seealso: TSRK 519 @*/ 520 PetscErrorCode TSRKGetTableau(TS ts, PetscInt *s, const PetscReal **A, const PetscReal **b, const PetscReal **c, const PetscReal **bembed, 521 PetscInt *p, const PetscReal **binterp, PetscBool *FSAL) 522 { 523 PetscFunctionBegin; 524 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 525 PetscUseMethod(ts,"TSRKGetTableau_C",(TS,PetscInt*,const PetscReal**,const PetscReal**,const PetscReal**,const PetscReal**, 526 PetscInt*,const PetscReal**,PetscBool*),(ts,s,A,b,c,bembed,p,binterp,FSAL)); 527 PetscFunctionReturn(0); 528 } 529 530 /* 531 This is for single-step RK method 532 The step completion formula is 533 534 x1 = x0 + h b^T YdotRHS 535 536 This function can be called before or after ts->vec_sol has been updated. 537 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 538 We can write 539 540 x1e = x0 + h be^T YdotRHS 541 = x1 - h b^T YdotRHS + h be^T YdotRHS 542 = x1 + h (be - b)^T YdotRHS 543 544 so we can evaluate the method with different order even after the step has been optimistically completed. 545 */ 546 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 547 { 548 TS_RK *rk = (TS_RK*)ts->data; 549 RKTableau tab = rk->tableau; 550 PetscScalar *w = rk->work; 551 PetscReal h; 552 PetscInt s = tab->s,j; 553 554 PetscFunctionBegin; 555 switch (rk->status) { 556 case TS_STEP_INCOMPLETE: 557 case TS_STEP_PENDING: 558 h = ts->time_step; break; 559 case TS_STEP_COMPLETE: 560 h = ts->ptime - ts->ptime_prev; break; 561 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 562 } 563 if (order == tab->order) { 564 if (rk->status == TS_STEP_INCOMPLETE) { 565 PetscCall(VecCopy(ts->vec_sol,X)); 566 for (j=0; j<s; j++) w[j] = h*tab->b[j]/rk->dtratio; 567 PetscCall(VecMAXPY(X,s,w,rk->YdotRHS)); 568 } else PetscCall(VecCopy(ts->vec_sol,X)); 569 PetscFunctionReturn(0); 570 } else if (order == tab->order-1) { 571 if (!tab->bembed) goto unavailable; 572 if (rk->status == TS_STEP_INCOMPLETE) { /*Complete with the embedded method (be)*/ 573 PetscCall(VecCopy(ts->vec_sol,X)); 574 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 575 PetscCall(VecMAXPY(X,s,w,rk->YdotRHS)); 576 } else { /*Rollback and re-complete using (be-b) */ 577 PetscCall(VecCopy(ts->vec_sol,X)); 578 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 579 PetscCall(VecMAXPY(X,s,w,rk->YdotRHS)); 580 } 581 if (done) *done = PETSC_TRUE; 582 PetscFunctionReturn(0); 583 } 584 unavailable: 585 if (done) *done = PETSC_FALSE; 586 else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"RK '%s' of order %D cannot evaluate step at order %D. Consider using -ts_adapt_type none or a different method that has an embedded estimate.",tab->name,tab->order,order); 587 PetscFunctionReturn(0); 588 } 589 590 static PetscErrorCode TSForwardCostIntegral_RK(TS ts) 591 { 592 TS_RK *rk = (TS_RK*)ts->data; 593 TS quadts = ts->quadraturets; 594 RKTableau tab = rk->tableau; 595 const PetscInt s = tab->s; 596 const PetscReal *b = tab->b,*c = tab->c; 597 Vec *Y = rk->Y; 598 PetscInt i; 599 600 PetscFunctionBegin; 601 /* No need to backup quadts->vec_sol since it can be reverted in TSRollBack_RK */ 602 for (i=s-1; i>=0; i--) { 603 /* Evolve quadrature TS solution to compute integrals */ 604 PetscCall(TSComputeRHSFunction(quadts,rk->ptime+rk->time_step*c[i],Y[i],ts->vec_costintegrand)); 605 PetscCall(VecAXPY(quadts->vec_sol,rk->time_step*b[i],ts->vec_costintegrand)); 606 } 607 PetscFunctionReturn(0); 608 } 609 610 static PetscErrorCode TSAdjointCostIntegral_RK(TS ts) 611 { 612 TS_RK *rk = (TS_RK*)ts->data; 613 RKTableau tab = rk->tableau; 614 TS quadts = ts->quadraturets; 615 const PetscInt s = tab->s; 616 const PetscReal *b = tab->b,*c = tab->c; 617 Vec *Y = rk->Y; 618 PetscInt i; 619 620 PetscFunctionBegin; 621 for (i=s-1; i>=0; i--) { 622 /* Evolve quadrature TS solution to compute integrals */ 623 PetscCall(TSComputeRHSFunction(quadts,ts->ptime+ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand)); 624 PetscCall(VecAXPY(quadts->vec_sol,-ts->time_step*b[i],ts->vec_costintegrand)); 625 } 626 PetscFunctionReturn(0); 627 } 628 629 static PetscErrorCode TSRollBack_RK(TS ts) 630 { 631 TS_RK *rk = (TS_RK*)ts->data; 632 TS quadts = ts->quadraturets; 633 RKTableau tab = rk->tableau; 634 const PetscInt s = tab->s; 635 const PetscReal *b = tab->b,*c = tab->c; 636 PetscScalar *w = rk->work; 637 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 638 PetscInt j; 639 PetscReal h; 640 641 PetscFunctionBegin; 642 switch (rk->status) { 643 case TS_STEP_INCOMPLETE: 644 case TS_STEP_PENDING: 645 h = ts->time_step; break; 646 case TS_STEP_COMPLETE: 647 h = ts->ptime - ts->ptime_prev; break; 648 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 649 } 650 for (j=0; j<s; j++) w[j] = -h*b[j]; 651 PetscCall(VecMAXPY(ts->vec_sol,s,w,YdotRHS)); 652 if (quadts && ts->costintegralfwd) { 653 for (j=0; j<s; j++) { 654 /* Revert the quadrature TS solution */ 655 PetscCall(TSComputeRHSFunction(quadts,rk->ptime+h*c[j],Y[j],ts->vec_costintegrand)); 656 PetscCall(VecAXPY(quadts->vec_sol,-h*b[j],ts->vec_costintegrand)); 657 } 658 } 659 PetscFunctionReturn(0); 660 } 661 662 static PetscErrorCode TSForwardStep_RK(TS ts) 663 { 664 TS_RK *rk = (TS_RK*)ts->data; 665 RKTableau tab = rk->tableau; 666 Mat J,*MatsFwdSensipTemp = rk->MatsFwdSensipTemp; 667 const PetscInt s = tab->s; 668 const PetscReal *A = tab->A,*c = tab->c,*b = tab->b; 669 Vec *Y = rk->Y; 670 PetscInt i,j; 671 PetscReal stage_time,h = ts->time_step; 672 PetscBool zero; 673 674 PetscFunctionBegin; 675 PetscCall(MatCopy(ts->mat_sensip,rk->MatFwdSensip0,SAME_NONZERO_PATTERN)); 676 PetscCall(TSGetRHSJacobian(ts,&J,NULL,NULL,NULL)); 677 678 for (i=0; i<s; i++) { 679 stage_time = ts->ptime + h*c[i]; 680 zero = PETSC_FALSE; 681 if (b[i] == 0 && i == s-1) zero = PETSC_TRUE; 682 /* TLM Stage values */ 683 if (!i) { 684 PetscCall(MatCopy(ts->mat_sensip,rk->MatsFwdStageSensip[i],SAME_NONZERO_PATTERN)); 685 } else if (!zero) { 686 PetscCall(MatZeroEntries(rk->MatsFwdStageSensip[i])); 687 for (j=0; j<i; j++) { 688 PetscCall(MatAXPY(rk->MatsFwdStageSensip[i],h*A[i*s+j],MatsFwdSensipTemp[j],SAME_NONZERO_PATTERN)); 689 } 690 PetscCall(MatAXPY(rk->MatsFwdStageSensip[i],1.,ts->mat_sensip,SAME_NONZERO_PATTERN)); 691 } else { 692 PetscCall(MatZeroEntries(rk->MatsFwdStageSensip[i])); 693 } 694 695 PetscCall(TSComputeRHSJacobian(ts,stage_time,Y[i],J,J)); 696 PetscCall(MatMatMult(J,rk->MatsFwdStageSensip[i],MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatsFwdSensipTemp[i])); 697 if (ts->Jacprhs) { 698 PetscCall(TSComputeRHSJacobianP(ts,stage_time,Y[i],ts->Jacprhs)); /* get f_p */ 699 if (ts->vecs_sensi2p) { /* TLM used for 2nd-order adjoint */ 700 PetscScalar *xarr; 701 PetscCall(MatDenseGetColumn(MatsFwdSensipTemp[i],0,&xarr)); 702 PetscCall(VecPlaceArray(rk->VecDeltaFwdSensipCol,xarr)); 703 PetscCall(MatMultAdd(ts->Jacprhs,ts->vec_dir,rk->VecDeltaFwdSensipCol,rk->VecDeltaFwdSensipCol)); 704 PetscCall(VecResetArray(rk->VecDeltaFwdSensipCol)); 705 PetscCall(MatDenseRestoreColumn(MatsFwdSensipTemp[i],&xarr)); 706 } else { 707 PetscCall(MatAXPY(MatsFwdSensipTemp[i],1.,ts->Jacprhs,SUBSET_NONZERO_PATTERN)); 708 } 709 } 710 } 711 712 for (i=0; i<s; i++) { 713 PetscCall(MatAXPY(ts->mat_sensip,h*b[i],rk->MatsFwdSensipTemp[i],SAME_NONZERO_PATTERN)); 714 } 715 rk->status = TS_STEP_COMPLETE; 716 PetscFunctionReturn(0); 717 } 718 719 static PetscErrorCode TSForwardGetStages_RK(TS ts,PetscInt *ns,Mat **stagesensip) 720 { 721 TS_RK *rk = (TS_RK*)ts->data; 722 RKTableau tab = rk->tableau; 723 724 PetscFunctionBegin; 725 if (ns) *ns = tab->s; 726 if (stagesensip) *stagesensip = rk->MatsFwdStageSensip; 727 PetscFunctionReturn(0); 728 } 729 730 static PetscErrorCode TSForwardSetUp_RK(TS ts) 731 { 732 TS_RK *rk = (TS_RK*)ts->data; 733 RKTableau tab = rk->tableau; 734 PetscInt i; 735 736 PetscFunctionBegin; 737 /* backup sensitivity results for roll-backs */ 738 PetscCall(MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatFwdSensip0)); 739 740 PetscCall(PetscMalloc1(tab->s,&rk->MatsFwdStageSensip)); 741 PetscCall(PetscMalloc1(tab->s,&rk->MatsFwdSensipTemp)); 742 for (i=0; i<tab->s; i++) { 743 PetscCall(MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdStageSensip[i])); 744 PetscCall(MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdSensipTemp[i])); 745 } 746 PetscCall(VecDuplicate(ts->vec_sol,&rk->VecDeltaFwdSensipCol)); 747 PetscFunctionReturn(0); 748 } 749 750 static PetscErrorCode TSForwardReset_RK(TS ts) 751 { 752 TS_RK *rk = (TS_RK*)ts->data; 753 RKTableau tab = rk->tableau; 754 PetscInt i; 755 756 PetscFunctionBegin; 757 PetscCall(MatDestroy(&rk->MatFwdSensip0)); 758 if (rk->MatsFwdStageSensip) { 759 for (i=0; i<tab->s; i++) { 760 PetscCall(MatDestroy(&rk->MatsFwdStageSensip[i])); 761 } 762 PetscCall(PetscFree(rk->MatsFwdStageSensip)); 763 } 764 if (rk->MatsFwdSensipTemp) { 765 for (i=0; i<tab->s; i++) { 766 PetscCall(MatDestroy(&rk->MatsFwdSensipTemp[i])); 767 } 768 PetscCall(PetscFree(rk->MatsFwdSensipTemp)); 769 } 770 PetscCall(VecDestroy(&rk->VecDeltaFwdSensipCol)); 771 PetscFunctionReturn(0); 772 } 773 774 static PetscErrorCode TSStep_RK(TS ts) 775 { 776 TS_RK *rk = (TS_RK*)ts->data; 777 RKTableau tab = rk->tableau; 778 const PetscInt s = tab->s; 779 const PetscReal *A = tab->A,*c = tab->c; 780 PetscScalar *w = rk->work; 781 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 782 PetscBool FSAL = tab->FSAL; 783 TSAdapt adapt; 784 PetscInt i,j; 785 PetscInt rejections = 0; 786 PetscBool stageok,accept = PETSC_TRUE; 787 PetscReal next_time_step = ts->time_step; 788 789 PetscFunctionBegin; 790 if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE; 791 if (FSAL) PetscCall(VecCopy(YdotRHS[s-1],YdotRHS[0])); 792 793 rk->status = TS_STEP_INCOMPLETE; 794 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 795 PetscReal t = ts->ptime; 796 PetscReal h = ts->time_step; 797 for (i=0; i<s; i++) { 798 rk->stage_time = t + h*c[i]; 799 PetscCall(TSPreStage(ts,rk->stage_time)); 800 PetscCall(VecCopy(ts->vec_sol,Y[i])); 801 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 802 PetscCall(VecMAXPY(Y[i],i,w,YdotRHS)); 803 PetscCall(TSPostStage(ts,rk->stage_time,i,Y)); 804 PetscCall(TSGetAdapt(ts,&adapt)); 805 PetscCall(TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok)); 806 if (!stageok) goto reject_step; 807 if (FSAL && !i) continue; 808 PetscCall(TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i])); 809 } 810 811 rk->status = TS_STEP_INCOMPLETE; 812 PetscCall(TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL)); 813 rk->status = TS_STEP_PENDING; 814 PetscCall(TSGetAdapt(ts,&adapt)); 815 PetscCall(TSAdaptCandidatesClear(adapt)); 816 PetscCall(TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE)); 817 PetscCall(TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept)); 818 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 819 if (!accept) { /* Roll back the current step */ 820 PetscCall(TSRollBack_RK(ts)); 821 ts->time_step = next_time_step; 822 goto reject_step; 823 } 824 825 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */ 826 rk->ptime = ts->ptime; 827 rk->time_step = ts->time_step; 828 } 829 830 ts->ptime += ts->time_step; 831 ts->time_step = next_time_step; 832 break; 833 834 reject_step: 835 ts->reject++; accept = PETSC_FALSE; 836 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 837 ts->reason = TS_DIVERGED_STEP_REJECTED; 838 PetscCall(PetscInfo(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections)); 839 } 840 } 841 PetscFunctionReturn(0); 842 } 843 844 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 845 { 846 TS_RK *rk = (TS_RK*)ts->data; 847 RKTableau tab = rk->tableau; 848 PetscInt s = tab->s; 849 850 PetscFunctionBegin; 851 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 852 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam)); 853 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecsSensiTemp)); 854 if (ts->vecs_sensip) { 855 PetscCall(VecDuplicate(ts->vecs_sensip[0],&rk->VecDeltaMu)); 856 } 857 if (ts->vecs_sensi2) { 858 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam2)); 859 PetscCall(VecDuplicateVecs(ts->vecs_sensi2[0],ts->numcost,&rk->VecsSensi2Temp)); 860 } 861 if (ts->vecs_sensi2p) { 862 PetscCall(VecDuplicate(ts->vecs_sensi2p[0],&rk->VecDeltaMu2)); 863 } 864 PetscFunctionReturn(0); 865 } 866 867 /* 868 Assumptions: 869 - TSStep_RK() always evaluates the step with b, not bembed. 870 */ 871 static PetscErrorCode TSAdjointStep_RK(TS ts) 872 { 873 TS_RK *rk = (TS_RK*)ts->data; 874 TS quadts = ts->quadraturets; 875 RKTableau tab = rk->tableau; 876 Mat J,Jpre,Jquad; 877 const PetscInt s = tab->s; 878 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 879 PetscScalar *w = rk->work,*xarr; 880 Vec *Y = rk->Y,*VecsDeltaLam = rk->VecsDeltaLam,VecDeltaMu = rk->VecDeltaMu,*VecsSensiTemp = rk->VecsSensiTemp; 881 Vec *VecsDeltaLam2 = rk->VecsDeltaLam2,VecDeltaMu2 = rk->VecDeltaMu2,*VecsSensi2Temp = rk->VecsSensi2Temp; 882 Vec VecDRDUTransCol = ts->vec_drdu_col,VecDRDPTransCol = ts->vec_drdp_col; 883 PetscInt i,j,nadj; 884 PetscReal t = ts->ptime; 885 PetscReal h = ts->time_step; 886 887 PetscFunctionBegin; 888 rk->status = TS_STEP_INCOMPLETE; 889 890 PetscCall(TSGetRHSJacobian(ts,&J,&Jpre,NULL,NULL)); 891 if (quadts) { 892 PetscCall(TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL)); 893 } 894 for (i=s-1; i>=0; i--) { 895 if (tab->FSAL && i == s-1) { 896 /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/ 897 continue; 898 } 899 rk->stage_time = t + h*(1.0-c[i]); 900 PetscCall(TSComputeSNESJacobian(ts,Y[i],J,Jpre)); 901 if (quadts) { 902 PetscCall(TSComputeRHSJacobian(quadts,rk->stage_time,Y[i],Jquad,Jquad)); /* get r_u^T */ 903 } 904 if (ts->vecs_sensip) { 905 PetscCall(TSComputeRHSJacobianP(ts,rk->stage_time,Y[i],ts->Jacprhs)); /* get f_p */ 906 if (quadts) { 907 PetscCall(TSComputeRHSJacobianP(quadts,rk->stage_time,Y[i],quadts->Jacprhs)); /* get f_p for the quadrature */ 908 } 909 } 910 911 if (b[i]) { 912 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]/b[i]; /* coefficients for computing VecsSensiTemp */ 913 } else { 914 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]; /* coefficients for computing VecsSensiTemp */ 915 } 916 917 for (nadj=0; nadj<ts->numcost; nadj++) { 918 /* Stage values of lambda */ 919 if (b[i]) { 920 /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */ 921 PetscCall(VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj])); /* VecDeltaLam is an vec array of size s by numcost */ 922 PetscCall(VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1])); 923 PetscCall(MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i])); /* VecsSensiTemp will be reused by 2nd-order adjoint */ 924 PetscCall(VecScale(VecsDeltaLam[nadj*s+i],-h*b[i])); 925 if (quadts) { 926 PetscCall(MatDenseGetColumn(Jquad,nadj,&xarr)); 927 PetscCall(VecPlaceArray(VecDRDUTransCol,xarr)); 928 PetscCall(VecAXPY(VecsDeltaLam[nadj*s+i],-h*b[i],VecDRDUTransCol)); 929 PetscCall(VecResetArray(VecDRDUTransCol)); 930 PetscCall(MatDenseRestoreColumn(Jquad,&xarr)); 931 } 932 } else { 933 /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */ 934 PetscCall(VecSet(VecsSensiTemp[nadj],0)); 935 PetscCall(VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1])); 936 PetscCall(MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i])); 937 PetscCall(VecScale(VecsDeltaLam[nadj*s+i],-h)); 938 } 939 940 /* Stage values of mu */ 941 if (ts->vecs_sensip) { 942 if (b[i]) { 943 PetscCall(MatMultTranspose(ts->Jacprhs,VecsSensiTemp[nadj],VecDeltaMu)); 944 PetscCall(VecScale(VecDeltaMu,-h*b[i])); 945 if (quadts) { 946 PetscCall(MatDenseGetColumn(quadts->Jacprhs,nadj,&xarr)); 947 PetscCall(VecPlaceArray(VecDRDPTransCol,xarr)); 948 PetscCall(VecAXPY(VecDeltaMu,-h*b[i],VecDRDPTransCol)); 949 PetscCall(VecResetArray(VecDRDPTransCol)); 950 PetscCall(MatDenseRestoreColumn(quadts->Jacprhs,&xarr)); 951 } 952 } else { 953 PetscCall(VecScale(VecDeltaMu,-h)); 954 } 955 PetscCall(VecAXPY(ts->vecs_sensip[nadj],1.,VecDeltaMu)); /* update sensip for each stage */ 956 } 957 } 958 959 if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */ 960 /* Get w1 at t_{n+1} from TLM matrix */ 961 PetscCall(MatDenseGetColumn(rk->MatsFwdStageSensip[i],0,&xarr)); 962 PetscCall(VecPlaceArray(ts->vec_sensip_col,xarr)); 963 /* lambda_s^T F_UU w_1 */ 964 PetscCall(TSComputeRHSHessianProductFunctionUU(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_guu)); 965 if (quadts) { 966 /* R_UU w_1 */ 967 PetscCall(TSComputeRHSHessianProductFunctionUU(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_guu)); 968 } 969 if (ts->vecs_sensip) { 970 /* lambda_s^T F_UP w_2 */ 971 PetscCall(TSComputeRHSHessianProductFunctionUP(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gup)); 972 if (quadts) { 973 /* R_UP w_2 */ 974 PetscCall(TSComputeRHSHessianProductFunctionUP(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gup)); 975 } 976 } 977 if (ts->vecs_sensi2p) { 978 /* lambda_s^T F_PU w_1 */ 979 PetscCall(TSComputeRHSHessianProductFunctionPU(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_gpu)); 980 /* lambda_s^T F_PP w_2 */ 981 PetscCall(TSComputeRHSHessianProductFunctionPP(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gpp)); 982 if (b[i] && quadts) { 983 /* R_PU w_1 */ 984 PetscCall(TSComputeRHSHessianProductFunctionPU(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gpu)); 985 /* R_PP w_2 */ 986 PetscCall(TSComputeRHSHessianProductFunctionPP(quadts,rk->stage_time,Y[i],NULL,ts->vec_dir,ts->vecs_gpp)); 987 } 988 } 989 PetscCall(VecResetArray(ts->vec_sensip_col)); 990 PetscCall(MatDenseRestoreColumn(rk->MatsFwdStageSensip[i],&xarr)); 991 992 for (nadj=0; nadj<ts->numcost; nadj++) { 993 /* Stage values of lambda */ 994 if (b[i]) { 995 /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */ 996 PetscCall(VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj])); 997 PetscCall(VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1])); 998 PetscCall(MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i])); 999 PetscCall(VecScale(VecsDeltaLam2[nadj*s+i],-h*b[i])); 1000 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_guu[nadj])); 1001 if (ts->vecs_sensip) { 1002 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_gup[nadj])); 1003 } 1004 } else { 1005 /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */ 1006 PetscCall(VecSet(VecsDeltaLam2[nadj*s+i],0)); 1007 PetscCall(VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1])); 1008 PetscCall(MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i])); 1009 PetscCall(VecScale(VecsDeltaLam2[nadj*s+i],-h)); 1010 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_guu[nadj])); 1011 if (ts->vecs_sensip) { 1012 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_gup[nadj])); 1013 } 1014 } 1015 if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */ 1016 PetscCall(MatMultTranspose(ts->Jacprhs,VecsSensi2Temp[nadj],VecDeltaMu2)); 1017 if (b[i]) { 1018 PetscCall(VecScale(VecDeltaMu2,-h*b[i])); 1019 PetscCall(VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpu[nadj])); 1020 PetscCall(VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpp[nadj])); 1021 } else { 1022 PetscCall(VecScale(VecDeltaMu2,-h)); 1023 PetscCall(VecAXPY(VecDeltaMu2,-h,ts->vecs_gpu[nadj])); 1024 PetscCall(VecAXPY(VecDeltaMu2,-h,ts->vecs_gpp[nadj])); 1025 } 1026 PetscCall(VecAXPY(ts->vecs_sensi2p[nadj],1,VecDeltaMu2)); /* update sensi2p for each stage */ 1027 } 1028 } 1029 } 1030 } 1031 1032 for (j=0; j<s; j++) w[j] = 1.0; 1033 for (nadj=0; nadj<ts->numcost; nadj++) { /* no need to do this for mu's */ 1034 PetscCall(VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecsDeltaLam[nadj*s])); 1035 if (ts->vecs_sensi2) { 1036 PetscCall(VecMAXPY(ts->vecs_sensi2[nadj],s,w,&VecsDeltaLam2[nadj*s])); 1037 } 1038 } 1039 rk->status = TS_STEP_COMPLETE; 1040 PetscFunctionReturn(0); 1041 } 1042 1043 static PetscErrorCode TSAdjointReset_RK(TS ts) 1044 { 1045 TS_RK *rk = (TS_RK*)ts->data; 1046 RKTableau tab = rk->tableau; 1047 1048 PetscFunctionBegin; 1049 PetscCall(VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam)); 1050 PetscCall(VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp)); 1051 PetscCall(VecDestroy(&rk->VecDeltaMu)); 1052 PetscCall(VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam2)); 1053 PetscCall(VecDestroy(&rk->VecDeltaMu2)); 1054 PetscCall(VecDestroyVecs(ts->numcost,&rk->VecsSensi2Temp)); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 1059 { 1060 TS_RK *rk = (TS_RK*)ts->data; 1061 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 1062 PetscReal h; 1063 PetscReal tt,t; 1064 PetscScalar *b; 1065 const PetscReal *B = rk->tableau->binterp; 1066 1067 PetscFunctionBegin; 1068 PetscCheck(B,PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 1069 1070 switch (rk->status) { 1071 case TS_STEP_INCOMPLETE: 1072 case TS_STEP_PENDING: 1073 h = ts->time_step; 1074 t = (itime - ts->ptime)/h; 1075 break; 1076 case TS_STEP_COMPLETE: 1077 h = ts->ptime - ts->ptime_prev; 1078 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1079 break; 1080 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1081 } 1082 PetscCall(PetscMalloc1(s,&b)); 1083 for (i=0; i<s; i++) b[i] = 0; 1084 for (j=0,tt=t; j<p; j++,tt*=t) { 1085 for (i=0; i<s; i++) { 1086 b[i] += h * B[i*p+j] * tt; 1087 } 1088 } 1089 PetscCall(VecCopy(rk->Y[0],X)); 1090 PetscCall(VecMAXPY(X,s,b,rk->YdotRHS)); 1091 PetscCall(PetscFree(b)); 1092 PetscFunctionReturn(0); 1093 } 1094 1095 /*------------------------------------------------------------*/ 1096 1097 static PetscErrorCode TSRKTableauReset(TS ts) 1098 { 1099 TS_RK *rk = (TS_RK*)ts->data; 1100 RKTableau tab = rk->tableau; 1101 1102 PetscFunctionBegin; 1103 if (!tab) PetscFunctionReturn(0); 1104 PetscCall(PetscFree(rk->work)); 1105 PetscCall(VecDestroyVecs(tab->s,&rk->Y)); 1106 PetscCall(VecDestroyVecs(tab->s,&rk->YdotRHS)); 1107 PetscFunctionReturn(0); 1108 } 1109 1110 static PetscErrorCode TSReset_RK(TS ts) 1111 { 1112 PetscFunctionBegin; 1113 PetscCall(TSRKTableauReset(ts)); 1114 if (ts->use_splitrhsfunction) { 1115 PetscTryMethod(ts,"TSReset_RK_MultirateSplit_C",(TS),(ts)); 1116 } else { 1117 PetscTryMethod(ts,"TSReset_RK_MultirateNonsplit_C",(TS),(ts)); 1118 } 1119 PetscFunctionReturn(0); 1120 } 1121 1122 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 1123 { 1124 PetscFunctionBegin; 1125 PetscFunctionReturn(0); 1126 } 1127 1128 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1129 { 1130 PetscFunctionBegin; 1131 PetscFunctionReturn(0); 1132 } 1133 1134 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 1135 { 1136 PetscFunctionBegin; 1137 PetscFunctionReturn(0); 1138 } 1139 1140 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1141 { 1142 PetscFunctionBegin; 1143 PetscFunctionReturn(0); 1144 } 1145 1146 static PetscErrorCode TSRKTableauSetUp(TS ts) 1147 { 1148 TS_RK *rk = (TS_RK*)ts->data; 1149 RKTableau tab = rk->tableau; 1150 1151 PetscFunctionBegin; 1152 PetscCall(PetscMalloc1(tab->s,&rk->work)); 1153 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y)); 1154 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS)); 1155 PetscFunctionReturn(0); 1156 } 1157 1158 static PetscErrorCode TSSetUp_RK(TS ts) 1159 { 1160 TS quadts = ts->quadraturets; 1161 DM dm; 1162 1163 PetscFunctionBegin; 1164 PetscCall(TSCheckImplicitTerm(ts)); 1165 PetscCall(TSRKTableauSetUp(ts)); 1166 if (quadts && ts->costintegralfwd) { 1167 Mat Jquad; 1168 PetscCall(TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL)); 1169 } 1170 PetscCall(TSGetDM(ts,&dm)); 1171 PetscCall(DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts)); 1172 PetscCall(DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts)); 1173 if (ts->use_splitrhsfunction) { 1174 PetscTryMethod(ts,"TSSetUp_RK_MultirateSplit_C",(TS),(ts)); 1175 } else { 1176 PetscTryMethod(ts,"TSSetUp_RK_MultirateNonsplit_C",(TS),(ts)); 1177 } 1178 PetscFunctionReturn(0); 1179 } 1180 1181 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 1182 { 1183 TS_RK *rk = (TS_RK*)ts->data; 1184 PetscErrorCode ierr; 1185 1186 PetscFunctionBegin; 1187 PetscCall(PetscOptionsHead(PetscOptionsObject,"RK ODE solver options")); 1188 { 1189 RKTableauLink link; 1190 PetscInt count,choice; 1191 PetscBool flg,use_multirate = PETSC_FALSE; 1192 const char **namelist; 1193 1194 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 1195 PetscCall(PetscMalloc1(count,(char***)&namelist)); 1196 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1197 PetscCall(PetscOptionsBool("-ts_rk_multirate","Use interpolation-based multirate RK method","TSRKSetMultirate",rk->use_multirate,&use_multirate,&flg)); 1198 if (flg) { 1199 PetscCall(TSRKSetMultirate(ts,use_multirate)); 1200 } 1201 PetscCall(PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg)); 1202 if (flg) PetscCall(TSRKSetType(ts,namelist[choice])); 1203 PetscCall(PetscFree(namelist)); 1204 } 1205 PetscCall(PetscOptionsTail()); 1206 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)ts),NULL,"Multirate methods options","");PetscCall(ierr); 1207 PetscCall(PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL)); 1208 ierr = PetscOptionsEnd();PetscCall(ierr); 1209 PetscFunctionReturn(0); 1210 } 1211 1212 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 1213 { 1214 TS_RK *rk = (TS_RK*)ts->data; 1215 PetscBool iascii; 1216 1217 PetscFunctionBegin; 1218 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1219 if (iascii) { 1220 RKTableau tab = rk->tableau; 1221 TSRKType rktype; 1222 const PetscReal *c; 1223 PetscInt s; 1224 char buf[512]; 1225 PetscBool FSAL; 1226 1227 PetscCall(TSRKGetType(ts,&rktype)); 1228 PetscCall(TSRKGetTableau(ts,&s,NULL,NULL,&c,NULL,NULL,NULL,&FSAL)); 1229 PetscCall(PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype)); 1230 PetscCall(PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order)); 1231 PetscCall(PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",FSAL ? "yes" : "no")); 1232 PetscCall(PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",s,c)); 1233 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf)); 1234 } 1235 PetscFunctionReturn(0); 1236 } 1237 1238 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 1239 { 1240 TSAdapt adapt; 1241 1242 PetscFunctionBegin; 1243 PetscCall(TSGetAdapt(ts,&adapt)); 1244 PetscCall(TSAdaptLoad(adapt,viewer)); 1245 PetscFunctionReturn(0); 1246 } 1247 1248 /*@ 1249 TSRKGetOrder - Get the order of RK scheme 1250 1251 Not collective 1252 1253 Input Parameter: 1254 . ts - timestepping context 1255 1256 Output Parameter: 1257 . order - order of RK-scheme 1258 1259 Level: intermediate 1260 1261 .seealso: TSRKGetType() 1262 @*/ 1263 PetscErrorCode TSRKGetOrder(TS ts,PetscInt *order) 1264 { 1265 PetscFunctionBegin; 1266 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1267 PetscValidIntPointer(order,2); 1268 PetscUseMethod(ts,"TSRKGetOrder_C",(TS,PetscInt*),(ts,order)); 1269 PetscFunctionReturn(0); 1270 } 1271 1272 /*@C 1273 TSRKSetType - Set the type of RK scheme 1274 1275 Logically collective 1276 1277 Input Parameters: 1278 + ts - timestepping context 1279 - rktype - type of RK-scheme 1280 1281 Options Database: 1282 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1283 1284 Level: intermediate 1285 1286 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK2B, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS, TSRK6VR, TSRK7VR, TSRK8VR 1287 @*/ 1288 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1289 { 1290 PetscFunctionBegin; 1291 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1292 PetscValidCharPointer(rktype,2); 1293 PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype)); 1294 PetscFunctionReturn(0); 1295 } 1296 1297 /*@C 1298 TSRKGetType - Get the type of RK scheme 1299 1300 Not collective 1301 1302 Input Parameter: 1303 . ts - timestepping context 1304 1305 Output Parameter: 1306 . rktype - type of RK-scheme 1307 1308 Level: intermediate 1309 1310 .seealso: TSRKSetType() 1311 @*/ 1312 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1313 { 1314 PetscFunctionBegin; 1315 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1316 PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype)); 1317 PetscFunctionReturn(0); 1318 } 1319 1320 static PetscErrorCode TSRKGetOrder_RK(TS ts,PetscInt *order) 1321 { 1322 TS_RK *rk = (TS_RK*)ts->data; 1323 1324 PetscFunctionBegin; 1325 *order = rk->tableau->order; 1326 PetscFunctionReturn(0); 1327 } 1328 1329 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1330 { 1331 TS_RK *rk = (TS_RK*)ts->data; 1332 1333 PetscFunctionBegin; 1334 *rktype = rk->tableau->name; 1335 PetscFunctionReturn(0); 1336 } 1337 1338 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1339 { 1340 TS_RK *rk = (TS_RK*)ts->data; 1341 PetscBool match; 1342 RKTableauLink link; 1343 1344 PetscFunctionBegin; 1345 if (rk->tableau) { 1346 PetscCall(PetscStrcmp(rk->tableau->name,rktype,&match)); 1347 if (match) PetscFunctionReturn(0); 1348 } 1349 for (link = RKTableauList; link; link=link->next) { 1350 PetscCall(PetscStrcmp(link->tab.name,rktype,&match)); 1351 if (match) { 1352 if (ts->setupcalled) PetscCall(TSRKTableauReset(ts)); 1353 rk->tableau = &link->tab; 1354 if (ts->setupcalled) PetscCall(TSRKTableauSetUp(ts)); 1355 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1356 PetscFunctionReturn(0); 1357 } 1358 } 1359 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1360 } 1361 1362 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1363 { 1364 TS_RK *rk = (TS_RK*)ts->data; 1365 1366 PetscFunctionBegin; 1367 if (ns) *ns = rk->tableau->s; 1368 if (Y) *Y = rk->Y; 1369 PetscFunctionReturn(0); 1370 } 1371 1372 static PetscErrorCode TSDestroy_RK(TS ts) 1373 { 1374 PetscFunctionBegin; 1375 PetscCall(TSReset_RK(ts)); 1376 if (ts->dm) { 1377 PetscCall(DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts)); 1378 PetscCall(DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts)); 1379 } 1380 PetscCall(PetscFree(ts->data)); 1381 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetOrder_C",NULL)); 1382 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL)); 1383 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL)); 1384 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetTableau_C",NULL)); 1385 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",NULL)); 1386 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",NULL)); 1387 PetscFunctionReturn(0); 1388 } 1389 1390 /* 1391 This defines the nonlinear equation that is to be solved with SNES 1392 We do not need to solve the equation; we just use SNES to approximate the Jacobian 1393 */ 1394 static PetscErrorCode SNESTSFormFunction_RK(SNES snes,Vec x,Vec y,TS ts) 1395 { 1396 TS_RK *rk = (TS_RK*)ts->data; 1397 DM dm,dmsave; 1398 1399 PetscFunctionBegin; 1400 PetscCall(SNESGetDM(snes,&dm)); 1401 /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */ 1402 dmsave = ts->dm; 1403 ts->dm = dm; 1404 PetscCall(TSComputeRHSFunction(ts,rk->stage_time,x,y)); 1405 ts->dm = dmsave; 1406 PetscFunctionReturn(0); 1407 } 1408 1409 static PetscErrorCode SNESTSFormJacobian_RK(SNES snes,Vec x,Mat A,Mat B,TS ts) 1410 { 1411 TS_RK *rk = (TS_RK*)ts->data; 1412 DM dm,dmsave; 1413 1414 PetscFunctionBegin; 1415 PetscCall(SNESGetDM(snes,&dm)); 1416 dmsave = ts->dm; 1417 ts->dm = dm; 1418 PetscCall(TSComputeRHSJacobian(ts,rk->stage_time,x,A,B)); 1419 ts->dm = dmsave; 1420 PetscFunctionReturn(0); 1421 } 1422 1423 /*@C 1424 TSRKSetMultirate - Use the interpolation-based multirate RK method 1425 1426 Logically collective 1427 1428 Input Parameters: 1429 + ts - timestepping context 1430 - use_multirate - PETSC_TRUE enables the multirate RK method, sets the basic method to be RK2A and sets the ratio between slow stepsize and fast stepsize to be 2 1431 1432 Options Database: 1433 . -ts_rk_multirate - <true,false> 1434 1435 Notes: 1436 The multirate method requires interpolation. The default interpolation works for 1st- and 2nd- order RK, but not for high-order RKs except TSRK5DP which comes with the interpolation coeffcients (binterp). 1437 1438 Level: intermediate 1439 1440 .seealso: TSRKGetMultirate() 1441 @*/ 1442 PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate) 1443 { 1444 PetscFunctionBegin; 1445 PetscTryMethod(ts,"TSRKSetMultirate_C",(TS,PetscBool),(ts,use_multirate)); 1446 PetscFunctionReturn(0); 1447 } 1448 1449 /*@C 1450 TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method 1451 1452 Not collective 1453 1454 Input Parameter: 1455 . ts - timestepping context 1456 1457 Output Parameter: 1458 . use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise 1459 1460 Level: intermediate 1461 1462 .seealso: TSRKSetMultirate() 1463 @*/ 1464 PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate) 1465 { 1466 PetscFunctionBegin; 1467 PetscUseMethod(ts,"TSRKGetMultirate_C",(TS,PetscBool*),(ts,use_multirate)); 1468 PetscFunctionReturn(0); 1469 } 1470 1471 /*MC 1472 TSRK - ODE and DAE solver using Runge-Kutta schemes 1473 1474 The user should provide the right hand side of the equation 1475 using TSSetRHSFunction(). 1476 1477 Notes: 1478 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1479 1480 Level: beginner 1481 1482 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRK2D, TTSRK2E, TSRK3, 1483 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirate(), TSRKGetMultirate() 1484 1485 M*/ 1486 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1487 { 1488 TS_RK *rk; 1489 1490 PetscFunctionBegin; 1491 PetscCall(TSRKInitializePackage()); 1492 1493 ts->ops->reset = TSReset_RK; 1494 ts->ops->destroy = TSDestroy_RK; 1495 ts->ops->view = TSView_RK; 1496 ts->ops->load = TSLoad_RK; 1497 ts->ops->setup = TSSetUp_RK; 1498 ts->ops->interpolate = TSInterpolate_RK; 1499 ts->ops->step = TSStep_RK; 1500 ts->ops->evaluatestep = TSEvaluateStep_RK; 1501 ts->ops->rollback = TSRollBack_RK; 1502 ts->ops->setfromoptions = TSSetFromOptions_RK; 1503 ts->ops->getstages = TSGetStages_RK; 1504 1505 ts->ops->snesfunction = SNESTSFormFunction_RK; 1506 ts->ops->snesjacobian = SNESTSFormJacobian_RK; 1507 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1508 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1509 ts->ops->adjointstep = TSAdjointStep_RK; 1510 ts->ops->adjointreset = TSAdjointReset_RK; 1511 1512 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1513 ts->ops->forwardsetup = TSForwardSetUp_RK; 1514 ts->ops->forwardreset = TSForwardReset_RK; 1515 ts->ops->forwardstep = TSForwardStep_RK; 1516 ts->ops->forwardgetstages= TSForwardGetStages_RK; 1517 1518 PetscCall(PetscNewLog(ts,&rk)); 1519 ts->data = (void*)rk; 1520 1521 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetOrder_C",TSRKGetOrder_RK)); 1522 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK)); 1523 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK)); 1524 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetTableau_C",TSRKGetTableau_RK)); 1525 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",TSRKSetMultirate_RK)); 1526 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",TSRKGetMultirate_RK)); 1527 1528 PetscCall(TSRKSetType(ts,TSRKDefault)); 1529 rk->dtratio = 1; 1530 PetscFunctionReturn(0); 1531 } 1532