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 14108c343cSJed Brown static const TSARKIMEXType TSARKIMEXDefault = TSARKIMEX3; 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 */ 26108c343cSJed Brown PetscReal *bembedt,*bembed; /* Embedded formula of order one less (order-1) */ 27cd652676SJed Brown PetscReal *binterpt,*binterp; /* Dense output formula */ 28108c343cSJed Brown PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 298a381b04SJed Brown }; 308a381b04SJed Brown typedef struct _ARKTableauLink *ARKTableauLink; 318a381b04SJed Brown struct _ARKTableauLink { 328a381b04SJed Brown struct _ARKTableau tab; 338a381b04SJed Brown ARKTableauLink next; 348a381b04SJed Brown }; 358a381b04SJed Brown static ARKTableauLink ARKTableauList; 368a381b04SJed Brown 378a381b04SJed Brown typedef struct { 388a381b04SJed Brown ARKTableau tableau; 398a381b04SJed Brown Vec *Y; /* States computed during the step */ 408a381b04SJed Brown Vec *YdotI; /* Time derivatives for the stiff part */ 418a381b04SJed Brown Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 428a381b04SJed Brown Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 438a381b04SJed Brown Vec Work; /* Generic work vector */ 448a381b04SJed Brown Vec Z; /* Ydot = shift(Y-Z) */ 458a381b04SJed Brown PetscScalar *work; /* Scalar work */ 468a381b04SJed Brown PetscReal shift; 478a381b04SJed Brown PetscReal stage_time; 484cc180ffSJed Brown PetscBool imex; 49108c343cSJed Brown TSStepStatus status; 508a381b04SJed Brown } TS_ARKIMEX; 518a381b04SJed Brown 5264f491ddSJed Brown /*MC 5364f491ddSJed Brown TSARKIMEX2D - Second order ARK IMEX scheme with L-stable implicit part. 5464f491ddSJed Brown 5564f491ddSJed Brown This method has one explicit stage and two implicit stages. 5664f491ddSJed Brown 57b330ce4dSSatish Balay Level: advanced 58b330ce4dSSatish Balay 5964f491ddSJed Brown .seealso: TSARKIMEX 6064f491ddSJed Brown M*/ 6164f491ddSJed Brown /*MC 6264f491ddSJed Brown TSARKIMEX2E - Second order ARK IMEX scheme with L-stable implicit part. 6364f491ddSJed Brown 6464f491ddSJed Brown This method has one explicit stage and two implicit stages. It is is an optimal method developed by Emil Constantinescu. 6564f491ddSJed Brown 66b330ce4dSSatish Balay Level: advanced 67b330ce4dSSatish Balay 6864f491ddSJed Brown .seealso: TSARKIMEX 6964f491ddSJed Brown M*/ 7064f491ddSJed Brown /*MC 716cf0794eSJed Brown TSARKIMEXPRSSP2 - Second order SSP ARK IMEX scheme. 726cf0794eSJed Brown 736cf0794eSJed Brown This method has three implicit stages. 746cf0794eSJed Brown 756cf0794eSJed Brown References: 766cf0794eSJed 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 776cf0794eSJed Brown 786cf0794eSJed Brown This method is referred to as SSP2-(3,3,2) in http://arxiv.org/abs/1110.4375 796cf0794eSJed Brown 806cf0794eSJed Brown Level: advanced 816cf0794eSJed Brown 826cf0794eSJed Brown .seealso: TSARKIMEX 836cf0794eSJed Brown M*/ 846cf0794eSJed Brown /*MC 8564f491ddSJed Brown TSARKIMEX3 - Third order ARK IMEX scheme with L-stable implicit part. 8664f491ddSJed Brown 8764f491ddSJed Brown This method has one explicit stage and three implicit stages. 8864f491ddSJed Brown 8964f491ddSJed Brown References: 9064f491ddSJed Brown Kennedy and Carpenter 2003. 9164f491ddSJed Brown 92b330ce4dSSatish Balay Level: advanced 93b330ce4dSSatish Balay 9464f491ddSJed Brown .seealso: TSARKIMEX 9564f491ddSJed Brown M*/ 9664f491ddSJed Brown /*MC 976cf0794eSJed Brown TSARKIMEXARS443 - Third order ARK IMEX scheme. 986cf0794eSJed Brown 996cf0794eSJed Brown This method has one explicit stage and four implicit stages. 1006cf0794eSJed Brown 1016cf0794eSJed Brown References: 1026cf0794eSJed 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. 1036cf0794eSJed Brown 1046cf0794eSJed Brown This method is referred to as ARS(4,4,3) in http://arxiv.org/abs/1110.4375 1056cf0794eSJed Brown 1066cf0794eSJed Brown Level: advanced 1076cf0794eSJed Brown 1086cf0794eSJed Brown .seealso: TSARKIMEX 1096cf0794eSJed Brown M*/ 1106cf0794eSJed Brown /*MC 1116cf0794eSJed Brown TSARKIMEXBPR3 - Third order ARK IMEX scheme. 1126cf0794eSJed Brown 1136cf0794eSJed Brown This method has one explicit stage and four implicit stages. 1146cf0794eSJed Brown 1156cf0794eSJed Brown References: 1166cf0794eSJed Brown This method is referred to as ARK3 in http://arxiv.org/abs/1110.4375 1176cf0794eSJed Brown 1186cf0794eSJed Brown Level: advanced 1196cf0794eSJed Brown 1206cf0794eSJed Brown .seealso: TSARKIMEX 1216cf0794eSJed Brown M*/ 1226cf0794eSJed Brown /*MC 12364f491ddSJed Brown TSARKIMEX4 - Fourth order ARK IMEX scheme with L-stable implicit part. 12464f491ddSJed Brown 12564f491ddSJed Brown This method has one explicit stage and four implicit stages. 12664f491ddSJed Brown 12764f491ddSJed Brown References: 12864f491ddSJed Brown Kennedy and Carpenter 2003. 12964f491ddSJed Brown 130b330ce4dSSatish Balay Level: advanced 131b330ce4dSSatish Balay 13264f491ddSJed Brown .seealso: TSARKIMEX 13364f491ddSJed Brown M*/ 13464f491ddSJed Brown /*MC 13564f491ddSJed Brown TSARKIMEX5 - Fifth order ARK IMEX scheme with L-stable implicit part. 13664f491ddSJed Brown 13764f491ddSJed Brown This method has one explicit stage and five implicit stages. 13864f491ddSJed Brown 13964f491ddSJed Brown References: 14064f491ddSJed Brown Kennedy and Carpenter 2003. 14164f491ddSJed Brown 142b330ce4dSSatish Balay Level: advanced 143b330ce4dSSatish Balay 14464f491ddSJed Brown .seealso: TSARKIMEX 14564f491ddSJed Brown M*/ 14664f491ddSJed Brown 1478a381b04SJed Brown #undef __FUNCT__ 1488a381b04SJed Brown #define __FUNCT__ "TSARKIMEXRegisterAll" 1498a381b04SJed Brown /*@C 1508a381b04SJed Brown TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSARKIMEX 1518a381b04SJed Brown 152fca742c7SJed Brown Not Collective, but should be called by all processes which will need the schemes to be registered 1538a381b04SJed Brown 1548a381b04SJed Brown Level: advanced 1558a381b04SJed Brown 1568a381b04SJed Brown .keywords: TS, TSARKIMEX, register, all 1578a381b04SJed Brown 1588a381b04SJed Brown .seealso: TSARKIMEXRegisterDestroy() 1598a381b04SJed Brown @*/ 1608a381b04SJed Brown PetscErrorCode TSARKIMEXRegisterAll(void) 1618a381b04SJed Brown { 1628a381b04SJed Brown PetscErrorCode ierr; 1638a381b04SJed Brown 1648a381b04SJed Brown PetscFunctionBegin; 1658a381b04SJed Brown if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(0); 1668a381b04SJed Brown TSARKIMEXRegisterAllCalled = PETSC_TRUE; 1678a381b04SJed Brown 1688a381b04SJed Brown { 1698a381b04SJed Brown const PetscReal 1708a381b04SJed Brown A[3][3] = {{0,0,0}, 1718a381b04SJed Brown {0.41421356237309504880,0,0}, 1728a381b04SJed Brown {0.75,0.25,0}}, 1738a381b04SJed Brown At[3][3] = {{0,0,0}, 1748a381b04SJed Brown {0.12132034355964257320,0.29289321881345247560,0}, 17506db7b1cSJed Brown {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}}, 176108c343cSJed Brown bembedt[3] = {0.29289321881345247560,0.50000000000000000000,0.20710678118654752440}, 17706db7b1cSJed Brown binterpt[3][2] = {{1,-0.5},{0,0},{0,0.5}}; 178108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEX2D,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,bembedt,bembedt,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 1798a381b04SJed Brown } 18006db7b1cSJed Brown { /* Optimal for linear implicit part */ 18103403c7fSJed Brown const PetscReal s2 = PetscSqrtReal((PetscReal)2.0), 182a3a57f36SJed Brown A[3][3] = {{0,0,0}, 183a3a57f36SJed Brown {2-s2,0,0}, 184a3a57f36SJed Brown {(3-2*s2)/6,(3+2*s2)/6,0}}, 185a3a57f36SJed Brown At[3][3] = {{0,0,0}, 186a3a57f36SJed Brown {1-1/s2,1-1/s2,0}, 187cd652676SJed Brown {1/(2*s2),1/(2*s2),1-1/s2}}, 188cd652676SJed Brown binterpt[3][2] = {{1,-0.5},{0,0},{0,0.5}}; 189108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEX2E,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 190a3a57f36SJed Brown } 1916cf0794eSJed Brown { /* Optimal for linear implicit part */ 1926cf0794eSJed Brown const PetscReal 1936cf0794eSJed Brown A[3][3] = {{0,0,0}, 1946cf0794eSJed Brown {0.5,0,0}, 1956cf0794eSJed Brown {0.5,0.5,0}}, 1966cf0794eSJed Brown At[3][3] = {{0.25,0,0}, 1976cf0794eSJed Brown {0,0.25,0}, 1986cf0794eSJed Brown {1./3,1./3,1./3}}; 199108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEXPRSSP2,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,0,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2006cf0794eSJed Brown } 201a3a57f36SJed Brown { 202a3a57f36SJed Brown const PetscReal 203a3a57f36SJed Brown A[4][4] = {{0,0,0,0}, 2044040e9f2SJed Brown {1767732205903./2027836641118.,0,0,0}, 2054040e9f2SJed Brown {5535828885825./10492691773637.,788022342437./10882634858940.,0,0}, 2064040e9f2SJed Brown {6485989280629./16251701735622.,-4246266847089./9704473918619.,10755448449292./10357097424841.,0}}, 207a3a57f36SJed Brown At[4][4] = {{0,0,0,0}, 2084040e9f2SJed Brown {1767732205903./4055673282236.,1767732205903./4055673282236.,0,0}, 2094040e9f2SJed Brown {2746238789719./10658868560708.,-640167445237./6845629431997.,1767732205903./4055673282236.,0}, 2104040e9f2SJed Brown {1471266399579./7840856788654.,-4482444167858./7529755066697.,11266239266428./11593286722821.,1767732205903./4055673282236.}}, 211*cc46b9d1SJed Brown bembedt[4] = {2756255671327./12835298489170.,-10771552573575./22201958757719.,9247589265047./10645013368117.,2193209047091./5459859503100.}, 2124040e9f2SJed Brown binterpt[4][2] = {{4655552711362./22874653954995., -215264564351./13552729205753.}, 2134040e9f2SJed Brown {-18682724506714./9892148508045.,17870216137069./13817060693119.}, 2144040e9f2SJed Brown {34259539580243./13192909600954.,-28141676662227./17317692491321.}, 2154040e9f2SJed Brown {584795268549./6622622206610., 2508943948391./7218656332882.}}; 216108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEX3,3,4,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,bembedt,bembedt,2,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 217a3a57f36SJed Brown } 218a3a57f36SJed Brown { 219a3a57f36SJed Brown const PetscReal 2206cf0794eSJed Brown A[5][5] = {{0,0,0,0}, 2216cf0794eSJed Brown {1./2,0,0,0,0}, 2226cf0794eSJed Brown {11./18,1./18,0,0,0}, 2236cf0794eSJed Brown {5./6,-5./6,.5,0,0}, 2246cf0794eSJed Brown {1./4,7./4,3./4,-7./4,0}}, 2256cf0794eSJed Brown At[5][5] = {{0,0,0,0,0}, 2266cf0794eSJed Brown {0,1./2,0,0,0}, 2276cf0794eSJed Brown {0,1./6,1./2,0,0}, 2286cf0794eSJed Brown {0,-1./2,1./2,1./2,0}, 229108c343cSJed Brown {0,3./2,-3./2,1./2,1./2}}, 230108c343cSJed Brown *bembedt = PETSC_NULL; 231108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEXARS443,3,5,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,bembedt,bembedt,0,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2326cf0794eSJed Brown } 2336cf0794eSJed Brown { 2346cf0794eSJed Brown const PetscReal 2356cf0794eSJed Brown A[5][5] = {{0,0,0,0}, 2366cf0794eSJed Brown {1,0,0,0,0}, 2376cf0794eSJed Brown {4./9,2./9,0,0,0}, 2386cf0794eSJed Brown {1./4,0,3./4,0,0}, 2396cf0794eSJed Brown {1./4,0,3./5,0,0}}, 2406cf0794eSJed Brown At[5][5] = {{0,0,0,0}, 2416cf0794eSJed Brown {.5,.5,0,0,0}, 2426cf0794eSJed Brown {5./18,-1./9,.5,0,0}, 2436cf0794eSJed Brown {.5,0,0,.5,0}, 244108c343cSJed Brown {.25,0,.75,-.5,.5}}, 245108c343cSJed Brown *bembedt = PETSC_NULL; 246108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEXBPR3,3,5,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,bembedt,bembedt,0,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2476cf0794eSJed Brown } 2486cf0794eSJed Brown { 2496cf0794eSJed Brown const PetscReal 250a3a57f36SJed Brown A[6][6] = {{0,0,0,0,0,0}, 251a3a57f36SJed Brown {1./2,0,0,0,0,0}, 2524040e9f2SJed Brown {13861./62500.,6889./62500.,0,0,0,0}, 2534040e9f2SJed Brown {-116923316275./2393684061468.,-2731218467317./15368042101831.,9408046702089./11113171139209.,0,0,0}, 2544040e9f2SJed Brown {-451086348788./2902428689909.,-2682348792572./7519795681897.,12662868775082./11960479115383.,3355817975965./11060851509271.,0,0}, 2554040e9f2SJed Brown {647845179188./3216320057751.,73281519250./8382639484533.,552539513391./3454668386233.,3354512671639./8306763924573.,4040./17871.,0}}, 256a3a57f36SJed Brown At[6][6] = {{0,0,0,0,0,0}, 257a3a57f36SJed Brown {1./4,1./4,0,0,0,0}, 2584040e9f2SJed Brown {8611./62500.,-1743./31250.,1./4,0,0,0}, 2594040e9f2SJed Brown {5012029./34652500.,-654441./2922500.,174375./388108.,1./4,0,0}, 2604040e9f2SJed Brown {15267082809./155376265600.,-71443401./120774400.,730878875./902184768.,2285395./8070912.,1./4,0}, 2614040e9f2SJed Brown {82889./524892.,0,15625./83664.,69875./102672.,-2260./8211,1./4}}, 262*cc46b9d1SJed Brown bembedt[6] = {4586570599./29645900160.,0,178811875./945068544.,814220225./1159782912.,-3700637./11593932.,61727./225920.}, 2634040e9f2SJed Brown binterpt[6][3] = {{6943876665148./7220017795957.,-54480133./30881146.,6818779379841./7100303317025.}, 264cd652676SJed Brown {0,0,0}, 2654040e9f2SJed Brown {7640104374378./9702883013639.,-11436875./14766696.,2173542590792./12501825683035.}, 2664040e9f2SJed Brown {-20649996744609./7521556579894.,174696575./18121608.,-31592104683404./5083833661969.}, 2674040e9f2SJed Brown {8854892464581./2390941311638.,-12120380./966161.,61146701046299./7138195549469.}, 2684040e9f2SJed Brown {-11397109935349./6675773540249.,3843./706.,-17219254887155./4939391667607.}}; 269108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEX4,4,6,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,bembedt,bembedt,3,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 270a3a57f36SJed Brown } 271a3a57f36SJed Brown { 272a3a57f36SJed Brown const PetscReal 273a3a57f36SJed Brown A[8][8] = {{0,0,0,0,0,0,0,0}, 274a3a57f36SJed Brown {41./100,0,0,0,0,0,0,0}, 2754040e9f2SJed Brown {367902744464./2072280473677.,677623207551./8224143866563.,0,0,0,0,0,0}, 2764040e9f2SJed Brown {1268023523408./10340822734521.,0,1029933939417./13636558850479.,0,0,0,0,0}, 2774040e9f2SJed Brown {14463281900351./6315353703477.,0,66114435211212./5879490589093.,-54053170152839./4284798021562.,0,0,0,0}, 2784040e9f2SJed Brown {14090043504691./34967701212078.,0,15191511035443./11219624916014.,-18461159152457./12425892160975.,-281667163811./9011619295870.,0,0,0}, 2794040e9f2SJed Brown {19230459214898./13134317526959.,0,21275331358303./2942455364971.,-38145345988419./4862620318723.,-1./8,-1./8,0,0}, 2804040e9f2SJed Brown {-19977161125411./11928030595625.,0,-40795976796054./6384907823539.,177454434618887./12078138498510.,782672205425./8267701900261.,-69563011059811./9646580694205.,7356628210526./4942186776405.,0}}, 281a3a57f36SJed Brown At[8][8] = {{0,0,0,0,0,0,0,0}, 2824040e9f2SJed Brown {41./200.,41./200.,0,0,0,0,0,0}, 2834040e9f2SJed Brown {41./400.,-567603406766./11931857230679.,41./200.,0,0,0,0,0}, 2844040e9f2SJed Brown {683785636431./9252920307686.,0,-110385047103./1367015193373.,41./200.,0,0,0,0}, 2854040e9f2SJed Brown {3016520224154./10081342136671.,0,30586259806659./12414158314087.,-22760509404356./11113319521817.,41./200.,0,0,0}, 2864040e9f2SJed Brown {218866479029./1489978393911.,0,638256894668./5436446318841.,-1179710474555./5321154724896.,-60928119172./8023461067671.,41./200.,0,0}, 2874040e9f2SJed Brown {1020004230633./5715676835656.,0,25762820946817./25263940353407.,-2161375909145./9755907335909.,-211217309593./5846859502534.,-4269925059573./7827059040749.,41./200,0}, 2884040e9f2SJed Brown {-872700587467./9133579230613.,0,0,22348218063261./9555858737531.,-1143369518992./8141816002931.,-39379526789629./19018526304540.,32727382324388./42900044865799.,41./200.}}, 289*cc46b9d1SJed Brown bembedt[8] = {-975461918565./9796059967033.,0,0,78070527104295./32432590147079.,-548382580838./3424219808633.,-33438840321285./15594753105479.,3629800801594./4656183773603.,4035322873751./18575991585200.}, 2904040e9f2SJed Brown binterpt[8][3] = {{-17674230611817./10670229744614. , 43486358583215./12773830924787. , -9257016797708./5021505065439.}, 291cd652676SJed Brown {0 , 0 , 0 }, 292cd652676SJed Brown {0 , 0 , 0 }, 2934040e9f2SJed Brown {65168852399939./7868540260826. , -91478233927265./11067650958493., 26096422576131./11239449250142.}, 2944040e9f2SJed Brown {15494834004392./5936557850923. , -79368583304911./10890268929626., 92396832856987./20362823103730.}, 2954040e9f2SJed Brown {-99329723586156./26959484932159., -12239297817655./9152339842473. , 30029262896817./10175596800299.}, 2964040e9f2SJed Brown {-19024464361622./5461577185407. , 115839755401235./10719374521269., -26136350496073./3983972220547.}, 2974040e9f2SJed Brown {-6511271360970./6095937251113. , 5843115559534./2180450260947. , -5289405421727./3760307252460. }}; 298108c343cSJed Brown ierr = TSARKIMEXRegister(TSARKIMEX5,5,8,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL,bembedt,bembedt,3,binterpt[0],PETSC_NULL);CHKERRQ(ierr); 299a3a57f36SJed Brown } 300a3a57f36SJed Brown 3018a381b04SJed Brown PetscFunctionReturn(0); 3028a381b04SJed Brown } 3038a381b04SJed Brown 3048a381b04SJed Brown #undef __FUNCT__ 3058a381b04SJed Brown #define __FUNCT__ "TSARKIMEXRegisterDestroy" 3068a381b04SJed Brown /*@C 3078a381b04SJed Brown TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by TSARKIMEXRegister(). 3088a381b04SJed Brown 3098a381b04SJed Brown Not Collective 3108a381b04SJed Brown 3118a381b04SJed Brown Level: advanced 3128a381b04SJed Brown 3138a381b04SJed Brown .keywords: TSARKIMEX, register, destroy 3148a381b04SJed Brown .seealso: TSARKIMEXRegister(), TSARKIMEXRegisterAll(), TSARKIMEXRegisterDynamic() 3158a381b04SJed Brown @*/ 3168a381b04SJed Brown PetscErrorCode TSARKIMEXRegisterDestroy(void) 3178a381b04SJed Brown { 3188a381b04SJed Brown PetscErrorCode ierr; 3198a381b04SJed Brown ARKTableauLink link; 3208a381b04SJed Brown 3218a381b04SJed Brown PetscFunctionBegin; 3228a381b04SJed Brown while ((link = ARKTableauList)) { 3238a381b04SJed Brown ARKTableau t = &link->tab; 3248a381b04SJed Brown ARKTableauList = link->next; 3258a381b04SJed Brown ierr = PetscFree6(t->At,t->bt,t->ct,t->A,t->b,t->c);CHKERRQ(ierr); 326108c343cSJed Brown ierr = PetscFree2(t->bembedt,t->bembed);CHKERRQ(ierr); 327cd652676SJed Brown ierr = PetscFree2(t->binterpt,t->binterp);CHKERRQ(ierr); 3288a381b04SJed Brown ierr = PetscFree(t->name);CHKERRQ(ierr); 3298a381b04SJed Brown ierr = PetscFree(link);CHKERRQ(ierr); 3308a381b04SJed Brown } 3318a381b04SJed Brown TSARKIMEXRegisterAllCalled = PETSC_FALSE; 3328a381b04SJed Brown PetscFunctionReturn(0); 3338a381b04SJed Brown } 3348a381b04SJed Brown 3358a381b04SJed Brown #undef __FUNCT__ 3368a381b04SJed Brown #define __FUNCT__ "TSARKIMEXInitializePackage" 3378a381b04SJed Brown /*@C 3388a381b04SJed Brown TSARKIMEXInitializePackage - This function initializes everything in the TSARKIMEX package. It is called 3398a381b04SJed Brown from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_ARKIMEX() 3408a381b04SJed Brown when using static libraries. 3418a381b04SJed Brown 3428a381b04SJed Brown Input Parameter: 3438a381b04SJed Brown path - The dynamic library path, or PETSC_NULL 3448a381b04SJed Brown 3458a381b04SJed Brown Level: developer 3468a381b04SJed Brown 3478a381b04SJed Brown .keywords: TS, TSARKIMEX, initialize, package 3488a381b04SJed Brown .seealso: PetscInitialize() 3498a381b04SJed Brown @*/ 3508a381b04SJed Brown PetscErrorCode TSARKIMEXInitializePackage(const char path[]) 3518a381b04SJed Brown { 3528a381b04SJed Brown PetscErrorCode ierr; 3538a381b04SJed Brown 3548a381b04SJed Brown PetscFunctionBegin; 3558a381b04SJed Brown if (TSARKIMEXPackageInitialized) PetscFunctionReturn(0); 3568a381b04SJed Brown TSARKIMEXPackageInitialized = PETSC_TRUE; 3578a381b04SJed Brown ierr = TSARKIMEXRegisterAll();CHKERRQ(ierr); 3588a381b04SJed Brown ierr = PetscRegisterFinalize(TSARKIMEXFinalizePackage);CHKERRQ(ierr); 3598a381b04SJed Brown PetscFunctionReturn(0); 3608a381b04SJed Brown } 3618a381b04SJed Brown 3628a381b04SJed Brown #undef __FUNCT__ 3638a381b04SJed Brown #define __FUNCT__ "TSARKIMEXFinalizePackage" 3648a381b04SJed Brown /*@C 3658a381b04SJed Brown TSARKIMEXFinalizePackage - This function destroys everything in the TSARKIMEX package. It is 3668a381b04SJed Brown called from PetscFinalize(). 3678a381b04SJed Brown 3688a381b04SJed Brown Level: developer 3698a381b04SJed Brown 3708a381b04SJed Brown .keywords: Petsc, destroy, package 3718a381b04SJed Brown .seealso: PetscFinalize() 3728a381b04SJed Brown @*/ 3738a381b04SJed Brown PetscErrorCode TSARKIMEXFinalizePackage(void) 3748a381b04SJed Brown { 3758a381b04SJed Brown PetscErrorCode ierr; 3768a381b04SJed Brown 3778a381b04SJed Brown PetscFunctionBegin; 3788a381b04SJed Brown TSARKIMEXPackageInitialized = PETSC_FALSE; 3798a381b04SJed Brown ierr = TSARKIMEXRegisterDestroy();CHKERRQ(ierr); 3808a381b04SJed Brown PetscFunctionReturn(0); 3818a381b04SJed Brown } 3828a381b04SJed Brown 3838a381b04SJed Brown #undef __FUNCT__ 3848a381b04SJed Brown #define __FUNCT__ "TSARKIMEXRegister" 385cd652676SJed Brown /*@C 386cd652676SJed Brown TSARKIMEXRegister - register an ARK IMEX scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 387cd652676SJed Brown 388cd652676SJed Brown Not Collective, but the same schemes should be registered on all processes on which they will be used 389cd652676SJed Brown 390cd652676SJed Brown Input Parameters: 391cd652676SJed Brown + name - identifier for method 392cd652676SJed Brown . order - approximation order of method 393cd652676SJed Brown . s - number of stages, this is the dimension of the matrices below 394cd652676SJed Brown . At - Butcher table of stage coefficients for stiff part (dimension s*s, row-major) 395cd652676SJed Brown . bt - Butcher table for completing the stiff part of the step (dimension s; PETSC_NULL to use the last row of At) 396cd652676SJed Brown . ct - Abscissa of each stiff stage (dimension s, PETSC_NULL to use row sums of At) 397cd652676SJed Brown . A - Non-stiff stage coefficients (dimension s*s, row-major) 398cd652676SJed Brown . b - Non-stiff step completion table (dimension s; PETSC_NULL to use last row of At) 399cd652676SJed Brown . c - Non-stiff abscissa (dimension s; PETSC_NULL to use row sums of A) 400108c343cSJed Brown . bembedt - Stiff part of completion table for embedded method (dimension s; PETSC_NULL if not available) 401108c343cSJed Brown . bembed - Non-stiff part of completion table for embedded method (dimension s; PETSC_NULL to use bembedt if provided) 402cd652676SJed Brown . pinterp - Order of the interpolation scheme, equal to the number of columns of binterpt and binterp 403cd652676SJed Brown . binterpt - Coefficients of the interpolation formula for the stiff part (dimension s*pinterp) 404cd652676SJed Brown - binterp - Coefficients of the interpolation formula for the non-stiff part (dimension s*pinterp; PETSC_NULL to reuse binterpt) 405cd652676SJed Brown 406cd652676SJed Brown Notes: 407cd652676SJed Brown Several ARK IMEX methods are provided, this function is only needed to create new methods. 408cd652676SJed Brown 409cd652676SJed Brown Level: advanced 410cd652676SJed Brown 411cd652676SJed Brown .keywords: TS, register 412cd652676SJed Brown 413cd652676SJed Brown .seealso: TSARKIMEX 414cd652676SJed Brown @*/ 4158a381b04SJed Brown PetscErrorCode TSARKIMEXRegister(const TSARKIMEXType name,PetscInt order,PetscInt s, 4168a381b04SJed Brown const PetscReal At[],const PetscReal bt[],const PetscReal ct[], 417cd652676SJed Brown const PetscReal A[],const PetscReal b[],const PetscReal c[], 418108c343cSJed Brown const PetscReal bembedt[],const PetscReal bembed[], 419cd652676SJed Brown PetscInt pinterp,const PetscReal binterpt[],const PetscReal binterp[]) 4208a381b04SJed Brown { 4218a381b04SJed Brown PetscErrorCode ierr; 4228a381b04SJed Brown ARKTableauLink link; 4238a381b04SJed Brown ARKTableau t; 4248a381b04SJed Brown PetscInt i,j; 4258a381b04SJed Brown 4268a381b04SJed Brown PetscFunctionBegin; 4278a381b04SJed Brown ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 428cd652676SJed Brown ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 4298a381b04SJed Brown t = &link->tab; 4308a381b04SJed Brown ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 4318a381b04SJed Brown t->order = order; 4328a381b04SJed Brown t->s = s; 4338a381b04SJed 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); 4348a381b04SJed Brown ierr = PetscMemcpy(t->At,At,s*s*sizeof(At[0]));CHKERRQ(ierr); 4358a381b04SJed Brown ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 4368a381b04SJed Brown if (bt) {ierr = PetscMemcpy(t->bt,bt,s*sizeof(bt[0]));CHKERRQ(ierr);} 4378a381b04SJed Brown else for (i=0; i<s; i++) t->bt[i] = At[(s-1)*s+i]; 4388a381b04SJed Brown if (b) {ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr);} 4398a381b04SJed Brown else for (i=0; i<s; i++) t->b[i] = At[(s-1)*s+i]; 4408a381b04SJed Brown if (ct) {ierr = PetscMemcpy(t->ct,ct,s*sizeof(ct[0]));CHKERRQ(ierr);} 4418a381b04SJed 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]; 4428a381b04SJed Brown if (c) {ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr);} 4438a381b04SJed 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]; 444108c343cSJed Brown if (bembedt) { 445108c343cSJed Brown ierr = PetscMalloc2(s,PetscReal,&t->bembedt,s,PetscReal,&t->bembed);CHKERRQ(ierr); 446108c343cSJed Brown ierr = PetscMemcpy(t->bembedt,bembedt,s*sizeof(bembedt[0]));CHKERRQ(ierr); 447108c343cSJed Brown ierr = PetscMemcpy(t->bembed,bembed?bembed:bembedt,s*sizeof(bembed[0]));CHKERRQ(ierr); 448108c343cSJed Brown } 449108c343cSJed Brown 4504f385281SJed Brown t->pinterp = pinterp; 451cd652676SJed Brown ierr = PetscMalloc2(s*pinterp,PetscReal,&t->binterpt,s*pinterp,PetscReal,&t->binterp);CHKERRQ(ierr); 452cd652676SJed Brown ierr = PetscMemcpy(t->binterpt,binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 453cd652676SJed Brown ierr = PetscMemcpy(t->binterp,binterp?binterp:binterpt,s*pinterp*sizeof(binterpt[0]));CHKERRQ(ierr); 4548a381b04SJed Brown link->next = ARKTableauList; 4558a381b04SJed Brown ARKTableauList = link; 4568a381b04SJed Brown PetscFunctionReturn(0); 4578a381b04SJed Brown } 4588a381b04SJed Brown 4598a381b04SJed Brown #undef __FUNCT__ 460108c343cSJed Brown #define __FUNCT__ "TSEvaluateStep_ARKIMEX" 461108c343cSJed Brown /* 462108c343cSJed Brown The step completion formula is 463108c343cSJed Brown 464108c343cSJed Brown x1 = x0 - h bt^T YdotI + h b^T YdotRHS 465108c343cSJed Brown 466108c343cSJed Brown This function can be called before or after ts->vec_sol has been updated. 467108c343cSJed Brown Suppose we have a completion formula (bt,b) and an embedded formula (bet,be) of different order. 468108c343cSJed Brown We can write 469108c343cSJed Brown 470108c343cSJed Brown x1e = x0 - h bet^T YdotI + h be^T YdotRHS 471108c343cSJed Brown = x1 + h bt^T YdotI - h b^T YdotRHS - h bet^T YdotI + h be^T YdotRHS 472108c343cSJed Brown = x1 - h (bet - bt)^T YdotI + h (be - b)^T YdotRHS 473108c343cSJed Brown 474108c343cSJed Brown so we can evaluate the method with different order even after the step has been optimistically completed. 475108c343cSJed Brown */ 476108c343cSJed Brown static PetscErrorCode TSEvaluateStep_ARKIMEX(TS ts,PetscInt order,Vec X,PetscBool *done) 477108c343cSJed Brown { 478108c343cSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 479108c343cSJed Brown ARKTableau tab = ark->tableau; 480108c343cSJed Brown PetscScalar *w = ark->work; 481108c343cSJed Brown PetscReal h; 482108c343cSJed Brown PetscInt s = tab->s,j; 483108c343cSJed Brown PetscErrorCode ierr; 484108c343cSJed Brown 485108c343cSJed Brown PetscFunctionBegin; 486108c343cSJed Brown switch (ark->status) { 487108c343cSJed Brown case TS_STEP_INCOMPLETE: 488108c343cSJed Brown case TS_STEP_PENDING: 489108c343cSJed Brown h = ts->time_step; break; 490108c343cSJed Brown case TS_STEP_COMPLETE: 491108c343cSJed Brown h = ts->time_step_prev; break; 492b9ce6d65SJed Brown default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 493108c343cSJed Brown } 494108c343cSJed Brown if (order == tab->order) { 495108c343cSJed Brown if (ark->status == TS_STEP_INCOMPLETE) { /* Use the standard completion formula (bt,b) */ 496108c343cSJed Brown ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 497108c343cSJed Brown for (j=0; j<s; j++) w[j] = -h*tab->bt[j]; 498108c343cSJed Brown ierr = VecMAXPY(X,s,w,ark->YdotI);CHKERRQ(ierr); 499108c343cSJed Brown for (j=0; j<s; j++) w[j] = h*tab->b[j]; 500108c343cSJed Brown ierr = VecMAXPY(X,s,w,ark->YdotRHS);CHKERRQ(ierr); 501108c343cSJed Brown } else {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 502108c343cSJed Brown if (done) *done = PETSC_TRUE; 503108c343cSJed Brown PetscFunctionReturn(0); 504108c343cSJed Brown } else if (order == tab->order-1) { 505108c343cSJed Brown if (!tab->bembedt) goto unavailable; 506108c343cSJed Brown if (ark->status == TS_STEP_INCOMPLETE) { /* Complete with the embedded method (bet,be) */ 507108c343cSJed Brown ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 508108c343cSJed Brown for (j=0; j<s; j++) w[j] = -h*tab->bembedt[j]; 509108c343cSJed Brown ierr = VecMAXPY(X,s,w,ark->YdotI);CHKERRQ(ierr); 510108c343cSJed Brown for (j=0; j<s; j++) w[j] = h*tab->bembed[j]; 511108c343cSJed Brown ierr = VecMAXPY(X,s,w,ark->YdotRHS);CHKERRQ(ierr); 512108c343cSJed Brown } else { /* Rollback and re-complete using (bet-be,be-b) */ 513108c343cSJed Brown ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 514108c343cSJed Brown for (j=0; j<s; j++) w[j] = -h*(tab->bembedt[j] - tab->bt[j]); 515108c343cSJed Brown ierr = VecMAXPY(X,tab->s,w,ark->YdotI);CHKERRQ(ierr); 516108c343cSJed Brown for (j=0; j<s; j++) w[j] = h*(tab->bembed[j] - tab->b[j]); 517108c343cSJed Brown ierr = VecMAXPY(X,s,w,ark->YdotRHS);CHKERRQ(ierr); 518108c343cSJed Brown } 519108c343cSJed Brown if (done) *done = PETSC_TRUE; 520108c343cSJed Brown PetscFunctionReturn(0); 521108c343cSJed Brown } 522108c343cSJed Brown unavailable: 523108c343cSJed Brown if (done) *done = PETSC_FALSE; 524108c343cSJed Brown else SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_SUP,"ARKIMEX '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 525108c343cSJed Brown PetscFunctionReturn(0); 526108c343cSJed Brown } 527108c343cSJed Brown 528108c343cSJed Brown #undef __FUNCT__ 5298a381b04SJed Brown #define __FUNCT__ "TSStep_ARKIMEX" 5308a381b04SJed Brown static PetscErrorCode TSStep_ARKIMEX(TS ts) 5318a381b04SJed Brown { 5328a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 5338a381b04SJed Brown ARKTableau tab = ark->tableau; 5348a381b04SJed Brown const PetscInt s = tab->s; 5358a381b04SJed Brown const PetscReal *At = tab->At,*A = tab->A,*bt = tab->bt,*b = tab->b,*ct = tab->ct,*c = tab->c; 536406d0ec2SJed Brown PetscScalar *w = ark->work; 5378a381b04SJed Brown Vec *Y = ark->Y,*YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS,Ydot = ark->Ydot,W = ark->Work,Z = ark->Z; 538108c343cSJed Brown TSAdapt adapt; 5398a381b04SJed Brown SNES snes; 540108c343cSJed Brown PetscInt i,j,its,lits,reject,next_scheme; 541cdbf8f93SLisandro Dalcin PetscReal next_time_step; 542108c343cSJed Brown PetscReal t; 543108c343cSJed Brown PetscBool accept; 5448a381b04SJed Brown PetscErrorCode ierr; 5458a381b04SJed Brown 5468a381b04SJed Brown PetscFunctionBegin; 5478a381b04SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 548cdbf8f93SLisandro Dalcin next_time_step = ts->time_step; 5498a381b04SJed Brown t = ts->ptime; 550108c343cSJed Brown accept = PETSC_TRUE; 551108c343cSJed Brown ark->status = TS_STEP_INCOMPLETE; 5528a381b04SJed Brown 55397335746SJed Brown for (reject=0; reject<ts->max_reject && !ts->reason; reject++,ts->reject++) { 554108c343cSJed Brown PetscReal h = ts->time_step; 5558a381b04SJed Brown for (i=0; i<s; i++) { 5568a381b04SJed Brown if (At[i*s+i] == 0) { /* This stage is explicit */ 5578a381b04SJed Brown ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 5588a381b04SJed Brown for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 5598a381b04SJed Brown ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 5608a381b04SJed Brown for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 5618a381b04SJed Brown ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 5628a381b04SJed Brown } else { 5638a381b04SJed Brown ark->stage_time = t + h*ct[i]; 5648a381b04SJed Brown ark->shift = 1./(h*At[i*s+i]); 5658a381b04SJed Brown /* Affine part */ 5668a381b04SJed Brown ierr = VecZeroEntries(W);CHKERRQ(ierr); 5678a381b04SJed Brown for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 5688a381b04SJed Brown ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 569f16577ceSEmil Constantinescu ierr = VecScale(W, ark->shift);CHKERRQ(ierr); 570f16577ceSEmil Constantinescu 5718a381b04SJed Brown /* Ydot = shift*(Y-Z) */ 5728a381b04SJed Brown ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 5734f385281SJed Brown for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 5744f385281SJed Brown ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr); 575f16577ceSEmil Constantinescu 5768a381b04SJed Brown /* Initial guess taken from last stage */ 5778a381b04SJed Brown ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr); 5788a381b04SJed Brown ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 5798a381b04SJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 5808a381b04SJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 5818a381b04SJed Brown ts->nonlinear_its += its; ts->linear_its += lits; 58297335746SJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 58397335746SJed Brown ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 58497335746SJed Brown if (!accept) goto reject_step; 5858a381b04SJed Brown } 5868a381b04SJed Brown ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 5874cc180ffSJed Brown ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr); 5884cc180ffSJed Brown if (ark->imex) { 5898a381b04SJed Brown ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 5904cc180ffSJed Brown } else { 5914cc180ffSJed Brown ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr); 5924cc180ffSJed Brown } 5938a381b04SJed Brown } 594108c343cSJed Brown ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 595108c343cSJed Brown ark->status = TS_STEP_PENDING; 5968a381b04SJed Brown 597108c343cSJed Brown /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 598108c343cSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 599108c343cSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 600108c343cSJed Brown ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 601108c343cSJed Brown ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 602108c343cSJed Brown if (accept) { 603108c343cSJed Brown /* ignore next_scheme for now */ 6048a381b04SJed Brown ts->ptime += ts->time_step; 605cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 6068a381b04SJed Brown ts->steps++; 607108c343cSJed Brown ark->status = TS_STEP_COMPLETE; 608108c343cSJed Brown break; 609108c343cSJed Brown } else { /* Roll back the current step */ 610108c343cSJed Brown for (j=0; j<s; j++) w[j] = h*bt[j]; 611108c343cSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotI);CHKERRQ(ierr); 612108c343cSJed Brown for (j=0; j<s; j++) w[j] = -h*b[j]; 613108c343cSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotRHS);CHKERRQ(ierr); 614108c343cSJed Brown ts->time_step = next_time_step; 615108c343cSJed Brown ark->status = TS_STEP_INCOMPLETE; 616108c343cSJed Brown } 617476b6736SJed Brown reject_step: continue; 618108c343cSJed Brown } 619b2ce242eSJed Brown if (ark->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 6208a381b04SJed Brown PetscFunctionReturn(0); 6218a381b04SJed Brown } 6228a381b04SJed Brown 623cd652676SJed Brown #undef __FUNCT__ 624cd652676SJed Brown #define __FUNCT__ "TSInterpolate_ARKIMEX" 625cd652676SJed Brown static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X) 626cd652676SJed Brown { 627cd652676SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 6284f385281SJed Brown PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j; 629108c343cSJed Brown PetscReal h; 630108c343cSJed Brown PetscReal tt,t; 631cd652676SJed Brown PetscScalar *bt,*b; 632cd652676SJed Brown const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp; 633cd652676SJed Brown PetscErrorCode ierr; 634cd652676SJed Brown 635cd652676SJed Brown PetscFunctionBegin; 636cd652676SJed Brown if (!Bt || !B) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name); 637108c343cSJed Brown switch (ark->status) { 638108c343cSJed Brown case TS_STEP_INCOMPLETE: 639108c343cSJed Brown case TS_STEP_PENDING: 640108c343cSJed Brown h = ts->time_step; 641108c343cSJed Brown t = (itime - ts->ptime)/h; 642108c343cSJed Brown break; 643108c343cSJed Brown case TS_STEP_COMPLETE: 644108c343cSJed Brown h = ts->time_step_prev; 645108c343cSJed Brown t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 646108c343cSJed Brown break; 647b9ce6d65SJed Brown default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 648108c343cSJed Brown } 649cd652676SJed Brown ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr); 650cd652676SJed Brown for (i=0; i<s; i++) bt[i] = b[i] = 0; 6514f385281SJed Brown for (j=0,tt=t; j<pinterp; j++,tt*=t) { 652cd652676SJed Brown for (i=0; i<s; i++) { 653108c343cSJed Brown bt[i] += h * Bt[i*pinterp+j] * tt * -1.0; 654108c343cSJed Brown b[i] += h * B[i*pinterp+j] * tt; 655cd652676SJed Brown } 656cd652676SJed Brown } 657cd652676SJed 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"); 658cd652676SJed Brown ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr); 659cd652676SJed Brown ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr); 660cd652676SJed Brown ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr); 661cd652676SJed Brown ierr = PetscFree2(bt,b);CHKERRQ(ierr); 662cd652676SJed Brown PetscFunctionReturn(0); 663cd652676SJed Brown } 664cd652676SJed Brown 6658a381b04SJed Brown /*------------------------------------------------------------*/ 6668a381b04SJed Brown #undef __FUNCT__ 6678a381b04SJed Brown #define __FUNCT__ "TSReset_ARKIMEX" 6688a381b04SJed Brown static PetscErrorCode TSReset_ARKIMEX(TS ts) 6698a381b04SJed Brown { 6708a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 6718a381b04SJed Brown PetscInt s; 6728a381b04SJed Brown PetscErrorCode ierr; 6738a381b04SJed Brown 6748a381b04SJed Brown PetscFunctionBegin; 6758a381b04SJed Brown if (!ark->tableau) PetscFunctionReturn(0); 6768a381b04SJed Brown s = ark->tableau->s; 6778a381b04SJed Brown ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 6788a381b04SJed Brown ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 6798a381b04SJed Brown ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 6808a381b04SJed Brown ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 6818a381b04SJed Brown ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 6828a381b04SJed Brown ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 6838a381b04SJed Brown ierr = PetscFree(ark->work);CHKERRQ(ierr); 6848a381b04SJed Brown PetscFunctionReturn(0); 6858a381b04SJed Brown } 6868a381b04SJed Brown 6878a381b04SJed Brown #undef __FUNCT__ 6888a381b04SJed Brown #define __FUNCT__ "TSDestroy_ARKIMEX" 6898a381b04SJed Brown static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 6908a381b04SJed Brown { 6918a381b04SJed Brown PetscErrorCode ierr; 6928a381b04SJed Brown 6938a381b04SJed Brown PetscFunctionBegin; 6948a381b04SJed Brown ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 6958a381b04SJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 6968a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","",PETSC_NULL);CHKERRQ(ierr); 6978a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","",PETSC_NULL);CHKERRQ(ierr); 698995b3938SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","",PETSC_NULL);CHKERRQ(ierr); 6998a381b04SJed Brown PetscFunctionReturn(0); 7008a381b04SJed Brown } 7018a381b04SJed Brown 7028a381b04SJed Brown /* 7038a381b04SJed Brown This defines the nonlinear equation that is to be solved with SNES 7048a381b04SJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 7058a381b04SJed Brown */ 7068a381b04SJed Brown #undef __FUNCT__ 7078a381b04SJed Brown #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 7088a381b04SJed Brown static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 7098a381b04SJed Brown { 7108a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7118a381b04SJed Brown PetscErrorCode ierr; 7128a381b04SJed Brown 7138a381b04SJed Brown PetscFunctionBegin; 7148a381b04SJed Brown ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 7154cc180ffSJed Brown ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,F,ark->imex);CHKERRQ(ierr); 7168a381b04SJed Brown PetscFunctionReturn(0); 7178a381b04SJed Brown } 7188a381b04SJed Brown 7198a381b04SJed Brown #undef __FUNCT__ 7208a381b04SJed Brown #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 7218a381b04SJed Brown static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 7228a381b04SJed Brown { 7238a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7248a381b04SJed Brown PetscErrorCode ierr; 7258a381b04SJed Brown 7268a381b04SJed Brown PetscFunctionBegin; 7278a381b04SJed Brown /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 7288a381b04SJed Brown ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 7298a381b04SJed Brown PetscFunctionReturn(0); 7308a381b04SJed Brown } 7318a381b04SJed Brown 7328a381b04SJed Brown #undef __FUNCT__ 7338a381b04SJed Brown #define __FUNCT__ "TSSetUp_ARKIMEX" 7348a381b04SJed Brown static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 7358a381b04SJed Brown { 7368a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7378a381b04SJed Brown ARKTableau tab = ark->tableau; 7388a381b04SJed Brown PetscInt s = tab->s; 7398a381b04SJed Brown PetscErrorCode ierr; 7408a381b04SJed Brown 7418a381b04SJed Brown PetscFunctionBegin; 7428a381b04SJed Brown if (!ark->tableau) { 743e24355feSJed Brown ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 7448a381b04SJed Brown } 7458a381b04SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 7468a381b04SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 7478a381b04SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 7488a381b04SJed Brown ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 7498a381b04SJed Brown ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 7508a381b04SJed Brown ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 7518a381b04SJed Brown ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 7528a381b04SJed Brown PetscFunctionReturn(0); 7538a381b04SJed Brown } 7548a381b04SJed Brown /*------------------------------------------------------------*/ 7558a381b04SJed Brown 7568a381b04SJed Brown #undef __FUNCT__ 7578a381b04SJed Brown #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 7588a381b04SJed Brown static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 7598a381b04SJed Brown { 7604cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7618a381b04SJed Brown PetscErrorCode ierr; 7628a381b04SJed Brown char arktype[256]; 7638a381b04SJed Brown 7648a381b04SJed Brown PetscFunctionBegin; 7658a381b04SJed Brown ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 7668a381b04SJed Brown { 7678a381b04SJed Brown ARKTableauLink link; 7688a381b04SJed Brown PetscInt count,choice; 7698a381b04SJed Brown PetscBool flg; 7708a381b04SJed Brown const char **namelist; 7718a381b04SJed Brown ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof arktype);CHKERRQ(ierr); 7728a381b04SJed Brown for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 7738a381b04SJed Brown ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 7748a381b04SJed Brown for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 7758a381b04SJed Brown ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 7768a381b04SJed Brown ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 7778a381b04SJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 7784cc180ffSJed Brown flg = (PetscBool)!ark->imex; 7794cc180ffSJed Brown ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 7804cc180ffSJed Brown ark->imex = (PetscBool)!flg; 781d52bd9f3SBarry Smith ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 7828a381b04SJed Brown } 7838a381b04SJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 7848a381b04SJed Brown PetscFunctionReturn(0); 7858a381b04SJed Brown } 7868a381b04SJed Brown 7878a381b04SJed Brown #undef __FUNCT__ 7888a381b04SJed Brown #define __FUNCT__ "PetscFormatRealArray" 7898a381b04SJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 7908a381b04SJed Brown { 791257d2499SJed Brown PetscErrorCode ierr; 792f1d86077SJed Brown PetscInt i; 793f1d86077SJed Brown size_t left,count; 7948a381b04SJed Brown char *p; 7958a381b04SJed Brown 7968a381b04SJed Brown PetscFunctionBegin; 797f1d86077SJed Brown for (i=0,p=buf,left=len; i<n; i++) { 798f1d86077SJed Brown ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 7998a381b04SJed Brown if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 8008a381b04SJed Brown left -= count; 8018a381b04SJed Brown p += count; 8028a381b04SJed Brown *p++ = ' '; 8038a381b04SJed Brown } 8048a381b04SJed Brown p[i ? 0 : -1] = 0; 8058a381b04SJed Brown PetscFunctionReturn(0); 8068a381b04SJed Brown } 8078a381b04SJed Brown 8088a381b04SJed Brown #undef __FUNCT__ 8098a381b04SJed Brown #define __FUNCT__ "TSView_ARKIMEX" 8108a381b04SJed Brown static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 8118a381b04SJed Brown { 8128a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 8138a381b04SJed Brown ARKTableau tab = ark->tableau; 8148a381b04SJed Brown PetscBool iascii; 8158a381b04SJed Brown PetscErrorCode ierr; 8168a381b04SJed Brown 8178a381b04SJed Brown PetscFunctionBegin; 8188a381b04SJed Brown ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 8198a381b04SJed Brown if (iascii) { 8208a381b04SJed Brown const TSARKIMEXType arktype; 8218a381b04SJed Brown char buf[512]; 8228a381b04SJed Brown ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 8238a381b04SJed Brown ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 8248a381b04SJed Brown ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 82531f6fcc0SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 82631f6fcc0SJed Brown ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 82731f6fcc0SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 8288a381b04SJed Brown } 829d52bd9f3SBarry Smith ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 8308a381b04SJed Brown PetscFunctionReturn(0); 8318a381b04SJed Brown } 8328a381b04SJed Brown 8338a381b04SJed Brown #undef __FUNCT__ 8348a381b04SJed Brown #define __FUNCT__ "TSARKIMEXSetType" 8358a381b04SJed Brown /*@C 8368a381b04SJed Brown TSARKIMEXSetType - Set the type of ARK IMEX scheme 8378a381b04SJed Brown 8388a381b04SJed Brown Logically collective 8398a381b04SJed Brown 8408a381b04SJed Brown Input Parameter: 8418a381b04SJed Brown + ts - timestepping context 8428a381b04SJed Brown - arktype - type of ARK-IMEX scheme 8438a381b04SJed Brown 8448a381b04SJed Brown Level: intermediate 8458a381b04SJed Brown 846020d8f30SJed Brown .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2, TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5 8478a381b04SJed Brown @*/ 8488a381b04SJed Brown PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype) 8498a381b04SJed Brown { 8508a381b04SJed Brown PetscErrorCode ierr; 8518a381b04SJed Brown 8528a381b04SJed Brown PetscFunctionBegin; 8538a381b04SJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 8548a381b04SJed Brown ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 8558a381b04SJed Brown PetscFunctionReturn(0); 8568a381b04SJed Brown } 8578a381b04SJed Brown 8588a381b04SJed Brown #undef __FUNCT__ 8598a381b04SJed Brown #define __FUNCT__ "TSARKIMEXGetType" 8608a381b04SJed Brown /*@C 8618a381b04SJed Brown TSARKIMEXGetType - Get the type of ARK IMEX scheme 8628a381b04SJed Brown 8638a381b04SJed Brown Logically collective 8648a381b04SJed Brown 8658a381b04SJed Brown Input Parameter: 8668a381b04SJed Brown . ts - timestepping context 8678a381b04SJed Brown 8688a381b04SJed Brown Output Parameter: 8698a381b04SJed Brown . arktype - type of ARK-IMEX scheme 8708a381b04SJed Brown 8718a381b04SJed Brown Level: intermediate 8728a381b04SJed Brown 8738a381b04SJed Brown .seealso: TSARKIMEXGetType() 8748a381b04SJed Brown @*/ 8758a381b04SJed Brown PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype) 8768a381b04SJed Brown { 8778a381b04SJed Brown PetscErrorCode ierr; 8788a381b04SJed Brown 8798a381b04SJed Brown PetscFunctionBegin; 8808a381b04SJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 8818a381b04SJed Brown ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 8828a381b04SJed Brown PetscFunctionReturn(0); 8838a381b04SJed Brown } 8848a381b04SJed Brown 8854cc180ffSJed Brown #undef __FUNCT__ 8864cc180ffSJed Brown #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 8874cc180ffSJed Brown /*@C 8884cc180ffSJed Brown TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 8894cc180ffSJed Brown 8904cc180ffSJed Brown Logically collective 8914cc180ffSJed Brown 8924cc180ffSJed Brown Input Parameter: 8934cc180ffSJed Brown + ts - timestepping context 8944cc180ffSJed Brown - flg - PETSC_TRUE for fully implicit 8954cc180ffSJed Brown 8964cc180ffSJed Brown Level: intermediate 8974cc180ffSJed Brown 8984cc180ffSJed Brown .seealso: TSARKIMEXGetType() 8994cc180ffSJed Brown @*/ 9004cc180ffSJed Brown PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 9014cc180ffSJed Brown { 9024cc180ffSJed Brown PetscErrorCode ierr; 9034cc180ffSJed Brown 9044cc180ffSJed Brown PetscFunctionBegin; 9054cc180ffSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 9064cc180ffSJed Brown ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 9074cc180ffSJed Brown PetscFunctionReturn(0); 9084cc180ffSJed Brown } 9094cc180ffSJed Brown 9108a381b04SJed Brown EXTERN_C_BEGIN 9118a381b04SJed Brown #undef __FUNCT__ 9128a381b04SJed Brown #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 9138a381b04SJed Brown PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype) 9148a381b04SJed Brown { 9158a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 9168a381b04SJed Brown PetscErrorCode ierr; 9178a381b04SJed Brown 9188a381b04SJed Brown PetscFunctionBegin; 9198a381b04SJed Brown if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 9208a381b04SJed Brown *arktype = ark->tableau->name; 9218a381b04SJed Brown PetscFunctionReturn(0); 9228a381b04SJed Brown } 9238a381b04SJed Brown #undef __FUNCT__ 9248a381b04SJed Brown #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 9258a381b04SJed Brown PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype) 9268a381b04SJed Brown { 9278a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 9288a381b04SJed Brown PetscErrorCode ierr; 9298a381b04SJed Brown PetscBool match; 9308a381b04SJed Brown ARKTableauLink link; 9318a381b04SJed Brown 9328a381b04SJed Brown PetscFunctionBegin; 9338a381b04SJed Brown if (ark->tableau) { 9348a381b04SJed Brown ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 9358a381b04SJed Brown if (match) PetscFunctionReturn(0); 9368a381b04SJed Brown } 9378a381b04SJed Brown for (link = ARKTableauList; link; link=link->next) { 9388a381b04SJed Brown ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 9398a381b04SJed Brown if (match) { 9408a381b04SJed Brown ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 9418a381b04SJed Brown ark->tableau = &link->tab; 9428a381b04SJed Brown PetscFunctionReturn(0); 9438a381b04SJed Brown } 9448a381b04SJed Brown } 9458a381b04SJed Brown SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 9468a381b04SJed Brown PetscFunctionReturn(0); 9478a381b04SJed Brown } 9484cc180ffSJed Brown #undef __FUNCT__ 9494cc180ffSJed Brown #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 9504cc180ffSJed Brown PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 9514cc180ffSJed Brown { 9524cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 9534cc180ffSJed Brown 9544cc180ffSJed Brown PetscFunctionBegin; 9554cc180ffSJed Brown ark->imex = (PetscBool)!flg; 9564cc180ffSJed Brown PetscFunctionReturn(0); 9574cc180ffSJed Brown } 9588a381b04SJed Brown EXTERN_C_END 9598a381b04SJed Brown 9608a381b04SJed Brown /* ------------------------------------------------------------ */ 9618a381b04SJed Brown /*MC 9628a381b04SJed Brown TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes 9638a381b04SJed Brown 964fca742c7SJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 965fca742c7SJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 966fca742c7SJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 967fca742c7SJed Brown 968fca742c7SJed Brown Notes: 969c8058688SBarry Smith The default is TSARKIMEX2E, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 970c8058688SBarry Smith 971fca742c7SJed Brown This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 972fca742c7SJed Brown 9738a381b04SJed Brown Level: beginner 9748a381b04SJed Brown 975c8058688SBarry Smith .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3, 97626885f59SBarry Smith TSARKIMEX4, TSARKIMEX5, TSARKIMEXType, TSARKIMEXRegister() 9778a381b04SJed Brown 9788a381b04SJed Brown M*/ 9798a381b04SJed Brown EXTERN_C_BEGIN 9808a381b04SJed Brown #undef __FUNCT__ 9818a381b04SJed Brown #define __FUNCT__ "TSCreate_ARKIMEX" 9828a381b04SJed Brown PetscErrorCode TSCreate_ARKIMEX(TS ts) 9838a381b04SJed Brown { 9848a381b04SJed Brown TS_ARKIMEX *th; 9858a381b04SJed Brown PetscErrorCode ierr; 9868a381b04SJed Brown 9878a381b04SJed Brown PetscFunctionBegin; 9888a381b04SJed Brown #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 9898a381b04SJed Brown ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 9908a381b04SJed Brown #endif 9918a381b04SJed Brown 9928a381b04SJed Brown ts->ops->reset = TSReset_ARKIMEX; 9938a381b04SJed Brown ts->ops->destroy = TSDestroy_ARKIMEX; 9948a381b04SJed Brown ts->ops->view = TSView_ARKIMEX; 9958a381b04SJed Brown ts->ops->setup = TSSetUp_ARKIMEX; 9968a381b04SJed Brown ts->ops->step = TSStep_ARKIMEX; 997cd652676SJed Brown ts->ops->interpolate = TSInterpolate_ARKIMEX; 998108c343cSJed Brown ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 9998a381b04SJed Brown ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 10008a381b04SJed Brown ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 10018a381b04SJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 10028a381b04SJed Brown 10038a381b04SJed Brown ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 10048a381b04SJed Brown ts->data = (void*)th; 10054cc180ffSJed Brown th->imex = PETSC_TRUE; 10068a381b04SJed Brown 10078a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 10088a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 10094cc180ffSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","TSARKIMEXSetFullyImplicit_ARKIMEX",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 10108a381b04SJed Brown PetscFunctionReturn(0); 10118a381b04SJed Brown } 10128a381b04SJed Brown EXTERN_C_END 1013