added explicit rungekutta with step size control of order 4

based on:
2016 9th EUROSIM Congress on Modelling and Simulation,
Solving Stiff Systems of ODEs by Explicit Methods with Conformed Stability
Domains,
A. E. Novikov, A. E. Novikov, A. E. Novikov, A. E. Novikov

ToDo:
- terms for stability error missing in the ste size control
- added internal event handling (??)

Author: vruge

SimulationRuntime/c/simulation/solver/solver_main.c

@@ -77,10 +77,12 @@ typedef struct RK4_DATA
   int work_states_ndims;
   double *b;
   double *c;
+	double h;
 }RK4_DATA;
 
 
 static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo);
+static int rungekutta_step_ssc(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
 static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
 static int sym_euler_im_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo);
 
@@ -154,6 +156,12 @@ int solver_main_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverIn
     TRACE_POP
     return retVal;
 #endif
+  case S_ERKSSC:
+    retVal = rungekutta_step_ssc(data, threadData, solverInfo);
+    if(omc_flag[FLAG_SOLVER_STEPS])
+      data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
+    TRACE_POP
+    return retVal;
   case S_SYM_EULER:
     retVal = sym_euler_im_step(data, threadData, solverInfo);
     if(omc_flag[FLAG_SOLVER_STEPS])
@@ -214,6 +222,7 @@ int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solv
     allocateSymEulerImp(solverInfo, data->modelData->nStates);
     break;
   }
+  case S_ERKSSC:
   case S_RUNGEKUTTA:
   case S_HEUN:
   {
@@ -237,6 +246,9 @@ int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solv
       rungeData->work_states_ndims = heun_s;
       rungeData->b = heun_b;
       rungeData->c = heun_c;
+    } else if (solverInfo->solverMethod==S_ERKSSC) {
+      rungeData->h = (omc_flag[FLAG_INITIAL_STEP_SIZE]) ? atof(omc_flagValue[FLAG_INITIAL_STEP_SIZE]) : solverInfo->currentStepSize;
+      rungeData->work_states_ndims = 5;
     } else {
       throwStreamPrint(threadData, "Unknown RK solver");
     }
@@ -357,7 +369,8 @@ int freeSolverData(DATA* data, SOLVER_INFO* solverInfo)
   {
     freeSymEulerImp(solverInfo);
   }
-  else if(solverInfo->solverMethod == S_RUNGEKUTTA || solverInfo->solverMethod == S_HEUN)
+  else if(solverInfo->solverMethod == S_RUNGEKUTTA || solverInfo->solverMethod == S_HEUN
+				 || solverInfo->solverMethod == S_ERKSSC)
   {
     /* free RK work arrays */
     for(i = 0; i < ((RK4_DATA*)(solverInfo->solverData))->work_states_ndims + 1; i++)
@@ -802,6 +815,127 @@ static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo)
   return 0;
 }
 
+/***************************************   EXP_RK_SSC     *********************************/
+static int rungekutta_step_ssc(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
+{
+  // see.
+	// Solving Stiff Systems of ODEs by Explicit Methods with Conformed Stability
+	// Domains
+	// 2016 9th EUROSIM Congress on Modelling and Simulation
+	// A. E. Novikov...
+
+  RK4_DATA *rk = ((RK4_DATA*)(solverInfo->solverData));
+	const double Atol = data->simulationInfo->tolerance;
+  const double Rtol = data->simulationInfo->tolerance;
+  const int nx = data->modelData->nStates;
+  modelica_real h = rk->h;
+  double** k = rk->work_states;
+  int j, i;
+	double sum;
+  SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
+  SIMULATION_DATA *sDataOld = (SIMULATION_DATA*)data->localData[1];
+  modelica_real* stateDer = sData->realVars + nx;
+  modelica_real* stateDerOld = sDataOld->realVars + nx;
+	double t = sDataOld->timeValue;
+  const double targetTime = t + solverInfo->currentStepSize;
+	const short isMaxStepSizeSet = (short) omc_flagValue[FLAG_MAX_STEP_SIZE];
+	const double maxStepSize = isMaxStepSizeSet ? atof(omc_flagValue[FLAG_MAX_STEP_SIZE]) : -1;
+	double x0[nx];
+
+  if(t + h > targetTime)
+	  h = solverInfo->currentStepSize;
+
+  memcpy(k[0], stateDerOld, nx*sizeof(double));
+  memcpy(x0, sDataOld->realVars, nx*sizeof(double));
+	while(t < targetTime && h > 0){
+		//printf("\nfrom %g to %g by %g\n", t, targetTime, h);
+	  for(j = 0; j < 5; ++j){
+
+			if(j > 0)
+        memcpy(k[j], stateDer, nx*sizeof(double));
+
+      switch(j){
+         case 0:
+	       //yn + k1/3
+	         for(i = 0; i < nx; ++i)
+             sData->realVars[i] = x0[i] + h * k[0][i]/3.0;
+           sData->timeValue = t + h/3.0;
+				   break;
+
+			   case 1:
+			   //yn + 1/6(k1 + k2)
+	         for(i = 0; i < nx; ++i)
+             sData->realVars[i] = x0[i] + h/6.0 * (k[0][i] + k[1][i]);
+           sData->timeValue = t + h/3.0;
+				   break;
+
+			   case 2:
+			   //yn + 1/8*(k1 + 3*k3)
+	         for(i = 0; i < nx; ++i)
+             sData->realVars[i] = x0[i] + h/8.0 * (k[0][i] + 3*k[2][i]);
+           sData->timeValue = t + h/2.0;
+				   break;
+
+			   case 3:
+			   //yn + 1/2*(k1 - 3*k3 + 4*k4)
+	         for(i = 0; i < nx; ++i)
+             sData->realVars[i] = x0[i] + h/2.0 * (k[0][i] - 3*k[2][i] + 4*k[3][i]);
+           sData->timeValue = t + h;
+				   break;
+
+			   case 4:
+			   //yn + 1/6*(k1 + 4*k3 + k5)
+	         for(i = 0; i < nx; ++i){
+             sData->realVars[i] = x0[i] + h/6.0 * (k[0][i] + 4*k[3][i] + k[4][i]);
+					 }
+           sData->timeValue = t + h;
+			  	 break;
+
+		  }
+	    //f(yn + ...)
+	    /* read input vars */
+      externalInputUpdate(data);
+      data->callback->input_function(data, threadData);
+      /* eval ode equations */
+      data->callback->functionODE(data, threadData);
+	  }
+		t += h;
+    sData->timeValue = t;
+	  solverInfo->currentTime = t;
+
+	  /* save stats */
+    /* steps */
+    solverInfo->solverStatsTmp[0] += 1;
+    /* function ODE evaluation is done directly after this */
+    solverInfo->solverStatsTmp[1] += 4;
+
+
+	  //stepsize
+	  for(i = 0, sum = 0.0; i < nx; ++i)
+	    sum = fmax(fabs(k[0][i] + 4*k[3][i] -(4.5*k[2][i] + k[4][i]))/(fabs(k[4][i]) + fabs(k[2][i]) + fabs(k[3][i])+ + fabs(k[0][i]) +  Atol), sum);
+	  sum *= 2.0/30.0;
+
+
+    h = fmin(0.9*fmax(pow(sum,1/4.0)/(Atol ), 1e-12)*h + 1e-12, (targetTime - h));
+		if(isMaxStepSizeSet && h > maxStepSize) h = maxStepSize;
+	  if (h > 0) rk->h = h;
+
+	  if(t  < targetTime){
+			 memcpy(x0, sData->realVars, nx*sizeof(double));
+			 memcpy(k[0], k[4], nx*sizeof(double));
+			 sim_result.emit(&sim_result, data, threadData);
+		}
+	}
+
+  //assert(sData->timeValue == targetTime);
+	//assert(solverInfo->currentTime == targetTime);
+
+  return 0;
+}
+
+
+
+
 /***************************************    SYM_EULER_IMP     *********************************/
 static int sym_euler_im_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo){
   int retVal,i,j;

SimulationRuntime/c/util/simulation_options.c

@@ -482,6 +482,7 @@ const char *SOLVER_METHOD_NAME[S_MAX] = {
   "symEulerSsc",
   "heun",
   "ida",
+  "rungekutta_ssc",
   "qss"
 };
 
@@ -501,6 +502,7 @@ const char *SOLVER_METHOD_DESC[S_MAX] = {
   "symEulerSsc - symbolic implicit euler with step-size control, [compiler flag +symEuler needed]",
   "heun - Heun's method (Runge-Kutta fixed step, order 2)",
   "ida - Sundials ida solver",
+  "rungekutta_ssc - Runge-Kutta (with step size control, see. Novikov (2016), Solving Stiff Systems of ODEs...)",
   "qss - A QSS solver [experimental]"
 };
 

SimulationRuntime/c/util/simulation_options.h

@@ -149,6 +149,7 @@ enum SOLVER_METHOD
   S_SYM_IMP_EULER, /* 12 */
   S_HEUN,          /* 13 */
   S_IDA,           /* 14 */
+  S_ERKSSC,        /* 15 */
   S_QSS,
 
   S_MAX