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