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