Skip to content

Commit

Permalink
Set appropriate method-specific defaults for gbode solver (#10330)
Browse files Browse the repository at this point in the history 
* Set default interpolation to dense_output, if available
* Use single-rate method if percentage = 1
* Set default interpolation for bi-rate mode to dense_output_errctrl
* set default error control to richardson extrapolation for gauss, radau and lobatto RK methods
* Use enumeration and strings for error estimation flags

---------

Co-authored-by: AnHeuermann <AnHeuermann@users.noreply.github.com>
  • Loading branch information
@bernhardbachmann @AnHeuermann
bernhardbachmann and AnHeuermann committed Mar 7, 2023
1 parent 7a6d3ac commit ce22cb1
Show file tree
Hide file tree
Showing 7 changed files with 124 additions and 23 deletions.
62 changes: 59 additions & 3 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.c
View file
Expand Up @@ -176,7 +176,13 @@ enum GB_CTRL_METHOD getControllerMethod(enum _FLAG flag) {
dumOptions(FLAG_NAME[flag], flag_value, GB_CTRL_METHOD_NAME, GB_CTRL_MAX);
return GB_CTRL_UNKNOWN;
} else {
return GB_CTRL_I;
// Default value
if (flag == FLAG_MR_CTRL) {
return getControllerMethod(FLAG_SR_CTRL);
} else {
infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode step size control: i [default]");
return GB_CTRL_I;
}
}
}

Expand Down Expand Up @@ -205,7 +211,7 @@ enum GB_INTERPOL_METHOD getInterpolationMethod(enum _FLAG flag) {
if (strcmp(flag_value, GB_INTERPOL_METHOD_NAME[method]) == 0) {
if (flag == FLAG_MR_INT && (method == GB_INTERPOL_HERMITE_ERRCTRL || method == GB_DENSE_OUTPUT_ERRCTRL)) {
warningStreamPrint(LOG_SOLVER, 0, "Chosen gbode interpolation method %s not supported for fast state integration", GB_INTERPOL_METHOD_NAME[method]);
method = GB_INTERPOL_HERMITE;
method = GB_DENSE_OUTPUT;
}
infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode interpolation method: %s", GB_INTERPOL_METHOD_NAME[method]);
return method;
Expand All @@ -214,7 +220,19 @@ enum GB_INTERPOL_METHOD getInterpolationMethod(enum _FLAG flag) {
dumOptions(FLAG_NAME[flag], flag_value, GB_INTERPOL_METHOD_NAME, GB_INTERPOL_MAX);
return GB_INTERPOL_UNKNOWN;
} else {
return GB_INTERPOL_HERMITE;
// Default value
if (flag == FLAG_MR_INT) {
method = getInterpolationMethod(FLAG_SR_INT);
if (method == GB_INTERPOL_HERMITE_ERRCTRL || method == GB_DENSE_OUTPUT_ERRCTRL) {
warningStreamPrint(LOG_SOLVER, 0, "Chosen gbode interpolation method %s not supported for fast state integration", GB_INTERPOL_METHOD_NAME[method]);
infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode interpolation method: dense_output [default]");
method = GB_DENSE_OUTPUT;
}
return method;
} else {
infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode interpolation method: dense_output [default]");
return GB_DENSE_OUTPUT;
}
}
}

Expand All @@ -241,6 +259,44 @@ double getGBRatio() {
return percentage;
}

/**
* @brief Get extrapolation method from user flag.
*
* Reads flag FLAG_SR_ERR, FLAG_MR_ERR.
* Defaults to GB_EXT_DEFAULT.
*
* @param flag Flag specifying error estimation.
* Allowed values: FLAG_SR_ERR, FLAG_MR_ERR
* @return enum GB_EXTRAPOL_METHOD Extrapolation method.
*/
enum GB_EXTRAPOL_METHOD getGBErr(enum _FLAG flag) {
assertStreamPrint(NULL, flag==FLAG_SR_ERR || flag==FLAG_MR_ERR, "Illegal input 'flag' to getGBErr!");

enum GB_EXTRAPOL_METHOD extrapolationMethod;
const char *flag_value = omc_flagValue[flag];

if (flag_value != NULL) {
if (strcmp(flag_value, "default")==0) {
extrapolationMethod = GB_EXT_DEFAULT;
} else if (strcmp(flag_value, "richardson")==0) {
extrapolationMethod = GB_EXT_RICHARDSON;
} else if (strcmp(flag_value, "embedded")==0) {
extrapolationMethod = GB_EXT_EMBEDDED;
} else {
errorStreamPrint(LOG_STDOUT, 0, "Illegal value '%s' for flag -%s", flag_value, FLAG_NAME[flag]);
infoStreamPrint(LOG_STDOUT, 1, "Allowed values are:");
infoStreamPrint(LOG_STDOUT, 0, "default");
infoStreamPrint(LOG_STDOUT, 0, "richardson");
infoStreamPrint(LOG_STDOUT, 0, "embedded");
messageClose(LOG_STDOUT);
omc_throw(NULL);
}
} else {
extrapolationMethod = GB_EXT_DEFAULT;
}
return extrapolationMethod;
}

/**
* @brief Dump available flag options to stdout.
*
Expand Down
12 changes: 12 additions & 0 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.h
View file
Expand Up @@ -40,11 +40,23 @@
extern "C" {
#endif

/**
* @brief Extrapolation method for single-rate / multi-rate error estimation.
*/
enum GB_EXTRAPOL_METHOD{
GB_EXT_UNKNOWN = 0, /* Unknown method */

GB_EXT_DEFAULT, /* Default, depending on the Runge-Kutta method */
GB_EXT_RICHARDSON, /* Richardson extrapolation*/
GB_EXT_EMBEDDED /* Embedded scheme */
};

enum GB_METHOD getGB_method(enum _FLAG flag);
enum GB_INTERPOL_METHOD getInterpolationMethod(enum _FLAG flag);
enum GB_CTRL_METHOD getControllerMethod(enum _FLAG flag);
enum GB_NLS_METHOD getGB_NLS_method(enum _FLAG flag);
double getGBRatio();
enum GB_EXTRAPOL_METHOD getGBErr(enum _FLAG flag);

#ifdef __cplusplus
};
Expand Down
14 changes: 10 additions & 4 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c
View file
Expand Up @@ -264,7 +264,7 @@ int gbodef_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solve
break;
case GB_DENSE_OUTPUT:
case GB_DENSE_OUTPUT_ERRCTRL:
infoStreamPrint(LOG_SOLVER, 0, "If available, dense output is used for emitting results, otherwise hermite");
infoStreamPrint(LOG_SOLVER, 0, "Dense output is used for emitting results");
break;
default:
throwStreamPrint(NULL, "Unhandled interpolation case.");
Expand Down Expand Up @@ -467,7 +467,7 @@ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solver
}

gbData->percentage = getGBRatio();
gbData->multi_rate = gbData->percentage > 0;
gbData->multi_rate = gbData->percentage > 0 && gbData->percentage < 1;

gbData->fastStatesIdx = malloc(sizeof(int) * gbData->nStates);
gbData->slowStatesIdx = malloc(sizeof(int) * gbData->nStates);
Expand All @@ -481,7 +481,13 @@ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solver
gbData->sortedStatesIdx[i] = i;
}

gbData->interpolation = getInterpolationMethod(FLAG_SR_INT);

if (gbData->multi_rate && omc_flagValue[FLAG_SR_INT]==NULL) {
gbData->interpolation = GB_DENSE_OUTPUT_ERRCTRL;
} else {
gbData->interpolation = getInterpolationMethod(FLAG_SR_INT);
}

if (!gbData->tableau->withDenseOutput) {
if (gbData->interpolation == GB_DENSE_OUTPUT) gbData->interpolation = GB_INTERPOL_HERMITE;
if (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL) gbData->interpolation = GB_INTERPOL_HERMITE_ERRCTRL;
Expand All @@ -507,7 +513,7 @@ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solver
break;
case GB_DENSE_OUTPUT:
case GB_DENSE_OUTPUT_ERRCTRL:
infoStreamPrint(LOG_SOLVER, 0, "If available, dense output is used for emitting results%s", buffer);
infoStreamPrint(LOG_SOLVER, 0, "Dense output is used for emitting results%s", buffer);
break;
default:
throwStreamPrint(NULL, "Unhandled interpolation case.");
Expand Down
37 changes: 29 additions & 8 deletions OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c
View file
Expand Up @@ -34,6 +34,7 @@
*/

#include "gbode_tableau.h"
#include "gbode_conf.h"

#include <string.h>

Expand Down Expand Up @@ -1393,23 +1394,26 @@ void analyseButcherTableau(BUTCHER_TABLEAU* tableau, int nStates, unsigned int*
/**
* @brief Allocate memory and initialize Butcher tableau for given method.
*
* @param GM_method Runge-Kutta method.
* @param method Runge-Kutta method.
* @param flag Flag specifying error estimation.
* Allowed values: FLAG_SR_ERR, FLAG_MR_ERR
* @return BUTCHER_TABLEAU* Return pointer to Butcher tableau on success, NULL on failure.
*/
BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ERR) {
BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD method, enum _FLAG flag) {
BUTCHER_TABLEAU* tableau = (BUTCHER_TABLEAU*) malloc(sizeof(BUTCHER_TABLEAU));
const char* flag_value = omc_flagValue[FLAG_ERR];
int richardson = 0;
enum GB_EXTRAPOL_METHOD extrapolMethod;

if (flag_value != NULL) richardson = atoi(flag_value);
if (richardson) {
assertStreamPrint(NULL, flag==FLAG_SR_ERR || flag==FLAG_MR_ERR, "Illegal input 'flag' to initButcherTableau!");

extrapolMethod = getGBErr(flag);
if (extrapolMethod == GB_EXT_RICHARDSON) {
tableau->richardson = TRUE;
infoStreamPrint(LOG_SOLVER, 0, "Richardson extrapolation is used for step size control");
} else {
tableau->richardson = FALSE;
}

switch(GM_method)
switch(method)
{
case MS_ADAMS_MOULTON:
getButcherTableau_MS(tableau);
Expand Down Expand Up @@ -1493,58 +1497,75 @@ BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ER
getButcherTableau_ESDIRK4(tableau);
break;
case RK_RADAU_IA_2:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_RADAU_IA_2(tableau);
break;
case RK_RADAU_IA_3:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_RADAU_IA_3(tableau);
break;
case RK_RADAU_IA_4:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_RADAU_IA_4(tableau);
break;
case RK_RADAU_IIA_2:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_RADAU_IIA_2(tableau);
break;
case RK_RADAU_IIA_3:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_RADAU_IIA_3(tableau);
break;
case RK_RADAU_IIA_4:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_RADAU_IIA_4(tableau);
break;
case RK_LOBA_IIIA_3:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_LOBATTO_IIIA_3(tableau);
break;
case RK_LOBA_IIIA_4:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_LOBATTO_IIIA_4(tableau);
break;
case RK_LOBA_IIIB_3:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_LOBATTO_IIIB_3(tableau);
break;
case RK_LOBA_IIIB_4:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_LOBATTO_IIIB_4(tableau);
break;
case RK_LOBA_IIIC_3:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_LOBATTO_IIIC_3(tableau);
break;
case RK_LOBA_IIIC_4:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_LOBATTO_IIIC_4(tableau);
break;
case RK_GAUSS2:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_GAUSS2(tableau);
break;
case RK_GAUSS3:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_GAUSS3(tableau);
break;
case RK_GAUSS4:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_GAUSS4(tableau);
break;
case RK_GAUSS5:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_GAUSS5(tableau);
break;
case RK_GAUSS6:
if (extrapolMethod == GB_EXT_DEFAULT) tableau->richardson = TRUE;
getButcherTableau_GAUSS6(tableau);
break;
default:
throwStreamPrint(NULL, "Error: Unknown Runge Kutta method %i.", GM_method);
throwStreamPrint(NULL, "Error: Unknown Runge Kutta method %i.", method);
}

return tableau;
Expand Down
2 changes: 1 addition & 1 deletion OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.h
View file
Expand Up @@ -99,7 +99,7 @@ enum GM_TYPE {

/* Function prototypes */

BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD GM_method, enum _FLAG FLAG_ERR);
BUTCHER_TABLEAU* initButcherTableau(enum GB_METHOD method, enum _FLAG flag);
void freeButcherTableau(BUTCHER_TABLEAU* tableau);

void analyseButcherTableau(BUTCHER_TABLEAU* tableau, int nStates, unsigned int* nlSystemSize, enum GM_TYPE* expl);
Expand Down
18 changes: 12 additions & 6 deletions OMCompiler/SimulationRuntime/c/util/simulation_options.c
View file
Expand Up @@ -273,12 +273,12 @@ const char *FLAG_DESC[FLAG_MAX+1] = {
/* FLAG_DATA_RECONCILE_STATE */ "Run the State Estimation numerical computation algorithm for constrained equations",
/* FLAG_SR */ "Value specifies the chosen solver of solver gbode (single-rate, slow states integrator)",
/* FLAG_SR_CTRL */ "Step size control of solver gbode (single-rate, slow states integrator)",
/* FLAG_SR_ERR */ "Error estimation done by Richardson extrapolation (-gberr=1) of solver gbode (single-rate, slow states integrator)",
/* FLAG_SR_ERR */ "Error estimation method for solver gbode (single-rate, slow states integrator).",
/* FLAG_SR_INT */ "Interpolation method of solver gbode (single-rate, slow states integrator)",
/* FLAG_SR_NLS */ "Non-linear solver method of solver gbode (single-rate, slow states integrator)",
/* FLAG_MR */ "Value specifies the chosen solver of solver gbode (multi-rate, fast states integrator)",
/* FLAG_MR_CTRL */ "Step size control of solver gbode (multi-rate, fast states integrator)",
/* FLAG_MR_ERR */ "Error estimation done by Richardson extrapolation (-gberr=1) of solver gbode (multi-rate, fast states integrator)",
/* FLAG_MR_ERR */ "Error estimation method for gbode solver (multi-rate, fast states integrator).",
/* FLAG_MR_INT */ "Interpolation method of solver gbode (multi-rate, fast states integrator)",
/* FLAG_MR_NLS */ "Non-linear solver method of solver gbode (multi-rate, fast states integrator)",
/* FLAG_MR_PAR */ "Define percentage of states for the fast states selection of solver gbode",
Expand Down Expand Up @@ -579,8 +579,11 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
/* FLAG_SR_CTRL */
" Step size control of solver gbode (single-rate, slow states integrator).",
/* FLAG_SR_ERR */
" Error estimation done by Richardson extrapolation (-gberr=1) of solver gbode\n"
" (single-rate, slow states integrator).",
" Error estimation method for solver gbode (single-rate, slow states integrator)\n"
" Possible values:\n\n"
" * default - depending on the Runge-Kutta method\n"
" * richardson - Richardson extrapolation\n"
" * embedded - Embedded scheme\n",
/* FLAG_SR_INT */
" Interpolation method of solver gbode (single-rate, slow states integrator).",
/* FLAG_SR_NLS */
Expand All @@ -591,8 +594,11 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
/* FLAG_MR_CTRL */
" Step size control of solver gbode (multi-rate, fast states integrator).",
/* FLAG_MR_ERR */
" Error estimation done by Richardson extrapolation (-gberr=1) of solver gbode\n"
" (multi-rate, fast states integrator).",
" Error estimation method for solver gbode (multi-rate, fast states integrator)\n"
" Possible values:\n\n"
" * default - depending on the Runge-Kutta method\n"
" * richardson - Richardson extrapolation\n"
" * embedded - Embedded scheme\n",
/* FLAG_MR_INT */
" Interpolation method of solver gbode (multi-rate, fast states integrator).",
/* FLAG_MR_NLS */
Expand Down
2 changes: 1 addition & 1 deletion doc/UsersGuide/source/solving.rst
View file
Expand Up @@ -138,7 +138,7 @@ adjust the behavior of the solver for specific simulation problems:
:ref:`gbfctrl <simflag-gbctrl>`,
:ref:`gbfnls <simflag-gbnls>`,
:ref:`gbfint <simflag-gbint>`,
:ref:`gbferr <simflag-gberr>`.
:ref:`gbferr <simflag-gbferr>`.

This solver will replace obsolete and no longer maintained solvers providing a
lot more using the following simulation flags:
Expand Down

0 comments on commit ce22cb1

Please sign in to comment.