From de7728679b5ccb666043b756a47c11599da69f46 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Thu, 8 Jan 2026 16:57:46 -0800 Subject: [PATCH 1/6] Add fixed time step for IDA forward simulation --- GridKit/Solver/Dynamic/Ida.cpp | 40 +++++++++++++++++++++++++++++----- GridKit/Solver/Dynamic/Ida.hpp | 4 +++- 2 files changed, 38 insertions(+), 6 deletions(-) diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index baa303a5..9b88bbea 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -16,8 +16,8 @@ namespace AnalysisManager { template - Ida::Ida(GridKit::Model::Evaluator* model) - : DynamicSolver(model) + Ida::Ida(GridKit::Model::Evaluator* model, RealT dt) + : DynamicSolver(model), dt_(dt) { int retval = 0; @@ -81,14 +81,44 @@ namespace AnalysisManager // Set pointer to model data 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"); + + if (dt_ > 0) { + retval = IDASetMinStep(solver_, dt_); + checkOutput(retval, "IDASetMinStep"); + retval = IDASetMaxStep(solver_, dt_); + checkOutput(retval, "IDASetMaxStep"); + + /* Since the starting procedure is first order, the maximum global order + * of convergence is two */ + IDASetMaxOrd(solver_, 2); + checkOutput(retval, "IDASetMaxOrd"); + + + /* Enable more nonlinear iterations because a failed nonlinear solve + * causes a failed integration with fixed steps */ + retval = IDASetMaxNonlinIters(solver_, 100); + checkOutput(retval, "IDASetMaxNonlinIters"); + + // A large number so that the error test will never fail + static constexpr RealT fac = 1e10; + retval = IDASStolerances(solver_, fac, fac * abs_tol / rel_tol); + checkOutput(retval, "IDASStolerances"); + + /* 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(solver_, rel_tol / fac); + checkOutput(retval, "IDASetNonlinConvCoef"); + } else { + retval = IDASStolerances(solver_, rel_tol, abs_tol); + checkOutput(retval, "IDASStolerances"); + } IdxT msa; model_->setMaxSteps(msa); diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 8e1715e5..dcc91adc 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -32,7 +32,7 @@ namespace AnalysisManager using RealT = typename GridKit::ScalarTraits::RealT; public: - Ida(GridKit::Model::Evaluator* model); + Ida(GridKit::Model::Evaluator* model, RealT dt=0); ~Ida(); int configureSimulation(); @@ -156,6 +156,8 @@ namespace AnalysisManager void* user_data); private: + RealT dt_; + void* solver_{}; SUNContext context_{}; SUNMatrix JacobianMat_{}; From d6e668c565e8043e404e24f9a15e4ac279e781d8 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 9 Jan 2026 13:40:25 -0800 Subject: [PATCH 2/6] Add fixed steps to adjoint --- GridKit/Solver/Dynamic/Ida.cpp | 112 ++++++++++++++++++--------------- GridKit/Solver/Dynamic/Ida.hpp | 2 + 2 files changed, 63 insertions(+), 51 deletions(-) diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index 9b88bbea..ea1900aa 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -81,51 +81,8 @@ namespace AnalysisManager // Set pointer to model data 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"! - - if (dt_ > 0) { - retval = IDASetMinStep(solver_, dt_); - checkOutput(retval, "IDASetMinStep"); - retval = IDASetMaxStep(solver_, dt_); - checkOutput(retval, "IDASetMaxStep"); - - /* Since the starting procedure is first order, the maximum global order - * of convergence is two */ - IDASetMaxOrd(solver_, 2); - checkOutput(retval, "IDASetMaxOrd"); - - - /* Enable more nonlinear iterations because a failed nonlinear solve - * causes a failed integration with fixed steps */ - retval = IDASetMaxNonlinIters(solver_, 100); - checkOutput(retval, "IDASetMaxNonlinIters"); - - // A large number so that the error test will never fail - static constexpr RealT fac = 1e10; - retval = IDASStolerances(solver_, fac, fac * abs_tol / rel_tol); - checkOutput(retval, "IDASStolerances"); - /* 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(solver_, rel_tol / fac); - checkOutput(retval, "IDASetNonlinConvCoef"); - } else { - 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"); + configureIDA(solver_); // Tag differential variables std::vector& tag = model_->tag(); @@ -185,6 +142,64 @@ namespace AnalysisManager return retval; } + /** + * @brief Configure IDA for forward or backward simulation + * + * @tparam ScalarT + * @tparam IdxT + */ + template + void Ida::configureIDA(void *ida) + { + int retval = 0; + + // Set tolerances + RealT rel_tol; + RealT abs_tol; + + model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"! + + if (dt_ > 0) { + retval = IDASetMinStep(ida, dt_); + checkOutput(retval, "IDASetMinStep"); + retval = IDASetMaxStep(ida, dt_); + checkOutput(retval, "IDASetMaxStep"); + + /* Since the starting procedure is first order, the maximum global order + * of convergence is two */ + IDASetMaxOrd(ida, 2); + checkOutput(retval, "IDASetMaxOrd"); + + + /* Enable more nonlinear iterations because a failed nonlinear solve + * causes a failed integration with fixed steps */ + retval = IDASetMaxNonlinIters(ida, 100); + checkOutput(retval, "IDASetMaxNonlinIters"); + + // Set a large tolerance so the error test will never fail + static constexpr RealT FIXED_STEP_TOL_FAC = 1e10; + retval = IDASStolerances(ida, FIXED_STEP_TOL_FAC, + FIXED_STEP_TOL_FAC * abs_tol / rel_tol); + checkOutput(retval, "IDASStolerances"); + + /* 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(ida, rel_tol / FIXED_STEP_TOL_FAC); + checkOutput(retval, "IDASetNonlinConvCoef"); + } else { + retval = IDASStolerances(ida, 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(ida, static_cast(msa)); + checkOutput(retval, "IDASetMaxNumSteps"); + } + #ifdef GRIDKIT_ENABLE_SUNDIALS_SPARSE /** * @brief Configure a sparse linear solver @@ -589,16 +604,10 @@ 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"); - 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"); + configureIDA(IDAGetAdjIDABmem(solver_, backwardID_)); // Allocate Jacobian matrix, if not already if (JacobianMatB_ == nullptr) @@ -625,6 +634,7 @@ namespace AnalysisManager checkOutput(retval, "IDAQuadInitB"); // retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol*1.1, abs_tol*1.1); + model_->setTolerances(rel_tol, abs_tol); retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol * 0.1, abs_tol * 0.1); checkOutput(retval, "IDAQuadSStolerancesB"); diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index dcc91adc..2b424048 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -184,6 +184,8 @@ namespace AnalysisManager int backwardID_{}; private: + void configureIDA(void *ida); + // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); static void copyVec(const N_Vector x, std::vector& y); static void copyVec(const std::vector& x, N_Vector y); From 391c6d0f903b269ef51eabb064481c0dcaee9f11 Mon Sep 17 00:00:00 2001 From: Steven-Roberts Date: Fri, 9 Jan 2026 22:16:39 +0000 Subject: [PATCH 3/6] Apply pre-commmit fixes --- GridKit/Solver/Dynamic/Ida.cpp | 19 ++++++++++--------- GridKit/Solver/Dynamic/Ida.hpp | 4 ++-- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index ea1900aa..f80b76d2 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -149,28 +149,28 @@ namespace AnalysisManager * @tparam IdxT */ template - void Ida::configureIDA(void *ida) + void Ida::configureIDA(void* ida) { int retval = 0; - + // Set tolerances RealT rel_tol; RealT abs_tol; model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"! - - if (dt_ > 0) { + + if (dt_ > 0) + { retval = IDASetMinStep(ida, dt_); checkOutput(retval, "IDASetMinStep"); retval = IDASetMaxStep(ida, dt_); checkOutput(retval, "IDASetMaxStep"); - + /* Since the starting procedure is first order, the maximum global order * of convergence is two */ IDASetMaxOrd(ida, 2); checkOutput(retval, "IDASetMaxOrd"); - /* Enable more nonlinear iterations because a failed nonlinear solve * causes a failed integration with fixed steps */ retval = IDASetMaxNonlinIters(ida, 100); @@ -178,8 +178,7 @@ namespace AnalysisManager // Set a large tolerance so the error test will never fail static constexpr RealT FIXED_STEP_TOL_FAC = 1e10; - retval = IDASStolerances(ida, FIXED_STEP_TOL_FAC, - FIXED_STEP_TOL_FAC * abs_tol / rel_tol); + retval = IDASStolerances(ida, FIXED_STEP_TOL_FAC, FIXED_STEP_TOL_FAC * abs_tol / rel_tol); checkOutput(retval, "IDASStolerances"); /* We want the nonlinear solver tolerance to be ~rel_tol, but the with @@ -187,7 +186,9 @@ namespace AnalysisManager * "undo" the fac scaling. */ retval = IDASetNonlinConvCoef(ida, rel_tol / FIXED_STEP_TOL_FAC); checkOutput(retval, "IDASetNonlinConvCoef"); - } else { + } + else + { retval = IDASStolerances(ida, rel_tol, abs_tol); checkOutput(retval, "IDASStolerances"); } diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 2b424048..722dcc3f 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -32,7 +32,7 @@ namespace AnalysisManager using RealT = typename GridKit::ScalarTraits::RealT; public: - Ida(GridKit::Model::Evaluator* model, RealT dt=0); + Ida(GridKit::Model::Evaluator* model, RealT dt = 0); ~Ida(); int configureSimulation(); @@ -184,7 +184,7 @@ namespace AnalysisManager int backwardID_{}; private: - void configureIDA(void *ida); + void configureIDA(void* ida); // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); static void copyVec(const N_Vector x, std::vector& y); From 92a40204bf05480bc9bc401c27e34705d53bcab5 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 30 Jan 2026 18:11:40 -0800 Subject: [PATCH 4/6] Add abstoluteTolerance function to Evaluator --- GridKit/Model/Evaluator.hpp | 5 +- GridKit/Model/PhasorDynamics/BusBase.hpp | 22 +- GridKit/Model/PhasorDynamics/Component.hpp | 12 +- .../PhasorDynamics/SignalNode/SignalNode.hpp | 22 +- GridKit/Model/PhasorDynamics/SystemModel.hpp | 14 +- .../PowerElectronics/CircuitComponent.hpp | 20 +- .../SystemModelPowerElectronics.hpp | 18 +- GridKit/Model/PowerFlow/Bus/BaseBus.hpp | 2 - GridKit/Model/PowerFlow/Bus/BusSlack.hpp | 2 - GridKit/Model/PowerFlow/MiniGrid/MiniGrid.cpp | 2 - GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp | 2 - .../Model/PowerFlow/ModelEvaluatorImpl.hpp | 20 +- GridKit/Model/PowerFlow/SystemModel.hpp | 9 +- .../Model/PowerFlow/SystemModelPowerFlow.hpp | 11 +- GridKit/Solver/Dynamic/DynamicSolver.hpp | 4 + GridKit/Solver/Dynamic/Ida.cpp | 218 +++++++++++------- GridKit/Solver/Dynamic/Ida.hpp | 40 +++- GridKit/Solver/SteadyState/Kinsol.cpp | 28 ++- GridKit/Solver/SteadyState/Kinsol.hpp | 13 +- .../Solver/SteadyState/SteadyStateSolver.hpp | 2 + .../PowerElectronics/Microgrid/Microgrid.cpp | 5 +- .../PowerElectronics/RLCircuit/RLCircuit.cpp | 4 +- .../ScaleMicrogrid/ScaleMicrogrid.cpp | 8 +- .../ScaleMicrogridArbitrary.cpp | 8 +- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 20 +- 25 files changed, 271 insertions(+), 240 deletions(-) diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index 9bc9f010..e21c838f 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -94,8 +94,9 @@ namespace GridKit 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 std::vector& absoluteTolerance() = 0; + 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/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 603f9699..a59f5f75 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -58,17 +58,6 @@ namespace GridKit // No time to update in bus models } - virtual void setTolerances(RealT& rtol, RealT& atol) const override - { - rtol = rtol_; - atol = atol_; - } - - virtual void setMaxSteps(IdxT& msa) const override - { - msa = max_steps_; - } - virtual ScalarT& Vr() = 0; virtual const ScalarT& Vr() const = 0; virtual ScalarT& Vi() = 0; @@ -108,6 +97,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() override + { + return absTol_; + } + + const std::vector& absoluteTolerance() const override + { + return absTol_; + } + MatrixT& getJacobian() override { return J_; @@ -175,6 +174,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector absTol_; std::vector f_; MatrixT J_; diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 7fdf3594..e0b2f610 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -49,15 +49,14 @@ namespace GridKit // std::cout << "updateTime: t = " << time_ << "\n"; // } - virtual void setTolerances(RealT& rtol, RealT& atol) const override + std::vector& absoluteTolerance() override { - rtol = rel_tol_; - atol = abs_tol_; + return absTol_; } - virtual void setMaxSteps(IdxT& msa) const override + const std::vector& absoluteTolerance() const override { - msa = max_steps_; + return absTol_; } std::vector& y() override @@ -161,6 +160,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector absTol_; std::vector f_; std::vector g_; std::vector wb_; @@ -174,8 +174,6 @@ namespace GridKit RealT rel_tol_; RealT abs_tol_; - IdxT max_steps_; - /* ------ WARNING: Temporary ------ diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp index 7bd6e91b..f2226df3 100644 --- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp +++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp @@ -70,17 +70,6 @@ namespace GridKit // No time to update in bus models } - virtual void setTolerances(RealT& rtol, RealT& atol) const override - { - rtol = rtol_; - atol = atol_; - } - - virtual void setMaxSteps(IdxT& msa) const override - { - msa = max_steps_; - } - std::vector& y() override { return y_; @@ -111,6 +100,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() override + { + return absTol_; + } + + const std::vector& absoluteTolerance() const override + { + return absTol_; + } + MatrixT& getJacobian() override { return J_; @@ -140,6 +139,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector absTol_; std::vector f_; MatrixT J_; diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 2d787a02..223cd837 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -49,20 +49,13 @@ namespace GridKit using PhasorDynamics::Component::tag_; using PhasorDynamics::Component::f_; using PhasorDynamics::Component::J_; - using PhasorDynamics::Component::rel_tol_; - using PhasorDynamics::Component::abs_tol_; public: /** * @brief Constructor for the system model */ SystemModel() - { - // Set system model tolerances - rel_tol_ = 1e-7; - abs_tol_ = 1e-9; - this->max_steps_ = 2000; - } + {} /** * @brief Construct a new System Model object @@ -81,11 +74,6 @@ namespace GridKit using namespace Governor; using namespace Exciter; - // Set system model tolerances - rel_tol_ = 1e-7; - abs_tol_ = 1e-9; - this->max_steps_ = 2000; - owns_components_ = true; // Add electrical buses diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index 9d7f53ee..12424307 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -97,12 +97,6 @@ namespace GridKit return size_opt_; } - 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_; @@ -138,6 +132,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() override + { + return absTol_; + } + + const std::vector& absoluteTolerance() const override + { + return absTol_; + } + std::vector& yB() { return yB_; @@ -260,6 +264,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector absTol_; std::vector f_; std::vector g_; @@ -277,9 +282,6 @@ namespace GridKit RealT time_; RealT alpha_; - RealT rel_tol_; - RealT abs_tol_; - IdxT max_steps_; IdxT idc_; diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 736f9d00..63e43a24 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -69,8 +69,6 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::f_; using CircuitComponent::jac_; - using CircuitComponent::rel_tol_; - using CircuitComponent::abs_tol_; public: /** @@ -80,10 +78,6 @@ 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; } @@ -91,22 +85,12 @@ namespace GridKit /** * @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; } diff --git a/GridKit/Model/PowerFlow/Bus/BaseBus.hpp b/GridKit/Model/PowerFlow/Bus/BaseBus.hpp index 75068f10..e5dadebb 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_; diff --git a/GridKit/Model/PowerFlow/Bus/BusSlack.hpp b/GridKit/Model/PowerFlow/Bus/BusSlack.hpp index d35b7dc3..65edd619 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/MiniGrid/MiniGrid.cpp b/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.cpp index 25ce8c10..8dc9e8e5 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 19d4a1da..18a201c1 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; diff --git a/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp b/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp index 2457b16b..8cbdb35a 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() override + { + return absTol_; + } + + const std::vector& absoluteTolerance() const override + { + return absTol_; + } + std::vector& yB() { return yB_; @@ -233,6 +237,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector absTol_; 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 d6253cd7..10b51c3c 100644 --- a/GridKit/Model/PowerFlow/SystemModel.hpp +++ b/GridKit/Model/PowerFlow/SystemModel.hpp @@ -43,8 +43,6 @@ namespace GridKit 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_; @@ -55,12 +53,7 @@ namespace GridKit */ SystemModel() : ModelEvaluatorImpl(0, 0, 0) - { - // Set system model tolerances - rel_tol_ = 1e-7; - abs_tol_ = 1e-9; - this->max_steps_ = 2000; - } + {} /** * @brief Destructor for the system model diff --git a/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp b/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp index 88ace8c3..d3d2bbf9 100644 --- a/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp +++ b/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp @@ -52,8 +52,6 @@ namespace GridKit // 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_; @@ -64,11 +62,7 @@ namespace GridKit */ SystemSteadyStateModel() : ModelEvaluatorImpl(0, 0, 0) - { - // Set system model tolerances - rel_tol_ = 1e-5; - abs_tol_ = 1e-5; - } + {} /** * @brief Construct a new System Steady State Model object. Allows for simple allocation. @@ -78,9 +72,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) { diff --git a/GridKit/Solver/Dynamic/DynamicSolver.hpp b/GridKit/Solver/Dynamic/DynamicSolver.hpp index 74489c22..a86aed55 100644 --- a/GridKit/Solver/Dynamic/DynamicSolver.hpp +++ b/GridKit/Solver/Dynamic/DynamicSolver.hpp @@ -23,6 +23,10 @@ namespace AnalysisManager return model_; } + virtual void setTimeStep(ScalarT timeStep) = 0; + virtual void setTolerance(ScalarT relTol, ScalarT absTolFac=1) = 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 ea1900aa..436d2257 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -16,8 +16,8 @@ namespace AnalysisManager { template - Ida::Ida(GridKit::Model::Evaluator* model, RealT dt) - : DynamicSolver(model), dt_(dt) + Ida::Ida(GridKit::Model::Evaluator* model) + : DynamicSolver(model) { int retval = 0; @@ -82,10 +82,8 @@ namespace AnalysisManager retval = IDASetUserData(solver_, model_); checkOutput(retval, "IDASetUserData"); - configureIDA(solver_); - // Tag differential variables - std::vector& tag = model_->tag(); + const std::vector& tag = model_->tag(); if (static_cast(tag.size()) == model_->size()) { tag_ = N_VClone(yy_); @@ -99,6 +97,12 @@ namespace AnalysisManager checkOutput(retval, "IDASetSuppressAlg"); } + const std::vector& absTol = model_->absoluteTolerance(); + absTol_ = N_VClone(yy_); + checkAllocation((void*) absTol_, "N_VClone"); + copyVec(absTol, absTol_); + setTolerance(DEFAULT_REL_TOL); + // Set up linear solver return this->configureLinearSolver(); } @@ -142,64 +146,6 @@ namespace AnalysisManager return retval; } - /** - * @brief Configure IDA for forward or backward simulation - * - * @tparam ScalarT - * @tparam IdxT - */ - template - void Ida::configureIDA(void *ida) - { - int retval = 0; - - // Set tolerances - RealT rel_tol; - RealT abs_tol; - - model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"! - - if (dt_ > 0) { - retval = IDASetMinStep(ida, dt_); - checkOutput(retval, "IDASetMinStep"); - retval = IDASetMaxStep(ida, dt_); - checkOutput(retval, "IDASetMaxStep"); - - /* Since the starting procedure is first order, the maximum global order - * of convergence is two */ - IDASetMaxOrd(ida, 2); - checkOutput(retval, "IDASetMaxOrd"); - - - /* Enable more nonlinear iterations because a failed nonlinear solve - * causes a failed integration with fixed steps */ - retval = IDASetMaxNonlinIters(ida, 100); - checkOutput(retval, "IDASetMaxNonlinIters"); - - // Set a large tolerance so the error test will never fail - static constexpr RealT FIXED_STEP_TOL_FAC = 1e10; - retval = IDASStolerances(ida, FIXED_STEP_TOL_FAC, - FIXED_STEP_TOL_FAC * abs_tol / rel_tol); - checkOutput(retval, "IDASStolerances"); - - /* 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(ida, rel_tol / FIXED_STEP_TOL_FAC); - checkOutput(retval, "IDASetNonlinConvCoef"); - } else { - retval = IDASStolerances(ida, 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(ida, static_cast(msa)); - checkOutput(retval, "IDASetMaxNumSteps"); - } - #ifdef GRIDKIT_ENABLE_SUNDIALS_SPARSE /** * @brief Configure a sparse linear solver @@ -418,6 +364,7 @@ namespace AnalysisManager N_VDestroy(yy_); N_VDestroy(yp_); N_VDestroy(tag_); + N_VDestroy(absTol_); N_VDestroy(yy0_); N_VDestroy(yp0_); SUNLinSolFree(linearSolver_); @@ -446,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, 0.1); // Include quadrature in eror checking retval = IDASetQuadErrCon(solver_, SUNTRUE); @@ -588,8 +530,6 @@ namespace AnalysisManager int Ida::initializeBackwardSimulation(RealT tf) { int retval = 0; - RealT rel_tol; - RealT abs_tol; model_->initializeAdjoint(); @@ -604,11 +544,11 @@ namespace AnalysisManager retval = IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_); checkOutput(retval, "IDAInitB"); + setBackwardTolerance(DEFAULT_REL_TOL); + retval = IDASetUserDataB(solver_, backwardID_, model_); checkOutput(retval, "IDASetUserDataB"); - configureIDA(IDAGetAdjIDABmem(solver_, backwardID_)); - // Allocate Jacobian matrix, if not already if (JacobianMatB_ == nullptr) { @@ -633,10 +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); - model_->setTolerances(rel_tol, abs_tol); - retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol * 0.1, abs_tol * 0.1); - checkOutput(retval, "IDAQuadSStolerancesB"); + setBackwardQuadratureTolerance(0.1 * DEFAULT_REL_TOL, 0.1); // Include quadratures in error control retval = IDASetQuadErrConB(solver_, backwardID_, SUNTRUE); @@ -958,7 +895,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_); @@ -982,7 +919,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)); @@ -1002,7 +939,7 @@ namespace AnalysisManager * @tparam IdxT */ template - void Ida::printFinalStats() + void Ida::printFinalStats() const { int retval = IDAPrintAllStats(solver_, stdout, SUN_OUTPUTFORMAT_TABLE); checkOutput(retval, "IDAPrintAllStats"); @@ -1040,6 +977,127 @@ namespace AnalysisManager } } + template + void Ida::setTimeStep(ScalarT timeStep) { + setTimeStep(solver_, timeStep); + } + + template + void Ida::setBackwardTimeStep(ScalarT timeStep) { + setTimeStep(IDAGetAdjIDABmem(solver_, backwardID_), timeStep); + } + + template + void Ida::setTimeStep(void *mem, ScalarT timeStep) { + int retval = 0; + retval = IDASetMinStep(mem, timeStep); + checkOutput(retval, "IDASetMinStep"); + retval = IDASetMaxStep(mem, timeStep); + checkOutput(retval, "IDASetMaxStep"); + } + + template + void Ida::setFixedStep(ScalarT timeStep, ScalarT relTol, ScalarT absTolFac) { + setFixedStep(solver_, timeStep, relTol, absTolFac); + } + + template + void Ida::setBackwardFixedStep(ScalarT timeStep, ScalarT relTol, ScalarT absTolFac) { + setFixedStep(IDAGetAdjIDABmem(solver_, backwardID_), timeStep, relTol, absTolFac); + } + + template + void Ida::setFixedStep(void *mem, ScalarT timeStep, ScalarT relTol, ScalarT absTolFac) { + setTimeStep(mem, timeStep); + + if (timeStep != 0) { + int retval = 0; + /* 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 * relTol, FIXED_STEP_TOL_FAC * absTolFac); + + /* 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"); + } + } + + template + void Ida::setTolerance(ScalarT relTol, ScalarT absTolFac) { + setTolerance(solver_, relTol, absTolFac); + } + + template + void Ida::setBackwardTolerance(ScalarT relTol, ScalarT absTolFac) { + setTolerance(IDAGetAdjIDABmem(solver_, backwardID_), relTol, absTolFac); + } + + template + void Ida::setTolerance(void *mem, ScalarT relTol, ScalarT absTolFac) { + int retval = 0; + if (absTolFac == 1) { + retval = IDASVtolerances(mem, relTol, absTol_); + } else { + N_Vector tmp = N_VClone(absTol_); + N_VScale(absTolFac, absTol_, tmp); + retval = IDASVtolerances(mem, relTol, tmp); + N_VDestroy(tmp); + } + checkOutput(retval, "IDASVtolerances"); + } + + template + void Ida::setMaxSteps(IdxT maxSteps) { + setMaxSteps(solver_, maxSteps); + } + + template + void Ida::setBackwardMaxSteps(IdxT maxSteps) { + setMaxSteps(IDAGetAdjIDABmem(solver_, backwardID_), maxSteps); + } + + template + void Ida::setMaxSteps(void *mem, IdxT maxSteps) { + int retval = IDASetMaxNumSteps(mem, maxSteps); + checkOutput(retval, "IDASetMaxNumSteps"); + } + + template + void Ida::setQuadratureTolerance(ScalarT relTol, ScalarT absTolFac) { + setQuadratureTolerance(solver_, relTol, absTolFac); + } + + template + void Ida::setBackwardQuadratureTolerance(ScalarT relTol, ScalarT absTolFac) { + setQuadratureTolerance(IDAGetAdjIDABmem(solver_, backwardID_), relTol, absTolFac); + } + + template + void Ida::setQuadratureTolerance(void *mem, ScalarT relTol, ScalarT absTolFac) { + int retval = 0; + if (absTolFac == 1) { + retval = IDAQuadSVtolerances(mem, relTol, absTol_); + } else { + N_Vector tmp = N_VClone(absTol_); + N_VScale(absTolFac, absTol_, tmp); + retval = IDAQuadSVtolerances(mem, relTol, tmp); + N_VDestroy(tmp); + } + 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 2b424048..d51c6c16 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -32,15 +32,15 @@ namespace AnalysisManager using RealT = typename GridKit::ScalarTraits::RealT; public: - Ida(GridKit::Model::Evaluator* model, RealT dt=0); + Ida(GridKit::Model::Evaluator* model); ~Ida(); int configureSimulation(); - int configureLinearSolver(); + int configureLinearSolver(); // TODO: make private #ifdef GRIDKIT_ENABLE_SUNDIALS_SPARSE - int configureLinearSolverSparse(); + int configureLinearSolverSparse(); // TODO: make private #endif - int configureLinearSolverDense(); + int configureLinearSolverDense(); // TODO: make private int getDefaultInitialCondition(); int setIntegrationTime(RealT t_init, RealT t_final, int nout); int initializeSimulation(RealT t0, bool findConsistent = false); @@ -111,9 +111,20 @@ 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 setTimeStep(ScalarT timeStep) override; + void setBackwardTimeStep(ScalarT timeStep); + void setFixedStep(ScalarT timeStep, ScalarT relTol, ScalarT absTolFac=1); + void setBackwardFixedStep(ScalarT timeStep, ScalarT relTol, ScalarT absTolFac=1); + void setTolerance(ScalarT relTol, ScalarT absTolFac=1) override; + void setBackwardTolerance(ScalarT relTol, ScalarT absTolFac=1); + void setQuadratureTolerance(ScalarT relTol, ScalarT); + void setBackwardQuadratureTolerance(ScalarT relTol, ScalarT); + void setMaxSteps(IdxT maxSteps) override; + void setBackwardMaxSteps(IdxT maxSteps); private: static int Residual(RealT t, @@ -156,7 +167,7 @@ namespace AnalysisManager void* user_data); private: - RealT dt_; + static constexpr ScalarT DEFAULT_REL_TOL = 1e-5; void* solver_{}; SUNContext context_{}; @@ -172,6 +183,7 @@ namespace AnalysisManager N_Vector yy_{}; ///< Solution vector N_Vector yp_{}; ///< Solution derivatives vector N_Vector tag_{}; ///< Tags differential variables + N_Vector absTol_{}; ///< Tags differential variables N_Vector q_{}; ///< Integrand vector N_Vector yy0_{}; ///< Storage for initial values @@ -184,16 +196,20 @@ namespace AnalysisManager int backwardID_{}; private: - void configureIDA(void *ida); - // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); static void copyVec(const N_Vector x, std::vector& y); static void copyVec(const std::vector& x, N_Vector y); 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 setTimeStep(void *mem, ScalarT timeStep); + void setFixedStep(void *mem, ScalarT timeStep, ScalarT relTol, ScalarT absTolFac); + void setTolerance(void *mem, ScalarT relTol, ScalarT absTolFac); + void setMaxSteps(void *mem, IdxT maxSteps); + void setQuadratureTolerance(void *mem, ScalarT relTol, ScalarT absTolFac); }; /// Simple exception to use within Ida class. diff --git a/GridKit/Solver/SteadyState/Kinsol.cpp b/GridKit/Solver/SteadyState/Kinsol.cpp index 7e305318..dc08485f 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,15 @@ 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 3b4b5c40..97ce474d 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,11 @@ 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 05fa9da5..19643e33 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/examples/PowerElectronics/Microgrid/Microgrid.cpp b/examples/PowerElectronics/Microgrid/Microgrid.cpp index fc2f379e..b1eee826 100644 --- a/examples/PowerElectronics/Microgrid/Microgrid.cpp +++ b/examples/PowerElectronics/Microgrid/Microgrid.cpp @@ -19,13 +19,12 @@ 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; // 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; @@ -337,6 +336,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 dce13668..b3f6326d 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -16,13 +16,12 @@ int main(int /* argc */, char const** /* argv */) { - double abs_tol = 1.0e-8; double rel_tol = 1.0e-8; bool use_jac = true; // 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; @@ -107,6 +106,7 @@ int main(int /* argc */, char const** /* argv */) // setup simulation idas.configureSimulation(); + idas.setTolerance(rel_tol); idas.getDefaultInitialCondition(); idas.initializeSimulation(t_init); diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp index b9e4ca2b..d76f272a 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -70,13 +70,9 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) 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 - 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; @@ -333,6 +329,8 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) // setup simulation idas->configureSimulation(); + idas->setTolerance(rel_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 ffde3e7c..fefb5ce6 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp @@ -66,13 +66,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) @@ -304,6 +300,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 da4d49ff..2695b487 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -61,15 +61,6 @@ namespace GridKit return 0; } - void setTolerances([[maybe_unused]] RealT& rel_tol, [[maybe_unused]] RealT& abs_tol) const override - { - } - - void setMaxSteps(IdxT& msa) const override - { - msa = 2000; - } - int tagDifferentiable() override { return 0; @@ -140,6 +131,16 @@ namespace GridKit return tag_; } + std::vector& absoluteTolerance() override + { + return absTol_; + } + + const std::vector& absoluteTolerance() const override + { + return absTol_; + } + std::vector& yB() override { return yB_; @@ -249,6 +250,7 @@ namespace GridKit std::vector y_; std::vector yp_; std::vector tag_; + std::vector absTol_; std::vector f_; std::vector g_; From d92935f60b41c044f3aa6f5c9a64576fce350899 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 30 Jan 2026 18:41:18 -0800 Subject: [PATCH 5/6] Add absTol resizing on allocate --- GridKit/Model/PhasorDynamics/Bus/Bus.hpp | 1 + GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp | 1 + GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 1 + GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp | 1 + GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp | 1 + GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp | 1 + .../PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp | 1 + .../PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp | 1 + .../SynchronousMachine/GenClassical/GenClassical.hpp | 1 + .../SynchronousMachine/GenClassical/GenClassicalImpl.hpp | 1 + .../Model/PowerElectronics/SystemModelPowerElectronics.hpp | 4 ++++ GridKit/Model/PowerFlow/Bus/BusPQ.cpp | 1 + GridKit/Model/PowerFlow/Bus/BusPQ.hpp | 1 + GridKit/Model/PowerFlow/Bus/BusPV.cpp | 1 + GridKit/Model/PowerFlow/Bus/BusPV.hpp | 1 + 15 files changed, 18 insertions(+) diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp index 59493689..22e007d4 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp @@ -25,6 +25,7 @@ namespace GridKit using BusBase::f_; using BusBase::J_; using BusBase::tag_; + using BusBase::absTol_; using BusBase::variable_indices_; using BusBase::residual_indices_; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 37b70c01..4b00e39a 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); + absTol_.resize(size); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index fffb5674..0b380aa2 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -72,6 +72,7 @@ namespace GridKit using Component::nnz_; using Component::size_; using Component::tag_; + using Component::absTol_; using Component::time_; using Component::y_; using Component::yp_; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 6e935124..f316b1a0 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); + absTol_.resize(size); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index 95bb6fda..bc59cd95 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::absTol_; using Component::time_; using Component::y_; using Component::yp_; diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index f6add0df..864b5a08 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); + absTol_.resize(size_); // Resize external variable data ws_.resize(1); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp index 8c9aaa61..fc680283 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::absTol_; using Component::time_; using Component::y_; using Component::yp_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index 7f3ae05f..3dc6afb2 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -314,6 +314,7 @@ namespace GridKit y_.resize(static_cast(size_)); yp_.resize(static_cast(size_)); tag_.resize(static_cast(size_)); + absTol_.resize(static_cast(size_)); // Resize coupling data wb_.resize(2); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index 7a6d69df..d47479aa 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::absTol_; using Component::time_; using Component::y_; using Component::yp_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index 58ba7039..4ab2db3d 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -170,6 +170,7 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + absTol_.resize(size); // Resize coupling data wb_.resize(2); diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 63e43a24..06e9267c 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -69,6 +69,8 @@ namespace GridKit using CircuitComponent::yp_; using CircuitComponent::f_; using CircuitComponent::jac_; + using CircuitComponent::tag_; + using CircuitComponent::absTol_; public: /** @@ -167,6 +169,8 @@ namespace GridKit y_.resize(size_); yp_.resize(size_); f_.resize(size_); + tag_.resize(size_); + absTol_.resize(size_); return 0; } diff --git a/GridKit/Model/PowerFlow/Bus/BusPQ.cpp b/GridKit/Model/PowerFlow/Bus/BusPQ.cpp index 7b53a9c8..0666eda6 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_)); + absTol_.resize(static_cast(size_)); fB_.resize(static_cast(size_)); yB_.resize(static_cast(size_)); diff --git a/GridKit/Model/PowerFlow/Bus/BusPQ.hpp b/GridKit/Model/PowerFlow/Bus/BusPQ.hpp index e247f593..780499c9 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::absTol_; public: using RealT = typename ModelEvaluatorImpl::RealT; diff --git a/GridKit/Model/PowerFlow/Bus/BusPV.cpp b/GridKit/Model/PowerFlow/Bus/BusPV.cpp index c3312c70..ea2e9f12 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_)); + absTol_.resize(static_cast(size_)); fB_.resize(static_cast(size_)); yB_.resize(static_cast(size_)); diff --git a/GridKit/Model/PowerFlow/Bus/BusPV.hpp b/GridKit/Model/PowerFlow/Bus/BusPV.hpp index e56ad69e..62e45294 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::absTol_; public: using RealT = typename ModelEvaluatorImpl::RealT; From 78ff43da971bdf54088226e5644e2a441eb175d1 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 30 Jan 2026 19:46:47 -0800 Subject: [PATCH 6/6] Add function to set abs tol --- GridKit/Model/Evaluator.hpp | 13 +++---- .../Model/PhasorDynamics/Branch/Branch.hpp | 1 + .../PhasorDynamics/Branch/BranchImpl.hpp | 9 +++++ GridKit/Model/PhasorDynamics/Bus/Bus.hpp | 1 + GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp | 9 +++++ .../Model/PhasorDynamics/Bus/BusInfinite.hpp | 1 + .../PhasorDynamics/Bus/BusInfiniteImpl.hpp | 9 +++++ .../PhasorDynamics/BusFault/BusFault.hpp | 1 + .../PhasorDynamics/BusFault/BusFaultImpl.hpp | 9 +++++ .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 1 + .../Exciter/IEEET1/Ieeet1Impl.hpp | 9 +++++ .../PhasorDynamics/Governor/Tgov1/Tgov1.hpp | 1 + .../Governor/Tgov1/Tgov1Impl.hpp | 9 +++++ GridKit/Model/PhasorDynamics/Load/Load.hpp | 1 + .../Model/PhasorDynamics/Load/LoadImpl.hpp | 9 +++++ .../PhasorDynamics/SignalNode/SignalNode.hpp | 1 + .../SignalNode/SignalNodeImpl.hpp | 6 ++++ .../SynchronousMachine/GENROUwS/Genrou.hpp | 1 + .../GENROUwS/GenrouImpl.hpp | 9 +++++ .../GenClassical/GenClassical.hpp | 1 + .../GenClassical/GenClassicalImpl.hpp | 9 +++++ GridKit/Model/PhasorDynamics/SystemModel.hpp | 35 +++++++++++++++++++ .../DistributedGenerator.cpp | 9 +++++ .../DistributedGenerator.hpp | 1 + .../PowerElectronics/Inductor/Inductor.cpp | 9 +++++ .../PowerElectronics/Inductor/Inductor.hpp | 1 + .../MicrogridBusDQ/MicrogridBusDQ.cpp | 9 +++++ .../MicrogridBusDQ/MicrogridBusDQ.hpp | 1 + .../MicrogridLine/MicrogridLine.cpp | 9 +++++ .../MicrogridLine/MicrogridLine.hpp | 1 + .../MicrogridLoad/MicrogridLoad.cpp | 9 +++++ .../MicrogridLoad/MicrogridLoad.hpp | 1 + .../PowerElectronics/Resistor/Resistor.cpp | 9 +++++ .../PowerElectronics/Resistor/Resistor.hpp | 1 + .../SystemModelPowerElectronics.hpp | 5 +++ .../VoltageSource/VoltageSource.cpp | 9 +++++ .../VoltageSource/VoltageSource.hpp | 1 + GridKit/Model/PowerFlow/Branch/Branch.cpp | 9 +++++ GridKit/Model/PowerFlow/Branch/Branch.hpp | 1 + GridKit/Model/PowerFlow/Bus/BaseBus.hpp | 5 +++ GridKit/Model/PowerFlow/Bus/BusPQ.cpp | 6 ++++ GridKit/Model/PowerFlow/Bus/BusPQ.hpp | 1 + GridKit/Model/PowerFlow/Bus/BusPV.cpp | 6 ++++ GridKit/Model/PowerFlow/Bus/BusPV.hpp | 1 + .../PowerFlow/Generator/GeneratorBase.hpp | 5 +++ .../Model/PowerFlow/Generator2/Generator2.cpp | 1 + .../Model/PowerFlow/Generator2/Generator2.hpp | 1 + .../Model/PowerFlow/Generator4/Generator4.cpp | 10 ++++++ .../Model/PowerFlow/Generator4/Generator4.hpp | 2 ++ .../Generator4Governor/Generator4Governor.cpp | 1 + .../Generator4Governor/Generator4Governor.hpp | 1 + .../Generator4Param/Generator4Param.cpp | 1 + .../Generator4Param/Generator4Param.hpp | 1 + GridKit/Model/PowerFlow/Load/Load.cpp | 9 +++++ GridKit/Model/PowerFlow/Load/Load.hpp | 1 + GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp | 5 +++ GridKit/Model/PowerFlow/SystemModel.hpp | 35 +++++++++++++++++++ .../Model/PowerFlow/SystemModelPowerFlow.hpp | 11 ++++++ tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 5 +++ 59 files changed, 332 insertions(+), 6 deletions(-) diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index e21c838f..88ec2ac9 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -31,12 +31,13 @@ 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; + virtual int setAbsoluteTolerance() = 0; + virtual int evaluateResidual() = 0; + virtual int evaluateJacobian() = 0; + virtual int evaluateIntegrand() = 0; virtual int initializeAdjoint() = 0; virtual int evaluateAdjointResidual() = 0; diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index c917e431..62029552 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -68,6 +68,7 @@ namespace GridKit virtual int allocate() override; virtual int initialize() override; virtual int tagDifferentiable() override; + virtual int setAbsoluteTolerance() override; virtual int evaluateResidual() override; virtual int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp index 4b6016e3..0660a751 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp @@ -140,6 +140,15 @@ namespace GridKit return 0; } + /** + * @brief Set absolute tolerance + */ + template + int Branch::setAbsoluteTolerance() + { + 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 22e007d4..79c860c7 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp @@ -42,6 +42,7 @@ namespace GridKit virtual int setBusID(IdxT) override; virtual int allocate() override; virtual int tagDifferentiable() override; + virtual int setAbsoluteTolerance() override; virtual int initialize() override; virtual int evaluateResidual() override; virtual int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 4b00e39a..0e42ce78 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -119,6 +119,15 @@ namespace GridKit return 0; } + /** + * @brief Set absolute tolerance + */ + template + int Bus::setAbsoluteTolerance() + { + 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 ad97593f..1d5b094f 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp @@ -36,6 +36,7 @@ namespace GridKit virtual int setBusID(IdxT) override; virtual int allocate() override; virtual int tagDifferentiable() override; + virtual int setAbsoluteTolerance() override; virtual int initialize() override; virtual int evaluateResidual() override; virtual int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp index 56d05d82..7b7ae784 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp @@ -95,6 +95,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int BusInfinite::setAbsoluteTolerance() + { + 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 77c85e42..262ac85d 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp @@ -49,6 +49,7 @@ namespace GridKit int allocate() override; int initialize() override; int tagDifferentiable() override; + int setAbsoluteTolerance() override; int evaluateResidual() override; int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 97761a53..b1013c97 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -141,6 +141,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int BusFault::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Bus residual * diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index 0b380aa2..42b26d1c 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -98,6 +98,7 @@ namespace GridKit int verify() const override; int initialize() override; int tagDifferentiable() override; + int setAbsoluteTolerance() override; int evaluateResidual() override; int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index f316b1a0..231b7efd 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -249,6 +249,15 @@ namespace GridKit return 0; } + /** + * @brief Set absolute tolerance + */ + template + int Ieeet1::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Residuals of system equations * diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index bc59cd95..7ef24700 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -87,6 +87,7 @@ namespace GridKit int verify() const override; int initialize() override; int tagDifferentiable() override; + int setAbsoluteTolerance() override; int evaluateResidual() override; // Still to be implemented diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 864b5a08..2b4b54e7 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -230,6 +230,15 @@ namespace GridKit return 0; } + /** + * @brief Set absolute tolerance + */ + template + int Tgov1::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Internal residuals * diff --git a/GridKit/Model/PhasorDynamics/Load/Load.hpp b/GridKit/Model/PhasorDynamics/Load/Load.hpp index 091853b5..23990d54 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.hpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.hpp @@ -55,6 +55,7 @@ namespace GridKit virtual int allocate() override; virtual int initialize() override; virtual int tagDifferentiable() override; + virtual int setAbsoluteTolerance() override; virtual int evaluateResidual() override; virtual int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp index 144d49eb..9447e749 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp @@ -110,6 +110,15 @@ namespace GridKit return 0; } + /** + * @brief Set absolute tolerance + */ + template + int Load::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Bus residual * diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp index f2226df3..739f7023 100644 --- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp +++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp @@ -37,6 +37,7 @@ namespace GridKit virtual int allocate() override; virtual int initialize() override; virtual int tagDifferentiable() override; + virtual int setAbsoluteTolerance() override; virtual int evaluateResidual() override; virtual int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeImpl.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeImpl.hpp index c4556d20..aa0b70e4 100644 --- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeImpl.hpp @@ -39,6 +39,12 @@ namespace GridKit return 0; } + template + int SignalNode::setAbsoluteTolerance() + { + return 0; + } + template int SignalNode::evaluateResidual() { diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp index fc680283..2964fcfd 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp @@ -129,6 +129,7 @@ namespace GridKit int verify() const override; int initialize() override; int tagDifferentiable() override; + int setAbsoluteTolerance() override; int evaluateResidual() override; // Still to be implemented diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index 3dc6afb2..f4006388 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -469,6 +469,15 @@ namespace GridKit return 0; } + /** + * @brief Set absolute tolerance + */ + template + int Genrou::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Internal residual * diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index d47479aa..321456cf 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -70,6 +70,7 @@ namespace GridKit int allocate() override; int initialize() override; int tagDifferentiable() override; + int setAbsoluteTolerance() override; int evaluateResidual() override; int verify() const override diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index 4ab2db3d..7add51f3 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -233,6 +233,15 @@ namespace GridKit return 0; } + /** + * @brief Set absolute tolerance + */ + template + int GenClassical::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Internal residual * diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 223cd837..733ea39c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -47,6 +47,7 @@ namespace GridKit using PhasorDynamics::Component::y_; using PhasorDynamics::Component::yp_; using PhasorDynamics::Component::tag_; + using PhasorDynamics::Component::absTol_; using PhasorDynamics::Component::f_; using PhasorDynamics::Component::J_; @@ -315,6 +316,7 @@ namespace GridKit yp_.resize(size_); f_.resize(size_); tag_.resize(size_); + absTol_.resize(size_); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) @@ -483,6 +485,39 @@ namespace GridKit return 0; } + /** + * @todo Specify absolute tolerance + * + * Specify a "noise" level close to zero for which pure relative error + * cannot be used. + */ + int setAbsoluteTolerance() + { + // Set initial values for global solution vectors + IdxT offset = 0; + for (const auto& bus : buses_) + { + bus->setAbsoluteTolerance(); + for (IdxT j = 0; j < bus->size(); ++j) + { + absTol_[offset + j] = bus->absoluteTolerance()[j]; + } + offset += bus->size(); + } + + for (const auto& component : components_) + { + component->setAbsoluteTolerance(); + for (IdxT j = 0; j < component->size(); ++j) + { + absTol_[offset + j] = component->absoluteTolerance()[j]; + } + offset += component->size(); + } + + return 0; + } + /** * @brief Compute system residual vector * diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index 5e1eee4d..dc2f3db4 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -82,6 +82,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance. + */ + template + int DistributedGenerator::setAbsoluteTolerance() + { + 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 d2c930ac..a952d7f1 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -72,6 +72,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp index ebf70788..c7c7acd3 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp @@ -63,6 +63,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int Inductor::setAbsoluteTolerance() + { + 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 a45042ec..354fe20b 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp @@ -50,6 +50,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp index fd5dbe8f..8796104a 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -67,6 +67,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int MicrogridBusDQ::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Evaluate residual of microgrid line * This model has "Virtual resistors". The voltage of the bus divided by its virtual resistance. diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index 8c0d8dad..49466f31 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -49,6 +49,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index 44f9073d..a167780e 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -71,6 +71,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int MicrogridLine::setAbsoluteTolerance() + { + 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 64992c5c..ca9cd5ca 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -50,6 +50,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index e16b70b6..6b20ddca 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -69,6 +69,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int MicrogridLoad::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Eval Micro Load */ diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 0a8f34bb..34b6a7fc 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -50,6 +50,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp index e42eefff..c8505432 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp @@ -62,6 +62,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int Resistor::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Computes the resistors resisdual * diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp index 1f2dac7f..1680f418 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp @@ -49,6 +49,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 06e9267c..09647811 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -225,6 +225,11 @@ namespace GridKit return 0; } + int setAbsoluteTolerance() + { + 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 1625d378..9c76419c 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp @@ -62,6 +62,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int VoltageSource::setAbsoluteTolerance() + { + return 0; + } + /** * @brief Evaluate resisdual of component */ diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp index ef4a9474..2270299f 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp @@ -49,6 +49,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Branch/Branch.cpp b/GridKit/Model/PowerFlow/Branch/Branch.cpp index 95048043..20a7b10d 100644 --- a/GridKit/Model/PowerFlow/Branch/Branch.cpp +++ b/GridKit/Model/PowerFlow/Branch/Branch.cpp @@ -96,6 +96,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance. + */ + template + int Branch::setAbsoluteTolerance() + { + 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 57c49813..fb4e324d 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(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Bus/BaseBus.hpp b/GridKit/Model/PowerFlow/Bus/BaseBus.hpp index e5dadebb..1fe69006 100644 --- a/GridKit/Model/PowerFlow/Bus/BaseBus.hpp +++ b/GridKit/Model/PowerFlow/Bus/BaseBus.hpp @@ -74,6 +74,11 @@ namespace GridKit return 0; } + virtual int setAbsoluteTolerance() + { + return 0; + } + virtual int evaluateResidual() { return 0; diff --git a/GridKit/Model/PowerFlow/Bus/BusPQ.cpp b/GridKit/Model/PowerFlow/Bus/BusPQ.cpp index 0666eda6..f0c30bcf 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPQ.cpp +++ b/GridKit/Model/PowerFlow/Bus/BusPQ.cpp @@ -91,6 +91,12 @@ namespace GridKit return 0; } + template + int BusPQ::setAbsoluteTolerance() + { + 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 780499c9..92002755 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPQ.hpp +++ b/GridKit/Model/PowerFlow/Bus/BusPQ.hpp @@ -38,6 +38,7 @@ namespace GridKit virtual int allocate(); virtual int tagDifferentiable(); + virtual int setAbsoluteTolerance(); 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 ea2e9f12..8db383f3 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPV.cpp +++ b/GridKit/Model/PowerFlow/Bus/BusPV.cpp @@ -88,6 +88,12 @@ namespace GridKit return 0; } + template + int BusPV::setAbsoluteTolerance() + { + 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 62e45294..80e7f45e 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPV.hpp +++ b/GridKit/Model/PowerFlow/Bus/BusPV.hpp @@ -40,6 +40,7 @@ namespace GridKit virtual int allocate(); virtual int tagDifferentiable(); + virtual int setAbsoluteTolerance(); virtual int initialize(); virtual int evaluateResidual(); virtual int initializeAdjoint(); diff --git a/GridKit/Model/PowerFlow/Generator/GeneratorBase.hpp b/GridKit/Model/PowerFlow/Generator/GeneratorBase.hpp index 5750c686..a4e796a8 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() + { + return 0; + } + virtual int evaluateResidual() { return 0; diff --git a/GridKit/Model/PowerFlow/Generator2/Generator2.cpp b/GridKit/Model/PowerFlow/Generator2/Generator2.cpp index 0e61a44d..85f8cd3a 100644 --- a/GridKit/Model/PowerFlow/Generator2/Generator2.cpp +++ b/GridKit/Model/PowerFlow/Generator2/Generator2.cpp @@ -49,6 +49,7 @@ namespace GridKit int Generator2::allocate() { tag_.resize(static_cast(size_)); + absTol_.resize(static_cast(size_)); return 0; } diff --git a/GridKit/Model/PowerFlow/Generator2/Generator2.hpp b/GridKit/Model/PowerFlow/Generator2/Generator2.hpp index c022279c..fd0832fe 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::absTol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; diff --git a/GridKit/Model/PowerFlow/Generator4/Generator4.cpp b/GridKit/Model/PowerFlow/Generator4/Generator4.cpp index edf91b3d..7a2dff7d 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_)); + absTol_.resize(static_cast(size_)); return 0; } @@ -149,6 +150,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance + */ + template + int Generator4::setAbsoluteTolerance() + { + 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 de17996a..8328c599 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::absTol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; @@ -45,6 +46,7 @@ namespace GridKit int allocate(); int initialize(); int tagDifferentiable(); + int setAbsoluteTolerance(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp index 284346f6..9670d716 100644 --- a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp +++ b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp @@ -67,6 +67,7 @@ namespace GridKit { // std::cout << "Allocate Gen2..." << std::endl; tag_.resize(static_cast(size_)); + absTol_.resize(static_cast(size_)); return 0; } diff --git a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.hpp b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.hpp index 5a34067a..fb157aeb 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::absTol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; diff --git a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp index 18b23551..23ae38f7 100644 --- a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp +++ b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp @@ -54,6 +54,7 @@ namespace GridKit { // std::cout << "Allocate Generator4Param..." << std::endl; tag_.resize(static_cast(size_)); + absTol_.resize(static_cast(size_)); return 0; } diff --git a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.hpp b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.hpp index 3ca69b50..b7a26f68 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::absTol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::yB_; diff --git a/GridKit/Model/PowerFlow/Load/Load.cpp b/GridKit/Model/PowerFlow/Load/Load.cpp index 366f5624..61afe9d0 100644 --- a/GridKit/Model/PowerFlow/Load/Load.cpp +++ b/GridKit/Model/PowerFlow/Load/Load.cpp @@ -70,6 +70,15 @@ namespace GridKit return 0; } + /** + * \brief Specify absolute tolerance. + */ + template + int Load::setAbsoluteTolerance() + { + 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 2326ba5b..8a21bc34 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(); int evaluateResidual(); int evaluateJacobian(); int evaluateIntegrand(); diff --git a/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp b/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp index 18a201c1..dda903d5 100644 --- a/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp +++ b/GridKit/Model/PowerFlow/MiniGrid/MiniGrid.hpp @@ -33,6 +33,11 @@ namespace GridKit return -1; } + int setAbsoluteTolerance() + { + return -1; + } + int evaluateResidual(); int evaluateJacobian(); diff --git a/GridKit/Model/PowerFlow/SystemModel.hpp b/GridKit/Model/PowerFlow/SystemModel.hpp index 10b51c3c..e48c772f 100644 --- a/GridKit/Model/PowerFlow/SystemModel.hpp +++ b/GridKit/Model/PowerFlow/SystemModel.hpp @@ -39,6 +39,7 @@ namespace GridKit using ModelEvaluatorImpl::yB_; using ModelEvaluatorImpl::ypB_; using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::absTol_; using ModelEvaluatorImpl::f_; using ModelEvaluatorImpl::fB_; using ModelEvaluatorImpl::g_; @@ -106,6 +107,7 @@ namespace GridKit f_.resize(size_); fB_.resize(size_); tag_.resize(size_); + absTol_.resize(size_); g_.resize(size_quad_); gB_.resize(size_quad_ * size_opt_); @@ -235,6 +237,39 @@ namespace GridKit return 0; } + /** + * @todo Specify absolute tolerance + * + * Specify a "noise" level close to zero for which pure relative error + * cannot be used. + */ + int setAbsoluteTolerance() + { + // Set initial values for global solution vectors + IdxT offset = 0; + for (const auto& bus : buses_) + { + bus->setAbsoluteTolerance(); + for (IdxT j = 0; j < bus->size(); ++j) + { + absTol_[offset + j] = bus->absoluteTolerance()[j]; + } + offset += bus->size(); + } + + for (const auto& component : components_) + { + component->setAbsoluteTolerance(); + for (IdxT j = 0; j < component->size(); ++j) + { + absTol_[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 d3d2bbf9..ffc04051 100644 --- a/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp +++ b/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp @@ -215,6 +215,17 @@ namespace GridKit return 0; } + /** + * @todo Specify absolute tolerance + * + * Specify a "noise" level close to zero for which pure relative error + * cannot be used. + */ + int setAbsoluteTolerance() + { + return 0; + } + /** * @brief Compute system residual vector * diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index 2695b487..bf92044f 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -66,6 +66,11 @@ namespace GridKit return 0; } + int setAbsoluteTolerance() override + { + return 0; + } + int evaluateResidual() override { f_ = y_;