@@ -77,7 +77,7 @@ typedef struct RK4_DATA
7777 int work_states_ndims ;
7878 double * b ;
7979 double * c ;
80- double h ;
80+ double h ;
8181}RK4_DATA ;
8282
8383
@@ -331,7 +331,7 @@ int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solv
331331 case S_IDA :
332332 {
333333 IDA_SOLVER * idaData = NULL ;
334- /* Allocate Lobatto6 IIIA work arrays */
334+ /* Allocate Lobatto6 IIIA work arrays */
335335 infoStreamPrint (LOG_SOLVER , 0 , "Initializing IDA DAE Solver" );
336336 idaData = (IDA_SOLVER * ) malloc (sizeof (IDA_SOLVER ));
337337 retValue = ida_solver_initial (data , threadData , solverInfo , idaData );
@@ -370,7 +370,7 @@ int freeSolverData(DATA* data, SOLVER_INFO* solverInfo)
370370 freeSymEulerImp (solverInfo );
371371 }
372372 else if (solverInfo -> solverMethod == S_RUNGEKUTTA || solverInfo -> solverMethod == S_HEUN
373- || solverInfo -> solverMethod == S_ERKSSC )
373+ || solverInfo -> solverMethod == S_ERKSSC )
374374 {
375375 /* free RK work arrays */
376376 for (i = 0 ; i < ((RK4_DATA * )(solverInfo -> solverData ))-> work_states_ndims + 1 ; i ++ )
@@ -819,123 +819,129 @@ static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo)
819819static int rungekutta_step_ssc (DATA * data , threadData_t * threadData , SOLVER_INFO * solverInfo )
820820{
821821 // see.
822- // Solving Stiff Systems of ODEs by Explicit Methods with Conformed Stability
823- // Domains
824- // 2016 9th EUROSIM Congress on Modelling and Simulation
825- // A. E. Novikov...
822+ // Solving Stiff Systems of ODEs by Explicit Methods with Conformed Stability
823+ // Domains
824+ // 2016 9th EUROSIM Congress on Modelling and Simulation
825+ // A. E. Novikov...
826826
827827 RK4_DATA * rk = ((RK4_DATA * )(solverInfo -> solverData ));
828- const double Atol = data -> simulationInfo -> tolerance ;
828+ const double Atol = data -> simulationInfo -> tolerance ;
829829 const double Rtol = data -> simulationInfo -> tolerance ;
830830 const int nx = data -> modelData -> nStates ;
831831 modelica_real h = rk -> h ;
832832 double * * k = rk -> work_states ;
833833 int j , i ;
834- double sum ;
834+ double sum ;
835835 SIMULATION_DATA * sData = (SIMULATION_DATA * )data -> localData [0 ];
836836 SIMULATION_DATA * sDataOld = (SIMULATION_DATA * )data -> localData [1 ];
837837 modelica_real * stateDer = sData -> realVars + nx ;
838838 modelica_real * stateDerOld = sDataOld -> realVars + nx ;
839- double t = sDataOld -> timeValue ;
839+ double t = sDataOld -> timeValue ;
840840 const double targetTime = t + solverInfo -> currentStepSize ;
841- const short isMaxStepSizeSet = (short ) omc_flagValue [FLAG_MAX_STEP_SIZE ];
842- const double maxStepSize = isMaxStepSizeSet ? atof (omc_flagValue [FLAG_MAX_STEP_SIZE ]) : -1 ;
843- double x0 [nx ];
841+ const short isMaxStepSizeSet = (short ) omc_flagValue [FLAG_MAX_STEP_SIZE ];
842+ const double maxStepSize = isMaxStepSizeSet ? atof (omc_flagValue [FLAG_MAX_STEP_SIZE ]) : -1 ;
843+ #if defined(_MSC_VER )
844+ /* handle stupid compilers */
845+ double * x0 = (double * )malloc (nx * sizeof (double ));
846+ #else
847+ double x0 [nx ];
848+ #endif
844849
845850 if (t + h > targetTime )
846- h = solverInfo -> currentStepSize ;
851+ h = solverInfo -> currentStepSize ;
847852
848853 memcpy (k [0 ], stateDerOld , nx * sizeof (double ));
849854 memcpy (x0 , sDataOld -> realVars , nx * sizeof (double ));
850- while (t < targetTime && h > 0 ){
851- //printf("\nfrom %g to %g by %g\n", t, targetTime, h);
852- for (j = 0 ; j < 5 ; ++ j ){
855+ while (t < targetTime && h > 0 ){
856+ //printf("\nfrom %g to %g by %g\n", t, targetTime, h);
857+ for (j = 0 ; j < 5 ; ++ j ){
853858
854- if (j > 0 )
859+ if (j > 0 )
855860 memcpy (k [j ], stateDer , nx * sizeof (double ));
856861
857862 switch (j ){
858863 case 0 :
859- //yn + k1/3
860- for (i = 0 ; i < nx ; ++ i )
864+ //yn + k1/3
865+ for (i = 0 ; i < nx ; ++ i )
861866 sData -> realVars [i ] = x0 [i ] + h * k [0 ][i ]/3.0 ;
862867 sData -> timeValue = t + h /3.0 ;
863- break ;
868+ break ;
864869
865- case 1 :
866- //yn + 1/6(k1 + k2)
867- for (i = 0 ; i < nx ; ++ i )
870+ case 1 :
871+ //yn + 1/6(k1 + k2)
872+ for (i = 0 ; i < nx ; ++ i )
868873 sData -> realVars [i ] = x0 [i ] + h /6.0 * (k [0 ][i ] + k [1 ][i ]);
869874 sData -> timeValue = t + h /3.0 ;
870- break ;
875+ break ;
871876
872- case 2 :
873- //yn + 1/8*(k1 + 3*k3)
874- for (i = 0 ; i < nx ; ++ i )
877+ case 2 :
878+ //yn + 1/8*(k1 + 3*k3)
879+ for (i = 0 ; i < nx ; ++ i )
875880 sData -> realVars [i ] = x0 [i ] + h /8.0 * (k [0 ][i ] + 3 * k [2 ][i ]);
876881 sData -> timeValue = t + h /2.0 ;
877- break ;
882+ break ;
878883
879- case 3 :
880- //yn + 1/2*(k1 - 3*k3 + 4*k4)
881- for (i = 0 ; i < nx ; ++ i )
884+ case 3 :
885+ //yn + 1/2*(k1 - 3*k3 + 4*k4)
886+ for (i = 0 ; i < nx ; ++ i )
882887 sData -> realVars [i ] = x0 [i ] + h /2.0 * (k [0 ][i ] - 3 * k [2 ][i ] + 4 * k [3 ][i ]);
883888 sData -> timeValue = t + h ;
884- break ;
889+ break ;
885890
886- case 4 :
887- //yn + 1/6*(k1 + 4*k3 + k5)
888- for (i = 0 ; i < nx ; ++ i ){
891+ case 4 :
892+ //yn + 1/6*(k1 + 4*k3 + k5)
893+ for (i = 0 ; i < nx ; ++ i ){
889894 sData -> realVars [i ] = x0 [i ] + h /6.0 * (k [0 ][i ] + 4 * k [3 ][i ] + k [4 ][i ]);
890- }
895+ }
891896 sData -> timeValue = t + h ;
892- break ;
897+ break ;
893898
894- }
895- //f(yn + ...)
896- /* read input vars */
899+ }
900+ //f(yn + ...)
901+ /* read input vars */
897902 externalInputUpdate (data );
898903 data -> callback -> input_function (data , threadData );
899904 /* eval ode equations */
900905 data -> callback -> functionODE (data , threadData );
901- }
902- t += h ;
906+ }
907+ t += h ;
903908 sData -> timeValue = t ;
904- solverInfo -> currentTime = t ;
909+ solverInfo -> currentTime = t ;
905910
906- /* save stats */
911+ /* save stats */
907912 /* steps */
908913 solverInfo -> solverStatsTmp [0 ] += 1 ;
909914 /* function ODE evaluation is done directly after this */
910915 solverInfo -> solverStatsTmp [1 ] += 4 ;
911916
912917
913- //stepsize
914- for (i = 0 , sum = 0.0 ; i < nx ; ++ i )
915- 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 );
916- sum *= 2.0 /30.0 ;
918+ //stepsize
919+ for (i = 0 , sum = 0.0 ; i < nx ; ++ i )
920+ 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 );
921+ sum *= 2.0 /30.0 ;
917922
918923
919924 h = fmin (0.9 * fmax (pow (sum ,1 /4.0 )/(Atol ), 1e-12 )* h + 1e-12 , (targetTime - h ));
920- if (isMaxStepSizeSet && h > maxStepSize ) h = maxStepSize ;
921- if (h > 0 ) rk -> h = h ;
925+ if (isMaxStepSizeSet && h > maxStepSize ) h = maxStepSize ;
926+ if (h > 0 ) rk -> h = h ;
922927
923- if (t < targetTime ){
924- memcpy (x0 , sData -> realVars , nx * sizeof (double ));
925- memcpy (k [0 ], k [4 ], nx * sizeof (double ));
926- sim_result .emit (& sim_result , data , threadData );
927- }
928- }
928+ if (t < targetTime ){
929+ memcpy (x0 , sData -> realVars , nx * sizeof (double ));
930+ memcpy (k [0 ], k [4 ], nx * sizeof (double ));
931+ sim_result .emit (& sim_result , data , threadData );
932+ }
933+ }
929934
930935 //assert(sData->timeValue == targetTime);
931- //assert(solverInfo->currentTime == targetTime);
936+ //assert(solverInfo->currentTime == targetTime);
937+ #if defined(_MSC_VER )
938+ /* handle stupid compilers */
939+ free (x0 );
940+ #endif
932941
933942 return 0 ;
934943}
935944
936-
937-
938-
939945/*************************************** SYM_EULER_IMP *********************************/
940946static int sym_euler_im_step (DATA * data , threadData_t * threadData , SOLVER_INFO * solverInfo ){
941947 int retVal ,i ,j ;
@@ -981,7 +987,7 @@ static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* so
981987
982988 /* We calculate k[0] before returning from this function.
983989 * We only want to calculate f() 4 times per call */
984- memcpy (k [0 ], stateDerOld , data -> modelData -> nStates * sizeof (modelica_real ));
990+ memcpy (k [0 ], stateDerOld , data -> modelData -> nStates * sizeof (modelica_real ));
985991
986992 for (j = 1 ; j < rk -> work_states_ndims ; j ++ )
987993 {
@@ -995,7 +1001,7 @@ static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* so
9951001 data -> callback -> input_function (data , threadData );
9961002 /* eval ode equations */
9971003 data -> callback -> functionODE (data , threadData );
998- memcpy (k [j ], stateDer , data -> modelData -> nStates * sizeof (modelica_real ));
1004+ memcpy (k [j ], stateDer , data -> modelData -> nStates * sizeof (modelica_real ));
9991005
10001006 }
10011007
0 commit comments