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