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