diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c index 8615de2d444..a5b6fcca6c1 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c @@ -185,6 +185,49 @@ void getButcherTableau_ESDIRK3(BUTCHER_TABLEAU* tableau) { tableau->isKRightAvailable = TRUE; } +// TODO: Describe me +void denseOutput_TSIT5(BUTCHER_TABLEAU* tableau, double* yOld, double* x, double* k, double dt, double stepSize, double* y, int nIdx, int* idx, int nStates) { + +tableau->b_dt[0] = (( -1.0530884977290216 * dt + 2.913255461821912743750681) * dt -2.763706197274825911336736) * dt + 1; +tableau->b_dt[1] = (( 0.1017 * dt -0.22339999999999999818) * dt + 0.13169999999999999727) * dt; +tableau->b_dt[2] = (( 2.490627285651252793 * dt -5.941033872131504734702492) * dt + 3.930296236894751528506874) * dt; +tableau->b_dt[3] = (( -16.54810288924490272 * dt + 30.33818863028232159817299) * dt -12.41107716693367698373438) * dt; +tableau->b_dt[4] = (( 47.37952196281928122 * dt -88.17890489476640110142767) * dt + 37.50931341651103919496904) * dt; +tableau->b_dt[5] = (( -34.87065786149660974 * dt + 65.09189467479367152629022) * dt -27.89652628919728780594826) * dt; +tableau->b_dt[6] = (( 2.5 * dt -4) * dt + 1.5) * dt; + + denseOutput(tableau, yOld, x, k, dt, stepSize, y, nIdx, idx, nStates); +} + +// TODO: Describe me +void getButcherTableau_TSIT5(BUTCHER_TABLEAU* tableau) { + + tableau->nStages = 7; + tableau->order_b = 5; + tableau->order_bt = 4; + tableau->fac = 1.0; + + const double c[] = { 0, 0.161, 0.327, 0.9, 0.980025540904509685729810286287, 1, 1}; + const double A[] = { + 0, 0, 0, 0, 0, 0, 0, + 0.161, 0, 0, 0, 0, 0, 0, + -0.00848065549235698854442687425023, 0.33548065549235698854442687425, 0, 0, 0, 0, 0, + 2.89715305710549343213043259419, -6.35944848997507484314815991238, 4.36229543286958141101772731819, 0, 0, 0, 0, + 5.32586482843925660442887792084, -11.7488835640628278777471703398, 7.49553934288983620830460478456, -0.0924950663617552492565020793321, 0, 0, 0, + 5.86145544294642002865925148698, -12.9209693178471092917061186818, 8.15936789857615864318040079454, -0.0715849732814009972245305425258, -0.0282690503940683829090030572127, 0, 0, + 0.0964607668180652295181673131651, 0.01, 0.479889650414499574775249532291, 1.37900857410374189319227482186, -3.29006951543608067990104758571, 2.3247105240997739824153559184, 0}; + const double b[] = {0.0964607668180652295181673131651, 0.01, 0.479889650414499574775249532291, 1.37900857410374189319227482186, -3.29006951543608067990104758571, 2.3247105240997739824153559184, 0}; + const double bt[] = {0.0982407778702910009615458637727, 0.0108164344596567469032236360634, 0.472008772404237578764934804618, 1.52371958127700480072943996983, -3.87242668088863590492098519636, 2.78279263002896092907699243723, -0.0151515151515151515151515151515}; + + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + + tableau->withDenseOutput = TRUE; + tableau->dense_output = denseOutput_TSIT5; + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = TRUE; +} + // TODO: Describe me void denseOutput_ESDIRK4(BUTCHER_TABLEAU* tableau, double* yOld, double* x, double* k, double dt, double stepSize, double* y, int nIdx, int* idx, int nStates) { @@ -343,6 +386,49 @@ void getButcherTableau_EXPLEULER(BUTCHER_TABLEAU* tableau) { } } +// TODO: Describe me +void getButcherTableau_RUNGEKUTTA(BUTCHER_TABLEAU* tableau) { + + if (tableau->richardson) { + tableau->nStages = 4; + tableau->order_b = 4; + + /* Butcher Tableau */ + const double c[] = { 0, 0.5, 0.5, 1}; + const double A[] = { + 0, 0, 0, 0, + 0.5, 0, 0, 0, + 0, 0.5, 0, 0, + 0, 0, 1, 0}; + const double b[] = {0.166666666666666666666666666667, 0.333333333333333333333333333333, 0.333333333333333333333333333333, 0.166666666666666666666666666667}; + const double* bt = NULL; + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + tableau->isKLeftAvailable = FALSE; + tableau->isKRightAvailable = FALSE; + } else { + tableau->nStages = 5; + tableau->order_b = 4; + tableau->order_bt = 3; + tableau->fac = 1.0; + + /* Butcher Tableau */ + const double c[] = { 0, 0.5, 0.5, 1, 1}; + const double A[] = { + 0, 0, 0, 0, 0, + 0.5, 0, 0, 0, 0, + 0, 0.5, 0, 0, 0, + 0, 0, 1, 0, 0, + 0.166666666666666666666666666667, 0.333333333333333333333333333333, 0.333333333333333333333333333333, 0.166666666666666666666666666667, 0}; + const double b[] = {0.166666666666666666666666666667, 0.333333333333333333333333333333, 0.333333333333333333333333333333, 0.166666666666666666666666666667, 0}; + const double bt[] = {0.166666666666666666666666666667, 0.333333333333333333333333333333, 0.333333333333333333333333333333, 0.0666666666666666666666666666667, 0.1}; + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = TRUE; + } +} + // TODO: Describe me void getButcherTableau_RADAU_IA_2(BUTCHER_TABLEAU* tableau) { @@ -1334,6 +1420,12 @@ BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ER case RK_TRAPEZOID: getButcherTableau_TRAPEZOID(tableau); break; + case RK_RUNGEKUTTA: + getButcherTableau_RUNGEKUTTA(tableau); + break; + case RK_TSIT5: + getButcherTableau_TSIT5(tableau); + break; case RK_DOPRI45: getButcherTableau_DOPRI45(tableau); break; diff --git a/OMCompiler/SimulationRuntime/c/util/simulation_options.c b/OMCompiler/SimulationRuntime/c/util/simulation_options.c index 0086e4dfa4c..0dfdba935bf 100644 --- a/OMCompiler/SimulationRuntime/c/util/simulation_options.c +++ b/OMCompiler/SimulationRuntime/c/util/simulation_options.c @@ -923,6 +923,8 @@ const char *GB_METHOD_NAME[RK_MAX] = { /* RK_DOPRI45 */ "dopri45", /* RK_DOPRISSC1 */ "dopriSsc1", /* RK_DOPRISSC2 */ "dopriSsc2", + /* RK_TSIT5 */ "tsit5", + /* RK_RUNGEKUTTA */ "rungekutta", /* RK_RKSSC */ "rungekuttaSsc" }; @@ -969,6 +971,8 @@ const char *GB_METHOD_DESC[RK_MAX] = { /* RK_DOPRI45 */ "Explicit Runge-Kutta method Dormand-Prince (order 5)", /* RK_DOPRISSC1 */ "Explicit Runge-Kutta method Dormand-Prince with large stability region (order 1)", /* RK_DOPRISSC2 */ "Explicit Runge-Kutta method Dormand-Prince with large stability region (order 2)", + /* RK_TSIT5 */ "Explicit Runge-Kutta method from Tsitouras (order 5)", + /* RK_RUNGEKUTTA */ "Explicit classical Runge-Kutta method (order 4)", /* RK_RKSSC */ "Explicit Runge-Kutta method with large stabiliy region (order 1)" }; diff --git a/OMCompiler/SimulationRuntime/c/util/simulation_options.h b/OMCompiler/SimulationRuntime/c/util/simulation_options.h index 8285a970098..8629fdf6006 100644 --- a/OMCompiler/SimulationRuntime/c/util/simulation_options.h +++ b/OMCompiler/SimulationRuntime/c/util/simulation_options.h @@ -245,6 +245,8 @@ enum GB_METHOD { RK_DOPRI45, /* dopri45*/ RK_DOPRISSC1, /* dopriSsc1*/ RK_DOPRISSC2, /* dopriSsc2*/ + RK_TSIT5, /* tsit5*/ + RK_RUNGEKUTTA, /* rungekutta*/ RK_RKSSC, /* rungekuttaSsc */