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