Skip to content

Commit

Permalink
Adding rungekuttaSsc, heun and trapezoid method to gbode framework (#…
Browse files Browse the repository at this point in the history 
…9232)

* Added Butcher tableau of the rungekuttaSsc method to the gbode framework
see:
 Solving Stiff Systems of ODEs by Explicit Methods with Conformed Stability Domains
 Anton E. Novikov1     Mikhail V. Rybkov2     Yury V. Shornikov3     Lyudmila V. Knaub4
 EUROSIM 2016 & SIMS 2016
* RK method heun added
* RK method trapezoid added
  • Loading branch information
@bernhardbachmann
bernhardbachmann committed Jul 13, 2022
1 parent e94216e commit e1e122d
Show file tree
Hide file tree
Showing 3 changed files with 97 additions and 8 deletions.
85 changes: 82 additions & 3 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c
View file
Expand Up @@ -288,6 +288,26 @@ void getButcherTableau_MS(BUTCHER_TABLEAU* tableau)
tableau->isKRightAvailable = TRUE;
}

// https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods
void getButcherTableau_HEUN(BUTCHER_TABLEAU* tableau)
{
tableau->nStages = 2;
tableau->order_b = 2;
tableau->order_bt = 1;
tableau->fac = 1.0;

/* Butcher Tableau */
const double c[] = {0.0, 1.0};
const double A[] = {0.0, 0.0,
1.0, 0.0};
const double b[] = {0.5, 0.5};
const double bt[] = {1.0, 0.0};

setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *)bt);
tableau->isKLeftAvailable = TRUE;
tableau->isKRightAvailable = FALSE;
}

// TODO: Describe me
void getButcherTableau_EXPLEULER(BUTCHER_TABLEAU* tableau) {

Expand Down Expand Up @@ -739,6 +759,27 @@ void getButcherTableau_IMPLEULER(BUTCHER_TABLEAU* tableau) {
tableau->isKRightAvailable = TRUE;
}

// https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods
void getButcherTableau_TRAPEZOID(BUTCHER_TABLEAU* tableau) {

tableau->nStages = 2;
tableau->order_b = 2;
tableau->order_bt = 1;
tableau->fac = 1.e0;

// /* Butcher Tableau */
const double c[] = {0.0, 1.0};
const double A[] = {0.0, 0.0,
0.5, 0.5};
const double b[] = {0.5, 0.5}; // trapezoidal rule
const double bt[] = {1.0, 0.0}; // implicit Euler step

setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt);

tableau->isKLeftAvailable = TRUE;
tableau->isKRightAvailable = TRUE;
}

// TODO: Describe me
void getButcherTableau_MERSON(BUTCHER_TABLEAU* tableau) {

Expand Down Expand Up @@ -1018,6 +1059,35 @@ void getButcherTableau_RK1214(BUTCHER_TABLEAU* tableau) {
tableau->isKRightAvailable = FALSE;
}

/** @brief Get Runge-Kutta Ssc Butcher tableau.
* Solving Stiff Systems of ODEs by Explicit Methods with Conformed Stability Domains
* From: Anton E. Novikov Mikhail V. Rybkov Yury V. Shornikov Lyudmila V. Knaub
* EUROSIM 2016 & SIMS 2016
*
* @param tableau Pointer to Butcher tableau to fill.
*/
void getButcherTableau_RKSSC(BUTCHER_TABLEAU* tableau) {

tableau->nStages = 5;
tableau->order_b = 1;
tableau->order_bt = 2;
tableau->fac = 1e0;

const double c[] = { 0, 0.041324301621055, 0.1611647763, 0.3608883044, 0.64049984};
const double A[] = {
0, 0, 0, 0, 0,
0.041324301621055, 0, 0, 0, 0,
0.0805823881610573, 0.0805823881610573, 0, 0, 0,
0.1191668151228434, 0.1597820013984078, 0.0819394878966193, 0, 0,
0.1570787892802991, 0.237958302195982, 0.1631711307360486, 0.0822916178203657, 0};
const double b[] = { 0.1945277188657676, 0.3151822878089125, 0.2437005934695969, 0.1641555613805598, 0.0824338384751631};
const double bt[] = { 0.12149281854707711872, 0.18276003767888932578, 0.14735209191059650291, -0.42005655782840578172, 0.96845160969184283432};

setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt);
tableau->isKLeftAvailable = FALSE;
tableau->isKRightAvailable = FALSE;
}

/**
* @brief Analyse Butcher tableau and return size and if the method is explicit.
*
Expand Down Expand Up @@ -1092,15 +1162,18 @@ BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ER
case MS_ADAMS_MOULTON:
getButcherTableau_MS(tableau);
break;
case RK_EXPL_EULER:
getButcherTableau_EXPLEULER(tableau);
break;
case RK_IMPL_EULER:
getButcherTableau_IMPLEULER(tableau);
break;
case RK_TRAPEZOID:
getButcherTableau_TRAPEZOID(tableau);
break;
case RK_DOPRI45:
getButcherTableau_DOPRI45(tableau);
break;
case RK_RKSSC:
getButcherTableau_RKSSC(tableau);
break;
case RK_RK810:
getButcherTableau_RK810(tableau);
break;
Expand All @@ -1113,6 +1186,12 @@ BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ER
case RK_MERSON:
getButcherTableau_MERSON(tableau);
break;
case RK_EXPL_EULER:
getButcherTableau_EXPLEULER(tableau);
break;
case RK_HEUN:
getButcherTableau_HEUN(tableau);
break;
case RK_FEHLBERG12:
getButcherTableau_FEHLBERG12(tableau);
break;
Expand Down
16 changes: 11 additions & 5 deletions OMCompiler/SimulationRuntime/c/util/simulation_options.c
View file
Expand Up @@ -876,6 +876,7 @@ const char *GB_METHOD_NAME[RK_MAX] = {
/* MS_ADAMS_MOULTON */ "adams",
/* RK_EXPL_EULER */ "expl_euler",
/* RK_IMPL_EULER */ "impl_euler",
/* RK_TRAPEZOID */ "trapezoid",
/* RK_SDIRK2 */ "sdirk2",
/* RK_SDIRK3 */ "sdirk3",
/* RK_ESDIRK2 */ "esdirk2",
Expand All @@ -899,20 +900,23 @@ const char *GB_METHOD_NAME[RK_MAX] = {
/* RK_GAUSS5 */ "gauss5",
/* RK_GAUSS6 */ "gauss6",
/* RK_MERSON */ "merson",
/* RK_HEUN */ "heun",
/* RK_FEHLBERG12 */ "fehlberg12",
/* RK_FEHLBERG45 */ "fehlberg45",
/* RK_FEHLBERG78 */ "fehlberg78",
/* RK_RK810 */ "rk810",
/* RK_RK1012 */ "rk1012",
/* RK_RK1214 */ "rk1214",
/* RK_DOPRI45 */ "dopri45"
/* RK_DOPRI45 */ "dopri45",
/* RK_RKSSC */ "rungekuttaSsc"
};

const char *GB_METHOD_DESC[RK_MAX] = {
/* GB_UNKNOWN = 0 */ "unknown",
/* MS_ADAMS_MOULTON */ "adams",
/* RK_EXPL_EULER */ "Explizit Euler method (order 1)",
/* RK_IMPL_EULER */ "Implizit Euler method (order 1)",
/* MS_ADAMS_MOULTON */ "Implicit multistep method of type Adams-Moulton (order 2)",
/* RK_EXPL_EULER */ "Explizit Runge-Kutta Euler method (order 1)",
/* RK_IMPL_EULER */ "Implizit Runge-Kutta Euler method (order 1)",
/* RK_TRAPEZOID */ "Implicit Runge-Kutta trapezoid method (order 2)",
/* RK_SDIRK2 */ "Singly-diagonal implicit Runge-Kutta (order 2)",
/* RK_SDIRK3 */ "Singly-diagonal implicit Runge-Kutta (order 3)",
/* RK_ESDIRK2 */ "Explicit singly-diagonal implicit Runge-Kutta (order 2)",
Expand All @@ -936,13 +940,15 @@ const char *GB_METHOD_DESC[RK_MAX] = {
/* RK_GAUSS5 */ "Implicit Runge-Kutta method of Gauss (order 10)",
/* RK_GAUSS6 */ "Implicit Runge-Kutta method of Gauss (order 12)",
/* RK_MERSON */ "Explicit Runge-Kutta Merson method (order 4)",
/* RK_HEUN */ "Explicit Runge-Kutta Heun method (order 2)",
/* RK_FEHLBERG12 */ "Explicit Runge-Kutta Fehlberg method (order 2)",
/* RK_FEHLBERG45 */ "Explicit Runge-Kutta Fehlberg method (order 5)",
/* RK_FEHLBERG78 */ "Explicit Runge-Kutta Fehlberg method (order 8)",
/* RK_RK810 */ "Explicit 8-10 Runge-Kutta method (order 10)",
/* RK_RK1012 */ "Explicit 10-12 Runge-Kutta method (order 12)",
/* RK_RK1214 */ "Explicit 12-14 Runge-Kutta method (order 14)",
/* RK_DOPRI45 */ "Explicit Runge-Kutta method Dormand-Prince (order 5)"
/* RK_DOPRI45 */ "Explicit Runge-Kutta method Dormand-Prince (order 5)",
/* RK_RKSSC */ "Explicit Runge-Kutta method with large stabiliy region (order 1)"
};

const char *GB_NLS_METHOD_NAME[GB_NLS_MAX] = {
Expand Down
4 changes: 4 additions & 0 deletions OMCompiler/SimulationRuntime/c/util/simulation_options.h
View file
Expand Up @@ -206,6 +206,7 @@ enum GB_METHOD {
MS_ADAMS_MOULTON, /* adams*/
RK_EXPL_EULER, /* expl_euler*/
RK_IMPL_EULER, /* impl_euler*/
RK_TRAPEZOID, /* trapezoid */
RK_SDIRK2, /* sdirk2*/
RK_SDIRK3, /* sdirk3*/
RK_ESDIRK2, /* esdirk2*/
Expand All @@ -229,13 +230,16 @@ enum GB_METHOD {
RK_GAUSS5, /* gauss5*/
RK_GAUSS6, /* gauss6*/
RK_MERSON, /* merson*/
RK_HEUN, /* heun */
RK_FEHLBERG12, /* fehlberg12*/
RK_FEHLBERG45, /* fehlberg45*/
RK_FEHLBERG78, /* fehlberg78*/
RK_RK810, /* rk810*/
RK_RK1012, /* rk1012*/
RK_RK1214, /* rk1214*/
RK_DOPRI45, /* dopri4*/
RK_RKSSC, /* rungekuttaSsc */


RK_MAX
};
Expand Down

0 comments on commit e1e122d

Please sign in to comment.