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