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; 540f1b97656SJed Brown SNESConvergedReason snesreason; 541108c343cSJed Brown PetscInt i,j,its,lits,reject,next_scheme; 542cdbf8f93SLisandro Dalcin PetscReal next_time_step; 543108c343cSJed Brown PetscReal t; 544108c343cSJed Brown PetscBool accept; 5458a381b04SJed Brown PetscErrorCode ierr; 5468a381b04SJed Brown 5478a381b04SJed Brown PetscFunctionBegin; 5488a381b04SJed Brown ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 549cdbf8f93SLisandro Dalcin next_time_step = ts->time_step; 5508a381b04SJed Brown t = ts->ptime; 551108c343cSJed Brown accept = PETSC_TRUE; 552108c343cSJed Brown ark->status = TS_STEP_INCOMPLETE; 5538a381b04SJed Brown 554108c343cSJed Brown for (reject=0; reject<ts->max_reject; reject++,ts->reject++) { 555108c343cSJed Brown PetscReal h = ts->time_step; 5568a381b04SJed Brown for (i=0; i<s; i++) { 5578a381b04SJed Brown if (At[i*s+i] == 0) { /* This stage is explicit */ 5588a381b04SJed Brown ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 5598a381b04SJed Brown for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 5608a381b04SJed Brown ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 5618a381b04SJed Brown for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 5628a381b04SJed Brown ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 5638a381b04SJed Brown } else { 5648a381b04SJed Brown ark->stage_time = t + h*ct[i]; 5658a381b04SJed Brown ark->shift = 1./(h*At[i*s+i]); 5668a381b04SJed Brown /* Affine part */ 5678a381b04SJed Brown ierr = VecZeroEntries(W);CHKERRQ(ierr); 5688a381b04SJed Brown for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 5698a381b04SJed Brown ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 570f16577ceSEmil Constantinescu ierr = VecScale(W, ark->shift);CHKERRQ(ierr); 571f16577ceSEmil Constantinescu 5728a381b04SJed Brown /* Ydot = shift*(Y-Z) */ 5738a381b04SJed Brown ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 5744f385281SJed Brown for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 5754f385281SJed Brown ierr = VecMAXPY(Z,i,w,YdotI);CHKERRQ(ierr); 576f16577ceSEmil Constantinescu 5778a381b04SJed Brown /* Initial guess taken from last stage */ 5788a381b04SJed Brown ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr); 5798a381b04SJed Brown ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 5808a381b04SJed Brown ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 5818a381b04SJed Brown ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 582f1b97656SJed Brown ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr); 5838a381b04SJed Brown ts->nonlinear_its += its; ts->linear_its += lits; 584476b6736SJed Brown if (snesreason < 0) { 585476b6736SJed Brown if (ts->max_snes_failures > 0 && ++ts->num_snes_failures >= ts->max_snes_failures) { 586f1b97656SJed Brown ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; 587f1b97656SJed 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); 588f1b97656SJed Brown PetscFunctionReturn(0); 589476b6736SJed Brown } else { 590476b6736SJed Brown ts->time_step *= ts->scale_solve_failed; 591476b6736SJed Brown goto reject_step; /* Do not try to solve any more stages */ 592476b6736SJed Brown } 593f1b97656SJed Brown } 5948a381b04SJed Brown } 5958a381b04SJed Brown ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 5964cc180ffSJed Brown ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],ark->imex);CHKERRQ(ierr); 5974cc180ffSJed Brown if (ark->imex) { 5988a381b04SJed Brown ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 5994cc180ffSJed Brown } else { 6004cc180ffSJed Brown ierr = VecZeroEntries(YdotRHS[i]);CHKERRQ(ierr); 6014cc180ffSJed Brown } 6028a381b04SJed Brown } 603108c343cSJed Brown ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 604108c343cSJed Brown ark->status = TS_STEP_PENDING; 6058a381b04SJed Brown 606108c343cSJed Brown /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 607108c343cSJed Brown ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 608108c343cSJed Brown ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 609108c343cSJed Brown ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 610108c343cSJed Brown ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 611108c343cSJed Brown if (accept) { 612108c343cSJed Brown /* ignore next_scheme for now */ 6138a381b04SJed Brown ts->ptime += ts->time_step; 614cdbf8f93SLisandro Dalcin ts->time_step = next_time_step; 6158a381b04SJed Brown ts->steps++; 616108c343cSJed Brown ark->status = TS_STEP_COMPLETE; 617108c343cSJed Brown break; 618108c343cSJed Brown } else { /* Roll back the current step */ 619108c343cSJed Brown for (j=0; j<s; j++) w[j] = h*bt[j]; 620108c343cSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotI);CHKERRQ(ierr); 621108c343cSJed Brown for (j=0; j<s; j++) w[j] = -h*b[j]; 622108c343cSJed Brown ierr = VecMAXPY(ts->vec_sol,s,w,ark->YdotRHS);CHKERRQ(ierr); 623108c343cSJed Brown ts->time_step = next_time_step; 624108c343cSJed Brown ark->status = TS_STEP_INCOMPLETE; 625108c343cSJed Brown } 626476b6736SJed Brown reject_step: continue; 627108c343cSJed Brown } 628476b6736SJed Brown ts->reason = TS_DIVERGED_STEP_REJECTED; 6298a381b04SJed Brown PetscFunctionReturn(0); 6308a381b04SJed Brown } 6318a381b04SJed Brown 632cd652676SJed Brown #undef __FUNCT__ 633cd652676SJed Brown #define __FUNCT__ "TSInterpolate_ARKIMEX" 634cd652676SJed Brown static PetscErrorCode TSInterpolate_ARKIMEX(TS ts,PetscReal itime,Vec X) 635cd652676SJed Brown { 636cd652676SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 6374f385281SJed Brown PetscInt s = ark->tableau->s,pinterp = ark->tableau->pinterp,i,j; 638108c343cSJed Brown PetscReal h; 639108c343cSJed Brown PetscReal tt,t; 640cd652676SJed Brown PetscScalar *bt,*b; 641cd652676SJed Brown const PetscReal *Bt = ark->tableau->binterpt,*B = ark->tableau->binterp; 642cd652676SJed Brown PetscErrorCode ierr; 643cd652676SJed Brown 644cd652676SJed Brown PetscFunctionBegin; 645cd652676SJed Brown if (!Bt || !B) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSARKIMEX %s does not have an interpolation formula",ark->tableau->name); 646108c343cSJed Brown switch (ark->status) { 647108c343cSJed Brown case TS_STEP_INCOMPLETE: 648108c343cSJed Brown case TS_STEP_PENDING: 649108c343cSJed Brown h = ts->time_step; 650108c343cSJed Brown t = (itime - ts->ptime)/h; 651108c343cSJed Brown break; 652108c343cSJed Brown case TS_STEP_COMPLETE: 653108c343cSJed Brown h = ts->time_step_prev; 654108c343cSJed Brown t = (itime - ts->ptime)/h + 1; /* In the interval [0,1] */ 655108c343cSJed Brown break; 656b9ce6d65SJed Brown default: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_PLIB,"Invalid TSStepStatus"); 657108c343cSJed Brown } 658cd652676SJed Brown ierr = PetscMalloc2(s,PetscScalar,&bt,s,PetscScalar,&b);CHKERRQ(ierr); 659cd652676SJed Brown for (i=0; i<s; i++) bt[i] = b[i] = 0; 6604f385281SJed Brown for (j=0,tt=t; j<pinterp; j++,tt*=t) { 661cd652676SJed Brown for (i=0; i<s; i++) { 662108c343cSJed Brown bt[i] += h * Bt[i*pinterp+j] * tt * -1.0; 663108c343cSJed Brown b[i] += h * B[i*pinterp+j] * tt; 664cd652676SJed Brown } 665cd652676SJed Brown } 666cd652676SJed 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"); 667cd652676SJed Brown ierr = VecCopy(ark->Y[0],X);CHKERRQ(ierr); 668cd652676SJed Brown ierr = VecMAXPY(X,s,bt,ark->YdotI);CHKERRQ(ierr); 669cd652676SJed Brown ierr = VecMAXPY(X,s,b,ark->YdotRHS);CHKERRQ(ierr); 670cd652676SJed Brown ierr = PetscFree2(bt,b);CHKERRQ(ierr); 671cd652676SJed Brown PetscFunctionReturn(0); 672cd652676SJed Brown } 673cd652676SJed Brown 6748a381b04SJed Brown /*------------------------------------------------------------*/ 6758a381b04SJed Brown #undef __FUNCT__ 6768a381b04SJed Brown #define __FUNCT__ "TSReset_ARKIMEX" 6778a381b04SJed Brown static PetscErrorCode TSReset_ARKIMEX(TS ts) 6788a381b04SJed Brown { 6798a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 6808a381b04SJed Brown PetscInt s; 6818a381b04SJed Brown PetscErrorCode ierr; 6828a381b04SJed Brown 6838a381b04SJed Brown PetscFunctionBegin; 6848a381b04SJed Brown if (!ark->tableau) PetscFunctionReturn(0); 6858a381b04SJed Brown s = ark->tableau->s; 6868a381b04SJed Brown ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 6878a381b04SJed Brown ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 6888a381b04SJed Brown ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 6898a381b04SJed Brown ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 6908a381b04SJed Brown ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 6918a381b04SJed Brown ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 6928a381b04SJed Brown ierr = PetscFree(ark->work);CHKERRQ(ierr); 6938a381b04SJed Brown PetscFunctionReturn(0); 6948a381b04SJed Brown } 6958a381b04SJed Brown 6968a381b04SJed Brown #undef __FUNCT__ 6978a381b04SJed Brown #define __FUNCT__ "TSDestroy_ARKIMEX" 6988a381b04SJed Brown static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 6998a381b04SJed Brown { 7008a381b04SJed Brown PetscErrorCode ierr; 7018a381b04SJed Brown 7028a381b04SJed Brown PetscFunctionBegin; 7038a381b04SJed Brown ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 7048a381b04SJed Brown ierr = PetscFree(ts->data);CHKERRQ(ierr); 7058a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","",PETSC_NULL);CHKERRQ(ierr); 7068a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","",PETSC_NULL);CHKERRQ(ierr); 707995b3938SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","",PETSC_NULL);CHKERRQ(ierr); 7088a381b04SJed Brown PetscFunctionReturn(0); 7098a381b04SJed Brown } 7108a381b04SJed Brown 7118a381b04SJed Brown /* 7128a381b04SJed Brown This defines the nonlinear equation that is to be solved with SNES 7138a381b04SJed Brown G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 7148a381b04SJed Brown */ 7158a381b04SJed Brown #undef __FUNCT__ 7168a381b04SJed Brown #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 7178a381b04SJed Brown static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 7188a381b04SJed Brown { 7198a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7208a381b04SJed Brown PetscErrorCode ierr; 7218a381b04SJed Brown 7228a381b04SJed Brown PetscFunctionBegin; 7238a381b04SJed Brown ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 7244cc180ffSJed Brown ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,F,ark->imex);CHKERRQ(ierr); 7258a381b04SJed Brown PetscFunctionReturn(0); 7268a381b04SJed Brown } 7278a381b04SJed Brown 7288a381b04SJed Brown #undef __FUNCT__ 7298a381b04SJed Brown #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 7308a381b04SJed Brown static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 7318a381b04SJed Brown { 7328a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7338a381b04SJed Brown PetscErrorCode ierr; 7348a381b04SJed Brown 7358a381b04SJed Brown PetscFunctionBegin; 7368a381b04SJed Brown /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 7378a381b04SJed Brown ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 7388a381b04SJed Brown PetscFunctionReturn(0); 7398a381b04SJed Brown } 7408a381b04SJed Brown 7418a381b04SJed Brown #undef __FUNCT__ 7428a381b04SJed Brown #define __FUNCT__ "TSSetUp_ARKIMEX" 7438a381b04SJed Brown static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 7448a381b04SJed Brown { 7458a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7468a381b04SJed Brown ARKTableau tab = ark->tableau; 7478a381b04SJed Brown PetscInt s = tab->s; 7488a381b04SJed Brown PetscErrorCode ierr; 7498a381b04SJed Brown 7508a381b04SJed Brown PetscFunctionBegin; 7518a381b04SJed Brown if (!ark->tableau) { 752e24355feSJed Brown ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 7538a381b04SJed Brown } 7548a381b04SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 7558a381b04SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 7568a381b04SJed Brown ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 7578a381b04SJed Brown ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 7588a381b04SJed Brown ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 7598a381b04SJed Brown ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 7608a381b04SJed Brown ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 7618a381b04SJed Brown PetscFunctionReturn(0); 7628a381b04SJed Brown } 7638a381b04SJed Brown /*------------------------------------------------------------*/ 7648a381b04SJed Brown 7658a381b04SJed Brown #undef __FUNCT__ 7668a381b04SJed Brown #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 7678a381b04SJed Brown static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 7688a381b04SJed Brown { 7694cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 7708a381b04SJed Brown PetscErrorCode ierr; 7718a381b04SJed Brown char arktype[256]; 7728a381b04SJed Brown 7738a381b04SJed Brown PetscFunctionBegin; 7748a381b04SJed Brown ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 7758a381b04SJed Brown { 7768a381b04SJed Brown ARKTableauLink link; 7778a381b04SJed Brown PetscInt count,choice; 7788a381b04SJed Brown PetscBool flg; 7798a381b04SJed Brown const char **namelist; 7808a381b04SJed Brown ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof arktype);CHKERRQ(ierr); 7818a381b04SJed Brown for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 7828a381b04SJed Brown ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 7838a381b04SJed Brown for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 7848a381b04SJed Brown ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 7858a381b04SJed Brown ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 7868a381b04SJed Brown ierr = PetscFree(namelist);CHKERRQ(ierr); 7874cc180ffSJed Brown flg = (PetscBool)!ark->imex; 7884cc180ffSJed Brown ierr = PetscOptionsBool("-ts_arkimex_fully_implicit","Solve the problem fully implicitly","TSARKIMEXSetFullyImplicit",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 7894cc180ffSJed Brown ark->imex = (PetscBool)!flg; 790d52bd9f3SBarry Smith ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 7918a381b04SJed Brown } 7928a381b04SJed Brown ierr = PetscOptionsTail();CHKERRQ(ierr); 7938a381b04SJed Brown PetscFunctionReturn(0); 7948a381b04SJed Brown } 7958a381b04SJed Brown 7968a381b04SJed Brown #undef __FUNCT__ 7978a381b04SJed Brown #define __FUNCT__ "PetscFormatRealArray" 7988a381b04SJed Brown static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 7998a381b04SJed Brown { 800257d2499SJed Brown PetscErrorCode ierr; 801f1d86077SJed Brown PetscInt i; 802f1d86077SJed Brown size_t left,count; 8038a381b04SJed Brown char *p; 8048a381b04SJed Brown 8058a381b04SJed Brown PetscFunctionBegin; 806f1d86077SJed Brown for (i=0,p=buf,left=len; i<n; i++) { 807f1d86077SJed Brown ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 8088a381b04SJed Brown if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 8098a381b04SJed Brown left -= count; 8108a381b04SJed Brown p += count; 8118a381b04SJed Brown *p++ = ' '; 8128a381b04SJed Brown } 8138a381b04SJed Brown p[i ? 0 : -1] = 0; 8148a381b04SJed Brown PetscFunctionReturn(0); 8158a381b04SJed Brown } 8168a381b04SJed Brown 8178a381b04SJed Brown #undef __FUNCT__ 8188a381b04SJed Brown #define __FUNCT__ "TSView_ARKIMEX" 8198a381b04SJed Brown static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 8208a381b04SJed Brown { 8218a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 8228a381b04SJed Brown ARKTableau tab = ark->tableau; 8238a381b04SJed Brown PetscBool iascii; 8248a381b04SJed Brown PetscErrorCode ierr; 8258a381b04SJed Brown 8268a381b04SJed Brown PetscFunctionBegin; 8278a381b04SJed Brown ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 8288a381b04SJed Brown if (iascii) { 8298a381b04SJed Brown const TSARKIMEXType arktype; 8308a381b04SJed Brown char buf[512]; 8318a381b04SJed Brown ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 8328a381b04SJed Brown ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 8338a381b04SJed Brown ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 83431f6fcc0SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 83531f6fcc0SJed Brown ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 83631f6fcc0SJed Brown ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 8378a381b04SJed Brown } 838d52bd9f3SBarry Smith ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 8398a381b04SJed Brown PetscFunctionReturn(0); 8408a381b04SJed Brown } 8418a381b04SJed Brown 8428a381b04SJed Brown #undef __FUNCT__ 8438a381b04SJed Brown #define __FUNCT__ "TSARKIMEXSetType" 8448a381b04SJed Brown /*@C 8458a381b04SJed Brown TSARKIMEXSetType - Set the type of ARK IMEX scheme 8468a381b04SJed Brown 8478a381b04SJed Brown Logically collective 8488a381b04SJed Brown 8498a381b04SJed Brown Input Parameter: 8508a381b04SJed Brown + ts - timestepping context 8518a381b04SJed Brown - arktype - type of ARK-IMEX scheme 8528a381b04SJed Brown 8538a381b04SJed Brown Level: intermediate 8548a381b04SJed Brown 855020d8f30SJed Brown .seealso: TSARKIMEXGetType(), TSARKIMEX, TSARKIMEX2D, TSARKIMEX2E, TSARKIMEXPRSSP2, TSARKIMEX3, TSARKIMEXBPR3, TSARKIMEXARS443, TSARKIMEX4, TSARKIMEX5 8568a381b04SJed Brown @*/ 8578a381b04SJed Brown PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype) 8588a381b04SJed Brown { 8598a381b04SJed Brown PetscErrorCode ierr; 8608a381b04SJed Brown 8618a381b04SJed Brown PetscFunctionBegin; 8628a381b04SJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 8638a381b04SJed Brown ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 8648a381b04SJed Brown PetscFunctionReturn(0); 8658a381b04SJed Brown } 8668a381b04SJed Brown 8678a381b04SJed Brown #undef __FUNCT__ 8688a381b04SJed Brown #define __FUNCT__ "TSARKIMEXGetType" 8698a381b04SJed Brown /*@C 8708a381b04SJed Brown TSARKIMEXGetType - Get the type of ARK IMEX scheme 8718a381b04SJed Brown 8728a381b04SJed Brown Logically collective 8738a381b04SJed Brown 8748a381b04SJed Brown Input Parameter: 8758a381b04SJed Brown . ts - timestepping context 8768a381b04SJed Brown 8778a381b04SJed Brown Output Parameter: 8788a381b04SJed Brown . arktype - type of ARK-IMEX scheme 8798a381b04SJed Brown 8808a381b04SJed Brown Level: intermediate 8818a381b04SJed Brown 8828a381b04SJed Brown .seealso: TSARKIMEXGetType() 8838a381b04SJed Brown @*/ 8848a381b04SJed Brown PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype) 8858a381b04SJed Brown { 8868a381b04SJed Brown PetscErrorCode ierr; 8878a381b04SJed Brown 8888a381b04SJed Brown PetscFunctionBegin; 8898a381b04SJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 8908a381b04SJed Brown ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 8918a381b04SJed Brown PetscFunctionReturn(0); 8928a381b04SJed Brown } 8938a381b04SJed Brown 8944cc180ffSJed Brown #undef __FUNCT__ 8954cc180ffSJed Brown #define __FUNCT__ "TSARKIMEXSetFullyImplicit" 8964cc180ffSJed Brown /*@C 8974cc180ffSJed Brown TSARKIMEXSetFullyImplicit - Solve both parts of the equation implicitly 8984cc180ffSJed Brown 8994cc180ffSJed Brown Logically collective 9004cc180ffSJed Brown 9014cc180ffSJed Brown Input Parameter: 9024cc180ffSJed Brown + ts - timestepping context 9034cc180ffSJed Brown - flg - PETSC_TRUE for fully implicit 9044cc180ffSJed Brown 9054cc180ffSJed Brown Level: intermediate 9064cc180ffSJed Brown 9074cc180ffSJed Brown .seealso: TSARKIMEXGetType() 9084cc180ffSJed Brown @*/ 9094cc180ffSJed Brown PetscErrorCode TSARKIMEXSetFullyImplicit(TS ts,PetscBool flg) 9104cc180ffSJed Brown { 9114cc180ffSJed Brown PetscErrorCode ierr; 9124cc180ffSJed Brown 9134cc180ffSJed Brown PetscFunctionBegin; 9144cc180ffSJed Brown PetscValidHeaderSpecific(ts,TS_CLASSID,1); 9154cc180ffSJed Brown ierr = PetscTryMethod(ts,"TSARKIMEXSetFullyImplicit_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 9164cc180ffSJed Brown PetscFunctionReturn(0); 9174cc180ffSJed Brown } 9184cc180ffSJed Brown 9198a381b04SJed Brown EXTERN_C_BEGIN 9208a381b04SJed Brown #undef __FUNCT__ 9218a381b04SJed Brown #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 9228a381b04SJed Brown PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype) 9238a381b04SJed Brown { 9248a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 9258a381b04SJed Brown PetscErrorCode ierr; 9268a381b04SJed Brown 9278a381b04SJed Brown PetscFunctionBegin; 9288a381b04SJed Brown if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 9298a381b04SJed Brown *arktype = ark->tableau->name; 9308a381b04SJed Brown PetscFunctionReturn(0); 9318a381b04SJed Brown } 9328a381b04SJed Brown #undef __FUNCT__ 9338a381b04SJed Brown #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 9348a381b04SJed Brown PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype) 9358a381b04SJed Brown { 9368a381b04SJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 9378a381b04SJed Brown PetscErrorCode ierr; 9388a381b04SJed Brown PetscBool match; 9398a381b04SJed Brown ARKTableauLink link; 9408a381b04SJed Brown 9418a381b04SJed Brown PetscFunctionBegin; 9428a381b04SJed Brown if (ark->tableau) { 9438a381b04SJed Brown ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 9448a381b04SJed Brown if (match) PetscFunctionReturn(0); 9458a381b04SJed Brown } 9468a381b04SJed Brown for (link = ARKTableauList; link; link=link->next) { 9478a381b04SJed Brown ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 9488a381b04SJed Brown if (match) { 9498a381b04SJed Brown ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 9508a381b04SJed Brown ark->tableau = &link->tab; 9518a381b04SJed Brown PetscFunctionReturn(0); 9528a381b04SJed Brown } 9538a381b04SJed Brown } 9548a381b04SJed Brown SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 9558a381b04SJed Brown PetscFunctionReturn(0); 9568a381b04SJed Brown } 9574cc180ffSJed Brown #undef __FUNCT__ 9584cc180ffSJed Brown #define __FUNCT__ "TSARKIMEXSetFullyImplicit_ARKIMEX" 9594cc180ffSJed Brown PetscErrorCode TSARKIMEXSetFullyImplicit_ARKIMEX(TS ts,PetscBool flg) 9604cc180ffSJed Brown { 9614cc180ffSJed Brown TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 9624cc180ffSJed Brown 9634cc180ffSJed Brown PetscFunctionBegin; 9644cc180ffSJed Brown ark->imex = (PetscBool)!flg; 9654cc180ffSJed Brown PetscFunctionReturn(0); 9664cc180ffSJed Brown } 9678a381b04SJed Brown EXTERN_C_END 9688a381b04SJed Brown 9698a381b04SJed Brown /* ------------------------------------------------------------ */ 9708a381b04SJed Brown /*MC 9718a381b04SJed Brown TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes 9728a381b04SJed Brown 973fca742c7SJed Brown These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 974fca742c7SJed Brown nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 975fca742c7SJed Brown of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 976fca742c7SJed Brown 977fca742c7SJed Brown Notes: 978c8058688SBarry Smith The default is TSARKIMEX2E, it can be changed with TSARKIMEXSetType() or -ts_arkimex_type 979c8058688SBarry Smith 980fca742c7SJed Brown This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 981fca742c7SJed Brown 9828a381b04SJed Brown Level: beginner 9838a381b04SJed Brown 984c8058688SBarry Smith .seealso: TSCreate(), TS, TSSetType(), TSARKIMEXSetType(), TSARKIMEXGetType(), TSARKIMEXSetFullyImplicit(), TSARKIMEX2D, TTSARKIMEX2E, TSARKIMEX3, 98526885f59SBarry Smith TSARKIMEX4, TSARKIMEX5, TSARKIMEXType, TSARKIMEXRegister() 9868a381b04SJed Brown 9878a381b04SJed Brown M*/ 9888a381b04SJed Brown EXTERN_C_BEGIN 9898a381b04SJed Brown #undef __FUNCT__ 9908a381b04SJed Brown #define __FUNCT__ "TSCreate_ARKIMEX" 9918a381b04SJed Brown PetscErrorCode TSCreate_ARKIMEX(TS ts) 9928a381b04SJed Brown { 9938a381b04SJed Brown TS_ARKIMEX *th; 9948a381b04SJed Brown PetscErrorCode ierr; 9958a381b04SJed Brown 9968a381b04SJed Brown PetscFunctionBegin; 9978a381b04SJed Brown #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 9988a381b04SJed Brown ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 9998a381b04SJed Brown #endif 10008a381b04SJed Brown 10018a381b04SJed Brown ts->ops->reset = TSReset_ARKIMEX; 10028a381b04SJed Brown ts->ops->destroy = TSDestroy_ARKIMEX; 10038a381b04SJed Brown ts->ops->view = TSView_ARKIMEX; 10048a381b04SJed Brown ts->ops->setup = TSSetUp_ARKIMEX; 10058a381b04SJed Brown ts->ops->step = TSStep_ARKIMEX; 1006cd652676SJed Brown ts->ops->interpolate = TSInterpolate_ARKIMEX; 1007108c343cSJed Brown ts->ops->evaluatestep = TSEvaluateStep_ARKIMEX; 10088a381b04SJed Brown ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 10098a381b04SJed Brown ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 10108a381b04SJed Brown ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 10118a381b04SJed Brown 10128a381b04SJed Brown ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 10138a381b04SJed Brown ts->data = (void*)th; 10144cc180ffSJed Brown th->imex = PETSC_TRUE; 10158a381b04SJed Brown 10168a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 10178a381b04SJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 10184cc180ffSJed Brown ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetFullyImplicit_C","TSARKIMEXSetFullyImplicit_ARKIMEX",TSARKIMEXSetFullyImplicit_ARKIMEX);CHKERRQ(ierr); 10198a381b04SJed Brown PetscFunctionReturn(0); 10208a381b04SJed Brown } 10218a381b04SJed Brown EXTERN_C_END 1022