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