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 PetscErrorCode ierr; 524 525 PetscFunctionBegin; 526 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 527 ierr = PetscUseMethod(ts,"TSRKGetTableau_C",(TS,PetscInt*,const PetscReal**,const PetscReal**,const PetscReal**,const PetscReal**, 528 PetscInt*,const PetscReal**,PetscBool*),(ts,s,A,b,c,bembed,p,binterp,FSAL));PetscCall(ierr); 529 PetscFunctionReturn(0); 530 } 531 532 /* 533 This is for single-step RK method 534 The step completion formula is 535 536 x1 = x0 + h b^T YdotRHS 537 538 This function can be called before or after ts->vec_sol has been updated. 539 Suppose we have a completion formula (b) and an embedded formula (be) of different order. 540 We can write 541 542 x1e = x0 + h be^T YdotRHS 543 = x1 - h b^T YdotRHS + h be^T YdotRHS 544 = x1 + h (be - b)^T YdotRHS 545 546 so we can evaluate the method with different order even after the step has been optimistically completed. 547 */ 548 static PetscErrorCode TSEvaluateStep_RK(TS ts,PetscInt order,Vec X,PetscBool *done) 549 { 550 TS_RK *rk = (TS_RK*)ts->data; 551 RKTableau tab = rk->tableau; 552 PetscScalar *w = rk->work; 553 PetscReal h; 554 PetscInt s = tab->s,j; 555 556 PetscFunctionBegin; 557 switch (rk->status) { 558 case TS_STEP_INCOMPLETE: 559 case TS_STEP_PENDING: 560 h = ts->time_step; break; 561 case TS_STEP_COMPLETE: 562 h = ts->ptime - ts->ptime_prev; break; 563 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 564 } 565 if (order == tab->order) { 566 if (rk->status == TS_STEP_INCOMPLETE) { 567 PetscCall(VecCopy(ts->vec_sol,X)); 568 for (j=0; j<s; j++) w[j] = h*tab->b[j]/rk->dtratio; 569 PetscCall(VecMAXPY(X,s,w,rk->YdotRHS)); 570 } else PetscCall(VecCopy(ts->vec_sol,X)); 571 PetscFunctionReturn(0); 572 } else if (order == tab->order-1) { 573 if (!tab->bembed) goto unavailable; 574 if (rk->status == TS_STEP_INCOMPLETE) { /*Complete with the embedded method (be)*/ 575 PetscCall(VecCopy(ts->vec_sol,X)); 576 for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 577 PetscCall(VecMAXPY(X,s,w,rk->YdotRHS)); 578 } else { /*Rollback and re-complete using (be-b) */ 579 PetscCall(VecCopy(ts->vec_sol,X)); 580 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 581 PetscCall(VecMAXPY(X,s,w,rk->YdotRHS)); 582 } 583 if (done) *done = PETSC_TRUE; 584 PetscFunctionReturn(0); 585 } 586 unavailable: 587 if (done) *done = PETSC_FALSE; 588 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); 589 PetscFunctionReturn(0); 590 } 591 592 static PetscErrorCode TSForwardCostIntegral_RK(TS ts) 593 { 594 TS_RK *rk = (TS_RK*)ts->data; 595 TS quadts = ts->quadraturets; 596 RKTableau tab = rk->tableau; 597 const PetscInt s = tab->s; 598 const PetscReal *b = tab->b,*c = tab->c; 599 Vec *Y = rk->Y; 600 PetscInt i; 601 602 PetscFunctionBegin; 603 /* No need to backup quadts->vec_sol since it can be reverted in TSRollBack_RK */ 604 for (i=s-1; i>=0; i--) { 605 /* Evolve quadrature TS solution to compute integrals */ 606 PetscCall(TSComputeRHSFunction(quadts,rk->ptime+rk->time_step*c[i],Y[i],ts->vec_costintegrand)); 607 PetscCall(VecAXPY(quadts->vec_sol,rk->time_step*b[i],ts->vec_costintegrand)); 608 } 609 PetscFunctionReturn(0); 610 } 611 612 static PetscErrorCode TSAdjointCostIntegral_RK(TS ts) 613 { 614 TS_RK *rk = (TS_RK*)ts->data; 615 RKTableau tab = rk->tableau; 616 TS quadts = ts->quadraturets; 617 const PetscInt s = tab->s; 618 const PetscReal *b = tab->b,*c = tab->c; 619 Vec *Y = rk->Y; 620 PetscInt i; 621 622 PetscFunctionBegin; 623 for (i=s-1; i>=0; i--) { 624 /* Evolve quadrature TS solution to compute integrals */ 625 PetscCall(TSComputeRHSFunction(quadts,ts->ptime+ts->time_step*(1.0-c[i]),Y[i],ts->vec_costintegrand)); 626 PetscCall(VecAXPY(quadts->vec_sol,-ts->time_step*b[i],ts->vec_costintegrand)); 627 } 628 PetscFunctionReturn(0); 629 } 630 631 static PetscErrorCode TSRollBack_RK(TS ts) 632 { 633 TS_RK *rk = (TS_RK*)ts->data; 634 TS quadts = ts->quadraturets; 635 RKTableau tab = rk->tableau; 636 const PetscInt s = tab->s; 637 const PetscReal *b = tab->b,*c = tab->c; 638 PetscScalar *w = rk->work; 639 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 640 PetscInt j; 641 PetscReal h; 642 643 PetscFunctionBegin; 644 switch (rk->status) { 645 case TS_STEP_INCOMPLETE: 646 case TS_STEP_PENDING: 647 h = ts->time_step; break; 648 case TS_STEP_COMPLETE: 649 h = ts->ptime - ts->ptime_prev; break; 650 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 651 } 652 for (j=0; j<s; j++) w[j] = -h*b[j]; 653 PetscCall(VecMAXPY(ts->vec_sol,s,w,YdotRHS)); 654 if (quadts && ts->costintegralfwd) { 655 for (j=0; j<s; j++) { 656 /* Revert the quadrature TS solution */ 657 PetscCall(TSComputeRHSFunction(quadts,rk->ptime+h*c[j],Y[j],ts->vec_costintegrand)); 658 PetscCall(VecAXPY(quadts->vec_sol,-h*b[j],ts->vec_costintegrand)); 659 } 660 } 661 PetscFunctionReturn(0); 662 } 663 664 static PetscErrorCode TSForwardStep_RK(TS ts) 665 { 666 TS_RK *rk = (TS_RK*)ts->data; 667 RKTableau tab = rk->tableau; 668 Mat J,*MatsFwdSensipTemp = rk->MatsFwdSensipTemp; 669 const PetscInt s = tab->s; 670 const PetscReal *A = tab->A,*c = tab->c,*b = tab->b; 671 Vec *Y = rk->Y; 672 PetscInt i,j; 673 PetscReal stage_time,h = ts->time_step; 674 PetscBool zero; 675 676 PetscFunctionBegin; 677 PetscCall(MatCopy(ts->mat_sensip,rk->MatFwdSensip0,SAME_NONZERO_PATTERN)); 678 PetscCall(TSGetRHSJacobian(ts,&J,NULL,NULL,NULL)); 679 680 for (i=0; i<s; i++) { 681 stage_time = ts->ptime + h*c[i]; 682 zero = PETSC_FALSE; 683 if (b[i] == 0 && i == s-1) zero = PETSC_TRUE; 684 /* TLM Stage values */ 685 if (!i) { 686 PetscCall(MatCopy(ts->mat_sensip,rk->MatsFwdStageSensip[i],SAME_NONZERO_PATTERN)); 687 } else if (!zero) { 688 PetscCall(MatZeroEntries(rk->MatsFwdStageSensip[i])); 689 for (j=0; j<i; j++) { 690 PetscCall(MatAXPY(rk->MatsFwdStageSensip[i],h*A[i*s+j],MatsFwdSensipTemp[j],SAME_NONZERO_PATTERN)); 691 } 692 PetscCall(MatAXPY(rk->MatsFwdStageSensip[i],1.,ts->mat_sensip,SAME_NONZERO_PATTERN)); 693 } else { 694 PetscCall(MatZeroEntries(rk->MatsFwdStageSensip[i])); 695 } 696 697 PetscCall(TSComputeRHSJacobian(ts,stage_time,Y[i],J,J)); 698 PetscCall(MatMatMult(J,rk->MatsFwdStageSensip[i],MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatsFwdSensipTemp[i])); 699 if (ts->Jacprhs) { 700 PetscCall(TSComputeRHSJacobianP(ts,stage_time,Y[i],ts->Jacprhs)); /* get f_p */ 701 if (ts->vecs_sensi2p) { /* TLM used for 2nd-order adjoint */ 702 PetscScalar *xarr; 703 PetscCall(MatDenseGetColumn(MatsFwdSensipTemp[i],0,&xarr)); 704 PetscCall(VecPlaceArray(rk->VecDeltaFwdSensipCol,xarr)); 705 PetscCall(MatMultAdd(ts->Jacprhs,ts->vec_dir,rk->VecDeltaFwdSensipCol,rk->VecDeltaFwdSensipCol)); 706 PetscCall(VecResetArray(rk->VecDeltaFwdSensipCol)); 707 PetscCall(MatDenseRestoreColumn(MatsFwdSensipTemp[i],&xarr)); 708 } else { 709 PetscCall(MatAXPY(MatsFwdSensipTemp[i],1.,ts->Jacprhs,SUBSET_NONZERO_PATTERN)); 710 } 711 } 712 } 713 714 for (i=0; i<s; i++) { 715 PetscCall(MatAXPY(ts->mat_sensip,h*b[i],rk->MatsFwdSensipTemp[i],SAME_NONZERO_PATTERN)); 716 } 717 rk->status = TS_STEP_COMPLETE; 718 PetscFunctionReturn(0); 719 } 720 721 static PetscErrorCode TSForwardGetStages_RK(TS ts,PetscInt *ns,Mat **stagesensip) 722 { 723 TS_RK *rk = (TS_RK*)ts->data; 724 RKTableau tab = rk->tableau; 725 726 PetscFunctionBegin; 727 if (ns) *ns = tab->s; 728 if (stagesensip) *stagesensip = rk->MatsFwdStageSensip; 729 PetscFunctionReturn(0); 730 } 731 732 static PetscErrorCode TSForwardSetUp_RK(TS ts) 733 { 734 TS_RK *rk = (TS_RK*)ts->data; 735 RKTableau tab = rk->tableau; 736 PetscInt i; 737 738 PetscFunctionBegin; 739 /* backup sensitivity results for roll-backs */ 740 PetscCall(MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatFwdSensip0)); 741 742 PetscCall(PetscMalloc1(tab->s,&rk->MatsFwdStageSensip)); 743 PetscCall(PetscMalloc1(tab->s,&rk->MatsFwdSensipTemp)); 744 for (i=0; i<tab->s; i++) { 745 PetscCall(MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdStageSensip[i])); 746 PetscCall(MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&rk->MatsFwdSensipTemp[i])); 747 } 748 PetscCall(VecDuplicate(ts->vec_sol,&rk->VecDeltaFwdSensipCol)); 749 PetscFunctionReturn(0); 750 } 751 752 static PetscErrorCode TSForwardReset_RK(TS ts) 753 { 754 TS_RK *rk = (TS_RK*)ts->data; 755 RKTableau tab = rk->tableau; 756 PetscInt i; 757 758 PetscFunctionBegin; 759 PetscCall(MatDestroy(&rk->MatFwdSensip0)); 760 if (rk->MatsFwdStageSensip) { 761 for (i=0; i<tab->s; i++) { 762 PetscCall(MatDestroy(&rk->MatsFwdStageSensip[i])); 763 } 764 PetscCall(PetscFree(rk->MatsFwdStageSensip)); 765 } 766 if (rk->MatsFwdSensipTemp) { 767 for (i=0; i<tab->s; i++) { 768 PetscCall(MatDestroy(&rk->MatsFwdSensipTemp[i])); 769 } 770 PetscCall(PetscFree(rk->MatsFwdSensipTemp)); 771 } 772 PetscCall(VecDestroy(&rk->VecDeltaFwdSensipCol)); 773 PetscFunctionReturn(0); 774 } 775 776 static PetscErrorCode TSStep_RK(TS ts) 777 { 778 TS_RK *rk = (TS_RK*)ts->data; 779 RKTableau tab = rk->tableau; 780 const PetscInt s = tab->s; 781 const PetscReal *A = tab->A,*c = tab->c; 782 PetscScalar *w = rk->work; 783 Vec *Y = rk->Y,*YdotRHS = rk->YdotRHS; 784 PetscBool FSAL = tab->FSAL; 785 TSAdapt adapt; 786 PetscInt i,j; 787 PetscInt rejections = 0; 788 PetscBool stageok,accept = PETSC_TRUE; 789 PetscReal next_time_step = ts->time_step; 790 791 PetscFunctionBegin; 792 if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE; 793 if (FSAL) PetscCall(VecCopy(YdotRHS[s-1],YdotRHS[0])); 794 795 rk->status = TS_STEP_INCOMPLETE; 796 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 797 PetscReal t = ts->ptime; 798 PetscReal h = ts->time_step; 799 for (i=0; i<s; i++) { 800 rk->stage_time = t + h*c[i]; 801 PetscCall(TSPreStage(ts,rk->stage_time)); 802 PetscCall(VecCopy(ts->vec_sol,Y[i])); 803 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 804 PetscCall(VecMAXPY(Y[i],i,w,YdotRHS)); 805 PetscCall(TSPostStage(ts,rk->stage_time,i,Y)); 806 PetscCall(TSGetAdapt(ts,&adapt)); 807 PetscCall(TSAdaptCheckStage(adapt,ts,rk->stage_time,Y[i],&stageok)); 808 if (!stageok) goto reject_step; 809 if (FSAL && !i) continue; 810 PetscCall(TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i])); 811 } 812 813 rk->status = TS_STEP_INCOMPLETE; 814 PetscCall(TSEvaluateStep(ts,tab->order,ts->vec_sol,NULL)); 815 rk->status = TS_STEP_PENDING; 816 PetscCall(TSGetAdapt(ts,&adapt)); 817 PetscCall(TSAdaptCandidatesClear(adapt)); 818 PetscCall(TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE)); 819 PetscCall(TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept)); 820 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 821 if (!accept) { /* Roll back the current step */ 822 PetscCall(TSRollBack_RK(ts)); 823 ts->time_step = next_time_step; 824 goto reject_step; 825 } 826 827 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */ 828 rk->ptime = ts->ptime; 829 rk->time_step = ts->time_step; 830 } 831 832 ts->ptime += ts->time_step; 833 ts->time_step = next_time_step; 834 break; 835 836 reject_step: 837 ts->reject++; accept = PETSC_FALSE; 838 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 839 ts->reason = TS_DIVERGED_STEP_REJECTED; 840 PetscCall(PetscInfo(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections)); 841 } 842 } 843 PetscFunctionReturn(0); 844 } 845 846 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 847 { 848 TS_RK *rk = (TS_RK*)ts->data; 849 RKTableau tab = rk->tableau; 850 PetscInt s = tab->s; 851 852 PetscFunctionBegin; 853 if (ts->adjointsetupcalled++) PetscFunctionReturn(0); 854 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam)); 855 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&rk->VecsSensiTemp)); 856 if (ts->vecs_sensip) { 857 PetscCall(VecDuplicate(ts->vecs_sensip[0],&rk->VecDeltaMu)); 858 } 859 if (ts->vecs_sensi2) { 860 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0],s*ts->numcost,&rk->VecsDeltaLam2)); 861 PetscCall(VecDuplicateVecs(ts->vecs_sensi2[0],ts->numcost,&rk->VecsSensi2Temp)); 862 } 863 if (ts->vecs_sensi2p) { 864 PetscCall(VecDuplicate(ts->vecs_sensi2p[0],&rk->VecDeltaMu2)); 865 } 866 PetscFunctionReturn(0); 867 } 868 869 /* 870 Assumptions: 871 - TSStep_RK() always evaluates the step with b, not bembed. 872 */ 873 static PetscErrorCode TSAdjointStep_RK(TS ts) 874 { 875 TS_RK *rk = (TS_RK*)ts->data; 876 TS quadts = ts->quadraturets; 877 RKTableau tab = rk->tableau; 878 Mat J,Jpre,Jquad; 879 const PetscInt s = tab->s; 880 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 881 PetscScalar *w = rk->work,*xarr; 882 Vec *Y = rk->Y,*VecsDeltaLam = rk->VecsDeltaLam,VecDeltaMu = rk->VecDeltaMu,*VecsSensiTemp = rk->VecsSensiTemp; 883 Vec *VecsDeltaLam2 = rk->VecsDeltaLam2,VecDeltaMu2 = rk->VecDeltaMu2,*VecsSensi2Temp = rk->VecsSensi2Temp; 884 Vec VecDRDUTransCol = ts->vec_drdu_col,VecDRDPTransCol = ts->vec_drdp_col; 885 PetscInt i,j,nadj; 886 PetscReal t = ts->ptime; 887 PetscReal h = ts->time_step; 888 889 PetscFunctionBegin; 890 rk->status = TS_STEP_INCOMPLETE; 891 892 PetscCall(TSGetRHSJacobian(ts,&J,&Jpre,NULL,NULL)); 893 if (quadts) { 894 PetscCall(TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL)); 895 } 896 for (i=s-1; i>=0; i--) { 897 if (tab->FSAL && i == s-1) { 898 /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/ 899 continue; 900 } 901 rk->stage_time = t + h*(1.0-c[i]); 902 PetscCall(TSComputeSNESJacobian(ts,Y[i],J,Jpre)); 903 if (quadts) { 904 PetscCall(TSComputeRHSJacobian(quadts,rk->stage_time,Y[i],Jquad,Jquad)); /* get r_u^T */ 905 } 906 if (ts->vecs_sensip) { 907 PetscCall(TSComputeRHSJacobianP(ts,rk->stage_time,Y[i],ts->Jacprhs)); /* get f_p */ 908 if (quadts) { 909 PetscCall(TSComputeRHSJacobianP(quadts,rk->stage_time,Y[i],quadts->Jacprhs)); /* get f_p for the quadrature */ 910 } 911 } 912 913 if (b[i]) { 914 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]/b[i]; /* coefficients for computing VecsSensiTemp */ 915 } else { 916 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]; /* coefficients for computing VecsSensiTemp */ 917 } 918 919 for (nadj=0; nadj<ts->numcost; nadj++) { 920 /* Stage values of lambda */ 921 if (b[i]) { 922 /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */ 923 PetscCall(VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj])); /* VecDeltaLam is an vec array of size s by numcost */ 924 PetscCall(VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1])); 925 PetscCall(MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i])); /* VecsSensiTemp will be reused by 2nd-order adjoint */ 926 PetscCall(VecScale(VecsDeltaLam[nadj*s+i],-h*b[i])); 927 if (quadts) { 928 PetscCall(MatDenseGetColumn(Jquad,nadj,&xarr)); 929 PetscCall(VecPlaceArray(VecDRDUTransCol,xarr)); 930 PetscCall(VecAXPY(VecsDeltaLam[nadj*s+i],-h*b[i],VecDRDUTransCol)); 931 PetscCall(VecResetArray(VecDRDUTransCol)); 932 PetscCall(MatDenseRestoreColumn(Jquad,&xarr)); 933 } 934 } else { 935 /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */ 936 PetscCall(VecSet(VecsSensiTemp[nadj],0)); 937 PetscCall(VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1])); 938 PetscCall(MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i])); 939 PetscCall(VecScale(VecsDeltaLam[nadj*s+i],-h)); 940 } 941 942 /* Stage values of mu */ 943 if (ts->vecs_sensip) { 944 if (b[i]) { 945 PetscCall(MatMultTranspose(ts->Jacprhs,VecsSensiTemp[nadj],VecDeltaMu)); 946 PetscCall(VecScale(VecDeltaMu,-h*b[i])); 947 if (quadts) { 948 PetscCall(MatDenseGetColumn(quadts->Jacprhs,nadj,&xarr)); 949 PetscCall(VecPlaceArray(VecDRDPTransCol,xarr)); 950 PetscCall(VecAXPY(VecDeltaMu,-h*b[i],VecDRDPTransCol)); 951 PetscCall(VecResetArray(VecDRDPTransCol)); 952 PetscCall(MatDenseRestoreColumn(quadts->Jacprhs,&xarr)); 953 } 954 } else { 955 PetscCall(VecScale(VecDeltaMu,-h)); 956 } 957 PetscCall(VecAXPY(ts->vecs_sensip[nadj],1.,VecDeltaMu)); /* update sensip for each stage */ 958 } 959 } 960 961 if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */ 962 /* Get w1 at t_{n+1} from TLM matrix */ 963 PetscCall(MatDenseGetColumn(rk->MatsFwdStageSensip[i],0,&xarr)); 964 PetscCall(VecPlaceArray(ts->vec_sensip_col,xarr)); 965 /* lambda_s^T F_UU w_1 */ 966 PetscCall(TSComputeRHSHessianProductFunctionUU(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_guu)); 967 if (quadts) { 968 /* R_UU w_1 */ 969 PetscCall(TSComputeRHSHessianProductFunctionUU(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_guu)); 970 } 971 if (ts->vecs_sensip) { 972 /* lambda_s^T F_UP w_2 */ 973 PetscCall(TSComputeRHSHessianProductFunctionUP(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gup)); 974 if (quadts) { 975 /* R_UP w_2 */ 976 PetscCall(TSComputeRHSHessianProductFunctionUP(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gup)); 977 } 978 } 979 if (ts->vecs_sensi2p) { 980 /* lambda_s^T F_PU w_1 */ 981 PetscCall(TSComputeRHSHessianProductFunctionPU(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_gpu)); 982 /* lambda_s^T F_PP w_2 */ 983 PetscCall(TSComputeRHSHessianProductFunctionPP(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gpp)); 984 if (b[i] && quadts) { 985 /* R_PU w_1 */ 986 PetscCall(TSComputeRHSHessianProductFunctionPU(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gpu)); 987 /* R_PP w_2 */ 988 PetscCall(TSComputeRHSHessianProductFunctionPP(quadts,rk->stage_time,Y[i],NULL,ts->vec_dir,ts->vecs_gpp)); 989 } 990 } 991 PetscCall(VecResetArray(ts->vec_sensip_col)); 992 PetscCall(MatDenseRestoreColumn(rk->MatsFwdStageSensip[i],&xarr)); 993 994 for (nadj=0; nadj<ts->numcost; nadj++) { 995 /* Stage values of lambda */ 996 if (b[i]) { 997 /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */ 998 PetscCall(VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj])); 999 PetscCall(VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1])); 1000 PetscCall(MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i])); 1001 PetscCall(VecScale(VecsDeltaLam2[nadj*s+i],-h*b[i])); 1002 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_guu[nadj])); 1003 if (ts->vecs_sensip) { 1004 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_gup[nadj])); 1005 } 1006 } else { 1007 /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */ 1008 PetscCall(VecSet(VecsDeltaLam2[nadj*s+i],0)); 1009 PetscCall(VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1])); 1010 PetscCall(MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i])); 1011 PetscCall(VecScale(VecsDeltaLam2[nadj*s+i],-h)); 1012 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_guu[nadj])); 1013 if (ts->vecs_sensip) { 1014 PetscCall(VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_gup[nadj])); 1015 } 1016 } 1017 if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */ 1018 PetscCall(MatMultTranspose(ts->Jacprhs,VecsSensi2Temp[nadj],VecDeltaMu2)); 1019 if (b[i]) { 1020 PetscCall(VecScale(VecDeltaMu2,-h*b[i])); 1021 PetscCall(VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpu[nadj])); 1022 PetscCall(VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpp[nadj])); 1023 } else { 1024 PetscCall(VecScale(VecDeltaMu2,-h)); 1025 PetscCall(VecAXPY(VecDeltaMu2,-h,ts->vecs_gpu[nadj])); 1026 PetscCall(VecAXPY(VecDeltaMu2,-h,ts->vecs_gpp[nadj])); 1027 } 1028 PetscCall(VecAXPY(ts->vecs_sensi2p[nadj],1,VecDeltaMu2)); /* update sensi2p for each stage */ 1029 } 1030 } 1031 } 1032 } 1033 1034 for (j=0; j<s; j++) w[j] = 1.0; 1035 for (nadj=0; nadj<ts->numcost; nadj++) { /* no need to do this for mu's */ 1036 PetscCall(VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecsDeltaLam[nadj*s])); 1037 if (ts->vecs_sensi2) { 1038 PetscCall(VecMAXPY(ts->vecs_sensi2[nadj],s,w,&VecsDeltaLam2[nadj*s])); 1039 } 1040 } 1041 rk->status = TS_STEP_COMPLETE; 1042 PetscFunctionReturn(0); 1043 } 1044 1045 static PetscErrorCode TSAdjointReset_RK(TS ts) 1046 { 1047 TS_RK *rk = (TS_RK*)ts->data; 1048 RKTableau tab = rk->tableau; 1049 1050 PetscFunctionBegin; 1051 PetscCall(VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam)); 1052 PetscCall(VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp)); 1053 PetscCall(VecDestroy(&rk->VecDeltaMu)); 1054 PetscCall(VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam2)); 1055 PetscCall(VecDestroy(&rk->VecDeltaMu2)); 1056 PetscCall(VecDestroyVecs(ts->numcost,&rk->VecsSensi2Temp)); 1057 PetscFunctionReturn(0); 1058 } 1059 1060 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 1061 { 1062 TS_RK *rk = (TS_RK*)ts->data; 1063 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 1064 PetscReal h; 1065 PetscReal tt,t; 1066 PetscScalar *b; 1067 const PetscReal *B = rk->tableau->binterp; 1068 1069 PetscFunctionBegin; 1070 PetscCheck(B,PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 1071 1072 switch (rk->status) { 1073 case TS_STEP_INCOMPLETE: 1074 case TS_STEP_PENDING: 1075 h = ts->time_step; 1076 t = (itime - ts->ptime)/h; 1077 break; 1078 case TS_STEP_COMPLETE: 1079 h = ts->ptime - ts->ptime_prev; 1080 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1081 break; 1082 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1083 } 1084 PetscCall(PetscMalloc1(s,&b)); 1085 for (i=0; i<s; i++) b[i] = 0; 1086 for (j=0,tt=t; j<p; j++,tt*=t) { 1087 for (i=0; i<s; i++) { 1088 b[i] += h * B[i*p+j] * tt; 1089 } 1090 } 1091 PetscCall(VecCopy(rk->Y[0],X)); 1092 PetscCall(VecMAXPY(X,s,b,rk->YdotRHS)); 1093 PetscCall(PetscFree(b)); 1094 PetscFunctionReturn(0); 1095 } 1096 1097 /*------------------------------------------------------------*/ 1098 1099 static PetscErrorCode TSRKTableauReset(TS ts) 1100 { 1101 TS_RK *rk = (TS_RK*)ts->data; 1102 RKTableau tab = rk->tableau; 1103 1104 PetscFunctionBegin; 1105 if (!tab) PetscFunctionReturn(0); 1106 PetscCall(PetscFree(rk->work)); 1107 PetscCall(VecDestroyVecs(tab->s,&rk->Y)); 1108 PetscCall(VecDestroyVecs(tab->s,&rk->YdotRHS)); 1109 PetscFunctionReturn(0); 1110 } 1111 1112 static PetscErrorCode TSReset_RK(TS ts) 1113 { 1114 PetscFunctionBegin; 1115 PetscCall(TSRKTableauReset(ts)); 1116 if (ts->use_splitrhsfunction) { 1117 PetscCall(PetscTryMethod(ts,"TSReset_RK_MultirateSplit_C",(TS),(ts))); 1118 } else { 1119 PetscCall(PetscTryMethod(ts,"TSReset_RK_MultirateNonsplit_C",(TS),(ts))); 1120 } 1121 PetscFunctionReturn(0); 1122 } 1123 1124 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 1125 { 1126 PetscFunctionBegin; 1127 PetscFunctionReturn(0); 1128 } 1129 1130 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1131 { 1132 PetscFunctionBegin; 1133 PetscFunctionReturn(0); 1134 } 1135 1136 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 1137 { 1138 PetscFunctionBegin; 1139 PetscFunctionReturn(0); 1140 } 1141 1142 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1143 { 1144 PetscFunctionBegin; 1145 PetscFunctionReturn(0); 1146 } 1147 1148 static PetscErrorCode TSRKTableauSetUp(TS ts) 1149 { 1150 TS_RK *rk = (TS_RK*)ts->data; 1151 RKTableau tab = rk->tableau; 1152 1153 PetscFunctionBegin; 1154 PetscCall(PetscMalloc1(tab->s,&rk->work)); 1155 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y)); 1156 PetscCall(VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS)); 1157 PetscFunctionReturn(0); 1158 } 1159 1160 static PetscErrorCode TSSetUp_RK(TS ts) 1161 { 1162 TS quadts = ts->quadraturets; 1163 DM dm; 1164 1165 PetscFunctionBegin; 1166 PetscCall(TSCheckImplicitTerm(ts)); 1167 PetscCall(TSRKTableauSetUp(ts)); 1168 if (quadts && ts->costintegralfwd) { 1169 Mat Jquad; 1170 PetscCall(TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL)); 1171 } 1172 PetscCall(TSGetDM(ts,&dm)); 1173 PetscCall(DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts)); 1174 PetscCall(DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts)); 1175 if (ts->use_splitrhsfunction) { 1176 PetscCall(PetscTryMethod(ts,"TSSetUp_RK_MultirateSplit_C",(TS),(ts))); 1177 } else { 1178 PetscCall(PetscTryMethod(ts,"TSSetUp_RK_MultirateNonsplit_C",(TS),(ts))); 1179 } 1180 PetscFunctionReturn(0); 1181 } 1182 1183 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 1184 { 1185 TS_RK *rk = (TS_RK*)ts->data; 1186 PetscErrorCode ierr; 1187 1188 PetscFunctionBegin; 1189 PetscCall(PetscOptionsHead(PetscOptionsObject,"RK ODE solver options")); 1190 { 1191 RKTableauLink link; 1192 PetscInt count,choice; 1193 PetscBool flg,use_multirate = PETSC_FALSE; 1194 const char **namelist; 1195 1196 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 1197 PetscCall(PetscMalloc1(count,(char***)&namelist)); 1198 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1199 PetscCall(PetscOptionsBool("-ts_rk_multirate","Use interpolation-based multirate RK method","TSRKSetMultirate",rk->use_multirate,&use_multirate,&flg)); 1200 if (flg) { 1201 PetscCall(TSRKSetMultirate(ts,use_multirate)); 1202 } 1203 PetscCall(PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg)); 1204 if (flg) PetscCall(TSRKSetType(ts,namelist[choice])); 1205 PetscCall(PetscFree(namelist)); 1206 } 1207 PetscCall(PetscOptionsTail()); 1208 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)ts),NULL,"Multirate methods options","");PetscCall(ierr); 1209 PetscCall(PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL)); 1210 ierr = PetscOptionsEnd();PetscCall(ierr); 1211 PetscFunctionReturn(0); 1212 } 1213 1214 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 1215 { 1216 TS_RK *rk = (TS_RK*)ts->data; 1217 PetscBool iascii; 1218 1219 PetscFunctionBegin; 1220 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1221 if (iascii) { 1222 RKTableau tab = rk->tableau; 1223 TSRKType rktype; 1224 const PetscReal *c; 1225 PetscInt s; 1226 char buf[512]; 1227 PetscBool FSAL; 1228 1229 PetscCall(TSRKGetType(ts,&rktype)); 1230 PetscCall(TSRKGetTableau(ts,&s,NULL,NULL,&c,NULL,NULL,NULL,&FSAL)); 1231 PetscCall(PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype)); 1232 PetscCall(PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order)); 1233 PetscCall(PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",FSAL ? "yes" : "no")); 1234 PetscCall(PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",s,c)); 1235 PetscCall(PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf)); 1236 } 1237 PetscFunctionReturn(0); 1238 } 1239 1240 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 1241 { 1242 TSAdapt adapt; 1243 1244 PetscFunctionBegin; 1245 PetscCall(TSGetAdapt(ts,&adapt)); 1246 PetscCall(TSAdaptLoad(adapt,viewer)); 1247 PetscFunctionReturn(0); 1248 } 1249 1250 /*@ 1251 TSRKGetOrder - Get the order of RK scheme 1252 1253 Not collective 1254 1255 Input Parameter: 1256 . ts - timestepping context 1257 1258 Output Parameter: 1259 . order - order of RK-scheme 1260 1261 Level: intermediate 1262 1263 .seealso: TSRKGetType() 1264 @*/ 1265 PetscErrorCode TSRKGetOrder(TS ts,PetscInt *order) 1266 { 1267 PetscFunctionBegin; 1268 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1269 PetscValidIntPointer(order,2); 1270 PetscCall(PetscUseMethod(ts,"TSRKGetOrder_C",(TS,PetscInt*),(ts,order))); 1271 PetscFunctionReturn(0); 1272 } 1273 1274 /*@C 1275 TSRKSetType - Set the type of RK scheme 1276 1277 Logically collective 1278 1279 Input Parameters: 1280 + ts - timestepping context 1281 - rktype - type of RK-scheme 1282 1283 Options Database: 1284 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1285 1286 Level: intermediate 1287 1288 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK2B, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS, TSRK6VR, TSRK7VR, TSRK8VR 1289 @*/ 1290 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1291 { 1292 PetscFunctionBegin; 1293 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1294 PetscValidCharPointer(rktype,2); 1295 PetscCall(PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype))); 1296 PetscFunctionReturn(0); 1297 } 1298 1299 /*@C 1300 TSRKGetType - Get the type of RK scheme 1301 1302 Not collective 1303 1304 Input Parameter: 1305 . ts - timestepping context 1306 1307 Output Parameter: 1308 . rktype - type of RK-scheme 1309 1310 Level: intermediate 1311 1312 .seealso: TSRKSetType() 1313 @*/ 1314 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1315 { 1316 PetscFunctionBegin; 1317 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1318 PetscCall(PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype))); 1319 PetscFunctionReturn(0); 1320 } 1321 1322 static PetscErrorCode TSRKGetOrder_RK(TS ts,PetscInt *order) 1323 { 1324 TS_RK *rk = (TS_RK*)ts->data; 1325 1326 PetscFunctionBegin; 1327 *order = rk->tableau->order; 1328 PetscFunctionReturn(0); 1329 } 1330 1331 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1332 { 1333 TS_RK *rk = (TS_RK*)ts->data; 1334 1335 PetscFunctionBegin; 1336 *rktype = rk->tableau->name; 1337 PetscFunctionReturn(0); 1338 } 1339 1340 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1341 { 1342 TS_RK *rk = (TS_RK*)ts->data; 1343 PetscBool match; 1344 RKTableauLink link; 1345 1346 PetscFunctionBegin; 1347 if (rk->tableau) { 1348 PetscCall(PetscStrcmp(rk->tableau->name,rktype,&match)); 1349 if (match) PetscFunctionReturn(0); 1350 } 1351 for (link = RKTableauList; link; link=link->next) { 1352 PetscCall(PetscStrcmp(link->tab.name,rktype,&match)); 1353 if (match) { 1354 if (ts->setupcalled) PetscCall(TSRKTableauReset(ts)); 1355 rk->tableau = &link->tab; 1356 if (ts->setupcalled) PetscCall(TSRKTableauSetUp(ts)); 1357 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1358 PetscFunctionReturn(0); 1359 } 1360 } 1361 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1362 } 1363 1364 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1365 { 1366 TS_RK *rk = (TS_RK*)ts->data; 1367 1368 PetscFunctionBegin; 1369 if (ns) *ns = rk->tableau->s; 1370 if (Y) *Y = rk->Y; 1371 PetscFunctionReturn(0); 1372 } 1373 1374 static PetscErrorCode TSDestroy_RK(TS ts) 1375 { 1376 PetscFunctionBegin; 1377 PetscCall(TSReset_RK(ts)); 1378 if (ts->dm) { 1379 PetscCall(DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts)); 1380 PetscCall(DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts)); 1381 } 1382 PetscCall(PetscFree(ts->data)); 1383 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetOrder_C",NULL)); 1384 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL)); 1385 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL)); 1386 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetTableau_C",NULL)); 1387 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",NULL)); 1388 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",NULL)); 1389 PetscFunctionReturn(0); 1390 } 1391 1392 /* 1393 This defines the nonlinear equation that is to be solved with SNES 1394 We do not need to solve the equation; we just use SNES to approximate the Jacobian 1395 */ 1396 static PetscErrorCode SNESTSFormFunction_RK(SNES snes,Vec x,Vec y,TS ts) 1397 { 1398 TS_RK *rk = (TS_RK*)ts->data; 1399 DM dm,dmsave; 1400 1401 PetscFunctionBegin; 1402 PetscCall(SNESGetDM(snes,&dm)); 1403 /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */ 1404 dmsave = ts->dm; 1405 ts->dm = dm; 1406 PetscCall(TSComputeRHSFunction(ts,rk->stage_time,x,y)); 1407 ts->dm = dmsave; 1408 PetscFunctionReturn(0); 1409 } 1410 1411 static PetscErrorCode SNESTSFormJacobian_RK(SNES snes,Vec x,Mat A,Mat B,TS ts) 1412 { 1413 TS_RK *rk = (TS_RK*)ts->data; 1414 DM dm,dmsave; 1415 1416 PetscFunctionBegin; 1417 PetscCall(SNESGetDM(snes,&dm)); 1418 dmsave = ts->dm; 1419 ts->dm = dm; 1420 PetscCall(TSComputeRHSJacobian(ts,rk->stage_time,x,A,B)); 1421 ts->dm = dmsave; 1422 PetscFunctionReturn(0); 1423 } 1424 1425 /*@C 1426 TSRKSetMultirate - Use the interpolation-based multirate RK method 1427 1428 Logically collective 1429 1430 Input Parameters: 1431 + ts - timestepping context 1432 - 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 1433 1434 Options Database: 1435 . -ts_rk_multirate - <true,false> 1436 1437 Notes: 1438 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). 1439 1440 Level: intermediate 1441 1442 .seealso: TSRKGetMultirate() 1443 @*/ 1444 PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate) 1445 { 1446 PetscFunctionBegin; 1447 PetscCall(PetscTryMethod(ts,"TSRKSetMultirate_C",(TS,PetscBool),(ts,use_multirate))); 1448 PetscFunctionReturn(0); 1449 } 1450 1451 /*@C 1452 TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method 1453 1454 Not collective 1455 1456 Input Parameter: 1457 . ts - timestepping context 1458 1459 Output Parameter: 1460 . use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise 1461 1462 Level: intermediate 1463 1464 .seealso: TSRKSetMultirate() 1465 @*/ 1466 PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate) 1467 { 1468 PetscFunctionBegin; 1469 PetscCall(PetscUseMethod(ts,"TSRKGetMultirate_C",(TS,PetscBool*),(ts,use_multirate))); 1470 PetscFunctionReturn(0); 1471 } 1472 1473 /*MC 1474 TSRK - ODE and DAE solver using Runge-Kutta schemes 1475 1476 The user should provide the right hand side of the equation 1477 using TSSetRHSFunction(). 1478 1479 Notes: 1480 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1481 1482 Level: beginner 1483 1484 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRK2D, TTSRK2E, TSRK3, 1485 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirate(), TSRKGetMultirate() 1486 1487 M*/ 1488 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1489 { 1490 TS_RK *rk; 1491 1492 PetscFunctionBegin; 1493 PetscCall(TSRKInitializePackage()); 1494 1495 ts->ops->reset = TSReset_RK; 1496 ts->ops->destroy = TSDestroy_RK; 1497 ts->ops->view = TSView_RK; 1498 ts->ops->load = TSLoad_RK; 1499 ts->ops->setup = TSSetUp_RK; 1500 ts->ops->interpolate = TSInterpolate_RK; 1501 ts->ops->step = TSStep_RK; 1502 ts->ops->evaluatestep = TSEvaluateStep_RK; 1503 ts->ops->rollback = TSRollBack_RK; 1504 ts->ops->setfromoptions = TSSetFromOptions_RK; 1505 ts->ops->getstages = TSGetStages_RK; 1506 1507 ts->ops->snesfunction = SNESTSFormFunction_RK; 1508 ts->ops->snesjacobian = SNESTSFormJacobian_RK; 1509 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1510 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1511 ts->ops->adjointstep = TSAdjointStep_RK; 1512 ts->ops->adjointreset = TSAdjointReset_RK; 1513 1514 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1515 ts->ops->forwardsetup = TSForwardSetUp_RK; 1516 ts->ops->forwardreset = TSForwardReset_RK; 1517 ts->ops->forwardstep = TSForwardStep_RK; 1518 ts->ops->forwardgetstages= TSForwardGetStages_RK; 1519 1520 PetscCall(PetscNewLog(ts,&rk)); 1521 ts->data = (void*)rk; 1522 1523 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetOrder_C",TSRKGetOrder_RK)); 1524 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK)); 1525 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK)); 1526 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetTableau_C",TSRKGetTableau_RK)); 1527 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",TSRKSetMultirate_RK)); 1528 PetscCall(PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",TSRKGetMultirate_RK)); 1529 1530 PetscCall(TSRKSetType(ts,TSRKDefault)); 1531 rk->dtratio = 1; 1532 PetscFunctionReturn(0); 1533 } 1534