Skip to content

Commit

Permalink
New rk methods added (#10160)
Browse files Browse the repository at this point in the history 
* Error in Dense output formular of esdirk methods corrected
* Tsit5 and classical RK scheme added
  • Loading branch information
@bernhardbachmann
bernhardbachmann committed Feb 6, 2023
1 parent eac10d9 commit c05a31d
Show file tree
Hide file tree
Showing 3 changed files with 98 additions and 0 deletions.
92 changes: 92 additions & 0 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c
View file
Expand Up @@ -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) {

Expand Down Expand Up @@ -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) {

Expand Down Expand Up @@ -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;
Expand Down
4 changes: 4 additions & 0 deletions OMCompiler/SimulationRuntime/c/util/simulation_options.c
View file
Expand Up @@ -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"
};

Expand Down Expand Up @@ -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)"
};

Expand Down
2 changes: 2 additions & 0 deletions OMCompiler/SimulationRuntime/c/util/simulation_options.h
View file
Expand Up @@ -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 */


Expand Down

0 comments on commit c05a31d

Please sign in to comment.