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 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_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 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 PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecDeltaMu)); 936 if (b[i]) { 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 PetscCall(PetscMalloc1(count, (char ***)&namelist)); 1182 for (link = RKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name; 1183 PetscCall(PetscOptionsBool("-ts_rk_multirate", "Use interpolation-based multirate RK method", "TSRKSetMultirate", rk->use_multirate, &use_multirate, &flg)); 1184 if (flg) PetscCall(TSRKSetMultirate(ts, use_multirate)); 1185 PetscCall(PetscOptionsEList("-ts_rk_type", "Family of RK method", "TSRKSetType", (const char *const *)namelist, count, rk->tableau->name, &choice, &flg)); 1186 if (flg) PetscCall(TSRKSetType(ts, namelist[choice])); 1187 PetscCall(PetscFree(namelist)); 1188 } 1189 PetscOptionsHeadEnd(); 1190 PetscOptionsBegin(PetscObjectComm((PetscObject)ts), NULL, "Multirate methods options", ""); 1191 PetscCall(PetscOptionsInt("-ts_rk_dtratio", "time step ratio between slow and fast", "", rk->dtratio, &rk->dtratio, NULL)); 1192 PetscOptionsEnd(); 1193 PetscFunctionReturn(PETSC_SUCCESS); 1194 } 1195 1196 static PetscErrorCode TSView_RK(TS ts, PetscViewer viewer) 1197 { 1198 TS_RK *rk = (TS_RK *)ts->data; 1199 PetscBool iascii; 1200 1201 PetscFunctionBegin; 1202 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1203 if (iascii) { 1204 RKTableau tab = rk->tableau; 1205 TSRKType rktype; 1206 const PetscReal *c; 1207 PetscInt s; 1208 char buf[512]; 1209 PetscBool FSAL; 1210 1211 PetscCall(TSRKGetType(ts, &rktype)); 1212 PetscCall(TSRKGetTableau(ts, &s, NULL, NULL, &c, NULL, NULL, NULL, &FSAL)); 1213 PetscCall(PetscViewerASCIIPrintf(viewer, " RK type %s\n", rktype)); 1214 PetscCall(PetscViewerASCIIPrintf(viewer, " Order: %" PetscInt_FMT "\n", tab->order)); 1215 PetscCall(PetscViewerASCIIPrintf(viewer, " FSAL property: %s\n", FSAL ? "yes" : "no")); 1216 PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", s, c)); 1217 PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf)); 1218 } 1219 PetscFunctionReturn(PETSC_SUCCESS); 1220 } 1221 1222 static PetscErrorCode TSLoad_RK(TS ts, PetscViewer viewer) 1223 { 1224 TSAdapt adapt; 1225 1226 PetscFunctionBegin; 1227 PetscCall(TSGetAdapt(ts, &adapt)); 1228 PetscCall(TSAdaptLoad(adapt, viewer)); 1229 PetscFunctionReturn(PETSC_SUCCESS); 1230 } 1231 1232 /*@ 1233 TSRKGetOrder - Get the order of the `TSRK` scheme 1234 1235 Not Collective 1236 1237 Input Parameter: 1238 . ts - timestepping context 1239 1240 Output Parameter: 1241 . order - order of `TSRK` scheme 1242 1243 Level: intermediate 1244 1245 .seealso: [](ch_ts), `TSRK`, `TSRKGetType()` 1246 @*/ 1247 PetscErrorCode TSRKGetOrder(TS ts, PetscInt *order) 1248 { 1249 PetscFunctionBegin; 1250 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1251 PetscAssertPointer(order, 2); 1252 PetscUseMethod(ts, "TSRKGetOrder_C", (TS, PetscInt *), (ts, order)); 1253 PetscFunctionReturn(PETSC_SUCCESS); 1254 } 1255 1256 /*@ 1257 TSRKSetType - Set the type of the `TSRK` scheme 1258 1259 Logically Collective 1260 1261 Input Parameters: 1262 + ts - timestepping context 1263 - rktype - type of `TSRK` scheme 1264 1265 Options Database Key: 1266 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1267 1268 Level: intermediate 1269 1270 .seealso: [](ch_ts), `TSRKGetType()`, `TSRK`, `TSRKType`, `TSRK1FE`, `TSRK2A`, `TSRK2B`, `TSRK3`, `TSRK3BS`, `TSRK4`, `TSRK5F`, `TSRK5DP`, `TSRK5BS`, `TSRK6VR`, `TSRK7VR`, `TSRK8VR` 1271 @*/ 1272 PetscErrorCode TSRKSetType(TS ts, TSRKType rktype) 1273 { 1274 PetscFunctionBegin; 1275 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1276 PetscAssertPointer(rktype, 2); 1277 PetscTryMethod(ts, "TSRKSetType_C", (TS, TSRKType), (ts, rktype)); 1278 PetscFunctionReturn(PETSC_SUCCESS); 1279 } 1280 1281 /*@ 1282 TSRKGetType - Get the type of `TSRK` scheme 1283 1284 Not Collective 1285 1286 Input Parameter: 1287 . ts - timestepping context 1288 1289 Output Parameter: 1290 . rktype - type of `TSRK`-scheme 1291 1292 Level: intermediate 1293 1294 .seealso: [](ch_ts), `TSRKSetType()` 1295 @*/ 1296 PetscErrorCode TSRKGetType(TS ts, TSRKType *rktype) 1297 { 1298 PetscFunctionBegin; 1299 PetscValidHeaderSpecific(ts, TS_CLASSID, 1); 1300 PetscUseMethod(ts, "TSRKGetType_C", (TS, TSRKType *), (ts, rktype)); 1301 PetscFunctionReturn(PETSC_SUCCESS); 1302 } 1303 1304 static PetscErrorCode TSRKGetOrder_RK(TS ts, PetscInt *order) 1305 { 1306 TS_RK *rk = (TS_RK *)ts->data; 1307 1308 PetscFunctionBegin; 1309 *order = rk->tableau->order; 1310 PetscFunctionReturn(PETSC_SUCCESS); 1311 } 1312 1313 static PetscErrorCode TSRKGetType_RK(TS ts, TSRKType *rktype) 1314 { 1315 TS_RK *rk = (TS_RK *)ts->data; 1316 1317 PetscFunctionBegin; 1318 *rktype = rk->tableau->name; 1319 PetscFunctionReturn(PETSC_SUCCESS); 1320 } 1321 1322 static PetscErrorCode TSRKSetType_RK(TS ts, TSRKType rktype) 1323 { 1324 TS_RK *rk = (TS_RK *)ts->data; 1325 PetscBool match; 1326 RKTableauLink link; 1327 1328 PetscFunctionBegin; 1329 if (rk->tableau) { 1330 PetscCall(PetscStrcmp(rk->tableau->name, rktype, &match)); 1331 if (match) PetscFunctionReturn(PETSC_SUCCESS); 1332 } 1333 for (link = RKTableauList; link; link = link->next) { 1334 PetscCall(PetscStrcmp(link->tab.name, rktype, &match)); 1335 if (match) { 1336 if (ts->setupcalled) PetscCall(TSRKTableauReset(ts)); 1337 rk->tableau = &link->tab; 1338 if (ts->setupcalled) PetscCall(TSRKTableauSetUp(ts)); 1339 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1340 PetscFunctionReturn(PETSC_SUCCESS); 1341 } 1342 } 1343 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rktype); 1344 } 1345 1346 static PetscErrorCode TSGetStages_RK(TS ts, PetscInt *ns, Vec **Y) 1347 { 1348 TS_RK *rk = (TS_RK *)ts->data; 1349 1350 PetscFunctionBegin; 1351 if (ns) *ns = rk->tableau->s; 1352 if (Y) *Y = rk->Y; 1353 PetscFunctionReturn(PETSC_SUCCESS); 1354 } 1355 1356 static PetscErrorCode TSDestroy_RK(TS ts) 1357 { 1358 PetscFunctionBegin; 1359 PetscCall(TSReset_RK(ts)); 1360 if (ts->dm) { 1361 PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRK, DMRestrictHook_TSRK, ts)); 1362 PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRK, DMSubDomainRestrictHook_TSRK, ts)); 1363 } 1364 PetscCall(PetscFree(ts->data)); 1365 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetOrder_C", NULL)); 1366 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetType_C", NULL)); 1367 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetType_C", NULL)); 1368 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetTableau_C", NULL)); 1369 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetMultirate_C", NULL)); 1370 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetMultirate_C", NULL)); 1371 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", NULL)); 1372 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", NULL)); 1373 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", NULL)); 1374 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", NULL)); 1375 PetscFunctionReturn(PETSC_SUCCESS); 1376 } 1377 1378 /* 1379 This defines the nonlinear equation that is to be solved with SNES 1380 We do not need to solve the equation; we just use SNES to approximate the Jacobian 1381 */ 1382 static PetscErrorCode SNESTSFormFunction_RK(SNES snes, Vec x, Vec y, TS ts) 1383 { 1384 TS_RK *rk = (TS_RK *)ts->data; 1385 DM dm, dmsave; 1386 1387 PetscFunctionBegin; 1388 PetscCall(SNESGetDM(snes, &dm)); 1389 /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */ 1390 dmsave = ts->dm; 1391 ts->dm = dm; 1392 PetscCall(TSComputeRHSFunction(ts, rk->stage_time, x, y)); 1393 ts->dm = dmsave; 1394 PetscFunctionReturn(PETSC_SUCCESS); 1395 } 1396 1397 static PetscErrorCode SNESTSFormJacobian_RK(SNES snes, Vec x, Mat A, Mat B, TS ts) 1398 { 1399 TS_RK *rk = (TS_RK *)ts->data; 1400 DM dm, dmsave; 1401 1402 PetscFunctionBegin; 1403 PetscCall(SNESGetDM(snes, &dm)); 1404 dmsave = ts->dm; 1405 ts->dm = dm; 1406 PetscCall(TSComputeRHSJacobian(ts, rk->stage_time, x, A, B)); 1407 ts->dm = dmsave; 1408 PetscFunctionReturn(PETSC_SUCCESS); 1409 } 1410 1411 /*@ 1412 TSRKSetMultirate - Use the interpolation-based multirate `TSRK` method 1413 1414 Logically Collective 1415 1416 Input Parameters: 1417 + ts - timestepping context 1418 - 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 1419 1420 Options Database Key: 1421 . -ts_rk_multirate - <true,false> 1422 1423 Level: intermediate 1424 1425 Note: 1426 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). 1427 1428 .seealso: [](ch_ts), `TSRK`, `TSRKGetMultirate()` 1429 @*/ 1430 PetscErrorCode TSRKSetMultirate(TS ts, PetscBool use_multirate) 1431 { 1432 PetscFunctionBegin; 1433 PetscTryMethod(ts, "TSRKSetMultirate_C", (TS, PetscBool), (ts, use_multirate)); 1434 PetscFunctionReturn(PETSC_SUCCESS); 1435 } 1436 1437 /*@ 1438 TSRKGetMultirate - Gets whether to use the interpolation-based multirate `TSRK` method 1439 1440 Not Collective 1441 1442 Input Parameter: 1443 . ts - timestepping context 1444 1445 Output Parameter: 1446 . use_multirate - `PETSC_TRUE` if the multirate RK method is enabled, `PETSC_FALSE` otherwise 1447 1448 Level: intermediate 1449 1450 .seealso: [](ch_ts), `TSRK`, `TSRKSetMultirate()` 1451 @*/ 1452 PetscErrorCode TSRKGetMultirate(TS ts, PetscBool *use_multirate) 1453 { 1454 PetscFunctionBegin; 1455 PetscUseMethod(ts, "TSRKGetMultirate_C", (TS, PetscBool *), (ts, use_multirate)); 1456 PetscFunctionReturn(PETSC_SUCCESS); 1457 } 1458 1459 /*MC 1460 TSRK - ODE and DAE solver using Runge-Kutta schemes 1461 1462 The user should provide the right-hand side of the equation 1463 using `TSSetRHSFunction()`. 1464 1465 Level: beginner 1466 1467 Notes: 1468 The default is `TSRK3BS`, it can be changed with `TSRKSetType()` or -ts_rk_type 1469 1470 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSRK`, `TSSetType()`, `TSRKSetType()`, `TSRKGetType()`, `TSRK2D`, `TSRK2E`, `TSRK3`, 1471 `TSRK4`, `TSRK5`, `TSRKPRSSP2`, `TSRKBPR3`, `TSRKType`, `TSRKRegister()`, `TSRKSetMultirate()`, `TSRKGetMultirate()`, `TSType` 1472 M*/ 1473 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1474 { 1475 TS_RK *rk; 1476 1477 PetscFunctionBegin; 1478 PetscCall(TSRKInitializePackage()); 1479 1480 ts->ops->reset = TSReset_RK; 1481 ts->ops->destroy = TSDestroy_RK; 1482 ts->ops->view = TSView_RK; 1483 ts->ops->load = TSLoad_RK; 1484 ts->ops->setup = TSSetUp_RK; 1485 ts->ops->interpolate = TSInterpolate_RK; 1486 ts->ops->step = TSStep_RK; 1487 ts->ops->evaluatestep = TSEvaluateStep_RK; 1488 ts->ops->rollback = TSRollBack_RK; 1489 ts->ops->setfromoptions = TSSetFromOptions_RK; 1490 ts->ops->getstages = TSGetStages_RK; 1491 1492 ts->ops->snesfunction = SNESTSFormFunction_RK; 1493 ts->ops->snesjacobian = SNESTSFormJacobian_RK; 1494 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1495 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1496 ts->ops->adjointstep = TSAdjointStep_RK; 1497 ts->ops->adjointreset = TSAdjointReset_RK; 1498 1499 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1500 ts->ops->forwardsetup = TSForwardSetUp_RK; 1501 ts->ops->forwardreset = TSForwardReset_RK; 1502 ts->ops->forwardstep = TSForwardStep_RK; 1503 ts->ops->forwardgetstages = TSForwardGetStages_RK; 1504 1505 PetscCall(PetscNew(&rk)); 1506 ts->data = (void *)rk; 1507 1508 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetOrder_C", TSRKGetOrder_RK)); 1509 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetType_C", TSRKGetType_RK)); 1510 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetType_C", TSRKSetType_RK)); 1511 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetTableau_C", TSRKGetTableau_RK)); 1512 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetMultirate_C", TSRKSetMultirate_RK)); 1513 PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetMultirate_C", TSRKGetMultirate_RK)); 1514 1515 PetscCall(TSRKSetType(ts, TSRKDefault)); 1516 rk->dtratio = 1; 1517 PetscFunctionReturn(PETSC_SUCCESS); 1518 } 1519