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 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: `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: `PetscInitialize()`
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: `PetscFinalize()`
409: @*/
410: PetscErrorCode TSRKFinalizePackage(void)
411: {
412:   TSRKPackageInitialized = PETSC_FALSE;
413:   TSRKRegisterDestroy();
414:   return 0;
415: }

417: /*@C
418:    TSRKRegister - register an RK 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:    Notes:
434:    Several RK methods are provided, this function is only needed to create new methods.

436:    Level: advanced

438: .seealso: `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 RK 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: `TSRK`
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 RK scheme

1214:   Not collective

1216:   Input Parameter:
1217: .  ts - timestepping context

1219:   Output Parameter:
1220: .  order - order of RK-scheme

1222:   Level: intermediate

1224: .seealso: `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 RK scheme

1237:   Logically collective

1239:   Input Parameters:
1240: +  ts - timestepping context
1241: -  rktype - type of RK-scheme

1243:   Options Database:
1244: .   -ts_rk_type - <1fe,2a,3,3bs,4,5f,5dp,5bs>

1246:   Level: intermediate

1248: .seealso: `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 RK scheme

1261:   Not collective

1263:   Input Parameter:
1264: .  ts - timestepping context

1266:   Output Parameter:
1267: .  rktype - type of RK-scheme

1269:   Level: intermediate

1271: .seealso: `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 RK method

1383:   Logically collective

1385:   Input Parameters:
1386: +  ts - timestepping context
1387: -  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

1389:   Options Database:
1390: .   -ts_rk_multirate - <true,false>

1392:   Notes:
1393:   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).

1395:   Level: intermediate

1397: .seealso: `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 RK 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: `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:   Notes:
1433:   The default is TSRK3BS, it can be changed with TSRKSetType() or -ts_rk_type

1435:   Level: beginner

1437: .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSRKSetType()`, `TSRKGetType()`, `TSRK2D`, `TTSRK2E`, `TSRK3`,
1438:           `TSRK4`, `TSRK5`, `TSRKPRSSP2`, `TSRKBPR3`, `TSRKType`, `TSRKRegister()`, `TSRKSetMultirate()`, `TSRKGetMultirate()`

1440: M*/
1441: PETSC_EXTERN PetscErrorCode TSCreate_RK(TS ts)
1442: {
1443:   TS_RK *rk;

1445:   TSRKInitializePackage();

1447:   ts->ops->reset          = TSReset_RK;
1448:   ts->ops->destroy        = TSDestroy_RK;
1449:   ts->ops->view           = TSView_RK;
1450:   ts->ops->load           = TSLoad_RK;
1451:   ts->ops->setup          = TSSetUp_RK;
1452:   ts->ops->interpolate    = TSInterpolate_RK;
1453:   ts->ops->step           = TSStep_RK;
1454:   ts->ops->evaluatestep   = TSEvaluateStep_RK;
1455:   ts->ops->rollback       = TSRollBack_RK;
1456:   ts->ops->setfromoptions = TSSetFromOptions_RK;
1457:   ts->ops->getstages      = TSGetStages_RK;

1459:   ts->ops->snesfunction    = SNESTSFormFunction_RK;
1460:   ts->ops->snesjacobian    = SNESTSFormJacobian_RK;
1461:   ts->ops->adjointintegral = TSAdjointCostIntegral_RK;
1462:   ts->ops->adjointsetup    = TSAdjointSetUp_RK;
1463:   ts->ops->adjointstep     = TSAdjointStep_RK;
1464:   ts->ops->adjointreset    = TSAdjointReset_RK;

1466:   ts->ops->forwardintegral  = TSForwardCostIntegral_RK;
1467:   ts->ops->forwardsetup     = TSForwardSetUp_RK;
1468:   ts->ops->forwardreset     = TSForwardReset_RK;
1469:   ts->ops->forwardstep      = TSForwardStep_RK;
1470:   ts->ops->forwardgetstages = TSForwardGetStages_RK;

1472:   PetscNew(&rk);
1473:   ts->data = (void *)rk;

1475:   PetscObjectComposeFunction((PetscObject)ts, "TSRKGetOrder_C", TSRKGetOrder_RK);
1476:   PetscObjectComposeFunction((PetscObject)ts, "TSRKGetType_C", TSRKGetType_RK);
1477:   PetscObjectComposeFunction((PetscObject)ts, "TSRKSetType_C", TSRKSetType_RK);
1478:   PetscObjectComposeFunction((PetscObject)ts, "TSRKGetTableau_C", TSRKGetTableau_RK);
1479:   PetscObjectComposeFunction((PetscObject)ts, "TSRKSetMultirate_C", TSRKSetMultirate_RK);
1480:   PetscObjectComposeFunction((PetscObject)ts, "TSRKGetMultirate_C", TSRKGetMultirate_RK);

1482:   TSRKSetType(ts, TSRKDefault);
1483:   rk->dtratio = 1;
1484:   return 0;
1485: }