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