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