diff --git a/include/OMSimulator/Types.h b/include/OMSimulator/Types.h index 33121e3f5..fe4fa0401 100644 --- a/include/OMSimulator/Types.h +++ b/include/OMSimulator/Types.h @@ -76,6 +76,7 @@ typedef enum { oms_solver_sc_min, oms_solver_sc_explicit_euler, oms_solver_sc_cvode, ///< default + oms_solver_sc_ida, ///< DAE mode: SUNDIALS IDA over the components' residuals (fmi-ls-dae) oms_solver_sc_max, oms_solver_wc_min, oms_solver_wc_ma, ///< Fixed stepsize (default) diff --git a/src/OMSimulatorLib/CMakeLists.txt b/src/OMSimulatorLib/CMakeLists.txt index 295f379f4..a126e5967 100644 --- a/src/OMSimulatorLib/CMakeLists.txt +++ b/src/OMSimulatorLib/CMakeLists.txt @@ -68,6 +68,7 @@ set(OMSIMULATORLIB_SOURCES Flags.cpp FMUInfo.cpp Logging.cpp + LsDae.cpp MatReader.cpp MatVer4.cpp MATWriter.cpp @@ -122,6 +123,8 @@ target_link_libraries(OMSimulatorLib oms::public_includes oms::3rd::sundials::kinsol oms::3rd::sundials::cvode + oms::3rd::sundials::ida + oms::3rd::sundials::sunlinsolklu oms::3rd::fmi4c oms::3rd::minizip oms::3rd::zlib @@ -155,6 +158,8 @@ target_link_libraries(OMSimulatorLib_static oms::public_includes oms::3rd::sundials::kinsol oms::3rd::sundials::cvode + oms::3rd::sundials::ida + oms::3rd::sundials::sunlinsolklu oms::3rd::fmi4c oms::3rd::minizip oms::3rd::zlib diff --git a/src/OMSimulatorLib/ComponentFMU3CS.cpp b/src/OMSimulatorLib/ComponentFMU3CS.cpp index c6492c5d8..54fee4b01 100644 --- a/src/OMSimulatorLib/ComponentFMU3CS.cpp +++ b/src/OMSimulatorLib/ComponentFMU3CS.cpp @@ -474,6 +474,22 @@ void oms::ComponentFMU3CS::dumpInitialUnknowns() logInfo("[" + std::string(getCref()) + ": " + getPath() + "] The FMU contains " + std::to_string(n) + " initial unknowns: " + str); } +/** + * \\brief The variable a ModelStructure dependency names. + * + * FMI 3.0 changed this from FMI 2.0: a `dependencies` list holds **value + * references**, not one-based indices into the variable list. Reading them as + * indices silently wires the dependency graph to the wrong variables, and fails + * outright on value reference 0, which is a perfectly ordinary one. + */ +const oms::Variable* oms::ComponentFMU3CS::variableByFMI3ValueReference(unsigned int vr) const +{ + for (const auto& v : allVariables) + if (v.getValueReferenceFMI3() == vr) + return &v; + return nullptr; +} + oms_status_enu_t oms::ComponentFMU3CS::initializeDependencyGraph_initialUnknowns() { if (initialUnknownsGraph.getEdges().connections.size() > 0) @@ -560,15 +576,16 @@ oms_status_enu_t oms::ComponentFMU3CS::initializeDependencyGraph_initialUnknowns else { //dependency exist - for (const auto &index : it.second) + for (const auto &vr : it.second) { - if (index < 1 || index > allVariables.size()) + const Variable* dependency = variableByFMI3ValueReference(vr); + if (!dependency) { - logWarning("Initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " has bad dependency on variable with index " + std::to_string(index) + " which couldn't be resolved"); + logWarning("Initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " has a dependency on value reference " + std::to_string(vr) + " which couldn't be resolved"); return logError(std::string(getCref()) + ": Erroneous initial unknowns detected in modelDescription.xml\nUse flag --ignoreInitialUnknowns=true to ignore all initial unknowns, but this can cause inflated loop size."); } - logDebug(std::string(getCref()) + ": " + getPath() + " initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " depends on " + std::string(allVariables[index - 1])); - initialUnknownsGraph.addEdge(allVariables[index - 1].makeConnector(this->getFullCref()), initialUnknownsGraph.getNodes()[i]); + logDebug(std::string(getCref()) + ": " + getPath() + " initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " depends on " + std::string(*dependency)); + initialUnknownsGraph.addEdge(dependency->makeConnector(this->getFullCref()), initialUnknownsGraph.getNodes()[i]); } } i = i + 1; @@ -610,15 +627,16 @@ oms_status_enu_t oms::ComponentFMU3CS::initializeDependencyGraph_outputs() } else { - for (const auto &index : it.second) + for (const auto &vr : it.second) { - if (index < 1 || index > allVariables.size()) + const Variable* dependency = variableByFMI3ValueReference(vr); + if (!dependency) { - logWarning("Output " + std::string(output) + " has bad dependency on variable with index " + std::to_string(index) + " which couldn't be resolved"); + logWarning("Output " + std::string(output) + " has a dependency on value reference " + std::to_string(vr) + " which couldn't be resolved"); return logError(std::string(getCref()) + ": erroneous dependencies detected in modelDescription.xml"); } - logDebug(std::string(getCref()) + ": " + getPath() + " output " + std::string(output) + " depends on " + std::string(allVariables[index - 1])); - outputsGraph.addEdge(allVariables[index - 1].makeConnector(this->getFullCref()), output.makeConnector(this->getFullCref())); + logDebug(std::string(getCref()) + ": " + getPath() + " output " + std::string(output) + " depends on " + std::string(*dependency)); + outputsGraph.addEdge(dependency->makeConnector(this->getFullCref()), output.makeConnector(this->getFullCref())); } } i = i + 1; diff --git a/src/OMSimulatorLib/ComponentFMU3CS.h b/src/OMSimulatorLib/ComponentFMU3CS.h index 5aa2317bc..e012749c6 100644 --- a/src/OMSimulatorLib/ComponentFMU3CS.h +++ b/src/OMSimulatorLib/ComponentFMU3CS.h @@ -77,6 +77,8 @@ namespace oms oms_status_enu_t stepUntil(double stopTime); + /// The variable a ModelStructure dependency names (FMI 3.0: a value reference). + const Variable* variableByFMI3ValueReference(unsigned int vr) const; oms_status_enu_t initializeDependencyGraph_initialUnknowns(); oms_status_enu_t initializeDependencyGraph_outputs(); diff --git a/src/OMSimulatorLib/ComponentFMU3ME.cpp b/src/OMSimulatorLib/ComponentFMU3ME.cpp index 53d0588dd..9dbb5ddc4 100644 --- a/src/OMSimulatorLib/ComponentFMU3ME.cpp +++ b/src/OMSimulatorLib/ComponentFMU3ME.cpp @@ -188,6 +188,23 @@ oms::Component* oms::ComponentFMU3ME::NewComponent(const oms::ComRef& cref, oms: return NULL; } + // fmi-ls-dae: a Model Exchange FMU may carry a DAE formulation beside its ODE + // face. Nothing changes until the system asks for DAE mode; here we only learn + // that the FMU has one, and remember the state derivatives it pairs with its + // states (the model description's own order, which + // is the order fmi3GetContinuousStates uses). + if (component->lsDae.load(tempDir.generic_string())) + { + const int n = fmi3_getNumberOfModelStructureContinuousStateDerivatives(component->fmu); + for (int i = 0; i < n; ++i) + component->derivativeVrs.push_back( + fmi3_getModelStructureValueReference(fmi3_getModelStructureContinuousStateDerivative(component->fmu, i))); + logInfo("FMU \"" + std::string(cref) + "\" carries a DAE formulation (fmi-ls-dae " + + component->lsDae.getVersion() + "): " + std::to_string(component->lsDae.getResiduals().size()) + + " residuals over " + std::to_string(component->derivativeVrs.size()) + " states and " + + std::to_string(component->lsDae.getAlgebraicVariables().size()) + " algebraic variables"); + } + // update FMU info component->fmuInfo.update(oms_component_fmu3, component->fmu); component->omsfmi3logger = oms::fmi3logger; @@ -472,6 +489,22 @@ void oms::ComponentFMU3ME::dumpInitialUnknowns() logInfo("[" + std::string(getCref()) + ": " + getPath() + "] The FMU contains " + std::to_string(n) + " initial unknowns: " + str); } +/** + * \\brief The variable a ModelStructure dependency names. + * + * FMI 3.0 changed this from FMI 2.0: a `dependencies` list holds **value + * references**, not one-based indices into the variable list. Reading them as + * indices silently wires the dependency graph to the wrong variables, and fails + * outright on value reference 0, which is a perfectly ordinary one. + */ +const oms::Variable* oms::ComponentFMU3ME::variableByFMI3ValueReference(unsigned int vr) const +{ + for (const auto& v : allVariables) + if (v.getValueReferenceFMI3() == vr) + return &v; + return nullptr; +} + oms_status_enu_t oms::ComponentFMU3ME::initializeDependencyGraph_initialUnknowns() { if (initialUnknownsGraph.getEdges().connections.size() > 0) @@ -558,15 +591,16 @@ oms_status_enu_t oms::ComponentFMU3ME::initializeDependencyGraph_initialUnknowns else { //dependency exist - for (const auto &index : it.second) + for (const auto &vr : it.second) { - if (index < 1 || index > allVariables.size()) + const Variable* dependency = variableByFMI3ValueReference(vr); + if (!dependency) { - logWarning("Initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " has bad dependency on variable with index " + std::to_string(index) + " which couldn't be resolved"); + logWarning("Initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " has a dependency on value reference " + std::to_string(vr) + " which couldn't be resolved"); return logError(std::string(getCref()) + ": Erroneous initial unknowns detected in modelDescription.xml\nUse flag --ignoreInitialUnknowns=true to ignore all initial unknowns, but this can cause inflated loop size."); } - logDebug(std::string(getCref()) + ": " + getPath() + " initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " depends on " + std::string(allVariables[index - 1])); - initialUnknownsGraph.addEdge(allVariables[index - 1].makeConnector(this->getFullCref()), initialUnknownsGraph.getNodes()[i]); + logDebug(std::string(getCref()) + ": " + getPath() + " initial unknown " + std::string(initialUnknownsGraph.getNodes()[i]) + " depends on " + std::string(*dependency)); + initialUnknownsGraph.addEdge(dependency->makeConnector(this->getFullCref()), initialUnknownsGraph.getNodes()[i]); } } i = i + 1; @@ -608,15 +642,16 @@ oms_status_enu_t oms::ComponentFMU3ME::initializeDependencyGraph_outputs() } else { - for (const auto &index : it.second) + for (const auto &vr : it.second) { - if (index < 1 || index > allVariables.size()) + const Variable* dependency = variableByFMI3ValueReference(vr); + if (!dependency) { - logWarning("Output " + std::string(output) + " has bad dependency on variable with index " + std::to_string(index) + " which couldn't be resolved"); + logWarning("Output " + std::string(output) + " has a dependency on value reference " + std::to_string(vr) + " which couldn't be resolved"); return logError(std::string(getCref()) + ": erroneous dependencies detected in modelDescription.xml"); } - logDebug(std::string(getCref()) + ": " + getPath() + " output " + std::string(output) + " depends on " + std::string(allVariables[index - 1])); - outputsGraph.addEdge(allVariables[index - 1].makeConnector(this->getFullCref()), output.makeConnector(this->getFullCref())); + logDebug(std::string(getCref()) + ": " + getPath() + " output " + std::string(output) + " depends on " + std::string(*dependency)); + outputsGraph.addEdge(dependency->makeConnector(this->getFullCref()), output.makeConnector(this->getFullCref())); } } i = i + 1; @@ -635,6 +670,13 @@ oms_status_enu_t oms::ComponentFMU3ME::instantiate() exit(1); } + // fmi-ls-dae: the switch is a structural parameter, so Configuration Mode is + // where it goes — before the start values and before initialization. An FMU + // that declares a DAE formulation is used as one: its residuals are what it + // has, and its ODE face need not even be implemented. + if (lsDae.isValid() && oms_status_ok != enableDaeMode()) + return oms_status_error; + // set start values from local resources if (values.hasResources()) { @@ -883,6 +925,10 @@ oms_status_enu_t oms::ComponentFMU3ME::initialize() fmistatus = fmi3_getNumberOfEventIndicators(instance, &nEventIndicators); if (fmi3OK != fmistatus) return logError_FMUCall("fmi3_getNumberOfEventIndicators", this); + // fmi-ls-dae: now that the state count is known, the DAE system can be checked + if (oms_status_ok != validateDaeMode()) + return oms_status_error; + // fmi3_exitInitialization_mode leaves FMU in event mode if (oms_status_ok != doEventIteration()) return oms_status_error; @@ -2055,6 +2101,147 @@ oms_status_enu_t oms::ComponentFMU3ME::getDerivatives(double* derivatives) return oms_status_ok; } +/** + * \brief fmi-ls-dae: switch the FMU into DAE mode. + * + * The switch is a structural parameter, so Configuration Mode is where it is + * set — before initialization, and only once. After this the FMU no longer owes + * the derivatives and the algebraic variables: the master sets them and reads + * the residuals of F(t, x, x', z) = 0 back. + */ +oms_status_enu_t oms::ComponentFMU3ME::enableDaeMode() +{ + CallClock callClock(clock); + + if (!lsDae.isValid()) + return logError("FMU \"" + std::string(getFullCref()) + "\" carries no fmi-ls-dae manifest, so it has no DAE mode"); + + if (daeMode) + return oms_status_ok; + + if (!lsDae.isFullyImplicit()) + return logError("fmi-ls-dae: FMU \"" + std::string(getFullCref()) + + "\" states a semi-explicit DAE, which is not supported yet"); + + // The number of continuous states is only known once the FMU has left + // Initialization Mode, which is long after Configuration Mode; the checks that + // need it are in validateDaeMode(), called from initialize(). + + if (fmi3OK != fmi3_enterConfigurationMode(instance)) + return logError_FMUCall("fmi3_enterConfigurationMode", this); + + const fmi3ValueReference vr = lsDae.getEnableValueReference(); + const fmi3Boolean on = fmi3True; + if (fmi3OK != fmi3_setBoolean(instance, &vr, 1, &on, 1)) + { + fmi3_exitConfigurationMode(instance); + return logError_FMUCall("fmi3_setBoolean (the fmi-ls-dae DAE-mode parameter)", this); + } + + if (fmi3OK != fmi3_exitConfigurationMode(instance)) + return logError_FMUCall("fmi3_exitConfigurationMode", this); + + daeMode = true; + return oms_status_ok; +} + +/** + * \brief The fmi-ls-dae checks that need the state count, which the FMU only + * reports once it has left Initialization Mode. + */ +oms_status_enu_t oms::ComponentFMU3ME::validateDaeMode() +{ + if (!daeMode) + return oms_status_ok; + + // The state each derivative belongs to: a residual's dependency on a state and + // on its derivative are the same column of the Jacobian, because one + // difference quotient carries dF/dx + cj*dF/dder(x) together. + // The `derivative` attribute is a value reference in FMI 3.0, where FMI 2.0 had + // an index — the same change as the ModelStructure dependencies. fmi4c hands + // the raw attribute through under its FMI 2.0 name, so what Variable calls the + // state index is the state's value reference here. + stateVrs.clear(); + for (const fmi3ValueReference vr : derivativeVrs) + { + const Variable* der = variableByFMI3ValueReference(vr); + const Variable* state = der ? variableByFMI3ValueReference(der->getStateIndex()) : nullptr; + if (!state) + { + logWarning("fmi-ls-dae: FMU \"" + std::string(getFullCref()) + "\" does not say which state the derivative at value reference " + + std::to_string(vr) + " belongs to; the DAE Jacobian will be treated as dense"); + stateVrs.clear(); + break; + } + stateVrs.push_back(state->getValueReferenceFMI3()); + } + + if (derivativeVrs.size() != getNumberOfContinuousStates()) + return logError("fmi-ls-dae: FMU \"" + std::string(getFullCref()) + "\" has " + + std::to_string(getNumberOfContinuousStates()) + " continuous states but its lists " + + std::to_string(derivativeVrs.size()) + " state derivatives"); + + // The implicit form's residuals cover the state rows too, so a square system + // wants one residual per state and per algebraic variable. + const size_t wanted = getNumberOfContinuousStates() + lsDae.getAlgebraicVariables().size(); + if (lsDae.getResiduals().size() != wanted) + return logError("fmi-ls-dae: FMU \"" + std::string(getFullCref()) + "\" states an implicit DAE with " + + std::to_string(getNumberOfContinuousStates()) + " states and " + + std::to_string(lsDae.getAlgebraicVariables().size()) + " algebraic variables, which wants " + + std::to_string(wanted) + " residuals, but lists " + std::to_string(lsDae.getResiduals().size()) + + "; only a square system can be integrated"); + + return oms_status_ok; +} + +oms_status_enu_t oms::ComponentFMU3ME::getDaeResiduals(double* residuals) +{ + CallClock callClock(clock); + const std::vector& vrs = lsDae.getResiduals(); + if (vrs.empty()) + return oms_status_ok; + fmi3Status fmistatus = fmi3_getFloat64(instance, vrs.data(), vrs.size(), residuals, vrs.size()); + if (fmi3OK != fmistatus) + return logError_FMUCall("fmi3_getFloat64 (fmi-ls-dae residuals)", this); + return oms_status_ok; +} + +oms_status_enu_t oms::ComponentFMU3ME::getAlgebraicVariables(double* values) +{ + CallClock callClock(clock); + const std::vector& vrs = lsDae.getAlgebraicVariables(); + if (vrs.empty()) + return oms_status_ok; + fmi3Status fmistatus = fmi3_getFloat64(instance, vrs.data(), vrs.size(), values, vrs.size()); + if (fmi3OK != fmistatus) + return logError_FMUCall("fmi3_getFloat64 (fmi-ls-dae algebraic variables)", this); + return oms_status_ok; +} + +oms_status_enu_t oms::ComponentFMU3ME::setAlgebraicVariables(const double* values) +{ + CallClock callClock(clock); + const std::vector& vrs = lsDae.getAlgebraicVariables(); + if (vrs.empty()) + return oms_status_ok; + fmi3Status fmistatus = fmi3_setFloat64(instance, vrs.data(), vrs.size(), values, vrs.size()); + if (fmi3OK != fmistatus) + return logError_FMUCall("fmi3_setFloat64 (fmi-ls-dae algebraic variables)", this); + return oms_status_ok; +} + +oms_status_enu_t oms::ComponentFMU3ME::setDerivatives(const double* derivativeValues) +{ + CallClock callClock(clock); + if (derivativeVrs.empty()) + return oms_status_ok; + fmi3Status fmistatus = fmi3_setFloat64(instance, derivativeVrs.data(), derivativeVrs.size(), + derivativeValues, derivativeVrs.size()); + if (fmi3OK != fmistatus) + return logError_FMUCall("fmi3_setFloat64 (state derivatives)", this); + return oms_status_ok; +} + oms_status_enu_t oms::ComponentFMU3ME::getNominalsOfContinuousStates(double* nominals) { CallClock callClock(clock); diff --git a/src/OMSimulatorLib/ComponentFMU3ME.h b/src/OMSimulatorLib/ComponentFMU3ME.h index 5bd51a9f6..9cbd6aeb5 100644 --- a/src/OMSimulatorLib/ComponentFMU3ME.h +++ b/src/OMSimulatorLib/ComponentFMU3ME.h @@ -38,6 +38,7 @@ #include "Component.h" #include "ComRef.h" +#include "LsDae.h" #include "ResultWriter.h" #include "Snapshot.h" #include "Values.h" @@ -73,6 +74,8 @@ namespace oms oms_status_enu_t terminate(); oms_status_enu_t reset(); + /// The variable a ModelStructure dependency names (FMI 3.0: a value reference). + const Variable* variableByFMI3ValueReference(unsigned int vr) const; oms_status_enu_t initializeDependencyGraph_initialUnknowns(); oms_status_enu_t initializeDependencyGraph_outputs(); @@ -111,6 +114,35 @@ namespace oms size_t getNumberOfContinuousStates() const {return nContinuousStates;} size_t getNumberOfEventIndicators() const {return nEventIndicators;} + + /// fmi-ls-dae: what the FMU's layered-standard manifest declares, if any. + const LsDaeManifest& getLsDae() const {return lsDae;} + bool hasDaeFormulation() const {return lsDae.isValid();} + /// Whether this instance was switched into DAE mode. + bool isInDaeMode() const {return daeMode;} + /** + * \brief Switch the FMU into fmi-ls-dae's DAE mode. + * + * The switch is a structural parameter, so it is set in Configuration Mode, + * before initialization. From then on the FMU no longer answers with the + * derivatives and the algebraic variables: the importer owns them and reads + * the residuals of F(t, x, x', z) = 0 back. + */ + oms_status_enu_t enableDaeMode(); + /// The DAE checks that need the state count, which is known after initialization. + oms_status_enu_t validateDaeMode(); + size_t getNumberOfDaeResiduals() const {return lsDae.getResiduals().size();} + size_t getNumberOfAlgebraicVariables() const {return lsDae.getAlgebraicVariables().size();} + oms_status_enu_t getDaeResiduals(double* residuals); + /// The value references of the continuous states, in the order + /// fmi3GetContinuousStates uses — that of . + const std::vector& getStateValueReferences() const {return stateVrs;} + /// The value references of the state derivatives, in the same order. + const std::vector& getDerivativeValueReferences() const {return derivativeVrs;} + oms_status_enu_t getAlgebraicVariables(double* values); + oms_status_enu_t setAlgebraicVariables(const double* values); + /// The state derivatives, by value reference: the implicit form's knowns. + oms_status_enu_t setDerivatives(const double* derivatives); oms_status_enu_t getContinuousStates(double* states); oms_status_enu_t setContinuousStates(double* states); oms_status_enu_t getDerivatives(double* derivatives); @@ -161,6 +193,16 @@ namespace oms size_t nEventIndicators; size_t nContinuousStates; + /// fmi-ls-dae: the manifest, and whether this instance runs in DAE mode. + LsDaeManifest lsDae; + bool daeMode = false; + /// The value references of the model + /// description, in the order fmi3GetContinuousStates uses for the states. + std::vector derivativeVrs; + /// The states those derivatives belong to, resolved from the derivative + /// links of the model description; same order. + std::vector stateVrs; + fmi3Boolean newDiscreteStatesNeeded; fmi3Boolean terminateSimulation; fmi3Boolean nominalsOfContinuousStatesChanged; diff --git a/src/OMSimulatorLib/LsDae.cpp b/src/OMSimulatorLib/LsDae.cpp new file mode 100644 index 000000000..f6bfafbe0 --- /dev/null +++ b/src/OMSimulatorLib/LsDae.cpp @@ -0,0 +1,114 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2026, 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 AGPL VERSION 3 LICENSE OR + * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GNU AGPL + * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: + * http://www.openmodelica.org or + * https://github.com/OpenModelica/ or + * http://www.ida.liu.se/projects/OpenModelica, + * and in the OpenModelica distribution. + * + * GNU AGPL version 3 is obtained from: + * https://www.gnu.org/licenses/licenses.html#GPL + * + * 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 "LsDae.h" + +#include "Logging.h" +#include "OMSFileSystem.h" + +#include +#include + +/** + * \brief The value references in a whitespace-separated `dependencies` list. + */ +static std::vector parseValueReferenceList(const char* text) +{ + std::vector vrs; + if (!text) + return vrs; + std::istringstream in(text); + unsigned long vr; + while (in >> vr) + vrs.push_back(static_cast(vr)); + return vrs; +} + +bool oms::LsDaeManifest::load(const std::string& unzippedFmuPath) +{ + const filesystem::path manifest = filesystem::path(unzippedFmuPath) / filesystem::path(path()); + if (!filesystem::exists(manifest)) + return false; // no DAE formulation, which is the ordinary case + + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_file(manifest.string().c_str()); + if (!result) + { + logWarning("fmi-ls-dae: could not read \"" + manifest.string() + "\": " + std::string(result.description())); + return false; + } + + const pugi::xml_node root = doc.child("fmiDAEManifest"); + if (!root) + { + logWarning("fmi-ls-dae: \"" + manifest.string() + "\" has no element"); + return false; + } + + version = root.attribute("fmi-ls:fmi-ls-version").as_string(); + + const pugi::xml_node enable = root.child("EnableDAEParameter"); + if (!enable) + { + logWarning("fmi-ls-dae: the manifest names no , so the FMU cannot be switched into DAE mode"); + return false; + } + enableVr = enable.attribute("valueReference").as_uint(); + + for (const pugi::xml_node& node : root.child("AlgebraicVariables").children("AlgebraicVariable")) + algebraicVrs.push_back(node.attribute("valueReference").as_uint()); + + // The manifest's replaces the model description's: no + // is the fully implicit form, where the residuals + // cover the state rows as well. + const pugi::xml_node structure = root.child("ModelStructure"); + for (const pugi::xml_node& node : structure.children("ContinuousStateDerivative")) + continuousStateDerivatives.push_back(node.attribute("valueReference").as_uint()); + + for (const pugi::xml_node& node : structure.children("Residual")) + { + residualVrs.push_back(node.attribute("valueReference").as_uint()); + residualDependencies.push_back(parseValueReferenceList(node.attribute("dependencies").value())); + } + + if (residualVrs.empty()) + { + logWarning("fmi-ls-dae: the manifest lists no , so there is no DAE system to solve"); + return false; + } + + valid = true; + return true; +} diff --git a/src/OMSimulatorLib/LsDae.h b/src/OMSimulatorLib/LsDae.h new file mode 100644 index 000000000..309d9d064 --- /dev/null +++ b/src/OMSimulatorLib/LsDae.h @@ -0,0 +1,114 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2026, 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 AGPL VERSION 3 LICENSE OR + * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GNU AGPL + * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: + * http://www.openmodelica.org or + * https://github.com/OpenModelica/ or + * http://www.ida.liu.se/projects/OpenModelica, + * and in the OpenModelica distribution. + * + * GNU AGPL version 3 is obtained from: + * https://www.gnu.org/licenses/licenses.html#GPL + * + * 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_LSDAE_H_ +#define _OMS_LSDAE_H_ + +#include +#include +#include + +namespace oms +{ + /** + * \brief fmi-ls-dae: what an FMU's layered-standard manifest says about its + * DAE formulation. + * + * A Model Exchange FMU that carries + * `extra/org.fmi-standard.fmi-ls-dae/fmi-ls-manifest.xml` is an ordinary ODE + * FMU until the importer sets the structural parameter the manifest names, + * and a DAE FMU after: it then exposes the residuals of + * `F(t, x, x', z) = 0` and lets the importer solve for the algebraic + * variables `z` beside the states, instead of iterating on them internally. + * + * The manifest names the switch, the algebraic variables and the residuals; + * the state derivatives are the model description's own + * `` list, which orders them with the states. + */ + class LsDaeManifest + { + public: + /// The layered standard this reader understands. + static const char* name() { return "org.fmi-standard.fmi-ls-dae"; } + /// Where the manifest sits inside an FMU. + static const char* path() { return "extra/org.fmi-standard.fmi-ls-dae/fmi-ls-manifest.xml"; } + + /** + * \brief Read the manifest out of an unzipped FMU. + * + * \param unzippedFmuPath Directory the FMU was unzipped into. + * \return true if a manifest is there and could be read. A missing file is + * not an error: the FMU simply has no DAE formulation. + */ + bool load(const std::string& unzippedFmuPath); + + bool isValid() const { return valid; } + /// The version the manifest declares, e.g. "1.0.0-alpha.1". + const std::string& getVersion() const { return version; } + /// The structural parameter that switches the FMU into DAE mode. + uint32_t getEnableValueReference() const { return enableVr; } + /// The algebraic variables the importer solves for, beside the states. + const std::vector& getAlgebraicVariables() const { return algebraicVrs; } + /// One value reference per residual of F. + const std::vector& getResiduals() const { return residualVrs; } + /** + * \brief Per residual, the value references it declares a dependency on; + * empty for a residual that named none, which the standard reads as + * a dependency on every known. + */ + const std::vector>& getResidualDependencies() const { return residualDependencies; } + /** + * \brief The state derivatives the manifest lists. + * + * Empty for the fully implicit form `F(t, x, x', z) = 0`, where the + * residuals cover the state rows too. A non-empty list is the semi-explicit + * form: the FMU answers with `der(x)` and the residuals constrain the + * algebraic variables alone. + */ + const std::vector& getContinuousStateDerivatives() const { return continuousStateDerivatives; } + bool isFullyImplicit() const { return continuousStateDerivatives.empty(); } + + private: + bool valid = false; + std::string version; + uint32_t enableVr = 0; + std::vector algebraicVrs; + std::vector residualVrs; + std::vector> residualDependencies; + std::vector continuousStateDerivatives; + }; +} + +#endif diff --git a/src/OMSimulatorLib/SystemSC.cpp b/src/OMSimulatorLib/SystemSC.cpp index 4376adf0d..4720250bd 100644 --- a/src/OMSimulatorLib/SystemSC.cpp +++ b/src/OMSimulatorLib/SystemSC.cpp @@ -44,6 +44,9 @@ #include "ssd/Tags.h" #include +#include +#include +#include #include #include #include @@ -152,6 +155,556 @@ int oms::cvode_roots(sunrealtype t, N_Vector y, sunrealtype *gout, void *user_da return 0; } +/* --------------------------------------------------------------------------- + * fmi-ls-dae: the system as one implicit DAE + * + * The unknowns are laid out per component, [x_i | z_i], and so are the rows: a + * component in DAE mode contributes the residuals it declares, an ODE component + * the explicit rows x'_i - f_i(x, u, t). Both blocks are square, so the system + * is, and a system may mix the two. + * + * The coupling is where it was for CVODE: updateInputs() propagates connected + * outputs into inputs before the rows are evaluated, so the residuals see a + * consistent coupling rather than the value left over from the last accepted + * step. + * ------------------------------------------------------------------------- */ + +/** + * \brief Which algebraic loops of the simulation graph become rows of the + * global system, and what their unknowns are. + * + * Every connection of a loop contributes one unknown, the value of the input it + * drives, and one row, output - input = 0 — the same system the inner Newton + * solved, handed to IDA instead. A loop is only lifted if all of its + * connections carry reals: an integer or boolean connection is not something a + * continuous unknown can stand for, so such a loop keeps its inner solver. + */ +std::vector oms::SystemSC::allLoopsOf(DirectedGraph& graph) const +{ + size_t loops = 0; + for (const scc_t& scc : graph.getSortedConnections()) + if (scc.thisIsALoop) + ++loops; + return std::vector(loops, true); +} + +bool oms::SystemSC::anyComponentInDaeMode() +{ + for (const auto& component : getComponents()) + { + const ComponentFMU3ME* fmu3 = dynamic_cast(component.second); + if (fmu3 && fmu3->isInDaeMode()) + return true; + } + return false; +} + +void oms::SystemSC::collectLiftedLoops() +{ + liftedConnections.clear(); + liftedLoops.clear(); + + const std::vector& sortedConnections = simulationGraph.getSortedConnections(); + for (const scc_t& scc : sortedConnections) + { + if (!scc.thisIsALoop) + continue; + + bool allReal = true; + for (const auto& connection : scc.connections) + if (simulationGraph.getNodes()[connection.second].getType() != oms_signal_type_real || + simulationGraph.getNodes()[connection.first].getType() != oms_signal_type_real) + allReal = false; + + liftedLoops.push_back(allReal); + if (!allReal) + { + logInfo("system \"" + std::string(getFullCref()) + "\": an algebraic loop over connections that are not all real keeps its own solver"); + continue; + } + + for (const auto& connection : scc.connections) + liftedConnections.push_back({connection.first, connection.second}); + } + + if (!liftedConnections.empty()) + logInfo("system \"" + std::string(getFullCref()) + "\": " + std::to_string(liftedConnections.size()) + + " connection(s) of " + std::to_string(std::count(liftedLoops.begin(), liftedLoops.end(), true)) + + " algebraic loop(s) are rows of the DAE system, not an inner solve"); +} + +/** + * \brief Does this component own the graph node, i.e. is the node one of its + * own variables? + */ +static bool componentOwnsNode(const oms::Component* component, const oms::ComRef& node, std::string& variableName) +{ + // The connection graph of a system names its nodes relative to that system + // ("A.u"), while a component knows itself by its full path + // ("model.root.A"). Try both, so this does not depend on which one a caller + // built the graph with. + const std::string name = std::string(node); + for (const std::string& prefix : {std::string(component->getFullCref()) + ".", + std::string(component->getCref()) + "."}) + { + if (name.rfind(prefix, 0) == 0 && name.size() > prefix.size()) + { + variableName = name.substr(prefix.size()); + return true; + } + } + return false; +} + +void oms::SystemSC::buildDaeJacobianSparsity() +{ + daeJacRowsByColumn.clear(); + daeJacColorOfColumn.clear(); + daeJacColumnsOfColor.clear(); + daeJacIncrement.clear(); + daeJacSparse = false; + + const size_t n = nDaeUnknowns; + if (0 == n) + return; + + // Where each component's unknowns sit, and which column a value reference of + // it stands for. A state and its derivative share a column: one difference + // quotient carries dF/dx + cj*dF/dder(x) together. + std::vector columnOfComponent(fmus.size(), 0); + std::vector> columnOfValueReference(fmus.size()); + const size_t loopBase = n - liftedConnections.size(); + + for (size_t i = 0, j = 0; i < fmus.size(); ++i) + { + columnOfComponent[i] = j; + if (daeFmus[i]) + { + const std::vector& stateVrs = daeFmus[i]->getStateValueReferences(); + const std::vector& derVrs = daeFmus[i]->getDerivativeValueReferences(); + if (stateVrs.size() != nStates[i]) + { + // The component could not say which state each derivative belongs to; + // without that a residual's dependency on a state cannot be placed. + logInfo("system \"" + std::string(getFullCref()) + "\": the DAE Jacobian of \"" + + std::string(fmus[i]->getFullCref()) + "\" has no usable structure; differencing it densely"); + return; + } + for (size_t k = 0; k < nStates[i]; ++k) + { + columnOfValueReference[i][stateVrs[k]] = static_cast(j + k); + columnOfValueReference[i][derVrs[k]] = static_cast(j + k); + } + const std::vector& algVrs = daeFmus[i]->getLsDae().getAlgebraicVariables(); + for (size_t k = 0; k < algVrs.size(); ++k) + columnOfValueReference[i][algVrs[k]] = static_cast(j + nStates[i] + k); + } + j += nStates[i] + nAlgebraic[i]; + } + + // A lifted loop's unknown is the value of an input of some component, so a + // residual that depends on that input depends on that column. + std::vector componentOfLoop(liftedConnections.size(), -1); + std::vector componentOfLoopOutput(liftedConnections.size(), -1); + for (size_t k = 0; k < liftedConnections.size(); ++k) + { + for (size_t i = 0; i < fmus.size(); ++i) + { + std::string variableName; + if (componentOwnsNode(fmus[i], simulationGraph.getNodes()[liftedConnections[k].input].getName(), variableName)) + { + componentOfLoop[k] = static_cast(i); + if (daeFmus[i]) + { + const Variable* v = daeFmus[i]->getVariable(ComRef(variableName)); + if (v) + columnOfValueReference[i][v->getValueReferenceFMI3()] = static_cast(loopBase + k); + } + } + if (componentOwnsNode(fmus[i], simulationGraph.getNodes()[liftedConnections[k].output].getName(), variableName)) + componentOfLoopOutput[k] = static_cast(i); + } + } + + // The loop unknowns that drive a component. A residual of it may depend on any + // of them, and the manifest does not say: the dependencies it declares are + // those of the FMU's own DAE Jacobian, whose columns are the states and the + // algebraic variables — an input is not one of them. It only becomes an + // unknown here, where the master carries it, so these columns are added to + // every row of the component whether the manifest names them or not. + auto incomingLoopColumns = [&](size_t i, std::vector& out) { + for (size_t k = 0; k < liftedConnections.size(); ++k) + if (componentOfLoop[k] == static_cast(i)) + out.push_back(static_cast(loopBase + k)); + }; + + // All the columns of one component, plus those: what a row is taken to reach + // when nothing more precise is known. + auto columnsAround = [&](size_t i, std::vector& out) { + for (size_t k = 0; k < nStates[i] + nAlgebraic[i]; ++k) + out.push_back(static_cast(columnOfComponent[i] + k)); + incomingLoopColumns(i, out); + }; + + std::vector> columnsByRow(n); + size_t row = 0; + for (size_t i = 0; i < fmus.size(); ++i) + { + if (daeFmus[i]) + { + const std::vector>& dependencies = daeFmus[i]->getLsDae().getResidualDependencies(); + for (size_t k = 0; k < nStates[i] + nAlgebraic[i]; ++k, ++row) + { + bool exact = k < dependencies.size() && !dependencies[k].empty(); + if (exact) + { + for (const uint32_t vr : dependencies[k]) + { + const auto it = columnOfValueReference[i].find(vr); + if (it == columnOfValueReference[i].end()) + { + // A value reference that is not one of the unknowns: a parameter or + // time, which contribute nothing, or a variable of the component + // that stands for the unknowns behind it. Cannot tell the two + // apart here, so take the whole component. + exact = false; + break; + } + columnsByRow[row].push_back(it->second); + } + } + if (exact) + incomingLoopColumns(i, columnsByRow[row]); + else + { + columnsByRow[row].clear(); + columnsAround(i, columnsByRow[row]); + } + } + } + else + { + // x' - f(x, u, t): its own states, and whatever drives its inputs. + for (size_t k = 0; k < nStates[i]; ++k, ++row) + columnsAround(i, columnsByRow[row]); + } + } + + // A loop row is output - input: the unknown it stands for, and whatever the + // component computing that output depends on. + for (size_t k = 0; k < liftedConnections.size(); ++k, ++row) + { + columnsByRow[row].push_back(static_cast(loopBase + k)); + if (componentOfLoopOutput[k] >= 0) + columnsAround(static_cast(componentOfLoopOutput[k]), columnsByRow[row]); + else + for (size_t c = 0; c < n; ++c) + columnsByRow[row].push_back(static_cast(c)); + } + + // Transpose, and drop the duplicates a conservative row may have collected. + daeJacRowsByColumn.assign(n, {}); + size_t nnz = 0; + for (size_t r = 0; r < n; ++r) + { + std::sort(columnsByRow[r].begin(), columnsByRow[r].end()); + columnsByRow[r].erase(std::unique(columnsByRow[r].begin(), columnsByRow[r].end()), columnsByRow[r].end()); + for (const int c : columnsByRow[r]) + daeJacRowsByColumn[c].push_back(static_cast(r)); + nnz += columnsByRow[r].size(); + } + + // Greedy distance-1 colouring of the columns: two columns may share a colour + // only when no row reaches both, which is what makes their contributions + // separable in a single differenced evaluation. + daeJacColorOfColumn.assign(n, -1); + for (size_t c = 0; c < n; ++c) + { + std::set forbidden; + for (const int r : daeJacRowsByColumn[c]) + for (const int other : columnsByRow[r]) + if (other != static_cast(c) && daeJacColorOfColumn[other] >= 0) + forbidden.insert(daeJacColorOfColumn[other]); + + int color = 0; + while (forbidden.count(color)) + ++color; + daeJacColorOfColumn[c] = color; + if (static_cast(color) >= daeJacColumnsOfColor.size()) + daeJacColumnsOfColor.resize(color + 1); + daeJacColumnsOfColor[color].push_back(static_cast(c)); + } + + daeJacIncrement.assign(n, 0.0); + logInfo("system \"" + std::string(getFullCref()) + "\": DAE Jacobian " + std::to_string(n) + "x" + std::to_string(n) + + ", " + std::to_string(nnz) + " structural nonzeros, " + std::to_string(daeJacColumnsOfColor.size()) + + " colour(s) — sparse (KLU), one residual evaluation per colour instead of one per unknown"); +} + +/** + * \brief The DAE Jacobian dF/dy + cj*dF/dy', differenced one colour at a time. + * + * The colouring makes every column of a colour distinguishable in a single + * evaluation, because no row reaches two of them, so the whole matrix costs as + * many residual evaluations as there are colours. IDA's own difference quotient + * would cost one per unknown, and each of those calls into every FMU of the + * system. + */ +int oms::ida_jac(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector rr, + SUNMatrix Jac, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) +{ + SystemSC* system = (SystemSC*)user_data; + const size_t n = system->nDaeUnknowns; + + // A sparse matrix carries its structure in the same arrays as its values, and + // SUNMatZero would wipe it; every entry of the pattern is written below + // anyway, structure and value together. + sunindextype* colptrs = nullptr; + sunindextype* rowvals = nullptr; + sunrealtype* data = nullptr; + if (system->daeJacSparse) + { + colptrs = SUNSparseMatrix_IndexPointers(Jac); + rowvals = SUNSparseMatrix_IndexValues(Jac); + data = SUNSparseMatrix_Data(Jac); + sunindextype at = 0; + for (size_t j = 0; j < n; ++j) + { + colptrs[j] = at; + at += static_cast(system->daeJacRowsByColumn[j].size()); + } + colptrs[n] = at; + } + else + SUNMatZero(Jac); + + sunrealtype hh = 0.0; + if (IDAGetCurrentStep(system->solverData.ida.mem, &hh) != IDA_SUCCESS) + return -1; + // Its own vector: tmp3 takes the perturbed residual below, which would + // otherwise overwrite the weights before the next colour reads them. + if (IDAGetErrWeights(system->solverData.ida.mem, system->solverData.ida.ewt) != IDA_SUCCESS) + return -1; + + const sunrealtype srur = std::sqrt(SUN_UNIT_ROUNDOFF); + + for (const std::vector& columns : system->daeJacColumnsOfColor) + { + N_VScale(1.0, yy, tmp1); + N_VScale(1.0, yp, tmp2); + + for (const int j : columns) + { + const sunrealtype yj = NV_Ith_S(yy, j); + const sunrealtype ypj = NV_Ith_S(yp, j); + const sunrealtype ewtj = NV_Ith_S(system->solverData.ida.ewt, j); + + // IDA's own difference quotient increment (IDADenseDQJac). + sunrealtype inc = std::max(srur * std::max(std::fabs(yj), std::fabs(hh * ypj)), 1.0 / ewtj); + if (hh * ypj < 0.0) + inc = -inc; + inc = (yj + inc) - yj; + if (0.0 == inc) + inc = srur; + + system->daeJacIncrement[j] = inc; + NV_Ith_S(tmp1, j) = yj + inc; + NV_Ith_S(tmp2, j) = ypj + cj * inc; + } + + if (0 != ida_res(t, tmp1, tmp2, tmp3, user_data)) + return -1; + + for (const int j : columns) + { + const sunrealtype scale = 1.0 / system->daeJacIncrement[j]; + const std::vector& rows = system->daeJacRowsByColumn[j]; + for (size_t k = 0; k < rows.size(); ++k) + { + const int i = rows[k]; + const sunrealtype value = (NV_Ith_S(tmp3, i) - NV_Ith_S(rr, i)) * scale; + if (system->daeJacSparse) + { + const sunindextype at = colptrs[j] + static_cast(k); + rowvals[at] = i; + data[at] = value; + } + else + SM_ELEMENT_D(Jac, i, j) = value; + } + } + } + + // The point the components hold is the perturbed one; put the original back so + // whatever reads them next does not see it. + if (oms_status_ok != system->setDaePoint(t, yy, yp)) + return -1; + + return 0; +} + +oms_status_enu_t oms::SystemSC::setDaePoint(double t, N_Vector yy, N_Vector yp) +{ + oms_status_enu_t status; + + for (size_t i = 0, j = 0; i < fmus.size(); ++i) + { + fmus[i]->setTime(t); + + const size_t xoff = j; // where this component's states sit in y and yp + for (size_t k = 0; k < nStates[i]; ++k, ++j) + states[i][k] = NV_Ith_S(yy, j); + for (size_t k = 0; k < nAlgebraic[i]; ++k, ++j) + algebraicVars[i][k] = NV_Ith_S(yy, j); + + if (nStates[i] > 0) + { + status = fmus[i]->setContinuousStates(states[i]); + if (oms_status_ok != status) return status; + } + + if (daeFmus[i]) + { + // The implicit form's residuals are a function of (t, x, x', z), so the + // derivatives are knowns the master hands over, not something the FMU owes. + if (nStates[i] > 0) + { + for (size_t k = 0; k < nStates[i]; ++k) + states_der[i][k] = NV_Ith_S(yp, xoff + k); + status = daeFmus[i]->setDerivatives(states_der[i]); + if (oms_status_ok != status) return status; + } + if (nAlgebraic[i] > 0) + { + status = daeFmus[i]->setAlgebraicVariables(algebraicVars[i]); + if (oms_status_ok != status) return status; + } + } + } + + // The loop connections IDA owns: the unknown is the value of the driven input. + // Set before the propagation below, which then leaves those loops alone. + for (size_t k = 0; k < liftedConnections.size(); ++k) + { + const double value = NV_Ith_S(yy, nDaeUnknowns - liftedConnections.size() + k); + if (oms_status_ok != setReal(simulationGraph.getNodes()[liftedConnections[k].input].getName(), value)) + return oms_status_error; + } + + // The rest of the coupling, before the rows below are asked for. + return updateInputsInternal(simulationGraph, &liftedLoops); +} + +oms_status_enu_t oms::SystemSC::getDaePoint(N_Vector yy, N_Vector yp) +{ + oms_status_enu_t status; + + for (size_t i = 0, j = 0; i < fmus.size(); ++i) + { + if (nStates[i] > 0) + { + status = fmus[i]->getContinuousStates(states[i]); + if (oms_status_ok != status) return status; + status = fmus[i]->getDerivatives(states_der[i]); + if (oms_status_ok != status) return status; + } + if (daeFmus[i] && nAlgebraic[i] > 0) + { + status = daeFmus[i]->getAlgebraicVariables(algebraicVars[i]); + if (oms_status_ok != status) return status; + } + + for (size_t k = 0; k < nStates[i]; ++k, ++j) + { + NV_Ith_S(yy, j) = states[i][k]; + NV_Ith_S(yp, j) = states_der[i][k]; + } + for (size_t k = 0; k < nAlgebraic[i]; ++k, ++j) + { + NV_Ith_S(yy, j) = algebraicVars[i][k]; + NV_Ith_S(yp, j) = 0.0; // an algebraic unknown has no derivative row + } + } + + // The lifted loop connections, seeded with the values the inputs hold — after + // initialization, or after an event, that is the best guess there is. + for (size_t k = 0, j = nDaeUnknowns - liftedConnections.size(); k < liftedConnections.size(); ++k, ++j) + { + double value = 0.0; + if (oms_status_ok != getReal(simulationGraph.getNodes()[liftedConnections[k].input].getName(), value)) + return oms_status_error; + NV_Ith_S(yy, j) = value; + NV_Ith_S(yp, j) = 0.0; + } + + return oms_status_ok; +} + +int oms::ida_res(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data) +{ + SystemSC* system = (SystemSC*)user_data; + + if (oms_status_ok != system->setDaePoint(t, yy, yp)) + return -1; + + for (size_t i = 0, j = 0; i < system->fmus.size(); ++i) + { + if (system->daeFmus[i]) + { + const size_t n = system->nStates[i] + system->nAlgebraic[i]; + if (n == 0) + continue; + if (oms_status_ok != system->daeFmus[i]->getDaeResiduals(system->residuals[i])) + return -1; + for (size_t k = 0; k < n; ++k, ++j) + NV_Ith_S(rr, j) = system->residuals[i][k]; + } + else + { + // An ODE component closes its own rows: x' - f(x, u, t). + if (0 == system->nStates[i]) + continue; + if (oms_status_ok != system->fmus[i]->getDerivatives(system->states_der[i])) + return -1; + for (size_t k = 0; k < system->nStates[i]; ++k, ++j) + NV_Ith_S(rr, j) = NV_Ith_S(yp, j) - system->states_der[i][k]; + } + } + + // The lifted loop rows: what the driving output now says, against the value + // IDA is carrying for the input it drives. + for (size_t k = 0, j = system->nDaeUnknowns - system->liftedConnections.size(); + k < system->liftedConnections.size(); ++k, ++j) + { + double value = 0.0; + if (oms_status_ok != system->getReal(system->simulationGraph.getNodes()[system->liftedConnections[k].output].getName(), value)) + return -1; + NV_Ith_S(rr, j) = value - NV_Ith_S(yy, j); + } + + return 0; +} + +int oms::ida_roots(sunrealtype t, N_Vector yy, N_Vector yp, sunrealtype *gout, void* user_data) +{ + SystemSC* system = (SystemSC*)user_data; + + if (oms_status_ok != system->setDaePoint(t, yy, yp)) + return -1; + + for (size_t i = 0, j_gout = 0; i < system->fmus.size(); ++i) + { + if (0 == system->nEventIndicators[i]) + continue; + if (oms_status_ok != system->fmus[i]->getEventindicators(system->event_indicators[i], system->nEventIndicators[i])) + return -1; + for (size_t k = 0; k < system->nEventIndicators[i]; ++k, ++j_gout) + gout[j_gout] = system->event_indicators[i][k]; + } + + return 0; +} + oms::SystemSC::SystemSC(const ComRef& cref, Model* parentModel, System* parentSystem) : oms::System(cref, oms_system_sc, parentModel, parentSystem, oms_solver_sc_cvode) { @@ -187,6 +740,8 @@ std::string oms::SystemSC::getSolverName() const return std::string("euler"); case oms_solver_sc_cvode: return std::string("cvode"); + case oms_solver_sc_ida: + return std::string("ida"); default: return std::string("unknown"); } @@ -198,6 +753,8 @@ oms_status_enu_t oms::SystemSC::setSolverMethod(std::string solver) solverMethod = oms_solver_sc_explicit_euler; else if (std::string("cvode") == solver) solverMethod = oms_solver_sc_cvode; + else if (std::string("ida") == solver) + solverMethod = oms_solver_sc_ida; else return oms_status_error; @@ -269,6 +826,8 @@ oms_status_enu_t oms::SystemSC::instantiate() ; else if (oms_solver_sc_cvode == solverMethod) solverData.cvode.mem = nullptr; + else if (oms_solver_sc_ida == solverMethod) + solverData.ida.mem = nullptr; else return logError_InternalError; @@ -283,7 +842,18 @@ oms_status_enu_t oms::SystemSC::initialize() if (oms_status_ok != updateDependencyGraphs()) return oms_status_error; - if (oms_status_ok != updateInputs(initializationGraph)) + // fmi-ls-dae: a component in DAE mode does not solve for its own algebraic + // outputs — the master does — so a loop through one cannot be closed by the + // inner Newton here. Leave every loop of the initialization graph to the + // consistent-initial-conditions solve at the end of this function. + daeMode = anyComponentInDaeMode(); + if (daeMode) + { + std::vector allLoops = allLoopsOf(initializationGraph); + if (oms_status_ok != updateInputsInternal(initializationGraph, &allLoops)) + return oms_status_error; + } + else if (oms_status_ok != updateInputs(initializationGraph)) return oms_status_error; for (const auto& subsystem : getSubSystems()) @@ -316,8 +886,41 @@ oms_status_enu_t oms::SystemSC::initialize() states_nominal.push_back((double*)calloc(nStates.back(), sizeof(double))); event_indicators.push_back((double*)calloc(nEventIndicators.back(), sizeof(double))); event_indicators_prev.push_back((double*)calloc(nEventIndicators.back(), sizeof(double))); + + // fmi-ls-dae: a component that was switched into DAE mode carries algebraic + // unknowns of its own and answers with residuals instead of derivatives. + ComponentFMU3ME* fmu3 = dynamic_cast(fmu); + const bool inDaeMode = fmu3 && fmu3->isInDaeMode(); + daeFmus.push_back(inDaeMode ? fmu3 : nullptr); + nAlgebraic.push_back(inDaeMode ? fmu3->getNumberOfAlgebraicVariables() : 0); + algebraicVars.push_back((double*)calloc(nAlgebraic.back(), sizeof(double))); + residuals.push_back((double*)calloc(inDaeMode ? fmu3->getNumberOfDaeResiduals() : 0, sizeof(double))); + if (inDaeMode) + daeMode = true; // already true from anyComponentInDaeMode(), kept for clarity + nDaeUnknowns += nStates.back() + nAlgebraic.back(); } + // fmi-ls-dae: with IDA driving the system, an algebraic loop over connections + // can be rows of it instead of an inner Newton of its own. + if (daeMode) + { + collectLiftedLoops(); + nDaeUnknowns += liftedConnections.size(); + } + + if (daeMode && oms_solver_sc_ida != solverMethod) + { + // DAE mode is not a choice the user makes: a component that declares a DAE + // formulation has no explicit ODE for CVODE to integrate. + logInfo("system \"" + std::string(getFullCref()) + "\" contains a component in DAE mode (fmi-ls-dae); integrating with IDA instead of " + getSolverName()); + solverMethod = oms_solver_sc_ida; + } + else if (!daeMode && oms_solver_sc_ida == solverMethod) + { + logInfo("system \"" + std::string(getFullCref()) + "\" contains no component in DAE mode; integrating with CVODE instead of ida"); + solverMethod = oms_solver_sc_cvode; + } + // Now that fmus is filled, the per-FMU flags the integrator steps write to // can be sized. make_unique value-initializes them, i.e. all false. callEventUpdate = std::make_unique(fmus.size()); @@ -472,6 +1075,135 @@ oms_status_enu_t oms::SystemSC::initialize() flag = CVodeSetMaxNumSteps(solverData.cvode.mem, Flags::CVODEMaxSteps()); // MAXIMUM NUMBER OF STEPS if (flag < 0) logError("SUNDIALS_ERROR: CVodeSetMaxNumSteps() failed with flag = " + std::to_string(flag)); } + else if (oms_solver_sc_ida == solverMethod) + { + // fmi-ls-dae: one implicit system over [states | algebraic variables]. + if (SUNContext_Create(SUN_COMM_NULL, &solverData.ida.sunctx) != SUN_SUCCESS) + logError("SUNDIALS_ERROR: SUNContext_Create() failed"); + + /* Mute SUNDIALS' own output, use OMSimulator's logger */ + { + SUNLogger logger = NULL; + if (SUNContext_GetLogger(solverData.ida.sunctx, &logger) == SUN_SUCCESS && logger != NULL) + { + SUNLogger_SetErrorFilename(logger, ""); + SUNLogger_SetWarningFilename(logger, ""); + SUNLogger_SetInfoFilename(logger, ""); + SUNLogger_SetDebugFilename(logger, ""); + } + } + + const long n = static_cast(nDaeUnknowns); + solverData.ida.y = N_VNew_Serial(n, solverData.ida.sunctx); + solverData.ida.yp = N_VNew_Serial(n, solverData.ida.sunctx); + solverData.ida.id = N_VNew_Serial(n, solverData.ida.sunctx); + solverData.ida.abstol = N_VNew_Serial(n, solverData.ida.sunctx); + solverData.ida.ewt = N_VNew_Serial(n, solverData.ida.sunctx); + if (!solverData.ida.y || !solverData.ida.yp || !solverData.ida.id || !solverData.ida.abstol || !solverData.ida.ewt) + return logError("SUNDIALS_ERROR: N_VNew_Serial() failed - returned NULL pointer"); + + // The point the components hold after their own initialization, and which + // of the unknowns are differential: IDACalcIC solves for the algebraic ones + // and for the derivatives of the differential ones. + if (oms_status_ok != getDaePoint(solverData.ida.y, solverData.ida.yp)) + return oms_status_error; + + for (size_t i = 0, j = 0; i < fmus.size(); ++i) + { + for (size_t k = 0; k < nStates[i]; ++k, ++j) + { + NV_Ith_S(solverData.ida.id, j) = 1.0; + NV_Ith_S(solverData.ida.abstol, j) = relativeTolerance * states_nominal[i][k]; + } + for (size_t k = 0; k < nAlgebraic[i]; ++k, ++j) + { + NV_Ith_S(solverData.ida.id, j) = 0.0; + // The FMU gives no nominal for an algebraic unknown of its own, so the + // relative tolerance stands in for it. + NV_Ith_S(solverData.ida.abstol, j) = relativeTolerance; + } + } + + for (size_t k = 0, j = nDaeUnknowns - liftedConnections.size(); k < liftedConnections.size(); ++k, ++j) + { + NV_Ith_S(solverData.ida.id, j) = 0.0; + NV_Ith_S(solverData.ida.abstol, j) = relativeTolerance; + } + + solverData.ida.mem = IDACreate(solverData.ida.sunctx); + if (!solverData.ida.mem) return logError("SUNDIALS_ERROR: IDACreate() failed - returned NULL pointer"); + + int flag = IDASetUserData(solverData.ida.mem, (void*)this); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASetUserData() failed with flag = " + std::to_string(flag)); + + flag = IDAInit(solverData.ida.mem, ida_res, time, solverData.ida.y, solverData.ida.yp); + if (flag < 0) return logError("SUNDIALS_ERROR: IDAInit() failed with flag = " + std::to_string(flag)); + + flag = IDASVtolerances(solverData.ida.mem, relativeTolerance, solverData.ida.abstol); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASVtolerances() failed with flag = " + std::to_string(flag)); + + flag = IDASetId(solverData.ida.mem, solverData.ida.id); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASetId() failed with flag = " + std::to_string(flag)); + + if (n_event_indicators > 0) + { + flag = IDARootInit(solverData.ida.mem, static_cast(n_event_indicators), ida_roots); + if (flag < 0) return logError("SUNDIALS_ERROR: IDARootInit() failed with flag = " + std::to_string(flag)); + } + + // The components' manifests say which unknowns each residual reaches. That + // structure does two things: it colours the columns, so the Jacobian is + // differenced a colour at a time instead of a column at a time — a handful + // of residual evaluations instead of one per unknown, each of which calls + // into every FMU — and it lets the matrix be held sparse and factorized by + // KLU. Without a usable structure the matrix is dense and IDA differences it + // itself. + buildDaeJacobianSparsity(); + daeJacSparse = !daeJacColumnsOfColor.empty(); + + if (daeJacSparse) + { + sunindextype nnz = 0; + for (const std::vector& rows : daeJacRowsByColumn) + nnz += static_cast(rows.size()); + + solverData.ida.J = SUNSparseMatrix(n, n, nnz, CSC_MAT, solverData.ida.sunctx); + if (!solverData.ida.J) return logError("SUNDIALS_ERROR: SUNSparseMatrix() failed"); + solverData.ida.linSol = SUNLinSol_KLU(solverData.ida.y, solverData.ida.J, solverData.ida.sunctx); + if (!solverData.ida.linSol) return logError("SUNDIALS_ERROR: SUNLinSol_KLU() failed"); + } + else + { + solverData.ida.J = SUNDenseMatrix(n, n, solverData.ida.sunctx); + if (!solverData.ida.J) return logError("SUNDIALS_ERROR: SUNDenseMatrix() failed"); + solverData.ida.linSol = SUNLinSol_Dense(solverData.ida.y, solverData.ida.J, solverData.ida.sunctx); + if (!solverData.ida.linSol) return logError("SUNDIALS_ERROR: SUNLinSol_Dense() failed"); + } + + flag = IDASetLinearSolver(solverData.ida.mem, solverData.ida.linSol, solverData.ida.J); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASetLinearSolver() failed with flag = " + std::to_string(flag)); + + if (daeJacSparse) + { + flag = IDASetJacFn(solverData.ida.mem, ida_jac); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASetJacFn() failed with flag = " + std::to_string(flag)); + } + + flag = IDASetMaxStep(solverData.ida.mem, maximumStepSize); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASetMaxStep() failed with flag = " + std::to_string(flag)); + if (initialStepSize > 0.0) + { + flag = IDASetInitStep(solverData.ida.mem, initialStepSize); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASetInitStep() failed with flag = " + std::to_string(flag)); + } + flag = IDASetMaxNumSteps(solverData.ida.mem, Flags::CVODEMaxSteps()); + if (flag < 0) return logError("SUNDIALS_ERROR: IDASetMaxNumSteps() failed with flag = " + std::to_string(flag)); + + // The start point the components initialized to need not satisfy the + // coupled system; IDA solves for a consistent one before the first step. + if (oms_status_ok != calcConsistentInitialConditions(time + maximumStepSize)) + return oms_status_error; + } // Mark algebraic loops to be updated on next call forceLoopsToBeUpdated(); @@ -489,6 +1221,33 @@ oms_status_enu_t oms::SystemSC::terminate() if (oms_status_ok != component.second->terminate()) return oms_status_error; + if (oms_solver_sc_ida == solverMethod && solverData.ida.mem) + { + long int nst = 0, nre = 0, nsetups = 0, nni = 0, ncfn = 0, netf = 0; + IDAGetNumSteps(solverData.ida.mem, &nst); + IDAGetNumResEvals(solverData.ida.mem, &nre); + IDAGetNumLinSolvSetups(solverData.ida.mem, &nsetups); + IDAGetNumErrTestFails(solverData.ida.mem, &netf); + IDAGetNumNonlinSolvIters(solverData.ida.mem, &nni); + IDAGetNumNonlinSolvConvFails(solverData.ida.mem, &ncfn); + + std::string msg = "Final Statistics for '" + std::string(getFullCref()) + "':\n"; + msg += "NumSteps = " + std::to_string(nst) + " NumResEvals = " + std::to_string(nre) + " NumLinSolvSetups = " + std::to_string(nsetups) + "\n"; + msg += "NumNonlinSolvIters = " + std::to_string(nni) + " NumNonlinSolvConvFails = " + std::to_string(ncfn) + " NumErrTestFails = " + std::to_string(netf); + logInfo(msg); + + SUNMatDestroy(solverData.ida.J); + SUNLinSolFree(solverData.ida.linSol); + N_VDestroy_Serial(solverData.ida.y); + N_VDestroy_Serial(solverData.ida.yp); + N_VDestroy_Serial(solverData.ida.id); + N_VDestroy_Serial(solverData.ida.abstol); + N_VDestroy_Serial(solverData.ida.ewt); + IDAFree(&(solverData.ida.mem)); + SUNContext_Free(&(solverData.ida.sunctx)); + solverData.ida.mem = nullptr; + } + if (oms_solver_sc_cvode == solverMethod && solverData.cvode.mem) { long int nst, nfe, nsetups, nni, ncfn, netf; @@ -529,13 +1288,28 @@ oms_status_enu_t oms::SystemSC::terminate() free(states_nominal[i]); free(event_indicators[i]); free(event_indicators_prev[i]); + free(algebraicVars[i]); + free(residuals[i]); } fmus.clear(); + daeFmus.clear(); callEventUpdate.reset(); terminateSimulation.reset(); nStates.clear(); nEventIndicators.clear(); + nAlgebraic.clear(); + algebraicVars.clear(); + residuals.clear(); + daeMode = false; + nDaeUnknowns = 0; + liftedConnections.clear(); + liftedLoops.clear(); + daeJacRowsByColumn.clear(); + daeJacColorOfColumn.clear(); + daeJacColumnsOfColor.clear(); + daeJacIncrement.clear(); + daeJacSparse = false; states.clear(); states_der.clear(); states_nominal.clear(); @@ -557,6 +1331,33 @@ oms_status_enu_t oms::SystemSC::reset() time = getModel().getStartTime(); + if (oms_solver_sc_ida == solverMethod && solverData.ida.mem) + { + long int nst = 0, nre = 0, nsetups = 0, nni = 0, ncfn = 0, netf = 0; + IDAGetNumSteps(solverData.ida.mem, &nst); + IDAGetNumResEvals(solverData.ida.mem, &nre); + IDAGetNumLinSolvSetups(solverData.ida.mem, &nsetups); + IDAGetNumErrTestFails(solverData.ida.mem, &netf); + IDAGetNumNonlinSolvIters(solverData.ida.mem, &nni); + IDAGetNumNonlinSolvConvFails(solverData.ida.mem, &ncfn); + + std::string msg = "Final Statistics for '" + std::string(getFullCref()) + "':\n"; + msg += "NumSteps = " + std::to_string(nst) + " NumResEvals = " + std::to_string(nre) + " NumLinSolvSetups = " + std::to_string(nsetups) + "\n"; + msg += "NumNonlinSolvIters = " + std::to_string(nni) + " NumNonlinSolvConvFails = " + std::to_string(ncfn) + " NumErrTestFails = " + std::to_string(netf); + logInfo(msg); + + SUNMatDestroy(solverData.ida.J); + SUNLinSolFree(solverData.ida.linSol); + N_VDestroy_Serial(solverData.ida.y); + N_VDestroy_Serial(solverData.ida.yp); + N_VDestroy_Serial(solverData.ida.id); + N_VDestroy_Serial(solverData.ida.abstol); + N_VDestroy_Serial(solverData.ida.ewt); + IDAFree(&(solverData.ida.mem)); + SUNContext_Free(&(solverData.ida.sunctx)); + solverData.ida.mem = nullptr; + } + if (oms_solver_sc_cvode == solverMethod && solverData.cvode.mem) { long int nst, nfe, nsetups, nni, ncfn, netf; @@ -595,16 +1396,31 @@ oms_status_enu_t oms::SystemSC::reset() // instead of appending duplicate entries for every FMU (which corrupts // solver setup and causes lifecycle calls like enterContinuousTimeMode() // to be issued twice for the same FMU). + for (double* ptr : algebraicVars) free(ptr); + for (double* ptr : residuals) free(ptr); for (double* ptr : states) free(ptr); for (double* ptr : states_der) free(ptr); for (double* ptr : states_nominal) free(ptr); for (double* ptr : event_indicators) free(ptr); for (double* ptr : event_indicators_prev) free(ptr); fmus.clear(); + daeFmus.clear(); callEventUpdate.reset(); terminateSimulation.reset(); nStates.clear(); nEventIndicators.clear(); + nAlgebraic.clear(); + algebraicVars.clear(); + residuals.clear(); + daeMode = false; + nDaeUnknowns = 0; + liftedConnections.clear(); + liftedLoops.clear(); + daeJacRowsByColumn.clear(); + daeJacColorOfColumn.clear(); + daeJacColumnsOfColor.clear(); + daeJacIncrement.clear(); + daeJacSparse = false; states.clear(); states_der.clear(); states_nominal.clear(); @@ -624,6 +1440,9 @@ oms_status_enu_t oms::SystemSC::doStep() case oms_solver_sc_cvode: return doStepCVODE(); + case oms_solver_sc_ida: + return doStepIDA(); + default: return logError_InternalError; } @@ -1026,6 +1845,165 @@ oms_status_enu_t oms::SystemSC::doStepCVODE() } +/** + * \brief IDACalcIC over the point the components hold, then read the consistent + * point back into them. + * + * The components initialize on their own, each consistent with itself but not + * necessarily with the coupled system; and after an event the same is true + * again. IDA_YA_YDP_INIT solves for the algebraic unknowns and the derivatives + * of the differential ones, keeping the states where they are. + */ +oms_status_enu_t oms::SystemSC::calcConsistentInitialConditions(double tout) +{ + // IDACalcIC wants a time it may integrate towards; anything but the current + // one will do, and the step it takes is thrown away. + if (tout <= time) + tout = time + 1.0; + + int flag = IDACalcIC(solverData.ida.mem, IDA_YA_YDP_INIT, tout); + if (flag < 0) + return logError("SUNDIALS_ERROR: IDACalcIC() failed with flag = " + std::to_string(flag) + + " at time " + std::to_string(time) + + " (the coupled DAE has no consistent point there, or its Jacobian is singular)"); + + flag = IDAGetConsistentIC(solverData.ida.mem, solverData.ida.y, solverData.ida.yp); + if (flag < 0) + return logError("SUNDIALS_ERROR: IDAGetConsistentIC() failed with flag = " + std::to_string(flag)); + + // Leave the components standing at the consistent point, so what is emitted + // and what the next step starts from agree with it. + return setDaePoint(time, solverData.ida.y, solverData.ida.yp); +} + +oms_status_enu_t oms::SystemSC::doStepIDA() +{ + oms_status_enu_t status; + int flag; + + double end_time = std::min(time + maximumStepSize, getModel().getStopTime()); + + // find next time event + double tnext = end_time + 1.0; + for (size_t i = 0; i < fmus.size(); ++i) + { + if (fmus[i]->getNextEventTimeDefined() && (tnext > fmus[i]->getNextEventTime())) + tnext = fmus[i]->getNextEventTime(); + + if (fmus[i]->getTerminateSimulation()) + { + logInfo("Simulation terminated by FMU " + std::string(fmus[i]->getFullCref()) + " at time " + std::to_string(time)); + getModel().setStopTime(time); + time = end_time; + } + } + snapToEventTime(end_time, tnext); + + while (time < end_time) + { + const double tout = std::min(tnext, end_time); + + // Nothing left to integrate in an interval IDA cannot tell apart from zero. + if (tout - time < 2.0 * SUN_UNIT_ROUNDOFF * std::max(std::fabs(time), std::fabs(tout))) + { + logDebug("IDA: skipping degenerate interval " + std::to_string(time) + " -> " + std::to_string(tout)); + time = end_time; + break; + } + + logDebug("IDA: " + std::to_string(time) + " -> " + std::to_string(tout)); + flag = IDASolve(solverData.ida.mem, tout, &time, solverData.ida.y, solverData.ida.yp, IDA_NORMAL); + if (flag < 0) + return logError("IDA failed with flag = " + std::to_string(flag) + " at time " + std::to_string(time)); + + // The components stand at the point IDA returned. + status = setDaePoint(time, solverData.ida.y, solverData.ida.yp); + if (oms_status_ok != status) return status; + + for (const auto& component : getComponents()) + component.second->setTime(time); + + for (size_t i = 0; i < fmus.size(); ++i) + { + status = fmus[i]->completedIntegratorStep(true, callEventUpdate[i], terminateSimulation[i]); + if (oms_status_ok != status) return status; + } + + if (IDA_ROOT_RETURN == flag || time == tnext) + { + logDebug("event found!!! " + std::to_string(time)); + + // emit the left limit of the event + if (isTopLevelSystem()) + getModel().emit(time, false); + + for (size_t i = 0; i < fmus.size(); ++i) + { + status = fmus[i]->enterEventMode(); + if (oms_status_ok != status) return status; + + fmus[i]->doEventIteration(); + } + + { + // As during initialization: the loops are the integrator's, and + // calcConsistentInitialConditions() below closes them. + std::vector allLoops = allLoopsOf(eventGraph); + updateInputsInternal(eventGraph, &allLoops); + } + + for (size_t i = 0; i < fmus.size(); ++i) + { + status = fmus[i]->enterContinuousTimeMode(); + if (oms_status_ok != status) return status; + } + + // find next time event + tnext = end_time + 1.0; + for (size_t i = 0; i < fmus.size(); ++i) + { + if (fmus[i]->getNextEventTimeDefined() && (tnext > fmus[i]->getNextEventTime())) + tnext = fmus[i]->getNextEventTime(); + + if (fmus[i]->getTerminateSimulation()) + { + logInfo("Simulation terminated by FMU " + std::string(fmus[i]->getFullCref()) + " at time " + std::to_string(time)); + getModel().setStopTime(time); + time = end_time; + } + } + + // The event moved the components; take the point they hold now, restart + // IDA there and let it find a consistent one again — the algebraic + // unknowns and the derivatives need not have survived the event. + status = getDaePoint(solverData.ida.y, solverData.ida.yp); + if (oms_status_ok != status) return status; + + flag = IDAReInit(solverData.ida.mem, time, solverData.ida.y, solverData.ida.yp); + if (flag < 0) return logError("SUNDIALS_ERROR: IDAReInit() failed with flag = " + std::to_string(flag)); + + status = calcConsistentInitialConditions(std::min(tnext, end_time)); + if (oms_status_ok != status) return status; + + // emit the right limit of the event + { + std::vector allLoops = allLoopsOf(eventGraph); + updateInputsInternal(eventGraph, &allLoops); + } + if (isTopLevelSystem()) + getModel().emit(time, true); + + continue; + } + + updateInputs(simulationGraph); + if (isTopLevelSystem()) + getModel().emit(time, false); + } + + return oms_status_ok; +} + oms_status_enu_t oms::SystemSC::stepUntil(double stopTime) { CallClock callClock(clock); @@ -1053,6 +2031,11 @@ oms_status_enu_t oms::SystemSC::stepUntil(double stopTime) } oms_status_enu_t oms::SystemSC::updateInputs(DirectedGraph& graph) +{ + return updateInputsInternal(graph, nullptr); +} + +oms_status_enu_t oms::SystemSC::updateInputsInternal(DirectedGraph& graph, const std::vector* liftedLoops) { CallClock callClock(clock); oms_status_enu_t status; @@ -1104,11 +2087,17 @@ oms_status_enu_t oms::SystemSC::updateInputs(DirectedGraph& graph) } else { - status = solveAlgLoop(graph, loopNum); - if (oms_status_ok != status) + // A loop IDA owns is already satisfied by the point it set; solving it + // again here would fight the integrator for the same unknowns. + const bool lifted = liftedLoops && loopNum < static_cast(liftedLoops->size()) && (*liftedLoops)[loopNum]; + if (!lifted) { - forceLoopsToBeUpdated(); - return status; + status = solveAlgLoop(graph, loopNum); + if (oms_status_ok != status) + { + forceLoopsToBeUpdated(); + return status; + } } loopNum++; } diff --git a/src/OMSimulatorLib/SystemSC.h b/src/OMSimulatorLib/SystemSC.h index b1246a299..af689da2b 100644 --- a/src/OMSimulatorLib/SystemSC.h +++ b/src/OMSimulatorLib/SystemSC.h @@ -46,16 +46,24 @@ #include /* SUNLogger */ #include /* prototypes for CVODE fcts., consts. */ #include /* serial N_Vector types, fcts., macros */ +#include /* prototypes for IDA fcts., consts. */ #include /* Default dense linear solver */ +#include /* KLU, for the structured DAE Jacobian */ +#include namespace oms { class Model; class ComponentFMUME; + class ComponentFMU3ME; class Component; int cvode_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); int cvode_rhs_algebraic(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); int cvode_roots(sunrealtype t, N_Vector y, sunrealtype *gout, void* user_data); + int ida_res(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data); + int ida_roots(sunrealtype t, N_Vector yy, N_Vector yp, sunrealtype *gout, void* user_data); + int ida_jac(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector rr, + SUNMatrix Jac, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); class SystemSC : public System { @@ -83,6 +91,41 @@ namespace oms private: oms_status_enu_t doStepEuler(); oms_status_enu_t doStepCVODE(); + oms_status_enu_t doStepIDA(); + /** + * \brief updateInputs, with the algebraic loops IDA owns left alone. + * + * \param liftedLoops one flag per loop of the graph, in the order the + * sorted connections hold them; a loop marked true is + * skipped, because its connections are rows of the + * global system and IDA has already set their inputs. + * nullptr solves every loop, which is what updateInputs + * itself does. + */ + oms_status_enu_t updateInputsInternal(DirectedGraph& graph, const std::vector* liftedLoops); + /// Which loops of the simulation graph can become rows: every connection real. + void collectLiftedLoops(); + /// A flag per loop of the graph, all true: leave every loop to the integrator. + std::vector allLoopsOf(DirectedGraph& graph) const; + /// Whether any component was switched into DAE mode; readable before initialize(). + bool anyComponentInDaeMode(); + /** + * \brief The structure of the DAE Jacobian, out of the components' manifests. + * + * Which unknowns each row reaches, and a colouring of the columns from it: two + * columns share a colour when no row reaches both, so one residual evaluation + * differences all of them at once. A dense difference quotient costs one + * residual evaluation per unknown — and every one of those calls into every + * FMU — where a coloured one costs as many as there are colours, which for a + * coupled system is a handful however many unknowns there are. + */ + void buildDaeJacobianSparsity(); + /// Distribute the DAE unknowns y (and the state derivatives yp) into the FMUs. + oms_status_enu_t setDaePoint(double t, N_Vector yy, N_Vector yp); + /// Read the point the FMUs hold back into y and yp, after a solve or an event. + oms_status_enu_t getDaePoint(N_Vector yy, N_Vector yp); + /// IDACalcIC over the current point, then read the consistent point back. + oms_status_enu_t calcConsistentInitialConditions(double tout); protected: SystemSC(const ComRef& cref, Model* parentModel, System* parentSystem); @@ -101,6 +144,59 @@ namespace oms std::vector nStates; std::vector nEventIndicators; + /** + * fmi-ls-dae. A component that declares a DAE formulation is integrated in + * DAE mode: it owns no explicit ODE, so the master carries its algebraic + * variables beside the states and drives the residuals of + * F(t, x, x', z) = 0 to zero with IDA. + * + * The unknowns are laid out per component, [x_i | z_i], and so are the + * rows: a DAE component contributes its residuals (one per state and per + * algebraic variable, checked square when DAE mode is switched on), an ODE + * component the explicit rows x'_i - f_i(x, u, t). So a system may mix the + * two, and one with no DAE component at all is not affected by any of this. + */ + /// Index-aligned with fmus: the component as an FMI 3.0 ME one when it runs + /// in DAE mode, else nullptr. + std::vector daeFmus; + std::vector nAlgebraic; ///< per component, 0 unless in DAE mode + std::vector algebraicVars; + std::vector residuals; + bool daeMode = false; ///< any component runs in DAE mode + size_t nDaeUnknowns = 0; ///< states, algebraic variables and lifted loop connections + + /** + * An algebraic loop lifted into the global system. The loop's connections + * are what the inner Newton solved for; here each becomes one unknown — the + * value of the driven input — and one row, output - input = 0. So a loop + * spanning components is resolved by the integrator, with its own error + * control, instead of being iterated to convergence inside every residual + * evaluation on top of whatever each FMU hides internally. + * + * These unknowns follow the components' in y, and their rows the components' + * rows. They are algebraic: their entry in IDASetId's vector is 0. + */ + struct LiftedConnection + { + int output; ///< graph node driving the connection + int input; ///< graph node it drives, and the unknown + }; + std::vector liftedConnections; + /// Per loop of the simulation graph, whether it was lifted; passed to + /// updateInputsInternal so the inner solver leaves those alone. + std::vector liftedLoops; + + /// Per unknown, the rows that reach it; empty when the structure is unknown + /// and IDA differences the Jacobian itself. + std::vector> daeJacRowsByColumn; + std::vector daeJacColorOfColumn; + std::vector> daeJacColumnsOfColor; + /// The increment each column was perturbed by, kept between the two passes. + std::vector daeJacIncrement; + /// Whether the Jacobian is held sparse and factorized by KLU, which needs a + /// structure; without one it is dense and IDA differences it itself. + bool daeJacSparse = false; + std::vector states; std::vector states_der; std::vector states_nominal; @@ -124,15 +220,33 @@ namespace oms N_Vector abstol; }; + struct SolverDataIDA_t + { + SUNContext sunctx; /* SUNDIALS simulation context */ + void *mem; + N_Vector y; /* [states | algebraic variables], per component */ + N_Vector yp; /* the state derivatives; unused rows for the algebraic variables */ + N_Vector id; /* 1.0 differential, 0.0 algebraic — what IDASetId takes */ + SUNLinearSolver linSol; + SUNMatrix J; + N_Vector abstol; + N_Vector ewt; /* the error weights, for the Jacobian's increments */ + }; + union SolverData_t { SolverDataEuler_t euler; SolverDataCVODE_t cvode; + SolverDataIDA_t ida; } solverData; friend int oms::cvode_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); friend int oms::cvode_rhs_algebraic(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); friend int oms::cvode_roots(sunrealtype t, N_Vector y, sunrealtype *gout, void* user_data); + friend int oms::ida_res(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data); + friend int oms::ida_roots(sunrealtype t, N_Vector yy, N_Vector yp, sunrealtype *gout, void* user_data); + friend int oms::ida_jac(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector rr, + SUNMatrix Jac, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); }; } diff --git a/src/OMSimulatorPython/cli.py b/src/OMSimulatorPython/cli.py index d21db1ce1..4acfdb32d 100644 --- a/src/OMSimulatorPython/cli.py +++ b/src/OMSimulatorPython/cli.py @@ -253,7 +253,7 @@ def main(argv=None) -> int: _addOption(parser, '--realTime', 'false', 'Enable experimental feature for (soft) real-time co-simulation', action=argparse.BooleanOptionalAction) _addOption(parser, '--resultFile', "the model name plus '_res.mat'", 'Specify the name of the output result file') _addOption(parser, '--skipCSVHeader', 'true', 'Skip the CSV delimiter row in the header of .csv result files', action=argparse.BooleanOptionalAction) - _addOption(parser, '--solver', 'cvode', 'Set the ODE solver for model-exchange FMUs (.fmu, mode=me only)', choices=['euler', 'cvode']) + _addOption(parser, '--solver', 'cvode', "Set the solver for model-exchange FMUs (.fmu, mode=me only). 'ida' is the DAE mode of fmi-ls-dae and is selected on its own for an FMU that declares a DAE formulation", choices=['euler', 'cvode', 'ida']) _addOption(parser, '--solverStats', 'false', 'Add solver stats to the result file, e.g., step size; not supported for all solvers', action=argparse.BooleanOptionalAction) _addOption(parser, '--startTime', 'from the model', 'Specify the start time', type=float) _addOption(parser, '--stepSize', 'from the model', 'Specify the (maximum) step size (.fmu only)', type=float) diff --git a/src/OMSimulatorPython/fmu.py b/src/OMSimulatorPython/fmu.py index 1a403713c..16775972a 100644 --- a/src/OMSimulatorPython/fmu.py +++ b/src/OMSimulatorPython/fmu.py @@ -56,11 +56,13 @@ # ODE solver methods usable for model-exchange ('me'/sc) FMUs. Kept local to # avoid importing instantiated_model.py, which imports system.py, which -# imports this module (circular import). +# imports this module (circular import). These are values of oms_solver_enu_t +# in include/OMSimulator/Types.h, as in instantiated_model.SolverType. _FMU_SOLVER = { 'euler': 2, 'cvode': 3, - 'ma':6 + 'ida': 4, + 'ma': 7 } ## default initial attribute table for FMI 2.0, based on variability and causality diff --git a/src/OMSimulatorPython/instantiated_model.py b/src/OMSimulatorPython/instantiated_model.py index 612167d0b..0161c4c60 100644 --- a/src/OMSimulatorPython/instantiated_model.py +++ b/src/OMSimulatorPython/instantiated_model.py @@ -44,12 +44,18 @@ import warnings class SolverType(Enum): - '''Enumeration for solver method to map with c api.''' + '''Enumeration for solver method to map with c api. + + These are the values of oms_solver_enu_t in include/OMSimulator/Types.h and + have to be kept in step with it: inserting a solver there renumbers the ones + after it, and setSolver() then quietly selects the wrong one. + ''' euler = 2 cvode = 3 - oms_ma = 6 - oms_mav = 7 - oms_mav2 = 8 + ida = 4 # DAE mode (fmi-ls-dae) + oms_ma = 7 + oms_mav = 8 + oms_mav2 = 9 class SystemType(Enum): '''Enumeration for system type to map with c api.''' wc = 1 diff --git a/src/OMSimulatorPython/variable.py b/src/OMSimulatorPython/variable.py index 73909c709..3111522dd 100644 --- a/src/OMSimulatorPython/variable.py +++ b/src/OMSimulatorPython/variable.py @@ -45,6 +45,10 @@ class Causality(Enum): calculatedParameter = 3 local = 4 independent = 5 + # FMI 3.0: a parameter that changes the structure of the FMU (the number of + # states, of variables, of an array's dimensions) and can only be set in + # Configuration Mode. fmi-ls-dae's DAE-mode switch is one. + structuralParameter = 6 class SignalType(Enum): @@ -136,6 +140,9 @@ def isParameter(self): def isCalculatedParameter(self): return self.causality == Causality.calculatedParameter + def isStructuralParameter(self): + return self.causality == Causality.structuralParameter + def isContinuous(self): return self.variability == "continuous" diff --git a/testsuite/tests/resources/CMakeLists.txt b/testsuite/tests/resources/CMakeLists.txt index 33b0927d1..731afdf25 100644 --- a/testsuite/tests/resources/CMakeLists.txt +++ b/testsuite/tests/resources/CMakeLists.txt @@ -15,6 +15,7 @@ set(OM_OMS_TESTSUITE_RESOURCES_FMUS ${CMAKE_CURRENT_SOURCE_DIR}/ClockedStepBoundary ${CMAKE_CURRENT_SOURCE_DIR}/Dahlquist ${CMAKE_CURRENT_SOURCE_DIR}/Dahlquist3 + ${CMAKE_CURRENT_SOURCE_DIR}/DaeLoopA ${CMAKE_CURRENT_SOURCE_DIR}/DualMassOscillator.System1 ${CMAKE_CURRENT_SOURCE_DIR}/DualMassOscillator.System2 ${CMAKE_CURRENT_SOURCE_DIR}/Enum1 @@ -28,6 +29,7 @@ set(OM_OMS_TESTSUITE_RESOURCES_FMUS ${CMAKE_CURRENT_SOURCE_DIR}/HelloWorldWithInput ${CMAKE_CURRENT_SOURCE_DIR}/Int1 ${CMAKE_CURRENT_SOURCE_DIR}/Lin2DimODE + ${CMAKE_CURRENT_SOURCE_DIR}/LoopB ${CMAKE_CURRENT_SOURCE_DIR}/Modelica.Blocks.Continuous.Integrator ${CMAKE_CURRENT_SOURCE_DIR}/Modelica.Blocks.Math.Add ${CMAKE_CURRENT_SOURCE_DIR}/Modelica.Blocks.Math.Add3 diff --git a/testsuite/tests/resources/DaeLoopA/binaries/x86_64-linux/DaeLoopA.so b/testsuite/tests/resources/DaeLoopA/binaries/x86_64-linux/DaeLoopA.so new file mode 100755 index 000000000..51a205b16 Binary files /dev/null and b/testsuite/tests/resources/DaeLoopA/binaries/x86_64-linux/DaeLoopA.so differ diff --git a/testsuite/tests/resources/DaeLoopA/extra/org.fmi-standard.fmi-ls-dae/fmi-ls-manifest.xml b/testsuite/tests/resources/DaeLoopA/extra/org.fmi-standard.fmi-ls-dae/fmi-ls-manifest.xml new file mode 100644 index 000000000..10a2fdcd4 --- /dev/null +++ b/testsuite/tests/resources/DaeLoopA/extra/org.fmi-standard.fmi-ls-dae/fmi-ls-manifest.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/testsuite/tests/resources/DaeLoopA/modelDescription.xml b/testsuite/tests/resources/DaeLoopA/modelDescription.xml new file mode 100644 index 000000000..057bb0cc4 --- /dev/null +++ b/testsuite/tests/resources/DaeLoopA/modelDescription.xml @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/testsuite/tests/resources/DaeLoopA/resources/DaeLoopA_JacLSJac0.bin b/testsuite/tests/resources/DaeLoopA/resources/DaeLoopA_JacLSJac0.bin new file mode 100644 index 000000000..eed5342f8 Binary files /dev/null and b/testsuite/tests/resources/DaeLoopA/resources/DaeLoopA_JacLSJac0.bin differ diff --git a/testsuite/tests/resources/DaeLoopA/terminalsAndIcons/terminalsAndIcons.xml b/testsuite/tests/resources/DaeLoopA/terminalsAndIcons/terminalsAndIcons.xml new file mode 100644 index 000000000..5d9a7ac11 --- /dev/null +++ b/testsuite/tests/resources/DaeLoopA/terminalsAndIcons/terminalsAndIcons.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + \ No newline at end of file diff --git a/testsuite/tests/resources/LoopB/binaries/x86_64-linux/LoopB.so b/testsuite/tests/resources/LoopB/binaries/x86_64-linux/LoopB.so new file mode 100755 index 000000000..c5517fbeb Binary files /dev/null and b/testsuite/tests/resources/LoopB/binaries/x86_64-linux/LoopB.so differ diff --git a/testsuite/tests/resources/LoopB/modelDescription.xml b/testsuite/tests/resources/LoopB/modelDescription.xml new file mode 100644 index 000000000..722259686 --- /dev/null +++ b/testsuite/tests/resources/LoopB/modelDescription.xml @@ -0,0 +1,46 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/testsuite/tests/resources/LoopB/terminalsAndIcons/terminalsAndIcons.xml b/testsuite/tests/resources/LoopB/terminalsAndIcons/terminalsAndIcons.xml new file mode 100644 index 000000000..f79f1973b --- /dev/null +++ b/testsuite/tests/resources/LoopB/terminalsAndIcons/terminalsAndIcons.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + \ No newline at end of file diff --git a/testsuite/tests/simulation/CMakeLists.txt b/testsuite/tests/simulation/CMakeLists.txt index bae7b435f..ae5bab293 100644 --- a/testsuite/tests/simulation/CMakeLists.txt +++ b/testsuite/tests/simulation/CMakeLists.txt @@ -1,5 +1,7 @@ oms_add_test( clockedStepBoundary.py + DaeMode1.py + DaeMode2.py DualMassOscillator1.py DualMassOscillator2.py exportJson1.py @@ -7,6 +9,7 @@ oms_add_test( exportJson4.py fmuSimulation1.py fmuSimulation2.py + fmuSimulation3.py LinearTransformation1.py LinearTransformation2.py LinearTransformation3.py diff --git a/testsuite/tests/simulation/DaeMode1.py b/testsuite/tests/simulation/DaeMode1.py new file mode 100644 index 000000000..925eb2e07 --- /dev/null +++ b/testsuite/tests/simulation/DaeMode1.py @@ -0,0 +1,82 @@ +## status: correct +## linux: yes +## ucrt64: no +## win: no +## mac: no +## asan: no + +from OMSimulator import SSP, CRef, Settings + +Settings.suppressPath = True + +# A component in DAE mode (fmi-ls-dae) inside an algebraic loop with an ordinary +# component: DaeLoopA declares its residuals in the fmi-ls-dae manifest, LoopB is +# an explicit FMU, and the two connections close the loop A.y -> B.p, B.q -> A.u. +# The system switches to IDA by itself, lifts the loop into rows of the DAE system +# and factorizes the coloured Jacobian with KLU. +# +# Both FMUs are exported by OpenModelica with --simCodeTarget=C+Rust, so they only +# carry x86_64-linux binaries: +# +# model DaeLoopA +# input Real u; +# output Real y; +# Real x(start = 1, fixed = true); +# Real z; +# equation +# der(x) = -x + u; +# z + 2*y = x; +# z - y = u; +# annotation(__OpenModelica_commandLineOptions = "--daeMode"); +# end DaeLoopA; +# +# model LoopB +# input Real p; +# output Real q; +# equation +# q = 0.3*p + 0.5; +# end LoopB; +# +# Closed loop: u = (0.1*x + 0.5)/1.1, so x(t) = 0.5 + 0.5*exp(-t/1.1) and +# x(1) = 0.701445. +model = SSP() +model.addResource('../resources/DaeLoopA.fmu', new_name='resources/DaeLoopA.fmu') +model.addResource('../resources/LoopB.fmu', new_name='resources/LoopB.fmu') + +model.addComponent(CRef('default', 'A'), 'resources/DaeLoopA.fmu') +model.addComponent(CRef('default', 'B'), 'resources/LoopB.fmu') + +model.addConnection(CRef('default', 'A', 'y'), CRef('default', 'B', 'p')) +model.addConnection(CRef('default', 'B', 'q'), CRef('default', 'A', 'u')) + +model.export('DaeMode1.ssp') + +model2 = SSP('DaeMode1.ssp') +instantiated_model = model2.instantiate() +instantiated_model.setResultFile('DaeMode1_res.mat') +instantiated_model.initialize() +instantiated_model.simulate() +for name in [('A', 'x'), ('A', 'y'), ('A', 'u'), ('B', 'q')]: + print(f"info: {'.'.join(name)} = {instantiated_model.getValue(CRef('default', *name)):.6f}", flush=True) +instantiated_model.terminate() +instantiated_model.delete() +## Result: +## info: FMU "A" carries a DAE formulation (fmi-ls-dae 1.0.0-alpha.1): 3 residuals over 1 states and 2 algebraic variables +## info: Alg. loop (size 2/4) +## model.root.B +## model.root.A +## info: Alg. loop (size 2/4) +## model.root.B +## model.root.A +## info: system "model.root": 2 connection(s) of 1 algebraic loop(s) are rows of the DAE system, not an inner solve +## info: system "model.root" contains a component in DAE mode (fmi-ls-dae); integrating with IDA instead of cvode +## info: system "model.root": DAE Jacobian 5x5, 16 structural nonzeros, 5 colour(s) — sparse (KLU), one residual evaluation per colour instead of one per unknown +## info: Result file: DaeMode1_res.mat (bufferSize=10) +## info: A.x = 0.701445 +## info: A.y = 0.061044 +## info: A.u = 0.518313 +## info: B.q = 0.518313 +## info: Final Statistics for 'model.root': +## NumSteps = 1010 NumResEvals = 1014 NumLinSolvSetups = 13 +## NumNonlinSolvIters = 1012 NumNonlinSolvConvFails = 0 NumErrTestFails = 0 +## endResult diff --git a/testsuite/tests/simulation/DaeMode2.py b/testsuite/tests/simulation/DaeMode2.py new file mode 100644 index 000000000..5ce4ade25 --- /dev/null +++ b/testsuite/tests/simulation/DaeMode2.py @@ -0,0 +1,45 @@ +## status: correct +## linux: yes +## ucrt64: no +## win: no +## mac: no +## asan: no + +from OMSimulator import FMU + +# A single model-exchange FMU in DAE mode (fmi-ls-dae), through the FMU API. +# DaeLoopA is described in DaeMode1.py; here its input is a constant, u = 0.5, +# so x(t) = 0.5 + 0.5*exp(-t), y = (x - u)/3 and z = u + y. +# Selecting IDA by name and leaving the default CVODE give the same result: +# a DAE-mode component has no explicit ODE, so the system switches to IDA itself. + +for solver in ['ida', 'default']: + fmu = FMU('../resources/DaeLoopA.fmu', instanceName=f'DaeLoopA_{solver}') + fmu.stopTime = 1.0 + fmu.instantiate() + if solver != 'default': + fmu.setSolver(solver) + fmu.setResultFile(f'DaeMode2_{solver}_res.mat') + fmu.setValue('u', 0.5) + fmu.initialize() + fmu.simulate() + print(f"info: {solver}: x = {fmu.getValue('x'):.6f}, y = {fmu.getValue('y'):.6f}, z = {fmu.getValue('z'):.6f}", flush=True) + fmu.terminate() + fmu.delete() +## Result: +## info: FMU "DaeLoopA_ida" carries a DAE formulation (fmi-ls-dae 1.0.0-alpha.1): 3 residuals over 1 states and 2 algebraic variables +## info: system "DaeLoopA_ida.root": DAE Jacobian 3x3, 6 structural nonzeros, 3 colour(s) — sparse (KLU), one residual evaluation per colour instead of one per unknown +## info: Result file: DaeMode2_ida_res.mat (bufferSize=10) +## info: ida: x = 0.683940, y = 0.061313, z = 0.561313 +## info: Final Statistics for 'DaeLoopA_ida.root': +## NumSteps = 511 NumResEvals = 514 NumLinSolvSetups = 13 +## NumNonlinSolvIters = 512 NumNonlinSolvConvFails = 0 NumErrTestFails = 0 +## info: FMU "DaeLoopA_default" carries a DAE formulation (fmi-ls-dae 1.0.0-alpha.1): 3 residuals over 1 states and 2 algebraic variables +## info: system "DaeLoopA_default.root" contains a component in DAE mode (fmi-ls-dae); integrating with IDA instead of cvode +## info: system "DaeLoopA_default.root": DAE Jacobian 3x3, 6 structural nonzeros, 3 colour(s) — sparse (KLU), one residual evaluation per colour instead of one per unknown +## info: Result file: DaeMode2_default_res.mat (bufferSize=10) +## info: default: x = 0.683940, y = 0.061313, z = 0.561313 +## info: Final Statistics for 'DaeLoopA_default.root': +## NumSteps = 511 NumResEvals = 514 NumLinSolvSetups = 13 +## NumNonlinSolvIters = 512 NumNonlinSolvConvFails = 0 NumErrTestFails = 0 +## endResult diff --git a/testsuite/tests/simulation/fmuSimulation3.py b/testsuite/tests/simulation/fmuSimulation3.py new file mode 100644 index 000000000..0c7108348 --- /dev/null +++ b/testsuite/tests/simulation/fmuSimulation3.py @@ -0,0 +1,40 @@ +## status: correct +## linux: yes +## ucrt64: yes +## win: yes +## mac: yes +## asan: yes + +from OMSimulator import FMU + +# This test case selects every solver the FMU API offers, by name, for both +# kinds of an FMU that exports Co-Simulation and Model Exchange. The names +# are mapped to the values of oms_solver_enu_t; a solver added to that enum +# renumbers the ones after it, and a stale value makes setSolver fail. + +for mode, solvers in [('cs', ['ma']), ('me', ['euler', 'cvode'])]: + for solver in solvers: + fmu = FMU('../resources/Dahlquist3.fmu', instanceName=f'Dahlquist_{mode}_{solver}') + fmu.mode = mode + fmu.stopTime = 1.0 + fmu.stepSize = 0.1 + fmu.instantiate() + fmu.setSolver(solver) + fmu.setResultFile(f'fmuSimulation3_{mode}_{solver}_res.mat') + fmu.initialize() + fmu.simulate() + print(f"info: {mode} {solver}: x = {fmu.getValue('x'):.4f}", flush=True) + fmu.terminate() + fmu.delete() +## Result: +## info: Result file: fmuSimulation3_cs_ma_res.mat (bufferSize=10) +## info: cs ma: x = 0.3487 +## info: Result file: fmuSimulation3_me_euler_res.mat (bufferSize=10) +## info: me euler: x = 0.3487 +## info: maximum step size for 'Dahlquist_me_cvode.root': 0.100000 +## info: Result file: fmuSimulation3_me_cvode_res.mat (bufferSize=10) +## info: me cvode: x = 0.3679 +## info: Final Statistics for 'Dahlquist_me_cvode.root': +## NumSteps = 29 NumRhsEvals = 39 NumLinSolvSetups = 9 +## NumNonlinSolvIters = 38 NumNonlinSolvConvFails = 0 NumErrTestFails = 1 +## endResult