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