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 PetscAssertPointer(name, 1); 452 PetscAssertPointer(A, 4); 453 if (b) PetscAssertPointer(b, 5); 454 if (c) PetscAssertPointer(c, 6); 455 if (bembed) PetscAssertPointer(bembed, 7); 456 if (binterp || p > 1) PetscAssertPointer(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 = (PetscBool)(tab->FSAL && !rk->newtableau); 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 rk->newtableau = PETSC_FALSE; 806 807 rk->status = TS_STEP_INCOMPLETE; 808 while (!ts->reason && rk->status != TS_STEP_COMPLETE) { 809 PetscReal t = ts->ptime; 810 PetscReal h = ts->time_step; 811 for (i = 0; i < s; i++) { 812 rk->stage_time = t + h * c[i]; 813 PetscCall(TSPreStage(ts, rk->stage_time)); 814 PetscCall(VecCopy(ts->vec_sol, Y[i])); 815 for (j = 0; j < i; j++) w[j] = h * A[i * s + j]; 816 PetscCall(VecMAXPY(Y[i], i, w, YdotRHS)); 817 PetscCall(TSPostStage(ts, rk->stage_time, i, Y)); 818 PetscCall(TSGetAdapt(ts, &adapt)); 819 PetscCall(TSAdaptCheckStage(adapt, ts, rk->stage_time, Y[i], &stageok)); 820 if (!stageok) goto reject_step; 821 if (FSAL && !i) continue; 822 PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i])); 823 } 824 825 rk->status = TS_STEP_INCOMPLETE; 826 PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL)); 827 rk->status = TS_STEP_PENDING; 828 PetscCall(TSGetAdapt(ts, &adapt)); 829 PetscCall(TSAdaptCandidatesClear(adapt)); 830 PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE)); 831 PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept)); 832 rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 833 if (!accept) { /* Roll back the current step */ 834 PetscCall(TSRollBack_RK(ts)); 835 ts->time_step = next_time_step; 836 goto reject_step; 837 } 838 839 if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */ 840 rk->ptime = ts->ptime; 841 rk->time_step = ts->time_step; 842 } 843 844 ts->ptime += ts->time_step; 845 ts->time_step = next_time_step; 846 break; 847 848 reject_step: 849 ts->reject++; 850 accept = PETSC_FALSE; 851 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 852 ts->reason = TS_DIVERGED_STEP_REJECTED; 853 PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections)); 854 } 855 } 856 PetscFunctionReturn(PETSC_SUCCESS); 857 } 858 859 static PetscErrorCode TSAdjointSetUp_RK(TS ts) 860 { 861 TS_RK *rk = (TS_RK *)ts->data; 862 RKTableau tab = rk->tableau; 863 PetscInt s = tab->s; 864 865 PetscFunctionBegin; 866 if (ts->adjointsetupcalled++) PetscFunctionReturn(PETSC_SUCCESS); 867 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], s * ts->numcost, &rk->VecsDeltaLam)); 868 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &rk->VecsSensiTemp)); 869 if (ts->vecs_sensip) PetscCall(VecDuplicate(ts->vecs_sensip[0], &rk->VecDeltaMu)); 870 if (ts->vecs_sensi2) { 871 PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], s * ts->numcost, &rk->VecsDeltaLam2)); 872 PetscCall(VecDuplicateVecs(ts->vecs_sensi2[0], ts->numcost, &rk->VecsSensi2Temp)); 873 } 874 if (ts->vecs_sensi2p) PetscCall(VecDuplicate(ts->vecs_sensi2p[0], &rk->VecDeltaMu2)); 875 PetscFunctionReturn(PETSC_SUCCESS); 876 } 877 878 /* 879 Assumptions: 880 - TSStep_RK() always evaluates the step with b, not bembed. 881 */ 882 static PetscErrorCode TSAdjointStep_RK(TS ts) 883 { 884 TS_RK *rk = (TS_RK *)ts->data; 885 TS quadts = ts->quadraturets; 886 RKTableau tab = rk->tableau; 887 Mat J, Jpre, Jquad; 888 const PetscInt s = tab->s; 889 const PetscReal *A = tab->A, *b = tab->b, *c = tab->c; 890 PetscScalar *w = rk->work, *xarr; 891 Vec *Y = rk->Y, *VecsDeltaLam = rk->VecsDeltaLam, VecDeltaMu = rk->VecDeltaMu, *VecsSensiTemp = rk->VecsSensiTemp; 892 Vec *VecsDeltaLam2 = rk->VecsDeltaLam2, VecDeltaMu2 = rk->VecDeltaMu2, *VecsSensi2Temp = rk->VecsSensi2Temp; 893 Vec VecDRDUTransCol = ts->vec_drdu_col, VecDRDPTransCol = ts->vec_drdp_col; 894 PetscInt i, j, nadj; 895 PetscReal t = ts->ptime; 896 PetscReal h = ts->time_step; 897 898 PetscFunctionBegin; 899 rk->status = TS_STEP_INCOMPLETE; 900 901 PetscCall(TSGetRHSJacobian(ts, &J, &Jpre, NULL, NULL)); 902 if (quadts) PetscCall(TSGetRHSJacobian(quadts, &Jquad, NULL, NULL, NULL)); 903 for (i = s - 1; i >= 0; i--) { 904 if (tab->FSAL && i == s - 1) { 905 /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/ 906 continue; 907 } 908 rk->stage_time = t + h * (1.0 - c[i]); 909 PetscCall(TSComputeSNESJacobian(ts, Y[i], J, Jpre)); 910 if (quadts) { PetscCall(TSComputeRHSJacobian(quadts, rk->stage_time, Y[i], Jquad, Jquad)); /* get r_u^T */ } 911 if (ts->vecs_sensip) { 912 PetscCall(TSComputeRHSJacobianP(ts, rk->stage_time, Y[i], ts->Jacprhs)); /* get f_p */ 913 if (quadts) { PetscCall(TSComputeRHSJacobianP(quadts, rk->stage_time, Y[i], quadts->Jacprhs)); /* get f_p for the quadrature */ } 914 } 915 916 if (b[i]) { 917 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i] / b[i]; /* coefficients for computing VecsSensiTemp */ 918 } else { 919 for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; /* coefficients for computing VecsSensiTemp */ 920 } 921 922 for (nadj = 0; nadj < ts->numcost; nadj++) { 923 /* Stage values of lambda */ 924 if (b[i]) { 925 /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */ 926 PetscCall(VecCopy(ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); /* VecDeltaLam is an vec array of size s by numcost */ 927 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 928 PetscCall(MatMultTranspose(J, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); /* VecsSensiTemp will be reused by 2nd-order adjoint */ 929 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], -h * b[i])); 930 if (quadts) { 931 PetscCall(MatDenseGetColumn(Jquad, nadj, &xarr)); 932 PetscCall(VecPlaceArray(VecDRDUTransCol, xarr)); 933 PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -h * b[i], VecDRDUTransCol)); 934 PetscCall(VecResetArray(VecDRDUTransCol)); 935 PetscCall(MatDenseRestoreColumn(Jquad, &xarr)); 936 } 937 } else { 938 /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */ 939 PetscCall(VecSet(VecsSensiTemp[nadj], 0)); 940 PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1])); 941 PetscCall(MatMultTranspose(J, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); 942 PetscCall(VecScale(VecsDeltaLam[nadj * s + i], -h)); 943 } 944 945 /* Stage values of mu */ 946 if (ts->vecs_sensip) { 947 if (b[i]) { 948 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecDeltaMu)); 949 PetscCall(VecScale(VecDeltaMu, -h * b[i])); 950 if (quadts) { 951 PetscCall(MatDenseGetColumn(quadts->Jacprhs, nadj, &xarr)); 952 PetscCall(VecPlaceArray(VecDRDPTransCol, xarr)); 953 PetscCall(VecAXPY(VecDeltaMu, -h * b[i], VecDRDPTransCol)); 954 PetscCall(VecResetArray(VecDRDPTransCol)); 955 PetscCall(MatDenseRestoreColumn(quadts->Jacprhs, &xarr)); 956 } 957 } else { 958 PetscCall(VecScale(VecDeltaMu, -h)); 959 } 960 PetscCall(VecAXPY(ts->vecs_sensip[nadj], 1., VecDeltaMu)); /* update sensip for each stage */ 961 } 962 } 963 964 if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */ 965 /* Get w1 at t_{n+1} from TLM matrix */ 966 PetscCall(MatDenseGetColumn(rk->MatsFwdStageSensip[i], 0, &xarr)); 967 PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr)); 968 /* lambda_s^T F_UU w_1 */ 969 PetscCall(TSComputeRHSHessianProductFunctionUU(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_sensip_col, ts->vecs_guu)); 970 if (quadts) { 971 /* R_UU w_1 */ 972 PetscCall(TSComputeRHSHessianProductFunctionUU(quadts, rk->stage_time, Y[i], NULL, ts->vec_sensip_col, ts->vecs_guu)); 973 } 974 if (ts->vecs_sensip) { 975 /* lambda_s^T F_UP w_2 */ 976 PetscCall(TSComputeRHSHessianProductFunctionUP(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_dir, ts->vecs_gup)); 977 if (quadts) { 978 /* R_UP w_2 */ 979 PetscCall(TSComputeRHSHessianProductFunctionUP(quadts, rk->stage_time, Y[i], NULL, ts->vec_sensip_col, ts->vecs_gup)); 980 } 981 } 982 if (ts->vecs_sensi2p) { 983 /* lambda_s^T F_PU w_1 */ 984 PetscCall(TSComputeRHSHessianProductFunctionPU(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_sensip_col, ts->vecs_gpu)); 985 /* lambda_s^T F_PP w_2 */ 986 PetscCall(TSComputeRHSHessianProductFunctionPP(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_dir, ts->vecs_gpp)); 987 if (b[i] && quadts) { 988 /* R_PU w_1 */ 989 PetscCall(TSComputeRHSHessianProductFunctionPU(quadts, rk->stage_time, Y[i], NULL, ts->vec_sensip_col, ts->vecs_gpu)); 990 /* R_PP w_2 */ 991 PetscCall(TSComputeRHSHessianProductFunctionPP(quadts, rk->stage_time, Y[i], NULL, ts->vec_dir, ts->vecs_gpp)); 992 } 993 } 994 PetscCall(VecResetArray(ts->vec_sensip_col)); 995 PetscCall(MatDenseRestoreColumn(rk->MatsFwdStageSensip[i], &xarr)); 996 997 for (nadj = 0; nadj < ts->numcost; nadj++) { 998 /* Stage values of lambda */ 999 if (b[i]) { 1000 /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */ 1001 PetscCall(VecCopy(ts->vecs_sensi2[nadj], VecsSensi2Temp[nadj])); 1002 PetscCall(VecMAXPY(VecsSensi2Temp[nadj], s - i - 1, w, &VecsDeltaLam2[nadj * s + i + 1])); 1003 PetscCall(MatMultTranspose(J, VecsSensi2Temp[nadj], VecsDeltaLam2[nadj * s + i])); 1004 PetscCall(VecScale(VecsDeltaLam2[nadj * s + i], -h * b[i])); 1005 PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h * b[i], ts->vecs_guu[nadj])); 1006 if (ts->vecs_sensip) PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h * b[i], ts->vecs_gup[nadj])); 1007 } else { 1008 /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */ 1009 PetscCall(VecSet(VecsDeltaLam2[nadj * s + i], 0)); 1010 PetscCall(VecMAXPY(VecsSensi2Temp[nadj], s - i - 1, w, &VecsDeltaLam2[nadj * s + i + 1])); 1011 PetscCall(MatMultTranspose(J, VecsSensi2Temp[nadj], VecsDeltaLam2[nadj * s + i])); 1012 PetscCall(VecScale(VecsDeltaLam2[nadj * s + i], -h)); 1013 PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h, ts->vecs_guu[nadj])); 1014 if (ts->vecs_sensip) PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h, ts->vecs_gup[nadj])); 1015 } 1016 if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */ 1017 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensi2Temp[nadj], VecDeltaMu2)); 1018 if (b[i]) { 1019 PetscCall(VecScale(VecDeltaMu2, -h * b[i])); 1020 PetscCall(VecAXPY(VecDeltaMu2, -h * b[i], ts->vecs_gpu[nadj])); 1021 PetscCall(VecAXPY(VecDeltaMu2, -h * b[i], ts->vecs_gpp[nadj])); 1022 } else { 1023 PetscCall(VecScale(VecDeltaMu2, -h)); 1024 PetscCall(VecAXPY(VecDeltaMu2, -h, ts->vecs_gpu[nadj])); 1025 PetscCall(VecAXPY(VecDeltaMu2, -h, ts->vecs_gpp[nadj])); 1026 } 1027 PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], 1, VecDeltaMu2)); /* update sensi2p for each stage */ 1028 } 1029 } 1030 } 1031 } 1032 1033 for (j = 0; j < s; j++) w[j] = 1.0; 1034 for (nadj = 0; nadj < ts->numcost; nadj++) { /* no need to do this for mu's */ 1035 PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s])); 1036 if (ts->vecs_sensi2) PetscCall(VecMAXPY(ts->vecs_sensi2[nadj], s, w, &VecsDeltaLam2[nadj * s])); 1037 } 1038 rk->status = TS_STEP_COMPLETE; 1039 PetscFunctionReturn(PETSC_SUCCESS); 1040 } 1041 1042 static PetscErrorCode TSAdjointReset_RK(TS ts) 1043 { 1044 TS_RK *rk = (TS_RK *)ts->data; 1045 RKTableau tab = rk->tableau; 1046 1047 PetscFunctionBegin; 1048 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &rk->VecsDeltaLam)); 1049 PetscCall(VecDestroyVecs(ts->numcost, &rk->VecsSensiTemp)); 1050 PetscCall(VecDestroy(&rk->VecDeltaMu)); 1051 PetscCall(VecDestroyVecs(tab->s * ts->numcost, &rk->VecsDeltaLam2)); 1052 PetscCall(VecDestroy(&rk->VecDeltaMu2)); 1053 PetscCall(VecDestroyVecs(ts->numcost, &rk->VecsSensi2Temp)); 1054 PetscFunctionReturn(PETSC_SUCCESS); 1055 } 1056 1057 static PetscErrorCode TSInterpolate_RK(TS ts, PetscReal itime, Vec X) 1058 { 1059 TS_RK *rk = (TS_RK *)ts->data; 1060 PetscInt s = rk->tableau->s, p = rk->tableau->p, i, j; 1061 PetscReal h; 1062 PetscReal tt, t; 1063 PetscScalar *b; 1064 const PetscReal *B = rk->tableau->binterp; 1065 1066 PetscFunctionBegin; 1067 PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRK %s does not have an interpolation formula", rk->tableau->name); 1068 1069 switch (rk->status) { 1070 case TS_STEP_INCOMPLETE: 1071 case TS_STEP_PENDING: 1072 h = ts->time_step; 1073 t = (itime - ts->ptime) / h; 1074 break; 1075 case TS_STEP_COMPLETE: 1076 h = ts->ptime - ts->ptime_prev; 1077 t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */ 1078 break; 1079 default: 1080 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus"); 1081 } 1082 PetscCall(PetscMalloc1(s, &b)); 1083 for (i = 0; i < s; i++) b[i] = 0; 1084 for (j = 0, tt = t; j < p; j++, tt *= t) { 1085 for (i = 0; i < s; i++) b[i] += h * B[i * p + j] * tt; 1086 } 1087 PetscCall(VecCopy(rk->Y[0], X)); 1088 PetscCall(VecMAXPY(X, s, b, rk->YdotRHS)); 1089 PetscCall(PetscFree(b)); 1090 PetscFunctionReturn(PETSC_SUCCESS); 1091 } 1092 1093 /*------------------------------------------------------------*/ 1094 1095 static PetscErrorCode TSRKTableauReset(TS ts) 1096 { 1097 TS_RK *rk = (TS_RK *)ts->data; 1098 RKTableau tab = rk->tableau; 1099 1100 PetscFunctionBegin; 1101 if (!tab) PetscFunctionReturn(PETSC_SUCCESS); 1102 PetscCall(PetscFree(rk->work)); 1103 PetscCall(VecDestroyVecs(tab->s, &rk->Y)); 1104 PetscCall(VecDestroyVecs(tab->s, &rk->YdotRHS)); 1105 PetscFunctionReturn(PETSC_SUCCESS); 1106 } 1107 1108 static PetscErrorCode TSReset_RK(TS ts) 1109 { 1110 PetscFunctionBegin; 1111 PetscCall(TSRKTableauReset(ts)); 1112 if (ts->use_splitrhsfunction) { 1113 PetscTryMethod(ts, "TSReset_RK_MultirateSplit_C", (TS), (ts)); 1114 } else { 1115 PetscTryMethod(ts, "TSReset_RK_MultirateNonsplit_C", (TS), (ts)); 1116 } 1117 PetscFunctionReturn(PETSC_SUCCESS); 1118 } 1119 1120 static PetscErrorCode DMCoarsenHook_TSRK(DM fine, DM coarse, void *ctx) 1121 { 1122 PetscFunctionBegin; 1123 PetscFunctionReturn(PETSC_SUCCESS); 1124 } 1125 1126 static PetscErrorCode DMRestrictHook_TSRK(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx) 1127 { 1128 PetscFunctionBegin; 1129 PetscFunctionReturn(PETSC_SUCCESS); 1130 } 1131 1132 static PetscErrorCode DMSubDomainHook_TSRK(DM dm, DM subdm, void *ctx) 1133 { 1134 PetscFunctionBegin; 1135 PetscFunctionReturn(PETSC_SUCCESS); 1136 } 1137 1138 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx) 1139 { 1140 PetscFunctionBegin; 1141 PetscFunctionReturn(PETSC_SUCCESS); 1142 } 1143 1144 static PetscErrorCode TSRKTableauSetUp(TS ts) 1145 { 1146 TS_RK *rk = (TS_RK *)ts->data; 1147 RKTableau tab = rk->tableau; 1148 1149 PetscFunctionBegin; 1150 PetscCall(PetscMalloc1(tab->s, &rk->work)); 1151 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &rk->Y)); 1152 PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &rk->YdotRHS)); 1153 rk->newtableau = PETSC_TRUE; 1154 PetscFunctionReturn(PETSC_SUCCESS); 1155 } 1156 1157 static PetscErrorCode TSSetUp_RK(TS ts) 1158 { 1159 TS quadts = ts->quadraturets; 1160 DM dm; 1161 1162 PetscFunctionBegin; 1163 PetscCall(TSCheckImplicitTerm(ts)); 1164 PetscCall(TSRKTableauSetUp(ts)); 1165 if (quadts && ts->costintegralfwd) { 1166 Mat Jquad; 1167 PetscCall(TSGetRHSJacobian(quadts, &Jquad, NULL, NULL, NULL)); 1168 } 1169 PetscCall(TSGetDM(ts, &dm)); 1170 PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRK, DMRestrictHook_TSRK, ts)); 1171 PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRK, DMSubDomainRestrictHook_TSRK, ts)); 1172 if (ts->use_splitrhsfunction) { 1173 PetscTryMethod(ts, "TSSetUp_RK_MultirateSplit_C", (TS), (ts)); 1174 } else { 1175 PetscTryMethod(ts, "TSSetUp_RK_MultirateNonsplit_C", (TS), (ts)); 1176 } 1177 PetscFunctionReturn(PETSC_SUCCESS); 1178 } 1179 1180 static PetscErrorCode TSSetFromOptions_RK(TS ts, PetscOptionItems *PetscOptionsObject) 1181 { 1182 TS_RK *rk = (TS_RK *)ts->data; 1183 1184 PetscFunctionBegin; 1185 PetscOptionsHeadBegin(PetscOptionsObject, "RK ODE solver options"); 1186 { 1187 RKTableauLink link; 1188 PetscInt count, choice; 1189 PetscBool flg, use_multirate = PETSC_FALSE; 1190 const char **namelist; 1191 1192 for (link = RKTableauList, count = 0; link; link = link->next, count++) 1193 ; 1194 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1195 for (link = RKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1196 PetscCall(PetscOptionsBool("-ts_rk_multirate", "Use interpolation-based multirate RK method", "TSRKSetMultirate", rk->use_multirate, &use_multirate, &flg)); 1197 if (flg) PetscCall(TSRKSetMultirate(ts, use_multirate)); 1198 PetscCall(PetscOptionsEList("-ts_rk_type", "Family of RK method", "TSRKSetType", (const char *const *)namelist, count, rk->tableau->name, &choice, &flg)); 1199 if (flg) PetscCall(TSRKSetType(ts, namelist[choice])); 1200 PetscCall(PetscFree(namelist)); 1201 } 1202 PetscOptionsHeadEnd(); 1203 PetscOptionsBegin(PetscObjectComm((PetscObject)ts), NULL, "Multirate methods options", ""); 1204 PetscCall(PetscOptionsInt("-ts_rk_dtratio", "time step ratio between slow and fast", "", rk->dtratio, &rk->dtratio, NULL)); 1205 PetscOptionsEnd(); 1206 PetscFunctionReturn(PETSC_SUCCESS); 1207 } 1208 1209 static PetscErrorCode TSView_RK(TS ts, PetscViewer viewer) 1210 { 1211 TS_RK *rk = (TS_RK *)ts->data; 1212 PetscBool iascii; 1213 1214 PetscFunctionBegin; 1215 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1216 if (iascii) { 1217 RKTableau tab = rk->tableau; 1218 TSRKType rktype; 1219 const PetscReal *c; 1220 PetscInt s; 1221 char buf[512]; 1222 PetscBool FSAL; 1223 1224 PetscCall(TSRKGetType(ts, &rktype)); 1225 PetscCall(TSRKGetTableau(ts, &s, NULL, NULL, &c, NULL, NULL, NULL, &FSAL)); 1226 PetscCall(PetscViewerASCIIPrintf(viewer, " RK type %s\n", rktype)); 1227 PetscCall(PetscViewerASCIIPrintf(viewer, " Order: %" PetscInt_FMT "\n", tab->order)); 1228 PetscCall(PetscViewerASCIIPrintf(viewer, " FSAL property: %s\n", FSAL ? "yes" : "no")); 1229 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", s, c)); 1230 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf)); 1231 } 1232 PetscFunctionReturn(PETSC_SUCCESS); 1233 } 1234 1235 static PetscErrorCode TSLoad_RK(TS ts, PetscViewer viewer) 1236 { 1237 TSAdapt adapt; 1238 1239 PetscFunctionBegin; 1240 PetscCall(TSGetAdapt(ts, &adapt)); 1241 PetscCall(TSAdaptLoad(adapt, viewer)); 1242 PetscFunctionReturn(PETSC_SUCCESS); 1243 } 1244 1245 /*@ 1246 TSRKGetOrder - Get the order of the `TSRK` scheme 1247 1248 Not Collective 1249 1250 Input Parameter: 1251 . ts - timestepping context 1252 1253 Output Parameter: 1254 . order - order of `TSRK` scheme 1255 1256 Level: intermediate 1257 1258 .seealso: [](ch_ts), `TSRK`, `TSRKGetType()` 1259 @*/ 1260 PetscErrorCode TSRKGetOrder(TS ts, PetscInt *order) 1261 { 1262 PetscFunctionBegin; 1263 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1264 PetscAssertPointer(order, 2); 1265 PetscUseMethod(ts, "TSRKGetOrder_C", (TS, PetscInt *), (ts, order)); 1266 PetscFunctionReturn(PETSC_SUCCESS); 1267 } 1268 1269 /*@C 1270 TSRKSetType - Set the type of the `TSRK` scheme 1271 1272 Logically Collective 1273 1274 Input Parameters: 1275 + ts - timestepping context 1276 - rktype - type of `TSRK` scheme 1277 1278 Options Database Key: 1279 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1280 1281 Level: intermediate 1282 1283 .seealso: [](ch_ts), `TSRKGetType()`, `TSRK`, `TSRKType`, `TSRK1FE`, `TSRK2A`, `TSRK2B`, `TSRK3`, `TSRK3BS`, `TSRK4`, `TSRK5F`, `TSRK5DP`, `TSRK5BS`, `TSRK6VR`, `TSRK7VR`, `TSRK8VR` 1284 @*/ 1285 PetscErrorCode TSRKSetType(TS ts, TSRKType rktype) 1286 { 1287 PetscFunctionBegin; 1288 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1289 PetscAssertPointer(rktype, 2); 1290 PetscTryMethod(ts, "TSRKSetType_C", (TS, TSRKType), (ts, rktype)); 1291 PetscFunctionReturn(PETSC_SUCCESS); 1292 } 1293 1294 /*@C 1295 TSRKGetType - Get the type of `TSRK` scheme 1296 1297 Not Collective 1298 1299 Input Parameter: 1300 . ts - timestepping context 1301 1302 Output Parameter: 1303 . rktype - type of `TSRK`-scheme 1304 1305 Level: intermediate 1306 1307 .seealso: [](ch_ts), `TSRKSetType()` 1308 @*/ 1309 PetscErrorCode TSRKGetType(TS ts, TSRKType *rktype) 1310 { 1311 PetscFunctionBegin; 1312 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1313 PetscUseMethod(ts, "TSRKGetType_C", (TS, TSRKType *), (ts, rktype)); 1314 PetscFunctionReturn(PETSC_SUCCESS); 1315 } 1316 1317 static PetscErrorCode TSRKGetOrder_RK(TS ts, PetscInt *order) 1318 { 1319 TS_RK *rk = (TS_RK *)ts->data; 1320 1321 PetscFunctionBegin; 1322 *order = rk->tableau->order; 1323 PetscFunctionReturn(PETSC_SUCCESS); 1324 } 1325 1326 static PetscErrorCode TSRKGetType_RK(TS ts, TSRKType *rktype) 1327 { 1328 TS_RK *rk = (TS_RK *)ts->data; 1329 1330 PetscFunctionBegin; 1331 *rktype = rk->tableau->name; 1332 PetscFunctionReturn(PETSC_SUCCESS); 1333 } 1334 1335 static PetscErrorCode TSRKSetType_RK(TS ts, TSRKType rktype) 1336 { 1337 TS_RK *rk = (TS_RK *)ts->data; 1338 PetscBool match; 1339 RKTableauLink link; 1340 1341 PetscFunctionBegin; 1342 if (rk->tableau) { 1343 PetscCall(PetscStrcmp(rk->tableau->name, rktype, &match)); 1344 if (match) PetscFunctionReturn(PETSC_SUCCESS); 1345 } 1346 for (link = RKTableauList; link; link = link->next) { 1347 PetscCall(PetscStrcmp(link->tab.name, rktype, &match)); 1348 if (match) { 1349 if (ts->setupcalled) PetscCall(TSRKTableauReset(ts)); 1350 rk->tableau = &link->tab; 1351 if (ts->setupcalled) PetscCall(TSRKTableauSetUp(ts)); 1352 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1353 PetscFunctionReturn(PETSC_SUCCESS); 1354 } 1355 } 1356 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rktype); 1357 } 1358 1359 static PetscErrorCode TSGetStages_RK(TS ts, PetscInt *ns, Vec **Y) 1360 { 1361 TS_RK *rk = (TS_RK *)ts->data; 1362 1363 PetscFunctionBegin; 1364 if (ns) *ns = rk->tableau->s; 1365 if (Y) *Y = rk->Y; 1366 PetscFunctionReturn(PETSC_SUCCESS); 1367 } 1368 1369 static PetscErrorCode TSDestroy_RK(TS ts) 1370 { 1371 PetscFunctionBegin; 1372 PetscCall(TSReset_RK(ts)); 1373 if (ts->dm) { 1374 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRK, DMRestrictHook_TSRK, ts)); 1375 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRK, DMSubDomainRestrictHook_TSRK, ts)); 1376 } 1377 PetscCall(PetscFree(ts->data)); 1378 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetOrder_C", NULL)); 1379 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetType_C", NULL)); 1380 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetType_C", NULL)); 1381 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetTableau_C", NULL)); 1382 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetMultirate_C", NULL)); 1383 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetMultirate_C", NULL)); 1384 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", NULL)); 1385 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", NULL)); 1386 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", NULL)); 1387 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", NULL)); 1388 PetscFunctionReturn(PETSC_SUCCESS); 1389 } 1390 1391 /* 1392 This defines the nonlinear equation that is to be solved with SNES 1393 We do not need to solve the equation; we just use SNES to approximate the Jacobian 1394 */ 1395 static PetscErrorCode SNESTSFormFunction_RK(SNES snes, Vec x, Vec y, TS ts) 1396 { 1397 TS_RK *rk = (TS_RK *)ts->data; 1398 DM dm, dmsave; 1399 1400 PetscFunctionBegin; 1401 PetscCall(SNESGetDM(snes, &dm)); 1402 /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */ 1403 dmsave = ts->dm; 1404 ts->dm = dm; 1405 PetscCall(TSComputeRHSFunction(ts, rk->stage_time, x, y)); 1406 ts->dm = dmsave; 1407 PetscFunctionReturn(PETSC_SUCCESS); 1408 } 1409 1410 static PetscErrorCode SNESTSFormJacobian_RK(SNES snes, Vec x, Mat A, Mat B, TS ts) 1411 { 1412 TS_RK *rk = (TS_RK *)ts->data; 1413 DM dm, dmsave; 1414 1415 PetscFunctionBegin; 1416 PetscCall(SNESGetDM(snes, &dm)); 1417 dmsave = ts->dm; 1418 ts->dm = dm; 1419 PetscCall(TSComputeRHSJacobian(ts, rk->stage_time, x, A, B)); 1420 ts->dm = dmsave; 1421 PetscFunctionReturn(PETSC_SUCCESS); 1422 } 1423 1424 /*@C 1425 TSRKSetMultirate - Use the interpolation-based multirate `TSRK` method 1426 1427 Logically Collective 1428 1429 Input Parameters: 1430 + ts - timestepping context 1431 - 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 1432 1433 Options Database Key: 1434 . -ts_rk_multirate - <true,false> 1435 1436 Level: intermediate 1437 1438 Note: 1439 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). 1440 1441 .seealso: [](ch_ts), `TSRK`, `TSRKGetMultirate()` 1442 @*/ 1443 PetscErrorCode TSRKSetMultirate(TS ts, PetscBool use_multirate) 1444 { 1445 PetscFunctionBegin; 1446 PetscTryMethod(ts, "TSRKSetMultirate_C", (TS, PetscBool), (ts, use_multirate)); 1447 PetscFunctionReturn(PETSC_SUCCESS); 1448 } 1449 1450 /*@C 1451 TSRKGetMultirate - Gets whether to use the interpolation-based multirate `TSRK` method 1452 1453 Not Collective 1454 1455 Input Parameter: 1456 . ts - timestepping context 1457 1458 Output Parameter: 1459 . use_multirate - `PETSC_TRUE` if the multirate RK method is enabled, `PETSC_FALSE` otherwise 1460 1461 Level: intermediate 1462 1463 .seealso: [](ch_ts), `TSRK`, `TSRKSetMultirate()` 1464 @*/ 1465 PetscErrorCode TSRKGetMultirate(TS ts, PetscBool *use_multirate) 1466 { 1467 PetscFunctionBegin; 1468 PetscUseMethod(ts, "TSRKGetMultirate_C", (TS, PetscBool *), (ts, use_multirate)); 1469 PetscFunctionReturn(PETSC_SUCCESS); 1470 } 1471 1472 /*MC 1473 TSRK - ODE and DAE solver using Runge-Kutta schemes 1474 1475 The user should provide the right hand side of the equation 1476 using `TSSetRHSFunction()`. 1477 1478 Level: beginner 1479 1480 Notes: 1481 The default is `TSRK3BS`, it can be changed with `TSRKSetType()` or -ts_rk_type 1482 1483 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSRK`, `TSSetType()`, `TSRKSetType()`, `TSRKGetType()`, `TSRK2D`, `TSRK2E`, `TSRK3`, 1484 `TSRK4`, `TSRK5`, `TSRKPRSSP2`, `TSRKBPR3`, `TSRKType`, `TSRKRegister()`, `TSRKSetMultirate()`, `TSRKGetMultirate()`, `TSType` 1485 M*/ 1486 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1487 { 1488 TS_RK *rk; 1489 1490 PetscFunctionBegin; 1491 PetscCall(TSRKInitializePackage()); 1492 1493 ts->ops->reset = TSReset_RK; 1494 ts->ops->destroy = TSDestroy_RK; 1495 ts->ops->view = TSView_RK; 1496 ts->ops->load = TSLoad_RK; 1497 ts->ops->setup = TSSetUp_RK; 1498 ts->ops->interpolate = TSInterpolate_RK; 1499 ts->ops->step = TSStep_RK; 1500 ts->ops->evaluatestep = TSEvaluateStep_RK; 1501 ts->ops->rollback = TSRollBack_RK; 1502 ts->ops->setfromoptions = TSSetFromOptions_RK; 1503 ts->ops->getstages = TSGetStages_RK; 1504 1505 ts->ops->snesfunction = SNESTSFormFunction_RK; 1506 ts->ops->snesjacobian = SNESTSFormJacobian_RK; 1507 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1508 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1509 ts->ops->adjointstep = TSAdjointStep_RK; 1510 ts->ops->adjointreset = TSAdjointReset_RK; 1511 1512 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1513 ts->ops->forwardsetup = TSForwardSetUp_RK; 1514 ts->ops->forwardreset = TSForwardReset_RK; 1515 ts->ops->forwardstep = TSForwardStep_RK; 1516 ts->ops->forwardgetstages = TSForwardGetStages_RK; 1517 1518 PetscCall(PetscNew(&rk)); 1519 ts->data = (void *)rk; 1520 1521 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetOrder_C", TSRKGetOrder_RK)); 1522 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetType_C", TSRKGetType_RK)); 1523 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetType_C", TSRKSetType_RK)); 1524 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetTableau_C", TSRKGetTableau_RK)); 1525 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetMultirate_C", TSRKSetMultirate_RK)); 1526 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetMultirate_C", TSRKGetMultirate_RK)); 1527 1528 PetscCall(TSRKSetType(ts, TSRKDefault)); 1529 rk->dtratio = 1; 1530 PetscFunctionReturn(PETSC_SUCCESS); 1531 } 1532