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