- fixed the ordering of the array elements in RefArray class

- added a counter for the number of ODE-evaluations to CVode
- included version check for Sundials, colored jacobians are now enabled if the version is new enough
- some typo fixes

Author: Marcus Walther

SimulationRuntime/cpp/CMakeLists.txt

@@ -297,10 +297,29 @@ MESSAGE(STATUS ${LAPACK_LIBRARIES})
 # Find Sundials solvers
 IF(USE_SUNDIALS)
   FIND_PATH(SUNDIALS_INCLUDE_DIR cvode/cvode.h PATHS ${SUNDIALS_INLCUDE_HOME} $ENV{SUNDIALS_ROOT}/include)
-  if (NOT SUNDIALS_INCLUDE_DIR)
-      MESSAGE(FATAL_ERROR "Could not find Sundials, specify environment variable SUNDIALS_ROOT")
+  IF (NOT SUNDIALS_INCLUDE_DIR)
+    MESSAGE(FATAL_ERROR "Could not find Sundials, specify environment variable SUNDIALS_ROOT")
+  ELSE(NOT SUNDIALS_INCLUDE_DIR)
+    FIND_PATH(SUNDIALS_CONFIG_FILE "sundials_config.h" ${SUNDIALS_INCLUDE_DIR} "${SUNDIALS_INCLUDE_DIR}/sundials")
+    SET(SUNDIALS_CONFIG_FILE "${SUNDIALS_CONFIG_FILE}/sundials_config.h")
+    FILE(READ "${SUNDIALS_CONFIG_FILE}" SUNDIALS_CONFIG_FILE_CONTENT)
+    STRING(REGEX MATCH "#define SUNDIALS_PACKAGE_VERSION .([0-9]+)\\.([0-9]+)\\.([0-9]+)." SUNDIALS_CONFIG_FILE_CONTENT ${SUNDIALS_CONFIG_FILE_CONTENT})
+    STRING(REGEX REPLACE "#define SUNDIALS_PACKAGE_VERSION .([0-9]+)\\.([0-9]+)\\.([0-9]+)." "\\1;\\2;\\3" SUNDIALS_CONFIG_FILE_VERSION ${SUNDIALS_CONFIG_FILE_CONTENT})
+    LIST(GET SUNDIALS_CONFIG_FILE_VERSION 0 SUNDIALS_MAJOR_VERSION)
+    IF(SUNDIALS_MAJOR_VERSION)
+      ADD_DEFINITIONS("-DSUNDIALS_MAJOR_VERSION=${SUNDIALS_MAJOR_VERSION}")
+    ELSE()
+      MESSAGE(FATAL_ERROR "Could not determine sundials version")
+    ENDIF()
+    LIST(GET SUNDIALS_CONFIG_FILE_VERSION 1 SUNDIALS_MINOR_VERSION)
+    IF(SUNDIALS_MINOR_VERSION)
+      ADD_DEFINITIONS("-DSUNDIALS_MINOR_VERSION=${SUNDIALS_MINOR_VERSION}")
+    ELSE()
+      MESSAGE(FATAL_ERROR "Could not determine sundials version")
+    ENDIF()
+    LIST(GET SUNDIALS_CONFIG_FILE_VERSION 2 SUNDIALS_PATCH_VERSION)
+    MESSAGE(STATUS "Using sundials ${SUNDIALS_MAJOR_VERSION}.${SUNDIALS_MINOR_VERSION}.${SUNDIALS_PATCH_VERSION}")
   ENDIF()
-  MESSAGE(STATUS ${SUNDIALS_INCLUDE_DIR})
 
   FIND_LIBRARY(SUNDIALS_NVECSERIAL_LIB "sundials_nvecserial" PATHS ${SUNDIALS_LIBRARY_RELEASE_HOME} $ENV{SUNDIALS_ROOT}/lib)
   IF(NOT SUNDIALS_NVECSERIAL_LIB)

SimulationRuntime/cpp/Core/SimController/Main.cpp

@@ -32,8 +32,9 @@ int main(int argc, const char* argv[])
             ("stop-time,e", po::value< double >()->default_value(1.0),  "simulation stop time")
             ("step-size,f", po::value< double >()->default_value(1e-2),  "simulation step size")
             ("solver,i", po::value< string >()->default_value("euler"),  "solver method")
-            ("number-of-intervalls,v", po::value< int >()->default_value(500),  "number of intervalls")
-            ("tollerance,y", po::value< double >()->default_value(1e-6),  "solver tollerance")
+            ("number-of-intervals,v", po::value< int >()->default_value(500),  "number of intervals")
+            ("tolerance,y", po::value< double >()->default_value(1e-6),  "solver tolerance")
+            ("solverLog,l", "print additional solver information after simulation")
            ;
         po::variables_map vm;
         po::store(po::parse_command_line(argc, argv, desc,
@@ -47,8 +48,8 @@ int main(int argc, const char* argv[])
         string runtime_lib_path;
         double starttime =  vm["start-time"].as<double>();
         double stoptime = vm["stop-time"].as<double>();
-        double stepsize =  stoptime/vm["number-of-intervalls"].as<int>();
-        double tollerance =vm["tollerance"].as<double>();
+        double stepsize =  stoptime/vm["number-of-intervals"].as<int>();
+        double tollerance =vm["tolerance"].as<double>();
         string solver =  vm["solver"].as<string>();
         if (vm.count("runtime-library"))
         {

SimulationRuntime/cpp/Include/Core/Math/Array.h

@@ -428,7 +428,7 @@ class RefArrayDim2 : public RefArray<T, size1*size2>
   virtual const T& operator()(const vector<size_t>& idx) const
   {
     return *(RefArray<T, size1*size2>::
-             _ref_array[idx[0]-1 + size1*(idx[1]-1)]);
+             _ref_array[(idx[0]-1)*size2 + (idx[1]-1)]);
   }
 
   /**
@@ -438,7 +438,7 @@ class RefArrayDim2 : public RefArray<T, size1*size2>
   virtual T& operator()(const vector<size_t>& idx)
   {
     return *(RefArray<T, size1*size2>::
-             _ref_array[idx[0]-1 + size1*(idx[1]-1)]);
+             _ref_array[(idx[0]-1)*size2 + (idx[1]-1)]);
   }
 
   /**
@@ -449,7 +449,7 @@ class RefArrayDim2 : public RefArray<T, size1*size2>
   inline virtual T& operator()(size_t i, size_t j)
   {
     return *(RefArray<T, size1*size2>::
-             _ref_array[i-1 + size1*(j-1)]);
+             _ref_array[(i-1)*size2 + (j-1)]);
   }
 
   /**
@@ -553,7 +553,7 @@ class RefArrayDim3 : public RefArray<T, size1*size2*size3>
   virtual const T& operator()(const vector<size_t>& idx) const
   {
     return *(RefArray<T, size1*size2*size3>::
-             _ref_array[idx[0]-1 + size1*(idx[1]-1 + size2*(idx[2]-1))]);
+             _ref_array[size3*(idx[0]-1 + size2*(idx[1]-1)) + idx[2]-1]);
   }
 
   /**
@@ -563,7 +563,7 @@ class RefArrayDim3 : public RefArray<T, size1*size2*size3>
   virtual T& operator()(const vector<size_t>& idx)
   {
     return *(RefArray<T, size1*size2*size3>::
-             _ref_array[idx[0]-1 + size1*(idx[1]-1 + size2*(idx[2]-1))]);
+             _ref_array[size3*(idx[0]-1 + size2*(idx[1]-1)) + idx[2]-1]);
   }
 
   /**
@@ -575,7 +575,7 @@ class RefArrayDim3 : public RefArray<T, size1*size2*size3>
   inline virtual T& operator()(size_t i, size_t j, size_t k)
   {
     return *(RefArray<T, size1*size2*size3>::
-             _ref_array[i-1 + size1*(j-1 + size2*(k-1))]);
+             _ref_array[size3*(i-1 + size2*(j-1)) + (k-1)]);
   }
 
   /**

SimulationRuntime/cpp/Include/Solver/CVode/CVode.h

@@ -234,6 +234,8 @@ double
    IMixedSystem* _mixed_system;
    ITime* _time_system;
 
+   int _numberOfOdeEvaluations;
+
    #ifdef RUNTIME_PROFILING
    std::vector<MeasureTimeData> measureTimeFunctionsArray;
    MeasureTimeValues *measuredFunctionStartValues, *measuredFunctionEndValues, *solveFunctionStartValues, *solveFunctionEndValues;

SimulationRuntime/cpp/Solver/CVode/CVode.cpp

@@ -31,6 +31,7 @@ Cvode::Cvode(IMixedSystem* system, ISolverSettings* settings)
       _event_system(NULL),
       _mixed_system(NULL),
       _time_system(NULL),
+      _numberOfOdeEvaluations(0),
     _delta(NULL),
   _deltaInv(NULL),
     _ysave(NULL),
@@ -275,7 +276,7 @@ void Cvode::initialize()
 
   // Use own jacobian matrix
   // Check if Colored Jacobians are worth to use
-   #if defined(_WIN32) || defined(__WIN32__) || defined(WIN32)
+   #if SUNDIALS_MAJOR_VERSION >= 2 || (SUNDIALS_MAJOR_VERSION == 2 && SUNDIALS_MINOR_VERSION >= 5)
     _maxColors = _system->getAMaxColors();
     if(_maxColors < _dimSys && _continuous_system->getDimContinuousStates() > 0)
     {
@@ -717,6 +718,7 @@ int Cvode::calcFunction(const double& time, const double* y, double* f)
     _continuous_system->setContinuousStates(y);
     _continuous_system->evaluateODE(IContinuous::CONTINUOUS);
     _continuous_system->getRHS(f);
+    _numberOfOdeEvaluations++;
   }      //workaround until exception can be catch from c- libraries
   catch (std::exception & ex )
   {
@@ -968,6 +970,7 @@ void Cvode::writeSimulationInfo()
   BOOST_LOG_SEV(slg, cvode_normal)<< " Linear solver setups " << "nsetups: " << nsetups;
   BOOST_LOG_SEV(slg, cvode_normal)<< " Nonlinear iterations " << "nni: " << nni;
   BOOST_LOG_SEV(slg, cvode_normal)<< " Convergence failures " << "ncfn: " << ncfn;
+  BOOST_LOG_SEV(slg, cvode_normal)<< " Number of evaluateODE calls " << "eODE: " << _numberOfOdeEvaluations;
 
 #endif
   //// Solver