From 7d596a735eccdd073dd22369b86012ff9f07546b Mon Sep 17 00:00:00 2001 From: Lennart Ochel Date: Fri, 22 Jan 2021 13:47:20 +0100 Subject: [PATCH] Revert "Revert "Implement Kinsol as non-linear solver for non-linear loops over components (#815)" (#847)" This reverts commit 839c272aba25544cc4c533153a1115ccce166126. --- src/OMSimulatorLib/AlgLoop.cpp | 541 +++++++++++++++++++++++++++ src/OMSimulatorLib/AlgLoop.h | 111 ++++++ src/OMSimulatorLib/CMakeLists.txt | 1 + src/OMSimulatorLib/DirectedGraph.cpp | 6 +- src/OMSimulatorLib/DirectedGraph.h | 17 +- src/OMSimulatorLib/Flags.cpp | 12 + src/OMSimulatorLib/Flags.h | 4 + src/OMSimulatorLib/System.cpp | 49 +++ src/OMSimulatorLib/System.h | 15 + src/OMSimulatorLib/SystemSC.cpp | 76 +--- src/OMSimulatorLib/SystemSC.h | 5 - src/OMSimulatorLib/SystemWC.cpp | 83 +--- src/OMSimulatorLib/SystemWC.h | 5 - src/OMSimulatorLib/Types.h | 6 + 14 files changed, 785 insertions(+), 146 deletions(-) create mode 100644 src/OMSimulatorLib/AlgLoop.cpp create mode 100644 src/OMSimulatorLib/AlgLoop.h diff --git a/src/OMSimulatorLib/AlgLoop.cpp b/src/OMSimulatorLib/AlgLoop.cpp new file mode 100644 index 000000000..0746bbf0c --- /dev/null +++ b/src/OMSimulatorLib/AlgLoop.cpp @@ -0,0 +1,541 @@ +/* + * 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. + * + */ + +#include "AlgLoop.h" + +#include "Component.h" +#include "DirectedGraph.h" +#include "Flags.h" +#include "System.h" + +/** + * @brief Check flag returned by KINSOL function and log error + * + * @param flag Flag with some error value retured by KINSOL function + * @param functionName Name of KINSOL function that returned the flag + * @return true The flag was no error code + * @return false The flag was an error code + */ +inline bool checkFlag(int flag, std::string functionName) +{ + if (flag < 0) + { + logError("SUNDIALS_ERROR: " + functionName + "() failed with flag = " + std::to_string(flag)); + return false; + } + return true; +} + +/** + * @brief Error handler function given to KINSOL. + * + * @param errorCode Error code from KINSOL + * @param module Name of the module reporting the error. + * @param function Name of the function in which the error occurred. + * @param msg Error Message. + * @param userData Pointer to user data. Unused. + */ +void oms::KinsolSolver::sundialsErrorHandlerFunction(int errorCode, const char *module, + const char *function, char *msg, + void *userData) +{ + KINSOL_USER_DATA* kinsolUserData; + std::string systNum = "unknown"; + std::string mod = module; + std::string func = function; + + if (userData != NULL) { + kinsolUserData = (KINSOL_USER_DATA *)userData; + systNum = std::to_string(kinsolUserData->algLoopNumber); + } + + logError("SUNDIALS_ERROR: [system] " + systNum + " [module] " + mod + " | [function] " + func + +" | [error_code] " + std::to_string(errorCode)); + logError(msg); +} + +/** + * @brief Info handler function given to KINSOL. + * + * Will only print information when debug loging is active. + * + * @param module Name of the module reporting the information. + * @param function Name of the function reporting the information. + * @param msg Message. + * @param userData Pointer to user data. Unused. + */ +void oms::KinsolSolver::sundialsInfoHandlerFunction(const char *module, const char *function, + char *msg, void *userData) +{ + KINSOL_USER_DATA* kinsolUserData; + std::string systNum = "unknown"; + std::string mod = module; + std::string func = function; + + if (userData != NULL) { + kinsolUserData = (KINSOL_USER_DATA *)userData; + systNum = std::to_string(kinsolUserData->algLoopNumber); + } + + logDebug("SUNDIALS_INFO: [system] " + systNum + " [module] " + mod + " | [function] " + func); + logDebug(msg); +} + +/** + * @brief Residual function for KINSOL + * + * @param uu Input value + * @param fval Contains residual at output + * @param userData Pointer to user data to access System, AlgLoop, Directed Graph and CSS as well as getReal and setReal + * @return int Return 0 on success, 1 if an recoverable error occured, -1 if a fatal error occured. + */ +int oms::KinsolSolver::nlsKinsolResiduals(N_Vector uu, N_Vector fval, void *userData) +{ + double *uu_data = NV_DATA_S(uu); + double *fval_data = NV_DATA_S(fval); + + KINSOL_USER_DATA* kinsoluserData = (KINSOL_USER_DATA*) userData; + System* syst = kinsoluserData->syst; + AlgLoop* algLoop = syst->getAlgLoop(kinsoluserData->algLoopNumber); + DirectedGraph* graph = kinsoluserData->graph; + const oms_ssc_t SCC = algLoop->getSCC(); + + const int size = SCC.size(); + double maxRes; + oms_status_enu_t status; + + // Set values from uu + for (int i=0; isetReal(graph->getNodes()[output].getName(), uu_data[i]); + if (status == oms_status_discard || status == oms_status_error || status == oms_status_warning) + { + return 1 /* recoverable error */; + } + else if (status == oms_status_fatal) + { + return -1 /* not recoverable error */; + } + } + + // Get updated values and calulate residual + for (int i=0; igetReal(graph->getNodes()[output].getName(), fval_data[i]); + if (status == oms_status_discard || status == oms_status_error || status == oms_status_warning) + { + return 1 /* recoverable error */; + } + else if (status == oms_status_fatal) + { + return -1 /* not recoverable error */; + } + fval_data[i] = fval_data[i] - uu_data[i]; + } + + return 0 /* success */; +} + +/** + * @brief Destroy the oms::KinsolSolver::KinsolSolver object + * + */ +oms::KinsolSolver::~KinsolSolver() +{ + KINFree(&(this->kinsolMemory)); + N_VDestroy_Serial(this->initialGuess); /* TODO: Or is N_VDestroy_Serial correct? It won't free internal data */ + N_VDestroy_Serial(this->uScale); + N_VDestroy_Serial(this->fScale); + N_VDestroy_Serial(this->fTmp); + + /* Free linear solver data */ + SUNLinSolFree(this->linSol); + SUNMatDestroy(this->J); + N_VDestroy_Serial(this->y); + + delete((KINSOL_USER_DATA*)(this->userData)); +} + +/** + * @brief Create new oms::KinsolSolver::KinsolSolver object + * + * @param algLoopNum Number of algebraic loop + * @param size Dimension of algebraic loop + * @param absoluteTolerance Tolerance used for solving the loop + * @return oms::KinsolSolver* Retruns pointer to KinsolSolver object + */ +oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance) +{ + int flag; + int printLevel; + KinsolSolver* kinsolSolver = new KinsolSolver(); + + logDebug("Create new KinsolSolver object for algebraic loop number " + std::to_string(algLoopNum)); + + kinsolSolver->size = size; + + /* Allocate memory */ + kinsolSolver->initialGuess = N_VNew_Serial(kinsolSolver->size); + kinsolSolver->uScale = N_VNew_Serial(kinsolSolver->size); + kinsolSolver->fScale = N_VNew_Serial(kinsolSolver->size); + kinsolSolver->fTmp = N_VNew_Serial(kinsolSolver->size); + kinsolSolver->y = N_VNew_Serial(kinsolSolver->size); + kinsolSolver->kinsolMemory = NULL; + + /* Create KINSOL memory block */ + kinsolSolver->kinsolMemory = KINCreate(); + if (kinsolSolver->kinsolMemory == NULL) + { + logError("SUNDIALS_ERROR: KINCreate() failed"); + return NULL; + } + + /* Set user data given to KINSOL */ + kinsolSolver->userData = new KINSOL_USER_DATA{/*.syst=*/NULL, /*.graph=*/NULL, /*.algLoopNumber=*/algLoopNum}; + flag = KINSetUserData(kinsolSolver->kinsolMemory, kinsolSolver->userData); + if (!checkFlag(flag, "KINSetUserData")) return NULL; + + /* Set error handler and print level */ + if (Log::DebugEnabled()) { + logDebug("SUNDIALS KINSOL: Set print level to maximum."); + printLevel = 3; + } else { + printLevel = 0; + } + flag = KINSetPrintLevel(kinsolSolver->kinsolMemory, printLevel); + if (!checkFlag(flag, "KINSetPrintLevel")) return NULL; + + flag = KINSetErrHandlerFn(kinsolSolver->kinsolMemory, sundialsErrorHandlerFunction, kinsolSolver->userData); + if (!checkFlag(flag, "KINSetErrHandlerFn")) return NULL; + + flag = KINSetInfoHandlerFn(kinsolSolver->kinsolMemory, sundialsInfoHandlerFunction, kinsolSolver->userData); + if (!checkFlag(flag, "KINSetInfoHandlerFn")) return NULL; + + /* Initialize KINSOL object */ + flag = KINInit(kinsolSolver->kinsolMemory, nlsKinsolResiduals, kinsolSolver->initialGuess); + if (!checkFlag(flag, "KINInit")) return NULL; + + /* Create matrix object */ + kinsolSolver->J = SUNDenseMatrix(kinsolSolver->size, kinsolSolver->size); + + /* Create linear solver object */ + kinsolSolver->linSol = SUNLinSol_Dense(kinsolSolver->y, kinsolSolver->J); + + /* Set linear solver */ + flag = KINSetLinearSolver(kinsolSolver->kinsolMemory, kinsolSolver->linSol, kinsolSolver->J); + if (!checkFlag(flag, "KINSetLinearSolver")) return NULL; + + /* Set Jacobian for linear solver */ + flag = KINSetJacFn(kinsolSolver->kinsolMemory, NULL); /* Use KINSOLs internal difference quotient function for Jacobian approximation */ + if (!checkFlag(flag, "KINSetJacFn")) return NULL; + + /* Set function-norm stopping tolerance */ + kinsolSolver->fnormtol = absoluteTolerance; + flag = KINSetFuncNormTol(kinsolSolver->kinsolMemory, kinsolSolver->fnormtol); + if (!checkFlag(flag, "KINSetFuncNormTol")) return NULL; + + /* Set scaled-step stopping tolerance */ + flag = KINSetScaledStepTol(kinsolSolver->kinsolMemory, 0.0); + if (!checkFlag(flag, "KINSetScaledStepTol")) return NULL; + + /* Set max. number of nonlinear iterations */ + flag = KINSetNumMaxIters(kinsolSolver->kinsolMemory, 100 * kinsolSolver->size); + if (!checkFlag(flag, "KINSetNumMaxIters")) return NULL; + + flag = KINSetNoInitSetup(kinsolSolver->kinsolMemory, SUNFALSE); + if (!checkFlag(flag, "KINSetNoInitSetup")) return NULL; + + /* Set default scaling */ + double *uScaleData = NV_DATA_S(kinsolSolver->uScale); + double *fScaleData = NV_DATA_S(kinsolSolver->fScale); + for (int i=0; i < kinsolSolver->size; i++) + { + uScaleData[i] = 1.0; + fScaleData[i] = 1.0; + } + + return kinsolSolver; +} + +/** + * @brief Solve algebraic system with KINSOL + * + * @param syst Reference to System object + * @param graph Reference to graph opbject + * @return oms_status_enu_t Return `oms_status_ok` on success, + * `oms_status_warning` if solving was computed, but solution is not within tolerance + * and `oms_status_error` if an error occured. + */ +oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& graph) +{ + /* Update user data */ + KINSOL_USER_DATA* kinsolUserData = (KINSOL_USER_DATA*) userData; + kinsolUserData->syst = &syst; + kinsolUserData->graph = &graph; + AlgLoop* algLoop = syst.getAlgLoop(kinsolUserData->algLoopNumber); + const oms_ssc_t SCC = algLoop->getSCC(); + + int flag; + double fNormValue; + + logDebug("Solving system " + std::to_string(kinsolUserData->algLoopNumber)); + + if (SCC.size() != size) + { + logError("The size of the loop changed! This shouldn't be possible..."); + throw("Serious problem encountered. Open a ticket!"); + } + + /* Set initial guess */ + double *initialGuess_data = NV_DATA_S(initialGuess); + for (int i=0; i < size; i++) + { + int output = SCC[i].first; + if (oms_status_ok != syst.getReal(graph.getNodes()[output].getName(), initialGuess_data[i])) + { + return oms_status_error; + } + } + + /* u and f scaling */ + // TODO: Add scaling that is not only constant ones + + /* Solve algebraic loop with KINSol() */ + flag = KINSol(kinsolMemory, /* KINSol memory block */ + initialGuess, /* initial guess on input; solution vector */ + KIN_NONE, /* global strategy choice: Basic newton iteration */ + uScale, /* scaling vector, for the variable uu */ + fScale); /* scaling vector for function values fval */ + if (flag < 0) + { + logError("SUNDIALS_ERROR: KINSol() failed with flag = " + std::to_string(flag)); + return oms_status_error; + } + else + { + logDebug("SUNDIALS_INFO: KINSol() succeded with flag = " + std::to_string(flag)); + } + + /* Check solution */ + flag = nlsKinsolResiduals(initialGuess, fTmp, userData); + fNormValue = N_VWL2Norm(fTmp, fTmp); + if ( fNormValue > fnormtol ) + { + logWarning("Solution of algebraic loop " + std::to_string(((KINSOL_USER_DATA *)userData)->algLoopNumber) + "not within precission given by fnormtol: " + std::to_string(fnormtol)); + logDebug("2-norm of residual of solution: " + std::to_string(fNormValue)); + return oms_status_warning; + } + + logDebug("Solved system " + std::to_string(kinsolUserData->algLoopNumber) + " successfully"); + return oms_status_ok; +} + +/** + * @brief Construct a new oms::AlgLoop::AlgLoop object + * + * @param method Specifies used solver for the loop. + * Can be `oms_alg_solver_fixedpoint` for fixed-point-iteration + * or `oms_alg_solver_kinsol` for SUNDIALS KINSOL + * @param absTol Tolerance used for the algebraic solver. + * @param scc Strong Connected Componten, a vector of connected + * @param number + */ +oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, oms_ssc_t scc, const int number): absoluteTolerance(absTol), SCC(scc), systNumber(number) +{ + switch (method) + { + case oms_alg_solver_fixedpoint: + case oms_alg_solver_kinsol: + algSolverMethod = method; + break; + default: + logError("Unknown alg. loop solver method"); + throw; + } + + if (method == oms_alg_solver_kinsol) + { + kinsolData = KinsolSolver::NewKinsolSolver(systNumber, SCC.size(), absoluteTolerance); + if (kinsolData==NULL) + { + logError("NewKinsolSolver() failed. Aborting!"); + throw("AlgLoop() failed!"); + } + } +} + +/** + * @brief Solve algebraic loop + * + * Using solver method saved during AlgLoop creation. + * Can use fixed-point-iteration and KINSOL. + * + * @param syst Reference to System + * @param graph Reference to directed graph + * @return oms_status_enu_t Returns `oms_status_ok` on success + */ +oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph) +{ + logDebug("Solving algebraic loop formed by connections\n" + dumpLoopVars(graph)); + logDebug("Using solver " + getAlgSolverName()); + + switch (algSolverMethod) + { + case oms_alg_solver_fixedpoint: + return fixPointIteration(syst, graph); + case oms_alg_solver_kinsol: + return kinsolData->kinsolSolve(syst, graph); + default: + logError("Invalid algebraic solver method!"); + return oms_status_error; + } +} + +/** + * @brief Fixed-point-iteration to solve algebraic loop. + * + * @param syst + * @param graph + * @return oms_status_enu_t + */ +oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& graph) +{ + const int size = SCC.size(); + const int maxIterations = Flags::MaxLoopIteration(); + int it=0; + double maxRes; + double *res = new double[size](); + + do + { + it++; + // get old values + for (int i=0; i maxRes) + maxRes = fabs(res[i]); + } + } while(maxRes > absoluteTolerance && it < maxIterations); + + delete[] res; + + if (it >= maxIterations) + { + return logError("max. number of iterations (" + std::to_string(maxIterations) + ") exceeded at time = " + std::to_string(syst.getTime())); + } + logDebug("CompositeModel::solveAlgLoop: maxRes: " + std::to_string(maxRes) + ", iterations: " + std::to_string(it) + " at time = " + std::to_string(syst.getTime())); + return oms_status_ok; +} + +/** + * @brief Return solver method + * + * @return std::string Name of solver method + */ +std::string oms::AlgLoop::getAlgSolverName(void) +{ + switch (algSolverMethod) + { + case oms_alg_solver_none: + return "None"; + case oms_alg_solver_fixedpoint: + return "Fixed-Point-Iteration"; + case oms_alg_solver_kinsol: + return "KINSOL"; + default: + logError("This should not be reachable!"); + return "Can't figure out what solver method is used!"; + } +} + +/** + * @brief Dump variables of algebraic loop. + * + * @param graph + * @return std::string + */ +std::string oms::AlgLoop::dumpLoopVars(DirectedGraph& graph) +{ + const int size = SCC.size(); + std::string varNames = "\t"; + int output; + + for (int i=0; i "); + output = SCC[i].second; + varNames.append(graph.getNodes()[output].getName().c_str()); + varNames.append("\n\t"); + } + output = SCC[size-1].first; + varNames.append(graph.getNodes()[output].getName().c_str()); + varNames.append(" -> "); + output = SCC[size-1].second; + varNames.append(graph.getNodes()[output].getName().c_str()); + + return varNames; +} diff --git a/src/OMSimulatorLib/AlgLoop.h b/src/OMSimulatorLib/AlgLoop.h new file mode 100644 index 000000000..b59a48d68 --- /dev/null +++ b/src/OMSimulatorLib/AlgLoop.h @@ -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 +#include +#include "Types.h" +#include "DirectedGraph.h" + +#include +#include +#include /* 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_ diff --git a/src/OMSimulatorLib/CMakeLists.txt b/src/OMSimulatorLib/CMakeLists.txt index be7752a81..a677486a3 100644 --- a/src/OMSimulatorLib/CMakeLists.txt +++ b/src/OMSimulatorLib/CMakeLists.txt @@ -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) diff --git a/src/OMSimulatorLib/DirectedGraph.cpp b/src/OMSimulatorLib/DirectedGraph.cpp index 9644177ef..a8d425a67 100644 --- a/src/OMSimulatorLib/DirectedGraph.cpp +++ b/src/OMSimulatorLib/DirectedGraph.cpp @@ -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 >& 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 > oms::DirectedGraph::getSCCs() return components; } -const std::vector< std::vector< std::pair > >& 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 > >& oms::DirectedGraph::get void oms::DirectedGraph::calculateSortedConnections() { std::deque< std::vector > components = getSCCs(); - std::vector< std::pair > SCC; + oms_ssc_t SCC; sortedConnections.clear(); for (int i = 0; i < components.size(); ++i) diff --git a/src/OMSimulatorLib/DirectedGraph.h b/src/OMSimulatorLib/DirectedGraph.h index 771d777f7..955022ba5 100644 --- a/src/OMSimulatorLib/DirectedGraph.h +++ b/src/OMSimulatorLib/DirectedGraph.h @@ -43,6 +43,13 @@ #include #include +/** + * @brief Strong connected components data type. + * + * A vector of pairs of connected components. + */ +typedef std::vector< std::pair > 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 > >& getSortedConnections(); + const std::vector< oms_ssc_t >& getSortedConnections(); const std::vector& getNodes() const {return nodes;} - const std::vector< std::pair >& getEdges() const {return edges;} + const oms_ssc_t& getEdges() const {return edges;} private: std::deque< std::vector > getSCCs(); void calculateSortedConnections(); void strongconnect(int v, std::vector< std::vector > G, int& index, int *d, int *low, std::stack& S, bool *stacked, std::deque< std::vector >& components); - static int getEdgeIndex(const std::vector< std::pair >& edges, int from, int to); + static int getEdgeIndex(const oms_ssc_t& edges, int from, int to); private: std::vector nodes; - std::vector< std::pair > edges; + oms_ssc_t edges; std::vector< std::vector > G; - std::vector< std::vector< std::pair > > sortedConnections; + std::vector< oms_ssc_t > sortedConnections; bool sortedConnectionsAreValid; }; } diff --git a/src/OMSimulatorLib/Flags.cpp b/src/OMSimulatorLib/Flags.cpp index cfd47089b..fc4a15d1c 100644 --- a/src/OMSimulatorLib/Flags.cpp +++ b/src/OMSimulatorLib/Flags.cpp @@ -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(); diff --git a/src/OMSimulatorLib/Flags.h b/src/OMSimulatorLib/Flags.h index 63e717d46..8bf8bf91a 100644 --- a/src/OMSimulatorLib/Flags.h +++ b/src/OMSimulatorLib/Flags.h @@ -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 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); diff --git a/src/OMSimulatorLib/System.cpp b/src/OMSimulatorLib/System.cpp index 503ca2793..50c4fbacb 100644 --- a/src/OMSimulatorLib/System.cpp +++ b/src/OMSimulatorLib/System.cpp @@ -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 1) + { + addAlgLoop(sortedConnections[i], systCount); + systCount++; + } + } + loopsNeedUpdate = false; + } + + return oms_status_ok; +} + oms_status_enu_t oms::System::importParameterMappingFromSSM(std::string fileName, std::multimap &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); +} diff --git a/src/OMSimulatorLib/System.h b/src/OMSimulatorLib/System.h index 8e4664bc3..b7c24d90f 100644 --- a/src/OMSimulatorLib/System.h +++ b/src/OMSimulatorLib/System.h @@ -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 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 resultFileMapping; std::unordered_map 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 &mappedEntry); oms_status_enu_t importParameterMappingFromSSM(std::string fileName, std::multimap &mappedEntry); + + std::vector algLoops; ///< Vector of algebraic loop objects }; } diff --git a/src/OMSimulatorLib/SystemSC.cpp b/src/OMSimulatorLib/SystemSC.cpp index 232259f92..cf1ffbda8 100644 --- a/src/OMSimulatorLib/SystemSC.cpp +++ b/src/OMSimulatorLib/SystemSC.cpp @@ -345,6 +345,9 @@ oms_status_enu_t oms::SystemSC::initialize() if (flag < 0) logError("SUNDIALS_ERROR: CVodeSetMaxNumSteps() failed with flag = " + std::to_string(flag)); } + // Mark algebraic loops to be updated on next call + loopsNeedUpdate = true; + return oms_status_ok; } @@ -651,6 +654,8 @@ oms_status_enu_t oms::SystemSC::stepUntil(double stopTime, void (*cb)(const char oms_status_enu_t oms::SystemSC::updateInputs(DirectedGraph& graph) { CallClock callClock(clock); + oms_status_enu_t status; + int loopNum = 0; if (getModel()->validState(oms_modelState_simulation)) { @@ -673,7 +678,9 @@ oms_status_enu_t oms::SystemSC::updateInputs(DirectedGraph& graph) } // input := output - const std::vector< std::vector< std::pair > >& sortedConnections = graph.getSortedConnections(); + const std::vector< oms_ssc_t >& sortedConnections = graph.getSortedConnections(); + updateAlgebraicLoops(sortedConnections); + for(int i=0; i >& SCC) -{ - CallClock callClock(clock); - - const int size = SCC.size(); - const int maxIterations = Flags::MaxLoopIteration(); - double maxRes; - double *res = new double[size](); - - int it=0; - do - { - it++; - // get old values - for (int i=0; i maxRes) - maxRes = fabs(res[i]); - } - } while(maxRes > absoluteTolerance && it < maxIterations); - - delete[] res; - - if (it >= maxIterations) - return logError("max. number of iterations (" + std::to_string(maxIterations) + ") exceeded at time = " + std::to_string(getTime())); - logDebug("CompositeModel::solveAlgLoop: maxRes: " + std::to_string(maxRes) + ", iterations: " + std::to_string(it) + " at time = " + std::to_string(getTime())); + } return oms_status_ok; } diff --git a/src/OMSimulatorLib/SystemSC.h b/src/OMSimulatorLib/SystemSC.h index 56472f025..b9ec5fbda 100644 --- a/src/OMSimulatorLib/SystemSC.h +++ b/src/OMSimulatorLib/SystemSC.h @@ -61,10 +61,7 @@ namespace oms oms_status_enu_t reset(); oms_status_enu_t stepUntil(double stopTime, void (*cb)(const char* ident, double time, oms_status_enu_t status)); - double getTime() const {return time;} - oms_status_enu_t updateInputs(DirectedGraph& graph); - oms_status_enu_t solveAlgLoop(DirectedGraph& graph, const std::vector< std::pair >& SCC); std::string getSolverName() const; oms_status_enu_t setSolverMethod(std::string); @@ -79,8 +76,6 @@ namespace oms SystemSC& operator=(SystemSC const& copy); ///< not implemented private: - double time; - std::vector fmus; std::vector callEventUpdate; diff --git a/src/OMSimulatorLib/SystemWC.cpp b/src/OMSimulatorLib/SystemWC.cpp index 02ada14f9..9f65db03d 100644 --- a/src/OMSimulatorLib/SystemWC.cpp +++ b/src/OMSimulatorLib/SystemWC.cpp @@ -216,6 +216,9 @@ oms_status_enu_t oms::SystemWC::initialize() if (solverMethod == oms_solver_wc_mav || solverMethod == oms_solver_wc_mav2) stepSize = initialStepSize; + // Mark algebraic loops to be updated on next call + loopsNeedUpdate = true; + return oms_status_ok; } @@ -916,7 +919,7 @@ oms_status_enu_t oms::SystemWC::setRealInputDerivative(const ComRef& cref, const oms_status_enu_t oms::SystemWC::getInputs(oms::DirectedGraph& graph, std::vector& inputs) { inputs.clear(); - const std::vector< std::vector< std::pair > >& sortedConnections = graph.getSortedConnections(); + const std::vector< oms_ssc_t >& sortedConnections = graph.getSortedConnections(); for(int i=0; i& inputsDer) { int derI = 0; - const std::vector< std::vector< std::pair > >& sortedConnections = graph.getSortedConnections(); + const std::vector< oms_ssc_t >& sortedConnections = graph.getSortedConnections(); for(int i=0; i& inputVect,std::vector& outputVect,std::map FMUcomponents) { // FMUcomponents in will be list of FMUs that CAN GET FMUs - const std::vector< std::vector< std::pair > >& sortedConnections = graph.getSortedConnections(); + const std::vector< oms_ssc_t >& sortedConnections = graph.getSortedConnections(); inputVect.clear(); int inCount = 0; outputVect.clear(); @@ -1011,9 +1014,13 @@ oms_status_enu_t oms::SystemWC::getInputAndOutput(oms::DirectedGraph& graph, std oms_status_enu_t oms::SystemWC::updateInputs(oms::DirectedGraph& graph) { CallClock callClock(clock); + oms_status_enu_t status; + int loopNum = 0; // input := output - const std::vector< std::vector< std::pair > >& sortedConnections = graph.getSortedConnections(); + const std::vector< oms_ssc_t >& sortedConnections = graph.getSortedConnections(); + updateAlgebraicLoops(sortedConnections); + for(int i=0; i >& SCC) -{ - CallClock callClock(clock); - - const int size = SCC.size(); - const int maxIterations = Flags::MaxLoopIteration(); - double maxRes; - double *res = new double[size](); - - int it=0; - do - { - it++; - // get old values - for (int i=0; i maxRes) - maxRes = fabs(res[i]); + loopNum++; } - } while(maxRes > absoluteTolerance && it < maxIterations); - - delete[] res; - - if (it >= maxIterations) - return logError("max. number of iterations (" + std::to_string(maxIterations) + ") exceeded at time = " + std::to_string(getTime())); - logDebug("CompositeModel::solveAlgLoop: maxRes: " + std::to_string(maxRes) + ", iterations: " + std::to_string(it) + " at time = " + std::to_string(getTime())); + } return oms_status_ok; } diff --git a/src/OMSimulatorLib/SystemWC.h b/src/OMSimulatorLib/SystemWC.h index 62a13017d..f0198b9e2 100644 --- a/src/OMSimulatorLib/SystemWC.h +++ b/src/OMSimulatorLib/SystemWC.h @@ -60,8 +60,6 @@ namespace oms oms_status_enu_t stepUntil(double stopTime, void (*cb)(const char* ident, double time, oms_status_enu_t status)); oms_status_enu_t stepUntilASSC(double stopTime, void (*cb)(const char* ident, double time, oms_status_enu_t status)); - double getTime() const {return time;} - std::string getSolverName() const; oms_status_enu_t setSolverMethod(std::string); oms_status_enu_t setSolver(oms_solver_enu_t solver) {if (solver > oms_solver_wc_min && solver < oms_solver_wc_max) {solverMethod=solver; return oms_status_ok;} return oms_status_error;} @@ -70,7 +68,6 @@ namespace oms oms_status_enu_t setInputsDer(oms::DirectedGraph& graph, const std::vector& inputsDer); oms_status_enu_t getInputAndOutput(DirectedGraph& graph, std::vector& inputVect,std::vector& outputVect,std::map FMUcomponents); oms_status_enu_t updateInputs(DirectedGraph& graph); - oms_status_enu_t solveAlgLoop(DirectedGraph& graph, const std::vector< std::pair >& SCC); oms_status_enu_t getRealOutputDerivative(const ComRef& cref, SignalDerivative& der); oms_status_enu_t setRealInputDerivative(const ComRef& cref, const SignalDerivative& der); @@ -93,8 +90,6 @@ namespace oms SystemWC& operator=(SystemWC const& copy); ///< not implemented private: - double time; - StepSizeConfiguration stepSizeConfiguration; ///< Configuration data structure for assc unsigned int h_id; diff --git a/src/OMSimulatorLib/Types.h b/src/OMSimulatorLib/Types.h index ab1a62950..0524a63f9 100644 --- a/src/OMSimulatorLib/Types.h +++ b/src/OMSimulatorLib/Types.h @@ -107,6 +107,12 @@ typedef enum { oms_solver_wc_max } oms_solver_enu_t; +typedef enum { + oms_alg_solver_none, + oms_alg_solver_fixedpoint, ///< Fixed-point-iteration (default) + oms_alg_solver_kinsol ///< Kinsol solver +} oms_alg_solver_enu_t; + typedef enum { oms_element_system, oms_element_component,