Merge branch 'master' of https://github.com/OpenModelica/OMCompiler

Compiler/BackEnd/ExpressionSolve.mo

@@ -668,8 +668,7 @@ algorithm
       orhs := DAE.IFEXP(Expression.makeNoEvent(DAE.RELATION(
                 dt,
                   DAE.EQUAL(tp),
-                DAE.RCONST(0.0),-1,NONE())),
-              Expression.makePureBuiltinCall("$_old", {X}, tp), rhs);
+                DAE.RCONST(0.0),-1,NONE())), X, rhs);
     end if;
   end if;
 

SimulationRuntime/c/Makefile.objs

@@ -33,11 +33,11 @@ MATH_HFILES = blaswrap.h
 
 SOLVER_OBJS_MINIMAL=delay$(OBJ_EXT) events$(OBJ_EXT) external_input$(OBJ_EXT) linearSystem$(OBJ_EXT) linearSolverLapack$(OBJ_EXT) linearSolverTotalPivot$(OBJ_EXT) mixedSystem$(OBJ_EXT) mixedSearchSolver$(OBJ_EXT) model_help$(OBJ_EXT) nonlinearSystem$(OBJ_EXT) nonlinearSolverHomotopy$(OBJ_EXT) omc_math$(OBJ_EXT) solver_main$(OBJ_EXT) stateset$(OBJ_EXT)
 ifeq ($(OMC_MINIMAL_RUNTIME),)
-SOLVER_OBJS=$(SOLVER_OBJS_MINIMAL) kinsolSolver$(OBJ_EXT) linearSolverLis$(OBJ_EXT) linearSolverUmfpack$(OBJ_EXT) dassl$(OBJ_EXT) radau$(OBJ_EXT) nonlinearSolverHybrd$(OBJ_EXT) nonlinearSolverNewton$(OBJ_EXT) newtonIteration$(OBJ_EXT)
+SOLVER_OBJS=$(SOLVER_OBJS_MINIMAL) kinsolSolver$(OBJ_EXT) linearSolverLis$(OBJ_EXT) linearSolverUmfpack$(OBJ_EXT) dassl$(OBJ_EXT) radau$(OBJ_EXT) sym_imp_euler$(OBJ_EXT) nonlinearSolverHybrd$(OBJ_EXT) nonlinearSolverNewton$(OBJ_EXT) newtonIteration$(OBJ_EXT)
 else
 SOLVER_OBJS=$(SOLVER_OBJS_MINIMAL)
 endif
-SOLVER_HFILES = dassl.h delay.h epsilon.h events.h external_input.h linearSystem.h mixedSystem.h model_help.h nonlinearSystem.h radau.h solver_main.h stateset.h
+SOLVER_HFILES = dassl.h delay.h epsilon.h events.h external_input.h linearSystem.h mixedSystem.h model_help.h nonlinearSystem.h radau.h sym_imp_euler.h solver_main.h stateset.h
 
 INITIALIZATION_OBJS = initialization$(OBJ_EXT)
 INITIALIZATION_HFILES = initialization.h

SimulationRuntime/c/simulation/solver/CMakeLists.txt

@@ -14,14 +14,14 @@ dassl.c           kinsolSolver.c            linearSystem.c             nonlinear
 delay.c           linearSolverLapack.c      mixedSearchSolver.c        nonlinearSolverNewton.c  newtonIteration.c solver_main.c
 linearSolverLis.c mixedSystem.c             nonlinearSystem.c          stateset.c
 events.c          linearSolverTotalPivot.c  model_help.c               omc_math.c
-external_input.c  linearSolverUmfpack.c     nonlinearSolverHomotopy.c)
+external_input.c  linearSolverUmfpack.c     nonlinearSolverHomotopy.c  sym_imp_euler.c)
 
 SET(solver_headers ../../../../3rdParty/Cdaskr/solver/ddaskr_types.h
 dassl.h    external_input.h          linearSolverUmfpack.h  nonlinearSolverHomotopy.h  radau.h
 delay.h    kinsolSolver.h            linearSystem.h         nonlinearSolverHybrd.h     solver_main.h
 linearSolverLapack.h      mixedSearchSolver.h    nonlinearSolverNewton.h newtonIteration.h   stateset.h
 epsilon.h  linearSolverLis.h         mixedSystem.h          nonlinearSystem.h
-events.h   linearSolverTotalPivot.h  model_help.h           omc_math.h)
+events.h   linearSolverTotalPivot.h  model_help.h           omc_math.h	       sym_imp_euler.h)
 
 # Library util
 ADD_LIBRARY(solver ${solver_sources} ${solver_headers})

SimulationRuntime/c/simulation/solver/perform_simulation.c

@@ -194,6 +194,8 @@ int prefixedName_performSimulation(DATA* data, SOLVER_INFO* solverInfo)
 
       updateContinuousSystem(data);
 
+      if (solverInfo->solverMethod == S_SYM_IMP_EULER) data->callback->symEulerUpdate(data, solverInfo->solverStepSize);
+
       saveZeroCrossings(data);
       messageClose(LOG_SOLVER);
 

SimulationRuntime/c/simulation/solver/solver_main.c

@@ -38,6 +38,7 @@
 #include "openmodelica_func.h"
 #include "initialization/initialization.h"
 #include "nonlinearSystem.h"
+#include "newtonIteration.h"
 #include "dassl.h"
 #include "delay.h"
 #include "events.h"
@@ -49,6 +50,7 @@
 #include "meta/meta_modelica.h"
 #include "simulation/solver/epsilon.h"
 #include "linearSystem.h"
+#include "sym_imp_euler.h"
 
 #include "optimization/OptimizerInterface.h"
 
@@ -136,6 +138,9 @@ int solver_main_step(DATA* data, SOLVER_INFO* solverInfo)
   case S_SYM_EULER:
     retVal = sym_euler_im_step(data, solverInfo);
     return retVal;
+  case S_SYM_IMP_EULER:
+    retVal = sym_euler_im_with_step_size_control_step(data, solverInfo);
+    return retVal;
   }
 
   TRACE_POP
@@ -173,13 +178,24 @@ int initializeSolverData(DATA* data, SOLVER_INFO* solverInfo)
   solverInfo->stateEvents = 0;
   solverInfo->sampleEvents = 0;
 
+  /* if FLAG_NOEQUIDISTANT_GRID is set, choose dassl step method */
+  if (omc_flag[FLAG_NOEQUIDISTANT_GRID])
+  {
+    solverInfo->integratorSteps = 1; /* TRUE */
+  }
+
   /* set tolerance for ZeroCrossings */
   setZCtol(fmin(simInfo->stepSize, simInfo->tolerance));
 
   switch (solverInfo->solverMethod)
   {
   case S_SYM_EULER:
   case S_EULER: break;
+  case S_SYM_IMP_EULER:
+  {
+    allocateSymEulerImp(solverInfo, data->modelData.nStates);
+    break;
+  }
   case S_RUNGEKUTTA:
   {
     /* Allocate RK work arrays */
@@ -290,7 +306,11 @@ int freeSolverData(DATA* data, SOLVER_INFO* solverInfo)
   int i;
 
   /* deintialize solver related workspace */
-  if(solverInfo->solverMethod == S_RUNGEKUTTA)
+  if (solverInfo->solverMethod == S_SYM_IMP_EULER)
+  {
+    freeSymEulerImp(solverInfo);
+  }
+  else if(solverInfo->solverMethod == S_RUNGEKUTTA)
   {
     /* free RK work arrays */
     for(i = 0; i < ((RK4_DATA*)(solverInfo->solverData))->work_states_ndims + 1; i++)
@@ -563,7 +583,7 @@ int finishSimulation(DATA* data, SOLVER_INFO* solverInfo, const char* outputVari
  *  \param [in]  [solverID] selects the ode solver
  *  \param [in]  [outputVariablesAtEnd] ???
  *
- *  This is the main function of the solver it perform the simulation.
+ *  This is the main function of the solver, it performs the simulation.
  */
 int solver_main(DATA* data, const char* init_initMethod, const char* init_file,
     double init_time, int lambda_steps, int solverID, const char* outputVariablesAtEnd)
@@ -662,6 +682,7 @@ int solver_main(DATA* data, const char* init_initMethod, const char* init_file,
       retVal = data->callback->performSimulation(data, &solverInfo);
       omc_alloc_interface.collect_a_little();
       /* terminate the simulation */
+      if (solverInfo.solverMethod == S_SYM_IMP_EULER) data->callback->symEulerUpdate(data, 0);
       finishSimulation(data, &solverInfo, outputVariablesAtEnd);
       omc_alloc_interface.collect_a_little();
     }
@@ -716,7 +737,6 @@ static int sym_euler_im_step(DATA* data, SOLVER_INFO* solverInfo){
   return retVal;
 }
 
-
 /***************************************    RK4      ***********************************/
 static int rungekutta_step(DATA* data, SOLVER_INFO* solverInfo)
 {

SimulationRuntime/c/simulation/solver/solver_main.h

@@ -48,6 +48,7 @@ typedef struct SOLVER_INFO
   double currentStepSize;
   double laststep;
   int solverMethod;
+  double solverStepSize; /* used by implicit radau solver */
 
   /* set by solver if an internal root finding method is activated  */
   modelica_boolean solverRootFinding;
@@ -56,10 +57,16 @@ typedef struct SOLVER_INFO
   LIST* eventLst;
   int didEventStep;
 
+  /* radau_new
+  void* userdata;
+*/
   /* stats */
   unsigned long stateEvents;
   unsigned long sampleEvents;
 
+  /* further options */
+  int integratorSteps;
+
   void* solverData;
 }SOLVER_INFO;
 
@@ -75,8 +82,6 @@ extern int solver_main(DATA* data, const char* init_initMethod,
 extern int initializeSolverData(DATA* data, SOLVER_INFO* solverInfo);
 extern int freeSolverData(DATA* data, SOLVER_INFO* solverInfo);
 
-extern int freeSolverData(DATA* data, SOLVER_INFO* solverInfo);
-
 extern int initializeModel(DATA* data, const char* init_initMethod,
     const char* init_file, double init_time, int lambda_steps);