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 Use SNES to compute the Jacobian so that finite differencing could be used when TS Jacobian is not available. 812 */ 813 static PetscErrorCode TSFormRHSJacobian_Private(Vec x,Mat J,Mat Jpre,TS ts) 814 { 815 SNES snes = ts->snes; 816 PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*) = NULL; 817 PetscErrorCode ierr; 818 819 PetscFunctionBegin; 820 /* 821 Unlike implicit methods, explicit methods do not have SNESMatFDColoring in the snes object 822 because SNESSolve() has not been called yet. Insteading of querying it, we check the 823 Jacobian compute function directly to determin if FD coloring is used. 824 */ 825 ierr = SNESGetJacobian(snes,NULL,NULL,&jac,NULL);CHKERRQ(ierr); 826 if (jac == SNESComputeJacobianDefaultColor) { 827 Vec f; 828 ierr = SNESSetSolution(snes,x);CHKERRQ(ierr); 829 ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr); 830 /* Force MatFDColoringApply to evaluate the SNES residual function for the base vector */ 831 ierr = SNESComputeFunction(snes,x,f);CHKERRQ(ierr); 832 } 833 ierr = SNESComputeJacobian(snes,x,J,Jpre);CHKERRQ(ierr); 834 PetscFunctionReturn(0); 835 } 836 837 /* 838 Assumptions: 839 - TSStep_RK() always evaluates the step with b, not bembed. 840 */ 841 static PetscErrorCode TSAdjointStep_RK(TS ts) 842 { 843 TS_RK *rk = (TS_RK*)ts->data; 844 TS quadts = ts->quadraturets; 845 RKTableau tab = rk->tableau; 846 Mat J,Jpre,Jquad; 847 const PetscInt s = tab->s; 848 const PetscReal *A = tab->A,*b = tab->b,*c = tab->c; 849 PetscScalar *w = rk->work,*xarr; 850 Vec *Y = rk->Y,*VecsDeltaLam = rk->VecsDeltaLam,VecDeltaMu = rk->VecDeltaMu,*VecsSensiTemp = rk->VecsSensiTemp; 851 Vec *VecsDeltaLam2 = rk->VecsDeltaLam2,VecDeltaMu2 = rk->VecDeltaMu2,*VecsSensi2Temp = rk->VecsSensi2Temp; 852 Vec VecDRDUTransCol = ts->vec_drdu_col,VecDRDPTransCol = ts->vec_drdp_col; 853 PetscInt i,j,nadj; 854 PetscReal t = ts->ptime; 855 PetscReal h = ts->time_step; 856 PetscErrorCode ierr; 857 858 PetscFunctionBegin; 859 rk->status = TS_STEP_INCOMPLETE; 860 861 ierr = TSGetRHSJacobian(ts,&J,&Jpre,NULL,NULL);CHKERRQ(ierr); 862 if (quadts) { 863 ierr = TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);CHKERRQ(ierr); 864 } 865 for (i=s-1; i>=0; i--) { 866 if (tab->FSAL && i == s-1) { 867 /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/ 868 continue; 869 } 870 rk->stage_time = t + h*(1.0-c[i]); 871 ierr = TSFormRHSJacobian_Private(Y[i],J,Jpre,ts);CHKERRQ(ierr); 872 if (quadts) { 873 ierr = TSComputeRHSJacobian(quadts,rk->stage_time,Y[i],Jquad,Jquad);CHKERRQ(ierr); /* get r_u^T */ 874 } 875 if (ts->vecs_sensip) { 876 ierr = TSComputeRHSJacobianP(ts,rk->stage_time,Y[i],ts->Jacprhs);CHKERRQ(ierr); /* get f_p */ 877 if (quadts) { 878 ierr = TSComputeRHSJacobianP(quadts,rk->stage_time,Y[i],quadts->Jacprhs);CHKERRQ(ierr); /* get f_p for the quadrature */ 879 } 880 } 881 882 if (b[i]) { 883 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]/b[i]; /* coefficients for computing VecsSensiTemp */ 884 } else { 885 for (j=i+1; j<s; j++) w[j-i-1] = A[j*s+i]; /* coefficients for computing VecsSensiTemp */ 886 } 887 888 for (nadj=0; nadj<ts->numcost; nadj++) { 889 /* Stage values of lambda */ 890 if (b[i]) { 891 /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */ 892 ierr = VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);CHKERRQ(ierr); /* VecDeltaLam is an vec array of size s by numcost */ 893 ierr = VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 894 ierr = MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]);CHKERRQ(ierr); /* VecsSensiTemp will be reused by 2nd-order adjoint */ 895 ierr = VecScale(VecsDeltaLam[nadj*s+i],-h*b[i]);CHKERRQ(ierr); 896 if (quadts) { 897 ierr = MatDenseGetColumn(Jquad,nadj,&xarr);CHKERRQ(ierr); 898 ierr = VecPlaceArray(VecDRDUTransCol,xarr);CHKERRQ(ierr); 899 ierr = VecAXPY(VecsDeltaLam[nadj*s+i],-h*b[i],VecDRDUTransCol);CHKERRQ(ierr); 900 ierr = VecResetArray(VecDRDUTransCol);CHKERRQ(ierr); 901 ierr = MatDenseRestoreColumn(Jquad,&xarr);CHKERRQ(ierr); 902 } 903 } else { 904 /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */ 905 ierr = VecSet(VecsSensiTemp[nadj],0);CHKERRQ(ierr); 906 ierr = VecMAXPY(VecsSensiTemp[nadj],s-i-1,w,&VecsDeltaLam[nadj*s+i+1]);CHKERRQ(ierr); 907 ierr = MatMultTranspose(J,VecsSensiTemp[nadj],VecsDeltaLam[nadj*s+i]);CHKERRQ(ierr); 908 ierr = VecScale(VecsDeltaLam[nadj*s+i],-h);CHKERRQ(ierr); 909 } 910 911 /* Stage values of mu */ 912 if (ts->vecs_sensip) { 913 ierr = MatMultTranspose(ts->Jacprhs,VecsSensiTemp[nadj],VecDeltaMu);CHKERRQ(ierr); 914 if (b[i]) { 915 ierr = VecScale(VecDeltaMu,-h*b[i]);CHKERRQ(ierr); 916 if (quadts) { 917 ierr = MatDenseGetColumn(quadts->Jacprhs,nadj,&xarr);CHKERRQ(ierr); 918 ierr = VecPlaceArray(VecDRDPTransCol,xarr);CHKERRQ(ierr); 919 ierr = VecAXPY(VecDeltaMu,-h*b[i],VecDRDPTransCol);CHKERRQ(ierr); 920 ierr = VecResetArray(VecDRDPTransCol);CHKERRQ(ierr); 921 ierr = MatDenseRestoreColumn(quadts->Jacprhs,&xarr);CHKERRQ(ierr); 922 } 923 } else { 924 ierr = VecScale(VecDeltaMu,-h);CHKERRQ(ierr); 925 } 926 ierr = VecAXPY(ts->vecs_sensip[nadj],1.,VecDeltaMu);CHKERRQ(ierr); /* update sensip for each stage */ 927 } 928 } 929 930 if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */ 931 /* Get w1 at t_{n+1} from TLM matrix */ 932 ierr = MatDenseGetColumn(rk->MatsFwdStageSensip[i],0,&xarr);CHKERRQ(ierr); 933 ierr = VecPlaceArray(ts->vec_sensip_col,xarr);CHKERRQ(ierr); 934 /* lambda_s^T F_UU w_1 */ 935 ierr = TSComputeRHSHessianProductFunctionUU(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_guu);CHKERRQ(ierr); 936 if (quadts) { 937 /* R_UU w_1 */ 938 ierr = TSComputeRHSHessianProductFunctionUU(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_guu);CHKERRQ(ierr); 939 } 940 if (ts->vecs_sensip) { 941 /* lambda_s^T F_UP w_2 */ 942 ierr = TSComputeRHSHessianProductFunctionUP(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gup);CHKERRQ(ierr); 943 if (quadts) { 944 /* R_UP w_2 */ 945 ierr = TSComputeRHSHessianProductFunctionUP(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gup);CHKERRQ(ierr); 946 } 947 } 948 if (ts->vecs_sensi2p) { 949 /* lambda_s^T F_PU w_1 */ 950 ierr = TSComputeRHSHessianProductFunctionPU(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_sensip_col,ts->vecs_gpu);CHKERRQ(ierr); 951 /* lambda_s^T F_PP w_2 */ 952 ierr = TSComputeRHSHessianProductFunctionPP(ts,rk->stage_time,Y[i],VecsSensiTemp,ts->vec_dir,ts->vecs_gpp);CHKERRQ(ierr); 953 if (b[i] && quadts) { 954 /* R_PU w_1 */ 955 ierr = TSComputeRHSHessianProductFunctionPU(quadts,rk->stage_time,Y[i],NULL,ts->vec_sensip_col,ts->vecs_gpu);CHKERRQ(ierr); 956 /* R_PP w_2 */ 957 ierr = TSComputeRHSHessianProductFunctionPP(quadts,rk->stage_time,Y[i],NULL,ts->vec_dir,ts->vecs_gpp);CHKERRQ(ierr); 958 } 959 } 960 ierr = VecResetArray(ts->vec_sensip_col);CHKERRQ(ierr); 961 ierr = MatDenseRestoreColumn(rk->MatsFwdStageSensip[i],&xarr);CHKERRQ(ierr); 962 963 for (nadj=0; nadj<ts->numcost; nadj++) { 964 /* Stage values of lambda */ 965 if (b[i]) { 966 /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */ 967 ierr = VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);CHKERRQ(ierr); 968 ierr = VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);CHKERRQ(ierr); 969 ierr = MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);CHKERRQ(ierr); 970 ierr = VecScale(VecsDeltaLam2[nadj*s+i],-h*b[i]);CHKERRQ(ierr); 971 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_guu[nadj]);CHKERRQ(ierr); 972 if (ts->vecs_sensip) { 973 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h*b[i],ts->vecs_gup[nadj]);CHKERRQ(ierr); 974 } 975 } else { 976 /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */ 977 ierr = VecSet(VecsDeltaLam2[nadj*s+i],0);CHKERRQ(ierr); 978 ierr = VecMAXPY(VecsSensi2Temp[nadj],s-i-1,w,&VecsDeltaLam2[nadj*s+i+1]);CHKERRQ(ierr); 979 ierr = MatMultTranspose(J,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj*s+i]);CHKERRQ(ierr); 980 ierr = VecScale(VecsDeltaLam2[nadj*s+i],-h);CHKERRQ(ierr); 981 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_guu[nadj]);CHKERRQ(ierr); 982 if (ts->vecs_sensip) { 983 ierr = VecAXPY(VecsDeltaLam2[nadj*s+i],-h,ts->vecs_gup[nadj]);CHKERRQ(ierr); 984 } 985 } 986 if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */ 987 ierr = MatMultTranspose(ts->Jacprhs,VecsSensi2Temp[nadj],VecDeltaMu2);CHKERRQ(ierr); 988 if (b[i]) { 989 ierr = VecScale(VecDeltaMu2,-h*b[i]);CHKERRQ(ierr); 990 ierr = VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpu[nadj]);CHKERRQ(ierr); 991 ierr = VecAXPY(VecDeltaMu2,-h*b[i],ts->vecs_gpp[nadj]);CHKERRQ(ierr); 992 } else { 993 ierr = VecScale(VecDeltaMu2,-h);CHKERRQ(ierr); 994 ierr = VecAXPY(VecDeltaMu2,-h,ts->vecs_gpu[nadj]);CHKERRQ(ierr); 995 ierr = VecAXPY(VecDeltaMu2,-h,ts->vecs_gpp[nadj]);CHKERRQ(ierr); 996 } 997 ierr = VecAXPY(ts->vecs_sensi2p[nadj],1,VecDeltaMu2);CHKERRQ(ierr); /* update sensi2p for each stage */ 998 } 999 } 1000 } 1001 } 1002 1003 for (j=0; j<s; j++) w[j] = 1.0; 1004 for (nadj=0; nadj<ts->numcost; nadj++) { /* no need to do this for mu's */ 1005 ierr = VecMAXPY(ts->vecs_sensi[nadj],s,w,&VecsDeltaLam[nadj*s]);CHKERRQ(ierr); 1006 if (ts->vecs_sensi2) { 1007 ierr = VecMAXPY(ts->vecs_sensi2[nadj],s,w,&VecsDeltaLam2[nadj*s]);CHKERRQ(ierr); 1008 } 1009 } 1010 rk->status = TS_STEP_COMPLETE; 1011 PetscFunctionReturn(0); 1012 } 1013 1014 static PetscErrorCode TSAdjointReset_RK(TS ts) 1015 { 1016 TS_RK *rk = (TS_RK*)ts->data; 1017 RKTableau tab = rk->tableau; 1018 PetscErrorCode ierr; 1019 1020 PetscFunctionBegin; 1021 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 1022 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 1023 ierr = VecDestroy(&rk->VecDeltaMu);CHKERRQ(ierr); 1024 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam2);CHKERRQ(ierr); 1025 ierr = VecDestroy(&rk->VecDeltaMu2);CHKERRQ(ierr); 1026 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensi2Temp);CHKERRQ(ierr); 1027 PetscFunctionReturn(0); 1028 } 1029 1030 static PetscErrorCode TSInterpolate_RK(TS ts,PetscReal itime,Vec X) 1031 { 1032 TS_RK *rk = (TS_RK*)ts->data; 1033 PetscInt s = rk->tableau->s,p = rk->tableau->p,i,j; 1034 PetscReal h; 1035 PetscReal tt,t; 1036 PetscScalar *b; 1037 const PetscReal *B = rk->tableau->binterp; 1038 PetscErrorCode ierr; 1039 1040 PetscFunctionBegin; 1041 if (!B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRK %s does not have an interpolation formula",rk->tableau->name); 1042 1043 switch (rk->status) { 1044 case TS_STEP_INCOMPLETE: 1045 case TS_STEP_PENDING: 1046 h = ts->time_step; 1047 t = (itime - ts->ptime)/h; 1048 break; 1049 case TS_STEP_COMPLETE: 1050 h = ts->ptime - ts->ptime_prev; 1051 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 1052 break; 1053 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus"); 1054 } 1055 ierr = PetscMalloc1(s,&b);CHKERRQ(ierr); 1056 for (i=0; i<s; i++) b[i] = 0; 1057 for (j=0,tt=t; j<p; j++,tt*=t) { 1058 for (i=0; i<s; i++) { 1059 b[i] += h * B[i*p+j] * tt; 1060 } 1061 } 1062 ierr = VecCopy(rk->Y[0],X);CHKERRQ(ierr); 1063 ierr = VecMAXPY(X,s,b,rk->YdotRHS);CHKERRQ(ierr); 1064 ierr = PetscFree(b);CHKERRQ(ierr); 1065 PetscFunctionReturn(0); 1066 } 1067 1068 /*------------------------------------------------------------*/ 1069 1070 static PetscErrorCode TSRKTableauReset(TS ts) 1071 { 1072 TS_RK *rk = (TS_RK*)ts->data; 1073 RKTableau tab = rk->tableau; 1074 PetscErrorCode ierr; 1075 1076 PetscFunctionBegin; 1077 if (!tab) PetscFunctionReturn(0); 1078 ierr = PetscFree(rk->work);CHKERRQ(ierr); 1079 ierr = VecDestroyVecs(tab->s,&rk->Y);CHKERRQ(ierr); 1080 ierr = VecDestroyVecs(tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1081 ierr = VecDestroyVecs(tab->s*ts->numcost,&rk->VecsDeltaLam);CHKERRQ(ierr); 1082 ierr = VecDestroyVecs(ts->numcost,&rk->VecsSensiTemp);CHKERRQ(ierr); 1083 ierr = VecDestroy(&rk->VecDeltaMu);CHKERRQ(ierr); 1084 PetscFunctionReturn(0); 1085 } 1086 1087 static PetscErrorCode TSReset_RK(TS ts) 1088 { 1089 PetscErrorCode ierr; 1090 1091 PetscFunctionBegin; 1092 ierr = TSRKTableauReset(ts);CHKERRQ(ierr); 1093 if (ts->use_splitrhsfunction) { 1094 ierr = PetscTryMethod(ts,"TSReset_RK_MultirateSplit_C",(TS),(ts));CHKERRQ(ierr); 1095 } else { 1096 ierr = PetscTryMethod(ts,"TSReset_RK_MultirateNonsplit_C",(TS),(ts));CHKERRQ(ierr); 1097 } 1098 PetscFunctionReturn(0); 1099 } 1100 1101 static PetscErrorCode DMCoarsenHook_TSRK(DM fine,DM coarse,void *ctx) 1102 { 1103 PetscFunctionBegin; 1104 PetscFunctionReturn(0); 1105 } 1106 1107 static PetscErrorCode DMRestrictHook_TSRK(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 1108 { 1109 PetscFunctionBegin; 1110 PetscFunctionReturn(0); 1111 } 1112 1113 1114 static PetscErrorCode DMSubDomainHook_TSRK(DM dm,DM subdm,void *ctx) 1115 { 1116 PetscFunctionBegin; 1117 PetscFunctionReturn(0); 1118 } 1119 1120 static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx) 1121 { 1122 1123 PetscFunctionBegin; 1124 PetscFunctionReturn(0); 1125 } 1126 1127 static PetscErrorCode TSRKTableauSetUp(TS ts) 1128 { 1129 TS_RK *rk = (TS_RK*)ts->data; 1130 RKTableau tab = rk->tableau; 1131 PetscErrorCode ierr; 1132 1133 PetscFunctionBegin; 1134 ierr = PetscMalloc1(tab->s,&rk->work);CHKERRQ(ierr); 1135 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->Y);CHKERRQ(ierr); 1136 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&rk->YdotRHS);CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 static PetscErrorCode TSSetUp_RK(TS ts) 1141 { 1142 TS quadts = ts->quadraturets; 1143 PetscErrorCode ierr; 1144 DM dm; 1145 1146 PetscFunctionBegin; 1147 ierr = TSCheckImplicitTerm(ts);CHKERRQ(ierr); 1148 ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr); 1149 if (quadts && ts->costintegralfwd) { 1150 Mat Jquad; 1151 ierr = TSGetRHSJacobian(quadts,&Jquad,NULL,NULL,NULL);CHKERRQ(ierr); 1152 } 1153 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1154 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1155 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1156 if (ts->use_splitrhsfunction) { 1157 ierr = PetscTryMethod(ts,"TSSetUp_RK_MultirateSplit_C",(TS),(ts));CHKERRQ(ierr); 1158 } else { 1159 ierr = PetscTryMethod(ts,"TSSetUp_RK_MultirateNonsplit_C",(TS),(ts));CHKERRQ(ierr); 1160 } 1161 PetscFunctionReturn(0); 1162 } 1163 1164 static PetscErrorCode TSSetFromOptions_RK(PetscOptionItems *PetscOptionsObject,TS ts) 1165 { 1166 TS_RK *rk = (TS_RK*)ts->data; 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 ierr = PetscOptionsHead(PetscOptionsObject,"RK ODE solver options");CHKERRQ(ierr); 1171 { 1172 RKTableauLink link; 1173 PetscInt count,choice; 1174 PetscBool flg,use_multirate = PETSC_FALSE; 1175 const char **namelist; 1176 1177 for (link=RKTableauList,count=0; link; link=link->next,count++) ; 1178 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr); 1179 for (link=RKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 1180 ierr = PetscOptionsBool("-ts_rk_multirate","Use interpolation-based multirate RK method","TSRKSetMultirate",rk->use_multirate,&use_multirate,&flg);CHKERRQ(ierr); 1181 if (flg) { 1182 ierr = TSRKSetMultirate(ts,use_multirate);CHKERRQ(ierr); 1183 } 1184 ierr = PetscOptionsEList("-ts_rk_type","Family of RK method","TSRKSetType",(const char*const*)namelist,count,rk->tableau->name,&choice,&flg);CHKERRQ(ierr); 1185 if (flg) {ierr = TSRKSetType(ts,namelist[choice]);CHKERRQ(ierr);} 1186 ierr = PetscFree(namelist);CHKERRQ(ierr); 1187 } 1188 ierr = PetscOptionsTail();CHKERRQ(ierr); 1189 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Multirate methods options","");CHKERRQ(ierr); 1190 ierr = PetscOptionsInt("-ts_rk_dtratio","time step ratio between slow and fast","",rk->dtratio,&rk->dtratio,NULL);CHKERRQ(ierr); 1191 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1192 PetscFunctionReturn(0); 1193 } 1194 1195 static PetscErrorCode TSView_RK(TS ts,PetscViewer viewer) 1196 { 1197 TS_RK *rk = (TS_RK*)ts->data; 1198 PetscBool iascii; 1199 PetscErrorCode ierr; 1200 1201 PetscFunctionBegin; 1202 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1203 if (iascii) { 1204 RKTableau tab = rk->tableau; 1205 TSRKType rktype; 1206 char buf[512]; 1207 ierr = TSRKGetType(ts,&rktype);CHKERRQ(ierr); 1208 ierr = PetscViewerASCIIPrintf(viewer," RK type %s\n",rktype);CHKERRQ(ierr); 1209 ierr = PetscViewerASCIIPrintf(viewer," Order: %D\n",tab->order);CHKERRQ(ierr); 1210 ierr = PetscViewerASCIIPrintf(viewer," FSAL property: %s\n",tab->FSAL ? "yes" : "no");CHKERRQ(ierr); 1211 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 1212 ierr = PetscViewerASCIIPrintf(viewer," Abscissa c = %s\n",buf);CHKERRQ(ierr); 1213 } 1214 PetscFunctionReturn(0); 1215 } 1216 1217 static PetscErrorCode TSLoad_RK(TS ts,PetscViewer viewer) 1218 { 1219 PetscErrorCode ierr; 1220 TSAdapt adapt; 1221 1222 PetscFunctionBegin; 1223 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 1224 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr); 1225 PetscFunctionReturn(0); 1226 } 1227 1228 /*@C 1229 TSRKSetType - Set the type of RK scheme 1230 1231 Logically collective 1232 1233 Input Parameter: 1234 + ts - timestepping context 1235 - rktype - type of RK-scheme 1236 1237 Options Database: 1238 . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs> 1239 1240 Level: intermediate 1241 1242 .seealso: TSRKGetType(), TSRK, TSRKType, TSRK1FE, TSRK2A, TSRK3, TSRK3BS, TSRK4, TSRK5F, TSRK5DP, TSRK5BS, TSRK6VR, TSRK7VR, TSRK8VR 1243 @*/ 1244 PetscErrorCode TSRKSetType(TS ts,TSRKType rktype) 1245 { 1246 PetscErrorCode ierr; 1247 1248 PetscFunctionBegin; 1249 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1250 PetscValidCharPointer(rktype,2); 1251 ierr = PetscTryMethod(ts,"TSRKSetType_C",(TS,TSRKType),(ts,rktype));CHKERRQ(ierr); 1252 PetscFunctionReturn(0); 1253 } 1254 1255 /*@C 1256 TSRKGetType - Get the type of RK scheme 1257 1258 Logically collective 1259 1260 Input Parameter: 1261 . ts - timestepping context 1262 1263 Output Parameter: 1264 . rktype - type of RK-scheme 1265 1266 Level: intermediate 1267 1268 .seealso: TSRKGetType() 1269 @*/ 1270 PetscErrorCode TSRKGetType(TS ts,TSRKType *rktype) 1271 { 1272 PetscErrorCode ierr; 1273 1274 PetscFunctionBegin; 1275 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1276 ierr = PetscUseMethod(ts,"TSRKGetType_C",(TS,TSRKType*),(ts,rktype));CHKERRQ(ierr); 1277 PetscFunctionReturn(0); 1278 } 1279 1280 static PetscErrorCode TSRKGetType_RK(TS ts,TSRKType *rktype) 1281 { 1282 TS_RK *rk = (TS_RK*)ts->data; 1283 1284 PetscFunctionBegin; 1285 *rktype = rk->tableau->name; 1286 PetscFunctionReturn(0); 1287 } 1288 1289 static PetscErrorCode TSRKSetType_RK(TS ts,TSRKType rktype) 1290 { 1291 TS_RK *rk = (TS_RK*)ts->data; 1292 PetscErrorCode ierr; 1293 PetscBool match; 1294 RKTableauLink link; 1295 1296 PetscFunctionBegin; 1297 if (rk->tableau) { 1298 ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr); 1299 if (match) PetscFunctionReturn(0); 1300 } 1301 for (link = RKTableauList; link; link=link->next) { 1302 ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr); 1303 if (match) { 1304 if (ts->setupcalled) {ierr = TSRKTableauReset(ts);CHKERRQ(ierr);} 1305 rk->tableau = &link->tab; 1306 if (ts->setupcalled) {ierr = TSRKTableauSetUp(ts);CHKERRQ(ierr);} 1307 ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE; 1308 PetscFunctionReturn(0); 1309 } 1310 } 1311 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype); 1312 PetscFunctionReturn(0); 1313 } 1314 1315 static PetscErrorCode TSGetStages_RK(TS ts,PetscInt *ns,Vec **Y) 1316 { 1317 TS_RK *rk = (TS_RK*)ts->data; 1318 1319 PetscFunctionBegin; 1320 if (ns) *ns = rk->tableau->s; 1321 if (Y) *Y = rk->Y; 1322 PetscFunctionReturn(0); 1323 } 1324 1325 static PetscErrorCode TSDestroy_RK(TS ts) 1326 { 1327 PetscErrorCode ierr; 1328 1329 PetscFunctionBegin; 1330 ierr = TSReset_RK(ts);CHKERRQ(ierr); 1331 if (ts->dm) { 1332 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSRK,DMRestrictHook_TSRK,ts);CHKERRQ(ierr); 1333 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSRK,DMSubDomainRestrictHook_TSRK,ts);CHKERRQ(ierr); 1334 } 1335 ierr = PetscFree(ts->data);CHKERRQ(ierr); 1336 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",NULL);CHKERRQ(ierr); 1337 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",NULL);CHKERRQ(ierr); 1338 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",NULL);CHKERRQ(ierr); 1339 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",NULL);CHKERRQ(ierr); 1340 PetscFunctionReturn(0); 1341 } 1342 1343 /* 1344 This defines the nonlinear equation that is to be solved with SNES 1345 We do not need to solve the equation; we just use SNES to approximate the Jacobian 1346 */ 1347 static PetscErrorCode SNESTSFormFunction_RK(SNES snes,Vec x,Vec y,TS ts) 1348 { 1349 TS_RK *rk = (TS_RK*)ts->data; 1350 DM dm,dmsave; 1351 PetscErrorCode ierr; 1352 1353 PetscFunctionBegin; 1354 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1355 /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */ 1356 dmsave = ts->dm; 1357 ts->dm = dm; 1358 ierr = TSComputeRHSFunction(ts,rk->stage_time,x,y);CHKERRQ(ierr); 1359 ts->dm = dmsave; 1360 PetscFunctionReturn(0); 1361 } 1362 1363 static PetscErrorCode SNESTSFormJacobian_RK(SNES snes,Vec x,Mat A,Mat B,TS ts) 1364 { 1365 TS_RK *rk = (TS_RK*)ts->data; 1366 DM dm,dmsave; 1367 PetscErrorCode ierr; 1368 1369 PetscFunctionBegin; 1370 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1371 dmsave = ts->dm; 1372 ts->dm = dm; 1373 ierr = TSComputeRHSJacobian(ts,rk->stage_time,x,A,B);CHKERRQ(ierr); 1374 ts->dm = dmsave; 1375 PetscFunctionReturn(0); 1376 } 1377 1378 /*@C 1379 TSRKSetMultirate - Use the interpolation-based multirate RK method 1380 1381 Logically collective 1382 1383 Input Parameter: 1384 + ts - timestepping context 1385 - 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 1386 1387 Options Database: 1388 . -ts_rk_multirate - <true,false> 1389 1390 Notes: 1391 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). 1392 1393 Level: intermediate 1394 1395 .seealso: TSRKGetMultirate() 1396 @*/ 1397 PetscErrorCode TSRKSetMultirate(TS ts,PetscBool use_multirate) 1398 { 1399 PetscErrorCode ierr; 1400 1401 PetscFunctionBegin; 1402 ierr = PetscTryMethod(ts,"TSRKSetMultirate_C",(TS,PetscBool),(ts,use_multirate));CHKERRQ(ierr); 1403 PetscFunctionReturn(0); 1404 } 1405 1406 /*@C 1407 TSRKGetMultirate - Gets whether to Use the interpolation-based multirate RK method 1408 1409 Not collective 1410 1411 Input Parameter: 1412 . ts - timestepping context 1413 1414 Output Parameter: 1415 . use_multirate - PETSC_TRUE if the multirate RK method is enabled, PETSC_FALSE otherwise 1416 1417 Level: intermediate 1418 1419 .seealso: TSRKSetMultirate() 1420 @*/ 1421 PetscErrorCode TSRKGetMultirate(TS ts,PetscBool *use_multirate) 1422 { 1423 PetscErrorCode ierr; 1424 1425 PetscFunctionBegin; 1426 ierr = PetscUseMethod(ts,"TSRKGetMultirate_C",(TS,PetscBool*),(ts,use_multirate));CHKERRQ(ierr); 1427 PetscFunctionReturn(0); 1428 } 1429 1430 /*MC 1431 TSRK - ODE and DAE solver using Runge-Kutta schemes 1432 1433 The user should provide the right hand side of the equation 1434 using TSSetRHSFunction(). 1435 1436 Notes: 1437 The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type 1438 1439 Level: beginner 1440 1441 .seealso: TSCreate(), TS, TSSetType(), TSRKSetType(), TSRKGetType(), TSRKSetFullyImplicit(), TSRK2D, TTSRK2E, TSRK3, 1442 TSRK4, TSRK5, TSRKPRSSP2, TSRKBPR3, TSRKType, TSRKRegister(), TSRKSetMultirate(), TSRKGetMultirate() 1443 1444 M*/ 1445 PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts) 1446 { 1447 TS_RK *rk; 1448 PetscErrorCode ierr; 1449 1450 PetscFunctionBegin; 1451 ierr = TSRKInitializePackage();CHKERRQ(ierr); 1452 1453 ts->ops->reset = TSReset_RK; 1454 ts->ops->destroy = TSDestroy_RK; 1455 ts->ops->view = TSView_RK; 1456 ts->ops->load = TSLoad_RK; 1457 ts->ops->setup = TSSetUp_RK; 1458 ts->ops->interpolate = TSInterpolate_RK; 1459 ts->ops->step = TSStep_RK; 1460 ts->ops->evaluatestep = TSEvaluateStep_RK; 1461 ts->ops->rollback = TSRollBack_RK; 1462 ts->ops->setfromoptions = TSSetFromOptions_RK; 1463 ts->ops->getstages = TSGetStages_RK; 1464 1465 ts->ops->snesfunction = SNESTSFormFunction_RK; 1466 ts->ops->snesjacobian = SNESTSFormJacobian_RK; 1467 ts->ops->adjointintegral = TSAdjointCostIntegral_RK; 1468 ts->ops->adjointsetup = TSAdjointSetUp_RK; 1469 ts->ops->adjointstep = TSAdjointStep_RK; 1470 ts->ops->adjointreset = TSAdjointReset_RK; 1471 1472 ts->ops->forwardintegral = TSForwardCostIntegral_RK; 1473 ts->ops->forwardsetup = TSForwardSetUp_RK; 1474 ts->ops->forwardreset = TSForwardReset_RK; 1475 ts->ops->forwardstep = TSForwardStep_RK; 1476 ts->ops->forwardgetstages= TSForwardGetStages_RK; 1477 1478 ierr = PetscNewLog(ts,&rk);CHKERRQ(ierr); 1479 ts->data = (void*)rk; 1480 1481 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetType_C",TSRKGetType_RK);CHKERRQ(ierr); 1482 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetType_C",TSRKSetType_RK);CHKERRQ(ierr); 1483 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKSetMultirate_C",TSRKSetMultirate_RK);CHKERRQ(ierr); 1484 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSRKGetMultirate_C",TSRKGetMultirate_RK);CHKERRQ(ierr); 1485 1486 ierr = TSRKSetType(ts,TSRKDefault);CHKERRQ(ierr); 1487 rk->dtratio = 1; 1488 PetscFunctionReturn(0); 1489 } 1490