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