From 2c53069f09c61e27500721cc69c4a0a47e8cf167 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 18 Aug 2026 11:32:33 -0700 Subject: [PATCH 1/7] Add option to pick IDA initialization type --- GridKit/Solver/Dynamic/Ida.cpp | 39 +++++++- GridKit/Solver/Dynamic/Ida.hpp | 10 ++ .../PhasorDynamics/AnalysisUtilities.hpp | 21 ++++ .../PhasorDynamics/ContingencyAnalysis.cpp | 1 + .../PhasorDynamics/DynamicSimulation.cpp | 1 + application/PhasorDynamics/README.md | 1 + tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 98 +++++++++++++++++++ .../UnitTests/Solver/Dynamic/runIdaTests.cpp | 1 + 8 files changed, 167 insertions(+), 5 deletions(-) diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index f2f2d5446..2a71b61fe 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -238,11 +238,7 @@ namespace AnalysisManager // Find a consistent set of initial conditions for DAE if (findConsistent) { - int initType = IDA_Y_INIT; - - if (tag_) - initType = IDA_YA_YDP_INIT; - + const int initType = getIDAInitType(); retval = IDACalcIC(solver_, initType, t0 + 0.1); checkOutput(retval, "IDACalcIC"); @@ -307,6 +303,26 @@ namespace AnalysisManager model_->updateTime(t, 0.0); } + template + int Ida::getIDAInitType() const + { + switch (init_type_) + { + case IdaInitType::Y: + return IDA_Y_INIT; + case IdaInitType::YA_YDP: + return IDA_YA_YDP_INIT; + case IdaInitType::AUTO: + return tag_ ? IDA_YA_YDP_INIT : IDA_Y_INIT; + default: + GridKit::Utilities::Logger::error() + << "Invalid IDA init type " + << static_cast(init_type_) + << ".\n"; + throw SundialsException(); + } + } + /** * @brief Run the IDA solver and optionally produce monitor output every `dt_monitor`. * @@ -1197,6 +1213,19 @@ namespace AnalysisManager backward_suppress_alg_ = suppress; } + /** + * @brief Set the IDA consistent-initial-condition calculation type + * + * @param init_type IDA init type. AUTO preserves the existing behavior. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setInitType(IdaInitType init_type) + { + init_type_ = init_type; + } + /** * @brief Set the maximum number of steps * diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index bbe87f8df..c214f8951 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -40,6 +40,13 @@ namespace AnalysisManager std::string report() const; }; + enum class IdaInitType + { + AUTO, + Y, + YA_YDP + }; + template class Ida : public DynamicSolver { @@ -133,6 +140,7 @@ namespace AnalysisManager ScalarT abs_tol_override = 0); void setSuppressAlgebraicErrors(bool suppress); void setBackwardSuppressAlgebraicErrors(bool suppress); + void setInitType(IdaInitType init_type); void setMaxSteps(IdxT maxSteps) override; void setBackwardMaxSteps(IdxT maxSteps); @@ -180,6 +188,7 @@ namespace AnalysisManager int getMonitorStepCount(RealT tf, RealT dt_monitor) const; RealT getMonitorTime(RealT tf, RealT dt_monitor, int step, int nsteps) const; + int getIDAInitType() const; void updateModelState(RealT t); private: @@ -213,6 +222,7 @@ namespace AnalysisManager RealT abs_tol_override_{}; IdxT max_steps_{}; bool suppress_alg_{false}; + IdaInitType init_type_{IdaInitType::AUTO}; RealT backward_time_step_{}; RealT backward_rel_tol_{DEFAULT_REL_TOL}; diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index bb30d2861..e4643cf65 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -12,6 +12,7 @@ #include #include +#include #include #include @@ -63,6 +64,8 @@ namespace GridKit double dt_fixed; /// maximum number of solver time steps, or 0 for the IDA default std::size_t max_steps; + /// IDA consistent initial condition calculation type + AnalysisManager::Sundials::IdaInitType init_type; /// set of system events std::vector events; /// path to output file @@ -100,6 +103,24 @@ namespace GridKit c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); c.dt_fixed = j.value("dt_fixed", 0.0); c.max_steps = j.value("max_steps", std::size_t{0}); + c.init_type = AnalysisManager::Sundials::IdaInitType::AUTO; + if (j.contains("init_type")) + { + const auto init_type_str = j.at("init_type").get(); + if (init_type_str == "y") + { + c.init_type = AnalysisManager::Sundials::IdaInitType::Y; + } + else if (init_type_str == "ya_ydp") + { + c.init_type = AnalysisManager::Sundials::IdaInitType::YA_YDP; + } + else + { + Log::error() << "Invalid IDA init type \"" << init_type_str << "\"; " + << "must be either \"y\" or \"ya_ydp\""; + } + } for (auto& raw_event : j.at("events")) { diff --git a/application/PhasorDynamics/ContingencyAnalysis.cpp b/application/PhasorDynamics/ContingencyAnalysis.cpp index 31415c6e2..18978d091 100644 --- a/application/PhasorDynamics/ContingencyAnalysis.cpp +++ b/application/PhasorDynamics/ContingencyAnalysis.cpp @@ -38,6 +38,7 @@ TestStatus runStudy(StudyData study_data) ida.setTolerance(study_data.rel_tol, study_data.abs_tol); ida.setFixedStep(study_data.dt_fixed); ida.setMaxSteps(study_data.max_steps); + ida.setInitType(study_data.init_type); ida.configureSimulation(); using EventType = SystemEvent::Type; diff --git a/application/PhasorDynamics/DynamicSimulation.cpp b/application/PhasorDynamics/DynamicSimulation.cpp index c9f8f32ea..a80613498 100644 --- a/application/PhasorDynamics/DynamicSimulation.cpp +++ b/application/PhasorDynamics/DynamicSimulation.cpp @@ -32,6 +32,7 @@ int main(int argc, const char* argv[]) ida.setTolerance(study.rel_tol, study.abs_tol); ida.setFixedStep(study.dt_fixed); ida.setMaxSteps(study.max_steps); + ida.setInitType(study.init_type); ida.configureSimulation(); // Start timer diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index 35ab99a5d..b5b24f8d6 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -11,6 +11,7 @@ `abs_tol` | Absolute solver tolerance override (default: 1.0e-9) `dt_fixed` | Fixed solver time step size, or 0 for adaptive stepping (default: 0) `max_steps` | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0) + `init_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: automatic) `events` | An array of event groups (see [Events](#events) below) `output_file` | Path to output (CSV) file (optional) `reference_file` | A string containing the name of the case (optional) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index 655f83afa..f998842e7 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -369,6 +369,74 @@ namespace GridKit private: RealT t_{}; }; + + template + class ConsistentInitTypeEvaluator : public NullEvaluator + { + protected: + using NullEvaluator::allocated_; + using NullEvaluator::y_; + using NullEvaluator::yp_; + using NullEvaluator::abs_tol_; + using NullEvaluator::tag_; + using NullEvaluator::f_; + + public: + ConsistentInitTypeEvaluator() = default; + + explicit ConsistentInitTypeEvaluator(bool consistent_derivative) + : consistent_derivative_(consistent_derivative) + { + } + + int initialize() override + { + if (!allocated_) + { + this->allocate(); + } + + auto* y = y_.getData(); + auto* yp = yp_.getData(); + auto* abs_tol = abs_tol_.getData(); + auto* f = f_.getData(); + + y[0] = 0.0; + y[1] = 10.0; + yp[0] = consistent_derivative_ ? 1.0 : 0.0; + yp[1] = 0.0; + tag_ = {true, false}; + abs_tol[0] = 0.0; + abs_tol[1] = 0.0; + f[0] = 0.0; + f[1] = 0.0; + y_.setDataUpdated(); + yp_.setDataUpdated(); + abs_tol_.setDataUpdated(); + f_.setDataUpdated(); + return 0; + } + + IdxT size() override + { + return 2; + } + + int evaluateResidual() override + { + auto* f = f_.getData(); + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + + f[0] = yp[0] - 1.0; + f[1] = y[1] - y[0]; + f_.setDataUpdated(); + return 0; + } + + private: + bool consistent_derivative_{false}; + }; } // namespace Model namespace Testing @@ -503,6 +571,36 @@ namespace GridKit return success.report(__func__); } + + TestOutcome initType() + { + TestStatus success = true; + + { + Model::ConsistentInitTypeEvaluator model; + + Ida ida(&model); + ida.setInitType(AnalysisManager::Sundials::IdaInitType::YA_YDP); + ida.configureSimulation(); + ida.initializeSimulation(0.0); + + success *= (std::abs(model.yp().getData()[0] - 1.0) < 1.0e-10); + success *= (std::abs(model.y().getData()[1]) < 1.0e-10); + } + + { + Model::NullEvaluator model; + + Ida ida(&model); + ida.setInitType(AnalysisManager::Sundials::IdaInitType::Y); + ida.configureSimulation(); + ida.initializeSimulation(0.0); + + success *= (std::abs(model.y().getData()[0]) < 1.0e-10); + } + + 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 e4f25a296..ee4723e84 100644 --- a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp +++ b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp @@ -13,6 +13,7 @@ int main() result += test.dtMonitorSuppressesEpsilonFinalStep(); result += test.fixedStep(); result += test.suppressAlgebraicErrors(); + result += test.initType(); return result.summary(); } From 386ba63b20601792d028d6b3a0b13a357b8fb156 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 18 Aug 2026 13:32:29 -0700 Subject: [PATCH 2/7] Remove functions to be more descriptive --- GridKit/Solver/Dynamic/Ida.cpp | 25 ++++++++++--------- GridKit/Solver/Dynamic/Ida.hpp | 8 +++--- .../PhasorDynamics/AnalysisUtilities.hpp | 19 +++++++------- .../PhasorDynamics/ContingencyAnalysis.cpp | 2 +- .../PhasorDynamics/DynamicSimulation.cpp | 2 +- application/PhasorDynamics/README.md | 2 +- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 14 +++++------ .../UnitTests/Solver/Dynamic/runIdaTests.cpp | 2 +- 8 files changed, 38 insertions(+), 36 deletions(-) diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index 2a71b61fe..48394cd7e 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -238,8 +238,8 @@ namespace AnalysisManager // Find a consistent set of initial conditions for DAE if (findConsistent) { - const int initType = getIDAInitType(); - retval = IDACalcIC(solver_, initType, t0 + 0.1); + const int consistentICType = getIDAConsistentICType(); + retval = IDACalcIC(solver_, consistentICType, t0 + 0.1); checkOutput(retval, "IDACalcIC"); retval = IDAGetConsistentIC(solver_, yy_, yp_); @@ -304,20 +304,20 @@ namespace AnalysisManager } template - int Ida::getIDAInitType() const + int Ida::getIDAConsistentICType() const { - switch (init_type_) + switch (consistent_ic_type_) { - case IdaInitType::Y: + case IdaConsistentICType::Y: return IDA_Y_INIT; - case IdaInitType::YA_YDP: + case IdaConsistentICType::YA_YDP: return IDA_YA_YDP_INIT; - case IdaInitType::AUTO: + case IdaConsistentICType::AUTO: return tag_ ? IDA_YA_YDP_INIT : IDA_Y_INIT; default: GridKit::Utilities::Logger::error() - << "Invalid IDA init type " - << static_cast(init_type_) + << "Invalid IDA consistent initial condition type " + << static_cast(consistent_ic_type_) << ".\n"; throw SundialsException(); } @@ -1216,14 +1216,15 @@ namespace AnalysisManager /** * @brief Set the IDA consistent-initial-condition calculation type * - * @param init_type IDA init type. AUTO preserves the existing behavior. + * @param consistent_ic_type IDA consistent initial condition type. AUTO + * preserves the existing behavior. * @tparam ScalarT Scalar data type * @tparam IdxT Index data type */ template - void Ida::setInitType(IdaInitType init_type) + void Ida::setConsistentICType(IdaConsistentICType consistent_ic_type) { - init_type_ = init_type; + consistent_ic_type_ = consistent_ic_type; } /** diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index c214f8951..e024682e8 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -40,7 +40,7 @@ namespace AnalysisManager std::string report() const; }; - enum class IdaInitType + enum class IdaConsistentICType { AUTO, Y, @@ -140,7 +140,7 @@ namespace AnalysisManager ScalarT abs_tol_override = 0); void setSuppressAlgebraicErrors(bool suppress); void setBackwardSuppressAlgebraicErrors(bool suppress); - void setInitType(IdaInitType init_type); + void setConsistentICType(IdaConsistentICType consistent_ic_type); void setMaxSteps(IdxT maxSteps) override; void setBackwardMaxSteps(IdxT maxSteps); @@ -188,7 +188,7 @@ namespace AnalysisManager int getMonitorStepCount(RealT tf, RealT dt_monitor) const; RealT getMonitorTime(RealT tf, RealT dt_monitor, int step, int nsteps) const; - int getIDAInitType() const; + int getIDAConsistentICType() const; void updateModelState(RealT t); private: @@ -222,7 +222,7 @@ namespace AnalysisManager RealT abs_tol_override_{}; IdxT max_steps_{}; bool suppress_alg_{false}; - IdaInitType init_type_{IdaInitType::AUTO}; + IdaConsistentICType consistent_ic_type_{IdaConsistentICType::AUTO}; RealT backward_time_step_{}; RealT backward_rel_tol_{DEFAULT_REL_TOL}; diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index e4643cf65..3b4761f9f 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -65,7 +65,7 @@ namespace GridKit /// maximum number of solver time steps, or 0 for the IDA default std::size_t max_steps; /// IDA consistent initial condition calculation type - AnalysisManager::Sundials::IdaInitType init_type; + AnalysisManager::Sundials::IdaConsistentICType consistent_ic_type; /// set of system events std::vector events; /// path to output file @@ -103,21 +103,22 @@ namespace GridKit c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); c.dt_fixed = j.value("dt_fixed", 0.0); c.max_steps = j.value("max_steps", std::size_t{0}); - c.init_type = AnalysisManager::Sundials::IdaInitType::AUTO; - if (j.contains("init_type")) + c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::AUTO; + if (j.contains("consistent_ic_type")) { - const auto init_type_str = j.at("init_type").get(); - if (init_type_str == "y") + const auto consistent_ic_type_str = j.at("consistent_ic_type").get(); + if (consistent_ic_type_str == "y") { - c.init_type = AnalysisManager::Sundials::IdaInitType::Y; + c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::Y; } - else if (init_type_str == "ya_ydp") + else if (consistent_ic_type_str == "ya_ydp") { - c.init_type = AnalysisManager::Sundials::IdaInitType::YA_YDP; + c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP; } else { - Log::error() << "Invalid IDA init type \"" << init_type_str << "\"; " + Log::error() << "Invalid IDA consistent initial condition type \"" + << consistent_ic_type_str << "\"; " << "must be either \"y\" or \"ya_ydp\""; } } diff --git a/application/PhasorDynamics/ContingencyAnalysis.cpp b/application/PhasorDynamics/ContingencyAnalysis.cpp index 18978d091..01d20821a 100644 --- a/application/PhasorDynamics/ContingencyAnalysis.cpp +++ b/application/PhasorDynamics/ContingencyAnalysis.cpp @@ -38,7 +38,7 @@ TestStatus runStudy(StudyData study_data) ida.setTolerance(study_data.rel_tol, study_data.abs_tol); ida.setFixedStep(study_data.dt_fixed); ida.setMaxSteps(study_data.max_steps); - ida.setInitType(study_data.init_type); + ida.setConsistentICType(study_data.consistent_ic_type); ida.configureSimulation(); using EventType = SystemEvent::Type; diff --git a/application/PhasorDynamics/DynamicSimulation.cpp b/application/PhasorDynamics/DynamicSimulation.cpp index a80613498..52c59e228 100644 --- a/application/PhasorDynamics/DynamicSimulation.cpp +++ b/application/PhasorDynamics/DynamicSimulation.cpp @@ -32,7 +32,7 @@ int main(int argc, const char* argv[]) ida.setTolerance(study.rel_tol, study.abs_tol); ida.setFixedStep(study.dt_fixed); ida.setMaxSteps(study.max_steps); - ida.setInitType(study.init_type); + ida.setConsistentICType(study.consistent_ic_type); ida.configureSimulation(); // Start timer diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index b5b24f8d6..4d3a3c8d8 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -11,7 +11,7 @@ `abs_tol` | Absolute solver tolerance override (default: 1.0e-9) `dt_fixed` | Fixed solver time step size, or 0 for adaptive stepping (default: 0) `max_steps` | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0) - `init_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: automatic) + `consistent_ic_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: automatic) `events` | An array of event groups (see [Events](#events) below) `output_file` | Path to output (CSV) file (optional) `reference_file` | A string containing the name of the case (optional) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index f998842e7..0e02f0800 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -371,7 +371,7 @@ namespace GridKit }; template - class ConsistentInitTypeEvaluator : public NullEvaluator + class ConsistentICTypeEvaluator : public NullEvaluator { protected: using NullEvaluator::allocated_; @@ -382,9 +382,9 @@ namespace GridKit using NullEvaluator::f_; public: - ConsistentInitTypeEvaluator() = default; + ConsistentICTypeEvaluator() = default; - explicit ConsistentInitTypeEvaluator(bool consistent_derivative) + explicit ConsistentICTypeEvaluator(bool consistent_derivative) : consistent_derivative_(consistent_derivative) { } @@ -572,15 +572,15 @@ namespace GridKit return success.report(__func__); } - TestOutcome initType() + TestOutcome consistentICType() { TestStatus success = true; { - Model::ConsistentInitTypeEvaluator model; + Model::ConsistentICTypeEvaluator model; Ida ida(&model); - ida.setInitType(AnalysisManager::Sundials::IdaInitType::YA_YDP); + ida.setConsistentICType(AnalysisManager::Sundials::IdaConsistentICType::YA_YDP); ida.configureSimulation(); ida.initializeSimulation(0.0); @@ -592,7 +592,7 @@ namespace GridKit Model::NullEvaluator model; Ida ida(&model); - ida.setInitType(AnalysisManager::Sundials::IdaInitType::Y); + ida.setConsistentICType(AnalysisManager::Sundials::IdaConsistentICType::Y); ida.configureSimulation(); ida.initializeSimulation(0.0); diff --git a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp index ee4723e84..6b127ce1a 100644 --- a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp +++ b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp @@ -13,7 +13,7 @@ int main() result += test.dtMonitorSuppressesEpsilonFinalStep(); result += test.fixedStep(); result += test.suppressAlgebraicErrors(); - result += test.initType(); + result += test.consistentICType(); return result.summary(); } From 9bffdbd8a282930721d59b487c24dd0f6d2c9a05 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 18 Aug 2026 13:39:55 -0700 Subject: [PATCH 3/7] Remove AUTO --- GridKit/Solver/Dynamic/Ida.cpp | 5 +---- GridKit/Solver/Dynamic/Ida.hpp | 3 +-- application/PhasorDynamics/AnalysisUtilities.hpp | 2 +- application/PhasorDynamics/README.md | 2 +- 4 files changed, 4 insertions(+), 8 deletions(-) diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index 48394cd7e..53c6851ac 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -312,8 +312,6 @@ namespace AnalysisManager return IDA_Y_INIT; case IdaConsistentICType::YA_YDP: return IDA_YA_YDP_INIT; - case IdaConsistentICType::AUTO: - return tag_ ? IDA_YA_YDP_INIT : IDA_Y_INIT; default: GridKit::Utilities::Logger::error() << "Invalid IDA consistent initial condition type " @@ -1216,8 +1214,7 @@ namespace AnalysisManager /** * @brief Set the IDA consistent-initial-condition calculation type * - * @param consistent_ic_type IDA consistent initial condition type. AUTO - * preserves the existing behavior. + * @param consistent_ic_type IDA consistent initial condition type. * @tparam ScalarT Scalar data type * @tparam IdxT Index data type */ diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index e024682e8..d57a308dc 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -42,7 +42,6 @@ namespace AnalysisManager enum class IdaConsistentICType { - AUTO, Y, YA_YDP }; @@ -222,7 +221,7 @@ namespace AnalysisManager RealT abs_tol_override_{}; IdxT max_steps_{}; bool suppress_alg_{false}; - IdaConsistentICType consistent_ic_type_{IdaConsistentICType::AUTO}; + IdaConsistentICType consistent_ic_type_{IdaConsistentICType::YA_YDP}; RealT backward_time_step_{}; RealT backward_rel_tol_{DEFAULT_REL_TOL}; diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index 3b4761f9f..88e2e274a 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -103,7 +103,7 @@ namespace GridKit c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); c.dt_fixed = j.value("dt_fixed", 0.0); c.max_steps = j.value("max_steps", std::size_t{0}); - c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::AUTO; + c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP; if (j.contains("consistent_ic_type")) { const auto consistent_ic_type_str = j.at("consistent_ic_type").get(); diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index 4d3a3c8d8..1d129eccd 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -11,7 +11,7 @@ `abs_tol` | Absolute solver tolerance override (default: 1.0e-9) `dt_fixed` | Fixed solver time step size, or 0 for adaptive stepping (default: 0) `max_steps` | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0) - `consistent_ic_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: automatic) + `consistent_ic_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: "ya_ydp") `events` | An array of event groups (see [Events](#events) below) `output_file` | Path to output (CSV) file (optional) `reference_file` | A string containing the name of the case (optional) From 12d491646e4bcfdb02f39f63149642203fdfd6fe Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 18 Aug 2026 13:50:00 -0700 Subject: [PATCH 4/7] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 59e9c884f..270ec7077 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -80,6 +80,7 @@ - Added `HYGOV` governor model implementation for PhasorDynamics. - Added `REPCA` controller model implementation for PhasorDynamics. - Added `REECB` electrical-control model implementation for PhasorDynamics. +- Added IDA option to choose the consistent initial condition calculation type. ## v0.1 From 53dd44a3beffb0fc11384827ef2a30d613fef65e Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 18 Aug 2026 13:57:31 -0700 Subject: [PATCH 5/7] Align table columns --- application/PhasorDynamics/README.md | 30 ++++++++++++++-------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index 1d129eccd..969979fe0 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -2,22 +2,22 @@ ## Root elements - Name | Value - ---------------------|------------------------------------------------------- - `system_model_file` | Path to the system model file[^1] - `dt_monitor` | Monitor output time interval for recorded simulation results (default: 0, no intermediate monitoring) - `tmax` | A floating-point value for max time - `rel_tol` | Relative solver tolerance (default: 1.0e-7) - `abs_tol` | Absolute solver tolerance override (default: 1.0e-9) - `dt_fixed` | Fixed solver time step size, or 0 for adaptive stepping (default: 0) - `max_steps` | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0) + Name | Value + ----------------------|------------------------------------------------------- + `system_model_file` | Path to the system model file[^1] + `dt_monitor` | Monitor output time interval for recorded simulation results (default: 0, no intermediate monitoring) + `tmax` | A floating-point value for max time + `rel_tol` | Relative solver tolerance (default: 1.0e-7) + `abs_tol` | Absolute solver tolerance override (default: 1.0e-9) + `dt_fixed` | Fixed solver time step size, or 0 for adaptive stepping (default: 0) + `max_steps` | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0) `consistent_ic_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: "ya_ydp") - `events` | An array of event groups (see [Events](#events) below) - `output_file` | Path to output (CSV) file (optional) - `reference_file` | A string containing the name of the case (optional) - `error_type` | One of { "relative" (default), "absolute" } - `error_tolerance` | A floating-point value for highest allowable total error (default: 1.0e-4) - `abs_err_threshold` | A floating-point value for the smallest value at which to scale relative error (default: machine epsilon for double-precision) + `events` | An array of event groups (see [Events](#events) below) + `output_file` | Path to output (CSV) file (optional) + `reference_file` | A string containing the name of the case (optional) + `error_type` | One of { "relative" (default), "absolute" } + `error_tolerance` | A floating-point value for highest allowable total error (default: 1.0e-4) + `abs_err_threshold` | A floating-point value for the smallest value at which to scale relative error (default: machine epsilon for double-precision) [^1]: See system model [case format](../../GridKit/Model/PhasorDynamics/INPUT_FORMAT.md) From bc9385f4d02ce302385f321c344318f20b1f43c8 Mon Sep 17 00:00:00 2001 From: Steven-Roberts Date: Tue, 18 Aug 2026 20:59:57 +0000 Subject: [PATCH 6/7] Apply pre-commit fixes --- GridKit/Solver/Dynamic/Ida.cpp | 2 +- GridKit/Solver/Dynamic/Ida.hpp | 10 +++--- .../PhasorDynamics/AnalysisUtilities.hpp | 36 +++++++++---------- 3 files changed, 24 insertions(+), 24 deletions(-) diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index 53c6851ac..ffeedebb9 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -239,7 +239,7 @@ namespace AnalysisManager if (findConsistent) { const int consistentICType = getIDAConsistentICType(); - retval = IDACalcIC(solver_, consistentICType, t0 + 0.1); + retval = IDACalcIC(solver_, consistentICType, t0 + 0.1); checkOutput(retval, "IDACalcIC"); retval = IDAGetConsistentIC(solver_, yy_, yp_); diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index d57a308dc..bbda253ee 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -216,11 +216,11 @@ namespace AnalysisManager int backwardID_{}; - RealT time_step_{}; - RealT rel_tol_{DEFAULT_REL_TOL}; - RealT abs_tol_override_{}; - IdxT max_steps_{}; - bool suppress_alg_{false}; + RealT time_step_{}; + RealT rel_tol_{DEFAULT_REL_TOL}; + RealT abs_tol_override_{}; + IdxT max_steps_{}; + bool suppress_alg_{false}; IdaConsistentICType consistent_ic_type_{IdaConsistentICType::YA_YDP}; RealT backward_time_step_{}; diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index 88e2e274a..867a73834 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -51,35 +51,35 @@ namespace GridKit struct StudyData { /// path to system model JSON file - fs::path system_model_file; + fs::path system_model_file; /// monitor output time step size, or 0 for no intermediate monitoring - double dt_monitor; + double dt_monitor; /// max time - double tmax; + double tmax; /// relative tolerance for the solver - double rel_tol; + double rel_tol; /// absolute tolerance for the solver - double abs_tol; + double abs_tol; /// fixed solver time step size, or 0 for adaptive stepping - double dt_fixed; + double dt_fixed; /// maximum number of solver time steps, or 0 for the IDA default - std::size_t max_steps; + std::size_t max_steps; /// IDA consistent initial condition calculation type AnalysisManager::Sundials::IdaConsistentICType consistent_ic_type; /// set of system events - std::vector events; + std::vector events; /// path to output file - fs::path output_file; + fs::path output_file; /// path to reference file for validation - fs::path reference_file; + fs::path reference_file; /// Error tolerance (between output file and reference file) - std::vector error_tol; + std::vector error_tol; /// Type of total error (relative or absolute) - Testing::ErrorType error_type; + Testing::ErrorType error_type; /// Smallest value at which to scale for relative error - double abs_err_threshold; + double abs_err_threshold; /// Instance of model data - SystemModelData<> model_data; + SystemModelData<> model_data; }; using json = ::nlohmann::json; @@ -99,10 +99,10 @@ namespace GridKit j.at("system_model_file").get_to(c.system_model_file); c.dt_monitor = j.value("dt_monitor", 0.0); j.at("tmax").get_to(c.tmax); - c.rel_tol = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL); - c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); - c.dt_fixed = j.value("dt_fixed", 0.0); - c.max_steps = j.value("max_steps", std::size_t{0}); + c.rel_tol = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL); + c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); + c.dt_fixed = j.value("dt_fixed", 0.0); + c.max_steps = j.value("max_steps", std::size_t{0}); c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP; if (j.contains("consistent_ic_type")) { From 7ff6b9e9b638dfb39f09371ccecab8d28cd75258 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan <72402802+nkoukpaizan@users.noreply.github.com> Date: Wed, 19 Aug 2026 17:42:04 -0400 Subject: [PATCH 7/7] Refine IDATest::consistentICType test for steady-state initialization (#543) * std::abs--> isEqual * Make all checks explicity in IDATest::consistentICType, including unchanged states. * Add test for IDATest::consistentICType with initially consistent derivative. * Change IDATest::consistentICType model to something that admits a steady-state solution. * Increasing consistentICType tolerances. * Update tolerance comment. --------- Co-authored-by: nkoukpaizan --- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 44 ++++++++++++++++----- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index 0e02f0800..aea551aa5 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -384,8 +384,8 @@ namespace GridKit public: ConsistentICTypeEvaluator() = default; - explicit ConsistentICTypeEvaluator(bool consistent_derivative) - : consistent_derivative_(consistent_derivative) + explicit ConsistentICTypeEvaluator(bool steady_state) + : steady_state_(steady_state) { } @@ -403,7 +403,7 @@ namespace GridKit y[0] = 0.0; y[1] = 10.0; - yp[0] = consistent_derivative_ ? 1.0 : 0.0; + yp[0] = steady_state_ ? 0.0 : 2.0; // Purposefully inconsistent guess for non-steady-state case. yp[1] = 0.0; tag_ = {true, false}; abs_tol[0] = 0.0; @@ -428,14 +428,14 @@ namespace GridKit const auto* y = y_.getData(); const auto* yp = yp_.getData(); - f[0] = yp[0] - 1.0; + f[0] = yp[0] + y[0] + y[1] - 1.0; f[1] = y[1] - y[0]; f_.setDataUpdated(); return 0; } private: - bool consistent_derivative_{false}; + bool steady_state_{false}; }; } // namespace Model @@ -576,16 +576,42 @@ namespace GridKit { TestStatus success = true; + using RealT = typename ScalarTraits::RealT; + + // If the tolerances are too tight, a finite difference approximation to the Jacobian + // cannot be generated, and initialization will fail with bad initial guesses. + static constexpr auto tol = 100.0 * std::numeric_limits::epsilon(); + + // IdaConsistentICType::YA_YDP with non-steady-state initial guess for the derivatives { - Model::ConsistentICTypeEvaluator model; + Model::ConsistentICTypeEvaluator model(false); Ida ida(&model); ida.setConsistentICType(AnalysisManager::Sundials::IdaConsistentICType::YA_YDP); + ida.setTolerance(tol); + ida.configureSimulation(); + ida.initializeSimulation(0.0); + + success *= isEqual(model.yp().getData()[0], 1.0, tol); + success *= isEqual(model.yp().getData()[1], 0.0, tol); + success *= isEqual(model.y().getData()[0], 0.0, tol); + success *= isEqual(model.y().getData()[1], 0.0, tol); + } + + // IdaConsistentICType::Y with steady-state initial guess for the derivatives + { + Model::ConsistentICTypeEvaluator model(true); + + Ida ida(&model); + ida.setConsistentICType(AnalysisManager::Sundials::IdaConsistentICType::Y); + ida.setTolerance(tol); ida.configureSimulation(); ida.initializeSimulation(0.0); - success *= (std::abs(model.yp().getData()[0] - 1.0) < 1.0e-10); - success *= (std::abs(model.y().getData()[1]) < 1.0e-10); + success *= isEqual(model.yp().getData()[0], 0.0, tol); + success *= isEqual(model.yp().getData()[1], 0.0, tol); + success *= isEqual(model.y().getData()[0], 0.5, tol); + success *= isEqual(model.y().getData()[1], 0.5, tol); } { @@ -596,7 +622,7 @@ namespace GridKit ida.configureSimulation(); ida.initializeSimulation(0.0); - success *= (std::abs(model.y().getData()[0]) < 1.0e-10); + success *= isEqual(model.y().getData()[0], 0.0); } return success.report(__func__);