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