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