diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp index 6e2f216d0..5d90ce0a2 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp @@ -13,7 +13,6 @@ #include #include #include -#include #include #include @@ -209,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; } /** @@ -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) { @@ -265,10 +269,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); @@ -305,8 +305,8 @@ namespace GridKit size_t variable_number_; ///< Independent variable ID bool is_fixed_; ///< Constant parameter flag. - mutable 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/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. */ 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 707902dfc..865864374 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. + */ int constructCsr() { - if (coo_jac_ == nullptr) + if constexpr (std::is_same_v) { - constructCoo(); + return constructCsrFromDependencies(); } - - if (csr_jac_ == nullptr) + else { - 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 constructCsrFromCoo(); } - - return 0; } protected: @@ -309,13 +301,17 @@ namespace GridKit */ void allocateVectors(IdxT n) { - y_.resize(n); yp_.resize(n); f_.resize(n); abs_tol_.resize(n); } + /** + * @brief COO construction from component-level raw buffers. + * + * @note the components retain ownership of the data in the raw buffers. + */ int constructCoo() { if (coo_jac_ == nullptr) @@ -340,6 +336,180 @@ namespace GridKit return 0; } + /** + * @brief CSR construction from COO. + * + * @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() + { + 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; + } + + /** + * @brief CSR construction from Dependency maps. + * + * @note Currently only used for testing. + * See \ref initializeDependencyTrackingVariableNumbers() + */ + int constructCsrFromDependencies() + requires std::is_same_v + { + 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 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)]; + 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 size_t 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 size_t 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; + } + + /** + * @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 size_{0}; IdxT nnz_{0}; /// Global (system-level) variable indices 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/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 0b8fd4615..2ed478014 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 @@ -319,6 +316,13 @@ namespace GridKit } commitInitialPoint(point); + + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return 0; } @@ -709,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..7fd338323 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 2645ed940..a4af51239 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 @@ -379,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()) { @@ -391,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 ada67c9e3..6be4cf1eb 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 @@ -301,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()) { @@ -315,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 d5a33b257..e35cc8a2a 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 @@ -388,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_); @@ -396,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/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/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 e85d05354..71ada99e0 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 @@ -401,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()) { @@ -421,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 229a58df3..01843aa0f 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 @@ -299,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_); @@ -308,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/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..f6105911e 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -7,15 +7,32 @@ 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. */ template bool SystemModel::hasJacobian() { + Log::warning() << "DependencyTracking::Variable Jacobians are only available for testing.\n" + << "Falling back to dense Jacobians for PhasorDyanmics simulations.\n"; + return false; } + /** + * @brief Evaluate system DependencyTracking::Variable Jacobian. + * + * @note Currently only used for testing. + */ + template + int SystemModel::evaluateJacobian() + { + this->constructCsr(); + + // 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..14300ccee 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp @@ -34,6 +34,192 @@ 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..598f1d23d 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -560,6 +560,12 @@ namespace GridKit y_.setDataUpdated(); yp_.setDataUpdated(); + // For DependencyTracking::Variable, set variable numbers + if constexpr (std::is_same_v) + { + this->initializeDependencyTrackingVariableNumbers(); + } + return status; } @@ -703,193 +709,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 * diff --git a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp index 03528173e..8d34e0480 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..e6f90ca2e 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..1288eee31 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,13 +694,7 @@ 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(), + const auto nrows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); for (size_t i = 0; i < nrows; ++i) @@ -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..05fc850b0 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..3297ceac6 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp @@ -303,94 +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 +336,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..283c1646d 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..164cbf038 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..c02c3263b 100644 --- a/tests/UnitTests/PhasorDynamics/GensalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GensalTests.hpp @@ -364,88 +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 +398,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..3231b0561 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..6ee89b191 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..24e959555 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 afc829795..b7a1b433d 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(); + auto* 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( @@ -538,8 +517,8 @@ namespace GridKit // Evaluate and get the system Jacobian system.evaluateResidual(); system.evaluateJacobian(); - GridKit::LinearAlgebra::CsrMatrix* system_jacobian = system.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: System Jacobian\n"; + auto* system_jacobian = system.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: System Jacobian with Enzyme\n"; system_jacobian->print(); return GridKit::Testing::MapFromCsr(system_jacobian);