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); 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,Jpre,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; 829 PetscErrorCode ierr; 830 831 PetscFunctionBegin; 832 rk->status = TS_STEP_INCOMPLETE; 833 834 ierr = TSGetRHSJacobian(ts,&J,&Jpre,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 rk->stage_time = t + h*(1.0-c[i]); 844 ierr = TSComputeSNESJacobian(ts,Y[i],J,Jpre);CHKERRQ(ierr); 845 if (quadts) { 846 ierr = TSComputeRHSJacobian(quadts,rk->stage_time,Y[i],Jquad,Jquad);CHKERRQ(ierr); /* get r_u^T */ 847 } 848 if (ts->vecs_sensip) { 849 ierr = TSComputeRHSJacobianP(ts,rk->stage_time,Y[i],ts->Jacprhs);CHKERRQ(ierr); /* get f_p */ 850 if (quadts) { 851 ierr = TSComputeRHSJacobianP(quadts,rk->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,rk->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,rk->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,rk->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,rk->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,rk->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,rk->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,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gpu);CHKERRQ(ierr); 929 /* R_PP w_2 */ 930 ierr = TSComputeRHSHessianProductFunctionPP(quadts,rk->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 /*@ 1202 TSRKGetOrder - Get the order of RK scheme 1203 1204 Not collective 1205 1206 Input Parameter: 1207 . ts - timestepping context 1208 1209 Output Parameter: 1210 . order - order of RK-scheme 1211 1212 Level: intermediate 1213 1214 .seealso: TSRKGetType() 1215 @*/ 1216 PetscErrorCode TSRKGetOrder(TS ts,PetscInt *order) 1217 { 1218 PetscErrorCode ierr; 1219 1220 PetscFunctionBegin; 1221 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1222 PetscValidIntPointer(order,2); 1223 ierr = PetscUseMethod(ts,"TSRKGetOrder_C",(TS,PetscInt*),(ts,order));CHKERRQ(ierr); 1224 PetscFunctionReturn(0); 1225 } 1226 1227 /*@C 1228 TSRKSetType - Set the type of RK scheme 1229 1230 Logically collective 1231 1232 Input Parameter: 1233 + ts - timestepping context 1234 - rktype - type of RK-scheme 1235 1236 Options Database: 1237 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1238 1239 Level: intermediate 1240 1241 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS, TSRK6VR, TSRK7VR, TSRK8VR 1242 @*/ 1243 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1244 { 1245 PetscErrorCode ierr; 1246 1247 PetscFunctionBegin; 1248 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1249 PetscValidCharPointer(rktype,2); 1250 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 1251 PetscFunctionReturn(0); 1252 } 1253 1254 /*@C 1255 TSRKGetType - Get the type of RK scheme 1256 1257 Not collective 1258 1259 Input Parameter: 1260 . ts - timestepping context 1261 1262 Output Parameter: 1263 . rktype - type of RK-scheme 1264 1265 Level: intermediate 1266 1267 .seealso: TSRKSetType() 1268 @*/ 1269 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1270 { 1271 PetscErrorCode ierr; 1272 1273 PetscFunctionBegin; 1274 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1275 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 1276 PetscFunctionReturn(0); 1277 } 1278 1279 static PetscErrorCode TSRKGetOrder_RK(TS ts,PetscInt *order) 1280 { 1281 TS_RK *rk = (TS_RK*)ts->data; 1282 1283 PetscFunctionBegin; 1284 *order = rk->tableau->order; 1285 PetscFunctionReturn(0); 1286 } 1287 1288 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1289 { 1290 TS_RK *rk = (TS_RK*)ts->data; 1291 1292 PetscFunctionBegin; 1293 *rktype = rk->tableau->name; 1294 PetscFunctionReturn(0); 1295 } 1296 1297 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1298 { 1299 TS_RK *rk = (TS_RK*)ts->data; 1300 PetscErrorCode ierr; 1301 PetscBool match; 1302 RKTableauLink link; 1303 1304 PetscFunctionBegin; 1305 if (rk->tableau) { 1306 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 1307 if (match) PetscFunctionReturn(0); 1308 } 1309 for (link = RKTableauList; link; link=link->next) { 1310 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 1311 if (match) { 1312 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 1313 rk->tableau = &link->tab; 1314 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 1315 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1316 PetscFunctionReturn(0); 1317 } 1318 } 1319 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1320 PetscFunctionReturn(0); 1321 } 1322 1323 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1324 { 1325 TS_RK *rk = (TS_RK*)ts->data; 1326 1327 PetscFunctionBegin; 1328 if (ns) *ns = rk->tableau->s; 1329 if (Y) *Y = rk->Y; 1330 PetscFunctionReturn(0); 1331 } 1332 1333 static PetscErrorCode TSDestroy_RK(TS ts) 1334 { 1335 PetscErrorCode ierr; 1336 1337 PetscFunctionBegin; 1338 ierr = TSReset_RK(ts);CHKERRQ(ierr); 1339 if (ts->dm) { 1340 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1341 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1342 } 1343 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1344 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetOrder_C",NULL);CHKERRQ(ierr); 1345 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 1346 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 1347 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",NULL);CHKERRQ(ierr); 1348 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",NULL);CHKERRQ(ierr); 1349 PetscFunctionReturn(0); 1350 } 1351 1352 /* 1353 This defines the nonlinear equation that is to be solved with SNES 1354 We do not need to solve the equation; we just use SNES to approximate the Jacobian 1355 */ 1356 static PetscErrorCode SNESTSFormFunction_RK(SNES snes,Vec x,Vec y,TS ts) 1357 { 1358 TS_RK *rk = (TS_RK*)ts->data; 1359 DM dm,dmsave; 1360 PetscErrorCode ierr; 1361 1362 PetscFunctionBegin; 1363 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1364 /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */ 1365 dmsave = ts->dm; 1366 ts->dm = dm; 1367 ierr = TSComputeRHSFunction(ts,rk->stage_time,x,y);CHKERRQ(ierr); 1368 ts->dm = dmsave; 1369 PetscFunctionReturn(0); 1370 } 1371 1372 static PetscErrorCode SNESTSFormJacobian_RK(SNES snes,Vec x,Mat A,Mat B,TS ts) 1373 { 1374 TS_RK *rk = (TS_RK*)ts->data; 1375 DM dm,dmsave; 1376 PetscErrorCode ierr; 1377 1378 PetscFunctionBegin; 1379 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1380 dmsave = ts->dm; 1381 ts->dm = dm; 1382 ierr = TSComputeRHSJacobian(ts,rk->stage_time,x,A,B);CHKERRQ(ierr); 1383 ts->dm = dmsave; 1384 PetscFunctionReturn(0); 1385 } 1386 1387 /*@C 1388 TSRKSetMultirate - Use the interpolation-based multirate RK method 1389 1390 Logically collective 1391 1392 Input Parameter: 1393 + ts - timestepping context 1394 - 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 1395 1396 Options Database: 1397 . -ts_rk_multirate - <true,false> 1398 1399 Notes: 1400 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). 1401 1402 Level: intermediate 1403 1404 .seealso: TSRKGetMultirate() 1405 @*/ 1406 PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate) 1407 { 1408 PetscErrorCode ierr; 1409 1410 PetscFunctionBegin; 1411 ierr = PetscTryMethod(ts,"TSRKSetMultirate_C",(TS,PetscBool),(ts,use_multirate));CHKERRQ(ierr); 1412 PetscFunctionReturn(0); 1413 } 1414 1415 /*@C 1416 TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method 1417 1418 Not collective 1419 1420 Input Parameter: 1421 . ts - timestepping context 1422 1423 Output Parameter: 1424 . use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise 1425 1426 Level: intermediate 1427 1428 .seealso: TSRKSetMultirate() 1429 @*/ 1430 PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate) 1431 { 1432 PetscErrorCode ierr; 1433 1434 PetscFunctionBegin; 1435 ierr = PetscUseMethod(ts,"TSRKGetMultirate_C",(TS,PetscBool*),(ts,use_multirate));CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 /*MC 1440 TSRK - ODE and DAE solver using Runge-Kutta schemes 1441 1442 The user should provide the right hand side of the equation 1443 using TSSetRHSFunction(). 1444 1445 Notes: 1446 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1447 1448 Level: beginner 1449 1450 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 1451 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirate(), TSRKGetMultirate() 1452 1453 M*/ 1454 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1455 { 1456 TS_RK *rk; 1457 PetscErrorCode ierr; 1458 1459 PetscFunctionBegin; 1460 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1461 1462 ts->ops->reset = TSReset_RK; 1463 ts->ops->destroy = TSDestroy_RK; 1464 ts->ops->view = TSView_RK; 1465 ts->ops->load = TSLoad_RK; 1466 ts->ops->setup = TSSetUp_RK; 1467 ts->ops->interpolate = TSInterpolate_RK; 1468 ts->ops->step = TSStep_RK; 1469 ts->ops->evaluatestep = TSEvaluateStep_RK; 1470 ts->ops->rollback = TSRollBack_RK; 1471 ts->ops->setfromoptions = TSSetFromOptions_RK; 1472 ts->ops->getstages = TSGetStages_RK; 1473 1474 ts->ops->snesfunction = SNESTSFormFunction_RK; 1475 ts->ops->snesjacobian = SNESTSFormJacobian_RK; 1476 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1477 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1478 ts->ops->adjointstep = TSAdjointStep_RK; 1479 ts->ops->adjointreset = TSAdjointReset_RK; 1480 1481 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1482 ts->ops->forwardsetup = TSForwardSetUp_RK; 1483 ts->ops->forwardreset = TSForwardReset_RK; 1484 ts->ops->forwardstep = TSForwardStep_RK; 1485 ts->ops->forwardgetstages= TSForwardGetStages_RK; 1486 1487 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1488 ts->data = (void*)rk; 1489 1490 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetOrder_C",TSRKGetOrder_RK);CHKERRQ(ierr); 1491 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1492 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1493 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",TSRKSetMultirate_RK);CHKERRQ(ierr); 1494 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",TSRKGetMultirate_RK);CHKERRQ(ierr); 1495 1496 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1497 rk->dtratio = 1; 1498 PetscFunctionReturn(0); 1499 } 1500