From 4c642c36ba2501073df7567c4b86ae8ef2e29b63 Mon Sep 17 00:00:00 2001 From: bernhardbachmann Date: Fri, 28 Oct 2022 16:26:07 +0200 Subject: [PATCH] New RK tableau with large stability regions added: (#9615) -gbm=dopriSsc1 (order 1, stages 7) -gbm=dopriSsc2 (order 2, stages 7) -gbm=mersonSsc1 (order 1, stages 5) -gbm=mersonSsc2 (order 2, stages 5) -gbm=fehlbergSsc1 (order 1, stages 13) -gbm=fehlbergSsc2 (order 2, stages 13) --- .../c/simulation/solver/gbode_tableau.c | 190 +++++++++++++++++- .../c/util/simulation_options.c | 12 ++ .../c/util/simulation_options.h | 8 +- 3 files changed, 206 insertions(+), 4 deletions(-) diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c index ac15c80d4b0..e592ebb7c8e 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c @@ -786,7 +786,7 @@ void getButcherTableau_MERSON(BUTCHER_TABLEAU* tableau) { tableau->nStages = 5; tableau->order_b = 4; tableau->order_bt = 3; - tableau->fac = 1.e5; + tableau->fac = 1e5; /* Butcher Tableau */ const double c[] = {0.0, 1./3, 1./3, 1./2, 1.0}; @@ -796,8 +796,54 @@ void getButcherTableau_MERSON(BUTCHER_TABLEAU* tableau) { 1./8, 0.0, 3./8, 0.0, 0.0, 1./2, 0.0, -3./2, 2.0, 0.0 }; - const double b[] = {1./6, 0.0, 0.0, 2./3, 1./6}; // 4th order??? - const double bt[] = {1./10, 0.0, 3./10, 2./5, 1./5}; // 3th order??? + const double b[] = {1./6, 0.0, 0.0, 2./3, 1./6}; // 4th order + const double bt[] = {1./10, 0.0, 3./10, 2./5, 1./5}; // 3th order + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = FALSE; +} + +void getButcherTableau_MERSONSSC1(BUTCHER_TABLEAU* tableau) { + + tableau->nStages = 5; + tableau->order_b = 1; + tableau->order_bt = 4; + tableau->fac = 1e0; + + /* Butcher Tableau */ + const double c[] = {0.0, 1./3, 1./3, 1./2, 1.0}; + const double A[] = { 0.0, 0.0, 0.0, 0.0, 0.0, + 1./3, 0.0, 0.0, 0.0, 0.0, + 1./6, 1./6, 0.0, 0.0, 0.0, + 1./8, 0.0, 3./8, 0.0, 0.0, + 1./2, 0.0, -3./2, 2.0, 0.0 + }; + const double b[] = {0.512782397120662718471749459233, 0.330103091995730873405418477521, 0.146713304970630735231072129528, 0.0103570041584038638446238467251, 4.42017545718090471360869927930e-05}; + const double bt[] = {1./6, 0.0, 0.0, 2./3, 1./6}; // 4th order + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = FALSE; +} + +void getButcherTableau_MERSONSSC2(BUTCHER_TABLEAU* tableau) { + + tableau->nStages = 5; + tableau->order_b = 2; + tableau->order_bt = 4; + tableau->fac = 1e0; + + /* Butcher Tableau */ + const double c[] = {0.0, 1./3, 1./3, 1./2, 1.0}; + const double A[] = { 0.0, 0.0, 0.0, 0.0, 0.0, + 1./3, 0.0, 0.0, 0.0, 0.0, + 1./6, 1./6, 0.0, 0.0, 0.0, + 1./8, 0.0, 3./8, 0.0, 0.0, + 1./2, 0.0, -3./2, 2.0, 0.0 + }; + const double b[] = {-0.35629337268078937564325457003, 0.146074652453948837652304806997, 0.934217301122925451486796787885, 0.272197473925746365767707013552, 0.00380394517816872073644596159715}; + const double bt[] = {1./6, 0.0, 0.0, 2./3, 1./6}; // 4th order setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); tableau->isKLeftAvailable = TRUE; @@ -846,6 +892,64 @@ void getButcherTableau_DOPRI45(BUTCHER_TABLEAU* tableau) { tableau->isKRightAvailable = TRUE; } +// TODO: Describe me +void getButcherTableau_DOPRISSC1(BUTCHER_TABLEAU* tableau) { + + tableau->nStages = 7; + tableau->order_b = 1; + tableau->order_bt = 5; + tableau->fac = 1e0; + + /* Butcher Tableau */ + const double c[] = {0.0, 1./5, 3./10, 4./5, 8./9, 1., 1.}; + const double A[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 1./5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 3./40, 9./40, 0.0, 0.0, 0.0, 0.0, 0.0, + 44./45, -56./15, 32./9, 0.0, 0.0, 0.0, 0.0, + 19372./6561, -25360./2187, 64448./6561, -212./729, 0.0, 0.0, 0.0, + 9017./3168, -355./33, 46732./5247, 49./176, -5103./18656, 0.0, 0.0, + 35./384, 0.0, 500./1113, 125./192, -2187./6784, 11./84, 0.0 + }; + const double b[] = {0.278585202707552297491652379451, 0.499359343897282505016199003177, 0.21994590092478885648226620836, 0.00221513041070919707891834807597, -0.000108554006807565712812909262366, 2.90820039199235629183419683848e-06, 6.78660827172874851360023304864e-08}; + const double bt[] = {35./384, 0.0, 500./1113, 125./192, -2187./6784, 11./84, 0.0}; + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + + tableau->withDenseOutput = TRUE; + tableau->dense_output = denseOutput_DOPRI45; + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = FALSE; +} + +// TODO: Describe me +void getButcherTableau_DOPRISSC2(BUTCHER_TABLEAU* tableau) { + + tableau->nStages = 7; + tableau->order_b = 2; + tableau->order_bt = 5; + tableau->fac = 1e0; + + /* Butcher Tableau */ + const double c[] = {0.0, 1./5, 3./10, 4./5, 8./9, 1., 1.}; + const double A[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 1./5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 3./40, 9./40, 0.0, 0.0, 0.0, 0.0, 0.0, + 44./45, -56./15, 32./9, 0.0, 0.0, 0.0, 0.0, + 19372./6561, -25360./2187, 64448./6561, -212./729, 0.0, 0.0, 0.0, + 9017./3168, -355./33, 46732./5247, 49./176, -5103./18656, 0.0, 0.0, + 35./384, 0.0, 500./1113, 125./192, -2187./6784, 11./84, 0.0 + }; + const double b[] = {-0.486346436598901828254513839047, -0.234874439261298143693150869933, 1.65868062062825029557032033231, 0.0708767352953961545635216713703, -0.00905214141822685142604628709823, 0.000667704407394041424018354604235, 4.79569473863318158506377905590e-05}; + const double bt[] = {35./384, 0.0, 500./1113, 125./192, -2187./6784, 11./84, 0.0}; + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + + tableau->withDenseOutput = TRUE; + tableau->dense_output = denseOutput_DOPRI45; + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = FALSE; +} + // TODO: Describe me void getButcherTableau_FEHLBERG12(BUTCHER_TABLEAU* tableau) { @@ -924,6 +1028,68 @@ void getButcherTableau_FEHLBERG78(BUTCHER_TABLEAU* tableau) { tableau->isKRightAvailable = FALSE; } +void getButcherTableau_FEHLBERGSSC1(BUTCHER_TABLEAU* tableau) { + + tableau->nStages = 13; + tableau->order_b = 1; + tableau->order_bt = 8; + tableau->fac = 1e0; + + /* Butcher Tableau */ + const double c[] = {0, 0.0740740740740740740740740740741, 0.111111111111111111111111111111, 0.166666666666666666666666666667, 0.416666666666666666666666666667, 0.5, 0.833333333333333333333333333333, 0.166666666666666666666666666667, 0.666666666666666666666666666667, 0.333333333333333333333333333333, 1, 0, 1}; + const double A[] = { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.0740740740740740740740740740741, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.0277777777777777777777777777778, 0.0833333333333333333333333333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.0416666666666666666666666666667, 0, 0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.416666666666666666666666666667, 0, -1.5625, 1.5625, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.05, 0, 0, 0.25, 0.2, 0, 0, 0, 0, 0, 0, 0, 0, + -0.231481481481481481481481481481, 0, 0, 1.15740740740740740740740740741, -2.40740740740740740740740740741, 2.31481481481481481481481481481, 0, 0, 0, 0, 0, 0, 0, + 0.103333333333333333333333333333, 0, 0, 0, 0.271111111111111111111111111111, -0.222222222222222222222222222222, 0.0144444444444444444444444444444, 0, 0, 0, 0, 0, 0, + 2, 0, 0, -8.83333333333333333333333333333, 15.6444444444444444444444444444, -11.8888888888888888888888888889, 0.744444444444444444444444444444, 3, 0, 0, 0, 0, 0, + -0.842592592592592592592592592593, 0, 0, 0.212962962962962962962962962963, -7.22962962962962962962962962963, 5.75925925925925925925925925926, -0.316666666666666666666666666667, 2.83333333333333333333333333333, -0.0833333333333333333333333333333, 0, 0, 0, 0, + 0.581219512195121951219512195122, 0, 0, -2.07926829268292682926829268293, 4.38634146341463414634146341463, -3.67073170731707317073170731707, 0.520243902439024390243902439024, 0.548780487804878048780487804878, 0.274390243902439024390243902439, 0.439024390243902439024390243902, 0, 0, 0, + 0.0146341463414634146341463414634, 0, 0, 0, 0, -0.146341463414634146341463414634, -0.0146341463414634146341463414634, -0.0731707317073170731707317073171, 0.0731707317073170731707317073171, 0.146341463414634146341463414634, 0, 0, 0, + -0.433414634146341463414634146341, 0, 0, -2.07926829268292682926829268293, 4.38634146341463414634146341463, -3.52439024390243902439024390244, 0.534878048780487804878048780488, 0.621951219512195121951219512195, 0.201219512195121951219512195122, 0.292682926829268292682926829268, 0, 1, 0}; + const double b[] = {-0.364851774598758815913721254923, 0.301330765163132508037187549082, 0.8325767564796758519952730731, 0, -0.0268487161606468393107786573864, 0.0959117920564658481117857153644, -0.0286739781279975214853581446154, 0.213522659333043888240291858263, -0.00908604512512093286675314270667, 0.00681453439485641598976617630506, 0.0103479967072222251389783702879, -0.0310439901218725048480067038039, -1.23088664838968206937943710922e-16}; + const double bt[] = { 0, 0, 0, 0, 0, 0.32380952380952380952380952381, 0.257142857142857142857142857143, 0.257142857142857142857142857143, 0.0321428571428571428571428571429, 0.0321428571428571428571428571429, 0, 0.0488095238095238095238095238095, 0.0488095238095238095238095238095}; + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = FALSE; +} + +void getButcherTableau_FEHLBERGSSC2(BUTCHER_TABLEAU* tableau) { + + tableau->nStages = 13; + tableau->order_b = 2; + tableau->order_bt = 8; + tableau->fac = 1e0; + + /* Butcher Tableau */ + const double c[] = { 0, 0.0740740740740740740740740740741, 0.111111111111111111111111111111, 0.166666666666666666666666666667, 0.416666666666666666666666666667, 0.5, 0.833333333333333333333333333333, 0.166666666666666666666666666667, 0.666666666666666666666666666667, 0.333333333333333333333333333333, 1, 0, 1}; + const double A[] = { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.0740740740740740740740740740741, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.0277777777777777777777777777778, 0.0833333333333333333333333333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.0416666666666666666666666666667, 0, 0.125, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.416666666666666666666666666667, 0, -1.5625, 1.5625, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0.05, 0, 0, 0.25, 0.2, 0, 0, 0, 0, 0, 0, 0, 0, + -0.231481481481481481481481481481, 0, 0, 1.15740740740740740740740740741, -2.40740740740740740740740740741, 2.31481481481481481481481481481, 0, 0, 0, 0, 0, 0, 0, + 0.103333333333333333333333333333, 0, 0, 0, 0.271111111111111111111111111111, -0.222222222222222222222222222222, 0.0144444444444444444444444444444, 0, 0, 0, 0, 0, 0, + 2, 0, 0, -8.83333333333333333333333333333, 15.6444444444444444444444444444, -11.8888888888888888888888888889, 0.744444444444444444444444444444, 3, 0, 0, 0, 0, 0, + -0.842592592592592592592592592593, 0, 0, 0.212962962962962962962962962963, -7.22962962962962962962962962963, 5.75925925925925925925925925926, -0.316666666666666666666666666667, 2.83333333333333333333333333333, -0.0833333333333333333333333333333, 0, 0, 0, 0, + 0.581219512195121951219512195122, 0, 0, -2.07926829268292682926829268293, 4.38634146341463414634146341463, -3.67073170731707317073170731707, 0.520243902439024390243902439024, 0.548780487804878048780487804878, 0.274390243902439024390243902439, 0.439024390243902439024390243902, 0, 0, 0, + 0.0146341463414634146341463414634, 0, 0, 0, 0, -0.146341463414634146341463414634, -0.0146341463414634146341463414634, -0.0731707317073170731707317073171, 0.0731707317073170731707317073171, 0.146341463414634146341463414634, 0, 0, 0, + -0.433414634146341463414634146341, 0, 0, -2.07926829268292682926829268293, 4.38634146341463414634146341463, -3.52439024390243902439024390244, 0.534878048780487804878048780488, 0.621951219512195121951219512195, 0.201219512195121951219512195122, 0.292682926829268292682926829268, 0, 1, 0}; + const double b[] = {1.36308696433418349654697814366, 0, -3.26459140897742935132414339163, -1.62051733050260843967905937243, -0.29368516237586603628166052316, 1.8370371692484130907915149173, -0.526247844078595836398397402697, 3.92724305645149417982124010038, -0.167071106248436357738384607217, 0.12530612140680376084342193382, 0.190280224596084454633821394807, -0.570840683836337962747432519784, -1.77049984678986730541092039109e-11}; + const double bt[] = { 0, 0, 0, 0, 0, 0.32380952380952380952380952381, 0.257142857142857142857142857143, 0.257142857142857142857142857143, 0.0321428571428571428571428571429, 0.0321428571428571428571428571429, 0, 0.0488095238095238095238095238095, 0.0488095238095238095238095238095}; + + setButcherTableau(tableau, (double *)c, (double *)A, (double *)b, (double *) bt); + tableau->isKLeftAvailable = TRUE; + tableau->isKRightAvailable = FALSE; +} + // TODO: Describe me void getButcherTableau_RK810(BUTCHER_TABLEAU* tableau) { @@ -1171,6 +1337,12 @@ BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ER case RK_DOPRI45: getButcherTableau_DOPRI45(tableau); break; + case RK_DOPRISSC1: + getButcherTableau_DOPRISSC1(tableau); + break; + case RK_DOPRISSC2: + getButcherTableau_DOPRISSC2(tableau); + break; case RK_RKSSC: getButcherTableau_RKSSC(tableau); break; @@ -1186,6 +1358,12 @@ BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ER case RK_MERSON: getButcherTableau_MERSON(tableau); break; + case RK_MERSONSSC1: + getButcherTableau_MERSONSSC1(tableau); + break; + case RK_MERSONSSC2: + getButcherTableau_MERSONSSC2(tableau); + break; case RK_EXPL_EULER: getButcherTableau_EXPLEULER(tableau); break; @@ -1201,6 +1379,12 @@ BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ER case RK_FEHLBERG78: getButcherTableau_FEHLBERG78(tableau); break; + case RK_FEHLBERGSSC1: + getButcherTableau_FEHLBERGSSC1(tableau); + break; + case RK_FEHLBERGSSC2: + getButcherTableau_FEHLBERGSSC2(tableau); + break; case RK_SDIRK2: getButcherTableau_SDIRK2(tableau); break; diff --git a/OMCompiler/SimulationRuntime/c/util/simulation_options.c b/OMCompiler/SimulationRuntime/c/util/simulation_options.c index d5f0726e3a7..0bddc5219b6 100644 --- a/OMCompiler/SimulationRuntime/c/util/simulation_options.c +++ b/OMCompiler/SimulationRuntime/c/util/simulation_options.c @@ -897,14 +897,20 @@ const char *GB_METHOD_NAME[RK_MAX] = { /* RK_GAUSS5 */ "gauss5", /* RK_GAUSS6 */ "gauss6", /* RK_MERSON */ "merson", + /* RK_MERSONSSC1 */ "mersonSsc1", + /* RK_MERSONSSC2 */ "mersonSsc2", /* RK_HEUN */ "heun", /* RK_FEHLBERG12 */ "fehlberg12", /* RK_FEHLBERG45 */ "fehlberg45", /* RK_FEHLBERG78 */ "fehlberg78", + /* RK_FEHLBERGSSC1 */ "fehlbergSsc1", + /* RK_FEHLBERGSSC2 */ "fehlbergSsc2", /* RK_RK810 */ "rk810", /* RK_RK1012 */ "rk1012", /* RK_RK1214 */ "rk1214", /* RK_DOPRI45 */ "dopri45", + /* RK_DOPRISSC1 */ "dopriSsc1", + /* RK_DOPRISSC2 */ "dopriSsc2", /* RK_RKSSC */ "rungekuttaSsc" }; @@ -937,14 +943,20 @@ 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_MERSONSSC1 */ "Explicit Runge-Kutta Merson method with large stability region (order 1)", + /* RK_MERSONSSC2 */ "Explicit Runge-Kutta Merson method with large stability region (order 2)", /* 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_FEHLBERGSSC1 */ "Explicit Runge-Kutta Fehlberg method with large stability region (order 1)", + /* RK_FEHLBERGSSC2 */ "Explicit Runge-Kutta Fehlberg method with large stability region (order 2)", /* 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_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_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 cc0eab4d4b4..fad60d61567 100644 --- a/OMCompiler/SimulationRuntime/c/util/simulation_options.h +++ b/OMCompiler/SimulationRuntime/c/util/simulation_options.h @@ -229,14 +229,20 @@ enum GB_METHOD { RK_GAUSS5, /* gauss5*/ RK_GAUSS6, /* gauss6*/ RK_MERSON, /* merson*/ + RK_MERSONSSC1, /* mersonSsc1*/ + RK_MERSONSSC2, /* mersonSsc2*/ RK_HEUN, /* heun */ RK_FEHLBERG12, /* fehlberg12*/ RK_FEHLBERG45, /* fehlberg45*/ RK_FEHLBERG78, /* fehlberg78*/ + RK_FEHLBERGSSC1, /* fehlbergSsc1*/ + RK_FEHLBERGSSC2, /* fehlbergSsc2*/ RK_RK810, /* rk810*/ RK_RK1012, /* rk1012*/ RK_RK1214, /* rk1214*/ - RK_DOPRI45, /* dopri4*/ + RK_DOPRI45, /* dopri45*/ + RK_DOPRISSC1, /* dopriSsc1*/ + RK_DOPRISSC2, /* dopriSsc2*/ RK_RKSSC, /* rungekuttaSsc */