xref: /petsc/src/ts/impls/arkimex/arkimex.c (revision 6cf0794eb0b2857d96a70b98984d398bcf9034cb) 
18a381b04SJed Brown /*
28a381b04SJed Brown   Code for timestepping with additive Runge-Kutta IMEX method
38a381b04SJed Brown 
48a381b04SJed Brown   Notes:
58a381b04SJed Brown   The general system is written as
68a381b04SJed Brown 
78a381b04SJed Brown   F(t,X,Xdot) = G(t,X)
88a381b04SJed Brown 
98a381b04SJed Brown   where F represents the stiff part of the physics and G represents the non-stiff part.
108a381b04SJed Brown 
118a381b04SJed Brown */
128a381b04SJed Brown #include <private/tsimpl.h>                /*I   "petscts.h"   I*/
138a381b04SJed Brown 
144f385281SJed Brown static const TSARKIMEXType TSARKIMEXDefault = TSARKIMEX2E;
158a381b04SJed Brown static PetscBool TSARKIMEXRegisterAllCalled;
168a381b04SJed Brown static PetscBool TSARKIMEXPackageInitialized;
178a381b04SJed Brown 
188a381b04SJed Brown typedef struct _ARKTableau *ARKTableau;
198a381b04SJed Brown struct _ARKTableau {
208a381b04SJed Brown   char *name;
214f385281SJed Brown   PetscInt order;               /* Classical approximation order of the method */
224f385281SJed Brown   PetscInt s;                   /* Number of stages */
234f385281SJed Brown   PetscInt pinterp;             /* Interpolation order */
244f385281SJed Brown   PetscReal *At,*bt,*ct;        /* Stiff tableau */
258a381b04SJed Brown   PetscReal *A,*b,*c;           /* Non-stiff tableau */
26cd652676SJed Brown   PetscReal *binterpt,*binterp; /* Dense output formula */
278a381b04SJed Brown };
288a381b04SJed Brown typedef struct _ARKTableauLink *ARKTableauLink;
298a381b04SJed Brown struct _ARKTableauLink {
308a381b04SJed Brown   struct _ARKTableau tab;
318a381b04SJed Brown   ARKTableauLink next;
328a381b04SJed Brown };
338a381b04SJed Brown static ARKTableauLink ARKTableauList;
348a381b04SJed Brown 
358a381b04SJed Brown typedef struct {
368a381b04SJed Brown   ARKTableau  tableau;
378a381b04SJed Brown   Vec         *Y;               /* States computed during the step */
388a381b04SJed Brown   Vec         *YdotI;           /* Time derivatives for the stiff part */
398a381b04SJed Brown   Vec         *YdotRHS;         /* Function evaluations for the non-stiff part */
408a381b04SJed Brown   Vec         Ydot;             /* Work vector holding Ydot during residual evaluation */
418a381b04SJed Brown   Vec         Work;             /* Generic work vector */
428a381b04SJed Brown   Vec         Z;                /* Ydot = shift(Y-Z) */
438a381b04SJed Brown   PetscScalar *work;            /* Scalar work */
448a381b04SJed Brown   PetscReal   shift;
458a381b04SJed Brown   PetscReal   stage_time;
464cc180ffSJed Brown   PetscBool   imex;
478a381b04SJed Brown } TS_ARKIMEX;
488a381b04SJed Brown 
4964f491ddSJed Brown /*MC
5064f491ddSJed Brown      TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part.
5164f491ddSJed Brown 
5264f491ddSJed Brown      This method has one explicit stage and two implicit stages.
5364f491ddSJed Brown 
54b330ce4dSSatish Balay      Level: advanced
55b330ce4dSSatish Balay 
5664f491ddSJed Brown .seealso: TSARKIMEX
5764f491ddSJed Brown M*/
5864f491ddSJed Brown /*MC
5964f491ddSJed Brown      TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part.
6064f491ddSJed Brown 
6164f491ddSJed Brown      This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu.
6264f491ddSJed Brown 
63b330ce4dSSatish Balay      Level: advanced
64b330ce4dSSatish Balay 
6564f491ddSJed Brown .seealso: TSARKIMEX
6664f491ddSJed Brown M*/
6764f491ddSJed Brown /*MC
68*6cf0794eSJed Brown      TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme.
69*6cf0794eSJed Brown 
70*6cf0794eSJed Brown      This method has three implicit stages.
71*6cf0794eSJed Brown 
72*6cf0794eSJed Brown      References:
73*6cf0794eSJed Brown      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
74*6cf0794eSJed Brown 
75*6cf0794eSJed Brown      This method is referred to as SSP2-(3,3,2) in http://arxiv.org/abs/1110.4375
76*6cf0794eSJed Brown 
77*6cf0794eSJed Brown      Level: advanced
78*6cf0794eSJed Brown 
79*6cf0794eSJed Brown .seealso: TSARKIMEX
80*6cf0794eSJed Brown M*/
81*6cf0794eSJed Brown /*MC
8264f491ddSJed Brown      TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part.
8364f491ddSJed Brown 
8464f491ddSJed Brown      This method has one explicit stage and three implicit stages.
8564f491ddSJed Brown 
8664f491ddSJed Brown      References:
8764f491ddSJed Brown      Kennedy and Carpenter 2003.
8864f491ddSJed Brown 
89b330ce4dSSatish Balay      Level: advanced
90b330ce4dSSatish Balay 
9164f491ddSJed Brown .seealso: TSARKIMEX
9264f491ddSJed Brown M*/
9364f491ddSJed Brown /*MC
94*6cf0794eSJed Brown      TSARKIMEXARS443 - Third order ARK IMEX scheme.
95*6cf0794eSJed Brown 
96*6cf0794eSJed Brown      This method has one explicit stage and four implicit stages.
97*6cf0794eSJed Brown 
98*6cf0794eSJed Brown      References:
99*6cf0794eSJed Brown      U. Ascher, S. Ruuth, R. J. Spitheri, Implicit-explicit Runge-Kutta methods for time dependent Partial Differential Equations. Appl. Numer. Math. 25, (1997), pp. 151–167.
100*6cf0794eSJed Brown 
101*6cf0794eSJed Brown      This method is referred to as ARS(4,4,3) in http://arxiv.org/abs/1110.4375
102*6cf0794eSJed Brown 
103*6cf0794eSJed Brown      Level: advanced
104*6cf0794eSJed Brown 
105*6cf0794eSJed Brown .seealso: TSARKIMEX
106*6cf0794eSJed Brown M*/
107*6cf0794eSJed Brown /*MC
108*6cf0794eSJed Brown      TSARKIMEXBPR3 - Third order ARK IMEX scheme.
109*6cf0794eSJed Brown 
110*6cf0794eSJed Brown      This method has one explicit stage and four implicit stages.
111*6cf0794eSJed Brown 
112*6cf0794eSJed Brown      References:
113*6cf0794eSJed Brown      This method is referred to as ARK3 in http://arxiv.org/abs/1110.4375
114*6cf0794eSJed Brown 
115*6cf0794eSJed Brown      Level: advanced
116*6cf0794eSJed Brown 
117*6cf0794eSJed Brown .seealso: TSARKIMEX
118*6cf0794eSJed Brown M*/
119*6cf0794eSJed Brown /*MC
12064f491ddSJed Brown      TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part.
12164f491ddSJed Brown 
12264f491ddSJed Brown      This method has one explicit stage and four implicit stages.
12364f491ddSJed Brown 
12464f491ddSJed Brown      References:
12564f491ddSJed Brown      Kennedy and Carpenter 2003.
12664f491ddSJed Brown 
127b330ce4dSSatish Balay      Level: advanced
128b330ce4dSSatish Balay 
12964f491ddSJed Brown .seealso: TSARKIMEX
13064f491ddSJed Brown M*/
13164f491ddSJed Brown /*MC
13264f491ddSJed Brown      TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part.
13364f491ddSJed Brown 
13464f491ddSJed Brown      This method has one explicit stage and five implicit stages.
13564f491ddSJed Brown 
13664f491ddSJed Brown      References:
13764f491ddSJed Brown      Kennedy and Carpenter 2003.
13864f491ddSJed Brown 
139b330ce4dSSatish Balay      Level: advanced
140b330ce4dSSatish Balay 
14164f491ddSJed Brown .seealso: TSARKIMEX
14264f491ddSJed Brown M*/
14364f491ddSJed Brown 
1448a381b04SJed Brown #undef __FUNCT__
1458a381b04SJed Brown #define __FUNCT__ "TSARKIMEXRegisterAll"
1468a381b04SJed Brown /*@C
1478a381b04SJed Brown   TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSARKIMEX
1488a381b04SJed Brown 
149fca742c7SJed Brown   Not Collective, but should be called by all processes which will need the schemes to be registered
1508a381b04SJed Brown 
1518a381b04SJed Brown   Level: advanced
1528a381b04SJed Brown 
1538a381b04SJed Brown .keywords: TS, TSARKIMEX, register, all
1548a381b04SJed Brown 
1558a381b04SJed Brown .seealso:  TSARKIMEXRegisterDestroy()
1568a381b04SJed Brown @*/
1578a381b04SJed Brown PetscErrorCode TSARKIMEXRegisterAll(void)
1588a381b04SJed Brown {
1598a381b04SJed Brown   PetscErrorCode ierr;
1608a381b04SJed Brown 
1618a381b04SJed Brown   PetscFunctionBegin;
1628a381b04SJed Brown   if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(0);
1638a381b04SJed Brown   TSARKIMEXRegisterAllCalled = PETSC_TRUE;
1648a381b04SJed Brown 
1658a381b04SJed Brown   {
1668a381b04SJed Brown     const PetscReal
1678a381b04SJed Brown       A[3][3] = {{0,0,0},
1688a381b04SJed Brown                  {0.41421356237309504880,0,0},
1698a381b04SJed Brown                  {0.75,0.25,0}},
1708a381b04SJed Brown       At[3][3] = {{0,0,0},
1718a381b04SJed Brown                   {0.12132034355964257320,0.29289321881345247560,0},
17206db7b1cSJed Brown                   {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}},
17306db7b1cSJed Brown       binterpt[3][2] = {{1,-0.5},{0,0},{0,0.5}};
17406db7b1cSJed Brown     ierr = TSARKIMEXRegister(TSARKIMEX2D,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr);
1758a381b04SJed Brown   }
17606db7b1cSJed Brown   {                             /* Optimal for linear implicit part */
17703403c7fSJed Brown     const PetscReal s2 = PetscSqrtReal((PetscReal)2.0),
178a3a57f36SJed Brown       A[3][3] = {{0,0,0},
179a3a57f36SJed Brown                  {2-s2,0,0},
180a3a57f36SJed Brown                  {(3-2*s2)/6,(3+2*s2)/6,0}},
181a3a57f36SJed Brown       At[3][3] = {{0,0,0},
182a3a57f36SJed Brown                   {1-1/s2,1-1/s2,0},
183cd652676SJed Brown                   {1/(2*s2),1/(2*s2),1-1/s2}},
184cd652676SJed Brown       binterpt[3][2] = {{1,-0.5},{0,0},{0,0.5}};
185cd652676SJed Brown     ierr = TSARKIMEXRegister(TSARKIMEX2E,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr);
186a3a57f36SJed Brown   }
187*6cf0794eSJed Brown   {                             /* Optimal for linear implicit part */
188*6cf0794eSJed Brown     const PetscReal
189*6cf0794eSJed Brown       A[3][3] = {{0,0,0},
190*6cf0794eSJed Brown                  {0.5,0,0},
191*6cf0794eSJed Brown                  {0.5,0.5,0}},
192*6cf0794eSJed Brown       At[3][3] = {{0.25,0,0},
193*6cf0794eSJed Brown                   {0,0.25,0},
194*6cf0794eSJed Brown                   {1./3,1./3,1./3}};
195*6cf0794eSJed Brown     ierr = TSARKIMEXRegister(TSARKIMEXPRSSP2,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,0,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
196*6cf0794eSJed Brown   }
197a3a57f36SJed Brown   {
198a3a57f36SJed Brown     const PetscReal
199a3a57f36SJed Brown       A[4][4] = {{0,0,0,0},
2004040e9f2SJed Brown                  {1767732205903./2027836641118.,0,0,0},
2014040e9f2SJed Brown                  {5535828885825./10492691773637.,788022342437./10882634858940.,0,0},
2024040e9f2SJed Brown                  {6485989280629./16251701735622.,-4246266847089./9704473918619.,10755448449292./10357097424841.,0}},
203a3a57f36SJed Brown       At[4][4] = {{0,0,0,0},
2044040e9f2SJed Brown                   {1767732205903./4055673282236.,1767732205903./4055673282236.,0,0},
2054040e9f2SJed Brown                   {2746238789719./10658868560708.,-640167445237./6845629431997.,1767732205903./4055673282236.,0},
2064040e9f2SJed Brown                   {1471266399579./7840856788654.,-4482444167858./7529755066697.,11266239266428./11593286722821.,1767732205903./4055673282236.}},
2074040e9f2SJed Brown       binterpt[4][2] = {{4655552711362./22874653954995., -215264564351./13552729205753.},
2084040e9f2SJed Brown                         {-18682724506714./9892148508045.,17870216137069./13817060693119.},
2094040e9f2SJed Brown                         {34259539580243./13192909600954.,-28141676662227./17317692491321.},
2104040e9f2SJed Brown                         {584795268549./6622622206610.,   2508943948391./7218656332882.}};
211cd652676SJed Brown     ierr = TSARKIMEXRegister(TSARKIMEX3,3,4,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr);
212a3a57f36SJed Brown   }
213a3a57f36SJed Brown   {
214a3a57f36SJed Brown     const PetscReal
215*6cf0794eSJed Brown       A[5][5] = {{0,0,0,0},
216*6cf0794eSJed Brown                  {1./2,0,0,0,0},
217*6cf0794eSJed Brown                  {11./18,1./18,0,0,0},
218*6cf0794eSJed Brown                  {5./6,-5./6,.5,0,0},
219*6cf0794eSJed Brown                  {1./4,7./4,3./4,-7./4,0}},
220*6cf0794eSJed Brown       At[5][5] = {{0,0,0,0,0},
221*6cf0794eSJed Brown                   {0,1./2,0,0,0},
222*6cf0794eSJed Brown                   {0,1./6,1./2,0,0},
223*6cf0794eSJed Brown                   {0,-1./2,1./2,1./2,0},
224*6cf0794eSJed Brown                   {0,3./2,-3./2,1./2,1./2}};
225*6cf0794eSJed Brown     ierr = TSARKIMEXRegister(TSARKIMEXARS443,3,5,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,0,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
226*6cf0794eSJed Brown   }
227*6cf0794eSJed Brown   {
228*6cf0794eSJed Brown     const PetscReal
229*6cf0794eSJed Brown       A[5][5] = {{0,0,0,0},
230*6cf0794eSJed Brown                  {1,0,0,0,0},
231*6cf0794eSJed Brown                  {4./9,2./9,0,0,0},
232*6cf0794eSJed Brown                  {1./4,0,3./4,0,0},
233*6cf0794eSJed Brown                  {1./4,0,3./5,0,0}},
234*6cf0794eSJed Brown       At[5][5] = {{0,0,0,0},
235*6cf0794eSJed Brown                   {.5,.5,0,0,0},
236*6cf0794eSJed Brown                   {5./18,-1./9,.5,0,0},
237*6cf0794eSJed Brown                   {.5,0,0,.5,0},
238*6cf0794eSJed Brown                   {.25,0,.75,-.5,.5}};
239*6cf0794eSJed Brown     ierr = TSARKIMEXRegister(TSARKIMEXBPR3,3,5,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,0,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
240*6cf0794eSJed Brown   }
241*6cf0794eSJed Brown   {
242*6cf0794eSJed Brown     const PetscReal
243a3a57f36SJed Brown       A[6][6] = {{0,0,0,0,0,0},
244a3a57f36SJed Brown                  {1./2,0,0,0,0,0},
2454040e9f2SJed Brown                  {13861./62500.,6889./62500.,0,0,0,0},
2464040e9f2SJed Brown                  {-116923316275./2393684061468.,-2731218467317./15368042101831.,9408046702089./11113171139209.,0,0,0},
2474040e9f2SJed Brown                  {-451086348788./2902428689909.,-2682348792572./7519795681897.,12662868775082./11960479115383.,3355817975965./11060851509271.,0,0},
2484040e9f2SJed Brown                  {647845179188./3216320057751.,73281519250./8382639484533.,552539513391./3454668386233.,3354512671639./8306763924573.,4040./17871.,0}},
249a3a57f36SJed Brown       At[6][6] = {{0,0,0,0,0,0},
250a3a57f36SJed Brown                   {1./4,1./4,0,0,0,0},
2514040e9f2SJed Brown                   {8611./62500.,-1743./31250.,1./4,0,0,0},
2524040e9f2SJed Brown                   {5012029./34652500.,-654441./2922500.,174375./388108.,1./4,0,0},
2534040e9f2SJed Brown                   {15267082809./155376265600.,-71443401./120774400.,730878875./902184768.,2285395./8070912.,1./4,0},
2544040e9f2SJed Brown                   {82889./524892.,0,15625./83664.,69875./102672.,-2260./8211,1./4}},
2554040e9f2SJed Brown       binterpt[6][3] = {{6943876665148./7220017795957.,-54480133./30881146.,6818779379841./7100303317025.},
256cd652676SJed Brown                         {0,0,0},
2574040e9f2SJed Brown                         {7640104374378./9702883013639.,-11436875./14766696.,2173542590792./12501825683035.},
2584040e9f2SJed Brown                         {-20649996744609./7521556579894.,174696575./18121608.,-31592104683404./5083833661969.},
2594040e9f2SJed Brown                         {8854892464581./2390941311638.,-12120380./966161.,61146701046299./7138195549469.},
2604040e9f2SJed Brown                         {-11397109935349./6675773540249.,3843./706.,-17219254887155./4939391667607.}};
261cd652676SJed Brown     ierr = TSARKIMEXRegister(TSARKIMEX4,4,6,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,3,binterpt[0],PETSC_NULL);CHKERRQ(ierr);
262a3a57f36SJed Brown   }
263a3a57f36SJed Brown   {
264a3a57f36SJed Brown     const PetscReal
265a3a57f36SJed Brown       A[8][8] = {{0,0,0,0,0,0,0,0},
266a3a57f36SJed Brown                  {41./100,0,0,0,0,0,0,0},
2674040e9f2SJed Brown                  {367902744464./2072280473677.,677623207551./8224143866563.,0,0,0,0,0,0},
2684040e9f2SJed Brown                  {1268023523408./10340822734521.,0,1029933939417./13636558850479.,0,0,0,0,0},
2694040e9f2SJed Brown                  {14463281900351./6315353703477.,0,66114435211212./5879490589093.,-54053170152839./4284798021562.,0,0,0,0},
2704040e9f2SJed Brown                  {14090043504691./34967701212078.,0,15191511035443./11219624916014.,-18461159152457./12425892160975.,-281667163811./9011619295870.,0,0,0},
2714040e9f2SJed Brown                  {19230459214898./13134317526959.,0,21275331358303./2942455364971.,-38145345988419./4862620318723.,-1./8,-1./8,0,0},
2724040e9f2SJed Brown                  {-19977161125411./11928030595625.,0,-40795976796054./6384907823539.,177454434618887./12078138498510.,782672205425./8267701900261.,-69563011059811./9646580694205.,7356628210526./4942186776405.,0}},
273a3a57f36SJed Brown       At[8][8] = {{0,0,0,0,0,0,0,0},
2744040e9f2SJed Brown                   {41./200.,41./200.,0,0,0,0,0,0},
2754040e9f2SJed Brown                   {41./400.,-567603406766./11931857230679.,41./200.,0,0,0,0,0},
2764040e9f2SJed Brown                   {683785636431./9252920307686.,0,-110385047103./1367015193373.,41./200.,0,0,0,0},
2774040e9f2SJed Brown                   {3016520224154./10081342136671.,0,30586259806659./12414158314087.,-22760509404356./11113319521817.,41./200.,0,0,0},
2784040e9f2SJed Brown                   {218866479029./1489978393911.,0,638256894668./5436446318841.,-1179710474555./5321154724896.,-60928119172./8023461067671.,41./200.,0,0},
2794040e9f2SJed Brown                   {1020004230633./5715676835656.,0,25762820946817./25263940353407.,-2161375909145./9755907335909.,-211217309593./5846859502534.,-4269925059573./7827059040749.,41./200,0},
2804040e9f2SJed Brown                   {-872700587467./9133579230613.,0,0,22348218063261./9555858737531.,-1143369518992./8141816002931.,-39379526789629./19018526304540.,32727382324388./42900044865799.,41./200.}},
2814040e9f2SJed Brown       binterpt[8][3] = {{-17674230611817./10670229744614. ,  43486358583215./12773830924787. , -9257016797708./5021505065439.},
282cd652676SJed Brown                         {0                               ,  0                              , 0                            },
283cd652676SJed Brown                         {0                               ,  0                              , 0                            },
2844040e9f2SJed Brown                         {65168852399939./7868540260826.  ,  -91478233927265./11067650958493., 26096422576131./11239449250142.},
2854040e9f2SJed Brown                         {15494834004392./5936557850923.  ,  -79368583304911./10890268929626., 92396832856987./20362823103730.},
2864040e9f2SJed Brown                         {-99329723586156./26959484932159.,  -12239297817655./9152339842473. , 30029262896817./10175596800299.},
2874040e9f2SJed Brown                         {-19024464361622./5461577185407. ,  115839755401235./10719374521269., -26136350496073./3983972220547.},
2884040e9f2SJed Brown                         {-6511271360970./6095937251113.  ,  5843115559534./2180450260947.   , -5289405421727./3760307252460. }};
289cd652676SJed Brown     ierr = TSARKIMEXRegister(TSARKIMEX5,5,8,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,3,binterpt[0],PETSC_NULL);CHKERRQ(ierr);
290a3a57f36SJed Brown   }
291a3a57f36SJed Brown 
2928a381b04SJed Brown   PetscFunctionReturn(0);
2938a381b04SJed Brown }
2948a381b04SJed Brown 
2958a381b04SJed Brown #undef __FUNCT__
2968a381b04SJed Brown #define __FUNCT__ "TSARKIMEXRegisterDestroy"
2978a381b04SJed Brown /*@C
2988a381b04SJed Brown    TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by TSARKIMEXRegister().
2998a381b04SJed Brown 
3008a381b04SJed Brown    Not Collective
3018a381b04SJed Brown 
3028a381b04SJed Brown    Level: advanced
3038a381b04SJed Brown 
3048a381b04SJed Brown .keywords: TSARKIMEX, register, destroy
3058a381b04SJed Brown .seealso: TSARKIMEXRegister(), TSARKIMEXRegisterAll(), TSARKIMEXRegisterDynamic()
3068a381b04SJed Brown @*/
3078a381b04SJed Brown PetscErrorCode TSARKIMEXRegisterDestroy(void)
3088a381b04SJed Brown {
3098a381b04SJed Brown   PetscErrorCode ierr;
3108a381b04SJed Brown   ARKTableauLink link;
3118a381b04SJed Brown 
3128a381b04SJed Brown   PetscFunctionBegin;
3138a381b04SJed Brown   while ((link = ARKTableauList)) {
3148a381b04SJed Brown     ARKTableau t = &link->tab;
3158a381b04SJed Brown     ARKTableauList = link->next;
3168a381b04SJed Brown     ierr = PetscFree6(t->At,t->bt,t->ct,t->A,t->b,t->c);CHKERRQ(ierr);
317cd652676SJed Brown     ierr = PetscFree2(t->binterpt,t->binterp);CHKERRQ(ierr);
3188a381b04SJed Brown     ierr = PetscFree(t->name);CHKERRQ(ierr);
3198a381b04SJed Brown     ierr = PetscFree(link);CHKERRQ(ierr);
3208a381b04SJed Brown   }
3218a381b04SJed Brown   TSARKIMEXRegisterAllCalled = PETSC_FALSE;
3228a381b04SJed Brown   PetscFunctionReturn(0);
3238a381b04SJed Brown }
3248a381b04SJed Brown 
3258a381b04SJed Brown #undef __FUNCT__
3268a381b04SJed Brown #define __FUNCT__ "TSARKIMEXInitializePackage"
3278a381b04SJed Brown /*@C
3288a381b04SJed Brown   TSARKIMEXInitializePackage - This function initializes everything in the TSARKIMEX package. It is called
3298a381b04SJed Brown   from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_ARKIMEX()
3308a381b04SJed Brown   when using static libraries.
3318a381b04SJed Brown 
3328a381b04SJed Brown   Input Parameter:
3338a381b04SJed Brown   path - The dynamic library path, or PETSC_NULL
3348a381b04SJed Brown 
3358a381b04SJed Brown   Level: developer
3368a381b04SJed Brown 
3378a381b04SJed Brown .keywords: TS, TSARKIMEX, initialize, package
3388a381b04SJed Brown .seealso: PetscInitialize()
3398a381b04SJed Brown @*/
3408a381b04SJed Brown PetscErrorCode TSARKIMEXInitializePackage(const char path[])
3418a381b04SJed Brown {
3428a381b04SJed Brown   PetscErrorCode ierr;
3438a381b04SJed Brown 
3448a381b04SJed Brown   PetscFunctionBegin;
3458a381b04SJed Brown   if (TSARKIMEXPackageInitialized) PetscFunctionReturn(0);
3468a381b04SJed Brown   TSARKIMEXPackageInitialized = PETSC_TRUE;
3478a381b04SJed Brown   ierr = TSARKIMEXRegisterAll();CHKERRQ(ierr);
3488a381b04SJed Brown   ierr = PetscRegisterFinalize(TSARKIMEXFinalizePackage);CHKERRQ(ierr);
3498a381b04SJed Brown   PetscFunctionReturn(0);
3508a381b04SJed Brown }
3518a381b04SJed Brown 
3528a381b04SJed Brown #undef __FUNCT__
3538a381b04SJed Brown #define __FUNCT__ "TSARKIMEXFinalizePackage"
3548a381b04SJed Brown /*@C
3558a381b04SJed Brown   TSARKIMEXFinalizePackage - This function destroys everything in the TSARKIMEX package. It is
3568a381b04SJed Brown   called from PetscFinalize().
3578a381b04SJed Brown 
3588a381b04SJed Brown   Level: developer
3598a381b04SJed Brown 
3608a381b04SJed Brown .keywords: Petsc, destroy, package
3618a381b04SJed Brown .seealso: PetscFinalize()
3628a381b04SJed Brown @*/
3638a381b04SJed Brown PetscErrorCode TSARKIMEXFinalizePackage(void)
3648a381b04SJed Brown {
3658a381b04SJed Brown   PetscErrorCode ierr;
3668a381b04SJed Brown 
3678a381b04SJed Brown   PetscFunctionBegin;
3688a381b04SJed Brown   TSARKIMEXPackageInitialized = PETSC_FALSE;
3698a381b04SJed Brown   ierr = TSARKIMEXRegisterDestroy();CHKERRQ(ierr);
3708a381b04SJed Brown   PetscFunctionReturn(0);
3718a381b04SJed Brown }
3728a381b04SJed Brown 
3738a381b04SJed Brown #undef __FUNCT__
3748a381b04SJed Brown #define __FUNCT__ "TSARKIMEXRegister"
375cd652676SJed Brown /*@C
376cd652676SJed Brown    TSARKIMEXRegister - register an ARK IMEX scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation
377cd652676SJed Brown 
378cd652676SJed Brown    Not Collective, but the same schemes should be registered on all processes on which they will be used
379cd652676SJed Brown 
380cd652676SJed Brown    Input Parameters:
381cd652676SJed Brown +  name - identifier for method
382cd652676SJed Brown .  order - approximation order of method
383cd652676SJed Brown .  s - number of stages, this is the dimension of the matrices below
384cd652676SJed Brown .  At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major)
385cd652676SJed Brown .  bt - Butcher table for completing the stiff part of the step (dimension s; PETSC_NULL to use the last row of At)
386cd652676SJed Brown .  ct - Abscissa of each stiff stage (dimension s, PETSC_NULL to use row sums of At)
387cd652676SJed Brown .  A - Non-stiff stage coefficients (dimension s*s, row-major)
388cd652676SJed Brown .  b - Non-stiff step completion table (dimension s; PETSC_NULL to use last row of At)
389cd652676SJed Brown .  c - Non-stiff abscissa (dimension s; PETSC_NULL to use row sums of A)
390cd652676SJed Brown .  pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp
391cd652676SJed Brown .  binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp)
392cd652676SJed Brown -  binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; PETSC_NULL to reuse binterpt)
393cd652676SJed Brown 
394cd652676SJed Brown    Notes:
395cd652676SJed Brown    Several ARK IMEX methods are provided, this function is only needed to create new methods.
396cd652676SJed Brown 
397cd652676SJed Brown    Level: advanced
398cd652676SJed Brown 
399cd652676SJed Brown .keywords: TS, register
400cd652676SJed Brown 
401cd652676SJed Brown .seealso: TSARKIMEX
402cd652676SJed Brown @*/
4038a381b04SJed Brown PetscErrorCode TSARKIMEXRegister(const TSARKIMEXType name,PetscInt order,PetscInt s,
4048a381b04SJed Brown                                  const PetscReal At[],const PetscReal bt[],const PetscReal ct[],
405cd652676SJed Brown                                  const PetscReal A[],const PetscReal b[],const PetscReal c[],
406cd652676SJed Brown                                  PetscInt pinterp,const PetscReal binterpt[],const PetscReal binterp[])
4078a381b04SJed Brown {
4088a381b04SJed Brown   PetscErrorCode ierr;
4098a381b04SJed Brown   ARKTableauLink link;
4108a381b04SJed Brown   ARKTableau t;
4118a381b04SJed Brown   PetscInt i,j;
4128a381b04SJed Brown 
4138a381b04SJed Brown   PetscFunctionBegin;
4148a381b04SJed Brown   ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr);
415cd652676SJed Brown   ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr);
4168a381b04SJed Brown   t = &link->tab;
4178a381b04SJed Brown   ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
4188a381b04SJed Brown   t->order = order;
4198a381b04SJed Brown   t->s = s;
4208a381b04SJed Brown   ierr = PetscMalloc6(s*s,PetscReal,&t->At,s,PetscReal,&t->bt,s,PetscReal,&t->ct,s*s,PetscReal,&t->A,s,PetscReal,&t->b,s,PetscReal,&t->c);CHKERRQ(ierr);
4218a381b04SJed Brown   ierr = PetscMemcpy(t->At,At,s*s*sizeof(At[0]));CHKERRQ(ierr);
4228a381b04SJed Brown   ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr);
4238a381b04SJed Brown   if (bt) {ierr = PetscMemcpy(t->bt,bt,s*sizeof(bt[0]));CHKERRQ(ierr);}
4248a381b04SJed Brown   else for (i=0; i<s; i++) t->bt[i] = At[(s-1)*s+i];
4258a381b04SJed Brown   if (b) {ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr);}
4268a381b04SJed Brown   else for (i=0; i<s; i++) t->b[i] = At[(s-1)*s+i];
4278a381b04SJed Brown   if (ct) {ierr = PetscMemcpy(t->ct,ct,s*sizeof(ct[0]));CHKERRQ(ierr);}
4288a381b04SJed Brown   else for (i=0; i<s; i++) for (j=0,t->ct[i]=0; j<s; j++) t->ct[i] += At[i*s+j];
4298a381b04SJed Brown   if (c) {ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr);}
4308a381b04SJed Brown   else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
4314f385281SJed Brown   t->pinterp = pinterp;
432cd652676SJed Brown   ierr = PetscMalloc2(s*pinterp,PetscReal,&t->binterpt,s*pinterp,PetscReal,&t->binterp);CHKERRQ(ierr);
433cd652676SJed Brown   ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr);
434cd652676SJed Brown   ierr = PetscMemcpy(t->binterp,binterp?binterp:binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr);
4358a381b04SJed Brown   link->next = ARKTableauList;
4368a381b04SJed Brown   ARKTableauList = link;
4378a381b04SJed Brown   PetscFunctionReturn(0);
4388a381b04SJed Brown }
4398a381b04SJed Brown 
4408a381b04SJed Brown #undef __FUNCT__
4418a381b04SJed Brown #define __FUNCT__ "TSStep_ARKIMEX"
4428a381b04SJed Brown static PetscErrorCode TSStep_ARKIMEX(TS ts)
4438a381b04SJed Brown {
4448a381b04SJed Brown   TS_ARKIMEX          *ark = (TS_ARKIMEX*)ts->data;
4458a381b04SJed Brown   ARKTableau          tab  = ark->tableau;
4468a381b04SJed Brown   const PetscInt      s    = tab->s;
4478a381b04SJed Brown   const PetscReal     *At  = tab->At,*A = tab->A,*bt = tab->bt,*b = tab->b,*ct = tab->ct,*c = tab->c;
448406d0ec2SJed Brown   PetscScalar         *w   = ark->work;
4498a381b04SJed Brown   Vec                 *Y   = ark->Y,*YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS,Ydot = ark->Ydot,W = ark->Work,Z = ark->Z;
4508a381b04SJed Brown   SNES                snes;
451f1b97656SJed Brown   SNESConvergedReason snesreason;
4528a381b04SJed Brown   PetscInt            i,j,its,lits;
453cdbf8f93SLisandro Dalcin   PetscReal           next_time_step;
4548a381b04SJed Brown   PetscReal           h,t;
4558a381b04SJed Brown   PetscErrorCode      ierr;
4568a381b04SJed Brown 
4578a381b04SJed Brown   PetscFunctionBegin;
4588a381b04SJed Brown   ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
459cdbf8f93SLisandro Dalcin   next_time_step = ts->time_step;
460cdbf8f93SLisandro Dalcin   h = ts->time_step;
4618a381b04SJed Brown   t = ts->ptime;
4628a381b04SJed Brown 
4638a381b04SJed Brown   for (i=0; i<s; i++) {
4648a381b04SJed Brown     if (At[i*s+i] == 0) {           /* This stage is explicit */
4658a381b04SJed Brown       ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr);
4668a381b04SJed Brown       for (j=0; j<i; j++) w[j] = -h*At[i*s+j];
4678a381b04SJed Brown       ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr);
4688a381b04SJed Brown       for (j=0; j<i; j++) w[j] = h*A[i*s+j];
4698a381b04SJed Brown       ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr);
4708a381b04SJed Brown     } else {
4718a381b04SJed Brown       ark->stage_time = t + h*ct[i];
4728a381b04SJed Brown       ark->shift = 1./(h*At[i*s+i]);
4738a381b04SJed Brown       /* Affine part */
4748a381b04SJed Brown       ierr = VecZeroEntries(W);CHKERRQ(ierr);
4758a381b04SJed Brown       for (j=0; j<i; j++) w[j] = h*A[i*s+j];
4768a381b04SJed Brown       ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr);
4778a381b04SJed Brown       /* Ydot = shift*(Y-Z) */
4788a381b04SJed Brown       ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr);
4794f385281SJed Brown       for (j=0; j<i; j++) w[j] = -h*At[i*s+j];
4804f385281SJed Brown       ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr);
4818a381b04SJed Brown       /* Initial guess taken from last stage */
4828a381b04SJed Brown       ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr);
4838a381b04SJed Brown       ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr);
4848a381b04SJed Brown       ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
4858a381b04SJed Brown       ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr);
486f1b97656SJed Brown       ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr);
4878a381b04SJed Brown       ts->nonlinear_its += its; ts->linear_its += lits;
488f1b97656SJed Brown       if (snesreason < 0 && ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) {
489f1b97656SJed Brown         ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
490f1b97656SJed Brown         ierr = PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);CHKERRQ(ierr);
491f1b97656SJed Brown         PetscFunctionReturn(0);
492f1b97656SJed Brown       }
4938a381b04SJed Brown     }
4948a381b04SJed Brown     ierr = VecZeroEntries(Ydot);CHKERRQ(ierr);
4954cc180ffSJed Brown     ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr);
4964cc180ffSJed Brown     if (ark->imex) {
4978a381b04SJed Brown       ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr);
4984cc180ffSJed Brown     } else {
4994cc180ffSJed Brown       ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr);
5004cc180ffSJed Brown     }
5018a381b04SJed Brown   }
502a17b68daSJed Brown   for (j=0; j<s; j++) w[j] = -h*bt[j];
5038a381b04SJed Brown   ierr = VecMAXPY(ts->vec_sol,s,w,YdotI);CHKERRQ(ierr);
5048a381b04SJed Brown   for (j=0; j<s; j++) w[j] = h*b[j];
5058a381b04SJed Brown   ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr);
5068a381b04SJed Brown 
5078a381b04SJed Brown   ts->ptime += ts->time_step;
508cdbf8f93SLisandro Dalcin   ts->time_step = next_time_step;
5098a381b04SJed Brown   ts->steps++;
5108a381b04SJed Brown   PetscFunctionReturn(0);
5118a381b04SJed Brown }
5128a381b04SJed Brown 
513cd652676SJed Brown #undef __FUNCT__
514cd652676SJed Brown #define __FUNCT__ "TSInterpolate_ARKIMEX"
515cd652676SJed Brown static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X)
516cd652676SJed Brown {
517cd652676SJed Brown   TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
5184f385281SJed Brown   PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j;
5194f385281SJed Brown   PetscReal tt,t = (itime - ts->ptime)/ts->time_step + 1; /* In the interval [0,1] */
520cd652676SJed Brown   PetscScalar *bt,*b;
521cd652676SJed Brown   const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp;
522cd652676SJed Brown   PetscErrorCode ierr;
523cd652676SJed Brown 
524cd652676SJed Brown   PetscFunctionBegin;
525cd652676SJed Brown   if (!Bt || !B) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name);
526cd652676SJed Brown   ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr);
527cd652676SJed Brown   for (i=0; i<s; i++) bt[i] = b[i] = 0;
5284f385281SJed Brown   for (j=0,tt=t; j<pinterp; j++,tt*=t) {
529cd652676SJed Brown     for (i=0; i<s; i++) {
5304f385281SJed Brown       bt[i] += ts->time_step * Bt[i*pinterp+j] * tt * -1.0;
5314f385281SJed Brown       b[i]  += ts->time_step * B[i*pinterp+j] * tt;
532cd652676SJed Brown     }
533cd652676SJed Brown   }
534cd652676SJed Brown   if (ark->tableau->At[0*s+0] != 0.0) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_SUP,"First stage not explicit so starting stage not saved");
535cd652676SJed Brown   ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr);
536cd652676SJed Brown   ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr);
537cd652676SJed Brown   ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr);
538cd652676SJed Brown   ierr = PetscFree2(bt,b);CHKERRQ(ierr);
539cd652676SJed Brown   PetscFunctionReturn(0);
540cd652676SJed Brown }
541cd652676SJed Brown 
5428a381b04SJed Brown /*------------------------------------------------------------*/
5438a381b04SJed Brown #undef __FUNCT__
5448a381b04SJed Brown #define __FUNCT__ "TSReset_ARKIMEX"
5458a381b04SJed Brown static PetscErrorCode TSReset_ARKIMEX(TS ts)
5468a381b04SJed Brown {
5478a381b04SJed Brown   TS_ARKIMEX      *ark = (TS_ARKIMEX*)ts->data;
5488a381b04SJed Brown   PetscInt        s;
5498a381b04SJed Brown   PetscErrorCode  ierr;
5508a381b04SJed Brown 
5518a381b04SJed Brown   PetscFunctionBegin;
5528a381b04SJed Brown   if (!ark->tableau) PetscFunctionReturn(0);
5538a381b04SJed Brown    s = ark->tableau->s;
5548a381b04SJed Brown   ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr);
5558a381b04SJed Brown   ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr);
5568a381b04SJed Brown   ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr);
5578a381b04SJed Brown   ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr);
5588a381b04SJed Brown   ierr = VecDestroy(&ark->Work);CHKERRQ(ierr);
5598a381b04SJed Brown   ierr = VecDestroy(&ark->Z);CHKERRQ(ierr);
5608a381b04SJed Brown   ierr = PetscFree(ark->work);CHKERRQ(ierr);
5618a381b04SJed Brown   PetscFunctionReturn(0);
5628a381b04SJed Brown }
5638a381b04SJed Brown 
5648a381b04SJed Brown #undef __FUNCT__
5658a381b04SJed Brown #define __FUNCT__ "TSDestroy_ARKIMEX"
5668a381b04SJed Brown static PetscErrorCode TSDestroy_ARKIMEX(TS ts)
5678a381b04SJed Brown {
5688a381b04SJed Brown   PetscErrorCode  ierr;
5698a381b04SJed Brown 
5708a381b04SJed Brown   PetscFunctionBegin;
5718a381b04SJed Brown   ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr);
5728a381b04SJed Brown   ierr = PetscFree(ts->data);CHKERRQ(ierr);
5738a381b04SJed Brown   ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","",PETSC_NULL);CHKERRQ(ierr);
5748a381b04SJed Brown   ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","",PETSC_NULL);CHKERRQ(ierr);
575995b3938SJed Brown   ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","",PETSC_NULL);CHKERRQ(ierr);
5768a381b04SJed Brown   PetscFunctionReturn(0);
5778a381b04SJed Brown }
5788a381b04SJed Brown 
5798a381b04SJed Brown /*
5808a381b04SJed Brown   This defines the nonlinear equation that is to be solved with SNES
5818a381b04SJed Brown   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
5828a381b04SJed Brown */
5838a381b04SJed Brown #undef __FUNCT__
5848a381b04SJed Brown #define __FUNCT__ "SNESTSFormFunction_ARKIMEX"
5858a381b04SJed Brown static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts)
5868a381b04SJed Brown {
5878a381b04SJed Brown   TS_ARKIMEX     *ark = (TS_ARKIMEX*)ts->data;
5888a381b04SJed Brown   PetscErrorCode ierr;
5898a381b04SJed Brown 
5908a381b04SJed Brown   PetscFunctionBegin;
5918a381b04SJed Brown   ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */
5924cc180ffSJed Brown   ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,F,ark->imex);CHKERRQ(ierr);
5938a381b04SJed Brown   PetscFunctionReturn(0);
5948a381b04SJed Brown }
5958a381b04SJed Brown 
5968a381b04SJed Brown #undef __FUNCT__
5978a381b04SJed Brown #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX"
5988a381b04SJed Brown static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts)
5998a381b04SJed Brown {
6008a381b04SJed Brown   TS_ARKIMEX       *ark = (TS_ARKIMEX*)ts->data;
6018a381b04SJed Brown   PetscErrorCode ierr;
6028a381b04SJed Brown 
6038a381b04SJed Brown   PetscFunctionBegin;
6048a381b04SJed Brown   /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */
6058a381b04SJed Brown   ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr);
6068a381b04SJed Brown   PetscFunctionReturn(0);
6078a381b04SJed Brown }
6088a381b04SJed Brown 
6098a381b04SJed Brown #undef __FUNCT__
6108a381b04SJed Brown #define __FUNCT__ "TSSetUp_ARKIMEX"
6118a381b04SJed Brown static PetscErrorCode TSSetUp_ARKIMEX(TS ts)
6128a381b04SJed Brown {
6138a381b04SJed Brown   TS_ARKIMEX     *ark = (TS_ARKIMEX*)ts->data;
6148a381b04SJed Brown   ARKTableau     tab  = ark->tableau;
6158a381b04SJed Brown   PetscInt       s = tab->s;
6168a381b04SJed Brown   PetscErrorCode ierr;
6178a381b04SJed Brown 
6188a381b04SJed Brown   PetscFunctionBegin;
6198a381b04SJed Brown   if (!ark->tableau) {
620e24355feSJed Brown     ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);
6218a381b04SJed Brown   }
6228a381b04SJed Brown   ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr);
6238a381b04SJed Brown   ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr);
6248a381b04SJed Brown   ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr);
6258a381b04SJed Brown   ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr);
6268a381b04SJed Brown   ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr);
6278a381b04SJed Brown   ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr);
6288a381b04SJed Brown   ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr);
6298a381b04SJed Brown   PetscFunctionReturn(0);
6308a381b04SJed Brown }
6318a381b04SJed Brown /*------------------------------------------------------------*/
6328a381b04SJed Brown 
6338a381b04SJed Brown #undef __FUNCT__
6348a381b04SJed Brown #define __FUNCT__ "TSSetFromOptions_ARKIMEX"
6358a381b04SJed Brown static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts)
6368a381b04SJed Brown {
6374cc180ffSJed Brown   TS_ARKIMEX     *ark = (TS_ARKIMEX*)ts->data;
6388a381b04SJed Brown   PetscErrorCode ierr;
6398a381b04SJed Brown   char           arktype[256];
6408a381b04SJed Brown 
6418a381b04SJed Brown   PetscFunctionBegin;
6428a381b04SJed Brown   ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr);
6438a381b04SJed Brown   {
6448a381b04SJed Brown     ARKTableauLink link;
6458a381b04SJed Brown     PetscInt count,choice;
6468a381b04SJed Brown     PetscBool flg;
6478a381b04SJed Brown     const char **namelist;
6488a381b04SJed Brown     ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof arktype);CHKERRQ(ierr);
6498a381b04SJed Brown     for (link=ARKTableauList,count=0; link; link=link->next,count++) ;
6508a381b04SJed Brown     ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr);
6518a381b04SJed Brown     for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name;
6528a381b04SJed Brown     ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr);
6538a381b04SJed Brown     ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr);
6548a381b04SJed Brown     ierr = PetscFree(namelist);CHKERRQ(ierr);
6554cc180ffSJed Brown     flg = (PetscBool)!ark->imex;
6564cc180ffSJed Brown     ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr);
6574cc180ffSJed Brown     ark->imex = (PetscBool)!flg;
658d52bd9f3SBarry Smith     ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr);
6598a381b04SJed Brown   }
6608a381b04SJed Brown   ierr = PetscOptionsTail();CHKERRQ(ierr);
6618a381b04SJed Brown   PetscFunctionReturn(0);
6628a381b04SJed Brown }
6638a381b04SJed Brown 
6648a381b04SJed Brown #undef __FUNCT__
6658a381b04SJed Brown #define __FUNCT__ "PetscFormatRealArray"
6668a381b04SJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[])
6678a381b04SJed Brown {
668257d2499SJed Brown   PetscErrorCode ierr;
669f1d86077SJed Brown   PetscInt i;
670f1d86077SJed Brown   size_t left,count;
6718a381b04SJed Brown   char *p;
6728a381b04SJed Brown 
6738a381b04SJed Brown   PetscFunctionBegin;
674f1d86077SJed Brown   for (i=0,p=buf,left=len; i<n; i++) {
675f1d86077SJed Brown     ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr);
6768a381b04SJed Brown     if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer");
6778a381b04SJed Brown     left -= count;
6788a381b04SJed Brown     p += count;
6798a381b04SJed Brown     *p++ = ' ';
6808a381b04SJed Brown   }
6818a381b04SJed Brown   p[i ? 0 : -1] = 0;
6828a381b04SJed Brown   PetscFunctionReturn(0);
6838a381b04SJed Brown }
6848a381b04SJed Brown 
6858a381b04SJed Brown #undef __FUNCT__
6868a381b04SJed Brown #define __FUNCT__ "TSView_ARKIMEX"
6878a381b04SJed Brown static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer)
6888a381b04SJed Brown {
6898a381b04SJed Brown   TS_ARKIMEX     *ark = (TS_ARKIMEX*)ts->data;
6908a381b04SJed Brown   ARKTableau     tab = ark->tableau;
6918a381b04SJed Brown   PetscBool      iascii;
6928a381b04SJed Brown   PetscErrorCode ierr;
6938a381b04SJed Brown 
6948a381b04SJed Brown   PetscFunctionBegin;
6958a381b04SJed Brown   ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
6968a381b04SJed Brown   if (iascii) {
6978a381b04SJed Brown     const TSARKIMEXType arktype;
6988a381b04SJed Brown     char buf[512];
6998a381b04SJed Brown     ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr);
7008a381b04SJed Brown     ierr = PetscViewerASCIIPrintf(viewer,"  ARK IMEX %s\n",arktype);CHKERRQ(ierr);
7018a381b04SJed Brown     ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr);
70231f6fcc0SJed Brown     ierr = PetscViewerASCIIPrintf(viewer,"  Stiff abscissa       ct = %s\n",buf);CHKERRQ(ierr);
70331f6fcc0SJed Brown     ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr);
70431f6fcc0SJed Brown     ierr = PetscViewerASCIIPrintf(viewer,"  Nonstiff abscissa     c = %s\n",buf);CHKERRQ(ierr);
7058a381b04SJed Brown   }
706d52bd9f3SBarry Smith   ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr);
7078a381b04SJed Brown   PetscFunctionReturn(0);
7088a381b04SJed Brown }
7098a381b04SJed Brown 
7108a381b04SJed Brown #undef __FUNCT__
7118a381b04SJed Brown #define __FUNCT__ "TSARKIMEXSetType"
7128a381b04SJed Brown /*@C
7138a381b04SJed Brown   TSARKIMEXSetType - Set the type of ARK IMEX scheme
7148a381b04SJed Brown 
7158a381b04SJed Brown   Logically collective
7168a381b04SJed Brown 
7178a381b04SJed Brown   Input Parameter:
7188a381b04SJed Brown +  ts - timestepping context
7198a381b04SJed Brown -  arktype - type of ARK-IMEX scheme
7208a381b04SJed Brown 
7218a381b04SJed Brown   Level: intermediate
7228a381b04SJed Brown 
7238a381b04SJed Brown .seealso: TSARKIMEXGetType()
7248a381b04SJed Brown @*/
7258a381b04SJed Brown PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype)
7268a381b04SJed Brown {
7278a381b04SJed Brown   PetscErrorCode ierr;
7288a381b04SJed Brown 
7298a381b04SJed Brown   PetscFunctionBegin;
7308a381b04SJed Brown   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
7318a381b04SJed Brown   ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr);
7328a381b04SJed Brown   PetscFunctionReturn(0);
7338a381b04SJed Brown }
7348a381b04SJed Brown 
7358a381b04SJed Brown #undef __FUNCT__
7368a381b04SJed Brown #define __FUNCT__ "TSARKIMEXGetType"
7378a381b04SJed Brown /*@C
7388a381b04SJed Brown   TSARKIMEXGetType - Get the type of ARK IMEX scheme
7398a381b04SJed Brown 
7408a381b04SJed Brown   Logically collective
7418a381b04SJed Brown 
7428a381b04SJed Brown   Input Parameter:
7438a381b04SJed Brown .  ts - timestepping context
7448a381b04SJed Brown 
7458a381b04SJed Brown   Output Parameter:
7468a381b04SJed Brown .  arktype - type of ARK-IMEX scheme
7478a381b04SJed Brown 
7488a381b04SJed Brown   Level: intermediate
7498a381b04SJed Brown 
7508a381b04SJed Brown .seealso: TSARKIMEXGetType()
7518a381b04SJed Brown @*/
7528a381b04SJed Brown PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype)
7538a381b04SJed Brown {
7548a381b04SJed Brown   PetscErrorCode ierr;
7558a381b04SJed Brown 
7568a381b04SJed Brown   PetscFunctionBegin;
7578a381b04SJed Brown   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
7588a381b04SJed Brown   ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr);
7598a381b04SJed Brown   PetscFunctionReturn(0);
7608a381b04SJed Brown }
7618a381b04SJed Brown 
7624cc180ffSJed Brown #undef __FUNCT__
7634cc180ffSJed Brown #define __FUNCT__ "TSARKIMEXSetFullyImplicit"
7644cc180ffSJed Brown /*@C
7654cc180ffSJed Brown   TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly
7664cc180ffSJed Brown 
7674cc180ffSJed Brown   Logically collective
7684cc180ffSJed Brown 
7694cc180ffSJed Brown   Input Parameter:
7704cc180ffSJed Brown +  ts - timestepping context
7714cc180ffSJed Brown -  flg - PETSC_TRUE for fully implicit
7724cc180ffSJed Brown 
7734cc180ffSJed Brown   Level: intermediate
7744cc180ffSJed Brown 
7754cc180ffSJed Brown .seealso: TSARKIMEXGetType()
7764cc180ffSJed Brown @*/
7774cc180ffSJed Brown PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg)
7784cc180ffSJed Brown {
7794cc180ffSJed Brown   PetscErrorCode ierr;
7804cc180ffSJed Brown 
7814cc180ffSJed Brown   PetscFunctionBegin;
7824cc180ffSJed Brown   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
7834cc180ffSJed Brown   ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr);
7844cc180ffSJed Brown   PetscFunctionReturn(0);
7854cc180ffSJed Brown }
7864cc180ffSJed Brown 
7878a381b04SJed Brown EXTERN_C_BEGIN
7888a381b04SJed Brown #undef __FUNCT__
7898a381b04SJed Brown #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX"
7908a381b04SJed Brown PetscErrorCode  TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype)
7918a381b04SJed Brown {
7928a381b04SJed Brown   TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
7938a381b04SJed Brown   PetscErrorCode ierr;
7948a381b04SJed Brown 
7958a381b04SJed Brown   PetscFunctionBegin;
7968a381b04SJed Brown   if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);}
7978a381b04SJed Brown   *arktype = ark->tableau->name;
7988a381b04SJed Brown   PetscFunctionReturn(0);
7998a381b04SJed Brown }
8008a381b04SJed Brown #undef __FUNCT__
8018a381b04SJed Brown #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX"
8028a381b04SJed Brown PetscErrorCode  TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype)
8038a381b04SJed Brown {
8048a381b04SJed Brown   TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
8058a381b04SJed Brown   PetscErrorCode ierr;
8068a381b04SJed Brown   PetscBool match;
8078a381b04SJed Brown   ARKTableauLink link;
8088a381b04SJed Brown 
8098a381b04SJed Brown   PetscFunctionBegin;
8108a381b04SJed Brown   if (ark->tableau) {
8118a381b04SJed Brown     ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr);
8128a381b04SJed Brown     if (match) PetscFunctionReturn(0);
8138a381b04SJed Brown   }
8148a381b04SJed Brown   for (link = ARKTableauList; link; link=link->next) {
8158a381b04SJed Brown     ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr);
8168a381b04SJed Brown     if (match) {
8178a381b04SJed Brown       ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr);
8188a381b04SJed Brown       ark->tableau = &link->tab;
8198a381b04SJed Brown       PetscFunctionReturn(0);
8208a381b04SJed Brown     }
8218a381b04SJed Brown   }
8228a381b04SJed Brown   SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype);
8238a381b04SJed Brown   PetscFunctionReturn(0);
8248a381b04SJed Brown }
8254cc180ffSJed Brown #undef __FUNCT__
8264cc180ffSJed Brown #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX"
8274cc180ffSJed Brown PetscErrorCode  TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg)
8284cc180ffSJed Brown {
8294cc180ffSJed Brown   TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data;
8304cc180ffSJed Brown 
8314cc180ffSJed Brown   PetscFunctionBegin;
8324cc180ffSJed Brown   ark->imex = (PetscBool)!flg;
8334cc180ffSJed Brown   PetscFunctionReturn(0);
8344cc180ffSJed Brown }
8358a381b04SJed Brown EXTERN_C_END
8368a381b04SJed Brown 
8378a381b04SJed Brown /* ------------------------------------------------------------ */
8388a381b04SJed Brown /*MC
8398a381b04SJed Brown       TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes
8408a381b04SJed Brown 
841fca742c7SJed Brown   These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly
842fca742c7SJed Brown   nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part
843fca742c7SJed Brown   of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction().
844fca742c7SJed Brown 
845fca742c7SJed Brown   Notes:
846c8058688SBarry Smith   The default is TSARKIMEX2E, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type
847c8058688SBarry Smith 
848fca742c7SJed Brown   This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X).
849fca742c7SJed Brown 
8508a381b04SJed Brown   Level: beginner
8518a381b04SJed Brown 
852c8058688SBarry Smith .seealso:  TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3,
85326885f59SBarry Smith            TSARKIMEX4, TSARKIMEX5, TSARKIMEXType, TSARKIMEXRegister()
8548a381b04SJed Brown 
8558a381b04SJed Brown M*/
8568a381b04SJed Brown EXTERN_C_BEGIN
8578a381b04SJed Brown #undef __FUNCT__
8588a381b04SJed Brown #define __FUNCT__ "TSCreate_ARKIMEX"
8598a381b04SJed Brown PetscErrorCode  TSCreate_ARKIMEX(TS ts)
8608a381b04SJed Brown {
8618a381b04SJed Brown   TS_ARKIMEX       *th;
8628a381b04SJed Brown   PetscErrorCode ierr;
8638a381b04SJed Brown 
8648a381b04SJed Brown   PetscFunctionBegin;
8658a381b04SJed Brown #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
8668a381b04SJed Brown   ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr);
8678a381b04SJed Brown #endif
8688a381b04SJed Brown 
8698a381b04SJed Brown   ts->ops->reset          = TSReset_ARKIMEX;
8708a381b04SJed Brown   ts->ops->destroy        = TSDestroy_ARKIMEX;
8718a381b04SJed Brown   ts->ops->view           = TSView_ARKIMEX;
8728a381b04SJed Brown   ts->ops->setup          = TSSetUp_ARKIMEX;
8738a381b04SJed Brown   ts->ops->step           = TSStep_ARKIMEX;
874cd652676SJed Brown   ts->ops->interpolate    = TSInterpolate_ARKIMEX;
8758a381b04SJed Brown   ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX;
8768a381b04SJed Brown   ts->ops->snesfunction   = SNESTSFormFunction_ARKIMEX;
8778a381b04SJed Brown   ts->ops->snesjacobian   = SNESTSFormJacobian_ARKIMEX;
8788a381b04SJed Brown 
8798a381b04SJed Brown   ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr);
8808a381b04SJed Brown   ts->data = (void*)th;
8814cc180ffSJed Brown   th->imex = PETSC_TRUE;
8828a381b04SJed Brown 
8838a381b04SJed Brown   ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr);
8848a381b04SJed Brown   ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr);
8854cc180ffSJed Brown   ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","TSARKIMEXSetFullyImplicit_ARKIMEX",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr);
8868a381b04SJed Brown   PetscFunctionReturn(0);
8878a381b04SJed Brown }
8888a381b04SJed Brown EXTERN_C_END
889