diff --git a/CHANGELOG.md b/CHANGELOG.md index ef67acf15..9c64b723d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -59,6 +59,7 @@ - Remove data copying between system and components in `PowerElectronics` models. - Added multi-contingency analysis application. - Added `BusToSignalAdapter` component for communicating bus voltages and injection currents. +- Added support for running IDA with fixed time steps ## v0.1 diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index d932f6565..c36b55e56 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -33,12 +33,23 @@ namespace GridKit { } - virtual int allocate() = 0; - virtual int initialize() = 0; - virtual int tagDifferentiable() = 0; - virtual int evaluateResidual() = 0; - virtual int evaluateJacobian() = 0; - virtual int evaluateIntegrand() = 0; + virtual int allocate() = 0; + virtual int initialize() = 0; + virtual int tagDifferentiable() = 0; + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + virtual int setAbsoluteTolerance(RealT rel_tol) = 0; + virtual int evaluateResidual() = 0; + virtual int evaluateJacobian() = 0; + virtual int evaluateIntegrand() = 0; virtual int initializeAdjoint() = 0; virtual int evaluateAdjointResidual() = 0; @@ -103,11 +114,26 @@ namespace GridKit */ virtual bool hasJacobian() = 0; - virtual IdxT sizeQuadrature() = 0; - virtual IdxT sizeParams() = 0; - virtual void updateTime(RealT t, RealT a) = 0; - virtual void setTolerances(RealT& rtol, RealT& atol) const = 0; - virtual void setMaxSteps(IdxT& msa) const = 0; + virtual IdxT sizeQuadrature() = 0; + virtual IdxT sizeParams() = 0; + virtual void updateTime(RealT t, RealT a) = 0; + + /** + * @brief Get the absolute tolerance for each variable in the model + * + * @return a reference to the absolute tolerance vector. + * + * @pre `setAbsoluteTolerance` must have been called first. + */ + virtual std::vector& absoluteTolerance() = 0; + /** + * @brief Get the absolute tolerance for each variable in the model + * + * @return a const reference to the absolute tolerance vector. + * + * @pre `setAbsoluteTolerance` must have been called first. + */ + virtual const std::vector& absoluteTolerance() const = 0; virtual std::vector& y() = 0; virtual const std::vector& y() const = 0; diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index 5355b0b8a..cea3fce7d 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -73,6 +73,7 @@ namespace GridKit virtual int allocate() override final; virtual int initialize() override final; virtual int tagDifferentiable() override final; + virtual int setAbsoluteTolerance(RealT rel_tol) override final; virtual int evaluateResidual() override final; virtual int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp index bcb0c268d..dbfc7b83b 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp @@ -140,6 +140,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Branch::setAbsoluteTolerance(RealT) + { + return 0; + } + /** * @brief Bus 1 residual contribution from bus 1 variables * diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp index a94afe863..c05370a36 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp @@ -28,6 +28,7 @@ namespace GridKit using BusBase::J_cols_buffer_; using BusBase::J_vals_buffer_; using BusBase::tag_; + using BusBase::abs_tol_; using BusBase::variable_indices_; using BusBase::residual_indices_; @@ -44,6 +45,7 @@ namespace GridKit virtual int setBusID(IdxT) override final; virtual int allocate() override final; virtual int tagDifferentiable() override final; + virtual int setAbsoluteTolerance(RealT rel_tol) override final; virtual int initialize() override final; virtual int evaluateResidual() override final; virtual int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 77fe6ecf0..c75a44c7c 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -86,6 +86,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -120,6 +121,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Bus::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /*! * @brief initialize method sets bus variables to stored initial values. */ diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp index 9f9866994..b1730b0f8 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp @@ -38,6 +38,7 @@ namespace GridKit virtual int setBusID(IdxT) override final; virtual int allocate() override final; virtual int tagDifferentiable() override final; + virtual int setAbsoluteTolerance(RealT rel_tol) override final; virtual int initialize() override final; virtual int evaluateResidual() override final; virtual int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp index 56d05d82a..60a244235 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp @@ -95,6 +95,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int BusInfinite::setAbsoluteTolerance(RealT) + { + return 0; + } + /*! * @brief initialize method sets bus variables to stored initial values. */ diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp index ba0d4784f..1ee7742ba 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp @@ -52,6 +52,7 @@ namespace GridKit int allocate() override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; int evaluateResidual() override final; int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 97761a533..e20b474f9 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -141,6 +141,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int BusFault::setAbsoluteTolerance(RealT) + { + return 0; + } + /** * @brief Bus residual * diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index 507ee52b1..b4033cdf4 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -73,6 +73,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::abs_tol_; using Component::time_; using Component::y_; using Component::yp_; @@ -105,6 +106,7 @@ namespace GridKit int verify() const override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; int evaluateResidual() override final; int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index f535862f0..7a3542a98 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -122,6 +122,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -283,6 +284,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Ieeet1::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Internal Residual * diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp index 87c143070..536da57d2 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp @@ -60,6 +60,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::abs_tol_; using Component::time_; using Component::y_; using Component::yp_; @@ -87,6 +88,7 @@ namespace GridKit int verify() const override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; int evaluateResidual() override final; int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index f09fd3b5f..cd541a5fe 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -62,6 +62,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -180,6 +181,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int SexsPti::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + template __attribute__((always_inline)) inline int SexsPti::evaluateInternalResidual( ScalarT* y, diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index 2b2a86ad7..b9368ba5b 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -64,6 +64,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::abs_tol_; using Component::time_; using Component::y_; using Component::yp_; @@ -91,6 +92,7 @@ namespace GridKit int verify() const override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT) override final; int evaluateResidual() override final; // Still to be implemented diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 836e9c3ff..93fa08944 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -140,6 +140,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -232,6 +233,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Tgov1::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Internal residuals * diff --git a/GridKit/Model/PhasorDynamics/GridElement.hpp b/GridKit/Model/PhasorDynamics/GridElement.hpp index aa961bfc7..e13069e22 100644 --- a/GridKit/Model/PhasorDynamics/GridElement.hpp +++ b/GridKit/Model/PhasorDynamics/GridElement.hpp @@ -54,17 +54,6 @@ namespace GridKit return nnz_; } - void setTolerances(RealT& rtol, RealT& atol) const override - { - rtol = rel_tol_; - atol = abs_tol_; - } - - void setMaxSteps(IdxT& msa) const override - { - msa = max_steps_; - } - std::vector& y() override { return y_; @@ -95,6 +84,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() override + { + return abs_tol_; + } + + const std::vector& absoluteTolerance() const override + { + return abs_tol_; + } + std::vector& getResidual() override { return f_; @@ -158,6 +157,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector abs_tol_; std::vector f_; std::vector g_; @@ -166,11 +166,6 @@ namespace GridKit IdxT* J_cols_buffer_{nullptr}; RealT* J_vals_buffer_{nullptr}; - RealT rel_tol_; - RealT abs_tol_; - - IdxT max_steps_; - // // Adjoint sensitivity members // diff --git a/GridKit/Model/PhasorDynamics/Load/Load.hpp b/GridKit/Model/PhasorDynamics/Load/Load.hpp index 6dd8390d5..bcf7d0ce3 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.hpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.hpp @@ -36,6 +36,7 @@ namespace GridKit using Component::alpha_; using Component::y_; using Component::yp_; + using Component::abs_tol_; using Component::tag_; using Component::wb_; using Component::h_; @@ -62,6 +63,7 @@ namespace GridKit virtual int allocate() override final; virtual int initialize() override final; virtual int tagDifferentiable() override final; + virtual int setAbsoluteTolerance(RealT) override final; virtual int evaluateResidual() override final; virtual int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp index e774074ea..ea289e47a 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -89,6 +90,7 @@ namespace GridKit f_.resize(size); y_.resize(size); yp_.resize(size); + abs_tol_.resize(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -140,6 +142,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Load::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Bus residual * diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp index 6ebfbd029..f90a9c6e0 100644 --- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp +++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp @@ -37,6 +37,7 @@ namespace GridKit using Component::y_; using Component::yp_; using Component::tag_; + using Component::abs_tol_; using Component::wb_; using Component::h_; using Component::f_; @@ -62,6 +63,7 @@ namespace GridKit int allocate() override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; int evaluateResidual() override final; int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp index 6896662d1..dfcc995db 100644 --- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp +++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp @@ -106,6 +106,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -158,6 +159,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int LoadZIP::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Bus residual * diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index 10a137d2f..f3f692343 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -66,6 +66,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::abs_tol_; using Component::time_; using Component::y_; using Component::yp_; @@ -93,6 +94,7 @@ namespace GridKit int verify() const override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; int evaluateResidual() override final; int evaluateJacobian() override final; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index e428c601a..f08e2a6aa 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -161,6 +161,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -243,6 +244,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Ieeest::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + template __attribute__((always_inline)) inline int Ieeest::evaluateInternalResidual( ScalarT* y, diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp index e0bccb2c0..afbf94745 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp @@ -75,6 +75,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::abs_tol_; using Component::time_; using Component::y_; using Component::yp_; @@ -134,6 +135,7 @@ namespace GridKit int verify() const override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT) override final; int evaluateResidual() override final; // Still to be implemented diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index b7ccac6d5..8993a55e3 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -333,6 +333,7 @@ namespace GridKit tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); + abs_tol_.resize(size); // Resize bus data wb_.resize(2); @@ -523,6 +524,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Genrou::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Internal residual * diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp index 93ba24dfd..de70a76c2 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp @@ -70,6 +70,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::abs_tol_; using Component::time_; using Component::y_; using Component::yp_; @@ -98,6 +99,7 @@ namespace GridKit int verify() const override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; int evaluateResidual() override final; // Still to be implemented diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp index 12227e69e..bdc1e7132 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp @@ -200,6 +200,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -346,6 +347,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Gensal::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Internal residual * diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index 2bf77e2ce..a25e71496 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -39,6 +39,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::abs_tol_; using Component::time_; using Component::y_; using Component::yp_; @@ -75,6 +76,7 @@ namespace GridKit int allocate() override final; int initialize() override final; int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT) override; int evaluateResidual() override final; int verify() const override final diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index e3836912f..f9461c133 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -184,6 +184,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -248,6 +249,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int GenClassical::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Internal residual * diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 25478738c..75f196cd5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -46,10 +46,9 @@ namespace GridKit using PhasorDynamics::Component::y_; using PhasorDynamics::Component::yp_; using PhasorDynamics::Component::tag_; + using PhasorDynamics::Component::abs_tol_; using PhasorDynamics::Component::f_; using PhasorDynamics::Component::J_; - using PhasorDynamics::Component::rel_tol_; - using PhasorDynamics::Component::abs_tol_; using PhasorDynamics::Component::variable_indices_; using PhasorDynamics::Component::residual_indices_; @@ -59,10 +58,6 @@ namespace GridKit */ SystemModel() { - // Set system model tolerances - rel_tol_ = 1e-7; - abs_tol_ = 1e-9; - this->max_steps_ = 2000; } /** @@ -83,11 +78,6 @@ namespace GridKit using namespace Exciter; using namespace Stabilizer; - // Set system model tolerances - rel_tol_ = 1e-7; - abs_tol_ = 1e-9; - this->max_steps_ = 2000; - owns_components_ = true; // Store parsed system bases before constructing data-driven components. @@ -481,6 +471,7 @@ namespace GridKit yp_.resize(size_); f_.resize(size_); tag_.resize(size_); + abs_tol_.resize(size_); variable_indices_.resize(size_); residual_indices_.resize(size_); @@ -681,6 +672,43 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) override + { + // Set initial values for global solution vectors + IdxT offset = 0; + for (const auto& bus : buses_) + { + bus->setAbsoluteTolerance(rel_tol); + for (IdxT j = 0; j < bus->size(); ++j) + { + abs_tol_[offset + j] = bus->absoluteTolerance()[j]; + } + offset += bus->size(); + } + + for (const auto& component : components_) + { + component->setAbsoluteTolerance(rel_tol); + for (IdxT j = 0; j < component->size(); ++j) + { + abs_tol_[offset + j] = component->absoluteTolerance()[j]; + } + offset += component->size(); + } + + return 0; + } + /** * @brief Compute system residual vector * diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index a99257454..e852dcfc3 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -99,6 +99,7 @@ namespace GridKit y_.resize(static_cast(size_)); yp_.resize(static_cast(size_)); + abs_tol_.resize(static_cast(size_)); f_.resize(static_cast(size_)); return 0; @@ -259,17 +260,6 @@ namespace GridKit return size_opt_; } - void setTolerances(RealT& rel_tol, RealT& abs_tol) const final - { - rel_tol = rel_tol_; - abs_tol = abs_tol_; - } - - void setMaxSteps(IdxT& msa) const final - { - msa = max_steps_; - } - std::vector& y() final { return y_; @@ -300,6 +290,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() final + { + return abs_tol_; + } + + const std::vector& absoluteTolerance() const final + { + return abs_tol_; + } + std::vector& yB() final { return yB_; @@ -437,6 +437,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector abs_tol_; std::vector f_; std::vector g_; @@ -452,9 +453,6 @@ namespace GridKit RealT time_; RealT alpha_; - RealT rel_tol_; - RealT abs_tol_; - IdxT max_steps_; IdxT idc_; diff --git a/GridKit/Model/PowerElectronics/CircuitNode.hpp b/GridKit/Model/PowerElectronics/CircuitNode.hpp index 5fd096f62..de0afd8a8 100644 --- a/GridKit/Model/PowerElectronics/CircuitNode.hpp +++ b/GridKit/Model/PowerElectronics/CircuitNode.hpp @@ -74,6 +74,7 @@ namespace GridKit yp_.resize(size); f_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_[0] = 0; residual_indices_[0] = 0; @@ -102,6 +103,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Node does not compute residuals, so here we just reset residual values. * @@ -163,6 +182,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector abs_tol_; std::vector f_; MatrixT J_; @@ -180,9 +200,6 @@ namespace GridKit RealT time_{0}; RealT alpha_{0}; - RealT rel_tol_{0}; - RealT abs_tol_{0}; - IdxT max_steps_{0}; public: @@ -211,17 +228,6 @@ namespace GridKit // No time to update in node models } - void setTolerances(RealT& rel_tol, RealT& abs_tol) const final - { - rel_tol = rel_tol_; - abs_tol = abs_tol_; - } - - void setMaxSteps(IdxT& msa) const final - { - msa = max_steps_; - } - std::vector& y() final { return y_; @@ -252,6 +258,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() final + { + return abs_tol_; + } + + const std::vector& absoluteTolerance() const final + { + return abs_tol_; + } + std::vector& yB() final { return yB_; diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index 1458c3db2..797022531 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -77,6 +77,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int DistributedGenerator::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Contributes to the resisdual of the Distributed Generator * diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp index 3632bde33..505311cea 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -51,6 +51,7 @@ namespace GridKit using CircuitComponent::y_int_; using CircuitComponent::yp_; using CircuitComponent::yp_int_; + using CircuitComponent::abs_tol_; using CircuitComponent::tag_; using CircuitComponent::f_; using CircuitComponent::f_int_; @@ -77,6 +78,7 @@ namespace GridKit int initialize(); int allocate() final; int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateInternalResidual() final; int evaluateExternalResidual() final; int evaluateJacobian(); diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp index 00248ad29..5a8527c55 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp @@ -51,6 +51,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Inductor::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Compute the resisdual of the component * diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp index ba02ff7a9..470a91741 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp @@ -32,6 +32,7 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::yp_int_; using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; using CircuitComponent::f_; using CircuitComponent::f_int_; using CircuitComponent::g_; @@ -53,6 +54,7 @@ namespace GridKit int initialize(); int allocate() final; int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateInternalResidual() final; int evaluateExternalResidual() final; int evaluateJacobian(); diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp index f9174ba15..c6e78d2b3 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -56,6 +56,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int MicrogridBusDQ::setAbsoluteTolerance(RealT) + { + return 0; + } + template int MicrogridBusDQ::evaluateInternalResidual() { diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index e44c131c4..05c02d65e 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -32,6 +32,7 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::yp_int_; using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; using CircuitComponent::f_; using CircuitComponent::f_int_; using CircuitComponent::g_; @@ -53,6 +54,7 @@ namespace GridKit int initialize(); int allocate() final; int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateInternalResidual() final; int evaluateExternalResidual() final; int evaluateJacobian(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index b5f6652dc..c0896889a 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -65,6 +65,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int MicrogridLine::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Evaluate residual of microgrid line * diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp index c36ce9632..f63826a43 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -32,6 +32,7 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::yp_int_; using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; using CircuitComponent::f_; using CircuitComponent::f_int_; using CircuitComponent::g_; @@ -53,6 +54,7 @@ namespace GridKit int initialize(); int allocate() final; int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateInternalResidual() final; int evaluateExternalResidual() final; int evaluateJacobian(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index caf203c44..d3e9a143d 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -61,6 +61,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int MicrogridLoad::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Eval Micro Load */ diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 99687d3da..dae888a4a 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -32,6 +32,7 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::yp_int_; using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; using CircuitComponent::f_; using CircuitComponent::f_int_; using CircuitComponent::g_; @@ -53,6 +54,7 @@ namespace GridKit int initialize(); int allocate() final; int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateInternalResidual() final; int evaluateExternalResidual() final; int evaluateJacobian(); diff --git a/GridKit/Model/PowerElectronics/NodeBase.hpp b/GridKit/Model/PowerElectronics/NodeBase.hpp index 9e9748529..b2fc80d85 100644 --- a/GridKit/Model/PowerElectronics/NodeBase.hpp +++ b/GridKit/Model/PowerElectronics/NodeBase.hpp @@ -49,17 +49,6 @@ namespace GridKit { } - void setTolerances(RealT& rtol, RealT& atol) const final - { - atol = atol_; - rtol = rtol_; - } - - void setMaxSteps(IdxT& msa) const final - { - msa = max_steps_; - } - std::vector& y() final { return y_; @@ -90,6 +79,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() final + { + return abs_tol_; + } + + const std::vector& absoluteTolerance() const final + { + return abs_tol_; + } + MatrixT& getJacobian() final { return J_; @@ -147,6 +146,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + abs_tol_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); @@ -160,6 +160,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) final + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + int initialize() override { // TODO: fill this in @@ -189,6 +207,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector abs_tol_; std::vector f_; MatrixT J_; @@ -196,11 +215,6 @@ namespace GridKit IdxT* J_cols_buffer_{nullptr}; RealT* J_vals_buffer_{nullptr}; - RealT rtol_; - RealT atol_; - - IdxT max_steps_; - // // Adjoint sensitivity members // diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp index 9c2243010..87d86cab4 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp @@ -52,6 +52,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Resistor::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Computes the resistors resisdual * diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp index b5ad84e80..a6a08b369 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp @@ -32,6 +32,7 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::yp_int_; using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; using CircuitComponent::f_; using CircuitComponent::f_int_; using CircuitComponent::g_; @@ -53,6 +54,7 @@ namespace GridKit int initialize(); int allocate() final; int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateInternalResidual() final; int evaluateExternalResidual() final; int evaluateJacobian(); diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 3b44f4ec6..725f182ef 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -78,7 +78,7 @@ namespace GridKit using CircuitComponent::yp_int_; using CircuitComponent::f_; using CircuitComponent::f_int_; - using CircuitComponent::rel_tol_; + using CircuitComponent::tag_; using CircuitComponent::abs_tol_; public: @@ -89,35 +89,21 @@ namespace GridKit */ PowerElectronicsModel() { - // Set system model parameters as default - rel_tol_ = 1e-4; - abs_tol_ = 1e-4; - this->max_steps_ = 2000; // By default don't use the jacobian - use_jac_ = false; + use_jac_ = false; } /** * @brief Constructor for the system model * - * @param[in] rel_tol Relative tolerance for the system model - * @param[in] abs_tol Absolute tolerance for the system model * @param[in] use_jac Boolean to choose if to use jacobian - * @param[in] max_steps Maximum number of steps for the system model * * @post System model parameters set as input */ - PowerElectronicsModel(double rel_tol = 1e-4, - double abs_tol = 1e-4, - bool use_jac = false, - IdxT max_steps = 2000) + PowerElectronicsModel(bool use_jac = false) { - // Set system model tolerances from input - rel_tol_ = rel_tol; - abs_tol_ = abs_tol; - this->max_steps_ = max_steps; // Can choose if to use jacobian - use_jac_ = use_jac; + use_jac_ = use_jac; } /** @@ -196,6 +182,8 @@ namespace GridKit y_.resize(size_); yp_.resize(size_); f_.resize(size_); + tag_.resize(size_); + abs_tol_.resize(size_); { // Start node internal indexing after all component internals for proper KLU ordering size_t node_internal_idx = component_internal_size; @@ -371,6 +359,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) final + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Evaluate Residuals at each component then collect them * diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp index 1f742bade..f2967c7bd 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp @@ -52,6 +52,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int VoltageSource::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Evaluate resisdual of component */ diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp index aa79ddd14..b1509468d 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp @@ -32,6 +32,7 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::yp_int_; using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; using CircuitComponent::f_; using CircuitComponent::f_int_; using CircuitComponent::g_; @@ -53,6 +54,7 @@ namespace GridKit int initialize(); int allocate() final; int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateInternalResidual() final; int evaluateExternalResidual() final; int evaluateJacobian() final; diff --git a/GridKit/Model/PowerFlow/Branch/Branch.cpp b/GridKit/Model/PowerFlow/Branch/Branch.cpp index 95048043c..ff24c0cb0 100644 --- a/GridKit/Model/PowerFlow/Branch/Branch.cpp +++ b/GridKit/Model/PowerFlow/Branch/Branch.cpp @@ -96,6 +96,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Branch::setAbsoluteTolerance(RealT) + { + return 0; + } + /** * \brief Residual contribution of the branch is pushed to the * two terminal buses. diff --git a/GridKit/Model/PowerFlow/Branch/Branch.hpp b/GridKit/Model/PowerFlow/Branch/Branch.hpp index 57c49813f..0e478f001 100644 --- a/GridKit/Model/PowerFlow/Branch/Branch.hpp +++ b/GridKit/Model/PowerFlow/Branch/Branch.hpp @@ -53,6 +53,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Bus/BaseBus.hpp b/GridKit/Model/PowerFlow/Bus/BaseBus.hpp index 75068f100..5002f320e 100644 --- a/GridKit/Model/PowerFlow/Bus/BaseBus.hpp +++ b/GridKit/Model/PowerFlow/Bus/BaseBus.hpp @@ -25,8 +25,6 @@ namespace GridKit using ModelEvaluatorImpl::nnz_; using ModelEvaluatorImpl::time_; using ModelEvaluatorImpl::alpha_; - using ModelEvaluatorImpl::rel_tol_; - using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::yp_; using ModelEvaluatorImpl::tag_; @@ -76,6 +74,11 @@ namespace GridKit return 0; } + virtual int setAbsoluteTolerance(RealT) + { + return 0; + } + virtual int evaluateResidual() { return 0; diff --git a/GridKit/Model/PowerFlow/Bus/BusPQ.cpp b/GridKit/Model/PowerFlow/Bus/BusPQ.cpp index 7b53a9c89..4cc8a450c 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPQ.cpp +++ b/GridKit/Model/PowerFlow/Bus/BusPQ.cpp @@ -74,6 +74,7 @@ namespace GridKit y_.resize(static_cast(size_)); yp_.resize(static_cast(size_)); tag_.resize(static_cast(size_)); + abs_tol_.resize(static_cast(size_)); fB_.resize(static_cast(size_)); yB_.resize(static_cast(size_)); @@ -90,6 +91,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int BusPQ::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /*! * @brief initialize method sets bus variables to stored initial values. */ diff --git a/GridKit/Model/PowerFlow/Bus/BusPQ.hpp b/GridKit/Model/PowerFlow/Bus/BusPQ.hpp index e247f593d..c3382740d 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPQ.hpp +++ b/GridKit/Model/PowerFlow/Bus/BusPQ.hpp @@ -25,6 +25,7 @@ namespace GridKit using BaseBus::f_; using BaseBus::fB_; using BaseBus::tag_; + using BaseBus::abs_tol_; public: using RealT = typename ModelEvaluatorImpl::RealT; @@ -37,6 +38,7 @@ namespace GridKit virtual int allocate(); virtual int tagDifferentiable(); + virtual int setAbsoluteTolerance(RealT); virtual int initialize(); virtual int evaluateResidual(); virtual int initializeAdjoint(); diff --git a/GridKit/Model/PowerFlow/Bus/BusPV.cpp b/GridKit/Model/PowerFlow/Bus/BusPV.cpp index c3312c70c..bde6dd49b 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPV.cpp +++ b/GridKit/Model/PowerFlow/Bus/BusPV.cpp @@ -72,6 +72,7 @@ namespace GridKit y_.resize(static_cast(size_)); yp_.resize(static_cast(size_)); tag_.resize(static_cast(size_)); + abs_tol_.resize(static_cast(size_)); fB_.resize(static_cast(size_)); yB_.resize(static_cast(size_)); @@ -87,6 +88,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int BusPV::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /*! * @brief initialize method sets bus variables to stored initial values. */ diff --git a/GridKit/Model/PowerFlow/Bus/BusPV.hpp b/GridKit/Model/PowerFlow/Bus/BusPV.hpp index e56ad69e6..4b7816254 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPV.hpp +++ b/GridKit/Model/PowerFlow/Bus/BusPV.hpp @@ -27,6 +27,7 @@ namespace GridKit using BaseBus::f_; using BaseBus::fB_; using BaseBus::tag_; + using BaseBus::abs_tol_; public: using RealT = typename ModelEvaluatorImpl::RealT; @@ -39,6 +40,7 @@ namespace GridKit virtual int allocate(); virtual int tagDifferentiable(); + virtual int setAbsoluteTolerance(RealT); virtual int initialize(); virtual int evaluateResidual(); virtual int initializeAdjoint(); diff --git a/GridKit/Model/PowerFlow/Bus/BusSlack.hpp b/GridKit/Model/PowerFlow/Bus/BusSlack.hpp index d35b7dc3e..65edd619b 100644 --- a/GridKit/Model/PowerFlow/Bus/BusSlack.hpp +++ b/GridKit/Model/PowerFlow/Bus/BusSlack.hpp @@ -23,8 +23,6 @@ namespace GridKit using BaseBus::yp_; using BaseBus::f_; using BaseBus::g_; - using BaseBus::abs_tol_; - using BaseBus::rel_tol_; public: using RealT = typename ModelEvaluatorImpl::RealT; diff --git a/GridKit/Model/PowerFlow/Generator/GeneratorBase.hpp b/GridKit/Model/PowerFlow/Generator/GeneratorBase.hpp index 5750c6863..142184ef2 100644 --- a/GridKit/Model/PowerFlow/Generator/GeneratorBase.hpp +++ b/GridKit/Model/PowerFlow/Generator/GeneratorBase.hpp @@ -65,6 +65,11 @@ namespace GridKit return 0; } + virtual int setAbsoluteTolerance(RealT) + { + return 0; + } + virtual int evaluateResidual() { return 0; diff --git a/GridKit/Model/PowerFlow/Generator2/Generator2.cpp b/GridKit/Model/PowerFlow/Generator2/Generator2.cpp index 0e61a44d5..eb8bd440e 100644 --- a/GridKit/Model/PowerFlow/Generator2/Generator2.cpp +++ b/GridKit/Model/PowerFlow/Generator2/Generator2.cpp @@ -2,6 +2,7 @@ #define _USE_MATH_DEFINES #include "Generator2.hpp" +#include #include #include @@ -49,6 +50,7 @@ namespace GridKit int Generator2::allocate() { tag_.resize(static_cast(size_)); + abs_tol_.resize(static_cast(size_)); return 0; } @@ -60,6 +62,13 @@ namespace GridKit return 0; } + template + int Generator2::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + template int Generator2::initialize() { diff --git a/GridKit/Model/PowerFlow/Generator2/Generator2.hpp b/GridKit/Model/PowerFlow/Generator2/Generator2.hpp index c022279c6..de429e06c 100644 --- a/GridKit/Model/PowerFlow/Generator2/Generator2.hpp +++ b/GridKit/Model/PowerFlow/Generator2/Generator2.hpp @@ -25,6 +25,7 @@ namespace GridKit using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::yp_; using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; @@ -45,6 +46,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Generator4/Generator4.cpp b/GridKit/Model/PowerFlow/Generator4/Generator4.cpp index edf91b3d7..296a56867 100644 --- a/GridKit/Model/PowerFlow/Generator4/Generator4.cpp +++ b/GridKit/Model/PowerFlow/Generator4/Generator4.cpp @@ -58,6 +58,7 @@ namespace GridKit { // std::cout << "Allocate Generator4..." << std::endl; tag_.resize(static_cast(size_)); + abs_tol_.resize(static_cast(size_)); return 0; } @@ -149,6 +150,25 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Generator4::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Computes residual vector for the generator model. * diff --git a/GridKit/Model/PowerFlow/Generator4/Generator4.hpp b/GridKit/Model/PowerFlow/Generator4/Generator4.hpp index de17996a0..d4b163e96 100644 --- a/GridKit/Model/PowerFlow/Generator4/Generator4.hpp +++ b/GridKit/Model/PowerFlow/Generator4/Generator4.hpp @@ -25,6 +25,7 @@ namespace GridKit using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::yp_; using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; @@ -45,6 +46,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp index 284346f6c..6caaa9674 100644 --- a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp +++ b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp @@ -2,6 +2,7 @@ #define _USE_MATH_DEFINES #include "Generator4Governor.hpp" +#include #include #include @@ -67,6 +68,7 @@ namespace GridKit { // std::cout << "Allocate Gen2..." << std::endl; tag_.resize(static_cast(size_)); + abs_tol_.resize(static_cast(size_)); return 0; } @@ -167,6 +169,13 @@ namespace GridKit return 0; } + template + int Generator4Governor::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Computes residual vector for the generator model. * diff --git a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.hpp b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.hpp index 5a34067a1..52f400157 100644 --- a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.hpp +++ b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.hpp @@ -26,6 +26,7 @@ namespace GridKit using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::yp_; using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; @@ -46,6 +47,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp index 18b235510..e08e33f73 100644 --- a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp +++ b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp @@ -2,6 +2,7 @@ #define _USE_MATH_DEFINES #include "Generator4Param.hpp" +#include #include #include @@ -54,6 +55,7 @@ namespace GridKit { // std::cout << "Allocate Generator4Param..." << std::endl; tag_.resize(static_cast(size_)); + abs_tol_.resize(static_cast(size_)); return 0; } @@ -149,6 +151,13 @@ namespace GridKit return 0; } + template + int Generator4Param::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Computes residual vector for the generator model. * diff --git a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.hpp b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.hpp index 3ca69b50f..bc1b94321 100644 --- a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.hpp +++ b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.hpp @@ -25,6 +25,7 @@ namespace GridKit using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::yp_; using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; @@ -45,6 +46,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Load/Load.cpp b/GridKit/Model/PowerFlow/Load/Load.cpp index 366f56242..087c7f6ec 100644 --- a/GridKit/Model/PowerFlow/Load/Load.cpp +++ b/GridKit/Model/PowerFlow/Load/Load.cpp @@ -70,6 +70,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + template + int Load::setAbsoluteTolerance(RealT) + { + return 0; + } + /** * @brief Contributes to the bus residual. * diff --git a/GridKit/Model/PowerFlow/Load/Load.hpp b/GridKit/Model/PowerFlow/Load/Load.hpp index 2326ba5b3..256ccda8e 100644 --- a/GridKit/Model/PowerFlow/Load/Load.hpp +++ b/GridKit/Model/PowerFlow/Load/Load.hpp @@ -46,6 +46,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(RealT); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.cpp b/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.cpp index 25ce8c10d..8dc9e8e5f 100644 --- a/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.cpp +++ b/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.cpp @@ -31,8 +31,6 @@ namespace GridKit B23_(12.0) { // std::cout << "Create a load model with " << size_ << " variables ...\n"; - rel_tol_ = 1e-5; - abs_tol_ = 1e-5; } template diff --git a/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp b/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp index 19d4a1daf..2f69e08bf 100644 --- a/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp +++ b/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp @@ -18,8 +18,6 @@ namespace GridKit using ModelEvaluatorImpl::time_; using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::f_; - using ModelEvaluatorImpl::rel_tol_; - using ModelEvaluatorImpl::abs_tol_; using RealT = typename ModelEvaluatorImpl::RealT; @@ -35,6 +33,11 @@ namespace GridKit return -1; } + int setAbsoluteTolerance(RealT) + { + return -1; + } + int evaluateResidual(); int evaluateJacobian(); diff --git a/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp b/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp index 2457b16b9..40762fcd3 100644 --- a/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp +++ b/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp @@ -77,12 +77,6 @@ namespace GridKit // std::cout << "updateTime: t = " << time_ << "\n"; // } - virtual void setTolerances(RealT& rel_tol, RealT& abs_tol) const - { - rel_tol = rel_tol_; - abs_tol = abs_tol_; - } - virtual void setMaxSteps(IdxT& msa) const { msa = max_steps_; @@ -118,6 +112,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() + { + return abs_tol_; + } + + const std::vector& absoluteTolerance() const + { + return abs_tol_; + } + std::vector& yB() { return yB_; @@ -233,6 +237,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector abs_tol_; std::vector f_; std::vector g_; @@ -250,9 +255,6 @@ namespace GridKit RealT time_; RealT alpha_; - RealT rel_tol_; - RealT abs_tol_; - IdxT max_steps_; IdxT idc_; diff --git a/GridKit/Model/PowerFlow/SystemModel.hpp b/GridKit/Model/PowerFlow/SystemModel.hpp index d6253cd7f..e174b4346 100644 --- a/GridKit/Model/PowerFlow/SystemModel.hpp +++ b/GridKit/Model/PowerFlow/SystemModel.hpp @@ -39,12 +39,11 @@ namespace GridKit using ModelEvaluatorImpl::yB_; using ModelEvaluatorImpl::ypB_; using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::fB_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::gB_; - using ModelEvaluatorImpl::rel_tol_; - using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::param_; using ModelEvaluatorImpl::param_up_; using ModelEvaluatorImpl::param_lo_; @@ -56,10 +55,6 @@ namespace GridKit SystemModel() : ModelEvaluatorImpl(0, 0, 0) { - // Set system model tolerances - rel_tol_ = 1e-7; - abs_tol_ = 1e-9; - this->max_steps_ = 2000; } /** @@ -113,6 +108,7 @@ namespace GridKit f_.resize(size_); fB_.resize(size_); tag_.resize(size_); + abs_tol_.resize(size_); g_.resize(size_quad_); gB_.resize(size_quad_ * size_opt_); @@ -242,6 +238,45 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) + { + // Set initial values for global solution vectors + IdxT offset = 0; + for (const auto& bus : buses_) + { + bus->setAbsoluteTolerance(rel_tol); + for (IdxT j = 0; j < bus->size(); ++j) + { + abs_tol_[offset + j] = bus->absoluteTolerance()[j]; + } + offset += bus->size(); + } + + for (const auto& component : components_) + { + component->setAbsoluteTolerance(rel_tol); + for (IdxT j = 0; j < component->size(); ++j) + { + abs_tol_[offset + j] = component->absoluteTolerance()[j]; + } + offset += component->size(); + } + + return 0; + } + /** * @brief Compute system residual vector * diff --git a/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp b/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp index 88ace8c3b..44b8a75fb 100644 --- a/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp +++ b/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp @@ -47,13 +47,12 @@ namespace GridKit // using ModelEvaluatorImpl::yp_; // using ModelEvaluatorImpl::yB_; // using ModelEvaluatorImpl::ypB_; - // using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::f_; // using ModelEvaluatorImpl::fB_; // using ModelEvaluatorImpl::g_; // using ModelEvaluatorImpl::gB_; - using ModelEvaluatorImpl::rel_tol_; - using ModelEvaluatorImpl::abs_tol_; // using ModelEvaluatorImpl::param_; // using ModelEvaluatorImpl::param_up_; // using ModelEvaluatorImpl::param_lo_; @@ -65,9 +64,6 @@ namespace GridKit SystemSteadyStateModel() : ModelEvaluatorImpl(0, 0, 0) { - // Set system model tolerances - rel_tol_ = 1e-5; - abs_tol_ = 1e-5; } /** @@ -78,9 +74,6 @@ namespace GridKit SystemSteadyStateModel(GridKit::PowerFlowData::SystemModelData mp) : ModelEvaluatorImpl(0, 0, 0) { - rel_tol_ = 1e-5; - abs_tol_ = 1e-5; - // add buses for (auto busdata : mp.bus) { @@ -157,6 +150,8 @@ namespace GridKit // Allocate global vectors y_.resize(size_); f_.resize(size_); + tag_.resize(size_); + abs_tol_.resize(size_); return 0; } @@ -224,6 +219,24 @@ namespace GridKit return 0; } + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + /** * @brief Compute system residual vector * diff --git a/GridKit/Solver/Dynamic/DynamicSolver.hpp b/GridKit/Solver/Dynamic/DynamicSolver.hpp index 74489c229..147f58f78 100644 --- a/GridKit/Solver/Dynamic/DynamicSolver.hpp +++ b/GridKit/Solver/Dynamic/DynamicSolver.hpp @@ -23,6 +23,14 @@ namespace AnalysisManager return model_; } + void setTolerance(ScalarT rel_tol) + { + setTolerance(rel_tol, 0); + } + + virtual void setTolerance(ScalarT rel_tol, ScalarT abs_tol_override) = 0; + virtual void setMaxSteps(IdxT msa) = 0; + protected: GridKit::Model::Evaluator* model_; }; diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index c61328e11..8f6aadbcb 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -83,35 +83,24 @@ namespace AnalysisManager retval = IDASetUserData(solver_, model_); checkOutput(retval, "IDASetUserData"); - // Set tolerances - RealT rel_tol; - RealT abs_tol; - - model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"! - retval = IDASStolerances(solver_, rel_tol, abs_tol); - checkOutput(retval, "IDASStolerances"); - - IdxT msa; - model_->setMaxSteps(msa); - - /// \todo Need to set max number of steps based on user input! - retval = IDASetMaxNumSteps(solver_, static_cast(msa)); - checkOutput(retval, "IDASetMaxNumSteps"); - // Tag differential variables - std::vector& tag = model_->tag(); - if (static_cast(tag.size()) == model_->size()) + const std::vector& tag = model_->tag(); + if (static_cast(tag.size()) != model_->size()) { - tag_ = N_VClone(yy_); - checkAllocation((void*) tag_, "N_VClone"); - model_->tagDifferentiable(); - copyVec(tag, tag_); - - retval = IDASetId(solver_, tag_); - checkOutput(retval, "IDASetId"); - retval = IDASetSuppressAlg(solver_, SUNTRUE); - checkOutput(retval, "IDASetSuppressAlg"); + std::cerr << "\nModel tag size does not match model size.\n\n"; + throw SundialsException(); } + tag_ = N_VClone(yy_); + checkAllocation((void*) tag_, "N_VClone"); + model_->tagDifferentiable(); + copyVec(tag, tag_); + + retval = IDASetId(solver_, tag_); + checkOutput(retval, "IDASetId"); + + abs_tol_ = N_VClone(yy_); + checkAllocation((void*) abs_tol_, "N_VClone"); + setTolerance(DEFAULT_REL_TOL); // Set up linear solver return this->configureLinearSolver(); @@ -375,6 +364,7 @@ namespace AnalysisManager N_VDestroy(yy_); N_VDestroy(yp_); N_VDestroy(tag_); + N_VDestroy(abs_tol_); N_VDestroy(yy0_); N_VDestroy(yp0_); SUNLinSolFree(linearSolver_); @@ -403,12 +393,7 @@ namespace AnalysisManager checkOutput(retval, "IDAQuadInit"); // Set tolerances and error control for quadratures - RealT rel_tol, abs_tol; - model_->setTolerances(rel_tol, abs_tol); - - // Set tolerances for quadrature stricter than for integration - retval = IDAQuadSStolerances(solver_, rel_tol * 0.1, abs_tol * 0.1); - checkOutput(retval, "IDAQuadSStolerances"); + setQuadratureTolerance(0.1 * DEFAULT_REL_TOL); // Include quadrature in eror checking retval = IDASetQuadErrCon(solver_, SUNTRUE); @@ -544,9 +529,7 @@ namespace AnalysisManager template int Ida::initializeBackwardSimulation(RealT tf) { - int retval = 0; - RealT rel_tol; - RealT abs_tol; + int retval = 0; model_->initializeAdjoint(); @@ -561,17 +544,11 @@ namespace AnalysisManager retval = IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_); checkOutput(retval, "IDAInitB"); - model_->setTolerances(rel_tol, abs_tol); - retval = IDASStolerancesB(solver_, backwardID_, rel_tol, abs_tol); - checkOutput(retval, "IDASStolerancesB"); + setBackwardTolerance(DEFAULT_REL_TOL); retval = IDASetUserDataB(solver_, backwardID_, model_); checkOutput(retval, "IDASetUserDataB"); - /// \todo Need to set max number of steps based on user input! - retval = IDASetMaxNumStepsB(solver_, backwardID_, 2000); - checkOutput(retval, "IDASetMaxNumSteps"); - // Allocate Jacobian matrix, if not already if (JacobianMatB_ == nullptr) { @@ -596,9 +573,7 @@ namespace AnalysisManager retval = IDAQuadInitB(solver_, backwardID_, this->adjointIntegrand, qB_); checkOutput(retval, "IDAQuadInitB"); - // retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol*1.1, abs_tol*1.1); - retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol * 0.1, abs_tol * 0.1); - checkOutput(retval, "IDAQuadSStolerancesB"); + setBackwardQuadratureTolerance(0.1 * DEFAULT_REL_TOL); // Include quadratures in error control retval = IDASetQuadErrConB(solver_, backwardID_, SUNTRUE); @@ -923,7 +898,7 @@ namespace AnalysisManager * @tparam IdxT */ template - void Ida::printOutput(RealT t) + void Ida::printOutput(RealT t) const { RealT* yval = N_VGetArrayPointer(yy_); RealT* ypval = N_VGetArrayPointer(yp_); @@ -947,7 +922,7 @@ namespace AnalysisManager * @tparam IdxT */ template - void Ida::printSpecial(RealT t, N_Vector y) + void Ida::printSpecial(RealT t, N_Vector y) const { RealT* yval = N_VGetArrayPointer(y); IdxT N = static_cast(N_VGetLength(y)); @@ -967,7 +942,7 @@ namespace AnalysisManager * @tparam IdxT */ template - void Ida::printFinalStats() + void Ida::printFinalStats() const { int retval = IDAPrintAllStats(solver_, stdout, SUN_OUTPUTFORMAT_TABLE); checkOutput(retval, "IDAPrintAllStats"); @@ -1075,6 +1050,289 @@ namespace AnalysisManager } } + /** + * @brief Set fixed step size and tolerances for the nonlinear solver + * + * @param time_step The fixed step size to use + * @param rel_tol The relative tolerance to use for the nonlinear solver + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance for the nonlinear solver rather than the + * model's default absolute tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setFixedStep(ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override) + { + setFixedStep(solver_, time_step, rel_tol, abs_tol_override); + } + + /** + * @brief Set fixed step size and tolerances for the nonlinear solver for + * the backward simulation + * + * @param time_step The fixed step size to use + * @param rel_tol The relative tolerance to use for the nonlinear solver + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance for the nonlinear solver rather than the + * model's default absolute tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setBackwardFixedStep(ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override) + { + setFixedStep(IDAGetAdjIDABmem(solver_, backwardID_), + time_step, + rel_tol, + abs_tol_override); + } + + /** + * @brief A helper function to set fixed step size and tolerances for the + * nonlinear solver + * + * @param mem The IDA memory (either forward or backward) + * @param time_step The fixed step size to use + * @param rel_tol The relative tolerance to use for the nonlinear solver + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance for the nonlinear solver rather than the + * model's default absolute tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setFixedStep(void* mem, + ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override) + { + int retval = 0; + retval = IDASetMinStep(mem, time_step); + checkOutput(retval, "IDASetMinStep"); + retval = IDASetMaxStep(mem, time_step); + checkOutput(retval, "IDASetMaxStep"); + + /* Since the starting procedure is first order, the maximum global order + * of convergence is two */ + retval = IDASetMaxOrd(mem, 2); + checkOutput(retval, "IDASetMaxOrd"); + + /* Enable more nonlinear iterations because a failed nonlinear solve + * causes a failed integration with fixed steps */ + retval = IDASetMaxNonlinIters(mem, 100); + checkOutput(retval, "IDASetMaxNonlinIters"); + + // Set a large tolerance so the error test will never fail + static constexpr RealT FIXED_STEP_TOL_FAC = 1e100; + setTolerance(mem, + FIXED_STEP_TOL_FAC * rel_tol, + FIXED_STEP_TOL_FAC * abs_tol_override, + FIXED_STEP_TOL_FAC); + + /* We want the nonlinear solver tolerance to be ~rel_tol, but the with + * the large tolerances set above, we need to choose this tolerance to + * "undo" the fac scaling. */ + retval = IDASetNonlinConvCoef(mem, 1 / FIXED_STEP_TOL_FAC); + checkOutput(retval, "IDASetNonlinConvCoef"); + } + + /** + * @brief Set the relative tolerance and optionally override the absolute + * tolerance + * + * @param rel_tol The relative tolerance to use + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance rather than the model's default absolute + * tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setTolerance(ScalarT rel_tol, + ScalarT abs_tol_override) + { + setTolerance(solver_, rel_tol, abs_tol_override); + } + + /** + * @brief Set the relative tolerance and optionally override the absolute + * tolerance for the backward simulation + * + * @param rel_tol The relative tolerance to use + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance rather than the model's default absolute + * tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setBackwardTolerance(ScalarT rel_tol, + ScalarT abs_tol_override) + { + setTolerance(IDAGetAdjIDABmem(solver_, backwardID_), + rel_tol, + abs_tol_override); + } + + /** + * @brief A helper function to set the relative tolerance and optionally + * override the absolute tolerance + * + * @param mem The IDA memory (either forward or backward) + * @param rel_tol The relative tolerance to use + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance rather than the model's default absolute + * tolerance + * @param abs_tol_fac A factor to apply to the absolute tolerance if not + * overridden + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setTolerance(void* mem, + ScalarT rel_tol, + ScalarT abs_tol_override, + ScalarT abs_tol_fac) + { + int retval = 0; + + if (abs_tol_override > 0) + { + retval = IDASStolerances(mem, rel_tol, abs_tol_override); + checkOutput(retval, "IDASStolerances"); + return; + } + + model_->setAbsoluteTolerance(rel_tol); + const std::vector& abs_tol = model_->absoluteTolerance(); + copyVec(abs_tol, abs_tol_); + N_VScale(abs_tol_fac, abs_tol_, abs_tol_); + retval = IDASVtolerances(mem, rel_tol, abs_tol_); + checkOutput(retval, "IDASVtolerances"); + } + + /** + * @brief Set the maximum number of steps + * + * @param max_steps The maximum number of steps + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setMaxSteps(IdxT max_steps) + { + setMaxSteps(solver_, max_steps); + } + + /** + * @brief Set the maximum number of steps for the backward simulation + * + * @param max_steps The maximum number of steps + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setBackwardMaxSteps(IdxT max_steps) + { + setMaxSteps(IDAGetAdjIDABmem(solver_, backwardID_), max_steps); + } + + /** + * @brief A helper function to set the maximum number of steps + * + * @param mem The IDA memory (either forward or backward) + * @param max_steps The maximum number of steps + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setMaxSteps(void* mem, IdxT max_steps) + { + int retval = IDASetMaxNumSteps(mem, static_cast(max_steps)); + checkOutput(retval, "IDASetMaxNumSteps"); + } + + /** + * @brief Set the quadrature relative tolerance and optionally override the + * absolute tolerance + * + * @param rel_tol The relative tolerance to use + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance rather than the model's default absolute + * tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setQuadratureTolerance(ScalarT rel_tol, + ScalarT abs_tol_override) + { + setQuadratureTolerance(solver_, rel_tol, abs_tol_override); + } + + /** + * @brief Set the quadrature relative tolerance and optionally override the + * absolute tolerance for the backward simulation + * + * @param rel_tol The relative tolerance to use + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance rather than the model's default absolute + * tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setBackwardQuadratureTolerance(ScalarT rel_tol, + ScalarT abs_tol_override) + { + setQuadratureTolerance(IDAGetAdjIDABmem(solver_, backwardID_), + rel_tol, + abs_tol_override); + } + + /** + * @brief A helper function to set the quadrature relative tolerance and + * optionally override the absolute tolerance + * + * @param mem The IDA memory (either forward or backward) + * @param rel_tol The relative tolerance to use + * @param abs_tol_override If positive, this value will be used as the + * absolute tolerance rather than the model's default absolute + * tolerance + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setQuadratureTolerance(void* mem, + ScalarT rel_tol, + ScalarT abs_tol_override) + { + int retval = 0; + + if (abs_tol_override > 0) + { + retval = IDAQuadSStolerances(mem, rel_tol, abs_tol_override); + checkOutput(retval, "IDAQuadSStolerances"); + return; + } + + model_->setAbsoluteTolerance(rel_tol); + const std::vector& abs_tol = model_->absoluteTolerance(); + if (static_cast(abs_tol.size()) != model_->size()) + { + std::cerr << "\nModel absolute tolerance does not match model size.\n\n"; + throw SundialsException(); + } + copyVec(abs_tol, abs_tol_); + retval = IDAQuadSVtolerances(mem, rel_tol, abs_tol_); + checkOutput(retval, "IDAQuadSVtolerances"); + } + // Compiler will prevent building modules with data type incompatible with sunrealtype template class Ida; template class Ida; diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 20ff56aa8..84600e680 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -59,7 +59,7 @@ namespace AnalysisManager int configureLinearSolverDense(); int getDefaultInitialCondition(); int setIntegrationTime(RealT t_init, RealT t_final, int nout); - int initializeSimulation(RealT t0, bool findConsistent = false); + int initializeSimulation(RealT t0, bool findConsistent = true); int runSimulation(RealT tf, int nout = 1, std::optional> step_callback = {}); int deleteSimulation(); @@ -127,9 +127,25 @@ namespace AnalysisManager return N_VGetArrayPointer(qB_); } - void printOutput(RealT t); - void printSpecial(RealT t, N_Vector x); - void printFinalStats(); + void printOutput(RealT t) const; + void printSpecial(RealT t, N_Vector x) const; + void printFinalStats() const; + + void setFixedStep(ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override = 0); + void setBackwardFixedStep(ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override = 0); + using DynamicSolver::setTolerance; + void setTolerance(ScalarT rel_tol, ScalarT abs_tol_override) override; + void setBackwardTolerance(ScalarT rel_tol, ScalarT abs_tol_override = 0); + void setQuadratureTolerance(ScalarT rel_tol, + ScalarT abs_tol_override = 0); + void setBackwardQuadratureTolerance(ScalarT rel_tol, + ScalarT abs_tol_override = 0); + void setMaxSteps(IdxT maxSteps) override; + void setBackwardMaxSteps(IdxT maxSteps); IdaStats getStats() const; @@ -174,6 +190,8 @@ namespace AnalysisManager void* user_data); private: + static constexpr ScalarT DEFAULT_REL_TOL = 1e-5; + void* solver_{}; SUNContext context_{}; SUNMatrix JacobianMat_{}; @@ -185,10 +203,11 @@ namespace AnalysisManager RealT t_final_{}; int nout_{}; ///< Number of integration outputs - N_Vector yy_{}; ///< Solution vector - N_Vector yp_{}; ///< Solution derivatives vector - N_Vector tag_{}; ///< Tags differential variables - N_Vector q_{}; ///< Integrand vector + N_Vector yy_{}; ///< Solution vector + N_Vector yp_{}; ///< Solution derivatives vector + N_Vector tag_{}; ///< Tags differential variables + N_Vector abs_tol_{}; ///< Absolute tolerance vector + N_Vector q_{}; ///< Integrand vector N_Vector yy0_{}; ///< Storage for initial values N_Vector yp0_{}; ///< Storage for initial derivatives @@ -208,6 +227,19 @@ namespace AnalysisManager // int check_flag(void *flagvalue, const char *funcname, int opt); static void checkAllocation(void* v, const char* functionName); static void checkOutput(int retval, const char* functionName); + + void setFixedStep(void* mem, + ScalarT time_step, + ScalarT rel_tol, + ScalarT abs_tol_override); + void setTolerance(void* mem, + ScalarT rel_tol, + ScalarT abs_tol_override, + ScalarT abs_tol_fac = 1); + void setMaxSteps(void* mem, IdxT max_steps); + void setQuadratureTolerance(void* mem, + ScalarT rel_tol, + ScalarT abs_tol_override); }; /// Simple exception to use within Ida class. diff --git a/GridKit/Solver/SteadyState/Kinsol.cpp b/GridKit/Solver/SteadyState/Kinsol.cpp index 7e305318f..aa04fe4e6 100644 --- a/GridKit/Solver/SteadyState/Kinsol.cpp +++ b/GridKit/Solver/SteadyState/Kinsol.cpp @@ -72,17 +72,6 @@ namespace AnalysisManager retval = KINSetUserData(solver_, model_); checkOutput(retval, "KINSetUserData"); - // Set tolerances - sunrealtype fnormtol; ///< Residual tolerance - sunrealtype scsteptol; ///< Scaled step tolerance - - model_->setTolerances(fnormtol, scsteptol); ///< \todo Function name should be "getTolerances"! - retval = KINSetFuncNormTol(solver_, fnormtol); - checkOutput(retval, "KINSetFuncNormTol"); - - retval = KINSetScaledStepTol(solver_, scsteptol); - checkOutput(retval, "KINSetScaledStepTol"); - // Set up linear solver return this->configureLinearSolver(); } @@ -171,7 +160,7 @@ namespace AnalysisManager } template - void Kinsol::printOutput() + void Kinsol::printOutput() const { sunrealtype* yval = N_VGetArrayPointer(yy_); @@ -184,7 +173,7 @@ namespace AnalysisManager } template - void Kinsol::printSpecial(sunrealtype t, N_Vector y) + void Kinsol::printSpecial(sunrealtype t, N_Vector y) const { sunrealtype* yval = N_VGetArrayPointer_Serial(y); IdxT N = static_cast(N_VGetLength_Serial(y)); @@ -198,10 +187,10 @@ namespace AnalysisManager } template - void Kinsol::printFinalStats() + void Kinsol::printFinalStats() const { int retval = KINPrintAllStats(solver_, stdout, SUN_OUTPUTFORMAT_TABLE); - checkOutput(retval, "IDAPrintAllStats"); + checkOutput(retval, "KINPrintAllStats"); } template @@ -224,6 +213,16 @@ namespace AnalysisManager } } + template + void Kinsol::setTolerance(ScalarT tol) + { + int retval = KINSetFuncNormTol(solver_, tol); + checkOutput(retval, "KINSetFuncNormTol"); + + retval = KINSetScaledStepTol(solver_, tol); + checkOutput(retval, "KINSetScaledStepTol"); + } + // Compiler will prevent building modules with data type incompatible with sunrealtype template class Kinsol; template class Kinsol; diff --git a/GridKit/Solver/SteadyState/Kinsol.hpp b/GridKit/Solver/SteadyState/Kinsol.hpp index 3b4b5c405..89fc6c720 100644 --- a/GridKit/Solver/SteadyState/Kinsol.hpp +++ b/GridKit/Solver/SteadyState/Kinsol.hpp @@ -105,9 +105,9 @@ namespace AnalysisManager // return N_VGetArrayPointer(qB_); // } - void printOutput(); - void printSpecial(sunrealtype t, N_Vector x); - void printFinalStats(); + void printOutput() const; + void printSpecial(sunrealtype t, N_Vector x) const; + void printFinalStats() const; private: static int Residual(N_Vector yy, N_Vector rr, void* user_data); @@ -145,8 +145,10 @@ namespace AnalysisManager static void copyVec(const std::vector& x, N_Vector y); // int check_flag(void *flagvalue, const char *funcname, int opt); - inline void checkAllocation(void* v, const char* functionName); - inline void checkOutput(int retval, const char* functionName); + static inline void checkAllocation(void* v, const char* functionName); + static inline void checkOutput(int retval, const char* functionName); + + void setTolerance(ScalarT tol) override; }; /// Simple exception to use within Kinsol class. diff --git a/GridKit/Solver/SteadyState/SteadyStateSolver.hpp b/GridKit/Solver/SteadyState/SteadyStateSolver.hpp index 05fa9da55..19643e331 100644 --- a/GridKit/Solver/SteadyState/SteadyStateSolver.hpp +++ b/GridKit/Solver/SteadyState/SteadyStateSolver.hpp @@ -22,6 +22,8 @@ namespace AnalysisManager return model_; } + virtual void setTolerance(ScalarT tol) = 0; + protected: GridKit::Model::Evaluator* model_; }; diff --git a/application/PhasorDynamics/DynamicSimulation.cpp b/application/PhasorDynamics/DynamicSimulation.cpp index 3aaadd21f..72861ea90 100644 --- a/application/PhasorDynamics/DynamicSimulation.cpp +++ b/application/PhasorDynamics/DynamicSimulation.cpp @@ -39,7 +39,7 @@ int main(int argc, const char* argv[]) real_type dt = study.dt; real_type final_time = study.tmax; real_type curr_time = 0.0; - ida.initializeSimulation(0.0, false); + ida.initializeSimulation(0.0); for (const auto& event : study.events) { // Run to event time @@ -58,7 +58,7 @@ int main(int argc, const char* argv[]) } // Re-initialize simulation at event time - ida.initializeSimulation(event.time, true); + ida.initializeSimulation(event.time); curr_time = event.time; } diff --git a/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp b/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp index ef4181079..1d67cf0ec 100644 --- a/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp +++ b/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp @@ -59,8 +59,8 @@ int main() // create initial condition after a fault { idas->getSavedInitialCondition(); - idas->initializeSimulation(t_init); gen->V() = 0.0; + idas->initializeSimulation(t_init); idas->runSimulation(0.1, 20); gen->V() = 1.0; idas->saveInitialCondition(); diff --git a/examples/Experimental/GenConstLoad/GenConstLoad.cpp b/examples/Experimental/GenConstLoad/GenConstLoad.cpp index d2d01d76a..9b3a7ed11 100644 --- a/examples/Experimental/GenConstLoad/GenConstLoad.cpp +++ b/examples/Experimental/GenConstLoad/GenConstLoad.cpp @@ -48,7 +48,7 @@ int main() idas->configureSimulation(); idas->configureAdjoint(); idas->getDefaultInitialCondition(); - idas->initializeSimulation(t_init, true); + idas->initializeSimulation(t_init); idas->configureQuadrature(); idas->initializeQuadrature(); @@ -80,7 +80,7 @@ int main() // Configure Ipopt application ipoptApp->Options()->SetStringValue("hessian_approximation", "limited-memory"); ipoptApp->Options()->SetNumericValue("tol", tol); - ipoptApp->Options()->SetIntegerValue("print_level", 5); + ipoptApp->Options()->SetIntegerValue("print_level", 0); // Create interface to Ipopt solver Ipopt::SmartPtr ipoptDynamicObjectiveInterface = diff --git a/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp b/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp index 48c2650d1..351af82de 100644 --- a/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp +++ b/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp @@ -55,8 +55,8 @@ int main() // create initial condition after a fault { idas.getSavedInitialCondition(); - idas.initializeSimulation(t_init); gen.V() = 0.0; + idas.initializeSimulation(t_init); idas.runSimulation(t_clear, 20); gen.V() = 1.0; idas.saveInitialCondition(); diff --git a/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp b/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp index cb7c2aff9..95a751a92 100644 --- a/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp +++ b/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp @@ -65,8 +65,8 @@ int main() // create initial condition after a fault { idas.getSavedInitialCondition(); - idas.initializeSimulation(t_init); gen.V() = 0.0; + idas.initializeSimulation(t_init); idas.runSimulation(t_clear, 20); gen.V() = 1.0; idas.saveInitialCondition(); diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp index 16e648a71..4f687a5bc 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp @@ -170,13 +170,13 @@ int main() // Introduce fault to ground and run for 0.1s fault.setStatus(1); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear fault and run until t = 10s. fault.setStatus(0); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp index 5f7b05e8b..e75edf180 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp @@ -156,13 +156,13 @@ int main() // Introduce fault to ground and run for 0.1s fault.setStatus(1); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear fault and run until t = 10s. fault.setStatus(0); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp index 326760e80..36acdbe0b 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp @@ -232,7 +232,7 @@ int main() // Run simulation, output each `dt` interval real_type start = static_cast(clock()); - ida.initializeSimulation(0.0, false); + ida.initializeSimulation(0.0); // Run for 1s int nout = static_cast(std::round((1.0 - 0.0) / dt)); @@ -240,13 +240,13 @@ int main() // Introduce fault to ground and run for 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp index d302d6e23..fe00fe68e 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp @@ -81,7 +81,7 @@ int main(int argc, const char* argv[]) // Run simulation, output each `dt` interval real_type start = static_cast(clock()); - ida.initializeSimulation(0.0, false); + ida.initializeSimulation(0.0); // Run for 1s int nout = static_cast(std::round((1.0 - 0.0) / dt)); @@ -89,13 +89,13 @@ int main(int argc, const char* argv[]) // Introduce fault to ground and run for 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout); // Clear fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp index 4373ef46c..cd3d12aa6 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp @@ -202,7 +202,7 @@ int main() // Run simulation, output each `dt` interval real_type start = static_cast(clock()); - ida.initializeSimulation(0.0, false); + ida.initializeSimulation(0.0); // Run for 1s int nout = static_cast(std::round((1.0 - 0.0) / dt)); @@ -210,13 +210,13 @@ int main() // Introduce fault to ground and run for 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp index aa7af92c4..4892d1b1b 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp @@ -149,7 +149,7 @@ int main(int argc, const char* argv[]) // Run simulation, output each `dt` interval real_type start = static_cast(clock()); - ida.initializeSimulation(0.0, false); + ida.initializeSimulation(0.0); // Run for 1s int nout = static_cast(std::round((1.0 - 0.0) / dt)); @@ -157,13 +157,13 @@ int main(int argc, const char* argv[]) // Introduce fault to ground and run for 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp index 43fe05c59..00d3fb551 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp @@ -155,13 +155,13 @@ int main() // Introduce fault and run for the next 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear the fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.solver.json b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.solver.json index ef55cd111..2f9b9d974 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.solver.json +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.solver.json @@ -7,5 +7,5 @@ {"time": 1.1, "type": "fault_off", "element_id": 0} ], "reference_file": "TwoBusBasic.ref.csv", - "error_tolerance": 2.01e-4 + "error_tolerance": 2.03e-4 } diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp index b1f5cbd06..1efb9ba3f 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp @@ -135,13 +135,13 @@ int main(int argc, const char* argv[]) // Introduce fault and run for the next 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear the fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp index 392a80dc6..3027075d5 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp @@ -211,13 +211,13 @@ int main() // Introduce fault and run for the next 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear the fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp index 623395822..692e9431a 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp @@ -145,14 +145,14 @@ int main(int argc, const char* argv[]) // Introduce fault and run for the next 0.1s // fault.setStatus(true); fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear the fault and run until t = 10s. // fault.setStatus(false); fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp index b50164656..a7dccb481 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp @@ -176,13 +176,13 @@ int main() // Introduce fault and run for the next 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear the fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp index 27aefa1fd..55d2ebd14 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp @@ -139,13 +139,13 @@ int main(int argc, const char* argv[]) // Introduce fault and run for the next 0.1s fault->setStatus(true); - ida.initializeSimulation(1.0, false); + ida.initializeSimulation(1.0); nout = static_cast(std::round((1.1 - 1.0) / dt)); ida.runSimulation(1.1, nout, output_cb); // Clear the fault and run until t = 10s. fault->setStatus(false); - ida.initializeSimulation(1.1, false); + ida.initializeSimulation(1.1); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); real_type stop = static_cast(clock()); diff --git a/examples/PowerElectronics/Microgrid/Microgrid.cpp b/examples/PowerElectronics/Microgrid/Microgrid.cpp index 3b20c14e0..7c6511c56 100644 --- a/examples/PowerElectronics/Microgrid/Microgrid.cpp +++ b/examples/PowerElectronics/Microgrid/Microgrid.cpp @@ -20,14 +20,13 @@ int main(int /* argc */, char const** /* argv */) { /// @todo Needs to be modified. Some components are small relative to others thus /// there error is high (or could be matlab vector issue) - double abs_tol = 1.0e-8; double rel_tol = 1.0e-8; size_t max_step_number = 3000; bool use_jac = true; bool debug_output = true; // Create model - auto* sysmodel = new GridKit::PowerElectronicsModel(rel_tol, abs_tol, use_jac, max_step_number); + auto* sysmodel = new GridKit::PowerElectronicsModel(use_jac); // Modeled after the problem in the paper double RN = 1.0e4; @@ -223,6 +222,8 @@ int main(int /* argc */, char const** /* argv */) // setup simulation idas->configureSimulation(); + idas->setTolerance(rel_tol); + idas->setMaxSteps(max_step_number); idas->getDefaultInitialCondition(); idas->initializeSimulation(t_init); diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index 5e21be567..924d99fb5 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -25,7 +25,7 @@ int main(int /* argc */, char const** /* argv */) // TODO:setup as named parameters // Create circuit model - GridKit::PowerElectronicsModel sysmodel(rel_tol, abs_tol, use_jac); + GridKit::PowerElectronicsModel sysmodel(use_jac); size_t idoff = 0; @@ -98,6 +98,7 @@ int main(int /* argc */, char const** /* argv */) // setup simulation idas.configureSimulation(); + idas.setTolerance(rel_tol, abs_tol); idas.getDefaultInitialCondition(); idas.initializeSimulation(t_init); diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp index 07246e83a..1befe7bc0 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -74,10 +74,7 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) real_type abs_tol = SCALE_MICROGRID_ABS_TOL; // Create circuit model - auto* sys_model = new PowerElectronicsModel(rel_tol, - abs_tol, - use_jac, - SCALE_MICROGRID_MAX_STEPS); + auto* sys_model = new PowerElectronicsModel(use_jac); const std::vector* true_vec = &answer_key_N8; @@ -290,6 +287,8 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) // setup simulation idas->configureSimulation(); + idas->setTolerance(rel_tol, abs_tol); + idas->setMaxSteps(SCALE_MICROGRID_MAX_STEPS); idas->getDefaultInitialCondition(); idas->initializeSimulation(t_init); diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp index b18846963..5c8d1d37f 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp @@ -68,13 +68,9 @@ int printMicrogridSystems(index_type N_size) real_type t_final = 1.0; real_type rel_tol = SCALE_MICROGRID_REL_TOL; - real_type abs_tol = SCALE_MICROGRID_ABS_TOL; // Create circuit model - PowerElectronicsModel sys_model(rel_tol, - abs_tol, - use_jac, - SCALE_MICROGRID_MAX_STEPS); + PowerElectronicsModel sys_model(use_jac); // Ensure minimum size requirement if (N_size < 1) @@ -260,6 +256,8 @@ int printMicrogridSystems(index_type N_size) // setup simulation idas.configureSimulation(); + idas.setTolerance(rel_tol); + idas.setMaxSteps(SCALE_MICROGRID_MAX_STEPS); idas.getDefaultInitialCondition(); idas.initializeSimulation(t_init); idas.runSimulation(t_final); diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index e3592bcf2..5a2f57a75 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -29,7 +29,8 @@ namespace GridKit y_ = {0}; yp_ = {0}; - tag_ = {false}; + tag_ = {false}; + abs_tol_ = {0}; f_ = {0}; g_ = {0}; @@ -61,17 +62,14 @@ namespace GridKit return 0; } - void setTolerances([[maybe_unused]] RealT& rel_tol, [[maybe_unused]] RealT& abs_tol) const override - { - } - - void setMaxSteps(IdxT& msa) const override + int tagDifferentiable() override { - msa = 2000; + return 0; } - int tagDifferentiable() override + int setAbsoluteTolerance(RealT rel_tol) override { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); return 0; } @@ -140,6 +138,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() override + { + return abs_tol_; + } + + const std::vector& absoluteTolerance() const override + { + return abs_tol_; + } + std::vector& yB() override { return yB_; @@ -249,6 +257,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector abs_tol_; std::vector f_; std::vector g_; @@ -271,7 +280,7 @@ namespace GridKit class IdaTests { public: - TestOutcome test() + TestOutcome callback() { const unsigned n_steps = 100; TestStatus success = true; @@ -294,6 +303,26 @@ namespace GridKit return success.report(__func__); } + + TestOutcome fixedStep() + { + const unsigned n_steps = 32; + TestStatus success = true; + + Model::NullEvaluator model; + + Ida ida(&model); + ida.configureSimulation(); + ida.setFixedStep(1.0 / n_steps, 1e-6); + + ida.initializeSimulation(0.0, false); + ida.runSimulation(1.0); + auto stats = ida.getStats(); + + success *= (stats.num_steps_ == n_steps); + + return success.report(__func__); + } }; } // namespace Testing } // namespace GridKit diff --git a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp index 5eebb666f..bb3853bd2 100644 --- a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp +++ b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp @@ -8,7 +8,8 @@ int main() GridKit::Testing::TestingResults result; GridKit::Testing::IdaTests test; - result += test.test(); + result += test.callback(); + result += test.fixedStep(); return result.summary(); }