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