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