Actual source code: rk.c
1: /*
2: Code for time stepping with the Runge-Kutta method
4: Notes:
5: The general system is written as
7: Udot = F(t,U)
9: */
11: #include <petsc/private/tsimpl.h>
12: #include <petscdm.h>
13: #include <../src/ts/impls/explicit/rk/rk.h>
14: #include <../src/ts/impls/explicit/rk/mrk.h>
16: static TSRKType TSRKDefault = TSRK3BS;
17: static PetscBool TSRKRegisterAllCalled;
18: static PetscBool TSRKPackageInitialized;
20: static RKTableauLink RKTableauList;
22: /*MC
23: TSRK1FE - First order forward Euler scheme.
25: This method has one stage.
27: Options Database Key:
28: . -ts_rk_type 1fe - use type 1fe
30: Level: advanced
32: .seealso: [](ch_ts), `TSRK`, `TSRKType`, `TSRKSetType()`
33: M*/
34: /*MC
35: TSRK2A - Second order RK scheme (Heun's method).
37: This method has two stages.
39: Options Database Key:
40: . -ts_rk_type 2a - use type 2a
42: Level: advanced
44: .seealso: [](ch_ts), `TSRK`, `TSRKType`, `TSRKSetType()`
45: M*/
46: /*MC
47: TSRK2B - Second order RK scheme (the midpoint method).
49: This method has two stages.
51: Options Database Key:
52: . -ts_rk_type 2b - use type 2b
54: Level: advanced
56: .seealso: [](ch_ts), `TSRK`, `TSRKType`, `TSRKSetType()`
57: M*/
58: /*MC
59: TSRK3 - Third order RK scheme.
61: This method has three stages.
63: Options Database Key:
64: . -ts_rk_type 3 - use type 3
66: Level: advanced
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>
73: This method has four stages with the First Same As Last (FSAL) property.
75: Options Database Key:
76: . -ts_rk_type 3bs - use type 3bs
78: Level: advanced
80: .seealso: [](ch_ts), `TSRK`, `TSRKType`, `TSRKSetType()`
81: M*/
82: /*MC
83: TSRK4 - Fourth order RK scheme.
85: This is the classical Runge-Kutta method with four stages.
87: Options Database Key:
88: . -ts_rk_type 4 - use type 4
90: Level: advanced
92: .seealso: [](ch_ts), `TSRK`, `TSRKType`, `TSRKSetType()`
93: M*/
94: /*MC
95: TSRK5F - Fifth order Fehlberg RK scheme with a 4th order embedded method.
97: This method has six stages.
99: Options Database Key:
100: . -ts_rk_type 5f - use type 5f
102: Level: advanced
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>
109: This method has seven stages with the First Same As Last (FSAL) property.
111: Options Database Key:
112: . -ts_rk_type 5dp - use type 5dp
114: Level: advanced
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>
121: This method has eight stages with the First Same As Last (FSAL) property.
123: Options Database Key:
124: . -ts_rk_type 5bs - use type 5bs
126: Level: advanced
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>
134: This method has nine stages with the First Same As Last (FSAL) property.
136: Options Database Key:
137: . -ts_rk_type 6vr - use type 6vr
139: Level: advanced
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>
147: This method has ten stages.
149: Options Database Key:
150: . -ts_rk_type 7vr - use type 7vr
152: Level: advanced
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>
160: This method has thirteen stages.
162: Options Database Key:
163: . -ts_rk_type 8vr - use type 8vr
165: Level: advanced
167: .seealso: [](ch_ts), `TSRK`, `TSRKType`, `TSRKSetType()`
168: M*/
170: /*@C
171: TSRKRegisterAll - Registers all of the Runge-Kutta explicit methods in `TSRK`
173: Not Collective, but should be called by all processes which will need the schemes to be registered
175: Level: advanced
177: .seealso: [](ch_ts), `TSRKRegisterDestroy()`, `TSRKRegister()`
178: @*/
179: PetscErrorCode TSRKRegisterAll(void)
180: {
181: PetscFunctionBegin;
182: if (TSRKRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
183: TSRKRegisterAllCalled = PETSC_TRUE;
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: }
348: /*@C
349: TSRKRegisterDestroy - Frees the list of schemes that were registered by `TSRKRegister()`.
351: Not Collective
353: Level: advanced
355: .seealso: [](ch_ts), `TSRK`, `TSRKRegister()`, `TSRKRegisterAll()`
356: @*/
357: PetscErrorCode TSRKRegisterDestroy(void)
358: {
359: RKTableauLink link;
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: }
375: /*@C
376: TSRKInitializePackage - This function initializes everything in the `TSRK` package. It is called
377: from `TSInitializePackage()`.
379: Level: developer
381: .seealso: [](ch_ts), `TSInitializePackage()`, `PetscInitialize()`, `TSRKFinalizePackage()`
382: @*/
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: }
393: /*@C
394: TSRKFinalizePackage - This function destroys everything in the `TSRK` package. It is
395: called from `PetscFinalize()`.
397: Level: developer
399: .seealso: [](ch_ts), `PetscFinalize()`, `TSRKInitializePackage()`
400: @*/
401: PetscErrorCode TSRKFinalizePackage(void)
402: {
403: PetscFunctionBegin;
404: TSRKPackageInitialized = PETSC_FALSE;
405: PetscCall(TSRKRegisterDestroy());
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: /*@C
410: TSRKRegister - register an `TSRK` scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
412: Not Collective, but the same schemes should be registered on all processes on which they will be used, No Fortran Support
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)
425: Level: advanced
427: Note:
428: Several `TSRK` methods are provided, this function is only needed to create new methods.
430: .seealso: [](ch_ts), `TSRK`
431: @*/
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;
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);
447: PetscCall(TSRKInitializePackage());
448: PetscCall(PetscNew(&link));
449: t = &link->tab;
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;
467: if (bembed) {
468: PetscCall(PetscMalloc1(s, &t->bembed));
469: PetscCall(PetscArraycpy(t->bembed, bembed, s));
470: }
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));
480: link->next = RKTableauList;
481: RKTableauList = link;
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
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;
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: }
502: /*@
503: TSRKGetTableau - Get info on the `TSRK` tableau
505: Not Collective
507: Input Parameter:
508: . ts - timestepping context
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
520: Level: developer
522: Fortran Note:
523: Call `TSRKRestoreTableau()` when you no longer need access to the tableau values.
525: .seealso: [](ch_ts), `TSRK`, `TSRKRegister()`, `TSRKSetType()`
526: @*/
527: 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)
528: {
529: PetscFunctionBegin;
531: 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));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*
536: This is for single-step RK method
537: The step completion formula is
539: x1 = x0 + h b^T YdotRHS
541: This function can be called before or after ts->vec_sol has been updated.
542: Suppose we have a completion formula (b) and an embedded formula (be) of different order.
543: We can write
545: x1e = x0 + h be^T YdotRHS
546: = x1 - h b^T YdotRHS + h be^T YdotRHS
547: = x1 + h (be - b)^T YdotRHS
549: so we can evaluate the method with different order even after the step has been optimistically completed.
550: */
551: static PetscErrorCode TSEvaluateStep_RK(TS ts, PetscInt order, Vec X, PetscBool *done)
552: {
553: TS_RK *rk = (TS_RK *)ts->data;
554: RKTableau tab = rk->tableau;
555: PetscScalar *w = rk->work;
556: PetscReal h;
557: PetscInt s = tab->s, j;
559: PetscFunctionBegin;
560: switch (rk->status) {
561: case TS_STEP_INCOMPLETE:
562: case TS_STEP_PENDING:
563: h = ts->time_step;
564: break;
565: case TS_STEP_COMPLETE:
566: h = ts->ptime - ts->ptime_prev;
567: break;
568: default:
569: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
570: }
571: if (order == tab->order) {
572: if (rk->status == TS_STEP_INCOMPLETE) {
573: PetscCall(VecCopy(ts->vec_sol, X));
574: for (j = 0; j < s; j++) w[j] = h * tab->b[j] / rk->dtratio;
575: PetscCall(VecMAXPY(X, s, w, rk->YdotRHS));
576: } else PetscCall(VecCopy(ts->vec_sol, X));
577: PetscFunctionReturn(PETSC_SUCCESS);
578: } else if (order == tab->order - 1) {
579: if (!tab->bembed) goto unavailable;
580: if (rk->status == TS_STEP_INCOMPLETE) { /*Complete with the embedded method (be)*/
581: PetscCall(VecCopy(ts->vec_sol, X));
582: for (j = 0; j < s; j++) w[j] = h * tab->bembed[j];
583: PetscCall(VecMAXPY(X, s, w, rk->YdotRHS));
584: } else { /*Rollback and re-complete using (be-b) */
585: PetscCall(VecCopy(ts->vec_sol, X));
586: for (j = 0; j < s; j++) w[j] = h * (tab->bembed[j] - tab->b[j]);
587: PetscCall(VecMAXPY(X, s, w, rk->YdotRHS));
588: }
589: if (done) *done = PETSC_TRUE;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
592: unavailable:
593: 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,
594: tab->order, order);
595: *done = PETSC_FALSE;
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: static PetscErrorCode TSForwardCostIntegral_RK(TS ts)
600: {
601: TS_RK *rk = (TS_RK *)ts->data;
602: TS quadts = ts->quadraturets;
603: RKTableau tab = rk->tableau;
604: const PetscInt s = tab->s;
605: const PetscReal *b = tab->b, *c = tab->c;
606: Vec *Y = rk->Y;
607: PetscInt i;
609: PetscFunctionBegin;
610: /* No need to backup quadts->vec_sol since it can be reverted in TSRollBack_RK */
611: for (i = s - 1; i >= 0; i--) {
612: /* Evolve quadrature TS solution to compute integrals */
613: PetscCall(TSComputeRHSFunction(quadts, rk->ptime + rk->time_step * c[i], Y[i], ts->vec_costintegrand));
614: PetscCall(VecAXPY(quadts->vec_sol, rk->time_step * b[i], ts->vec_costintegrand));
615: }
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: static PetscErrorCode TSAdjointCostIntegral_RK(TS ts)
620: {
621: TS_RK *rk = (TS_RK *)ts->data;
622: RKTableau tab = rk->tableau;
623: TS quadts = ts->quadraturets;
624: const PetscInt s = tab->s;
625: const PetscReal *b = tab->b, *c = tab->c;
626: Vec *Y = rk->Y;
627: PetscInt i;
629: PetscFunctionBegin;
630: for (i = s - 1; i >= 0; i--) {
631: /* Evolve quadrature TS solution to compute integrals */
632: PetscCall(TSComputeRHSFunction(quadts, ts->ptime + ts->time_step * (1.0 - c[i]), Y[i], ts->vec_costintegrand));
633: PetscCall(VecAXPY(quadts->vec_sol, -ts->time_step * b[i], ts->vec_costintegrand));
634: }
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: static PetscErrorCode TSRollBack_RK(TS ts)
639: {
640: TS_RK *rk = (TS_RK *)ts->data;
641: TS quadts = ts->quadraturets;
642: RKTableau tab = rk->tableau;
643: const PetscInt s = tab->s;
644: const PetscReal *b = tab->b, *c = tab->c;
645: PetscScalar *w = rk->work;
646: Vec *Y = rk->Y, *YdotRHS = rk->YdotRHS;
647: PetscInt j;
648: PetscReal h;
650: PetscFunctionBegin;
651: switch (rk->status) {
652: case TS_STEP_INCOMPLETE:
653: case TS_STEP_PENDING:
654: h = ts->time_step;
655: break;
656: case TS_STEP_COMPLETE:
657: h = ts->ptime - ts->ptime_prev;
658: break;
659: default:
660: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
661: }
662: for (j = 0; j < s; j++) w[j] = -h * b[j];
663: PetscCall(VecMAXPY(ts->vec_sol, s, w, YdotRHS));
664: if (quadts && ts->costintegralfwd) {
665: for (j = 0; j < s; j++) {
666: /* Revert the quadrature TS solution */
667: PetscCall(TSComputeRHSFunction(quadts, rk->ptime + h * c[j], Y[j], ts->vec_costintegrand));
668: PetscCall(VecAXPY(quadts->vec_sol, -h * b[j], ts->vec_costintegrand));
669: }
670: }
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: static PetscErrorCode TSForwardStep_RK(TS ts)
675: {
676: TS_RK *rk = (TS_RK *)ts->data;
677: RKTableau tab = rk->tableau;
678: Mat J, *MatsFwdSensipTemp = rk->MatsFwdSensipTemp;
679: const PetscInt s = tab->s;
680: const PetscReal *A = tab->A, *c = tab->c, *b = tab->b;
681: Vec *Y = rk->Y;
682: PetscInt i, j;
683: PetscReal stage_time, h = ts->time_step;
684: PetscBool zero;
686: PetscFunctionBegin;
687: PetscCall(MatCopy(ts->mat_sensip, rk->MatFwdSensip0, SAME_NONZERO_PATTERN));
688: PetscCall(TSGetRHSJacobian(ts, &J, NULL, NULL, NULL));
690: for (i = 0; i < s; i++) {
691: stage_time = ts->ptime + h * c[i];
692: zero = PETSC_FALSE;
693: if (b[i] == 0 && i == s - 1) zero = PETSC_TRUE;
694: /* TLM Stage values */
695: if (!i) {
696: PetscCall(MatCopy(ts->mat_sensip, rk->MatsFwdStageSensip[i], SAME_NONZERO_PATTERN));
697: } else if (!zero) {
698: PetscCall(MatZeroEntries(rk->MatsFwdStageSensip[i]));
699: for (j = 0; j < i; j++) PetscCall(MatAXPY(rk->MatsFwdStageSensip[i], h * A[i * s + j], MatsFwdSensipTemp[j], SAME_NONZERO_PATTERN));
700: PetscCall(MatAXPY(rk->MatsFwdStageSensip[i], 1., ts->mat_sensip, SAME_NONZERO_PATTERN));
701: } else {
702: PetscCall(MatZeroEntries(rk->MatsFwdStageSensip[i]));
703: }
705: PetscCall(TSComputeRHSJacobian(ts, stage_time, Y[i], J, J));
706: PetscCall(MatMatMult(J, rk->MatsFwdStageSensip[i], MAT_REUSE_MATRIX, PETSC_DETERMINE, &MatsFwdSensipTemp[i]));
707: if (ts->Jacprhs) {
708: PetscCall(TSComputeRHSJacobianP(ts, stage_time, Y[i], ts->Jacprhs)); /* get f_p */
709: if (ts->vecs_sensi2p) { /* TLM used for 2nd-order adjoint */
710: PetscScalar *xarr;
711: PetscCall(MatDenseGetColumn(MatsFwdSensipTemp[i], 0, &xarr));
712: PetscCall(VecPlaceArray(rk->VecDeltaFwdSensipCol, xarr));
713: PetscCall(MatMultAdd(ts->Jacprhs, ts->vec_dir, rk->VecDeltaFwdSensipCol, rk->VecDeltaFwdSensipCol));
714: PetscCall(VecResetArray(rk->VecDeltaFwdSensipCol));
715: PetscCall(MatDenseRestoreColumn(MatsFwdSensipTemp[i], &xarr));
716: } else {
717: PetscCall(MatAXPY(MatsFwdSensipTemp[i], 1., ts->Jacprhs, SUBSET_NONZERO_PATTERN));
718: }
719: }
720: }
722: for (i = 0; i < s; i++) PetscCall(MatAXPY(ts->mat_sensip, h * b[i], rk->MatsFwdSensipTemp[i], SAME_NONZERO_PATTERN));
723: rk->status = TS_STEP_COMPLETE;
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: static PetscErrorCode TSForwardGetStages_RK(TS ts, PetscInt *ns, Mat **stagesensip)
728: {
729: TS_RK *rk = (TS_RK *)ts->data;
730: RKTableau tab = rk->tableau;
732: PetscFunctionBegin;
733: if (ns) *ns = tab->s;
734: if (stagesensip) *stagesensip = rk->MatsFwdStageSensip;
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: static PetscErrorCode TSForwardSetUp_RK(TS ts)
739: {
740: TS_RK *rk = (TS_RK *)ts->data;
741: RKTableau tab = rk->tableau;
742: PetscInt i;
744: PetscFunctionBegin;
745: /* backup sensitivity results for roll-backs */
746: PetscCall(MatDuplicate(ts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &rk->MatFwdSensip0));
748: PetscCall(PetscMalloc1(tab->s, &rk->MatsFwdStageSensip));
749: PetscCall(PetscMalloc1(tab->s, &rk->MatsFwdSensipTemp));
750: for (i = 0; i < tab->s; i++) {
751: PetscCall(MatDuplicate(ts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &rk->MatsFwdStageSensip[i]));
752: PetscCall(MatDuplicate(ts->mat_sensip, MAT_DO_NOT_COPY_VALUES, &rk->MatsFwdSensipTemp[i]));
753: }
754: PetscCall(VecDuplicate(ts->vec_sol, &rk->VecDeltaFwdSensipCol));
755: PetscFunctionReturn(PETSC_SUCCESS);
756: }
758: static PetscErrorCode TSForwardReset_RK(TS ts)
759: {
760: TS_RK *rk = (TS_RK *)ts->data;
761: RKTableau tab = rk->tableau;
762: PetscInt i;
764: PetscFunctionBegin;
765: PetscCall(MatDestroy(&rk->MatFwdSensip0));
766: if (rk->MatsFwdStageSensip) {
767: for (i = 0; i < tab->s; i++) PetscCall(MatDestroy(&rk->MatsFwdStageSensip[i]));
768: PetscCall(PetscFree(rk->MatsFwdStageSensip));
769: }
770: if (rk->MatsFwdSensipTemp) {
771: for (i = 0; i < tab->s; i++) PetscCall(MatDestroy(&rk->MatsFwdSensipTemp[i]));
772: PetscCall(PetscFree(rk->MatsFwdSensipTemp));
773: }
774: PetscCall(VecDestroy(&rk->VecDeltaFwdSensipCol));
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: static PetscErrorCode TSStep_RK(TS ts)
779: {
780: TS_RK *rk = (TS_RK *)ts->data;
781: RKTableau tab = rk->tableau;
782: const PetscInt s = tab->s;
783: const PetscReal *A = tab->A, *c = tab->c;
784: PetscScalar *w = rk->work;
785: Vec *Y = rk->Y, *YdotRHS = rk->YdotRHS;
786: PetscBool FSAL = (PetscBool)(tab->FSAL && !rk->newtableau);
787: TSAdapt adapt;
788: PetscInt i, j;
789: PetscInt rejections = 0;
790: PetscBool stageok, accept = PETSC_TRUE;
791: PetscReal next_time_step = ts->time_step;
793: PetscFunctionBegin;
794: if (ts->steprollback || ts->steprestart) FSAL = PETSC_FALSE;
795: if (FSAL) PetscCall(VecCopy(YdotRHS[s - 1], YdotRHS[0]));
796: rk->newtableau = PETSC_FALSE;
798: rk->status = TS_STEP_INCOMPLETE;
799: while (!ts->reason && rk->status != TS_STEP_COMPLETE) {
800: PetscReal t = ts->ptime;
801: PetscReal h = ts->time_step;
802: for (i = 0; i < s; i++) {
803: rk->stage_time = t + h * c[i];
804: PetscCall(TSPreStage(ts, rk->stage_time));
805: PetscCall(VecCopy(ts->vec_sol, Y[i]));
806: for (j = 0; j < i; j++) w[j] = h * A[i * s + j];
807: PetscCall(VecMAXPY(Y[i], i, w, YdotRHS));
808: PetscCall(TSPostStage(ts, rk->stage_time, i, Y));
809: PetscCall(TSGetAdapt(ts, &adapt));
810: PetscCall(TSAdaptCheckStage(adapt, ts, rk->stage_time, Y[i], &stageok));
811: if (!stageok) goto reject_step;
812: if (FSAL && !i) continue;
813: PetscCall(TSComputeRHSFunction(ts, t + h * c[i], Y[i], YdotRHS[i]));
814: }
816: rk->status = TS_STEP_INCOMPLETE;
817: PetscCall(TSEvaluateStep(ts, tab->order, ts->vec_sol, NULL));
818: rk->status = TS_STEP_PENDING;
819: PetscCall(TSGetAdapt(ts, &adapt));
820: PetscCall(TSAdaptCandidatesClear(adapt));
821: PetscCall(TSAdaptCandidateAdd(adapt, tab->name, tab->order, 1, tab->ccfl, (PetscReal)tab->s, PETSC_TRUE));
822: PetscCall(TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept));
823: rk->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
824: if (!accept) { /* Roll back the current step */
825: PetscCall(TSRollBack_RK(ts));
826: ts->time_step = next_time_step;
827: goto reject_step;
828: }
830: if (ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */
831: rk->ptime = ts->ptime;
832: rk->time_step = ts->time_step;
833: }
835: ts->ptime += ts->time_step;
836: ts->time_step = next_time_step;
837: break;
839: reject_step:
840: ts->reject++;
841: accept = PETSC_FALSE;
842: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
843: ts->reason = TS_DIVERGED_STEP_REJECTED;
844: PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", step rejections %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, rejections));
845: }
846: }
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: static PetscErrorCode TSAdjointSetUp_RK(TS ts)
851: {
852: TS_RK *rk = (TS_RK *)ts->data;
853: RKTableau tab = rk->tableau;
854: PetscInt s = tab->s;
856: PetscFunctionBegin;
857: if (ts->adjointsetupcalled) PetscFunctionReturn(PETSC_SUCCESS);
858: ts->adjointsetupcalled = PETSC_TRUE;
859: PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], s * ts->numcost, &rk->VecsDeltaLam));
860: PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], ts->numcost, &rk->VecsSensiTemp));
861: if (ts->vecs_sensip) PetscCall(VecDuplicate(ts->vecs_sensip[0], &rk->VecDeltaMu));
862: if (ts->vecs_sensi2) {
863: PetscCall(VecDuplicateVecs(ts->vecs_sensi[0], s * ts->numcost, &rk->VecsDeltaLam2));
864: PetscCall(VecDuplicateVecs(ts->vecs_sensi2[0], ts->numcost, &rk->VecsSensi2Temp));
865: }
866: if (ts->vecs_sensi2p) PetscCall(VecDuplicate(ts->vecs_sensi2p[0], &rk->VecDeltaMu2));
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: /*
871: Assumptions:
872: - TSStep_RK() always evaluates the step with b, not bembed.
873: */
874: static PetscErrorCode TSAdjointStep_RK(TS ts)
875: {
876: TS_RK *rk = (TS_RK *)ts->data;
877: TS quadts = ts->quadraturets;
878: RKTableau tab = rk->tableau;
879: Mat J, Jpre, Jquad;
880: const PetscInt s = tab->s;
881: const PetscReal *A = tab->A, *b = tab->b, *c = tab->c;
882: PetscScalar *w = rk->work, *xarr;
883: Vec *Y = rk->Y, *VecsDeltaLam = rk->VecsDeltaLam, VecDeltaMu = rk->VecDeltaMu, *VecsSensiTemp = rk->VecsSensiTemp;
884: Vec *VecsDeltaLam2 = rk->VecsDeltaLam2, VecDeltaMu2 = rk->VecDeltaMu2, *VecsSensi2Temp = rk->VecsSensi2Temp;
885: Vec VecDRDUTransCol = ts->vec_drdu_col, VecDRDPTransCol = ts->vec_drdp_col;
886: PetscInt i, j, nadj;
887: PetscReal t = ts->ptime;
888: PetscReal h = ts->time_step;
890: PetscFunctionBegin;
891: rk->status = TS_STEP_INCOMPLETE;
893: PetscCall(TSGetRHSJacobian(ts, &J, &Jpre, NULL, NULL));
894: if (quadts) PetscCall(TSGetRHSJacobian(quadts, &Jquad, NULL, NULL, NULL));
895: for (i = s - 1; i >= 0; i--) {
896: if (tab->FSAL && i == s - 1) {
897: /* VecsDeltaLam[nadj*s+s-1] are initialized with zeros and the values never change.*/
898: continue;
899: }
900: rk->stage_time = t + h * (1.0 - c[i]);
901: PetscCall(TSComputeSNESJacobian(ts, Y[i], J, Jpre));
902: if (quadts) PetscCall(TSComputeRHSJacobian(quadts, rk->stage_time, Y[i], Jquad, Jquad)); /* get r_u^T */
903: if (ts->vecs_sensip) {
904: PetscCall(TSComputeRHSJacobianP(ts, rk->stage_time, Y[i], ts->Jacprhs)); /* get f_p */
905: if (quadts) PetscCall(TSComputeRHSJacobianP(quadts, rk->stage_time, Y[i], quadts->Jacprhs)); /* get f_p for the quadrature */
906: }
908: if (b[i]) {
909: for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i] / b[i]; /* coefficients for computing VecsSensiTemp */
910: } else {
911: for (j = i + 1; j < s; j++) w[j - i - 1] = A[j * s + i]; /* coefficients for computing VecsSensiTemp */
912: }
914: for (nadj = 0; nadj < ts->numcost; nadj++) {
915: /* Stage values of lambda */
916: if (b[i]) {
917: /* lambda_{n+1} + \sum_{j=i+1}^s a_{ji}/b[i]*lambda_{s,j} */
918: PetscCall(VecCopy(ts->vecs_sensi[nadj], VecsSensiTemp[nadj])); /* VecDeltaLam is an vec array of size s by numcost */
919: PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
920: PetscCall(MatMultTranspose(J, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i])); /* VecsSensiTemp will be reused by 2nd-order adjoint */
921: PetscCall(VecScale(VecsDeltaLam[nadj * s + i], -h * b[i]));
922: if (quadts) {
923: PetscCall(MatDenseGetColumn(Jquad, nadj, &xarr));
924: PetscCall(VecPlaceArray(VecDRDUTransCol, xarr));
925: PetscCall(VecAXPY(VecsDeltaLam[nadj * s + i], -h * b[i], VecDRDUTransCol));
926: PetscCall(VecResetArray(VecDRDUTransCol));
927: PetscCall(MatDenseRestoreColumn(Jquad, &xarr));
928: }
929: } else {
930: /* \sum_{j=i+1}^s a_{ji}*lambda_{s,j} */
931: PetscCall(VecSet(VecsSensiTemp[nadj], 0));
932: PetscCall(VecMAXPY(VecsSensiTemp[nadj], s - i - 1, w, &VecsDeltaLam[nadj * s + i + 1]));
933: PetscCall(MatMultTranspose(J, VecsSensiTemp[nadj], VecsDeltaLam[nadj * s + i]));
934: PetscCall(VecScale(VecsDeltaLam[nadj * s + i], -h));
935: }
937: /* Stage values of mu */
938: if (ts->vecs_sensip) {
939: PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensiTemp[nadj], VecDeltaMu));
940: if (b[i]) {
941: PetscCall(VecScale(VecDeltaMu, -h * b[i]));
942: if (quadts) {
943: PetscCall(MatDenseGetColumn(quadts->Jacprhs, nadj, &xarr));
944: PetscCall(VecPlaceArray(VecDRDPTransCol, xarr));
945: PetscCall(VecAXPY(VecDeltaMu, -h * b[i], VecDRDPTransCol));
946: PetscCall(VecResetArray(VecDRDPTransCol));
947: PetscCall(MatDenseRestoreColumn(quadts->Jacprhs, &xarr));
948: }
949: } else {
950: PetscCall(VecScale(VecDeltaMu, -h));
951: }
952: PetscCall(VecAXPY(ts->vecs_sensip[nadj], 1., VecDeltaMu)); /* update sensip for each stage */
953: }
954: }
956: if (ts->vecs_sensi2 && ts->forward_solve) { /* 2nd-order adjoint, TLM mode has to be turned on */
957: /* Get w1 at t_{n+1} from TLM matrix */
958: PetscCall(MatDenseGetColumn(rk->MatsFwdStageSensip[i], 0, &xarr));
959: PetscCall(VecPlaceArray(ts->vec_sensip_col, xarr));
960: /* lambda_s^T F_UU w_1 */
961: PetscCall(TSComputeRHSHessianProductFunctionUU(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_sensip_col, ts->vecs_guu));
962: if (quadts) {
963: /* R_UU w_1 */
964: PetscCall(TSComputeRHSHessianProductFunctionUU(quadts, rk->stage_time, Y[i], NULL, ts->vec_sensip_col, ts->vecs_guu));
965: }
966: if (ts->vecs_sensip) {
967: /* lambda_s^T F_UP w_2 */
968: PetscCall(TSComputeRHSHessianProductFunctionUP(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_dir, ts->vecs_gup));
969: if (quadts) {
970: /* R_UP w_2 */
971: PetscCall(TSComputeRHSHessianProductFunctionUP(quadts, rk->stage_time, Y[i], NULL, ts->vec_sensip_col, ts->vecs_gup));
972: }
973: }
974: if (ts->vecs_sensi2p) {
975: /* lambda_s^T F_PU w_1 */
976: PetscCall(TSComputeRHSHessianProductFunctionPU(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_sensip_col, ts->vecs_gpu));
977: /* lambda_s^T F_PP w_2 */
978: PetscCall(TSComputeRHSHessianProductFunctionPP(ts, rk->stage_time, Y[i], VecsSensiTemp, ts->vec_dir, ts->vecs_gpp));
979: if (b[i] && quadts) {
980: /* R_PU w_1 */
981: PetscCall(TSComputeRHSHessianProductFunctionPU(quadts, rk->stage_time, Y[i], NULL, ts->vec_sensip_col, ts->vecs_gpu));
982: /* R_PP w_2 */
983: PetscCall(TSComputeRHSHessianProductFunctionPP(quadts, rk->stage_time, Y[i], NULL, ts->vec_dir, ts->vecs_gpp));
984: }
985: }
986: PetscCall(VecResetArray(ts->vec_sensip_col));
987: PetscCall(MatDenseRestoreColumn(rk->MatsFwdStageSensip[i], &xarr));
989: for (nadj = 0; nadj < ts->numcost; nadj++) {
990: /* Stage values of lambda */
991: if (b[i]) {
992: /* J_i^T*(Lambda_{n+1}+\sum_{j=i+1}^s a_{ji}/b_i*Lambda_{s,j} */
993: PetscCall(VecCopy(ts->vecs_sensi2[nadj], VecsSensi2Temp[nadj]));
994: PetscCall(VecMAXPY(VecsSensi2Temp[nadj], s - i - 1, w, &VecsDeltaLam2[nadj * s + i + 1]));
995: PetscCall(MatMultTranspose(J, VecsSensi2Temp[nadj], VecsDeltaLam2[nadj * s + i]));
996: PetscCall(VecScale(VecsDeltaLam2[nadj * s + i], -h * b[i]));
997: PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h * b[i], ts->vecs_guu[nadj]));
998: if (ts->vecs_sensip) PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h * b[i], ts->vecs_gup[nadj]));
999: } else {
1000: /* \sum_{j=i+1}^s a_{ji}*Lambda_{s,j} */
1001: PetscCall(VecSet(VecsDeltaLam2[nadj * s + i], 0));
1002: PetscCall(VecMAXPY(VecsSensi2Temp[nadj], s - i - 1, w, &VecsDeltaLam2[nadj * s + i + 1]));
1003: PetscCall(MatMultTranspose(J, VecsSensi2Temp[nadj], VecsDeltaLam2[nadj * s + i]));
1004: PetscCall(VecScale(VecsDeltaLam2[nadj * s + i], -h));
1005: PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h, ts->vecs_guu[nadj]));
1006: if (ts->vecs_sensip) PetscCall(VecAXPY(VecsDeltaLam2[nadj * s + i], -h, ts->vecs_gup[nadj]));
1007: }
1008: if (ts->vecs_sensi2p) { /* 2nd-order adjoint for parameters */
1009: PetscCall(MatMultTranspose(ts->Jacprhs, VecsSensi2Temp[nadj], VecDeltaMu2));
1010: if (b[i]) {
1011: PetscCall(VecScale(VecDeltaMu2, -h * b[i]));
1012: PetscCall(VecAXPY(VecDeltaMu2, -h * b[i], ts->vecs_gpu[nadj]));
1013: PetscCall(VecAXPY(VecDeltaMu2, -h * b[i], ts->vecs_gpp[nadj]));
1014: } else {
1015: PetscCall(VecScale(VecDeltaMu2, -h));
1016: PetscCall(VecAXPY(VecDeltaMu2, -h, ts->vecs_gpu[nadj]));
1017: PetscCall(VecAXPY(VecDeltaMu2, -h, ts->vecs_gpp[nadj]));
1018: }
1019: PetscCall(VecAXPY(ts->vecs_sensi2p[nadj], 1, VecDeltaMu2)); /* update sensi2p for each stage */
1020: }
1021: }
1022: }
1023: }
1025: for (j = 0; j < s; j++) w[j] = 1.0;
1026: for (nadj = 0; nadj < ts->numcost; nadj++) { /* no need to do this for mu's */
1027: PetscCall(VecMAXPY(ts->vecs_sensi[nadj], s, w, &VecsDeltaLam[nadj * s]));
1028: if (ts->vecs_sensi2) PetscCall(VecMAXPY(ts->vecs_sensi2[nadj], s, w, &VecsDeltaLam2[nadj * s]));
1029: }
1030: rk->status = TS_STEP_COMPLETE;
1031: PetscFunctionReturn(PETSC_SUCCESS);
1032: }
1034: static PetscErrorCode TSAdjointReset_RK(TS ts)
1035: {
1036: TS_RK *rk = (TS_RK *)ts->data;
1037: RKTableau tab = rk->tableau;
1039: PetscFunctionBegin;
1040: PetscCall(VecDestroyVecs(tab->s * ts->numcost, &rk->VecsDeltaLam));
1041: PetscCall(VecDestroyVecs(ts->numcost, &rk->VecsSensiTemp));
1042: PetscCall(VecDestroy(&rk->VecDeltaMu));
1043: PetscCall(VecDestroyVecs(tab->s * ts->numcost, &rk->VecsDeltaLam2));
1044: PetscCall(VecDestroy(&rk->VecDeltaMu2));
1045: PetscCall(VecDestroyVecs(ts->numcost, &rk->VecsSensi2Temp));
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: static PetscErrorCode TSInterpolate_RK(TS ts, PetscReal itime, Vec X)
1050: {
1051: TS_RK *rk = (TS_RK *)ts->data;
1052: PetscInt s = rk->tableau->s, p = rk->tableau->p, i, j;
1053: PetscReal h;
1054: PetscReal tt, t;
1055: PetscScalar *b;
1056: const PetscReal *B = rk->tableau->binterp;
1058: PetscFunctionBegin;
1059: PetscCheck(B, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "TSRK %s does not have an interpolation formula", rk->tableau->name);
1061: switch (rk->status) {
1062: case TS_STEP_INCOMPLETE:
1063: case TS_STEP_PENDING:
1064: h = ts->time_step;
1065: t = (itime - ts->ptime) / h;
1066: break;
1067: case TS_STEP_COMPLETE:
1068: h = ts->ptime - ts->ptime_prev;
1069: t = (itime - ts->ptime) / h + 1; /* In the interval [0,1] */
1070: break;
1071: default:
1072: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_PLIB, "Invalid TSStepStatus");
1073: }
1074: PetscCall(PetscMalloc1(s, &b));
1075: for (i = 0; i < s; i++) b[i] = 0;
1076: for (j = 0, tt = t; j < p; j++, tt *= t) {
1077: for (i = 0; i < s; i++) b[i] += h * B[i * p + j] * tt;
1078: }
1079: PetscCall(VecCopy(rk->Y[0], X));
1080: PetscCall(VecMAXPY(X, s, b, rk->YdotRHS));
1081: PetscCall(PetscFree(b));
1082: PetscFunctionReturn(PETSC_SUCCESS);
1083: }
1085: static PetscErrorCode TSRKTableauReset(TS ts)
1086: {
1087: TS_RK *rk = (TS_RK *)ts->data;
1088: RKTableau tab = rk->tableau;
1090: PetscFunctionBegin;
1091: if (!tab) PetscFunctionReturn(PETSC_SUCCESS);
1092: PetscCall(PetscFree(rk->work));
1093: PetscCall(VecDestroyVecs(tab->s, &rk->Y));
1094: PetscCall(VecDestroyVecs(tab->s, &rk->YdotRHS));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: static PetscErrorCode TSReset_RK(TS ts)
1099: {
1100: PetscFunctionBegin;
1101: PetscCall(TSRKTableauReset(ts));
1102: if (ts->use_splitrhsfunction) {
1103: PetscTryMethod(ts, "TSReset_RK_MultirateSplit_C", (TS), (ts));
1104: } else {
1105: PetscTryMethod(ts, "TSReset_RK_MultirateNonsplit_C", (TS), (ts));
1106: }
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1110: static PetscErrorCode DMCoarsenHook_TSRK(DM fine, DM coarse, PetscCtx ctx)
1111: {
1112: PetscFunctionBegin;
1113: PetscFunctionReturn(PETSC_SUCCESS);
1114: }
1116: static PetscErrorCode DMRestrictHook_TSRK(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
1117: {
1118: PetscFunctionBegin;
1119: PetscFunctionReturn(PETSC_SUCCESS);
1120: }
1122: static PetscErrorCode DMSubDomainHook_TSRK(DM dm, DM subdm, PetscCtx ctx)
1123: {
1124: PetscFunctionBegin;
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: static PetscErrorCode DMSubDomainRestrictHook_TSRK(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
1129: {
1130: PetscFunctionBegin;
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: }
1134: static PetscErrorCode TSRKTableauSetUp(TS ts)
1135: {
1136: TS_RK *rk = (TS_RK *)ts->data;
1137: RKTableau tab = rk->tableau;
1139: PetscFunctionBegin;
1140: PetscCall(PetscMalloc1(tab->s, &rk->work));
1141: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &rk->Y));
1142: PetscCall(VecDuplicateVecs(ts->vec_sol, tab->s, &rk->YdotRHS));
1143: rk->newtableau = PETSC_TRUE;
1144: PetscFunctionReturn(PETSC_SUCCESS);
1145: }
1147: static PetscErrorCode TSSetUp_RK(TS ts)
1148: {
1149: TS quadts = ts->quadraturets;
1150: DM dm;
1152: PetscFunctionBegin;
1153: PetscCall(TSCheckImplicitTerm(ts));
1154: PetscCall(TSRKTableauSetUp(ts));
1155: if (quadts && ts->costintegralfwd) {
1156: Mat Jquad;
1157: PetscCall(TSGetRHSJacobian(quadts, &Jquad, NULL, NULL, NULL));
1158: }
1159: PetscCall(TSGetDM(ts, &dm));
1160: PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSRK, DMRestrictHook_TSRK, ts));
1161: PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSRK, DMSubDomainRestrictHook_TSRK, ts));
1162: if (ts->use_splitrhsfunction) {
1163: PetscTryMethod(ts, "TSSetUp_RK_MultirateSplit_C", (TS), (ts));
1164: } else {
1165: PetscTryMethod(ts, "TSSetUp_RK_MultirateNonsplit_C", (TS), (ts));
1166: }
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: static PetscErrorCode TSSetFromOptions_RK(TS ts, PetscOptionItems PetscOptionsObject)
1171: {
1172: TS_RK *rk = (TS_RK *)ts->data;
1174: PetscFunctionBegin;
1175: PetscOptionsHeadBegin(PetscOptionsObject, "RK ODE solver options");
1176: {
1177: RKTableauLink link;
1178: PetscInt count, choice;
1179: PetscBool flg, use_multirate = PETSC_FALSE;
1180: const char **namelist;
1182: for (link = RKTableauList, count = 0; link; link = link->next, count++);
1183: PetscCall(PetscMalloc1(count, (char ***)&namelist));
1184: for (link = RKTableauList, count = 0; link; link = link->next, count++) namelist[count] = link->tab.name;
1185: PetscCall(PetscOptionsBool("-ts_rk_multirate", "Use interpolation-based multirate RK method", "TSRKSetMultirate", rk->use_multirate, &use_multirate, &flg));
1186: if (flg) PetscCall(TSRKSetMultirate(ts, use_multirate));
1187: PetscCall(PetscOptionsEList("-ts_rk_type", "Family of RK method", "TSRKSetType", (const char *const *)namelist, count, rk->tableau->name, &choice, &flg));
1188: if (flg) PetscCall(TSRKSetType(ts, namelist[choice]));
1189: PetscCall(PetscFree(namelist));
1190: }
1191: PetscOptionsHeadEnd();
1192: PetscOptionsBegin(PetscObjectComm((PetscObject)ts), NULL, "Multirate methods options", "");
1193: PetscCall(PetscOptionsInt("-ts_rk_dtratio", "time step ratio between slow and fast", "", rk->dtratio, &rk->dtratio, NULL));
1194: PetscOptionsEnd();
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: static PetscErrorCode TSView_RK(TS ts, PetscViewer viewer)
1199: {
1200: TS_RK *rk = (TS_RK *)ts->data;
1201: PetscBool isascii;
1203: PetscFunctionBegin;
1204: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1205: if (isascii) {
1206: RKTableau tab = rk->tableau;
1207: TSRKType rktype;
1208: const PetscReal *c;
1209: PetscInt s;
1210: char buf[512];
1211: PetscBool FSAL;
1213: PetscCall(TSRKGetType(ts, &rktype));
1214: PetscCall(TSRKGetTableau(ts, &s, NULL, NULL, &c, NULL, NULL, NULL, &FSAL));
1215: PetscCall(PetscViewerASCIIPrintf(viewer, " RK type %s\n", rktype));
1216: PetscCall(PetscViewerASCIIPrintf(viewer, " Order: %" PetscInt_FMT "\n", tab->order));
1217: PetscCall(PetscViewerASCIIPrintf(viewer, " FSAL property: %s\n", FSAL ? "yes" : "no"));
1218: PetscCall(PetscFormatRealArray(buf, sizeof(buf), "% 8.6f", s, c));
1219: PetscCall(PetscViewerASCIIPrintf(viewer, " Abscissa c = %s\n", buf));
1220: }
1221: PetscFunctionReturn(PETSC_SUCCESS);
1222: }
1224: static PetscErrorCode TSLoad_RK(TS ts, PetscViewer viewer)
1225: {
1226: TSAdapt adapt;
1228: PetscFunctionBegin;
1229: PetscCall(TSGetAdapt(ts, &adapt));
1230: PetscCall(TSAdaptLoad(adapt, viewer));
1231: PetscFunctionReturn(PETSC_SUCCESS);
1232: }
1234: /*@
1235: TSRKGetOrder - Get the order of the `TSRK` scheme
1237: Not Collective
1239: Input Parameter:
1240: . ts - timestepping context
1242: Output Parameter:
1243: . order - order of `TSRK` scheme
1245: Level: intermediate
1247: .seealso: [](ch_ts), `TSRK`, `TSRKGetType()`
1248: @*/
1249: PetscErrorCode TSRKGetOrder(TS ts, PetscInt *order)
1250: {
1251: PetscFunctionBegin;
1253: PetscAssertPointer(order, 2);
1254: PetscUseMethod(ts, "TSRKGetOrder_C", (TS, PetscInt *), (ts, order));
1255: PetscFunctionReturn(PETSC_SUCCESS);
1256: }
1258: /*@
1259: TSRKSetType - Set the type of the `TSRK` scheme
1261: Logically Collective
1263: Input Parameters:
1264: + ts - timestepping context
1265: - rktype - type of `TSRK` scheme
1267: Options Database Key:
1268: . -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs>
1270: Level: intermediate
1272: .seealso: [](ch_ts), `TSRKGetType()`, `TSRK`, `TSRKType`, `TSRK1FE`, `TSRK2A`, `TSRK2B`, `TSRK3`, `TSRK3BS`, `TSRK4`, `TSRK5F`, `TSRK5DP`, `TSRK5BS`, `TSRK6VR`, `TSRK7VR`, `TSRK8VR`
1273: @*/
1274: PetscErrorCode TSRKSetType(TS ts, TSRKType rktype)
1275: {
1276: PetscFunctionBegin;
1278: PetscAssertPointer(rktype, 2);
1279: PetscTryMethod(ts, "TSRKSetType_C", (TS, TSRKType), (ts, rktype));
1280: PetscFunctionReturn(PETSC_SUCCESS);
1281: }
1283: /*@
1284: TSRKGetType - Get the type of `TSRK` scheme
1286: Not Collective
1288: Input Parameter:
1289: . ts - timestepping context
1291: Output Parameter:
1292: . rktype - type of `TSRK`-scheme
1294: Level: intermediate
1296: .seealso: [](ch_ts), `TSRKSetType()`
1297: @*/
1298: PetscErrorCode TSRKGetType(TS ts, TSRKType *rktype)
1299: {
1300: PetscFunctionBegin;
1302: PetscUseMethod(ts, "TSRKGetType_C", (TS, TSRKType *), (ts, rktype));
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: static PetscErrorCode TSRKGetOrder_RK(TS ts, PetscInt *order)
1307: {
1308: TS_RK *rk = (TS_RK *)ts->data;
1310: PetscFunctionBegin;
1311: *order = rk->tableau->order;
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: static PetscErrorCode TSRKGetType_RK(TS ts, TSRKType *rktype)
1316: {
1317: TS_RK *rk = (TS_RK *)ts->data;
1319: PetscFunctionBegin;
1320: *rktype = rk->tableau->name;
1321: PetscFunctionReturn(PETSC_SUCCESS);
1322: }
1324: static PetscErrorCode TSRKSetType_RK(TS ts, TSRKType rktype)
1325: {
1326: TS_RK *rk = (TS_RK *)ts->data;
1327: PetscBool match;
1328: RKTableauLink link;
1330: PetscFunctionBegin;
1331: if (rk->tableau) {
1332: PetscCall(PetscStrcmp(rk->tableau->name, rktype, &match));
1333: if (match) PetscFunctionReturn(PETSC_SUCCESS);
1334: }
1335: for (link = RKTableauList; link; link = link->next) {
1336: PetscCall(PetscStrcmp(link->tab.name, rktype, &match));
1337: if (match) {
1338: if (ts->setupcalled) PetscCall(TSRKTableauReset(ts));
1339: rk->tableau = &link->tab;
1340: if (ts->setupcalled) PetscCall(TSRKTableauSetUp(ts));
1341: ts->default_adapt_type = rk->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1342: PetscFunctionReturn(PETSC_SUCCESS);
1343: }
1344: }
1345: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Could not find '%s'", rktype);
1346: }
1348: static PetscErrorCode TSGetStages_RK(TS ts, PetscInt *ns, Vec **Y)
1349: {
1350: TS_RK *rk = (TS_RK *)ts->data;
1352: PetscFunctionBegin;
1353: if (ns) *ns = rk->tableau->s;
1354: if (Y) *Y = rk->Y;
1355: PetscFunctionReturn(PETSC_SUCCESS);
1356: }
1358: static PetscErrorCode TSDestroy_RK(TS ts)
1359: {
1360: PetscFunctionBegin;
1361: PetscCall(TSReset_RK(ts));
1362: if (ts->dm) {
1363: PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSRK, DMRestrictHook_TSRK, ts));
1364: PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSRK, DMSubDomainRestrictHook_TSRK, ts));
1365: }
1366: PetscCall(PetscFree(ts->data));
1367: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetOrder_C", NULL));
1368: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetType_C", NULL));
1369: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetType_C", NULL));
1370: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetTableau_C", NULL));
1371: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetMultirate_C", NULL));
1372: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetMultirate_C", NULL));
1373: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateSplit_C", NULL));
1374: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateSplit_C", NULL));
1375: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSetUp_RK_MultirateNonsplit_C", NULL));
1376: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSReset_RK_MultirateNonsplit_C", NULL));
1377: PetscFunctionReturn(PETSC_SUCCESS);
1378: }
1380: /*
1381: This defines the nonlinear equation that is to be solved with SNES
1382: We do not need to solve the equation; we just use SNES to approximate the Jacobian
1383: */
1384: static PetscErrorCode SNESTSFormFunction_RK(SNES snes, Vec x, Vec y, TS ts)
1385: {
1386: TS_RK *rk = (TS_RK *)ts->data;
1387: DM dm, dmsave;
1389: PetscFunctionBegin;
1390: PetscCall(SNESGetDM(snes, &dm));
1391: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
1392: dmsave = ts->dm;
1393: ts->dm = dm;
1394: PetscCall(TSComputeRHSFunction(ts, rk->stage_time, x, y));
1395: ts->dm = dmsave;
1396: PetscFunctionReturn(PETSC_SUCCESS);
1397: }
1399: static PetscErrorCode SNESTSFormJacobian_RK(SNES snes, Vec x, Mat A, Mat B, TS ts)
1400: {
1401: TS_RK *rk = (TS_RK *)ts->data;
1402: DM dm, dmsave;
1404: PetscFunctionBegin;
1405: PetscCall(SNESGetDM(snes, &dm));
1406: dmsave = ts->dm;
1407: ts->dm = dm;
1408: PetscCall(TSComputeRHSJacobian(ts, rk->stage_time, x, A, B));
1409: ts->dm = dmsave;
1410: PetscFunctionReturn(PETSC_SUCCESS);
1411: }
1413: /*@
1414: TSRKSetMultirate - Use the interpolation-based multirate `TSRK` method
1416: Logically Collective
1418: Input Parameters:
1419: + ts - timestepping context
1420: - 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
1422: Options Database Key:
1423: . -ts_rk_multirate (true|false) - enable the multirate RK method
1425: Level: intermediate
1427: Note:
1428: 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).
1430: .seealso: [](ch_ts), `TSRK`, `TSRKGetMultirate()`
1431: @*/
1432: PetscErrorCode TSRKSetMultirate(TS ts, PetscBool use_multirate)
1433: {
1434: PetscFunctionBegin;
1435: PetscTryMethod(ts, "TSRKSetMultirate_C", (TS, PetscBool), (ts, use_multirate));
1436: PetscFunctionReturn(PETSC_SUCCESS);
1437: }
1439: /*@
1440: TSRKGetMultirate - Gets whether to use the interpolation-based multirate `TSRK` method
1442: Not Collective
1444: Input Parameter:
1445: . ts - timestepping context
1447: Output Parameter:
1448: . use_multirate - `PETSC_TRUE` if the multirate RK method is enabled, `PETSC_FALSE` otherwise
1450: Level: intermediate
1452: .seealso: [](ch_ts), `TSRK`, `TSRKSetMultirate()`
1453: @*/
1454: PetscErrorCode TSRKGetMultirate(TS ts, PetscBool *use_multirate)
1455: {
1456: PetscFunctionBegin;
1457: PetscUseMethod(ts, "TSRKGetMultirate_C", (TS, PetscBool *), (ts, use_multirate));
1458: PetscFunctionReturn(PETSC_SUCCESS);
1459: }
1461: /*MC
1462: TSRK - ODE and DAE solver using Runge-Kutta schemes
1464: The user should provide the right-hand side of the equation
1465: using `TSSetRHSFunction()`.
1467: Level: beginner
1469: Notes:
1470: The default is `TSRK3BS`, it can be changed with `TSRKSetType()` or -ts_rk_type
1472: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSRK`, `TSSetType()`, `TSRKSetType()`, `TSRKGetType()`, `TSRK2D`, `TSRK2E`, `TSRK3`,
1473: `TSRK4`, `TSRK5`, `TSRKPRSSP2`, `TSRKBPR3`, `TSRKType`, `TSRKRegister()`, `TSRKSetMultirate()`, `TSRKGetMultirate()`, `TSType`
1474: M*/
1475: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
1476: {
1477: TS_RK *rk;
1479: PetscFunctionBegin;
1480: PetscCall(TSRKInitializePackage());
1482: ts->ops->reset = TSReset_RK;
1483: ts->ops->destroy = TSDestroy_RK;
1484: ts->ops->view = TSView_RK;
1485: ts->ops->load = TSLoad_RK;
1486: ts->ops->setup = TSSetUp_RK;
1487: ts->ops->interpolate = TSInterpolate_RK;
1488: ts->ops->step = TSStep_RK;
1489: ts->ops->evaluatestep = TSEvaluateStep_RK;
1490: ts->ops->rollback = TSRollBack_RK;
1491: ts->ops->setfromoptions = TSSetFromOptions_RK;
1492: ts->ops->getstages = TSGetStages_RK;
1494: ts->ops->snesfunction = SNESTSFormFunction_RK;
1495: ts->ops->snesjacobian = SNESTSFormJacobian_RK;
1496: ts->ops->adjointintegral = TSAdjointCostIntegral_RK;
1497: ts->ops->adjointsetup = TSAdjointSetUp_RK;
1498: ts->ops->adjointstep = TSAdjointStep_RK;
1499: ts->ops->adjointreset = TSAdjointReset_RK;
1501: ts->ops->forwardintegral = TSForwardCostIntegral_RK;
1502: ts->ops->forwardsetup = TSForwardSetUp_RK;
1503: ts->ops->forwardreset = TSForwardReset_RK;
1504: ts->ops->forwardstep = TSForwardStep_RK;
1505: ts->ops->forwardgetstages = TSForwardGetStages_RK;
1507: PetscCall(PetscNew(&rk));
1508: ts->data = (void *)rk;
1510: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetOrder_C", TSRKGetOrder_RK));
1511: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetType_C", TSRKGetType_RK));
1512: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetType_C", TSRKSetType_RK));
1513: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetTableau_C", TSRKGetTableau_RK));
1514: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKSetMultirate_C", TSRKSetMultirate_RK));
1515: PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSRKGetMultirate_C", TSRKGetMultirate_RK));
1517: PetscCall(TSRKSetType(ts, TSRKDefault));
1518: rk->dtratio = 1;
1519: PetscFunctionReturn(PETSC_SUCCESS);
1520: }