diff --git a/SimulationRuntime/c/simulation/solver/radau.c b/SimulationRuntime/c/simulation/solver/radau.c index e5ae3a1e411..c3bb95f0cf4 100644 --- a/SimulationRuntime/c/simulation/solver/radau.c +++ b/SimulationRuntime/c/simulation/solver/radau.c @@ -29,10 +29,14 @@ */ /*! \file radau.c - * author: team Bielefeld :) + * author: Ru(n)ge + * description: This file contains implicit Runge-Kutta methods of different order, + * based on Radau IIA and Lobatto IIA methods. The method with order one is + * corresponding to the implicit euler method. Further orders 2 to 6. */ #include +#include #include "radau.h" #include "external_input.h" @@ -55,7 +59,7 @@ extern int KINSetInfoHandlerFn(void *kinmem, KINInfoHandlerFn ihfun, void *ih_da } #endif -static int allocateNlpOde(KINODE *kinOde); +static int allocateNlpOde(KINODE *kinOde, int order); static int allocateKINSOLODE(KINODE *kinOde); @@ -81,16 +85,16 @@ static int lobatto2Res(N_Vector z, N_Vector f, void* user_data); static int lobatto4Res(N_Vector z, N_Vector f, void* user_data); static int lobatto6Res(N_Vector z, N_Vector f, void* user_data); -int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, int flag, int N) +int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, int order) { KINODE *kinOde = (KINODE*) solverInfo->solverData; kinOde->kData = (KDATAODE*) malloc(sizeof(KDATAODE)); kinOde->nlp = (NLPODE*) malloc(sizeof(NLPODE)); - kinOde->N = N; - kinOde->flag = flag; + kinOde->order = order; + kinOde->N = ceil((double)order/2.0); kinOde->data = data; kinOde->threadData = threadData; - allocateNlpOde(kinOde); + allocateNlpOde(kinOde, order); allocateKINSOLODE(kinOde); kinOde->solverInfo = solverInfo; kinOde->kData->mset = 50; @@ -104,29 +108,29 @@ int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo if (ACTIVE_STREAM(LOG_SOLVER)) { KINSetPrintLevel(kinOde->kData->kmem,2); } - //KINSetEtaForm(kinOde->kData->kmem, KIN_ETACHOICE2); + /* KINSetEtaForm(kinOde->kData->kmem, KIN_ETACHOICE2); */ KINSetMaxSetupCalls(kinOde->kData->kmem, kinOde->kData->mset); kinOde->nlp->currentStep = &kinOde->solverInfo->currentStepSize; KINSetErrHandlerFn(kinOde->kData->kmem, kinsol_errorHandler, NULL); KINSetInfoHandlerFn(kinOde->kData->kmem, kinsol_infoHandler, NULL); - switch(kinOde->flag) + switch(kinOde->order) { - case S_RADAU5: + case 5: KINInit(kinOde->kData->kmem, radau5Res, kinOde->kData->x); break; - case S_RADAU3: + case 3: KINInit(kinOde->kData->kmem, radau3Res, kinOde->kData->x); break; - case S_RADAU1: + case 1: KINInit(kinOde->kData->kmem, radau1Res, kinOde->kData->x); break; - case S_LOBATTO2: + case 2: KINInit(kinOde->kData->kmem, lobatto2Res, kinOde->kData->x); break; - case S_LOBATTO4: + case 4: KINInit(kinOde->kData->kmem, lobatto4Res, kinOde->kData->x); break; - case S_LOBATTO6: + case 6: KINInit(kinOde->kData->kmem, lobatto6Res, kinOde->kData->x); break; @@ -136,37 +140,37 @@ int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo /* Call KINDense to specify the linear solver */ /*KINSpbcg*/ if(kinOde->nlp->nStates < 10) - KINSpgmr(kinOde->kData->kmem, N*kinOde->nlp->nStates+1); + KINSpgmr(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates+1); else - KINSpbcg(kinOde->kData->kmem, N*kinOde->nlp->nStates+1); + KINSpbcg(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates+1); kinOde->kData->glstr = KIN_LINESEARCH; return 0; } -static int allocateNlpOde(KINODE *kinOde) +static int allocateNlpOde(KINODE *kinOde, int order) { NLPODE * nlp = (NLPODE*) kinOde->nlp; nlp->nStates = kinOde->data->modelData->nStates; - switch(kinOde->flag) + switch(order) { - case S_RADAU5: + case 5: radau5Coeff(kinOde); break; - case S_RADAU3: + case 3: radau3Coeff(kinOde); break; - case S_RADAU1: + case 1: radau1Coeff(kinOde); break; - case S_LOBATTO2: + case 2: radau1Coeff(kinOde); /* TODO: Is this right? */ break; - case S_LOBATTO4: + case 4: lobatto4Coeff(kinOde); break; - case S_LOBATTO6: + case 6: lobatto6Coeff(kinOde); break; default: @@ -214,14 +218,13 @@ static int boundsVars(KINODE *kinOde) static int radau5Coeff(KINODE *kinOde) { int i; - const int N = 3; NLPODE * nlp = (NLPODE*) kinOde->nlp; - nlp->c = (long double**) malloc(N * sizeof(long double*)); - for(i = 0; i < N; i++) - nlp->c[i] = (long double*) calloc(N+1, sizeof(long double)); + nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*)); + for(i = 0; i < kinOde->N; i++) + nlp->c[i] = (long double*) calloc(kinOde->N+1, sizeof(long double)); - nlp->a = (double*) malloc(N * sizeof(double)); + nlp->a = (double*) malloc(kinOde->N * sizeof(double)); nlp->c[0][0] = 4.1393876913398137178367408896470696703591369767880; nlp->c[0][1] = 3.2247448713915890490986420373529456959829737403284; @@ -247,14 +250,13 @@ static int radau5Coeff(KINODE *kinOde) static int radau3Coeff(KINODE *kinOde) { int i; - const int N = 2; NLPODE * nlp = (NLPODE*) kinOde->nlp; - nlp->c = (long double**) malloc(N * sizeof(long double*)); - for(i = 0; i < N; i++) - nlp->c[i] = (long double*) calloc(N+1, sizeof(long double)); + nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*)); + for(i = 0; i < kinOde->N; i++) + nlp->c[i] = (long double*) calloc(kinOde->N+1, sizeof(long double)); - nlp->a = (double*) malloc(N * sizeof(double)); + nlp->a = (double*) malloc(kinOde->N * sizeof(double)); nlp->c[0][0] = 2.0; nlp->c[0][1] = 1.50; @@ -294,15 +296,14 @@ static int lobatto4Coeff(KINODE *kinOde) static int lobatto6Coeff(KINODE *kinOde) { int i; - const int N = 3; NLPODE * nlp = (NLPODE*) kinOde->nlp; - nlp->c = (long double**) malloc(N * sizeof(long double*)); - for(i = 0; i < N; ++i) - nlp->c[i] = (long double*) malloc((N+2)* sizeof(long double)); + nlp->c = (long double**) malloc(kinOde->N * sizeof(long double*)); + for(i = 0; i < kinOde->N; ++i) + nlp->c[i] = (long double*) malloc((kinOde->N+2)* sizeof(long double)); - nlp->a = (double*) malloc(N * sizeof(double)); + nlp->a = (double*) malloc(kinOde->N * sizeof(double)); nlp->c[0][0] = 4.3013155617496424838955952368431696002490512113396; nlp->c[0][1] = 3.6180339887498948482045868343656381177203091798058; @@ -353,10 +354,10 @@ static int allocateKINSOLODE(KINODE *kinOde) return 0; } -int freeKinOde(DATA* data, SOLVER_INFO* solverInfo, int N) +int freeKinOde(DATA* data, SOLVER_INFO* solverInfo) { KINODE *kinOde = (KINODE*) solverInfo->solverData; - freeImOde((void*) kinOde->nlp ,N); + freeImOde((void*) kinOde->nlp, kinOde->N); freeKinsol((void*) kinOde->kData); free(kinOde); return 0; @@ -392,7 +393,7 @@ static int freeKinsol(void * kOde) static int initKinsol(KINODE *kinOde) { - int i,j,k, n; + int i,j,k,n; double *scal_eq, *scal_var, *x, *f2; long double tmp, h, hf, hf_min; DATA *data; diff --git a/SimulationRuntime/c/simulation/solver/radau.h b/SimulationRuntime/c/simulation/solver/radau.h index d87ac9f2258..9393b80718f 100644 --- a/SimulationRuntime/c/simulation/solver/radau.h +++ b/SimulationRuntime/c/simulation/solver/radau.h @@ -40,6 +40,9 @@ #include "omc_config.h" #ifdef WITH_SUNDIALS + + #define DEFAULT_IMPRK_ORDER 5 + #include #include @@ -84,7 +87,7 @@ threadData_t *threadData; SOLVER_INFO *solverInfo; int N; - int flag; + int order; }KINODE; #else @@ -94,12 +97,12 @@ DATA *data; SOLVER_INFO *solverInfo; int N; - int flag; + int order; }KINODE; #endif /* SUNDIALS */ - int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, int flag, int N); - int freeKinOde(DATA* data, SOLVER_INFO* solverInfo, int N); + int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, int order); + int freeKinOde(DATA* data, SOLVER_INFO* solverInfo); int kinsolOde(SOLVER_INFO* solverInfo); #ifdef __cplusplus }; diff --git a/SimulationRuntime/c/simulation/solver/solver_main.c b/SimulationRuntime/c/simulation/solver/solver_main.c index 67cce6be188..b296917c2b7 100644 --- a/SimulationRuntime/c/simulation/solver/solver_main.c +++ b/SimulationRuntime/c/simulation/solver/solver_main.c @@ -139,12 +139,9 @@ int solver_main_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverIn return retVal; #endif #ifdef WITH_SUNDIALS - case S_RADAU5: - case S_RADAU3: case S_RADAU1: case S_LOBATTO2: - case S_LOBATTO4: - case S_LOBATTO6: + case S_IMPRUNGEKUTTA: retVal = radau_lobatto_step(data, solverInfo); if(omc_flag[FLAG_SOLVER_STEPS]) data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0]; @@ -281,58 +278,37 @@ int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solv } #endif #ifdef WITH_SUNDIALS - case S_RADAU5: - { - /* Allocate Radau5 IIA work arrays */ - infoStreamPrint(LOG_SOLVER, 0, "Initializing Radau IIA of order 5"); - solverInfo->solverData = calloc(1, sizeof(KINODE)); - allocateKinOde(data, threadData, solverInfo, solverInfo->solverMethod, 3); - break; - } - case S_RADAU3: - { - /* Allocate Radau3 IIA work arrays */ - infoStreamPrint(LOG_SOLVER, 0, "Initializing Radau IIA of order 3"); - solverInfo->solverData = calloc(1, sizeof(KINODE)); - allocateKinOde(data, threadData, solverInfo, solverInfo->solverMethod, 2); - break; - } case S_RADAU1: - { - /* Allocate Radau1 IIA work arrays */ - infoStreamPrint(LOG_SOLVER, 0, "Initializing Radau IIA of order 1 (implicit euler) "); - solverInfo->solverData = calloc(1, sizeof(KINODE)); - allocateKinOde(data, threadData, solverInfo, solverInfo->solverMethod, 1); - break; - } - case S_LOBATTO6: - { - /* Allocate Lobatto2 IIIA work arrays */ - infoStreamPrint(LOG_SOLVER, 0, "Initializing Lobatto IIIA of order 6"); - solverInfo->solverData = calloc(1, sizeof(KINODE)); - allocateKinOde(data, threadData, solverInfo, solverInfo->solverMethod, 3); - break; - } - case S_LOBATTO4: - { - /* Allocate Lobatto4 IIIA work arrays */ - infoStreamPrint(LOG_SOLVER, 0, "Initializing Lobatto IIIA of order 4"); - solverInfo->solverData = calloc(1, sizeof(KINODE)); - allocateKinOde(data, threadData, solverInfo, solverInfo->solverMethod, 2); - break; - } case S_LOBATTO2: + case S_IMPRUNGEKUTTA: { - /* Allocate Lobatto6 IIIA work arrays */ - infoStreamPrint(LOG_SOLVER, 0, "Initializing Lobatto IIIA of order 2 (trapeze rule)"); + int usedImpRKOrder = DEFAULT_IMPRK_ORDER; + if (solverInfo->solverMethod == S_RADAU1) + usedImpRKOrder = 1; + if (solverInfo->solverMethod == S_LOBATTO2) + usedImpRKOrder = 2; + + /* Check the order if set */ + if (omc_flag[FLAG_IMPRK_ORDER]) + { + usedImpRKOrder = atoi(omc_flagValue[FLAG_IMPRK_ORDER]); + if (usedImpRKOrder>6 || usedImpRKOrder<1) + { + warningStreamPrint(LOG_STDOUT, 0, "Selected order %d is out of range[1-6]. Use default order %d", usedImpRKOrder, DEFAULT_IMPRK_ORDER); + usedImpRKOrder = DEFAULT_IMPRK_ORDER; + } + } + + /* Allocate implicit Runge-Kutta methods */ + infoStreamPrint(LOG_SOLVER, 0, "Initializing Runge-Kutta method with order %d", usedImpRKOrder); solverInfo->solverData = calloc(1, sizeof(KINODE)); - allocateKinOde(data, threadData, solverInfo, solverInfo->solverMethod, 1); + allocateKinOde(data, threadData, solverInfo, usedImpRKOrder); break; } case S_IDA: { IDA_SOLVER* idaData = NULL; - /* Allocate Lobatto6 IIIA work arrays */ + /* Allocate ida working data */ infoStreamPrint(LOG_SOLVER, 0, "Initializing IDA DAE Solver"); idaData = (IDA_SOLVER*) malloc(sizeof(IDA_SOLVER)); retValue = ida_solver_initial(data, threadData, solverInfo, idaData); @@ -394,35 +370,12 @@ int freeSolverData(DATA* data, SOLVER_INFO* solverInfo) } #endif #ifdef WITH_SUNDIALS - else if(solverInfo->solverMethod == S_RADAU5) - { - /* free work arrays */ - freeKinOde(data, solverInfo, 3); - } - else if(solverInfo->solverMethod == S_RADAU3) - { - /* free work arrays */ - freeKinOde(data, solverInfo, 2); - } - else if(solverInfo->solverMethod == S_RADAU1) - { - /* free work arrays */ - freeKinOde(data, solverInfo, 1); - } - else if(solverInfo->solverMethod == S_LOBATTO6) - { - /* free work arrays */ - freeKinOde(data, solverInfo, 3); - } - else if(solverInfo->solverMethod == S_LOBATTO4) - { - /* free work arrays */ - freeKinOde(data, solverInfo, 2); - } - else if(solverInfo->solverMethod == S_LOBATTO2) + else if(solverInfo->solverMethod == S_RADAU1 || + solverInfo->solverMethod == S_LOBATTO2 || + solverInfo->solverMethod == S_IMPRUNGEKUTTA) { /* free work arrays */ - freeKinOde(data, solverInfo, 1); + freeKinOde(data, solverInfo); } else if(solverInfo->solverMethod == S_IDA) { diff --git a/SimulationRuntime/c/util/simulation_options.c b/SimulationRuntime/c/util/simulation_options.c index 020b5ff2756..e24b191fb37 100644 --- a/SimulationRuntime/c/util/simulation_options.c +++ b/SimulationRuntime/c/util/simulation_options.c @@ -56,6 +56,7 @@ const char *FLAG_NAME[FLAG_MAX+1] = { /* FLAG_IIM */ "iim", /* FLAG_IIT */ "iit", /* FLAG_ILS */ "ils", + /* FLAG_IMPRK_ORDER */ "impRKOrder", /* FLAG_INITIAL_STEP_SIZE */ "initialStepSize", /* FLAG_INPUT_CSV */ "csvInput", /* FLAG_INPUT_FILE */ "exInputFile", @@ -139,6 +140,7 @@ const char *FLAG_DESC[FLAG_MAX+1] = { /* FLAG_IIM */ "value specifies the initialization method", /* FLAG_IIT */ "[double] value specifies a time for the initialization of the model", /* FLAG_ILS */ "[int] default: 1", + /* FLAG_IMPRK_ORDER */ "[int (default 5)] value specifies the integration order of the implicit Runge-Kutta method. Valid values: 1-6", /* FLAG_INITIAL_STEP_SIZE */ "value specifies an initial step size for supported solver", /* FLAG_INPUT_CSV */ "value specifies an csv-file with inputs for the simulation/optimization of the model", /* FLAG_INPUT_FILE */ "value specifies an external file with inputs for the simulation/optimization of the model", @@ -258,6 +260,8 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = { /* FLAG_ILS */ " Value specifies the number of steps for homotopy method (required: -iim=symbolic).\n" " The value is an Integer with default value 1.", + /* FLAG_IMPRK_ORDER */ + " Value specifies the integration order of the implicit Runge-Kutta method. Valid values: 1 to 6. Default order is 5.", /* FLAG_INITIAL_STEP_SIZE */ " Value specifies an initial step size, used by the methods: dassl, ida", /* FLAG_INPUT_CSV */ @@ -440,6 +444,7 @@ const int FLAG_TYPE[FLAG_MAX] = { /* FLAG_IIM */ FLAG_TYPE_OPTION, /* FLAG_IIT */ FLAG_TYPE_OPTION, /* FLAG_ILS */ FLAG_TYPE_OPTION, + /* FLAG_IMPRK_ORDER */ FLAG_TYPE_OPTION, /* FLAG_INITIAL_STEP_SIZE */ FLAG_TYPE_OPTION, /* FLAG_INPUT_CSV */ FLAG_TYPE_OPTION, /* FLAG_INPUT_FILE */ FLAG_TYPE_OPTION, @@ -501,12 +506,9 @@ const char *SOLVER_METHOD_NAME[S_MAX] = { "rungekutta", "dassl", "optimization", - "radau5", - "radau3", "impeuler", "trapezoid", - "lobatto4", - "lobatto6", + "imprungekutta", "symEuler", "symEulerSsc", "heun", @@ -521,12 +523,9 @@ const char *SOLVER_METHOD_DESC[S_MAX] = { "rungekutta - Runge-Kutta (fixed step, order 4)", "dassl - BDF solver with colored numerical Jacobian, with interval root finding - default", "optimization - Special solver for dynamic optimization", - "radau5 - Radau IIA with 3 points, \"Implicit Runge-Kutta\", order 5 [sundial/kinsol needed]", - "radau3 - Radau IIA with 2 points, \"Implicit Runge-Kutta\", order 3 [sundial/kinsol needed]", "impeuler - Implicit Euler (actually Radau IIA, order 1) [sundial/kinsol needed]", "trapezoid - Trapezoidal rule (actually Lobatto IIA with 2 points) [sundial/kinsol needed]", - "lobatto4 - Lobatto IIA with 3 points, order 4 [sundial/kinsol needed]", - "lobatto6 - Lobatto IIA with 4 points, order 6 [sundial/kinsol needed]", + "imprungekutta - Implicit Runge-Kutta bases on Radau and Lobatto IIA methods. Order based on flag 1-6 [sundial/kinsol needed]", "symEuler - symbolic implicit euler, [compiler flag +symEuler needed]", "symEulerSsc - symbolic implicit euler with step-size control, [compiler flag +symEuler needed]", "heun - Heun's method (Runge-Kutta fixed step, order 2)", diff --git a/SimulationRuntime/c/util/simulation_options.h b/SimulationRuntime/c/util/simulation_options.h index 345fa611fa6..1fe97929b29 100644 --- a/SimulationRuntime/c/util/simulation_options.h +++ b/SimulationRuntime/c/util/simulation_options.h @@ -64,6 +64,7 @@ enum _FLAG FLAG_IIM, FLAG_IIT, FLAG_ILS, + FLAG_IMPRK_ORDER, FLAG_INITIAL_STEP_SIZE, FLAG_INPUT_CSV, FLAG_INPUT_FILE, @@ -140,21 +141,18 @@ enum SOLVER_METHOD { S_UNKNOWN = 0, - S_EULER, /* 1 */ - S_RUNGEKUTTA, /* 2 */ - S_DASSL, /* 3 */ - S_OPTIMIZATION, /* 4 */ - S_RADAU5, /* 5 */ - S_RADAU3, /* 6 */ - S_RADAU1, /* 7 */ - S_LOBATTO2, /* 8 */ - S_LOBATTO4, /* 9 */ - S_LOBATTO6, /* 10 */ - S_SYM_EULER, /* 11 */ - S_SYM_IMP_EULER, /* 12 */ - S_HEUN, /* 13 */ - S_IDA, /* 14 */ - S_ERKSSC, /* 15 */ + S_EULER, + S_RUNGEKUTTA, + S_DASSL, + S_OPTIMIZATION, + S_RADAU1, + S_LOBATTO2, + S_IMPRUNGEKUTTA, + S_SYM_EULER, + S_SYM_IMP_EULER, + S_HEUN, + S_IDA, + S_ERKSSC, S_QSS, S_MAX