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