1 /*
2 Code for timestepping with additive Runge-Kutta IMEX method
3
4 Notes:
5 The general system is written as
6
7 F(t,U,Udot) = G(t,U)
8
9 where F represents the stiff part of the physics and G represents the non-stiff part.
10
11 */
12 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
13 #include <petscdm.h>
14
15 static TSARKIMEXType TSARKIMEXDefault = TSARKIMEX3;
16 static PetscBool TSARKIMEXRegisterAllCalled;
17 static PetscBool TSARKIMEXPackageInitialized;
18 static PetscErrorCode TSExtrapolate_ARKIMEX(TS,PetscReal,Vec);
19
20 typedef struct _ARKTableau *ARKTableau;
21 struct _ARKTableau {
22 char *name;
23 PetscInt order; /* Classical approximation order of the method */
24 PetscInt s; /* Number of stages */
25 PetscBool stiffly_accurate; /* The implicit part is stiffly accurate*/
26 PetscBool FSAL_implicit; /* The implicit part is FSAL*/
27 PetscBool explicit_first_stage; /* The implicit part has an explicit first stage*/
28 PetscInt pinterp; /* Interpolation order */
29 PetscReal *At,*bt,*ct; /* Stiff tableau */
30 PetscReal *A,*b,*c; /* Non-stiff tableau */
31 PetscReal *bembedt,*bembed; /* Embedded formula of order one less (order-1) */
32 PetscReal *binterpt,*binterp; /* Dense output formula */
33 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */
34 };
35 typedef struct _ARKTableauLink *ARKTableauLink;
36 struct _ARKTableauLink {
37 struct _ARKTableau tab;
38 ARKTableauLink next;
39 };
40 static ARKTableauLink ARKTableauList;
41
42 typedef struct {
43 ARKTableau tableau;
44 Vec *Y; /* States computed during the step */
45 Vec *YdotI; /* Time derivatives for the stiff part */
46 Vec *YdotRHS; /* Function evaluations for the non-stiff part */
47 Vec *Y_prev; /* States computed during the previous time step */
48 Vec *YdotI_prev; /* Time derivatives for the stiff part for the previous time step*/
49 Vec *YdotRHS_prev; /* Function evaluations for the non-stiff part for the previous time step*/
50 Vec Ydot0; /* Holds the slope from the previous step in FSAL case */
51 Vec Ydot; /* Work vector holding Ydot during residual evaluation */
52 Vec Z; /* Ydot = shift(Y-Z) */
53 PetscScalar *work; /* Scalar work */
54 PetscReal scoeff; /* shift = scoeff/dt */
55 PetscReal stage_time;
56 PetscBool imex;
57 PetscBool extrapolate; /* Extrapolate initial guess from previous time-step stage values */
58 TSStepStatus status;
59 } TS_ARKIMEX;
60 /*MC
61 TSARKIMEXARS122 - Second order ARK IMEX scheme.
62
63 This method has one explicit stage and one implicit stage.
64
65 Options Database:
66 . -ts_arkimex_type ars122
67
68 References:
69 . 1. - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997).
70
71 Level: advanced
72
73 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
74 M*/
75 /*MC
76 TSARKIMEXA2 - Second order ARK IMEX scheme with A-stable implicit part.
77
78 This method has an explicit stage and one implicit stage, and has an A-stable implicit scheme. This method was provided by Emil Constantinescu.
79
80 Options Database:
81 . -ts_arkimex_type a2
82
83 Level: advanced
84
85 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
86 M*/
87 /*MC
88 TSARKIMEXL2 - Second order ARK IMEX scheme with L-stable implicit part.
89
90 This method has two implicit stages, and L-stable implicit scheme.
91
92 Options Database:
93 . -ts_arkimex_type l2
94
95 References:
96 . 1. - L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005.
97
98 Level: advanced
99
100 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
101 M*/
102 /*MC
103 TSARKIMEX1BEE - First order Backward Euler represented as an ARK IMEX scheme with extrapolation as error estimator. This is a 3-stage method.
104
105 This method is aimed at starting the integration of implicit DAEs when explicit first-stage ARK methods are used.
106
107 Options Database:
108 . -ts_arkimex_type 1bee
109
110 Level: advanced
111
112 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
113 M*/
114 /*MC
115 TSARKIMEX2C - Second order ARK IMEX scheme with L-stable implicit part.
116
117 This method has one explicit stage and two implicit stages. The implicit part is the same as in TSARKIMEX2D and TSARKIMEX2E, but the explicit part has a larger stability region on the negative real axis. This method was provided by Emil Constantinescu.
118
119 Options Database:
120 . -ts_arkimex_type 2c
121
122 Level: advanced
123
124 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
125 M*/
126 /*MC
127 TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part.
128
129 This method has one explicit stage and two implicit stages. The stability function is independent of the explicit part in the infinity limit of the implict component. This method was provided by Emil Constantinescu.
130
131 Options Database:
132 . -ts_arkimex_type 2d
133
134 Level: advanced
135
136 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
137 M*/
138 /*MC
139 TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part.
140
141 This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu.
142
143 Options Database:
144 . -ts_arkimex_type 2e
145
146 Level: advanced
147
148 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
149 M*/
150 /*MC
151 TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme.
152
153 This method has three implicit stages.
154
155 References:
156 . 1. - L. Pareschi, G. Russo, Implicit Explicit Runge Kutta schemes and applications to hyperbolic systems with relaxations. Journal of Scientific Computing Volume: 25, Issue: 1, October, 2005.
157
158 This method is referred to as SSP2-(3,3,2) in https://arxiv.org/abs/1110.4375
159
160 Options Database:
161 . -ts_arkimex_type prssp2
162
163 Level: advanced
164
165 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
166 M*/
167 /*MC
168 TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part.
169
170 This method has one explicit stage and three implicit stages.
171
172 Options Database:
173 . -ts_arkimex_type 3
174
175 References:
176 . 1. - Kennedy and Carpenter 2003.
177
178 Level: advanced
179
180 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
181 M*/
182 /*MC
183 TSARKIMEXARS443 - Third order ARK IMEX scheme.
184
185 This method has one explicit stage and four implicit stages.
186
187 Options Database:
188 . -ts_arkimex_type ars443
189
190 References:
191 + 1. - U. Ascher, S. Ruuth, R. J. Spiteri, Implicit explicit Runge Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997).
192 - 2. - This method is referred to as ARS(4,4,3) in https://arxiv.org/abs/1110.4375
193
194 Level: advanced
195
196 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
197 M*/
198 /*MC
199 TSARKIMEXBPR3 - Third order ARK IMEX scheme.
200
201 This method has one explicit stage and four implicit stages.
202
203 Options Database:
204 . -ts_arkimex_type bpr3
205
206 References:
207 . This method is referred to as ARK3 in https://arxiv.org/abs/1110.4375
208
209 Level: advanced
210
211 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
212 M*/
213 /*MC
214 TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part.
215
216 This method has one explicit stage and four implicit stages.
217
218 Options Database:
219 . -ts_arkimex_type 4
220
221 References:
222 . Kennedy and Carpenter 2003.
223
224 Level: advanced
225
226 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
227 M*/
228 /*MC
229 TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part.
230
231 This method has one explicit stage and five implicit stages.
232
233 Options Database:
234 . -ts_arkimex_type 5
235
236 References:
237 . Kennedy and Carpenter 2003.
238
239 Level: advanced
240
241 .seealso: TSARKIMEX, TSARKIMEXType, TSARKIMEXSetType()
242 M*/
243
244 /*@C
245 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSARKIMEX
246
247 Not Collective, but should be called by all processes which will need the schemes to be registered
248
249 Level: advanced
250
251 .seealso: TSARKIMEXRegisterDestroy()
252 @*/
TSARKIMEXRegisterAll(void)253 PetscErrorCode TSARKIMEXRegisterAll(void)
254 {
255 PetscErrorCode ierr;
256
257 PetscFunctionBegin;
258 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(0);
259 TSARKIMEXRegisterAllCalled = PETSC_TRUE;
260
261 {
262 const PetscReal
263 A[3][3] = {{0.0,0.0,0.0},
264 {0.0,0.0,0.0},
265 {0.0,0.5,0.0}},
266 At[3][3] = {{1.0,0.0,0.0},
267 {0.0,0.5,0.0},
268 {0.0,0.5,0.5}},
269 b[3] = {0.0,0.5,0.5},
270 bembedt[3] = {1.0,0.0,0.0};
271 ierr = TSARKIMEXRegister(TSARKIMEX1BEE,2,3,&At[0][0],b,NULL,&A[0][0],b,NULL,bembedt,bembedt,1,b,NULL);CHKERRQ(ierr);
272 }
273 {
274 const PetscReal
275 A[2][2] = {{0.0,0.0},
276 {0.5,0.0}},
277 At[2][2] = {{0.0,0.0},
278 {0.0,0.5}},
279 b[2] = {0.0,1.0},
280 bembedt[2] = {0.5,0.5};
281 /* binterpt[2][2] = {{1.0,-1.0},{0.0,1.0}}; second order dense output has poor stability properties and hence it is not currently in use*/
282 ierr = TSARKIMEXRegister(TSARKIMEXARS122,2,2,&At[0][0],b,NULL,&A[0][0],b,NULL,bembedt,bembedt,1,b,NULL);CHKERRQ(ierr);
283 }
284 {
285 const PetscReal
286 A[2][2] = {{0.0,0.0},
287 {1.0,0.0}},
288 At[2][2] = {{0.0,0.0},
289 {0.5,0.5}},
290 b[2] = {0.5,0.5},
291 bembedt[2] = {0.0,1.0};
292 /* binterpt[2][2] = {{1.0,-0.5},{0.0,0.5}} second order dense output has poor stability properties and hence it is not currently in use*/
293 ierr = TSARKIMEXRegister(TSARKIMEXA2,2,2,&At[0][0],b,NULL,&A[0][0],b,NULL,bembedt,bembedt,1,b,NULL);CHKERRQ(ierr);
294 }
295 {
296 /* const PetscReal us2 = 1.0-1.0/PetscSqrtReal((PetscReal)2.0); Direct evaluation: 0.2928932188134524755992. Used below to ensure all values are available at compile time */
297 const PetscReal
298 A[2][2] = {{0.0,0.0},
299 {1.0,0.0}},
300 At[2][2] = {{0.2928932188134524755992,0.0},
301 {1.0-2.0*0.2928932188134524755992,0.2928932188134524755992}},
302 b[2] = {0.5,0.5},
303 bembedt[2] = {0.0,1.0},
304 binterpt[2][2] = {{ (0.2928932188134524755992-1.0)/(2.0*0.2928932188134524755992-1.0),-1/(2.0*(1.0-2.0*0.2928932188134524755992))},
305 {1-(0.2928932188134524755992-1.0)/(2.0*0.2928932188134524755992-1.0),-1/(2.0*(1.0-2.0*0.2928932188134524755992))}},
306 binterp[2][2] = {{1.0,-0.5},{0.0,0.5}};
307 ierr = TSARKIMEXRegister(TSARKIMEXL2,2,2,&At[0][0],b,NULL,&A[0][0],b,NULL,bembedt,bembedt,2,binterpt[0],binterp[0]);CHKERRQ(ierr);
308 }
309 {
310 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */
311 const PetscReal
312 A[3][3] = {{0,0,0},
313 {2-1.414213562373095048802,0,0},
314 {0.5,0.5,0}},
315 At[3][3] = {{0,0,0},
316 {1-1/1.414213562373095048802,1-1/1.414213562373095048802,0},
317 {1/(2*1.414213562373095048802),1/(2*1.414213562373095048802),1-1/1.414213562373095048802}},
318 bembedt[3] = {(4.-1.414213562373095048802)/8.,(4.-1.414213562373095048802)/8.,1/(2.*1.414213562373095048802)},
319 binterpt[3][2] = {{1.0/1.414213562373095048802,-1.0/(2.0*1.414213562373095048802)},
320 {1.0/1.414213562373095048802,-1.0/(2.0*1.414213562373095048802)},
321 {1.0-1.414213562373095048802,1.0/1.414213562373095048802}};
322 ierr = TSARKIMEXRegister(TSARKIMEX2C,2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,2,binterpt[0],NULL);CHKERRQ(ierr);
323 }
324 {
325 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */
326 const PetscReal
327 A[3][3] = {{0,0,0},
328 {2-1.414213562373095048802,0,0},
329 {0.75,0.25,0}},
330 At[3][3] = {{0,0,0},
331 {1-1/1.414213562373095048802,1-1/1.414213562373095048802,0},
332 {1/(2*1.414213562373095048802),1/(2*1.414213562373095048802),1-1/1.414213562373095048802}},
333 bembedt[3] = {(4.-1.414213562373095048802)/8.,(4.-1.414213562373095048802)/8.,1/(2.*1.414213562373095048802)},
334 binterpt[3][2] = {{1.0/1.414213562373095048802,-1.0/(2.0*1.414213562373095048802)},
335 {1.0/1.414213562373095048802,-1.0/(2.0*1.414213562373095048802)},
336 {1.0-1.414213562373095048802,1.0/1.414213562373095048802}};
337 ierr = TSARKIMEXRegister(TSARKIMEX2D,2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,2,binterpt[0],NULL);CHKERRQ(ierr);
338 }
339 { /* Optimal for linear implicit part */
340 /* const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), Direct evaluation: 1.414213562373095048802. Used below to ensure all values are available at compile time */
341 const PetscReal
342 A[3][3] = {{0,0,0},
343 {2-1.414213562373095048802,0,0},
344 {(3-2*1.414213562373095048802)/6,(3+2*1.414213562373095048802)/6,0}},
345 At[3][3] = {{0,0,0},
346 {1-1/1.414213562373095048802,1-1/1.414213562373095048802,0},
347 {1/(2*1.414213562373095048802),1/(2*1.414213562373095048802),1-1/1.414213562373095048802}},
348 bembedt[3] = {(4.-1.414213562373095048802)/8.,(4.-1.414213562373095048802)/8.,1/(2.*1.414213562373095048802)},
349 binterpt[3][2] = {{1.0/1.414213562373095048802,-1.0/(2.0*1.414213562373095048802)},
350 {1.0/1.414213562373095048802,-1.0/(2.0*1.414213562373095048802)},
351 {1.0-1.414213562373095048802,1.0/1.414213562373095048802}};
352 ierr = TSARKIMEXRegister(TSARKIMEX2E,2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,2,binterpt[0],NULL);CHKERRQ(ierr);
353 }
354 { /* Optimal for linear implicit part */
355 const PetscReal
356 A[3][3] = {{0,0,0},
357 {0.5,0,0},
358 {0.5,0.5,0}},
359 At[3][3] = {{0.25,0,0},
360 {0,0.25,0},
361 {1./3,1./3,1./3}};
362 ierr = TSARKIMEXRegister(TSARKIMEXPRSSP2,2,3,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,NULL,NULL,0,NULL,NULL);CHKERRQ(ierr);
363 }
364 {
365 const PetscReal
366 A[4][4] = {{0,0,0,0},
367 {1767732205903./2027836641118.,0,0,0},
368 {5535828885825./10492691773637.,788022342437./10882634858940.,0,0},
369 {6485989280629./16251701735622.,-4246266847089./9704473918619.,10755448449292./10357097424841.,0}},
370 At[4][4] = {{0,0,0,0},
371 {1767732205903./4055673282236.,1767732205903./4055673282236.,0,0},
372 {2746238789719./10658868560708.,-640167445237./6845629431997.,1767732205903./4055673282236.,0},
373 {1471266399579./7840856788654.,-4482444167858./7529755066697.,11266239266428./11593286722821.,1767732205903./4055673282236.}},
374 bembedt[4] = {2756255671327./12835298489170.,-10771552573575./22201958757719.,9247589265047./10645013368117.,2193209047091./5459859503100.},
375 binterpt[4][2] = {{4655552711362./22874653954995., -215264564351./13552729205753.},
376 {-18682724506714./9892148508045.,17870216137069./13817060693119.},
377 {34259539580243./13192909600954.,-28141676662227./17317692491321.},
378 {584795268549./6622622206610., 2508943948391./7218656332882.}};
379 ierr = TSARKIMEXRegister(TSARKIMEX3,3,4,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,2,binterpt[0],NULL);CHKERRQ(ierr);
380 }
381 {
382 const PetscReal
383 A[5][5] = {{0,0,0,0,0},
384 {1./2,0,0,0,0},
385 {11./18,1./18,0,0,0},
386 {5./6,-5./6,.5,0,0},
387 {1./4,7./4,3./4,-7./4,0}},
388 At[5][5] = {{0,0,0,0,0},
389 {0,1./2,0,0,0},
390 {0,1./6,1./2,0,0},
391 {0,-1./2,1./2,1./2,0},
392 {0,3./2,-3./2,1./2,1./2}},
393 *bembedt = NULL;
394 ierr = TSARKIMEXRegister(TSARKIMEXARS443,3,5,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,0,NULL,NULL);CHKERRQ(ierr);
395 }
396 {
397 const PetscReal
398 A[5][5] = {{0,0,0,0,0},
399 {1,0,0,0,0},
400 {4./9,2./9,0,0,0},
401 {1./4,0,3./4,0,0},
402 {1./4,0,3./5,0,0}},
403 At[5][5] = {{0,0,0,0,0},
404 {.5,.5,0,0,0},
405 {5./18,-1./9,.5,0,0},
406 {.5,0,0,.5,0},
407 {.25,0,.75,-.5,.5}},
408 *bembedt = NULL;
409 ierr = TSARKIMEXRegister(TSARKIMEXBPR3,3,5,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,0,NULL,NULL);CHKERRQ(ierr);
410 }
411 {
412 const PetscReal
413 A[6][6] = {{0,0,0,0,0,0},
414 {1./2,0,0,0,0,0},
415 {13861./62500.,6889./62500.,0,0,0,0},
416 {-116923316275./2393684061468.,-2731218467317./15368042101831.,9408046702089./11113171139209.,0,0,0},
417 {-451086348788./2902428689909.,-2682348792572./7519795681897.,12662868775082./11960479115383.,3355817975965./11060851509271.,0,0},
418 {647845179188./3216320057751.,73281519250./8382639484533.,552539513391./3454668386233.,3354512671639./8306763924573.,4040./17871.,0}},
419 At[6][6] = {{0,0,0,0,0,0},
420 {1./4,1./4,0,0,0,0},
421 {8611./62500.,-1743./31250.,1./4,0,0,0},
422 {5012029./34652500.,-654441./2922500.,174375./388108.,1./4,0,0},
423 {15267082809./155376265600.,-71443401./120774400.,730878875./902184768.,2285395./8070912.,1./4,0},
424 {82889./524892.,0,15625./83664.,69875./102672.,-2260./8211,1./4}},
425 bembedt[6] = {4586570599./29645900160.,0,178811875./945068544.,814220225./1159782912.,-3700637./11593932.,61727./225920.},
426 binterpt[6][3] = {{6943876665148./7220017795957.,-54480133./30881146.,6818779379841./7100303317025.},
427 {0,0,0},
428 {7640104374378./9702883013639.,-11436875./14766696.,2173542590792./12501825683035.},
429 {-20649996744609./7521556579894.,174696575./18121608.,-31592104683404./5083833661969.},
430 {8854892464581./2390941311638.,-12120380./966161.,61146701046299./7138195549469.},
431 {-11397109935349./6675773540249.,3843./706.,-17219254887155./4939391667607.}};
432 ierr = TSARKIMEXRegister(TSARKIMEX4,4,6,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,3,binterpt[0],NULL);CHKERRQ(ierr);
433 }
434 {
435 const PetscReal
436 A[8][8] = {{0,0,0,0,0,0,0,0},
437 {41./100,0,0,0,0,0,0,0},
438 {367902744464./2072280473677.,677623207551./8224143866563.,0,0,0,0,0,0},
439 {1268023523408./10340822734521.,0,1029933939417./13636558850479.,0,0,0,0,0},
440 {14463281900351./6315353703477.,0,66114435211212./5879490589093.,-54053170152839./4284798021562.,0,0,0,0},
441 {14090043504691./34967701212078.,0,15191511035443./11219624916014.,-18461159152457./12425892160975.,-281667163811./9011619295870.,0,0,0},
442 {19230459214898./13134317526959.,0,21275331358303./2942455364971.,-38145345988419./4862620318723.,-1./8,-1./8,0,0},
443 {-19977161125411./11928030595625.,0,-40795976796054./6384907823539.,177454434618887./12078138498510.,782672205425./8267701900261.,-69563011059811./9646580694205.,7356628210526./4942186776405.,0}},
444 At[8][8] = {{0,0,0,0,0,0,0,0},
445 {41./200.,41./200.,0,0,0,0,0,0},
446 {41./400.,-567603406766./11931857230679.,41./200.,0,0,0,0,0},
447 {683785636431./9252920307686.,0,-110385047103./1367015193373.,41./200.,0,0,0,0},
448 {3016520224154./10081342136671.,0,30586259806659./12414158314087.,-22760509404356./11113319521817.,41./200.,0,0,0},
449 {218866479029./1489978393911.,0,638256894668./5436446318841.,-1179710474555./5321154724896.,-60928119172./8023461067671.,41./200.,0,0},
450 {1020004230633./5715676835656.,0,25762820946817./25263940353407.,-2161375909145./9755907335909.,-211217309593./5846859502534.,-4269925059573./7827059040749.,41./200,0},
451 {-872700587467./9133579230613.,0,0,22348218063261./9555858737531.,-1143369518992./8141816002931.,-39379526789629./19018526304540.,32727382324388./42900044865799.,41./200.}},
452 bembedt[8] = {-975461918565./9796059967033.,0,0,78070527104295./32432590147079.,-548382580838./3424219808633.,-33438840321285./15594753105479.,3629800801594./4656183773603.,4035322873751./18575991585200.},
453 binterpt[8][3] = {{-17674230611817./10670229744614., 43486358583215./12773830924787., -9257016797708./5021505065439.},
454 {0, 0, 0 },
455 {0, 0, 0 },
456 {65168852399939./7868540260826., -91478233927265./11067650958493., 26096422576131./11239449250142.},
457 {15494834004392./5936557850923., -79368583304911./10890268929626., 92396832856987./20362823103730.},
458 {-99329723586156./26959484932159., -12239297817655./9152339842473., 30029262896817./10175596800299.},
459 {-19024464361622./5461577185407., 115839755401235./10719374521269., -26136350496073./3983972220547.},
460 {-6511271360970./6095937251113., 5843115559534./2180450260947., -5289405421727./3760307252460. }};
461 ierr = TSARKIMEXRegister(TSARKIMEX5,5,8,&At[0][0],NULL,NULL,&A[0][0],NULL,NULL,bembedt,bembedt,3,binterpt[0],NULL);CHKERRQ(ierr);
462 }
463 PetscFunctionReturn(0);
464 }
465
466 /*@C
467 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by TSARKIMEXRegister().
468
469 Not Collective
470
471 Level: advanced
472
473 .seealso: TSARKIMEXRegister(), TSARKIMEXRegisterAll()
474 @*/
TSARKIMEXRegisterDestroy(void)475 PetscErrorCode TSARKIMEXRegisterDestroy(void)
476 {
477 PetscErrorCode ierr;
478 ARKTableauLink link;
479
480 PetscFunctionBegin;
481 while ((link = ARKTableauList)) {
482 ARKTableau t = &link->tab;
483 ARKTableauList = link->next;
484 ierr = PetscFree6(t->At,t->bt,t->ct,t->A,t->b,t->c);CHKERRQ(ierr);
485 ierr = PetscFree2(t->bembedt,t->bembed);CHKERRQ(ierr);
486 ierr = PetscFree2(t->binterpt,t->binterp);CHKERRQ(ierr);
487 ierr = PetscFree(t->name);CHKERRQ(ierr);
488 ierr = PetscFree(link);CHKERRQ(ierr);
489 }
490 TSARKIMEXRegisterAllCalled = PETSC_FALSE;
491 PetscFunctionReturn(0);
492 }
493
494 /*@C
495 TSARKIMEXInitializePackage - This function initializes everything in the TSARKIMEX package. It is called
496 from TSInitializePackage().
497
498 Level: developer
499
500 .seealso: PetscInitialize()
501 @*/
TSARKIMEXInitializePackage(void)502 PetscErrorCode TSARKIMEXInitializePackage(void)
503 {
504 PetscErrorCode ierr;
505
506 PetscFunctionBegin;
507 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(0);
508 TSARKIMEXPackageInitialized = PETSC_TRUE;
509 ierr = TSARKIMEXRegisterAll();CHKERRQ(ierr);
510 ierr = PetscRegisterFinalize(TSARKIMEXFinalizePackage);CHKERRQ(ierr);
511 PetscFunctionReturn(0);
512 }
513
514 /*@C
515 TSARKIMEXFinalizePackage - This function destroys everything in the TSARKIMEX package. It is
516 called from PetscFinalize().
517
518 Level: developer
519
520 .seealso: PetscFinalize()
521 @*/
TSARKIMEXFinalizePackage(void)522 PetscErrorCode TSARKIMEXFinalizePackage(void)
523 {
524 PetscErrorCode ierr;
525
526 PetscFunctionBegin;
527 TSARKIMEXPackageInitialized = PETSC_FALSE;
528 ierr = TSARKIMEXRegisterDestroy();CHKERRQ(ierr);
529 PetscFunctionReturn(0);
530 }
531
532 /*@C
533 TSARKIMEXRegister - register an ARK IMEX scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
534
535 Not Collective, but the same schemes should be registered on all processes on which they will be used
536
537 Input Parameters:
538 + name - identifier for method
539 . order - approximation order of method
540 . s - number of stages, this is the dimension of the matrices below
541 . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major)
542 . bt - Butcher table for completing the stiff part of the step (dimension s; NULL to use the last row of At)
543 . ct - Abscissa of each stiff stage (dimension s, NULL to use row sums of At)
544 . A - Non-stiff stage coefficients (dimension s*s, row-major)
545 . b - Non-stiff step completion table (dimension s; NULL to use last row of At)
546 . c - Non-stiff abscissa (dimension s; NULL to use row sums of A)
547 . bembedt - Stiff part of completion table for embedded method (dimension s; NULL if not available)
548 . bembed - Non-stiff part of completion table for embedded method (dimension s; NULL to use bembedt if provided)
549 . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp
550 . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp)
551 - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; NULL to reuse binterpt)
552
553 Notes:
554 Several ARK IMEX methods are provided, this function is only needed to create new methods.
555
556 Level: advanced
557
558 .seealso: TSARKIMEX
559 @*/
TSARKIMEXRegister(TSARKIMEXType name,PetscInt order,PetscInt s,const PetscReal At[],const PetscReal bt[],const PetscReal ct[],const PetscReal A[],const PetscReal b[],const PetscReal c[],const PetscReal bembedt[],const PetscReal bembed[],PetscInt pinterp,const PetscReal binterpt[],const PetscReal binterp[])560 PetscErrorCode TSARKIMEXRegister(TSARKIMEXType name,PetscInt order,PetscInt s,
561 const PetscReal At[],const PetscReal bt[],const PetscReal ct[],
562 const PetscReal A[],const PetscReal b[],const PetscReal c[],
563 const PetscReal bembedt[],const PetscReal bembed[],
564 PetscInt pinterp,const PetscReal binterpt[],const PetscReal binterp[])
565 {
566 PetscErrorCode ierr;
567 ARKTableauLink link;
568 ARKTableau t;
569 PetscInt i,j;
570
571 PetscFunctionBegin;
572 ierr = TSARKIMEXInitializePackage();CHKERRQ(ierr);
573 ierr = PetscNew(&link);CHKERRQ(ierr);
574 t = &link->tab;
575 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
576 t->order = order;
577 t->s = s;
578 ierr = PetscMalloc6(s*s,&t->At,s,&t->bt,s,&t->ct,s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr);
579 ierr = PetscArraycpy(t->At,At,s*s);CHKERRQ(ierr);
580 ierr = PetscArraycpy(t->A,A,s*s);CHKERRQ(ierr);
581 if (bt) { ierr = PetscArraycpy(t->bt,bt,s);CHKERRQ(ierr); }
582 else for (i=0; i<s; i++) t->bt[i] = At[(s-1)*s+i];
583 if (b) { ierr = PetscArraycpy(t->b,b,s);CHKERRQ(ierr); }
584 else for (i=0; i<s; i++) t->b[i] = t->bt[i];
585 if (ct) { ierr = PetscArraycpy(t->ct,ct,s);CHKERRQ(ierr); }
586 else for (i=0; i<s; i++) for (j=0,t->ct[i]=0; j<s; j++) t->ct[i] += At[i*s+j];
587 if (c) { ierr = PetscArraycpy(t->c,c,s);CHKERRQ(ierr); }
588 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
589 t->stiffly_accurate = PETSC_TRUE;
590 for (i=0; i<s; i++) if (t->At[(s-1)*s+i] != t->bt[i]) t->stiffly_accurate = PETSC_FALSE;
591 t->explicit_first_stage = PETSC_TRUE;
592 for (i=0; i<s; i++) if (t->At[i] != 0.0) t->explicit_first_stage = PETSC_FALSE;
593 /*def of FSAL can be made more precise*/
594 t->FSAL_implicit = (PetscBool)(t->explicit_first_stage && t->stiffly_accurate);
595 if (bembedt) {
596 ierr = PetscMalloc2(s,&t->bembedt,s,&t->bembed);CHKERRQ(ierr);
597 ierr = PetscArraycpy(t->bembedt,bembedt,s);CHKERRQ(ierr);
598 ierr = PetscArraycpy(t->bembed,bembed ? bembed : bembedt,s);CHKERRQ(ierr);
599 }
600
601 t->pinterp = pinterp;
602 ierr = PetscMalloc2(s*pinterp,&t->binterpt,s*pinterp,&t->binterp);CHKERRQ(ierr);
603 ierr = PetscArraycpy(t->binterpt,binterpt,s*pinterp);CHKERRQ(ierr);
604 ierr = PetscArraycpy(t->binterp,binterp ? binterp : binterpt,s*pinterp);CHKERRQ(ierr);
605 link->next = ARKTableauList;
606 ARKTableauList = link;
607 PetscFunctionReturn(0);
608 }
609
610 /*
611 The step completion formula is
612
613 x1 = x0 - h bt^T YdotI + h b^T YdotRHS
614
615 This function can be called before or after ts->vec_sol has been updated.
616 Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order.
617 We can write
618
619 x1e = x0 - h bet^T YdotI + h be^T YdotRHS
620 = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS
621 = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS
622
623 so we can evaluate the method with different order even after the step has been optimistically completed.
624 */
TSEvaluateStep_ARKIMEX(TS ts,PetscInt order,Vec X,PetscBool * done)625 static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts,PetscInt order,Vec X,PetscBool *done)
626 {
627 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
628 ARKTableau tab = ark->tableau;
629 PetscScalar *w = ark->work;
630 PetscReal h;
631 PetscInt s = tab->s,j;
632 PetscErrorCode ierr;
633
634 PetscFunctionBegin;
635 switch (ark->status) {
636 case TS_STEP_INCOMPLETE:
637 case TS_STEP_PENDING:
638 h = ts->time_step; break;
639 case TS_STEP_COMPLETE:
640 h = ts->ptime - ts->ptime_prev; break;
641 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
642 }
643 if (order == tab->order) {
644 if (ark->status == TS_STEP_INCOMPLETE) {
645 if (!ark->imex && tab->stiffly_accurate) { /* Only the stiffly accurate implicit formula is used */
646 ierr = VecCopy(ark->Y[s-1],X);CHKERRQ(ierr);
647 } else { /* Use the standard completion formula (bt,b) */
648 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);
649 for (j=0; j<s; j++) w[j] = h*tab->bt[j];
650 ierr = VecMAXPY(X,s,w,ark->YdotI);CHKERRQ(ierr);
651 if (ark->imex) { /* Method is IMEX, complete the explicit formula */
652 for (j=0; j<s; j++) w[j] = h*tab->b[j];
653 ierr = VecMAXPY(X,s,w,ark->YdotRHS);CHKERRQ(ierr);
654 }
655 }
656 } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);}
657 if (done) *done = PETSC_TRUE;
658 PetscFunctionReturn(0);
659 } else if (order == tab->order-1) {
660 if (!tab->bembedt) goto unavailable;
661 if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */
662 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);
663 for (j=0; j<s; j++) w[j] = h*tab->bembedt[j];
664 ierr = VecMAXPY(X,s,w,ark->YdotI);CHKERRQ(ierr);
665 for (j=0; j<s; j++) w[j] = h*tab->bembed[j];
666 ierr = VecMAXPY(X,s,w,ark->YdotRHS);CHKERRQ(ierr);
667 } else { /* Rollback and re-complete using (bet-be,be-b) */
668 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);
669 for (j=0; j<s; j++) w[j] = h*(tab->bembedt[j] - tab->bt[j]);
670 ierr = VecMAXPY(X,tab->s,w,ark->YdotI);CHKERRQ(ierr);
671 for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]);
672 ierr = VecMAXPY(X,s,w,ark->YdotRHS);CHKERRQ(ierr);
673 }
674 if (done) *done = PETSC_TRUE;
675 PetscFunctionReturn(0);
676 }
677 unavailable:
678 if (done) *done = PETSC_FALSE;
679 else SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"ARKIMEX '%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);
680 PetscFunctionReturn(0);
681 }
682
TSRollBack_ARKIMEX(TS ts)683 static PetscErrorCode TSRollBack_ARKIMEX(TS ts)
684 {
685 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
686 ARKTableau tab = ark->tableau;
687 const PetscInt s = tab->s;
688 const PetscReal *bt = tab->bt,*b = tab->b;
689 PetscScalar *w = ark->work;
690 Vec *YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS;
691 PetscInt j;
692 PetscReal h;
693 PetscErrorCode ierr;
694
695 PetscFunctionBegin;
696 switch (ark->status) {
697 case TS_STEP_INCOMPLETE:
698 case TS_STEP_PENDING:
699 h = ts->time_step; break;
700 case TS_STEP_COMPLETE:
701 h = ts->ptime - ts->ptime_prev; break;
702 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
703 }
704 for (j=0; j<s; j++) w[j] = -h*bt[j];
705 ierr = VecMAXPY(ts->vec_sol,s,w,YdotI);CHKERRQ(ierr);
706 for (j=0; j<s; j++) w[j] = -h*b[j];
707 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr);
708 PetscFunctionReturn(0);
709 }
710
TSStep_ARKIMEX(TS ts)711 static PetscErrorCode TSStep_ARKIMEX(TS ts)
712 {
713 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
714 ARKTableau tab = ark->tableau;
715 const PetscInt s = tab->s;
716 const PetscReal *At = tab->At,*A = tab->A,*ct = tab->ct,*c = tab->c;
717 PetscScalar *w = ark->work;
718 Vec *Y = ark->Y,*YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS,Ydot = ark->Ydot,Ydot0 = ark->Ydot0,Z = ark->Z;
719 PetscBool extrapolate = ark->extrapolate;
720 TSAdapt adapt;
721 SNES snes;
722 PetscInt i,j,its,lits;
723 PetscInt rejections = 0;
724 PetscBool stageok,accept = PETSC_TRUE;
725 PetscReal next_time_step = ts->time_step;
726 PetscErrorCode ierr;
727
728 PetscFunctionBegin;
729 if (ark->extrapolate && !ark->Y_prev) {
730 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->Y_prev);CHKERRQ(ierr);
731 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->YdotI_prev);CHKERRQ(ierr);
732 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->YdotRHS_prev);CHKERRQ(ierr);
733 }
734
735 if (!ts->steprollback) {
736 if (ts->equation_type >= TS_EQ_IMPLICIT) { /* Save the initial slope for the next step */
737 ierr = VecCopy(YdotI[s-1],Ydot0);CHKERRQ(ierr);
738 }
739 if (ark->extrapolate && !ts->steprestart) { /* Save the Y, YdotI, YdotRHS for extrapolation initial guess */
740 for (i = 0; i<s; i++) {
741 ierr = VecCopy(Y[i],ark->Y_prev[i]);CHKERRQ(ierr);
742 ierr = VecCopy(YdotRHS[i],ark->YdotRHS_prev[i]);CHKERRQ(ierr);
743 ierr = VecCopy(YdotI[i],ark->YdotI_prev[i]);CHKERRQ(ierr);
744 }
745 }
746 }
747
748 if (ts->equation_type >= TS_EQ_IMPLICIT && tab->explicit_first_stage && ts->steprestart) {
749 TS ts_start;
750 ierr = TSClone(ts,&ts_start);CHKERRQ(ierr);
751 ierr = TSSetSolution(ts_start,ts->vec_sol);CHKERRQ(ierr);
752 ierr = TSSetTime(ts_start,ts->ptime);CHKERRQ(ierr);
753 ierr = TSSetMaxSteps(ts_start,ts->steps+1);CHKERRQ(ierr);
754 ierr = TSSetMaxTime(ts_start,ts->ptime+ts->time_step);CHKERRQ(ierr);
755 ierr = TSSetExactFinalTime(ts_start,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
756 ierr = TSSetTimeStep(ts_start,ts->time_step);CHKERRQ(ierr);
757 ierr = TSSetType(ts_start,TSARKIMEX);CHKERRQ(ierr);
758 ierr = TSARKIMEXSetFullyImplicit(ts_start,PETSC_TRUE);CHKERRQ(ierr);
759 ierr = TSARKIMEXSetType(ts_start,TSARKIMEX1BEE);CHKERRQ(ierr);
760
761 ierr = TSRestartStep(ts_start);CHKERRQ(ierr);
762 ierr = TSSolve(ts_start,ts->vec_sol);CHKERRQ(ierr);
763 ierr = TSGetTime(ts_start,&ts->ptime);CHKERRQ(ierr);
764 ierr = TSGetTimeStep(ts_start,&ts->time_step);CHKERRQ(ierr);
765
766 { /* Save the initial slope for the next step */
767 TS_ARKIMEX *ark_start = (TS_ARKIMEX*)ts_start->data;
768 ierr = VecCopy(ark_start->YdotI[ark_start->tableau->s-1],Ydot0);CHKERRQ(ierr);
769 }
770 ts->steps++;
771
772 /* Set the correct TS in SNES */
773 /* We'll try to bypass this by changing the method on the fly */
774 {
775 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
776 ierr = TSSetSNES(ts,snes);CHKERRQ(ierr);
777 }
778 ierr = TSDestroy(&ts_start);CHKERRQ(ierr);
779 }
780
781 ark->status = TS_STEP_INCOMPLETE;
782 while (!ts->reason && ark->status != TS_STEP_COMPLETE) {
783 PetscReal t = ts->ptime;
784 PetscReal h = ts->time_step;
785 for (i=0; i<s; i++) {
786 ark->stage_time = t + h*ct[i];
787 ierr = TSPreStage(ts,ark->stage_time);CHKERRQ(ierr);
788 if (At[i*s+i] == 0) { /* This stage is explicit */
789 if (i!=0 && ts->equation_type >= TS_EQ_IMPLICIT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Explicit stages other than the first one are not supported for implicit problems");
790 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr);
791 for (j=0; j<i; j++) w[j] = h*At[i*s+j];
792 ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr);
793 for (j=0; j<i; j++) w[j] = h*A[i*s+j];
794 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr);
795 } else {
796 ark->scoeff = 1./At[i*s+i];
797 /* Ydot = shift*(Y-Z) */
798 ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr);
799 for (j=0; j<i; j++) w[j] = h*At[i*s+j];
800 ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr);
801 for (j=0; j<i; j++) w[j] = h*A[i*s+j];
802 ierr = VecMAXPY(Z,i,w,YdotRHS);CHKERRQ(ierr);
803 if (extrapolate && !ts->steprestart) {
804 /* Initial guess extrapolated from previous time step stage values */
805 ierr = TSExtrapolate_ARKIMEX(ts,c[i],Y[i]);CHKERRQ(ierr);
806 } else {
807 /* Initial guess taken from last stage */
808 ierr = VecCopy(i>0 ? Y[i-1] : ts->vec_sol,Y[i]);CHKERRQ(ierr);
809 }
810 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
811 ierr = SNESSolve(snes,NULL,Y[i]);CHKERRQ(ierr);
812 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
813 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr);
814 ts->snes_its += its; ts->ksp_its += lits;
815 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr);
816 ierr = TSAdaptCheckStage(adapt,ts,ark->stage_time,Y[i],&stageok);CHKERRQ(ierr);
817 if (!stageok) {
818 /* We are likely rejecting the step because of solver or function domain problems so we should not attempt to
819 * use extrapolation to initialize the solves on the next attempt. */
820 extrapolate = PETSC_FALSE;
821 goto reject_step;
822 }
823 }
824 if (ts->equation_type >= TS_EQ_IMPLICIT) {
825 if (i==0 && tab->explicit_first_stage) {
826 if (!tab->stiffly_accurate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSARKIMEX %s is not stiffly accurate and therefore explicit-first stage methods cannot be used if the equation is implicit because the slope cannot be evaluated",ark->tableau->name);
827 ierr = VecCopy(Ydot0,YdotI[0]);CHKERRQ(ierr); /* YdotI = YdotI(tn-1) */
828 } else {
829 ierr = VecAXPBYPCZ(YdotI[i],-ark->scoeff/h,ark->scoeff/h,0,Z,Y[i]);CHKERRQ(ierr); /* YdotI = shift*(X-Z) */
830 }
831 } else {
832 if (i==0 && tab->explicit_first_stage) {
833 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr);
834 ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr);/* YdotI = -G(t,Y,0) */
835 ierr = VecScale(YdotI[i],-1.0);CHKERRQ(ierr);
836 } else {
837 ierr = VecAXPBYPCZ(YdotI[i],-ark->scoeff/h,ark->scoeff/h,0,Z,Y[i]);CHKERRQ(ierr); /* YdotI = shift*(X-Z) */
838 }
839 if (ark->imex) {
840 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr);
841 } else {
842 ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr);
843 }
844 }
845 ierr = TSPostStage(ts,ark->stage_time,i,Y);CHKERRQ(ierr);
846 }
847
848 ark->status = TS_STEP_INCOMPLETE;
849 ierr = TSEvaluateStep_ARKIMEX(ts,tab->order,ts->vec_sol,NULL);CHKERRQ(ierr);
850 ark->status = TS_STEP_PENDING;
851 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr);
852 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr);
853 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,(PetscReal)tab->s,PETSC_TRUE);CHKERRQ(ierr);
854 ierr = TSAdaptChoose(adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr);
855 ark->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
856 if (!accept) { /* Roll back the current step */
857 ierr = TSRollBack_ARKIMEX(ts);CHKERRQ(ierr);
858 ts->time_step = next_time_step;
859 goto reject_step;
860 }
861
862 ts->ptime += ts->time_step;
863 ts->time_step = next_time_step;
864 break;
865
866 reject_step:
867 ts->reject++; accept = PETSC_FALSE;
868 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
869 ts->reason = TS_DIVERGED_STEP_REJECTED;
870 ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr);
871 }
872 }
873 PetscFunctionReturn(0);
874 }
875
TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X)876 static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X)
877 {
878 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
879 PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j;
880 PetscReal h;
881 PetscReal tt,t;
882 PetscScalar *bt,*b;
883 const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp;
884 PetscErrorCode ierr;
885
886 PetscFunctionBegin;
887 if (!Bt || !B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name);
888 switch (ark->status) {
889 case TS_STEP_INCOMPLETE:
890 case TS_STEP_PENDING:
891 h = ts->time_step;
892 t = (itime - ts->ptime)/h;
893 break;
894 case TS_STEP_COMPLETE:
895 h = ts->ptime - ts->ptime_prev;
896 t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */
897 break;
898 default: SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Invalid TSStepStatus");
899 }
900 ierr = PetscMalloc2(s,&bt,s,&b);CHKERRQ(ierr);
901 for (i=0; i<s; i++) bt[i] = b[i] = 0;
902 for (j=0,tt=t; j<pinterp; j++,tt*=t) {
903 for (i=0; i<s; i++) {
904 bt[i] += h * Bt[i*pinterp+j] * tt;
905 b[i] += h * B[i*pinterp+j] * tt;
906 }
907 }
908 ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr);
909 ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr);
910 ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr);
911 ierr = PetscFree2(bt,b);CHKERRQ(ierr);
912 PetscFunctionReturn(0);
913 }
914
TSExtrapolate_ARKIMEX(TS ts,PetscReal c,Vec X)915 static PetscErrorCode TSExtrapolate_ARKIMEX(TS ts,PetscReal c,Vec X)
916 {
917 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
918 PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j;
919 PetscReal h,h_prev,t,tt;
920 PetscScalar *bt,*b;
921 const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp;
922 PetscErrorCode ierr;
923
924 PetscFunctionBegin;
925 if (!Bt || !B) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name);
926 ierr = PetscCalloc2(s,&bt,s,&b);CHKERRQ(ierr);
927 h = ts->time_step;
928 h_prev = ts->ptime - ts->ptime_prev;
929 t = 1 + h/h_prev*c;
930 for (j=0,tt=t; j<pinterp; j++,tt*=t) {
931 for (i=0; i<s; i++) {
932 bt[i] += h * Bt[i*pinterp+j] * tt;
933 b[i] += h * B[i*pinterp+j] * tt;
934 }
935 }
936 if (!ark->Y_prev) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Stages from previous step have not been stored");
937 ierr = VecCopy(ark->Y_prev[0],X);CHKERRQ(ierr);
938 ierr = VecMAXPY(X,s,bt,ark->YdotI_prev);CHKERRQ(ierr);
939 ierr = VecMAXPY(X,s,b,ark->YdotRHS_prev);CHKERRQ(ierr);
940 ierr = PetscFree2(bt,b);CHKERRQ(ierr);
941 PetscFunctionReturn(0);
942 }
943
944 /*------------------------------------------------------------*/
945
TSARKIMEXTableauReset(TS ts)946 static PetscErrorCode TSARKIMEXTableauReset(TS ts)
947 {
948 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
949 ARKTableau tab = ark->tableau;
950 PetscErrorCode ierr;
951
952 PetscFunctionBegin;
953 if (!tab) PetscFunctionReturn(0);
954 ierr = PetscFree(ark->work);CHKERRQ(ierr);
955 ierr = VecDestroyVecs(tab->s,&ark->Y);CHKERRQ(ierr);
956 ierr = VecDestroyVecs(tab->s,&ark->YdotI);CHKERRQ(ierr);
957 ierr = VecDestroyVecs(tab->s,&ark->YdotRHS);CHKERRQ(ierr);
958 ierr = VecDestroyVecs(tab->s,&ark->Y_prev);CHKERRQ(ierr);
959 ierr = VecDestroyVecs(tab->s,&ark->YdotI_prev);CHKERRQ(ierr);
960 ierr = VecDestroyVecs(tab->s,&ark->YdotRHS_prev);CHKERRQ(ierr);
961 PetscFunctionReturn(0);
962 }
963
TSReset_ARKIMEX(TS ts)964 static PetscErrorCode TSReset_ARKIMEX(TS ts)
965 {
966 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
967 PetscErrorCode ierr;
968
969 PetscFunctionBegin;
970 ierr = TSARKIMEXTableauReset(ts);CHKERRQ(ierr);
971 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr);
972 ierr = VecDestroy(&ark->Ydot0);CHKERRQ(ierr);
973 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr);
974 PetscFunctionReturn(0);
975 }
976
TSARKIMEXGetVecs(TS ts,DM dm,Vec * Z,Vec * Ydot)977 static PetscErrorCode TSARKIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot)
978 {
979 TS_ARKIMEX *ax = (TS_ARKIMEX*)ts->data;
980 PetscErrorCode ierr;
981
982 PetscFunctionBegin;
983 if (Z) {
984 if (dm && dm != ts->dm) {
985 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr);
986 } else *Z = ax->Z;
987 }
988 if (Ydot) {
989 if (dm && dm != ts->dm) {
990 ierr = DMGetNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr);
991 } else *Ydot = ax->Ydot;
992 }
993 PetscFunctionReturn(0);
994 }
995
996
TSARKIMEXRestoreVecs(TS ts,DM dm,Vec * Z,Vec * Ydot)997 static PetscErrorCode TSARKIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot)
998 {
999 PetscErrorCode ierr;
1000
1001 PetscFunctionBegin;
1002 if (Z) {
1003 if (dm && dm != ts->dm) {
1004 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Z",Z);CHKERRQ(ierr);
1005 }
1006 }
1007 if (Ydot) {
1008 if (dm && dm != ts->dm) {
1009 ierr = DMRestoreNamedGlobalVector(dm,"TSARKIMEX_Ydot",Ydot);CHKERRQ(ierr);
1010 }
1011 }
1012 PetscFunctionReturn(0);
1013 }
1014
1015 /*
1016 This defines the nonlinear equation that is to be solved with SNES
1017 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
1018 */
SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts)1019 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts)
1020 {
1021 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1022 DM dm,dmsave;
1023 Vec Z,Ydot;
1024 PetscReal shift = ark->scoeff / ts->time_step;
1025 PetscErrorCode ierr;
1026
1027 PetscFunctionBegin;
1028 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
1029 ierr = TSARKIMEXGetVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr);
1030 ierr = VecAXPBYPCZ(Ydot,-shift,shift,0,Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */
1031 dmsave = ts->dm;
1032 ts->dm = dm;
1033
1034 ierr = TSComputeIFunction(ts,ark->stage_time,X,Ydot,F,ark->imex);CHKERRQ(ierr);
1035
1036 ts->dm = dmsave;
1037 ierr = TSARKIMEXRestoreVecs(ts,dm,&Z,&Ydot);CHKERRQ(ierr);
1038 PetscFunctionReturn(0);
1039 }
1040
SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)1041 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat A,Mat B,TS ts)
1042 {
1043 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1044 DM dm,dmsave;
1045 Vec Ydot;
1046 PetscReal shift = ark->scoeff / ts->time_step;
1047 PetscErrorCode ierr;
1048
1049 PetscFunctionBegin;
1050 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
1051 ierr = TSARKIMEXGetVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr);
1052 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */
1053 dmsave = ts->dm;
1054 ts->dm = dm;
1055
1056 ierr = TSComputeIJacobian(ts,ark->stage_time,X,Ydot,shift,A,B,ark->imex);CHKERRQ(ierr);
1057
1058 ts->dm = dmsave;
1059 ierr = TSARKIMEXRestoreVecs(ts,dm,NULL,&Ydot);CHKERRQ(ierr);
1060 PetscFunctionReturn(0);
1061 }
1062
DMCoarsenHook_TSARKIMEX(DM fine,DM coarse,void * ctx)1063 static PetscErrorCode DMCoarsenHook_TSARKIMEX(DM fine,DM coarse,void *ctx)
1064 {
1065 PetscFunctionBegin;
1066 PetscFunctionReturn(0);
1067 }
1068
DMRestrictHook_TSARKIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void * ctx)1069 static PetscErrorCode DMRestrictHook_TSARKIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
1070 {
1071 TS ts = (TS)ctx;
1072 PetscErrorCode ierr;
1073 Vec Z,Z_c;
1074
1075 PetscFunctionBegin;
1076 ierr = TSARKIMEXGetVecs(ts,fine,&Z,NULL);CHKERRQ(ierr);
1077 ierr = TSARKIMEXGetVecs(ts,coarse,&Z_c,NULL);CHKERRQ(ierr);
1078 ierr = MatRestrict(restrct,Z,Z_c);CHKERRQ(ierr);
1079 ierr = VecPointwiseMult(Z_c,rscale,Z_c);CHKERRQ(ierr);
1080 ierr = TSARKIMEXRestoreVecs(ts,fine,&Z,NULL);CHKERRQ(ierr);
1081 ierr = TSARKIMEXRestoreVecs(ts,coarse,&Z_c,NULL);CHKERRQ(ierr);
1082 PetscFunctionReturn(0);
1083 }
1084
1085
DMSubDomainHook_TSARKIMEX(DM dm,DM subdm,void * ctx)1086 static PetscErrorCode DMSubDomainHook_TSARKIMEX(DM dm,DM subdm,void *ctx)
1087 {
1088 PetscFunctionBegin;
1089 PetscFunctionReturn(0);
1090 }
1091
DMSubDomainRestrictHook_TSARKIMEX(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void * ctx)1092 static PetscErrorCode DMSubDomainRestrictHook_TSARKIMEX(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
1093 {
1094 TS ts = (TS)ctx;
1095 PetscErrorCode ierr;
1096 Vec Z,Z_c;
1097
1098 PetscFunctionBegin;
1099 ierr = TSARKIMEXGetVecs(ts,dm,&Z,NULL);CHKERRQ(ierr);
1100 ierr = TSARKIMEXGetVecs(ts,subdm,&Z_c,NULL);CHKERRQ(ierr);
1101
1102 ierr = VecScatterBegin(gscat,Z,Z_c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1103 ierr = VecScatterEnd(gscat,Z,Z_c,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1104
1105 ierr = TSARKIMEXRestoreVecs(ts,dm,&Z,NULL);CHKERRQ(ierr);
1106 ierr = TSARKIMEXRestoreVecs(ts,subdm,&Z_c,NULL);CHKERRQ(ierr);
1107 PetscFunctionReturn(0);
1108 }
1109
TSARKIMEXTableauSetUp(TS ts)1110 static PetscErrorCode TSARKIMEXTableauSetUp(TS ts)
1111 {
1112 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1113 ARKTableau tab = ark->tableau;
1114 PetscErrorCode ierr;
1115
1116 PetscFunctionBegin;
1117 ierr = PetscMalloc1(tab->s,&ark->work);CHKERRQ(ierr);
1118 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->Y);CHKERRQ(ierr);
1119 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->YdotI);CHKERRQ(ierr);
1120 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->YdotRHS);CHKERRQ(ierr);
1121 if (ark->extrapolate) {
1122 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->Y_prev);CHKERRQ(ierr);
1123 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->YdotI_prev);CHKERRQ(ierr);
1124 ierr = VecDuplicateVecs(ts->vec_sol,tab->s,&ark->YdotRHS_prev);CHKERRQ(ierr);
1125 }
1126 PetscFunctionReturn(0);
1127 }
1128
TSSetUp_ARKIMEX(TS ts)1129 static PetscErrorCode TSSetUp_ARKIMEX(TS ts)
1130 {
1131 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1132 PetscErrorCode ierr;
1133 DM dm;
1134 SNES snes;
1135
1136 PetscFunctionBegin;
1137 ierr = TSARKIMEXTableauSetUp(ts);CHKERRQ(ierr);
1138 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr);
1139 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot0);CHKERRQ(ierr);
1140 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr);
1141 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
1142 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSARKIMEX,DMRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr);
1143 ierr = DMSubDomainHookAdd(dm,DMSubDomainHook_TSARKIMEX,DMSubDomainRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr);
1144 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
1145 PetscFunctionReturn(0);
1146 }
1147 /*------------------------------------------------------------*/
1148
TSSetFromOptions_ARKIMEX(PetscOptionItems * PetscOptionsObject,TS ts)1149 static PetscErrorCode TSSetFromOptions_ARKIMEX(PetscOptionItems *PetscOptionsObject,TS ts)
1150 {
1151 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1152 PetscErrorCode ierr;
1153
1154 PetscFunctionBegin;
1155 ierr = PetscOptionsHead(PetscOptionsObject,"ARKIMEX ODE solver options");CHKERRQ(ierr);
1156 {
1157 ARKTableauLink link;
1158 PetscInt count,choice;
1159 PetscBool flg;
1160 const char **namelist;
1161 for (link=ARKTableauList,count=0; link; link=link->next,count++) ;
1162 ierr = PetscMalloc1(count,(char***)&namelist);CHKERRQ(ierr);
1163 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
1164 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,ark->tableau->name,&choice,&flg);CHKERRQ(ierr);
1165 if (flg) {ierr = TSARKIMEXSetType(ts,namelist[choice]);CHKERRQ(ierr);}
1166 ierr = PetscFree(namelist);CHKERRQ(ierr);
1167
1168 flg = (PetscBool) !ark->imex;
1169 ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,NULL);CHKERRQ(ierr);
1170 ark->imex = (PetscBool) !flg;
1171 ierr = PetscOptionsBool("-ts_arkimex_initial_guess_extrapolate","Extrapolate the initial guess for the stage solution from stage values of the previous time step","",ark->extrapolate,&ark->extrapolate,NULL);CHKERRQ(ierr);
1172 }
1173 ierr = PetscOptionsTail();CHKERRQ(ierr);
1174 PetscFunctionReturn(0);
1175 }
1176
TSView_ARKIMEX(TS ts,PetscViewer viewer)1177 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer)
1178 {
1179 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1180 PetscBool iascii;
1181 PetscErrorCode ierr;
1182
1183 PetscFunctionBegin;
1184 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1185 if (iascii) {
1186 ARKTableau tab = ark->tableau;
1187 TSARKIMEXType arktype;
1188 char buf[512];
1189 PetscBool flg;
1190
1191 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr);
1192 ierr = TSARKIMEXGetFullyImplicit(ts,&flg);CHKERRQ(ierr);
1193 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr);
1194 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr);
1195 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr);
1196 ierr = PetscFormatRealArray(buf,sizeof(buf),"% 8.6f",tab->s,tab->c);CHKERRQ(ierr);
1197 ierr = PetscViewerASCIIPrintf(viewer,"Fully implicit: %s\n",flg ? "yes" : "no");CHKERRQ(ierr);
1198 ierr = PetscViewerASCIIPrintf(viewer,"Stiffly accurate: %s\n",tab->stiffly_accurate ? "yes" : "no");CHKERRQ(ierr);
1199 ierr = PetscViewerASCIIPrintf(viewer,"Explicit first stage: %s\n",tab->explicit_first_stage ? "yes" : "no");CHKERRQ(ierr);
1200 ierr = PetscViewerASCIIPrintf(viewer,"FSAL property: %s\n",tab->FSAL_implicit ? "yes" : "no");CHKERRQ(ierr);
1201 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr);
1202 }
1203 PetscFunctionReturn(0);
1204 }
1205
TSLoad_ARKIMEX(TS ts,PetscViewer viewer)1206 static PetscErrorCode TSLoad_ARKIMEX(TS ts,PetscViewer viewer)
1207 {
1208 PetscErrorCode ierr;
1209 SNES snes;
1210 TSAdapt adapt;
1211
1212 PetscFunctionBegin;
1213 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr);
1214 ierr = TSAdaptLoad(adapt,viewer);CHKERRQ(ierr);
1215 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
1216 ierr = SNESLoad(snes,viewer);CHKERRQ(ierr);
1217 /* function and Jacobian context for SNES when used with TS is always ts object */
1218 ierr = SNESSetFunction(snes,NULL,NULL,ts);CHKERRQ(ierr);
1219 ierr = SNESSetJacobian(snes,NULL,NULL,NULL,ts);CHKERRQ(ierr);
1220 PetscFunctionReturn(0);
1221 }
1222
1223 /*@C
1224 TSARKIMEXSetType - Set the type of ARK IMEX scheme
1225
1226 Logically collective
1227
1228 Input Parameter:
1229 + ts - timestepping context
1230 - arktype - type of ARK-IMEX scheme
1231
1232 Options Database:
1233 . -ts_arkimex_type <1bee,a2,l2,ars122,2c,2d,2e,prssp2,3,bpr3,ars443,4,5>
1234
1235 Level: intermediate
1236
1237 .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEXType, TSARKIMEX1BEE, TSARKIMEXA2, TSARKIMEXL2, TSARKIMEXARS122, TSARKIMEX2C, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2,
1238 TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5
1239 @*/
TSARKIMEXSetType(TS ts,TSARKIMEXType arktype)1240 PetscErrorCode TSARKIMEXSetType(TS ts,TSARKIMEXType arktype)
1241 {
1242 PetscErrorCode ierr;
1243
1244 PetscFunctionBegin;
1245 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1246 PetscValidCharPointer(arktype,2);
1247 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,TSARKIMEXType),(ts,arktype));CHKERRQ(ierr);
1248 PetscFunctionReturn(0);
1249 }
1250
1251 /*@C
1252 TSARKIMEXGetType - Get the type of ARK IMEX scheme
1253
1254 Logically collective
1255
1256 Input Parameter:
1257 . ts - timestepping context
1258
1259 Output Parameter:
1260 . arktype - type of ARK-IMEX scheme
1261
1262 Level: intermediate
1263
1264 .seealso: TSARKIMEXGetType()
1265 @*/
TSARKIMEXGetType(TS ts,TSARKIMEXType * arktype)1266 PetscErrorCode TSARKIMEXGetType(TS ts,TSARKIMEXType *arktype)
1267 {
1268 PetscErrorCode ierr;
1269
1270 PetscFunctionBegin;
1271 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1272 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr);
1273 PetscFunctionReturn(0);
1274 }
1275
1276 /*@
1277 TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly
1278
1279 Logically collective
1280
1281 Input Parameter:
1282 + ts - timestepping context
1283 - flg - PETSC_TRUE for fully implicit
1284
1285 Level: intermediate
1286
1287 .seealso: TSARKIMEXGetType(), TSARKIMEXGetFullyImplicit()
1288 @*/
TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg)1289 PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg)
1290 {
1291 PetscErrorCode ierr;
1292
1293 PetscFunctionBegin;
1294 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1295 PetscValidLogicalCollectiveBool(ts,flg,2);
1296 ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr);
1297 PetscFunctionReturn(0);
1298 }
1299
1300 /*@
1301 TSARKIMEXGetFullyImplicit - Inquires if both parts of the equation are solved implicitly
1302
1303 Logically collective
1304
1305 Input Parameter:
1306 . ts - timestepping context
1307
1308 Output Parameter:
1309 . flg - PETSC_TRUE for fully implicit
1310
1311 Level: intermediate
1312
1313 .seealso: TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit()
1314 @*/
TSARKIMEXGetFullyImplicit(TS ts,PetscBool * flg)1315 PetscErrorCode TSARKIMEXGetFullyImplicit(TS ts,PetscBool *flg)
1316 {
1317 PetscErrorCode ierr;
1318
1319 PetscFunctionBegin;
1320 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
1321 PetscValidPointer(flg,2);
1322 ierr = PetscUseMethod(ts,"TSARKIMEXGetFullyImplicit_C",(TS,PetscBool*),(ts,flg));CHKERRQ(ierr);
1323 PetscFunctionReturn(0);
1324 }
1325
TSARKIMEXGetType_ARKIMEX(TS ts,TSARKIMEXType * arktype)1326 static PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,TSARKIMEXType *arktype)
1327 {
1328 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1329
1330 PetscFunctionBegin;
1331 *arktype = ark->tableau->name;
1332 PetscFunctionReturn(0);
1333 }
TSARKIMEXSetType_ARKIMEX(TS ts,TSARKIMEXType arktype)1334 static PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,TSARKIMEXType arktype)
1335 {
1336 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1337 PetscErrorCode ierr;
1338 PetscBool match;
1339 ARKTableauLink link;
1340
1341 PetscFunctionBegin;
1342 if (ark->tableau) {
1343 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr);
1344 if (match) PetscFunctionReturn(0);
1345 }
1346 for (link = ARKTableauList; link; link=link->next) {
1347 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr);
1348 if (match) {
1349 if (ts->setupcalled) {ierr = TSARKIMEXTableauReset(ts);CHKERRQ(ierr);}
1350 ark->tableau = &link->tab;
1351 if (ts->setupcalled) {ierr = TSARKIMEXTableauSetUp(ts);CHKERRQ(ierr);}
1352 ts->default_adapt_type = ark->tableau->bembed ? TSADAPTBASIC : TSADAPTNONE;
1353 PetscFunctionReturn(0);
1354 }
1355 }
1356 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype);
1357 PetscFunctionReturn(0);
1358 }
1359
TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg)1360 static PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg)
1361 {
1362 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1363
1364 PetscFunctionBegin;
1365 ark->imex = (PetscBool)!flg;
1366 PetscFunctionReturn(0);
1367 }
1368
TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts,PetscBool * flg)1369 static PetscErrorCode TSARKIMEXGetFullyImplicit_ARKIMEX(TS ts,PetscBool *flg)
1370 {
1371 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
1372
1373 PetscFunctionBegin;
1374 *flg = (PetscBool)!ark->imex;
1375 PetscFunctionReturn(0);
1376 }
1377
TSDestroy_ARKIMEX(TS ts)1378 static PetscErrorCode TSDestroy_ARKIMEX(TS ts)
1379 {
1380 PetscErrorCode ierr;
1381
1382 PetscFunctionBegin;
1383 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr);
1384 if (ts->dm) {
1385 ierr = DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSARKIMEX,DMRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr);
1386 ierr = DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSARKIMEX,DMSubDomainRestrictHook_TSARKIMEX,ts);CHKERRQ(ierr);
1387 }
1388 ierr = PetscFree(ts->data);CHKERRQ(ierr);
1389 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXGetType_C",NULL);CHKERRQ(ierr);
1390 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetType_C",NULL);CHKERRQ(ierr);
1391 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C",NULL);CHKERRQ(ierr);
1392 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C",NULL);CHKERRQ(ierr);
1393 PetscFunctionReturn(0);
1394 }
1395
1396 /* ------------------------------------------------------------ */
1397 /*MC
1398 TSARKIMEX - ODE and DAE solver using Additive Runge-Kutta IMEX schemes
1399
1400 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
1401 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
1402 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction().
1403
1404 Notes:
1405 The default is TSARKIMEX3, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type
1406
1407 If the equation is implicit or a DAE, then TSSetEquationType() needs to be set accordingly. Refer to the manual for further information.
1408
1409 Methods with an explicit stage can only be used with ODE in which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).
1410
1411 Consider trying TSROSW if the stiff part is linear or weakly nonlinear.
1412
1413 Level: beginner
1414
1415 .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEXGetFullyImplicit(),
1416 TSARKIMEX1BEE, TSARKIMEX2C, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEX3, TSARKIMEXL2, TSARKIMEXA2, TSARKIMEXARS122,
1417 TSARKIMEX4, TSARKIMEX5, TSARKIMEXPRSSP2, TSARKIMEXARS443, TSARKIMEXBPR3, TSARKIMEXType, TSARKIMEXRegister()
1418
1419 M*/
TSCreate_ARKIMEX(TS ts)1420 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS ts)
1421 {
1422 TS_ARKIMEX *th;
1423 PetscErrorCode ierr;
1424
1425 PetscFunctionBegin;
1426 ierr = TSARKIMEXInitializePackage();CHKERRQ(ierr);
1427
1428 ts->ops->reset = TSReset_ARKIMEX;
1429 ts->ops->destroy = TSDestroy_ARKIMEX;
1430 ts->ops->view = TSView_ARKIMEX;
1431 ts->ops->load = TSLoad_ARKIMEX;
1432 ts->ops->setup = TSSetUp_ARKIMEX;
1433 ts->ops->step = TSStep_ARKIMEX;
1434 ts->ops->interpolate = TSInterpolate_ARKIMEX;
1435 ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX;
1436 ts->ops->rollback = TSRollBack_ARKIMEX;
1437 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX;
1438 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX;
1439 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX;
1440
1441 ts->usessnes = PETSC_TRUE;
1442
1443 ierr = PetscNewLog(ts,&th);CHKERRQ(ierr);
1444 ts->data = (void*)th;
1445 th->imex = PETSC_TRUE;
1446
1447 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXGetType_C",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr);
1448 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetType_C",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr);
1449 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr);
1450 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSARKIMEXGetFullyImplicit_C",TSARKIMEXGetFullyImplicit_ARKIMEX);CHKERRQ(ierr);
1451
1452 ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);
1453 PetscFunctionReturn(0);
1454 }
1455