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