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