From c9225e0f52e508a3f707faec5042030dc367bf82 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sat, 5 Sep 2026 17:48:53 -0400 Subject: [PATCH 01/12] SystemModel DependencyTracking Jacobian. --- .../PhasorDynamics/Controller/REECB/Reecb.hpp | 1 - .../Controller/REECB/ReecbImpl.hpp | 3 - .../Exciter/ESDC1A/Esdc1aImpl.hpp | 1 - .../Exciter/IEEET1/Ieeet1Impl.hpp | 1 - .../Governor/GASTPTI/GastPtiImpl.hpp | 1 - .../PhasorDynamics/Governor/HYGOV/Hygov.hpp | 1 - .../Governor/HYGOV/HygovImpl.hpp | 2 - .../Governor/Tgov1/Tgov1Impl.hpp | 1 - GridKit/Model/PhasorDynamics/SystemModel.cpp | 10 + GridKit/Model/PhasorDynamics/SystemModel.hpp | 1 - .../SystemModelDependencyTracking.cpp | 131 +++++++++++- .../PhasorDynamics/SystemModelEnzyme.cpp | 187 ++++++++++++++++++ .../Model/PhasorDynamics/SystemModelImpl.hpp | 187 ------------------ 13 files changed, 326 insertions(+), 201 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp index 302361cbe..6ec4bb77c 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp @@ -7,7 +7,6 @@ #pragma once #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp index 0b8fd4615..705cc5d78 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -8,10 +8,7 @@ #include #include -#include -#include #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 2645ed940..f2ed50915 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -7,7 +7,6 @@ #pragma once #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index ada67c9e3..a67560295 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -8,7 +8,6 @@ */ #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp index d5a33b257..a94410c75 100644 --- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp @@ -8,7 +8,6 @@ #include #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp index b11641318..10ed07c26 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -8,7 +8,6 @@ #include #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp index e85d05354..84a82d356 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -7,8 +7,6 @@ #pragma once #include -#include -#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 229a58df3..298ce785f 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -9,7 +9,6 @@ */ #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/SystemModel.cpp b/GridKit/Model/PhasorDynamics/SystemModel.cpp index 6cbce2002..c1dde94f7 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.cpp @@ -17,6 +17,16 @@ namespace GridKit return false; } + /** + * @brief By default, Jacobians are not available + * + */ + template + int SystemModel::evaluateJacobian() + { + return 0; + } + // Available template instantiations // template class SystemModel; template class SystemModel; diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index c39e0894b..dc0ed0b7b 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -2,7 +2,6 @@ #include #include -#include #include #include diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index 063b86d3f..af23e1069 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -7,15 +7,142 @@ namespace GridKit /** * @brief By default, Jacobians are not available * - * DependencyTracking::Variable stores the Jacobian as dependency maps. - * @todo Construct a Jacobian based on the dependency maps. + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). */ template bool SystemModel::hasJacobian() { + Log::warning() << "GridKit was not built with Enzyme. " + << "DependencyTracking::Variable Jacobians are only available for testing in PhasorDynamics.\n"; + return false; } + /** + * @brief Evaluate system DependencyTracking::Variable Jacobian. + * + */ + template + int SystemModel::evaluateJacobian() + { + using DependencyMap = typename ScalarT::DependencyMap; + + const auto* f = f_.getData(); + + if (csr_jac_ == nullptr) + { + IdxT* row_ptrs = new IdxT[static_cast(size_) + 1]; + row_ptrs[0] = 0; + + // Count the number of non-zeros + IdxT nnz = 0; + for (IdxT row = 0; row < size_; ++row) + { + DependencyMap row_map; + + for (const auto& dep : f[row].getDependencies()) + { + const IdxT col = static_cast(dep.first); + + // Merge-count y and yp dependencies + IdxT jac_col; + if (col < size_) + { + jac_col = col; + } + else + { + jac_col = col - size_; + } + + if (row_map.insert({jac_col, RealT{}}).second) + { + ++nnz; + } + } + + row_ptrs[static_cast(row) + 1] = nnz; + } + + // Allocate column and value pointers + IdxT* cols = new IdxT[static_cast(nnz)]; + RealT* vals = new RealT[static_cast(nnz)]; + + // Store column and values + IdxT i = 0; + for (IdxT row = 0; row < size_; ++row) + { + DependencyMap row_map; + + for (const auto& dep : f[row].getDependencies()) + { + const IdxT col = static_cast(dep.first); + + IdxT jac_col; + if (col < size_) + { + jac_col = col; + row_map[jac_col] += static_cast(dep.second); + } + else + { + jac_col = col - size_; + row_map[jac_col] += alpha_ * static_cast(dep.second); + } + } + + for (const auto& entry : row_map) + { + cols[i] = static_cast(entry.first); + vals[i] = static_cast(entry.second); + ++i; + } + } + + nnz_ = nnz; + csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); + } + else + { + RealT* vals = csr_jac_->getValues(); + + IdxT i = 0; + for (IdxT row = 0; row < size_; ++row) + { + DependencyMap row_map; + + for (const auto& dep : f[row].getDependencies()) + { + const IdxT col = static_cast(dep.first); + + IdxT jac_col; + if (col < size_) + { + jac_col = col; + row_map[jac_col] += static_cast(dep.second); + } + else + { + jac_col = col - size_; + row_map[jac_col] += alpha_ * static_cast(dep.second); + } + } + + for (const auto& entry : row_map) + { + vals[i] = static_cast(entry.second); + ++i; + } + } + } + + //Log::misc() << "System DependencyTracking Jacobian\n"; + //csr_jac_->print(Log::misc()); + + return 0; + } + // Available template instantiations // template class SystemModel; template class SystemModel; diff --git a/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp b/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp index 6d20c9c67..08c434a8c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp @@ -34,6 +34,193 @@ namespace GridKit return has_jacobian; } + /** + * @brief Evaluate system Jacobian using component-level Jacobians. + * + * - Evaluate component-level Jacobians (stored as CooMatrixT). + * - If not already constructed, construct system-level Jacobian. + * This used a system-level CooMatrix to deduplicate and sort the sparsity pattern. + * - If already constructed, reset Jacobian values to zero and accumulate contributions. + * + */ + template + int SystemModel::evaluateJacobian() + { + // Initialize bus Jacobians + for (const auto& bus : buses_) + { + bus->evaluateJacobian(); + } + + // Evaluate component Jacobians, including contribution to the bus Jacobians + for (const auto& component : components_) + { + component->evaluateJacobian(); + } + + // Build or update system CSR Jacobian + if (csr_jac_ == nullptr) + { + // Count the number of non-zeros + IdxT nnz_dup = 0; + for (const auto& component : components_) + { + auto component_jacobian = component->getCooJacobian(); + + if (component_jacobian != nullptr) + { + nnz_dup += component_jacobian->getNnz(); + } + else + { + Log::warning() << "A component has returned a nullptr Jacobian.\n"; + } + } + + for (const auto& bus : buses_) + { + auto bus_jacobian = bus->getCooJacobian(); + + if (bus_jacobian != nullptr) + { + nnz_dup += bus_jacobian->getNnz(); + } + else + { + Log::warning() << "A bus has returned a nullptr Jacobian.\n"; + } + } + + // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix) + IdxT* rows_dup = new IdxT[static_cast(nnz_dup)]; + IdxT* cols_dup = new IdxT[static_cast(nnz_dup)]; + RealT* vals_dup = new RealT[static_cast(nnz_dup)]; + + IdxT counter = 0; + for (const auto& component : components_) + { + auto component_jacobian = component->getCooJacobian(); + + if (component_jacobian != nullptr) + { + const IdxT* rows = component_jacobian->getRowData(); + const IdxT* columns = component_jacobian->getColData(); + const RealT* values = component_jacobian->getValues(); + for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + { + rows_dup[counter] = rows[i]; + cols_dup[counter] = columns[i]; + vals_dup[counter] = values[i]; + counter++; + } + } + else + { + Log::warning() << "A component has returned a nullptr Jacobian.\n"; + } + } + + for (const auto& bus : buses_) + { + auto bus_jacobian = bus->getCooJacobian(); + + if (bus_jacobian != nullptr) + { + const IdxT* rows = bus_jacobian->getRowData(); + const IdxT* columns = bus_jacobian->getColData(); + const RealT* values = bus_jacobian->getValues(); + for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + { + rows_dup[counter] = rows[i]; + cols_dup[counter] = columns[i]; + vals_dup[counter] = values[i]; + counter++; + } + } + else + { + Log::warning() << "A bus has returned a nullptr Jacobian.\n"; + } + } + + // Build the system COO Jacobian + CooMatrixT jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); + + // Populate CSR data with sort and deduplicate + IdxT* row_ptrs = jac.getCsrRowData(); + + // Deduplicated nnz + nnz_ = jac.getNnz(); + + // Allocate cols/vals with deduplicated nnz + IdxT* cols = new IdxT[static_cast(nnz_)]; + RealT* vals = new RealT[static_cast(nnz_)]; + + std::copy(jac.getColData(), jac.getColData() + nnz_, cols); + std::copy(jac.getValues(), jac.getValues() + nnz_, vals); + + // Create the CSR Jacobian + csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); + + const IdxT* map_to_sorted = jac.getMapToSorted(); + const IdxT* map_to_dedup = jac.getMapToDeduplicated(); + + // Build a mappping from original COO index to CSR index + map_to_csr_ = new IdxT[static_cast(nnz_dup)]; + for (IdxT i = 0; i < nnz_dup; ++i) + { + map_to_csr_[map_to_sorted[i]] = map_to_dedup[i]; + } + } + else + { + // Zero out values + RealT* vals = csr_jac_->getValues(); + for (IdxT i = 0; i < csr_jac_->getNnz(); ++i) + { + vals[i] = 0.0; + } + + // Update CSR values from component and bus Jacobians + IdxT counter = 0; + for (const auto& component : components_) + { + auto component_jacobian = component->getCooJacobian(); + + if (component_jacobian != nullptr) + { + const RealT* values = component_jacobian->getValues(); + for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + { + vals[map_to_csr_[counter]] += values[i]; + counter++; + } + } + } + + for (const auto& bus : buses_) + { + auto bus_jacobian = bus->getCooJacobian(); + + if (bus_jacobian != nullptr) + { + const RealT* values = bus_jacobian->getValues(); + for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + { + vals[map_to_csr_[counter]] += values[i]; + counter++; + } + } + } + } + + //Log::misc() << "System DependencyTracking Jacobian\n"; + //csr_jac_->print(Log::misc()); + + return 0; + } + + // Available template instantiations // template class SystemModel; template class SystemModel; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 2abc4f153..f0834cd88 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -703,193 +703,6 @@ namespace GridKit return 0; } - /** - * @brief Evaluate system Jacobian. - * - * First, initialize bus Jacobians to 0. - * Then, evaluate component Jacobians (internal block and bus Jacobian contributions). - * Once component Jacobians are evaluated, store the result in the system Jacobian. - * Finally, store bus Jacobians into the system Jacobian after all component have added their - * contributions. - * - */ - template - int SystemModel::evaluateJacobian() - { - // Initialize bus Jacobians - for (const auto& bus : buses_) - { - bus->evaluateJacobian(); - } - - // Evaluate component Jacobians, including contribution to the bus Jacobians - for (const auto& component : components_) - { - component->evaluateJacobian(); - } - - // Build or update system CSR Jacobian - if (csr_jac_ == nullptr) - { - // Count the number of non-zeros - IdxT nnz_dup = 0; - for (const auto& component : components_) - { - auto component_jacobian = component->getCooJacobian(); - - if (component_jacobian != nullptr) - { - nnz_dup += component_jacobian->getNnz(); - } - else - { - Log::warning() << "A component has returned a nullptr Jacobian.\n"; - } - } - - for (const auto& bus : buses_) - { - auto bus_jacobian = bus->getCooJacobian(); - - if (bus_jacobian != nullptr) - { - nnz_dup += bus_jacobian->getNnz(); - } - else - { - Log::warning() << "A bus has returned a nullptr Jacobian.\n"; - } - } - - // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix) - IdxT* rows_dup = new IdxT[static_cast(nnz_dup)]; - IdxT* cols_dup = new IdxT[static_cast(nnz_dup)]; - RealT* vals_dup = new RealT[static_cast(nnz_dup)]; - - IdxT counter = 0; - for (const auto& component : components_) - { - auto component_jacobian = component->getCooJacobian(); - - if (component_jacobian != nullptr) - { - const IdxT* rows = component_jacobian->getRowData(); - const IdxT* columns = component_jacobian->getColData(); - const RealT* values = component_jacobian->getValues(); - for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) - { - rows_dup[counter] = rows[i]; - cols_dup[counter] = columns[i]; - vals_dup[counter] = values[i]; - counter++; - } - } - else - { - Log::warning() << "A component has returned a nullptr Jacobian.\n"; - } - } - - for (const auto& bus : buses_) - { - auto bus_jacobian = bus->getCooJacobian(); - - if (bus_jacobian != nullptr) - { - const IdxT* rows = bus_jacobian->getRowData(); - const IdxT* columns = bus_jacobian->getColData(); - const RealT* values = bus_jacobian->getValues(); - for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) - { - rows_dup[counter] = rows[i]; - cols_dup[counter] = columns[i]; - vals_dup[counter] = values[i]; - counter++; - } - } - else - { - Log::warning() << "A bus has returned a nullptr Jacobian.\n"; - } - } - - // Build the system COO Jacobian - CooMatrixT jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); - - // Populate CSR data with sort and deduplicate - IdxT* row_ptrs = jac.getCsrRowData(); - - // Deduplicated nnz - nnz_ = jac.getNnz(); - - // Allocate cols/vals with deduplicated nnz - IdxT* cols = new IdxT[static_cast(nnz_)]; - RealT* vals = new RealT[static_cast(nnz_)]; - - std::copy(jac.getColData(), jac.getColData() + nnz_, cols); - std::copy(jac.getValues(), jac.getValues() + nnz_, vals); - - // Create the CSR Jacobian - csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); - - const IdxT* map_to_sorted = jac.getMapToSorted(); - const IdxT* map_to_dedup = jac.getMapToDeduplicated(); - - // Build a mappping from original COO index to CSR index - map_to_csr_ = new IdxT[static_cast(nnz_dup)]; - for (IdxT i = 0; i < nnz_dup; ++i) - { - map_to_csr_[map_to_sorted[i]] = map_to_dedup[i]; - } - } - else - { - // Zero out values - RealT* vals = csr_jac_->getValues(); - for (IdxT i = 0; i < csr_jac_->getNnz(); ++i) - { - vals[i] = 0.0; - } - - // Update CSR values from component and bus Jacobians - IdxT counter = 0; - for (const auto& component : components_) - { - auto component_jacobian = component->getCooJacobian(); - - if (component_jacobian != nullptr) - { - const RealT* values = component_jacobian->getValues(); - for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) - { - vals[map_to_csr_[counter]] += values[i]; - counter++; - } - } - } - - for (const auto& bus : buses_) - { - auto bus_jacobian = bus->getCooJacobian(); - - if (bus_jacobian != nullptr) - { - const RealT* values = bus_jacobian->getValues(); - for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) - { - vals[map_to_csr_[counter]] += values[i]; - counter++; - } - } - } - } - - // std::cout << "System Jacobian\n"; - // csr_jac_->print(std::cout); - - return 0; - } - /** * @brief Update time * From 8c9f8aeb6f9ac0bf1a4e068af6709d9632e571a8 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sat, 5 Sep 2026 19:17:01 -0400 Subject: [PATCH 02/12] Set variableNumber in SystemModel's initialize(). --- GridKit/Model/PhasorDynamics/SystemModel.hpp | 3 +++ .../SystemModelDependencyTracking.cpp | 14 ++++++------- .../Model/PhasorDynamics/SystemModelImpl.hpp | 20 +++++++++++++++++++ 3 files changed, 30 insertions(+), 7 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index dc0ed0b7b..cc4cbdb3c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -117,6 +117,9 @@ namespace GridKit bool owns_components_{false}; + /// Offset between y and yp for DependencyTracking::Variable numbers + IdxT y_yp_offset_; + /// Variable monitor std::unique_ptr monitor_; }; // class SystemModel diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index af23e1069..3eb26e968 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -29,7 +29,7 @@ namespace GridKit using DependencyMap = typename ScalarT::DependencyMap; const auto* f = f_.getData(); - + if (csr_jac_ == nullptr) { IdxT* row_ptrs = new IdxT[static_cast(size_) + 1]; @@ -47,13 +47,13 @@ namespace GridKit // Merge-count y and yp dependencies IdxT jac_col; - if (col < size_) + if (col < y_yp_offset_) { jac_col = col; } else { - jac_col = col - size_; + jac_col = col - y_yp_offset_; } if (row_map.insert({jac_col, RealT{}}).second) @@ -80,14 +80,14 @@ namespace GridKit const IdxT col = static_cast(dep.first); IdxT jac_col; - if (col < size_) + if (col < y_yp_offset_) { jac_col = col; row_map[jac_col] += static_cast(dep.second); } else { - jac_col = col - size_; + jac_col = col - y_yp_offset_; row_map[jac_col] += alpha_ * static_cast(dep.second); } } @@ -117,14 +117,14 @@ namespace GridKit const IdxT col = static_cast(dep.first); IdxT jac_col; - if (col < size_) + if (col < y_yp_offset_) { jac_col = col; row_map[jac_col] += static_cast(dep.second); } else { - jac_col = col - size_; + jac_col = col - y_yp_offset_; row_map[jac_col] += alpha_ * static_cast(dep.second); } } diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index f0834cd88..61883cf32 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -557,6 +557,26 @@ namespace GridKit status += component->initialize(); } + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + auto* y = y_.getData(); + auto* yp = yp_.getData(); + + // Offset the number for yp variables to track y and yp simultaneously. + // @todo For a hierarchical system, the offsets must be provided by a higher level + y_yp_offset_ = size_; + for (IdxT j = 0; j < size_; ++j) + { + const IdxT var_idx = this->getVariableIndex(j); + if (var_idx != INVALID_INDEX) + { + y[j].setVariableNumber(static_cast(var_idx)); + yp[j].setVariableNumber(static_cast(var_idx) + static_cast(y_yp_offset_)); + } + } + } + y_.setDataUpdated(); yp_.setDataUpdated(); From 3c7539cb20d5c896b37c81f44d99942606e3d652 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sat, 5 Sep 2026 19:39:09 -0400 Subject: [PATCH 03/12] Update SystemTest to get the DependencyTracking Jacobian in CSC form directly. --- .../UnitTests/PhasorDynamics/SystemTests.hpp | 35 ++++--------------- 1 file changed, 7 insertions(+), 28 deletions(-) diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index afc829795..21145ff90 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -494,35 +494,14 @@ namespace GridKit system.allocate(); system.initialize(); - // Set independent variables - auto* y = system.y().getData(); - for (size_t i = 0; i < system.size(); ++i) - { - y[i].setVariableNumber(i); - } - system.y().setDataUpdated(); - - // Evaluate and get the system residuals + // Evaluate and get the system Jacobian system.evaluateResidual(); - auto& residual = system.getResidual(); - const auto* residual_data = residual.getData(); - - // Print the dependencies - for (size_t i = 0; i < residual.getSize(); ++i) - { - std::cout << i << "th residual: "; - residual_data[i].print(std::cout); - std::cout << "\n"; - } - - // Extract the dependencies - std::vector dependencies(residual.getSize()); - for (IdxT i = 0; i < residual.getSize(); ++i) - { - dependencies[i] = residual_data[i].getDependencies(); - } + system.evaluateJacobian(); + GridKit::LinearAlgebra::CsrMatrix* system_jacobian = system.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: System Jacobian with DependencyTracking\n"; + system_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(system_jacobian); } std::vector EnzymeJacobian( @@ -539,7 +518,7 @@ namespace GridKit system.evaluateResidual(); system.evaluateJacobian(); GridKit::LinearAlgebra::CsrMatrix* system_jacobian = system.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: System Jacobian\n"; + std::cout << "Sparse Csr Matrix: System Jacobian with Enzyme\n"; system_jacobian->print(); return GridKit::Testing::MapFromCsr(system_jacobian); From a7ccb2ccc22544bf7797b04adb0b01f3ef33f8db Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sat, 5 Sep 2026 19:44:23 -0400 Subject: [PATCH 04/12] Minor variable name change. --- GridKit/Model/PhasorDynamics/SystemModel.hpp | 2 +- .../PhasorDynamics/SystemModelDependencyTracking.cpp | 12 ++++++------ GridKit/Model/PhasorDynamics/SystemModelImpl.hpp | 4 ++-- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index cc4cbdb3c..22133be79 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -118,7 +118,7 @@ namespace GridKit bool owns_components_{false}; /// Offset between y and yp for DependencyTracking::Variable numbers - IdxT y_yp_offset_; + IdxT y_yp_tracking_offset_; /// Variable monitor std::unique_ptr monitor_; diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index 3eb26e968..ce4f55ffe 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -47,13 +47,13 @@ namespace GridKit // Merge-count y and yp dependencies IdxT jac_col; - if (col < y_yp_offset_) + if (col < y_yp_tracking_offset_) { jac_col = col; } else { - jac_col = col - y_yp_offset_; + jac_col = col - y_yp_tracking_offset_; } if (row_map.insert({jac_col, RealT{}}).second) @@ -80,14 +80,14 @@ namespace GridKit const IdxT col = static_cast(dep.first); IdxT jac_col; - if (col < y_yp_offset_) + if (col < y_yp_tracking_offset_) { jac_col = col; row_map[jac_col] += static_cast(dep.second); } else { - jac_col = col - y_yp_offset_; + jac_col = col - y_yp_tracking_offset_; row_map[jac_col] += alpha_ * static_cast(dep.second); } } @@ -117,14 +117,14 @@ namespace GridKit const IdxT col = static_cast(dep.first); IdxT jac_col; - if (col < y_yp_offset_) + if (col < y_yp_tracking_offset_) { jac_col = col; row_map[jac_col] += static_cast(dep.second); } else { - jac_col = col - y_yp_offset_; + jac_col = col - y_yp_tracking_offset_; row_map[jac_col] += alpha_ * static_cast(dep.second); } } diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 61883cf32..fe8024165 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -565,14 +565,14 @@ namespace GridKit // Offset the number for yp variables to track y and yp simultaneously. // @todo For a hierarchical system, the offsets must be provided by a higher level - y_yp_offset_ = size_; + y_yp_tracking_offset_ = size_; for (IdxT j = 0; j < size_; ++j) { const IdxT var_idx = this->getVariableIndex(j); if (var_idx != INVALID_INDEX) { y[j].setVariableNumber(static_cast(var_idx)); - yp[j].setVariableNumber(static_cast(var_idx) + static_cast(y_yp_offset_)); + yp[j].setVariableNumber(static_cast(var_idx) + static_cast(y_yp_tracking_offset_)); } } } From 11db87ec3b21a7b1f0e66f1565c12e06afafb646 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 7 Sep 2026 16:14:27 -0400 Subject: [PATCH 05/12] Switch from y_yp_tracking_offset to odd/even split. --- GridKit/Model/PhasorDynamics/SystemModel.hpp | 3 -- .../SystemModelDependencyTracking.cpp | 32 +++++++------------ .../Model/PhasorDynamics/SystemModelImpl.hpp | 9 ++---- 3 files changed, 14 insertions(+), 30 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 22133be79..dc0ed0b7b 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -117,9 +117,6 @@ namespace GridKit bool owns_components_{false}; - /// Offset between y and yp for DependencyTracking::Variable numbers - IdxT y_yp_tracking_offset_; - /// Variable monitor std::unique_ptr monitor_; }; // class SystemModel diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index ce4f55ffe..b044cfb19 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -29,7 +29,7 @@ namespace GridKit using DependencyMap = typename ScalarT::DependencyMap; const auto* f = f_.getData(); - + if (csr_jac_ == nullptr) { IdxT* row_ptrs = new IdxT[static_cast(size_) + 1]; @@ -43,18 +43,10 @@ namespace GridKit for (const auto& dep : f[row].getDependencies()) { - const IdxT col = static_cast(dep.first); + const auto col = dep.first; // Merge-count y and yp dependencies - IdxT jac_col; - if (col < y_yp_tracking_offset_) - { - jac_col = col; - } - else - { - jac_col = col - y_yp_tracking_offset_; - } + const IdxT jac_col = static_cast(col / 2); if (row_map.insert({jac_col, RealT{}}).second) { @@ -77,17 +69,16 @@ namespace GridKit for (const auto& dep : f[row].getDependencies()) { - const IdxT col = static_cast(dep.first); + const auto col = dep.first; - IdxT jac_col; - if (col < y_yp_tracking_offset_) + const IdxT jac_col = static_cast(col / 2); + // Even indices for y and odd indices for yp + if (col % 2 == 0) { - jac_col = col; row_map[jac_col] += static_cast(dep.second); } else { - jac_col = col - y_yp_tracking_offset_; row_map[jac_col] += alpha_ * static_cast(dep.second); } } @@ -114,17 +105,16 @@ namespace GridKit for (const auto& dep : f[row].getDependencies()) { - const IdxT col = static_cast(dep.first); + const auto col = dep.first; - IdxT jac_col; - if (col < y_yp_tracking_offset_) + const IdxT jac_col = static_cast(col / 2); + // Even indices for y and odd indices for yp + if (col % 2 == 0) { - jac_col = col; row_map[jac_col] += static_cast(dep.second); } else { - jac_col = col - y_yp_tracking_offset_; row_map[jac_col] += alpha_ * static_cast(dep.second); } } diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index fe8024165..fda98a3e5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -10,7 +10,6 @@ // Include all components #include - namespace GridKit { namespace PhasorDynamics @@ -563,16 +562,14 @@ namespace GridKit auto* y = y_.getData(); auto* yp = yp_.getData(); - // Offset the number for yp variables to track y and yp simultaneously. - // @todo For a hierarchical system, the offsets must be provided by a higher level - y_yp_tracking_offset_ = size_; for (IdxT j = 0; j < size_; ++j) { const IdxT var_idx = this->getVariableIndex(j); if (var_idx != INVALID_INDEX) { - y[j].setVariableNumber(static_cast(var_idx)); - yp[j].setVariableNumber(static_cast(var_idx) + static_cast(y_yp_tracking_offset_)); + // Even indices for y and odd indices for yp + y[j].setVariableNumber(static_cast(2 * var_idx)); + yp[j].setVariableNumber(static_cast(2 * var_idx + 1)); } } } From 054e0b55ec3adaada7e069bc3b0b432f9f0c5b6d Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 7 Sep 2026 19:24:27 -0400 Subject: [PATCH 06/12] Made constructCsrFromCoo protected and added a public constructCsr that will dispatch based on ScalarT. --- GridKit/Model/PhasorDynamics/Component.hpp | 67 +++++++++++++++------- 1 file changed, 46 insertions(+), 21 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 707902dfc..4cccfd095 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -261,29 +261,21 @@ namespace GridKit return gridkit_component_id_; } + /** + * @brief CSR construction dispatch depending on ScalarT + * + * @note Currently only used for testing purposes. + */ int constructCsr() { - if (coo_jac_ == nullptr) - { - constructCoo(); - } - - if (csr_jac_ == nullptr) - { - IdxT* row_ptrs = coo_jac_->getCsrRowData(); - - nnz_ = coo_jac_->getNnz(); - - IdxT* cols = new IdxT[static_cast(nnz_)]; - RealT* vals = new RealT[static_cast(nnz_)]; - - std::copy(coo_jac_->getColData(), coo_jac_->getColData() + nnz_, cols); - std::copy(coo_jac_->getValues(), coo_jac_->getValues() + nnz_, vals); - - csr_jac_ = new CsrMatrixT(coo_jac_->getNumRows(), coo_jac_->getNumColumns(), nnz_, &row_ptrs, &cols, &vals); - } - - return 0; + //if constexpr (std::is_same_v) + //{ + // return constructCsrFromDependencies(); + //} + //else + //{ + return constructCsrFromCoo(); + //} } protected: @@ -316,6 +308,9 @@ namespace GridKit abs_tol_.resize(n); } + /** + * @brief COO construction from raw buffers. + */ int constructCoo() { if (coo_jac_ == nullptr) @@ -340,6 +335,36 @@ namespace GridKit return 0; } + /** + * @brief CSR construction from COO. + * + * @note Currently only used for testing purposes. + */ + int constructCsrFromCoo() + { + if (coo_jac_ == nullptr) + { + constructCoo(); + } + + if (csr_jac_ == nullptr) + { + IdxT* row_ptrs = coo_jac_->getCsrRowData(); + + nnz_ = coo_jac_->getNnz(); + + IdxT* cols = new IdxT[static_cast(nnz_)]; + RealT* vals = new RealT[static_cast(nnz_)]; + + std::copy(coo_jac_->getColData(), coo_jac_->getColData() + nnz_, cols); + std::copy(coo_jac_->getValues(), coo_jac_->getValues() + nnz_, vals); + + csr_jac_ = new CsrMatrixT(coo_jac_->getNumRows(), coo_jac_->getNumColumns(), nnz_, &row_ptrs, &cols, &vals); + } + + return 0; + } + IdxT size_{0}; IdxT nnz_{0}; /// Global (system-level) variable indices From 617dc1ab42c6df553fb2f1a8c5c206e000d61876 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 7 Sep 2026 19:40:50 -0400 Subject: [PATCH 07/12] Move constructCsrFromDependencies to base class. --- GridKit/Model/PhasorDynamics/Component.hpp | 133 ++++++++++++++++-- .../SystemModelDependencyTracking.cpp | 103 +------------- 2 files changed, 126 insertions(+), 110 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 4cccfd095..1a50d13c1 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -268,14 +268,14 @@ namespace GridKit */ int constructCsr() { - //if constexpr (std::is_same_v) - //{ - // return constructCsrFromDependencies(); - //} - //else - //{ + if constexpr (std::is_same_v) + { + return constructCsrFromDependencies(); + } + else + { return constructCsrFromCoo(); - //} + } } protected: @@ -301,7 +301,6 @@ namespace GridKit */ void allocateVectors(IdxT n) { - y_.resize(n); yp_.resize(n); f_.resize(n); @@ -309,7 +308,9 @@ namespace GridKit } /** - * @brief COO construction from raw buffers. + * @brief COO construction from component-level raw buffers. + * + * @note the components retain ownership of the data in the raw buffers. */ int constructCoo() { @@ -365,6 +366,120 @@ namespace GridKit return 0; } + /** + * @brief CSR construction from Dependency maps. + * + * @note Currently only used for testing purposes. + */ + int constructCsrFromDependencies() + { + static_assert(std::is_same_v, + "constructCsrFromDependencies() requires ScalarT = DependencyTracking::Variable"); + + using DependencyMap = typename ScalarT::DependencyMap; + + const auto* f = f_.getData(); + + if (csr_jac_ == nullptr) + { + IdxT* row_ptrs = new IdxT[static_cast(size_) + 1]; + row_ptrs[0] = 0; + + // Count the number of non-zeros + IdxT nnz = 0; + for (IdxT row = 0; row < size_; ++row) + { + DependencyMap row_map; + + for (const auto& dep : f[row].getDependencies()) + { + const auto col = dep.first; + + // Merge-count y and yp dependencies + const IdxT jac_col = static_cast(col / 2); + + if (row_map.insert({jac_col, RealT{}}).second) + { + ++nnz; + } + } + + row_ptrs[static_cast(row) + 1] = nnz; + } + + // Allocate column and value pointers + IdxT* cols = new IdxT[static_cast(nnz)]; + RealT* vals = new RealT[static_cast(nnz)]; + + // Store column and values + IdxT i = 0; + for (IdxT row = 0; row < size_; ++row) + { + DependencyMap row_map; + + for (const auto& dep : f[row].getDependencies()) + { + const auto col = dep.first; + + const IdxT jac_col = static_cast(col / 2); + // Even indices for y and odd indices for yp + if (col % 2 == 0) + { + row_map[jac_col] += static_cast(dep.second); + } + else + { + row_map[jac_col] += alpha_ * static_cast(dep.second); + } + } + + for (const auto& entry : row_map) + { + cols[i] = static_cast(entry.first); + vals[i] = static_cast(entry.second); + ++i; + } + } + + nnz_ = nnz; + csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); + } + else + { + RealT* vals = csr_jac_->getValues(); + + IdxT i = 0; + for (IdxT row = 0; row < size_; ++row) + { + DependencyMap row_map; + + for (const auto& dep : f[row].getDependencies()) + { + const auto col = dep.first; + + const IdxT jac_col = static_cast(col / 2); + // Even indices for y and odd indices for yp + if (col % 2 == 0) + { + row_map[jac_col] += static_cast(dep.second); + } + else + { + row_map[jac_col] += alpha_ * static_cast(dep.second); + } + } + + for (const auto& entry : row_map) + { + vals[i] = static_cast(entry.second); + ++i; + } + } + } + + return 0; + } + IdxT size_{0}; IdxT nnz_{0}; /// Global (system-level) variable indices diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index b044cfb19..447fd61cb 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -25,107 +25,8 @@ namespace GridKit */ template int SystemModel::evaluateJacobian() - { - using DependencyMap = typename ScalarT::DependencyMap; - - const auto* f = f_.getData(); - - if (csr_jac_ == nullptr) - { - IdxT* row_ptrs = new IdxT[static_cast(size_) + 1]; - row_ptrs[0] = 0; - - // Count the number of non-zeros - IdxT nnz = 0; - for (IdxT row = 0; row < size_; ++row) - { - DependencyMap row_map; - - for (const auto& dep : f[row].getDependencies()) - { - const auto col = dep.first; - - // Merge-count y and yp dependencies - const IdxT jac_col = static_cast(col / 2); - - if (row_map.insert({jac_col, RealT{}}).second) - { - ++nnz; - } - } - - row_ptrs[static_cast(row) + 1] = nnz; - } - - // Allocate column and value pointers - IdxT* cols = new IdxT[static_cast(nnz)]; - RealT* vals = new RealT[static_cast(nnz)]; - - // Store column and values - IdxT i = 0; - for (IdxT row = 0; row < size_; ++row) - { - DependencyMap row_map; - - for (const auto& dep : f[row].getDependencies()) - { - const auto col = dep.first; - - const IdxT jac_col = static_cast(col / 2); - // Even indices for y and odd indices for yp - if (col % 2 == 0) - { - row_map[jac_col] += static_cast(dep.second); - } - else - { - row_map[jac_col] += alpha_ * static_cast(dep.second); - } - } - - for (const auto& entry : row_map) - { - cols[i] = static_cast(entry.first); - vals[i] = static_cast(entry.second); - ++i; - } - } - - nnz_ = nnz; - csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); - } - else - { - RealT* vals = csr_jac_->getValues(); - - IdxT i = 0; - for (IdxT row = 0; row < size_; ++row) - { - DependencyMap row_map; - - for (const auto& dep : f[row].getDependencies()) - { - const auto col = dep.first; - - const IdxT jac_col = static_cast(col / 2); - // Even indices for y and odd indices for yp - if (col % 2 == 0) - { - row_map[jac_col] += static_cast(dep.second); - } - else - { - row_map[jac_col] += alpha_ * static_cast(dep.second); - } - } - - for (const auto& entry : row_map) - { - vals[i] = static_cast(entry.second); - ++i; - } - } - } + { + this->constructCsr(); //Log::misc() << "System DependencyTracking Jacobian\n"; //csr_jac_->print(Log::misc()); From 02d9430caa3d33708241e349bb540a4ea76a1f51 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 7 Sep 2026 20:15:15 -0400 Subject: [PATCH 08/12] Move initializeDependencyTrackingVariableNumbers to base class. --- GridKit/Model/PhasorDynamics/Component.hpp | 79 +++++++++++++------ .../SystemModelDependencyTracking.cpp | 12 +-- .../PhasorDynamics/SystemModelEnzyme.cpp | 5 +- .../Model/PhasorDynamics/SystemModelImpl.hpp | 25 ++---- 4 files changed, 70 insertions(+), 51 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 1a50d13c1..69b94db8d 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -263,7 +263,7 @@ namespace GridKit /** * @brief CSR construction dispatch depending on ScalarT - * + * * @note Currently only used for testing purposes. */ int constructCsr() @@ -309,7 +309,7 @@ namespace GridKit /** * @brief COO construction from component-level raw buffers. - * + * * @note the components retain ownership of the data in the raw buffers. */ int constructCoo() @@ -370,6 +370,7 @@ namespace GridKit * @brief CSR construction from Dependency maps. * * @note Currently only used for testing purposes. + * See \ref initializeDependencyTrackingVariableNumbers() */ int constructCsrFromDependencies() { @@ -379,49 +380,49 @@ namespace GridKit using DependencyMap = typename ScalarT::DependencyMap; const auto* f = f_.getData(); - + if (csr_jac_ == nullptr) { IdxT* row_ptrs = new IdxT[static_cast(size_) + 1]; - row_ptrs[0] = 0; - + row_ptrs[0] = 0; + // Count the number of non-zeros IdxT nnz = 0; for (IdxT row = 0; row < size_; ++row) { DependencyMap row_map; - + for (const auto& dep : f[row].getDependencies()) { const auto col = dep.first; - + // Merge-count y and yp dependencies - const IdxT jac_col = static_cast(col / 2); - + const size_t jac_col = static_cast(col / 2); + if (row_map.insert({jac_col, RealT{}}).second) { ++nnz; } } - + row_ptrs[static_cast(row) + 1] = nnz; } - + // Allocate column and value pointers - IdxT* cols = new IdxT[static_cast(nnz)]; + IdxT* cols = new IdxT[static_cast(nnz)]; RealT* vals = new RealT[static_cast(nnz)]; - + // Store column and values IdxT i = 0; for (IdxT row = 0; row < size_; ++row) { DependencyMap row_map; - + for (const auto& dep : f[row].getDependencies()) { const auto col = dep.first; - - const IdxT jac_col = static_cast(col / 2); + + const size_t jac_col = static_cast(col / 2); // Even indices for y and odd indices for yp if (col % 2 == 0) { @@ -432,7 +433,7 @@ namespace GridKit row_map[jac_col] += alpha_ * static_cast(dep.second); } } - + for (const auto& entry : row_map) { cols[i] = static_cast(entry.first); @@ -440,24 +441,24 @@ namespace GridKit ++i; } } - - nnz_ = nnz; + + nnz_ = nnz; csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); } else { RealT* vals = csr_jac_->getValues(); - + IdxT i = 0; for (IdxT row = 0; row < size_; ++row) { DependencyMap row_map; - + for (const auto& dep : f[row].getDependencies()) { const auto col = dep.first; - - const IdxT jac_col = static_cast(col / 2); + + const size_t jac_col = static_cast(col / 2); // Even indices for y and odd indices for yp if (col % 2 == 0) { @@ -468,7 +469,7 @@ namespace GridKit row_map[jac_col] += alpha_ * static_cast(dep.second); } } - + for (const auto& entry : row_map) { vals[i] = static_cast(entry.second); @@ -480,6 +481,36 @@ namespace GridKit return 0; } + /** + * @brief Initialize DependencyTracking variable numbers. + * + * @note Assigns even indices to y and odd indices to yp. + */ + int initializeDependencyTrackingVariableNumbers() + { + static_assert(std::is_same_v, + "initializeDependencyTrackingVariableNumbers() requires ScalarT = DependencyTracking::Variable"); + + auto* y = y_.getData(); + auto* yp = yp_.getData(); + + for (IdxT j = 0; j < size_; ++j) + { + const IdxT var_idx = this->getVariableIndex(j); + if (var_idx != INVALID_INDEX) + { + // Even indices for y and odd indices for yp + y[j].setVariableNumber(static_cast(2 * var_idx)); + yp[j].setVariableNumber(static_cast(2 * var_idx + 1)); + } + } + + y_.setDataUpdated(); + yp_.setDataUpdated(); + + return 0; + } + IdxT size_{0}; IdxT nnz_{0}; /// Global (system-level) variable indices diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index 447fd61cb..e1c7ba7f1 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -13,8 +13,8 @@ namespace GridKit template bool SystemModel::hasJacobian() { - Log::warning() << "GridKit was not built with Enzyme. " - << "DependencyTracking::Variable Jacobians are only available for testing in PhasorDynamics.\n"; + Log::warning() << "DependencyTracking::Variable Jacobians are only available for testing.\n" + << "Falling back to dense Jacobians for PhasorDyanmics simulations.\n"; return false; } @@ -25,12 +25,12 @@ namespace GridKit */ template int SystemModel::evaluateJacobian() - { + { this->constructCsr(); - //Log::misc() << "System DependencyTracking Jacobian\n"; - //csr_jac_->print(Log::misc()); - + // Log::misc() << "System DependencyTracking Jacobian\n"; + // csr_jac_->print(Log::misc()); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp b/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp index 08c434a8c..14300ccee 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp @@ -214,13 +214,12 @@ namespace GridKit } } - //Log::misc() << "System DependencyTracking Jacobian\n"; - //csr_jac_->print(Log::misc()); + // Log::misc() << "System DependencyTracking Jacobian\n"; + // csr_jac_->print(Log::misc()); return 0; } - // Available template instantiations // template class SystemModel; template class SystemModel; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index fda98a3e5..598f1d23d 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -10,6 +10,7 @@ // Include all components #include + namespace GridKit { namespace PhasorDynamics @@ -556,27 +557,15 @@ namespace GridKit status += component->initialize(); } - // For DependencyTracking::Variable, set variable numbers - if constexpr (std::is_same_v) - { - auto* y = y_.getData(); - auto* yp = yp_.getData(); - - for (IdxT j = 0; j < size_; ++j) - { - const IdxT var_idx = this->getVariableIndex(j); - if (var_idx != INVALID_INDEX) - { - // Even indices for y and odd indices for yp - y[j].setVariableNumber(static_cast(2 * var_idx)); - yp[j].setVariableNumber(static_cast(2 * var_idx + 1)); - } - } - } - y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return status; } From 88ec226c1c1bcb2fe08caa5b17635862b73f3f5d Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Tue, 8 Sep 2026 10:55:00 -0400 Subject: [PATCH 09/12] Remove std::vector from DependencyTracking::Variable. --- .../DependencyTracking/Variable.hpp | 5 ----- .../DependencyTracking/VariableImplementation.hpp | 13 ------------- 2 files changed, 18 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp index 6e2f216d0..3e874e986 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp @@ -13,7 +13,6 @@ #include #include #include -#include #include #include @@ -265,10 +264,6 @@ namespace GridKit using DependencyMap = std::map; inline const DependencyMap& getDependencies() const; - // set as the independent state variable and assign ID to it - inline void registerVariable(std::vector& x, - const size_t& offset); - // adds all dependencies of v to *this inline void addDependencies(const Variable& v); diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/VariableImplementation.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/VariableImplementation.hpp index 96b3d12f4..43288eb82 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/VariableImplementation.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/VariableImplementation.hpp @@ -16,19 +16,6 @@ namespace GridKit return dependencies_; } - /** - @brief Registers a variable as an unknown of the system and - adds a pointer to the global @a x vector. - */ - void Variable::registerVariable(std::vector& x, - const size_t& offset) - { - setVariableNumber(offset); // define global variable number - setFixed(false); // not a constant - - x[offset] = this; - } - /** @brief Adds all dependencies of v to *this. */ From fde40bb3f40f19d3094939b5f50e235361361044 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Tue, 8 Sep 2026 12:00:40 -0400 Subject: [PATCH 10/12] Update Variable::der(size_t i) implementation... no insertion if not found --> no need for mutable. --- .../AutomaticDifferentiation/DependencyTracking/Variable.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp index 3e874e986..6e23c8211 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp @@ -208,7 +208,8 @@ namespace GridKit */ double der(size_t i) const { - return dependencies_[i]; + auto it = dependencies_.find(i); + return it != dependencies_.end() ? it->second : 0.0; } /** @@ -300,7 +301,7 @@ namespace GridKit size_t variable_number_; ///< Independent variable ID bool is_fixed_; ///< Constant parameter flag. - mutable DependencyMap dependencies_; + DependencyMap dependencies_; static const size_t INVALID_VAR_NUMBER = INVALID_INDEX; }; From a76935eff35341d968423e2bfaae8140763c8141 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Tue, 8 Sep 2026 18:01:04 -0400 Subject: [PATCH 11/12] WIP: Component-level DependencyTracking Jacobians. --- .../DependencyTracking/Variable.hpp | 4 + .../Branch/BranchDependencyTracking.cpp | 10 +- GridKit/Model/PhasorDynamics/Bus/Bus.hpp | 29 +++++ .../Bus/BusDependencyTracking.cpp | 17 ++- GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp | 6 + .../Bus/BusInfiniteDependencyTracking.cpp | 5 + GridKit/Model/PhasorDynamics/BusBase.hpp | 1 - .../BusFault/BusFaultDependencyTracking.cpp | 22 +++- .../PhasorDynamics/BusFault/BusFaultImpl.hpp | 6 + .../BusToSignalAdapterDependencyTracking.cpp | 2 - GridKit/Model/PhasorDynamics/Component.hpp | 17 ++- .../REECB/ReecbDependencyTracking.cpp | 17 ++- .../Controller/REECB/ReecbImpl.hpp | 17 ++- .../REPCA/RepcaDependencyTracking.cpp | 17 ++- .../Controller/REPCA/RepcaImpl.hpp | 15 ++- .../REGCA/RegcaDependencyTracking.cpp | 16 ++- .../Converter/REGCA/RegcaImpl.hpp | 11 +- .../ESDC1A/Esdc1aDependencyTracking.cpp | 19 ++- .../Exciter/ESDC1A/Esdc1aImpl.hpp | 15 ++- .../IEEET1/Ieeet1DependencyTracking.cpp | 14 ++- .../Exciter/IEEET1/Ieeet1Impl.hpp | 16 ++- .../SEXS-PTI/SexsPtiDependencyTracking.cpp | 18 ++- .../Exciter/SEXS-PTI/SexsPtiImpl.hpp | 14 ++- .../GASTPTI/GastPtiDependencyTracking.cpp | 19 ++- .../Governor/GASTPTI/GastPtiImpl.hpp | 9 +- .../HYGOV/HygovDependencyTracking.cpp | 19 ++- .../Governor/HYGOV/HygovImpl.hpp | 11 +- .../Tgov1/Tgov1DependencyTracking.cpp | 15 ++- .../Governor/Tgov1/Tgov1Impl.hpp | 8 +- .../Load/LoadZ/LoadZDependencyTracking.cpp | 13 +- .../PhasorDynamics/Load/LoadZ/LoadZImpl.hpp | 6 + .../LoadZIP/LoadZIPDependencyTracking.cpp | 13 +- .../Load/LoadZIP/LoadZIPImpl.hpp | 7 ++ .../SignalNodeDependencyTracking.cpp | 4 +- ...ConstantSignalSourceDependencyTracking.cpp | 5 + .../IEEEST/IeeestDependencyTracking.cpp | 15 ++- .../Stabilizer/IEEEST/IeeestImpl.hpp | 6 + .../GENROU/GenrouDependencyTracking.cpp | 15 ++- .../SynchronousMachine/GENROU/GenrouImpl.hpp | 10 +- .../GENSAL/GensalDependencyTracking.cpp | 15 ++- .../SynchronousMachine/GENSAL/GensalImpl.hpp | 10 +- .../GenClassicalDependencyTracking.cpp | 13 +- .../GenClassical/GenClassicalImpl.hpp | 10 +- .../SystemModelDependencyTracking.cpp | 3 +- .../PhasorDynamics/BusFaultTests.hpp | 114 +++--------------- .../PhasorDynamics/ControllerReecbTests.hpp | 46 +++---- .../PhasorDynamics/ControllerRepcaTests.hpp | 38 +++--- .../PhasorDynamics/ConverterRegcaTests.hpp | 42 ++----- .../PhasorDynamics/ExciterEsdc1aTests.hpp | 28 ++--- .../PhasorDynamics/ExciterIeeet1Tests.hpp | 102 +++------------- .../PhasorDynamics/ExciterSexsPtiTests.hpp | 103 +++------------- .../PhasorDynamics/GenClassicalTests.hpp | 92 +++----------- .../UnitTests/PhasorDynamics/GenrouTests.hpp | 101 +++------------- .../UnitTests/PhasorDynamics/GensalTests.hpp | 94 +++------------ .../PhasorDynamics/GovernorGastPtiTests.hpp | 44 +++---- .../PhasorDynamics/GovernorHygovTests.hpp | 23 ++-- .../PhasorDynamics/GovernorTgov1Tests.hpp | 91 ++------------ .../UnitTests/PhasorDynamics/LoadZIPTests.hpp | 84 +++---------- tests/UnitTests/PhasorDynamics/LoadZTests.hpp | 55 ++++----- .../PhasorDynamics/StabilizerIeeestTests.hpp | 76 ++---------- .../UnitTests/PhasorDynamics/SystemTests.hpp | 4 +- 61 files changed, 668 insertions(+), 1003 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp index 6e23c8211..3f9bc2cfe 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp @@ -255,6 +255,10 @@ namespace GridKit /** @brief Turns variable into parameter, or vice versa. + + @todo is_fixed_ is currently not contributing to the semantics of + the derivatives. Leaving as-is for now, as it is not used + for anything other than printed diagnostics. */ void setFixed(bool b = false) { diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp index c06da9396..547e8526c 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp @@ -11,15 +11,19 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * No-op, because branch currently does not own residual equations. * * @return int - error code, 0 = success */ template int Branch::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Branch...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp index 7a3d375da..2c1aa7617 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp @@ -124,6 +124,35 @@ namespace GridKit return 0; } + /** + * @brief Initialize DependencyTracking variable numbers. + * + * @note Assigns even indices to y and odd indices to yp. + * Should be called in intialize(), after variables have been set (and updated as needed). + */ + int initializeDependencyTrackingVariableNumbers() + requires std::is_same_v + { + auto* y = y_.getData(); + auto* yp = yp_.getData(); + + for (IdxT j = 0; j < size_; ++j) + { + const IdxT var_idx = this->getVariableIndex(j); + if (var_idx != INVALID_INDEX) + { + // Even indices for y and odd indices for yp + y[j].setVariableNumber(static_cast(2 * var_idx)); + yp[j].setVariableNumber(static_cast(2 * var_idx + 1)); + } + } + + y_.setDataUpdated(); + yp_.setDataUpdated(); + + return 0; + } + IdxT* J_rows_buffer_{nullptr}; IdxT* J_cols_buffer_{nullptr}; RealT* J_vals_buffer_{nullptr}; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp index 355ad5b05..a24ca20b8 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp @@ -1,3 +1,8 @@ +/** + * @file BusDependencyTracking.cpp + * @author Slaven Peles (peless@ornl.gov) + * + */ #include "BusImpl.hpp" @@ -6,15 +11,21 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). + * Not yet implemented for bus residuals. * * @return int - error code */ template int Bus::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Bus..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Bus...\n"; + Log::misc() << "Jacobian evaluation is not implemented!\n"; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 78104751f..7fbd8b9d2 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -185,6 +185,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteDependencyTracking.cpp index d1c41fe84..61a12689f 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteDependencyTracking.cpp @@ -1,3 +1,8 @@ +/** + * @file BusInfiniteDependencyTracking.cpp + * @author Slaven Peles (peless@ornl.gov) + * + */ #include "BusInfiniteImpl.hpp" diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 815ee83ed..7fed07b6c 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -252,7 +252,6 @@ namespace GridKit */ void allocateVectors(IdxT n) { - y_.resize(n); yp_.resize(n); f_.resize(n); diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp index 097a97d14..b50cc68d9 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp @@ -1,3 +1,9 @@ +/** + * @file BusFaultDependencyTracking.cpp + * @author Slaven Peles (peless@ornl.gov) + * + */ + #include "BusFaultImpl.hpp" namespace GridKit @@ -5,19 +11,27 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented yet + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int BusFault::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; - Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for BusFault...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } - // Additional template instantiations + // Available template instantiations template class BusFault; template class BusFault; } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 39d825a8b..89d121246 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -170,6 +170,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterDependencyTracking.cpp index 68e5b1c95..b8ef837ee 100644 --- a/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterDependencyTracking.cpp @@ -2,8 +2,6 @@ * @file BusToSignalAdapterDependencyTracking.cpp * @author Philip Fackler (facklerpw@ornl.gov) * - * @brief Deinition of BusToSignalAdapter connector interface. - * */ #include "BusToSignalAdapterImpl.hpp" diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 69b94db8d..865864374 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -264,7 +264,7 @@ namespace GridKit /** * @brief CSR construction dispatch depending on ScalarT * - * @note Currently only used for testing purposes. + * @note Currently only used for testing. */ int constructCsr() { @@ -339,7 +339,9 @@ namespace GridKit /** * @brief CSR construction from COO. * - * @note Currently only used for testing purposes. + * @note Currently only used for testing. + * @todo The matrix is only computed on the first call, and the data is stale on subsequent calls. + * @todo Unify with system-level construction that retains map_to_csr_. */ int constructCsrFromCoo() { @@ -369,14 +371,12 @@ namespace GridKit /** * @brief CSR construction from Dependency maps. * - * @note Currently only used for testing purposes. + * @note Currently only used for testing. * See \ref initializeDependencyTrackingVariableNumbers() */ int constructCsrFromDependencies() + requires std::is_same_v { - static_assert(std::is_same_v, - "constructCsrFromDependencies() requires ScalarT = DependencyTracking::Variable"); - using DependencyMap = typename ScalarT::DependencyMap; const auto* f = f_.getData(); @@ -485,12 +485,11 @@ namespace GridKit * @brief Initialize DependencyTracking variable numbers. * * @note Assigns even indices to y and odd indices to yp. + * Should be called in intialize(), after variables have been set (and updated as needed). */ int initializeDependencyTrackingVariableNumbers() + requires std::is_same_v { - static_assert(std::is_same_v, - "initializeDependencyTrackingVariableNumbers() requires ScalarT = DependencyTracking::Variable"); - auto* y = y_.getData(); auto* yp = yp_.getData(); diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp index 1b2f97f32..b7dd60c5f 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp @@ -13,19 +13,28 @@ namespace GridKit namespace Controller { /** - * @brief Report that DependencyTracking exposes structure through the - * residual rather than a separately assembled Jacobian. + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). */ template int Reecb::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Reecb...\n"; - Log::misc() << "Jacobian evaluation is not implemented!\n"; + Log::misc() << "Evaluate DependencyTracking Jacobian for Reecb...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } + // Available template instantiations template class Reecb; template class Reecb; + } // namespace Controller } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp index 705cc5d78..2ed478014 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -316,6 +316,13 @@ namespace GridKit } commitInitialPoint(point); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } @@ -706,11 +713,11 @@ namespace GridKit Vref0_ = point.vref; } - pe_set_ = static_cast(point.signals[PE]); - qgen_set_ = static_cast(point.signals[QGEN]); - qext_set_ = static_cast(point.signals[QEXT]); - pfaref_set_ = static_cast(point.signals[PFAREF]); - pref_set_ = static_cast(point.signals[PREF]); + pe_set_ = static_cast(point.signals[PE]); + qgen_set_ = static_cast(point.signals[QGEN]); + qext_set_ = static_cast(point.signals[QEXT]); + pfaref_set_ = static_cast(point.signals[PFAREF]); + pref_set_ = static_cast(point.signals[PREF]); if (ports_.in.template port()) { diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp index e3d8095a7..185097dda 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp @@ -13,19 +13,28 @@ namespace GridKit namespace Controller { /** - * @brief Report that DependencyTracking exposes structure through the - * residual rather than a separately assembled Jacobian. + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). */ template int Repca::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Repca...\n"; - Log::misc() << "Jacobian evaluation is not implemented!\n"; + Log::misc() << "Evaluate DependencyTracking Jacobian for Repca...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } + // Available template instantiations template class Repca; template class Repca; + } // namespace Controller } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp index dbfd8c6bc..df9da208a 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp @@ -495,10 +495,10 @@ namespace GridKit Pmin_ = pmin; Pmax_ = pmax; - freqref_set_ = freqref0; - vref_set_ = vref0; - qref_set_ = qref0_system; - pref_set_ = pref0_system; + freqref_set_ = static_cast(freqref0); + vref_set_ = static_cast(vref0); + qref_set_ = static_cast(qref0_system); + pref_set_ = static_cast(pref0_system); if (auto vref_port = ports_.in.template port()) { @@ -537,6 +537,13 @@ namespace GridKit y_.setDataUpdated(); yp_.setToConst(static_cast(ZERO)); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp index 1cb6e1b82..6945845d9 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp @@ -13,20 +13,28 @@ namespace GridKit namespace Converter { /** - * @brief Report that dependency tracking does not assemble a separate Jacobian. + * @brief Evaluate DependencyTracking::Variable Jacobian. * - * Dependency tracking recovers the sparsity pattern from the residual. + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). */ template int Regca::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Regca...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } + // Available template instantiations template class Regca; template class Regca; + } // namespace Converter } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp index 591fd614c..233a5be86 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -495,8 +495,8 @@ namespace GridKit y[PBR] = vr * y[IR] + vi * y[II]; y[QBR] = vi * y[IR] - vr * y[II]; - ipcmd_set_ = this->toSystemBase(ipcmd0); - iqcmd_set_ = this->toSystemBase(iqcmd0); + ipcmd_set_ = static_cast(this->toSystemBase(ipcmd0)); + iqcmd_set_ = static_cast(this->toSystemBase(iqcmd0)); // Publish the resolved system-base commands for downstream controller // initialization. Unattached ports retain these values as constant @@ -512,6 +512,13 @@ namespace GridKit y_.setDataUpdated(); yp_.setToConst(static_cast(ZERO)); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp index dced64f5b..387677c72 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp @@ -13,21 +13,28 @@ namespace GridKit namespace Exciter { /** - * @brief Report that dependency tracking does not assemble a separate Jacobian. + * @brief Evaluate DependencyTracking::Variable Jacobian. * - * Dependency tracking exposes the Jacobian structure through the - * residual rather than a separately assembled matrix. + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). */ template int Esdc1a::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Esdc1a...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } - + + // Available template instantiations template class Esdc1a; template class Esdc1a; + } // namespace Exciter } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index f2ed50915..425a29f91 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -378,10 +378,10 @@ namespace GridKit y[VFE] = vfe0; y[EFD] = efd0; - omega_set_ = omega0; - vref_set_ = vref0; - vs_set_ = vs0; - vuel_set_ = vuel0; + omega_set_ = static_cast(omega0); + vref_set_ = static_cast(vref0); + vs_set_ = static_cast(vs0); + vuel_set_ = static_cast(vuel0); if (auto vref_port = ports_.in.template port()) { @@ -390,6 +390,13 @@ namespace GridKit y_.setDataUpdated(); yp_.setToConst(static_cast(ZERO)); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp index 6638f952b..43df43a0d 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp @@ -16,15 +16,23 @@ namespace GridKit namespace Exciter { /** - * @brief Jacobian evaluation not implemented yet + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int Ieeet1::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; - Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Ieeet1...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index a67560295..6be4cf1eb 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -300,11 +300,11 @@ namespace GridKit yp[i] = 0.0; } - omega_set_ = omega; - vref_set_ = vref; - vs_set_ = vs; - vuel_set_ = vuel; - voel_set_ = voel; + omega_set_ = static_cast(omega); + vref_set_ = static_cast(vref); + vs_set_ = static_cast(vs); + vuel_set_ = static_cast(vuel); + voel_set_ = static_cast(voel); if (auto vref_port = ports_.in.template port()) { @@ -314,6 +314,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp index 705734f27..e79ad9715 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp @@ -12,14 +12,28 @@ namespace GridKit { namespace Exciter { + /** + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). + * + * @return int - error code, 0 = success + */ template int SexsPti::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for SexsPti..." << std::endl; - Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for SexsPti...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } + // Available template instantiations template class SexsPti; template class SexsPti; diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index 16c1e0cce..5c9b035d2 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -206,10 +206,10 @@ namespace GridKit yp[static_cast(i)] = 0.0; } - vref_set_ = vref; - vs_set_ = vs; - vuel_set_ = vuel; - voel_set_ = voel; + vref_set_ = static_cast(vref); + vs_set_ = static_cast(vs); + vuel_set_ = static_cast(vuel); + voel_set_ = static_cast(voel); if (auto vref_port = ports_.in.template port()) { @@ -219,6 +219,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp index 3e681c007..6dde71a2f 100644 --- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp @@ -13,19 +13,30 @@ namespace GridKit namespace Governor { /** - * @brief Report that DependencyTracking exposes structure through the - * residual rather than a separately assembled Jacobian. + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). + * + * @return int - error code, 0 = success */ template int GastPti::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for GastPti..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for GastPti...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } + // Available template instantiations template class GastPti; template class GastPti; + } // namespace Governor } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp index a94410c75..e35cc8a2a 100644 --- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp @@ -387,7 +387,7 @@ namespace GridKit y[VTEMP] = static_cast(vtemp0); y[VLV] = static_cast(vlv0); - pref_set_ = static_cast(pref0); + pref_set_ = static_cast(pref0); if (auto pref_port = ports_.in.template port()) { pref_port.writeValue(pref_set_); @@ -395,6 +395,13 @@ namespace GridKit y_.setDataUpdated(); yp_.setToConst(static_cast(ZERO)); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp index 760bac957..406faf306 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp @@ -12,16 +12,31 @@ namespace GridKit { namespace Governor { + /** + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). + * + * @return int - error code, 0 = success + */ template int Hygov::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Hygov...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } + // Available template instantiations template class Hygov; template class Hygov; + } // namespace Governor } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp index 84a82d356..71ada99e0 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -399,8 +399,8 @@ namespace GridKit Gmin_response_ = Gmin_response; Gmax_response_ = Gmax_response; Hdam_eff_ = Hdam0; - pref_set_ = pref0; - paux_set_ = paux0_system; + pref_set_ = static_cast(pref0); + paux_set_ = static_cast(paux0_system); if (auto pref_port = ports_.in.template port()) { @@ -419,6 +419,13 @@ namespace GridKit y_.setDataUpdated(); yp_.setToConst(static_cast(ZERO)); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp index 278641e18..7f6178231 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp @@ -12,21 +12,30 @@ namespace GridKit namespace Governor { /** - * @brief Jacobian evaluation not implemented yet + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int Tgov1::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; - Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Tgov1...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } // Available template instantiations template class Tgov1; template class Tgov1; + } // namespace Governor } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 298ce785f..01843aa0f 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -298,7 +298,7 @@ namespace GridKit y[PTX] = pturb0; y[PV] = pv0; - pref_set_ = pref0; + pref_set_ = static_cast(pref0); if (auto pref_port = ports_.in.template port()) { pref_port.writeValue(pref_set_); @@ -307,6 +307,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setToConst(static_cast(ZERO)); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp index 2682dba88..b02b683a2 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp @@ -6,15 +6,22 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int LoadZ::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for LoadZ..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for LoadZ...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp index cea110760..b3e24fc8a 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp @@ -132,6 +132,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp index 514ead2e9..2c84ee6b1 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp @@ -6,15 +6,22 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int LoadZIP::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for LoadZIP..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for LoadZIP...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp index 19bd6ed6c..929cfb746 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp @@ -154,6 +154,13 @@ namespace GridKit yp[1] = 0.0; y_.setDataUpdated(); yp_.setDataUpdated(); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp index df00cf7a9..33dc3810d 100644 --- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp @@ -1,6 +1,7 @@ /** - * @file SignalNode model implementation. + * @file SignalNodeDependencyTracking.cpp */ + #include #include "SignalNodeImpl.hpp" @@ -9,6 +10,7 @@ namespace GridKit { namespace PhasorDynamics { + // Available template instantiations template class SignalNode; template class SignalNode; diff --git a/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceDependencyTracking.cpp index d3f1df113..dafdd1a71 100644 --- a/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceDependencyTracking.cpp @@ -1,3 +1,8 @@ +/** + * @file ConstantSignalSourceDependencyTracking.cpp + * @author Philip Fackler (facklerpw@ornl.gov) + * + */ #include "ConstantSignalSourceImpl.hpp" diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp index 083b5c7e2..eed3e9bd9 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp @@ -13,21 +13,30 @@ namespace GridKit namespace Stabilizer { /** - * @brief Jacobian evaluation not implemented yet + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int Ieeest::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Ieeest..." << std::endl; - Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Ieeest...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } // Available template instantiations template class Ieeest; template class Ieeest; + } // namespace Stabilizer } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index 1d7d6d9b0..2bc836363 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -249,6 +249,12 @@ namespace GridKit y[10] = bypass_T6_block_ * Ks_ * u; y[11] = Math::clamp(y[10], Lsmin_, Lsmax_); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + y_.setDataUpdated(); yp_.setDataUpdated(); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp index 7d5001b5e..fdf9b1805 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp @@ -5,20 +5,29 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented yet + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int Genrou::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; - Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Genrou...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } // Available template instantiations template class Genrou; template class Genrou; + } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp index 2051a9214..c052fec4f 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp @@ -437,13 +437,13 @@ namespace GridKit // Convert Te to system base for governor PM signal. ScalarT Te = y[12]; - pmech_set_ = this->toSystemBase(Te); + pmech_set_ = static_cast(this->toSystemBase(Te)); if (auto pmech_port = ports_.in.template port()) { pmech_port.writeValue(pmech_set_); } - efd_set_ = Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + psidpp * ksat; + efd_set_ = static_cast(Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + psidpp * ksat); if (auto efd_port = ports_.in.template port()) { efd_port.writeValue(efd_set_); @@ -457,6 +457,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp index 741453ce6..f366f8029 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp @@ -5,20 +5,29 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented yet + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int Gensal::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for Gensal..." << std::endl; - Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for Gensal...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); + return 0; } // Available template instantiations template class Gensal; template class Gensal; + } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp index e656a2f6b..8435bdb71 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp @@ -297,13 +297,13 @@ namespace GridKit y[13] = ii; // Convert Te to system base for governor PM signal. - pmech_set_ = this->toSystemBase(Te); + pmech_set_ = static_cast(this->toSystemBase(Te)); if (auto pmech_port = ports_.in.template port()) { pmech_port.writeValue(pmech_set_); } - efd_set_ = Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + Eqp * ksat; + efd_set_ = static_cast(Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + Eqp * ksat); if (auto efd_port = ports_.in.template port()) { efd_port.writeValue(efd_set_); @@ -317,6 +317,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp index 7a173534a..e4eebae9e 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp @@ -6,15 +6,22 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation not implemented + * @brief Evaluate DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + * + * DependencyTracking::Variable stores the Jacobian as dependency maps, + * updated during calls to evaluateResidual(). * * @return int - error code, 0 = success */ template int GenClassical::evaluateJacobian() { - Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; - Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate DependencyTracking Jacobian for GenClassical...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + this->constructCsr(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index cd61c29ea..0c3b12cf1 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -220,13 +220,13 @@ namespace GridKit y[4] = ii; // Convert Te to system base for governor PM signal. - pmech_set_ = this->toSystemBase(Te); + pmech_set_ = static_cast(this->toSystemBase(Te)); if (auto pmech_port = ports_.in.template port()) { pmech_port.writeValue(pmech_set_); } - efd_set_ = efd; + efd_set_ = static_cast(efd); if (auto efd_port = ports_.in.template port()) { efd_port.writeValue(efd_set_); @@ -240,6 +240,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index e1c7ba7f1..f6105911e 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -7,8 +7,6 @@ namespace GridKit /** * @brief By default, Jacobians are not available * - * DependencyTracking::Variable stores the Jacobian as dependency maps, - * updated during calls to evaluateResidual(). */ template bool SystemModel::hasJacobian() @@ -22,6 +20,7 @@ namespace GridKit /** * @brief Evaluate system DependencyTracking::Variable Jacobian. * + * @note Currently only used for testing. */ template int SystemModel::evaluateJacobian() diff --git a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp index 03528173e..9fa45e096 100644 --- a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp @@ -133,108 +133,26 @@ namespace GridKit bus.allocate(); fault.allocate(); - // Get d/dy - bus.initialize(); - fault.initialize(); - - auto* fault_y = fault.y().getData(); - for (size_t i = 0; i < fault.size(); ++i) - { - fault_y[i].setVariableNumber(i); ///< fault independent variables - } - fault.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + fault.size()); // Bus independent variables + bus.setVariableIndex(i, i + fault.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + fault.size()); // Reset bus residual indices } - bus.y().setDataUpdated(); - bus.evaluateResidual(); - fault.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_y_view = fault.getResidual(); - std::vector residual_y( - residual_y_view.getData(), - residual_y_view.getData() + residual_y_view.getSize()); - auto& bus_residual_y_view = bus.getResidual(); - std::vector bus_residual_y( - bus_residual_y_view.getData(), - bus_residual_y_view.getData() + bus_residual_y_view.getSize()); - - // Get d/dy' bus.initialize(); fault.initialize(); - auto* fault_yp = fault.yp().getData(); - for (size_t i = 0; i < fault.size(); ++i) - { - fault_yp[i].setVariableNumber(i); ///< fault independent variables - } - fault.yp().setDataUpdated(); + fault.updateTime(0.0, 1.0); bus.evaluateResidual(); - fault.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_yp_view = fault.getResidual(); - std::vector residual_yp( - residual_yp_view.getData(), - residual_yp_view.getData() + residual_yp_view.getSize()); - - // Print the dependencies - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - // Extract the dependencies and add d/dy' to d/dy - std::vector dependencies( - residual_y.size() + bus_residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); - DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - - for (const auto& pair_y : dependency_y) - { - auto index_y = pair_y.first; - auto value_y = pair_y.second; - auto it_yp = dependency_yp.find(index_y); - if (it_yp != dependency_yp.end()) - { - auto value_yp = it_yp->second; - dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); - } - else - { - dependencies[i].insert(std::make_pair(index_y, value_y)); - } - } - - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) - { - auto index_yp = pair_yp.first; - auto value_yp = pair_yp.second; - auto it_y = dependency_y.find(index_yp); - if (it_y == dependency_y.end()) - { - dependencies[i].insert(std::make_pair(index_yp, value_yp)); - } - } - } + fault.evaluateResidual(); - for (size_t i = 0; i < bus_residual_y.size(); ++i) - { - dependencies[residual_y.size() + i] = bus_residual_y[i].getDependencies(); - } + fault.evaluateJacobian(); + auto* model_jacobian = fault.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: BusFault DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian( @@ -249,25 +167,25 @@ namespace GridKit bus.allocate(); fault.allocate(); - bus.initialize(); - fault.initialize(); - - fault.updateTime(0.0, 1.0); - for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + fault.size()); // Reset bus variable indices bus.setResidualIndex(i, i + fault.size()); // Reset bus residual indices } + bus.initialize(); + fault.initialize(); + + fault.updateTime(0.0, 1.0); + bus.evaluateResidual(); fault.evaluateResidual(); bus.evaluateJacobian(); fault.evaluateJacobian(); fault.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = fault.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: BusFault Jacobian\n"; + auto* model_jacobian = fault.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: BusFault Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp index e25665382..65568cecb 100644 --- a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp @@ -355,21 +355,6 @@ namespace GridKit success *= scalarPreserved(latched.ipcmd(), kInitialIpcmd, "unassigned ipcmd"); success *= allResidualsWithinInitTolerance(latched.reecb); - // Initialization preserves AD metadata on the owned command inputs. - Fixture tracked_commands( - makeData(), 1.0, 0.0, kSystemBaseVa, false); - success *= tracked_commands.prepare(kInitialIqcmd, kInitialIpcmd); - auto* tracked_state = tracked_commands.reecb.y().getData(); - const auto iqcmd_index = index(Vars::IQCMD); - const auto ipcmd_index = index(Vars::IPCMD); - tracked_state[iqcmd_index].setVariableNumber(iqcmd_index); - tracked_state[ipcmd_index].setVariableNumber(ipcmd_index); - const auto iqcmd_dependencies = tracked_state[iqcmd_index].getDependencies(); - const auto ipcmd_dependencies = tracked_state[ipcmd_index].getDependencies(); - success *= (tracked_commands.reecb.initialize() == 0); - success *= isEqual(tracked_state[iqcmd_index].getDependencies(), iqcmd_dependencies); - success *= isEqual(tracked_state[ipcmd_index].getDependencies(), ipcmd_dependencies); - return success.report(__func__); } @@ -2753,7 +2738,8 @@ namespace GridKit return success; } - void numberVariables(Fixture& fixture, RealT alpha) const + /// @todo Remove and setup the test to not rely on explicit variable numbering + void numberVariables(Fixture& fixture) const { auto* y = fixture.reecb.y().getData(); auto* yp = fixture.reecb.yp().getData(); @@ -2761,17 +2747,16 @@ namespace GridKit for (size_t row = 0; row < Utilities::enum_size(); ++row) { - y[row].setVariableNumber(row); - yp[row].setVariableNumber(row); - yp[row].scaleDependencies(alpha); + y[row].setVariableNumber(2 * row); + yp[row].setVariableNumber(2 * row + 1); } for (size_t row = 0; row < static_cast(fixture.bus.size()); ++row) { - bus_y[row].setVariableNumber(kBusVrColumn + row); + bus_y[row].setVariableNumber(2 * (kBusVrColumn + row)); } for (auto variable : Utilities::enum_values()) { - fixture.input(variable).setVariableNumber(fixture.inputIndex(variable)); + fixture.input(variable).setVariableNumber(2 * fixture.inputIndex(variable)); } fixture.reecb.y().setDataUpdated(); @@ -2793,16 +2778,12 @@ namespace GridKit fixture.attachAllInputs(); success *= fixture.prepare(0.0, 0.2); setJacobianState(fixture, capacity, epiv, ipcmd); - numberVariables(fixture, alpha); - success *= (fixture.evaluate() == 0); + numberVariables(fixture); + fixture.reecb.updateTime(0.0, alpha); + success *= (fixture.reecb.evaluateResidual() == 0); + success *= (fixture.reecb.evaluateJacobian() == 0); - std::vector rows(Utilities::enum_size()); - const auto* f = fixture.reecb.getResidual().getData(); - for (size_t row = 0; row < rows.size(); ++row) - { - rows[row] = f[row].getDependencies(); - } - return rows; + return GridKit::Testing::MapFromCsr(fixture.reecb.getCsrJacobian()); } #ifdef GRIDKIT_ENABLE_ENZYME @@ -2825,10 +2806,11 @@ namespace GridKit setJacobianState(fixture, capacity, epiv, ipcmd); fixture.reecb.updateTime(0.0, alpha); - success *= (fixture.evaluate() == 0); + success *= (fixture.reecb.evaluateResidual() == 0); success *= (fixture.reecb.evaluateJacobian() == 0); success *= (fixture.reecb.constructCsr() == 0); - return MapFromCsr(fixture.reecb.getCsrJacobian()); + + return GridKit::Testing::MapFromCsr(fixture.reecb.getCsrJacobian()); } bool jacobiansMatch( diff --git a/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp index bb6688cde..7643d30a6 100644 --- a/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp @@ -1014,15 +1014,17 @@ namespace GridKit setState(blocked.repca, {{Vars::QPI, 1.7}, {Vars::ERQLIM, 0.4}, {Vars::SFRZ, 1.0}}); setDerivative(blocked.repca, {{Vars::XQPI, 0.0}}); - numberVariables(blocked, 1.0); + numberVariables(blocked); + blocked.repca.updateTime(0.0, 1.0); success *= (blocked.repca.evaluateResidual() == 0); const DependencyTracking::Variable::DependencyMap expected{ - {index(Vars::XQPI), -1.0}, - {index(Vars::SFRZ), 0.0}, - {index(Vars::ERQLIM), 0.0}, - {index(Vars::QPI), 0.0}, + {2 * index(Vars::XQPI) + 1, -1.0}, // @todo Remove these + {2 * index(Vars::SFRZ), 0.0}, // @todo Remove these + {2 * index(Vars::ERQLIM), 0.0}, // @todo Remove these + {2 * index(Vars::QPI), 0.0}, // @todo Remove these }; + success *= jacobianRowMatches( blocked.repca.getResidual().getData()[index(Vars::XQPI)].getDependencies(), expected, @@ -2282,8 +2284,8 @@ namespace GridKit return success; } - void numberVariables(Fixture& fixture, - RealT alpha) const + /// @todo Remove and setup the test to not rely on explicit variable numbering + void numberVariables(Fixture& fixture) const { auto* y = fixture.repca.y().getData(); auto* yp = fixture.repca.yp().getData(); @@ -2291,17 +2293,16 @@ namespace GridKit for (size_t row = 0; row < Utilities::enum_size(); ++row) { - y[row].setVariableNumber(row); - yp[row].setVariableNumber(row); - yp[row].scaleDependencies(alpha); + y[row].setVariableNumber(2 * row); + yp[row].setVariableNumber(2 * row + 1); } for (size_t row = 0; row < static_cast(fixture.bus.size()); ++row) { - bus_y[row].setVariableNumber(kBusVrColumn + row); + bus_y[row].setVariableNumber(2 * (kBusVrColumn + row)); } for (auto variable : Utilities::enum_values()) { - fixture.input(variable).setVariableNumber(fixture.inputIndex(variable)); + fixture.input(variable).setVariableNumber(2 * fixture.inputIndex(variable)); } fixture.repca.y().setDataUpdated(); @@ -2321,16 +2322,12 @@ namespace GridKit setAnswerKeyInputs(fixture); success *= fixture.prepare(0.0, 0.0); setAnswerKeyState(fixture.repca); - numberVariables(fixture, alpha); + numberVariables(fixture); + fixture.repca.updateTime(0.0, alpha); success *= (fixture.repca.evaluateResidual() == 0); + success *= (fixture.repca.evaluateJacobian() == 0); - std::vector rows(Utilities::enum_size()); - const auto* f = fixture.repca.getResidual().getData(); - for (size_t row = 0; row < Utilities::enum_size(); ++row) - { - rows[row] = f[row].getDependencies(); - } - return rows; + return MapFromCsr(fixture.repca.getCsrJacobian()); } #ifdef GRIDKIT_ENABLE_ENZYME @@ -2354,6 +2351,7 @@ namespace GridKit success *= (fixture.repca.evaluateResidual() == 0); success *= (fixture.repca.evaluateJacobian() == 0); success *= (fixture.repca.constructCsr() == 0); + return MapFromCsr(fixture.repca.getCsrJacobian()); } #endif diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index 47504fdca..cb7ba7f87 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -660,8 +660,8 @@ namespace GridKit fixture.regca.getResidual().getData()[index(Vars::IQEXTRA)].getDependencies(); const DepVar::DependencyMap expected{{ - {index(Vars::VT), 0.5 * kHvrcmGain}, - {index(Vars::IQEXTRA), -1.0}, + {2 * index(Vars::VT), 0.5 * kHvrcmGain}, // @todo remove these + {2 * index(Vars::IQEXTRA), -1.0}, // @todo remove these }}; success *= isEqual(dependencies, expected, kTol); } @@ -694,12 +694,6 @@ namespace GridKit dependencyTrackingJacobian(data, current, success); const auto enzyme_jacobian = enzymeJacobian(data, current, success); - const auto ip_row = index(Vars::IP); - const auto il_col = index(Vars::IL); - success *= dependency_tracking_jacobian[ip_row].contains(il_col); - success *= enzyme_jacobian[ip_row].contains(il_col); - - success *= (dependency_tracking_jacobian.size() == enzyme_jacobian.size()); const auto nrows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); @@ -930,11 +924,7 @@ namespace GridKit regca.yp().setDataUpdated(); } - /// Numbers regca y and yp together, the bus after the regca block, - /// and the command signals at their port indices, matching the - /// Jacobian layout. Write state values first; numbering resets each - /// dependency map. Numbering an unattached command is a no-op for the - /// model. + /// @todo Remove and setup the test to not rely on explicit variable numbering void numberVariables(Fixture& fixture) { auto* y = fixture.regca.y().getData(); @@ -944,15 +934,15 @@ namespace GridKit const auto regca_size = static_cast(fixture.regca.size()); for (size_t i = 0; i < regca_size; ++i) { - y[i].setVariableNumber(i); - yp[i].setVariableNumber(i); + y[i].setVariableNumber(2 * i); + yp[i].setVariableNumber(2 * i + 1); } for (size_t i = 0; i < static_cast(fixture.bus.size()); ++i) { - bus_y[i].setVariableNumber(i + regca_size); + bus_y[i].setVariableNumber(2 * (i + regca_size)); } - fixture.ipcmd.setVariableNumber(fixture.ipcmd_index); - fixture.iqcmd.setVariableNumber(fixture.iqcmd_index); + fixture.ipcmd.setVariableNumber(2 * fixture.ipcmd_index); + fixture.iqcmd.setVariableNumber(2 * fixture.iqcmd_index); fixture.regca.y().setDataUpdated(); fixture.regca.yp().setDataUpdated(); @@ -1037,22 +1027,12 @@ namespace GridKit fixture.iqcmd = kStateIqcmd; setJacobianState(fixture.regca, current); numberVariables(fixture); + fixture.regca.updateTime(0.0, 1.0); success *= (fixture.evaluate() == 0); + success *= (fixture.regca.evaluateJacobian() == 0); - const auto regca_size = static_cast(fixture.regca.size()); - const auto* f = fixture.regca.getResidual().getData(); - - std::vector dependencies( - regca_size + static_cast(fixture.bus.size())); - for (size_t i = 0; i < regca_size; ++i) - { - dependencies[i] = f[i].getDependencies(); - } - dependencies[regca_size] = fixture.bus.Ir().getDependencies(); - dependencies[regca_size + 1] = fixture.bus.Ii().getDependencies(); - - return dependencies; + return MapFromCsr(fixture.regca.getCsrJacobian()); } std::vector enzymeJacobian( diff --git a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp index ba1112e52..3efa7d742 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp @@ -804,9 +804,10 @@ namespace GridKit const auto& dependencies = transition.esdc1a.getResidual().getData()[static_cast(Internal::EFDP)].getDependencies(); const DepVar::DependencyMap expected{{ - {static_cast(Internal::EFDP), -13.0}, - {static_cast(Internal::VR), 1.0}, - {static_cast(Internal::VFE), -1.0}, + {2 * static_cast(Internal::EFDP), -12.0}, // @todo Remove these + {2 * static_cast(Internal::EFDP) + 1, -1.0}, // @todo Remove these + {2 * static_cast(Internal::VR), 1.0}, // @todo Remove these + {2 * static_cast(Internal::VFE), -1.0}, // @todo Remove these }}; success *= isEqual(dependencies, expected, kTol); } @@ -1438,6 +1439,7 @@ namespace GridKit return false; } + /// @todo Remove and setup the test to not rely on explicit variable numbering void numberVariables(Fixture& fixture) const { auto* y = fixture.esdc1a.y().getData(); @@ -1447,16 +1449,16 @@ namespace GridKit const auto model_size = static_cast(fixture.esdc1a.size()); for (size_t i = 0; i < model_size; ++i) { - y[i].setVariableNumber(i); - yp[i].setVariableNumber(i); + y[i].setVariableNumber(2 * i); + yp[i].setVariableNumber(2 * i + 1); } for (size_t i = 0; i < static_cast(fixture.bus.size()); ++i) { - bus_y[i].setVariableNumber(model_size + i); + bus_y[i].setVariableNumber(2 * (model_size + i)); } for (auto port : Utilities::enum_values()) { - fixture.input(port).setVariableNumber(fixture.inputIndex(port)); + fixture.input(port).setVariableNumber(2 * fixture.inputIndex(port)); } fixture.esdc1a.y().setDataUpdated(); @@ -1478,16 +1480,11 @@ namespace GridKit fixture.input(External::vuel) = kJacobianVuel; setAnswerKeyState(fixture.esdc1a); numberVariables(fixture); + fixture.esdc1a.updateTime(0.0, 1.0); success *= (fixture.evaluate() == 0); + success *= (fixture.esdc1a.evaluateJacobian() == 0); - const auto model_size = static_cast(fixture.esdc1a.size()); - std::vector rows(model_size); - const auto* f = fixture.esdc1a.getResidual().getData(); - for (size_t i = 0; i < model_size; ++i) - { - rows[i] = f[i].getDependencies(); - } - return rows; + return MapFromCsr(fixture.esdc1a.getCsrJacobian()); } std::vector enzymeJacobian( @@ -1510,6 +1507,7 @@ namespace GridKit success *= (fixture.evaluate() == 0); success *= (fixture.esdc1a.evaluateJacobian() == 0); success *= (fixture.esdc1a.constructCsr() == 0); + return MapFromCsr(fixture.esdc1a.getCsrJacobian()); } #endif diff --git a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp index c674cb0fe..b8d15e58a 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp @@ -303,95 +303,27 @@ namespace GridKit bus.allocate(); exciter.allocate(); - // Get d/dy - bus.initialize(); - exciter.initialize(); - - auto* exciter_y = exciter.y().getData(); - for (size_t i = 0; i < exciter.size(); ++i) - { - exciter_y[i].setVariableNumber(i); ///< Exciter independent variables - } - exciter.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + exciter.size()); // Bus independent variables + bus.setVariableIndex(i, i + exciter.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + exciter.size()); // Reset bus residual indices } - bus.y().setDataUpdated(); - - bus.evaluateResidual(); - exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_y_view = exciter.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); - // Get d/dy' bus.initialize(); exciter.initialize(); - auto* exciter_yp = exciter.yp().getData(); - for (size_t i = 0; i < exciter.size(); ++i) - { - exciter_yp[i].setVariableNumber(i); ///< Exciter independent variables - } - exciter.yp().setDataUpdated(); + exciter.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term bus.evaluateResidual(); - exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_yp_view = exciter.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); - - // Print the dependencies - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - // Extract the dependencies and add d/dy' to d/dy - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); - DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - - for (const auto& pair_y : dependency_y) - { - auto index_y = pair_y.first; - auto value_y = pair_y.second; - auto it_yp = dependency_yp.find(index_y); - if (it_yp != dependency_yp.end()) - { - auto value_yp = it_yp->second; - dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); - } - else - { - dependencies[i].insert(std::make_pair(index_y, value_y)); - } - } + exciter.evaluateResidual(); - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) - { - auto index_yp = pair_yp.first; - auto value_yp = pair_yp.second; - auto it_y = dependency_y.find(index_yp); - if (it_y == dependency_y.end()) - { - dependencies[i].insert(std::make_pair(index_yp, value_yp)); - } - } - } + bus.evaluateJacobian(); + exciter.evaluateJacobian(); + auto* model_jacobian = exciter.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Ieeet1 DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; - } + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian() { @@ -403,25 +335,25 @@ namespace GridKit bus.allocate(); exciter.allocate(); - bus.initialize(); - exciter.initialize(); - - exciter.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term - for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + exciter.size()); // Reset bus variable indices bus.setResidualIndex(i, i + exciter.size()); // Reset bus residual indices } + bus.initialize(); + exciter.initialize(); + + exciter.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term + bus.evaluateResidual(); exciter.evaluateResidual(); bus.evaluateJacobian(); exciter.evaluateJacobian(); exciter.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = exciter.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Ieeet1 Jacobian\n"; + auto* model_jacobian = exciter.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Ieeet1 Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp index 0de084ff5..0388d5742 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp @@ -32,7 +32,7 @@ namespace GridKit // Init and saturated-limiter checks are exact algebraic identities // (no iteration, sigmoid underflowed to 0/1 at the depths tested here). - static constexpr ScalarT kTol = static_cast(1.0e-14); + static constexpr ScalarT kTol = 100 * std::numeric_limits::epsilon(); TestOutcome constructor() { @@ -430,94 +430,27 @@ namespace GridKit bus.allocate(); exciter.allocate(); - // Get d/dy - bus.initialize(); - exciter.initialize(); - - auto* exciter_y = exciter.y().getData(); - for (size_t i = 0; i < exciter.size(); ++i) - { - exciter_y[i].setVariableNumber(i); ///< Exciter independent variables - } - exciter.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + exciter.size()); // Bus independent variables + bus.setVariableIndex(i, i + exciter.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + exciter.size()); // Reset bus residual indices } - bus.y().setDataUpdated(); - - bus.evaluateResidual(); - exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_y_view = exciter.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); - // Get d/dy' bus.initialize(); exciter.initialize(); - auto* exciter_yp = exciter.yp().getData(); - for (size_t i = 0; i < exciter.size(); ++i) - { - exciter_yp[i].setVariableNumber(i); ///< Exciter independent variables - } - exciter.yp().setDataUpdated(); + exciter.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term bus.evaluateResidual(); - exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_yp_view = exciter.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); - - // Print the dependencies - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - // Extract the dependencies and add d/dy' to d/dy - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); - DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - - for (const auto& pair_y : dependency_y) - { - auto index_y = pair_y.first; - auto value_y = pair_y.second; - auto it_yp = dependency_yp.find(index_y); - if (it_yp != dependency_yp.end()) - { - auto value_yp = it_yp->second; - dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); - } - else - { - dependencies[i].insert(std::make_pair(index_y, value_y)); - } - } + exciter.evaluateResidual(); - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) - { - auto index_yp = pair_yp.first; - auto value_yp = pair_yp.second; - auto it_y = dependency_y.find(index_yp); - if (it_y == dependency_y.end()) - { - dependencies[i].insert(std::make_pair(index_yp, value_yp)); - } - } - } + bus.evaluateJacobian(); + exciter.evaluateJacobian(); + auto* model_jacobian = exciter.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: SexsPti DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian() @@ -530,25 +463,25 @@ namespace GridKit bus.allocate(); exciter.allocate(); - bus.initialize(); - exciter.initialize(); - - exciter.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term - for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + exciter.size()); // Reset bus variable indices bus.setResidualIndex(i, i + exciter.size()); // Reset bus residual indices } + bus.initialize(); + exciter.initialize(); + + exciter.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term + bus.evaluateResidual(); exciter.evaluateResidual(); bus.evaluateJacobian(); exciter.evaluateJacobian(); exciter.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = exciter.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: SexsPti Jacobian\n"; + auto* model_jacobian = exciter.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: SexsPti Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index e81833e25..aa2ebff67 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -458,88 +458,26 @@ namespace GridKit bus.allocate(); gen.allocate(); - bus.initialize(); - gen.initialize(); - - auto* gen_y = gen.y().getData(); - for (size_t i = 0; i < gen.size(); ++i) - { - gen_y[i].setVariableNumber(i); - } - gen.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + gen.size()); + bus.setVariableIndex(i, i + gen.size()); + bus.setResidualIndex(i, i + gen.size()); } - bus.y().setDataUpdated(); - - bus.evaluateResidual(); - gen.evaluateResidual(); - auto& residual_y_view = gen.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); bus.initialize(); gen.initialize(); - auto* gen_yp = gen.yp().getData(); - for (size_t i = 0; i < gen.size(); ++i) - { - gen_yp[i].setVariableNumber(i); - } - gen.yp().setDataUpdated(); + gen.updateTime(0.0, 1.0); bus.evaluateResidual(); gen.evaluateResidual(); - auto& residual_yp_view = gen.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); - - // Print the dependencies - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); - DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - for (const auto& pair_y : dependency_y) - { - auto index_y = pair_y.first; - auto value_y = pair_y.second; - auto it_yp = dependency_yp.find(index_y); - if (it_yp != dependency_yp.end()) - { - auto value_yp = it_yp->second; - dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); - } - else - { - dependencies[i].insert(std::make_pair(index_y, value_y)); - } - } - - for (const auto& pair_yp : dependency_yp) - { - auto index_yp = pair_yp.first; - auto value_yp = pair_yp.second; - auto it_y = dependency_y.find(index_yp); - if (it_y == dependency_y.end()) - { - dependencies[i].insert(std::make_pair(index_yp, value_yp)); - } - } - } + gen.evaluateJacobian(); + auto* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: GenClassical DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian() @@ -553,25 +491,25 @@ namespace GridKit bus.allocate(); gen.allocate(); - bus.initialize(); - gen.initialize(); - - gen.updateTime(0.0, 1.0); - for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + gen.size()); bus.setResidualIndex(i, i + gen.size()); } + bus.initialize(); + gen.initialize(); + + gen.updateTime(0.0, 1.0); + bus.evaluateResidual(); gen.evaluateResidual(); bus.evaluateJacobian(); gen.evaluateJacobian(); gen.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gen.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: GenClassical Jacobian\n"; + auto* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: GenClassical Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 1d81ae2f6..42ba9feff 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -362,94 +362,27 @@ namespace GridKit bus.allocate(); gen.allocate(); - // Get d/dy - bus.initialize(); - gen.initialize(); - - auto* gen_y = gen.y().getData(); - for (size_t i = 0; i < gen.size(); ++i) - { - gen_y[i].setVariableNumber(i); ///< Generator independent variables - } - gen.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + gen.size()); // Bus independent variables + bus.setVariableIndex(i, i + gen.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + gen.size()); // Reset bus residual indices } - bus.y().setDataUpdated(); - - bus.evaluateResidual(); - gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_y_view = gen.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); - // Get d/dy' bus.initialize(); gen.initialize(); - auto* gen_yp = gen.yp().getData(); - for (size_t i = 0; i < gen.size(); ++i) - { - gen_yp[i].setVariableNumber(i); ///< Generator independent variables - } - gen.yp().setDataUpdated(); + gen.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term bus.evaluateResidual(); - gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_yp_view = gen.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); - - // Print the dependencies - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - // Extract the dependencies and add d/dy' to d/dy - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); - DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - - for (const auto& pair_y : dependency_y) - { - auto index_y = pair_y.first; - auto value_y = pair_y.second; - auto it_yp = dependency_yp.find(index_y); - if (it_yp != dependency_yp.end()) - { - auto value_yp = it_yp->second; - dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); - } - else - { - dependencies[i].insert(std::make_pair(index_y, value_y)); - } - } + gen.evaluateResidual(); - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) - { - auto index_yp = pair_yp.first; - auto value_yp = pair_yp.second; - auto it_y = dependency_y.find(index_yp); - if (it_y == dependency_y.end()) - { - dependencies[i].insert(std::make_pair(index_yp, value_yp)); - } - } - } + bus.evaluateJacobian(); + gen.evaluateJacobian(); + auto* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Genrou DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian() @@ -480,25 +413,25 @@ namespace GridKit bus.allocate(); gen.allocate(); - bus.initialize(); - gen.initialize(); - - gen.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term - for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + gen.size()); // Reset bus variable indices bus.setResidualIndex(i, i + gen.size()); // Reset bus residual indices } + bus.initialize(); + gen.initialize(); + + gen.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term + bus.evaluateResidual(); gen.evaluateResidual(); bus.evaluateJacobian(); gen.evaluateJacobian(); gen.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gen.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Genrou Jacobian\n"; + auto* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Genrou Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/GensalTests.hpp b/tests/UnitTests/PhasorDynamics/GensalTests.hpp index ba62c8f65..bf676ef8e 100644 --- a/tests/UnitTests/PhasorDynamics/GensalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GensalTests.hpp @@ -364,89 +364,27 @@ namespace GridKit bus.allocate(); gen.allocate(); - bus.initialize(); - gen.initialize(); - - auto* gen_y = gen.y().getData(); - for (size_t i = 0; i < gen.size(); ++i) - { - gen_y[i].setVariableNumber(i); - } - gen.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + gen.size()); + bus.setVariableIndex(i, i + gen.size()); + bus.setResidualIndex(i, i + gen.size()); } - bus.y().setDataUpdated(); - - bus.evaluateResidual(); - gen.evaluateResidual(); - auto& residual_y_view = gen.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); bus.initialize(); gen.initialize(); - auto* gen_yp = gen.yp().getData(); - for (size_t i = 0; i < gen.size(); ++i) - { - gen_yp[i].setVariableNumber(i); - } - gen.yp().setDataUpdated(); + gen.updateTime(0.0, 1.0); bus.evaluateResidual(); gen.evaluateResidual(); - auto& residual_yp_view = gen.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); - - // Print the dependencies - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); - DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - for (const auto& pair_y : dependency_y) - { - auto index_y = pair_y.first; - auto value_y = pair_y.second; - auto it_yp = dependency_yp.find(index_y); - if (it_yp != dependency_yp.end()) - { - auto value_yp = it_yp->second; - dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); - } - else - { - dependencies[i].insert(std::make_pair(index_y, value_y)); - } - } - - for (const auto& pair_yp : dependency_yp) - { - auto index_yp = pair_yp.first; - auto value_yp = pair_yp.second; - auto it_y = dependency_y.find(index_yp); - if (it_y == dependency_y.end()) - { - dependencies[i].insert(std::make_pair(index_yp, value_yp)); - } - } - } + bus.evaluateJacobian(); + gen.evaluateJacobian(); + auto* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Gensal DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; - } + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian() { @@ -459,25 +397,25 @@ namespace GridKit bus.allocate(); gen.allocate(); - bus.initialize(); - gen.initialize(); - - gen.updateTime(0.0, 1.0); - for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + gen.size()); bus.setResidualIndex(i, i + gen.size()); } + bus.initialize(); + gen.initialize(); + + gen.updateTime(0.0, 1.0); + bus.evaluateResidual(); gen.evaluateResidual(); bus.evaluateJacobian(); gen.evaluateJacobian(); gen.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gen.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Gensal Jacobian\n"; + auto* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Gensal Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp index e2f87c43c..d1cf326a0 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp @@ -823,9 +823,9 @@ namespace GridKit success *= (selector.evaluate() == 0); const DependencyTracking::Variable::DependencyMap expected{ - {index(Internal::VLOAD), 0.5}, - {index(Internal::VTEMP), 0.5}, - {index(Internal::VLV), -1.0}, + {2 * index(Internal::VLOAD), 0.5}, // @todo Remove these + {2 * index(Internal::VTEMP), 0.5}, // @todo Remove these + {2 * index(Internal::VLV), -1.0}, // @todo Remove these }; success *= jacobianRowMatches( selector.gastpti.getResidual().getData()[index(Internal::VLV)].getDependencies(), @@ -858,7 +858,7 @@ namespace GridKit const auto dependency_jacobian = dependencyTrackingJacobian(case_data, pmech, success, overrides); const auto enzyme_jacobian = - enzymeJacobian(case_data, pmech, success, overrides, context); + enzymeJacobian(case_data, pmech, success, overrides); success *= jacobianMatches(enzyme_jacobian, dependency_jacobian, @@ -1607,6 +1607,7 @@ namespace GridKit return success; } + /// @todo Remove and setup the test to not rely on explicit variable numbering void numberVariables(Fixture& fixture) const { auto* y = fixture.gastpti.y().getData(); @@ -1615,13 +1616,13 @@ namespace GridKit const auto model_size = static_cast(fixture.gastpti.size()); for (size_t i = 0; i < model_size; ++i) { - y[i].setVariableNumber(i); - yp[i].setVariableNumber(i); + y[i].setVariableNumber(2 * i); + yp[i].setVariableNumber(2 * i + 1); } for (auto variant : Utilities::enum_values()) { const auto port = static_cast(variant); - fixture.input(port).setVariableNumber(fixture.inputIndex(port)); + fixture.input(port).setVariableNumber(2 * fixture.inputIndex(port)); } fixture.gastpti.y().setDataUpdated(); @@ -1643,16 +1644,11 @@ namespace GridKit setAnswerKeyState(fixture.gastpti); setState(fixture.gastpti, overrides); numberVariables(fixture); + fixture.gastpti.updateTime(0.0, 1.0); success *= (fixture.evaluate() == 0); + success *= (fixture.gastpti.evaluateJacobian() == 0); - const auto model_size = static_cast(fixture.gastpti.size()); - std::vector rows(model_size); - const auto* f = fixture.gastpti.getResidual().getData(); - for (size_t i = 0; i < model_size; ++i) - { - rows[i] = f[i].getDependencies(); - } - return rows; + return MapFromCsr(fixture.gastpti.getCsrJacobian()); } #ifdef GRIDKIT_ENABLE_ENZYME @@ -1660,8 +1656,7 @@ namespace GridKit const Data& data, RealT pmech, TestStatus& success, - std::initializer_list overrides, - const char* context) const + std::initializer_list overrides) const { Fixture fixture(data); fixture.attachAllInputs(); @@ -1670,19 +1665,10 @@ namespace GridKit setAnswerKeyState(fixture.gastpti); setState(fixture.gastpti, overrides); fixture.gastpti.updateTime(0.0, 1.0); - if (fixture.evaluate() != 0 || fixture.gastpti.evaluateJacobian() != 0) - { - std::cout << "GASTPTI Jacobian evaluation failed for " << context << '\n'; - success = false; - return {}; - } + success *= (fixture.evaluate() == 0); + success *= (fixture.gastpti.evaluateJacobian() == 0); + success *= (fixture.gastpti.constructCsr() == 0); - if (fixture.gastpti.constructCsr() != 0) - { - std::cout << "GASTPTI CSR construction failed for " << context << '\n'; - success = false; - return {}; - } return MapFromCsr(fixture.gastpti.getCsrJacobian()); } diff --git a/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp index 07b9b5d39..8eb99c2c4 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp @@ -770,8 +770,8 @@ namespace GridKit const auto& dependencies = blocked.hygov.getResidual().getData()[static_cast(Internal::C)].getDependencies(); const DepVar::DependencyMap expected{{ - {static_cast(Internal::C), -1.0}, - {static_cast(Internal::RC), 0.0}, + {2 * static_cast(Internal::C), -1.0}, // @todo Remove these + {2 * static_cast(Internal::RC), 0.0}, // @todo Remove these }}; success *= isEqual(dependencies, expected, kTol); } @@ -1556,6 +1556,7 @@ namespace GridKit return false; } + /// @todo Remove and setup the test to not rely on explicit variable numbering void numberVariables(Fixture& fixture) const { auto* y = fixture.hygov.y().getData(); @@ -1564,12 +1565,12 @@ namespace GridKit const auto model_size = static_cast(fixture.hygov.size()); for (size_t i = 0; i < model_size; ++i) { - y[i].setVariableNumber(i); - yp[i].setVariableNumber(i); + y[i].setVariableNumber(2 * i); + yp[i].setVariableNumber(2 * i); } for (auto port : Utilities::enum_values()) { - fixture.input(port).setVariableNumber(fixture.inputIndex(port)); + fixture.input(port).setVariableNumber(2 * fixture.inputIndex(port)); } fixture.hygov.y().setDataUpdated(); @@ -1608,16 +1609,11 @@ namespace GridKit setAnswerKeyState(fixture.hygov); setState(fixture.hygov, {{Internal::G, gate}}); numberVariables(fixture); + fixture.hygov.updateTime(0.0, 1.0); success *= (fixture.evaluate() == 0); + success *= (fixture.hygov.evaluateJacobian() == 0); - const auto model_size = static_cast(fixture.hygov.size()); - std::vector rows(model_size); - const auto* f = fixture.hygov.getResidual().getData(); - for (size_t i = 0; i < model_size; ++i) - { - rows[i] = f[i].getDependencies(); - } - return rows; + return MapFromCsr(fixture.hygov.getCsrJacobian()); } std::vector enzymeJacobian( @@ -1635,6 +1631,7 @@ namespace GridKit success *= (fixture.evaluate() == 0); success *= (fixture.hygov.evaluateJacobian() == 0); success *= (fixture.hygov.constructCsr() == 0); + return MapFromCsr(fixture.hygov.getCsrJacobian()); } #endif diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 532a21d81..1c401eda2 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -310,95 +310,24 @@ namespace GridKit gov.allocate(); gen.allocate(); - // Get d/dy - bus.initialize(); - gen.initialize(); - gov.initialize(); - - auto* gov_y = gov.y().getData(); - for (size_t i = 0; i < gov.size(); ++i) - { - gov_y[i].setVariableNumber(i); // Governor independent variables - } - gov.y().setDataUpdated(); - auto* gen_y = gen.y().getData(); - gen_y[1].setVariableNumber(gov.size()); // omega as an additional independent variable - gen.y().setDataUpdated(); - - bus.evaluateResidual(); - gen.evaluateResidual(); - gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking - // the dependencies - auto& residual_y_view = gov.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); + gen.setVariableIndex(1, gov.size()); // Reset omega index - // Get d/dy' bus.initialize(); gen.initialize(); gov.initialize(); - auto* gov_yp = gov.yp().getData(); - for (size_t i = 0; i < gov.size(); ++i) - { - gov_yp[i].setVariableNumber(i); ///< Governor independent variables - } - gov.yp().setDataUpdated(); + gov.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term bus.evaluateResidual(); gen.evaluateResidual(); - gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking - // the dependencies - auto& residual_yp = gov.getResidual(); - const auto* residual_yp_data = residual_yp.getData(); - - // Print the dependencies - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - residual_yp_data[i].print(std::cout); - std::cout << "\n"; - } - - // Extract the dependencies and add d/dy' to d/dy - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); - DependencyTracking::Variable::DependencyMap dependency_yp = residual_yp_data[i].getDependencies(); - - for (const auto& pair_y : dependency_y) - { - auto index_y = pair_y.first; - auto value_y = pair_y.second; - auto it_yp = dependency_yp.find(index_y); - if (it_yp != dependency_yp.end()) - { - auto value_yp = it_yp->second; - dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); - } - else - { - dependencies[i].insert(std::make_pair(index_y, value_y)); - } - } + gov.evaluateResidual(); - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) - { - auto index_yp = pair_yp.first; - auto value_yp = pair_yp.second; - auto it_y = dependency_y.find(index_yp); - if (it_y == dependency_y.end()) - { - dependencies[i].insert(std::make_pair(index_yp, value_yp)); - } - } - } + gov.evaluateJacobian(); + auto* model_jacobian = gov.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Tgov1 DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian( @@ -430,8 +359,8 @@ namespace GridKit gov.evaluateJacobian(); gov.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gov.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Tgov1 Jacobian\n"; + auto* model_jacobian = gov.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Tgov1 Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp b/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp index b914a601b..4178cc24f 100644 --- a/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp +++ b/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp @@ -224,73 +224,25 @@ namespace GridKit bus.allocate(); load.allocate(); - bus.initialize(); - load.initialize(); - - auto* load_y = load.y().getData(); - for (size_t i = 0; i < load.size(); ++i) - { - load_y[i].setVariableNumber(i); - } - load.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + load.size()); + bus.setVariableIndex(i, i + load.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + load.size()); // Reset bus residual indices } - bus.y().setDataUpdated(); - - bus.evaluateResidual(); - load.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_y_view = load.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); bus.initialize(); load.initialize(); - auto* load_yp = load.yp().getData(); - for (size_t i = 0; i < load.size(); ++i) - { - load_yp[i].setVariableNumber(i); - } - load.yp().setDataUpdated(); + load.updateTime(0.0, 1.0); bus.evaluateResidual(); - load.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies - auto& residual_yp_view = load.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); - - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - auto dependency_y = residual_y[i].getDependencies(); - auto dependency_yp = residual_yp[i].getDependencies(); - - for (const auto& pair_y : dependency_y) - { - auto it_yp = dependency_yp.find(pair_y.first); - if (it_yp != dependency_yp.end()) - { - dependencies[i].insert(std::make_pair(pair_y.first, pair_y.second + it_yp->second)); - } - else - { - dependencies[i].insert(std::make_pair(pair_y.first, pair_y.second)); - } - } - - for (const auto& pair_yp : dependency_yp) - { - if (!dependency_y.contains(pair_yp.first)) - { - dependencies[i].insert(std::make_pair(pair_yp.first, pair_yp.second)); - } - } - } + load.evaluateResidual(); //< Tracks dependencies + load.evaluateJacobian(); //< Converts dependencies to CSR + auto* model_jacobian = load.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: LoadZIP DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian( @@ -305,25 +257,25 @@ namespace GridKit bus.allocate(); load.allocate(); + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + load.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + load.size()); // Reset bus residual indices + } + bus.initialize(); load.initialize(); load.updateTime(0.0, 1.0); - for (size_t i = 0; i < bus.size(); ++i) - { - bus.setVariableIndex(i, i + load.size()); - bus.setResidualIndex(i, i + load.size()); - } - bus.evaluateResidual(); load.evaluateResidual(); bus.evaluateJacobian(); load.evaluateJacobian(); load.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = load.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: LoadZIP Jacobian\n"; + auto* model_jacobian = load.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: LoadZIP Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); @@ -331,7 +283,7 @@ namespace GridKit #endif private: - static constexpr RealT tol_ = 1.0e-10; + static constexpr RealT tol_ = 10 * std::numeric_limits::epsilon(); auto makeData() -> DataT { diff --git a/tests/UnitTests/PhasorDynamics/LoadZTests.hpp b/tests/UnitTests/PhasorDynamics/LoadZTests.hpp index c0ac99d78..20c7d9c73 100644 --- a/tests/UnitTests/PhasorDynamics/LoadZTests.hpp +++ b/tests/UnitTests/PhasorDynamics/LoadZTests.hpp @@ -101,36 +101,28 @@ namespace GridKit bus.allocate(); load.allocate(); - bus.initialize(); - load.initialize(); - - auto* load_y = load.y().getData(); - for (size_t i = 0; i < load.size(); ++i) - { - load_y[i].setVariableNumber(i); ///< load independent variables - } - load.y().setDataUpdated(); - auto* bus_y = bus.y().getData(); for (size_t i = 0; i < bus.size(); ++i) { - bus_y[i].setVariableNumber(i + load.size()); // Bus independent variables + bus.setVariableIndex(i, i + load.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + load.size()); // Reset bus residual indices } - bus.y().setDataUpdated(); - bus.evaluateResidual(); - load.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking - ///< the dependencies + bus.initialize(); + load.initialize(); - auto& residuals = load.getResidual(); - const auto* residual_data = residuals.getData(); - std::vector ref = analyticalJacobian(R, X); + bus.evaluateResidual(); + load.evaluateResidual(); //< Tracks dependencies + load.evaluateJacobian(); //< Converts dependencies to CSR + auto* model_jacobian = load.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Load DependencyTracking Jacobian\n"; + model_jacobian->print(); - /// Compare dependencies computed automatically to the ones computed analytically - for (size_t i = 0; i < residuals.getSize(); ++i) + // Compare model Jacobian wih dependencies computed analytically + auto ref = analyticalJacobian(R, X); + auto model_dependencies = GridKit::Testing::MapFromCsr(model_jacobian); + for (size_t i = 0; i < ref.size(); ++i) { - DependencyTracking::Variable res = residual_data[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, ref[i])); + success *= (GridKit::Testing::isEqual(model_dependencies[i], ref[i])); } return success.report(__func__); @@ -190,23 +182,26 @@ namespace GridKit bus.allocate(); load.allocate(); - bus.initialize(); - load.initialize(); - for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + load.size()); // Reset bus variable indices bus.setResidualIndex(i, i + load.size()); // Reset bus residual indices } + bus.initialize(); + load.initialize(); + + bus.evaluateResidual(); + load.evaluateResidual(); + bus.evaluateJacobian(); load.evaluateJacobian(); load.constructCsr(); - GridKit::LinearAlgebra::CsrMatrix* model_jacobian = load.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Load Jacobian\n"; + auto* model_jacobian = load.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Load Enzyme Jacobian\n"; model_jacobian->print(); - /// Compare model Jacobian wih dependencies computed analytically + // Compare model Jacobian wih dependencies computed analytically std::vector ref = analyticalJacobian(R, X); std::vector model_dependencies = GridKit::Testing::MapFromCsr(model_jacobian); for (size_t i = 0; i < ref.size(); ++i) @@ -219,7 +214,7 @@ namespace GridKit #endif private: - static constexpr RealT tol_ = 1.0e-10; + static constexpr RealT tol_ = 10 * std::numeric_limits::epsilon(); auto makeData() -> DataT { diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index c26e24a35..99ba5b4d0 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -219,79 +219,19 @@ namespace GridKit stab.allocate(); stab.initialize(); - // --- d/dy: tag internal variables as independent --- - auto* y = stab.y().getData(); - for (size_t i = 0; i < stab.size(); ++i) - { - y[i].setVariableNumber(i); - } // Tag external signal u as an additional independent variable - u_value.setVariableNumber(stab.size()); - u_value.setValue(0.5); - + u_value.setVariableNumber(2 * stab.size()); // @todo avoid requiring knowledge of the numbering setStatePointDep(stab); - stab.evaluateResidual(); - auto& residual_y_view = stab.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); - - // --- d/dy': tag derivatives as independent --- - u_value = 0.5; - stab.initialize(); - auto* yp = stab.yp().getData(); - for (size_t i = 0; i < stab.size(); ++i) - { - yp[i].setVariableNumber(i); - } - - setStatePointDep(stab); + stab.updateTime(0.0, 1.0); // alpha = 1.0 to verify d/dy' term stab.evaluateResidual(); - auto& residual_yp_view = stab.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); - - // Print dependencies for debugging - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - // Merge d/dy and d/dy' into a single dependency map - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - auto dependency_y = (residual_y[i]).getDependencies(); - auto dependency_yp = (residual_yp[i]).getDependencies(); - - for (const auto& pair_y : dependency_y) - { - auto it_yp = dependency_yp.find(pair_y.first); - if (it_yp != dependency_yp.end()) - { - dependencies[i].insert(std::make_pair(pair_y.first, pair_y.second + it_yp->second)); - } - else - { - dependencies[i].insert(std::make_pair(pair_y.first, pair_y.second)); - } - } - - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) - { - if (!dependency_y.contains(pair_yp.first)) - { - dependencies[i].insert(std::make_pair(pair_yp.first, pair_yp.second)); - } - } - } + stab.evaluateJacobian(); + auto model_jacobian = stab.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Ieeest DependencyTracking Jacobian\n"; + model_jacobian->print(); - return dependencies; + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian( @@ -323,7 +263,7 @@ namespace GridKit stab.evaluateJacobian(); stab.constructCsr(); auto model_jacobian = stab.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Ieeest Jacobian\n"; + std::cout << "Sparse Csr Matrix: Ieeest Enzyme Jacobian\n"; model_jacobian->print(); return GridKit::Testing::MapFromCsr(model_jacobian); diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 21145ff90..b7a1b433d 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -497,7 +497,7 @@ namespace GridKit // Evaluate and get the system Jacobian system.evaluateResidual(); system.evaluateJacobian(); - GridKit::LinearAlgebra::CsrMatrix* system_jacobian = system.getCsrJacobian(); + auto* system_jacobian = system.getCsrJacobian(); std::cout << "Sparse Csr Matrix: System Jacobian with DependencyTracking\n"; system_jacobian->print(); @@ -517,7 +517,7 @@ namespace GridKit // Evaluate and get the system Jacobian system.evaluateResidual(); system.evaluateJacobian(); - GridKit::LinearAlgebra::CsrMatrix* system_jacobian = system.getCsrJacobian(); + auto* system_jacobian = system.getCsrJacobian(); std::cout << "Sparse Csr Matrix: System Jacobian with Enzyme\n"; system_jacobian->print(); From d2f982e4a41aafe4ebebe6544a8dabbf3f9afbd0 Mon Sep 17 00:00:00 2001 From: nkoukpaizan Date: Wed, 9 Sep 2026 20:27:48 +0000 Subject: [PATCH 12/12] Apply pre-commit fixes --- .../DependencyTracking/Variable.hpp | 6 +++--- .../Exciter/ESDC1A/Esdc1aDependencyTracking.cpp | 2 +- GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp | 2 +- tests/UnitTests/PhasorDynamics/BusFaultTests.hpp | 2 +- tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp | 6 +++--- tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp | 4 ++-- tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp | 6 +++--- tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp | 3 ++- tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp | 2 +- tests/UnitTests/PhasorDynamics/GenrouTests.hpp | 2 +- tests/UnitTests/PhasorDynamics/GensalTests.hpp | 3 ++- tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp | 2 +- tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp | 2 +- tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp | 2 +- 14 files changed, 23 insertions(+), 21 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp index 3f9bc2cfe..5d90ce0a2 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp @@ -257,7 +257,7 @@ namespace GridKit @brief Turns variable into parameter, or vice versa. @todo is_fixed_ is currently not contributing to the semantics of - the derivatives. Leaving as-is for now, as it is not used + the derivatives. Leaving as-is for now, as it is not used for anything other than printed diagnostics. */ void setFixed(bool b = false) @@ -305,8 +305,8 @@ namespace GridKit size_t variable_number_; ///< Independent variable ID bool is_fixed_; ///< Constant parameter flag. - DependencyMap dependencies_; - static const size_t INVALID_VAR_NUMBER = INVALID_INDEX; + DependencyMap dependencies_; + static const size_t INVALID_VAR_NUMBER = INVALID_INDEX; }; //------------------------------------ diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp index 387677c72..7fd338323 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp @@ -30,7 +30,7 @@ namespace GridKit return 0; } - + // Available template instantiations template class Esdc1a; template class Esdc1a; diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 425a29f91..a4af51239 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -390,7 +390,7 @@ namespace GridKit y_.setDataUpdated(); yp_.setToConst(static_cast(ZERO)); - + // For DependencyTracking::Variable, set variable numbers if constexpr (std::is_same_v) { diff --git a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp index 9fa45e096..8d34e0480 100644 --- a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp @@ -152,7 +152,7 @@ namespace GridKit std::cout << "Sparse Csr Matrix: BusFault DependencyTracking Jacobian\n"; model_jacobian->print(); - return GridKit::Testing::MapFromCsr(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian( diff --git a/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp index 7643d30a6..e6f90ca2e 100644 --- a/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp @@ -1020,9 +1020,9 @@ namespace GridKit const DependencyTracking::Variable::DependencyMap expected{ {2 * index(Vars::XQPI) + 1, -1.0}, // @todo Remove these - {2 * index(Vars::SFRZ), 0.0}, // @todo Remove these - {2 * index(Vars::ERQLIM), 0.0}, // @todo Remove these - {2 * index(Vars::QPI), 0.0}, // @todo Remove these + {2 * index(Vars::SFRZ), 0.0}, // @todo Remove these + {2 * index(Vars::ERQLIM), 0.0}, // @todo Remove these + {2 * index(Vars::QPI), 0.0}, // @todo Remove these }; success *= jacobianRowMatches( diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index cb7ba7f87..1288eee31 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -661,7 +661,7 @@ namespace GridKit const DepVar::DependencyMap expected{{ {2 * index(Vars::VT), 0.5 * kHvrcmGain}, // @todo remove these - {2 * index(Vars::IQEXTRA), -1.0}, // @todo remove these + {2 * index(Vars::IQEXTRA), -1.0}, // @todo remove these }}; success *= isEqual(dependencies, expected, kTol); } @@ -694,7 +694,7 @@ namespace GridKit dependencyTrackingJacobian(data, current, success); const auto enzyme_jacobian = enzymeJacobian(data, current, success); - const auto nrows = std::min(dependency_tracking_jacobian.size(), + const auto nrows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); for (size_t i = 0; i < nrows; ++i) diff --git a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp index 3efa7d742..05fc850b0 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp @@ -804,10 +804,10 @@ namespace GridKit const auto& dependencies = transition.esdc1a.getResidual().getData()[static_cast(Internal::EFDP)].getDependencies(); const DepVar::DependencyMap expected{{ - {2 * static_cast(Internal::EFDP), -12.0}, // @todo Remove these + {2 * static_cast(Internal::EFDP), -12.0}, // @todo Remove these {2 * static_cast(Internal::EFDP) + 1, -1.0}, // @todo Remove these - {2 * static_cast(Internal::VR), 1.0}, // @todo Remove these - {2 * static_cast(Internal::VFE), -1.0}, // @todo Remove these + {2 * static_cast(Internal::VR), 1.0}, // @todo Remove these + {2 * static_cast(Internal::VFE), -1.0}, // @todo Remove these }}; success *= isEqual(dependencies, expected, kTol); } diff --git a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp index b8d15e58a..3297ceac6 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp @@ -323,7 +323,8 @@ namespace GridKit std::cout << "Sparse Csr Matrix: Ieeet1 DependencyTracking Jacobian\n"; model_jacobian->print(); - return GridKit::Testing::MapFromCsr(model_jacobian); } + return GridKit::Testing::MapFromCsr(model_jacobian); + } std::vector EnzymeJacobian() { diff --git a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp index 0388d5742..283c1646d 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp @@ -450,7 +450,7 @@ namespace GridKit std::cout << "Sparse Csr Matrix: SexsPti DependencyTracking Jacobian\n"; model_jacobian->print(); - return GridKit::Testing::MapFromCsr(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian() diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 42ba9feff..164cbf038 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -382,7 +382,7 @@ namespace GridKit std::cout << "Sparse Csr Matrix: Genrou DependencyTracking Jacobian\n"; model_jacobian->print(); - return GridKit::Testing::MapFromCsr(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian() diff --git a/tests/UnitTests/PhasorDynamics/GensalTests.hpp b/tests/UnitTests/PhasorDynamics/GensalTests.hpp index bf676ef8e..c02c3263b 100644 --- a/tests/UnitTests/PhasorDynamics/GensalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GensalTests.hpp @@ -384,7 +384,8 @@ namespace GridKit std::cout << "Sparse Csr Matrix: Gensal DependencyTracking Jacobian\n"; model_jacobian->print(); - return GridKit::Testing::MapFromCsr(model_jacobian); } + return GridKit::Testing::MapFromCsr(model_jacobian); + } std::vector EnzymeJacobian() { diff --git a/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp index d1cf326a0..3231b0561 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp @@ -825,7 +825,7 @@ namespace GridKit const DependencyTracking::Variable::DependencyMap expected{ {2 * index(Internal::VLOAD), 0.5}, // @todo Remove these {2 * index(Internal::VTEMP), 0.5}, // @todo Remove these - {2 * index(Internal::VLV), -1.0}, // @todo Remove these + {2 * index(Internal::VLV), -1.0}, // @todo Remove these }; success *= jacobianRowMatches( selector.gastpti.getResidual().getData()[index(Internal::VLV)].getDependencies(), diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 1c401eda2..6ee89b191 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -327,7 +327,7 @@ namespace GridKit std::cout << "Sparse Csr Matrix: Tgov1 DependencyTracking Jacobian\n"; model_jacobian->print(); - return GridKit::Testing::MapFromCsr(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian( diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 99ba5b4d0..24e959555 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -231,7 +231,7 @@ namespace GridKit std::cout << "Sparse Csr Matrix: Ieeest DependencyTracking Jacobian\n"; model_jacobian->print(); - return GridKit::Testing::MapFromCsr(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } std::vector EnzymeJacobian(