Implement Kinsol as non-linear solver for non-linear loops over components (#815)

This reverts commit 839c272.

Author: lochel

src/OMSimulatorLib/AlgLoop.cpp

@@ -0,0 +1,111 @@
+/*
+ * This file is part of OpenModelica.
+ *
+ * Copyright (c) 1998-CurrentYear, Open Source Modelica Consortium (OSMC),
+ * c/o Linköpings universitet, Department of Computer and Information Science,
+ * SE-58183 Linköping, Sweden.
+ *
+ * All rights reserved.
+ *
+ * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF GPL VERSION 3 LICENSE OR
+ * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
+ * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
+ * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
+ * ACCORDING TO RECIPIENTS CHOICE.
+ *
+ * The OpenModelica software and the Open Source Modelica
+ * Consortium (OSMC) Public License (OSMC-PL) are obtained
+ * from OSMC, either from the above address,
+ * from the URLs: http://www.ida.liu.se/projects/OpenModelica or
+ * http://www.openmodelica.org, and in the OpenModelica distribution.
+ * GNU version 3 is obtained from: http://www.gnu.org/copyleft/gpl.html.
+ *
+ * This program is distributed WITHOUT ANY WARRANTY; without
+ * even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
+ * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL.
+ *
+ * See the full OSMC Public License conditions for more details.
+ *
+ */
+
+#ifndef _OMS_ALGLOOP_H_
+#define _OMS_ALGLOOP_H_
+
+#include <string>
+#include <vector>
+#include "Types.h"
+#include "DirectedGraph.h"
+
+#include <kinsol/kinsol.h>
+#include <nvector/nvector_serial.h>
+#include <sunlinsol/sunlinsol_dense.h>       /* Default dense linear solver */
+
+namespace oms
+{
+  class System;
+  class DirectedGraph;
+
+  typedef struct KINSOL_USER_DATA {
+    System*         syst;
+    DirectedGraph*  graph;
+    const int       algLoopNumber;
+  }KINSOL_USER_DATA;
+
+  class KinsolSolver
+  {
+    public:
+    ~KinsolSolver();
+    static KinsolSolver* NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance);
+    oms_status_enu_t kinsolSolve(System& syst, DirectedGraph& graph);
+
+    private:
+    /* tolerances */
+    double fnormtol;        /* function tolerance */
+
+    /* work arrays */
+    N_Vector initialGuess;
+    N_Vector uScale;        /* Scaling vector for u */
+    N_Vector fScale;        /* Scaling vector for f(u) */
+    N_Vector fTmp;          /* Vector used for tmp computations */
+
+    /* kinsol internal data */
+    void* kinsolMemory;
+    void* userData;
+    int size;
+
+    /* linear solver data */
+    SUNLinearSolver linSol; /* Linear solver object used by KINSOL */
+    N_Vector y;             /* Template for cloning vectors needed inside linear solver */
+    SUNMatrix J;            /* (Non-)Sparse matrix template for cloning matrices needed within linear solver */
+
+    /* member function */
+    static int nlsKinsolResiduals(N_Vector uu, N_Vector fval, void *userData);
+    static void sundialsErrorHandlerFunction(int error_code, const char *module, const char *function, char *msg, void *userData);
+    static void sundialsInfoHandlerFunction(const char *module, const char *function, char *msg, void *userData);
+  };
+
+  class AlgLoop
+  {
+    public:
+    AlgLoop(oms_alg_solver_enu_t method, double absTol, oms_ssc_t SCC, const int systNumber);
+
+    oms_ssc_t getSCC() {return SCC;}
+    oms_status_enu_t solveAlgLoop(System& syst, DirectedGraph& graph);
+    std::string getAlgSolverName();
+    std::string dumpLoopVars(DirectedGraph& graph);
+
+    private:
+    oms_alg_solver_enu_t algSolverMethod;
+    oms_status_enu_t fixPointIteration(System& syst, DirectedGraph& graph);
+
+    KinsolSolver* kinsolData;
+
+    /* Loop data */
+    const oms_ssc_t SCC;            ///< Strong connected components
+    const int systNumber;
+    double absoluteTolerance;
+  };
+}
+
+#endif // _OMS_ALGLOOP_H_

src/OMSimulatorLib/AlgLoop.h

@@ -28,6 +28,7 @@ ELSE ()
   set(TLM_STRING "-notlm")
 ENDIF ()
 
+list(APPEND OMSIMULATORLIB_SOURCES AlgLoop.cpp)
 list(APPEND OMSIMULATORLIB_SOURCES BusConnector.cpp)
 list(APPEND OMSIMULATORLIB_SOURCES Clock.cpp)
 list(APPEND OMSIMULATORLIB_SOURCES Clocks.cpp)

src/OMSimulatorLib/CMakeLists.txt

@@ -145,7 +145,7 @@ void oms::DirectedGraph::includeGraph(const oms::DirectedGraph& graph, const oms
     addEdge(graph.nodes[graph.edges[i].first].addPrefix(prefix), graph.nodes[graph.edges[i].second].addPrefix(prefix));
 }
 
-int oms::DirectedGraph::getEdgeIndex(const std::vector< std::pair<int, int> >& edges, int from, int to)
+int oms::DirectedGraph::getEdgeIndex(const oms_ssc_t& edges, int from, int to)
 {
   for (int i = 0; i < edges.size(); ++i)
     if (edges[i].first == from && edges[i].second == to)
@@ -243,7 +243,7 @@ std::deque< std::vector<int> > oms::DirectedGraph::getSCCs()
   return components;
 }
 
-const std::vector< std::vector< std::pair<int, int> > >& oms::DirectedGraph::getSortedConnections()
+const std::vector< oms_ssc_t >& oms::DirectedGraph::getSortedConnections()
 {
   if (!sortedConnectionsAreValid)
     calculateSortedConnections();
@@ -253,7 +253,7 @@ const std::vector< std::vector< std::pair<int, int> > >& oms::DirectedGraph::get
 void oms::DirectedGraph::calculateSortedConnections()
 {
   std::deque< std::vector<int> > components = getSCCs();
-  std::vector< std::pair<int, int> > SCC;
+  oms_ssc_t SCC;
   sortedConnections.clear();
 
   for (int i = 0; i < components.size(); ++i)

src/OMSimulatorLib/DirectedGraph.cpp

@@ -43,6 +43,13 @@
 #include <string>
 #include <vector>
 
+/**
+ * @brief Strong connected components data type.
+ *
+ * A vector of pairs of connected components.
+ */
+typedef std::vector< std::pair<int, int> > oms_ssc_t;
+
 namespace oms
 {
   class DirectedGraph
@@ -60,24 +67,24 @@ namespace oms
 
     void includeGraph(const DirectedGraph& graph, const ComRef& prefix);
 
-    const std::vector< std::vector< std::pair<int, int> > >& getSortedConnections();
+    const std::vector< oms_ssc_t >& getSortedConnections();
 
     const std::vector<Connector>& getNodes() const {return nodes;}
-    const std::vector< std::pair<int, int> >& getEdges() const {return edges;}
+    const oms_ssc_t& getEdges() const {return edges;}
 
   private:
     std::deque< std::vector<int> > getSCCs();
     void calculateSortedConnections();
     void strongconnect(int v, std::vector< std::vector<int> > G, int& index, int *d, int *low, std::stack<int>& S, bool *stacked, std::deque< std::vector<int> >& components);
 
-    static int getEdgeIndex(const std::vector< std::pair<int, int> >& edges, int from, int to);
+    static int getEdgeIndex(const oms_ssc_t& edges, int from, int to);
 
   private:
     std::vector<Connector> nodes;
-    std::vector< std::pair<int, int> > edges;
+    oms_ssc_t edges;
 
     std::vector< std::vector<int> > G;
-    std::vector< std::vector< std::pair<int, int> > > sortedConnections;
+    std::vector< oms_ssc_t > sortedConnections;
     bool sortedConnectionsAreValid;
   };
 }

src/OMSimulatorLib/DirectedGraph.h

@@ -62,6 +62,7 @@ oms::Flags::~Flags()
 void oms::Flags::setDefaults()
 {
   addParametersToCSV = false;
+  algLoopSolver = oms_alg_solver_fixedpoint;
   defaultModeIsCS = false;
   deleteTempFiles = true;
   emitEvents = true;
@@ -178,6 +179,17 @@ oms_status_enu_t oms::Flags::AddParametersToCSV(const std::string& value)
   return oms_status_ok;
 }
 
+oms_status_enu_t oms::Flags::AlgLoopSolver(const std::string& value)
+{
+  if (value == "fixedpoint")
+    GetInstance().algLoopSolver = oms_alg_solver_fixedpoint;
+  else if (value == "kinsol")
+    GetInstance().algLoopSolver = oms_alg_solver_kinsol;
+  else
+    return logError("Invalid solver method");
+  return oms_status_ok;
+  }
+
 oms_status_enu_t oms::Flags::ClearAllOptions(const std::string& value)
 {
   GetInstance().setDefaults();

src/OMSimulatorLib/Flags.cpp

@@ -55,6 +55,7 @@ namespace oms
   public:
     static oms_status_enu_t SetCommandLineOption(const std::string& cmd);
 
+    static oms_alg_solver_enu_t AlgLoopSolver() {return GetInstance().algLoopSolver;}
     static bool AddParametersToCSV() {return GetInstance().addParametersToCSV;}
     static bool DefaultModeIsCS() {return GetInstance().defaultModeIsCS;}
     static bool DeleteTempFiles() {return GetInstance().deleteTempFiles;}
@@ -83,6 +84,7 @@ namespace oms
 
   private:
     bool addParametersToCSV;
+    oms_alg_solver_enu_t algLoopSolver;
     bool defaultModeIsCS;
     bool deleteTempFiles;
     bool emitEvents;
@@ -133,6 +135,7 @@ namespace oms
     const std::vector<Flag> flags = {
       {"", "", "FMU or SSP file", re_filename, Flags::Filename, false},
       {"--addParametersToCSV", "", "Export parameters to .csv file (true, [false])", re_default, Flags::AddParametersToCSV, false},
+      {"--algLoopSolver", "", "Specifies the alg. loop solver method ([fixedpoint], kinsol) used for algebraic loops spanning over multiple components.", re_default, Flags::AlgLoopSolver, false},
       {"--clearAllOptions", "", "Reset all flags to default values", re_void, Flags::ClearAllOptions, false},
       {"--deleteTempFiles", "", "Deletes temp files as soon as they are no longer needed ([true], false)", re_bool, Flags::DeleteTempFiles, false},
       {"--emitEvents", "", "Specifies whether events should be emitted or not ([true], false)", re_bool, Flags::EmitEvents, false},
@@ -168,6 +171,7 @@ namespace oms
     };
 
     static oms_status_enu_t AddParametersToCSV(const std::string& value);
+    static oms_status_enu_t AlgLoopSolver(const std::string& value);
     static oms_status_enu_t ClearAllOptions(const std::string& value);
     static oms_status_enu_t DeleteTempFiles(const std::string& value);
     static oms_status_enu_t EmitEvents(const std::string& value);

src/OMSimulatorLib/Flags.h

@@ -31,6 +31,7 @@
 
 #include "System.h"
 
+#include "AlgLoop.h"
 #include "Component.h"
 #include "ComponentFMUCS.h"
 #include "ComponentFMUME.h"
@@ -2365,6 +2366,29 @@ oms_status_enu_t oms::System::importBusConnectorGeometry(const pugi::xml_node& n
   return oms_status_ok;
 }
 
+oms::AlgLoop* oms::System::getAlgLoop(const int systemNumber)
+{
+  if (systemNumber > algLoops.size()-1 || systemNumber < 0)
+  {
+    logError("Invalid system number for algebraic loop.");
+    return NULL;
+  }
+
+  return &algLoops[systemNumber];
+}
+
+oms_status_enu_t oms::System::addAlgLoop(oms_ssc_t SCC, const int algLoopNum)
+{
+  if (loopsNeedUpdate)
+  {
+    algLoops.clear();
+    loopsNeedUpdate = false;
+  }
+  algLoops.push_back( AlgLoop( Flags::AlgLoopSolver(), absoluteTolerance, SCC, algLoopNum ));
+
+  return oms_status_ok;
+}
+
 oms_status_enu_t oms::System::importStartValuesFromSSV()
 {
   for (const auto& file : startValuesFileSources)
@@ -2489,6 +2513,26 @@ oms_status_enu_t oms::System::importStartValuesFromSSVHelper(std::string fileNam
   return oms_status_ok;
 }
 
+oms_status_enu_t oms::System::updateAlgebraicLoops(const std::vector< oms_ssc_t >& sortedConnections)
+{
+  // Instantiate loops
+  if (loopsNeedUpdate)
+  {
+    int systCount = 0;
+    for(int i=0; i<sortedConnections.size(); i++)
+    {
+      if (sortedConnections[i].size() > 1)
+      {
+        addAlgLoop(sortedConnections[i], systCount);
+        systCount++;
+      }
+    }
+    loopsNeedUpdate = false;
+  }
+
+  return oms_status_ok;
+}
+
 oms_status_enu_t oms::System::importParameterMappingFromSSM(std::string fileName, std::multimap<ComRef, ComRef> &mappedEntry)
 {
   std::string tempdir = getModel()->getTempDirectory();
@@ -2528,3 +2572,8 @@ oms_status_enu_t oms::System::importParameterMappingFromSSM(std::string fileName
 
   return oms_status_ok;
 }
+
+oms_status_enu_t oms::System::solveAlgLoop(DirectedGraph& graph, int loopNumber)
+{
+  return algLoops[loopNumber].solveAlgLoop(*this, graph);
+}

src/OMSimulatorLib/System.cpp

@@ -32,6 +32,7 @@
 #ifndef _OMS_SYSTEM_H_
 #define _OMS_SYSTEM_H_
 
+#include "AlgLoop.h"
 #include "BusConnector.h"
 #include "Clock.h"
 #include "ComRef.h"
@@ -58,6 +59,7 @@
 
 namespace oms
 {
+  class AlgLoop;
   class Component;
   class Model;
   class Variable;
@@ -132,6 +134,8 @@ namespace oms
     virtual oms_status_enu_t reset() = 0;
     virtual oms_status_enu_t stepUntil(double stopTime, void (*cb)(const char* ident, double time, oms_status_enu_t status)) = 0;
 
+    double getTime() const {return time;}
+
     oms_status_enu_t getBoolean(const ComRef& cref, bool& value);
     oms_status_enu_t getInteger(const ComRef& cref, int& value);
     oms_status_enu_t getReal(const ComRef& cref, double& value);
@@ -160,12 +164,18 @@ namespace oms
     double getMaximumStepSize() {return maximumStepSize;}
     oms_solver_enu_t getSolver() {return solverMethod;}
 
+    AlgLoop* getAlgLoop(const int systemNumber);
+    oms_status_enu_t addAlgLoop(oms_ssc_t SCC, const int algLoopNum);
+    oms_status_enu_t updateAlgebraicLoops(const std::vector< oms_ssc_t >& sortedConnections);
+    oms_status_enu_t solveAlgLoop(DirectedGraph& graph, int loopNumber);
+
     bool useThreadPool();
     ctpl::thread_pool& getThreadPool();
 
     std::string getUniqueID() const;
     std::map<std::string, std::string> startValuesFileSources;  ///< ssvFileSource mapped with ssmFilesource if mapping is provided, otherwise only ssvFilesource entry is made
   protected:
+    double time;
     System(const ComRef& cref, oms_system_enu_t type, Model* parentModel, System* parentSystem, oms_solver_enu_t solverMethod);
 
     // stop the compiler generating methods copying the object
@@ -190,6 +200,9 @@ namespace oms
     std::unordered_map<unsigned int /*result file var ID*/, unsigned int /*allVariables ID*/> resultFileMapping;
     std::unordered_map<ComRef, bool> exportConnectors;
 
+  protected:
+    bool loopsNeedUpdate = true;
+
   private:
     ComRef cref;
     oms_system_enu_t type;
@@ -217,6 +230,8 @@ namespace oms
     oms_status_enu_t importStartValuesFromSSV();
     oms_status_enu_t importStartValuesFromSSVHelper(std::string fileName, std::multimap<ComRef, ComRef> &mappedEntry);
     oms_status_enu_t importParameterMappingFromSSM(std::string fileName, std::multimap<ComRef, ComRef> &mappedEntry);
+
+    std::vector<AlgLoop> algLoops;    ///< Vector of algebraic loop objects
   };
 }