From 10a92546e76baf9973b459e839b22888f2a7c00c Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 21 May 2026 22:05:19 -0500 Subject: [PATCH 1/3] Use global state views in PhasorDynamics SystemModel --- .../PhasorDynamics/Branch/BranchEnzyme.cpp | 16 +- GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp | 4 +- .../BusFault/BusFaultEnzyme.cpp | 2 +- .../Exciter/IEEET1/Ieeet1Enzyme.cpp | 2 +- .../Exciter/IEEET1/Ieeet1Impl.hpp | 4 +- .../Exciter/SEXS-PTI/SexsPtiEnzyme.cpp | 2 +- .../Exciter/SEXS-PTI/SexsPtiImpl.hpp | 4 +- .../Governor/Tgov1/Tgov1Impl.hpp | 4 +- GridKit/Model/PhasorDynamics/GridElement.hpp | 130 +++++++++++++- .../Model/PhasorDynamics/Load/LoadEnzyme.cpp | 2 +- .../Model/PhasorDynamics/Load/LoadImpl.hpp | 4 +- .../PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp | 2 +- .../PhasorDynamics/LoadZIP/LoadZIPImpl.hpp | 4 +- .../Stabilizer/IEEEST/IeeestImpl.hpp | 4 +- .../GENROUwS/GenrouEnzyme.cpp | 2 +- .../GENROUwS/GenrouImpl.hpp | 4 +- .../GENSALwS/GensalEnzyme.cpp | 2 +- .../GENSALwS/GensalImpl.hpp | 4 +- .../GenClassical/GenClassicalEnzyme.cpp | 2 +- .../GenClassical/GenClassicalImpl.hpp | 4 +- GridKit/Model/PhasorDynamics/SystemModel.hpp | 164 +++++++----------- .../UnitTests/PhasorDynamics/SystemTests.hpp | 50 ++++++ .../PhasorDynamics/runSystemTests.cpp | 1 + 23 files changed, 260 insertions(+), 157 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp index 535d944cd..523adf6a6 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp @@ -39,12 +39,12 @@ namespace GridKit ScalarT, IdxT>::eval(this, static_cast(bus1_->size()), - static_cast((bus1_->y()).size()), + static_cast(bus1_->size()), (bus1_->getResidualIndices()).data(), (bus1_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus1_->y()).data(), + bus1_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, @@ -56,12 +56,12 @@ namespace GridKit ScalarT, IdxT>::eval(this, static_cast(bus2_->size()), - static_cast((bus2_->y()).size()), + static_cast(bus2_->size()), (bus2_->getResidualIndices()).data(), (bus2_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus2_->y()).data(), + bus2_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, @@ -73,12 +73,12 @@ namespace GridKit ScalarT, IdxT>::eval(this, static_cast(bus1_->size()), - static_cast((bus2_->y()).size()), + static_cast(bus2_->size()), (bus1_->getResidualIndices()).data(), (bus2_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus2_->y()).data(), + bus2_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, @@ -91,12 +91,12 @@ namespace GridKit ScalarT, IdxT>::eval(this, static_cast(bus2_->size()), - static_cast((bus1_->y()).size()), + static_cast(bus1_->size()), (bus2_->getResidualIndices()).data(), (bus1_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus1_->y()).data(), + bus1_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 77fe6ecf0..3731dfdb7 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -82,9 +82,7 @@ namespace GridKit size_t size = static_cast(size_); // Resize component model data - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp index d607ead65..c58282b76 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp @@ -43,7 +43,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp index 9bb4496d6..87da765b2 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp @@ -65,7 +65,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), ws_.data(), J_rows_buffer_, J_cols_buffer_, diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index f535862f0..6995671af 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -118,9 +118,7 @@ namespace GridKit { // Resize component model data auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp index 41610a891..b470e8725 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp @@ -56,7 +56,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), ws_.data(), J_rows_buffer_, J_cols_buffer_, diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index f09fd3b5f..c032dec21 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -58,9 +58,7 @@ namespace GridKit int SexsPti::allocate() { auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 836e9c3ff..8eea53fc2 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -136,9 +136,7 @@ namespace GridKit { // Allocate local component data auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/GridElement.hpp b/GridKit/Model/PhasorDynamics/GridElement.hpp index aa961bfc7..652fe01e7 100644 --- a/GridKit/Model/PhasorDynamics/GridElement.hpp +++ b/GridKit/Model/PhasorDynamics/GridElement.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -14,6 +15,89 @@ namespace GridKit { using Log = ::GridKit::Utilities::Logger; + /// Owns state storage for standalone elements, or aliases a SystemModel-owned global slice. + /// The owned storage remains the std::vector returned through Evaluator; the data accessors + /// are the live state used by component implementations. + template + class StateVector + { + public: + using StorageT = std::vector; + using size_type = typename StorageT::size_type; + + StateVector() = default; + + // Non-copyable/movable: data_ aliases either storage_ or a global slice, and a GridElement + // (which owns these) is never copied or moved -- it is held by raw pointer in SystemModel. + StateVector(const StateVector&) = delete; + StateVector& operator=(const StateVector&) = delete; + StateVector(StateVector&&) = delete; + StateVector& operator=(StateVector&&) = delete; + + template + void resize(SizeT size) + { + const auto sz = static_cast(size); + if (viewing_) + { + assert(sz == size_); + return; + } + + storage_.resize(sz); + data_ = storage_.data(); + size_ = storage_.size(); + } + + void view(ScalarT* data, size_type size) + { + data_ = data; + size_ = size; + viewing_ = true; + } + + ScalarT* data() + { + return data_; + } + + const ScalarT* data() const + { + return data_; + } + + size_type size() const + { + return size_; + } + + ScalarT& operator[](size_type i) + { + return data_[i]; + } + + const ScalarT& operator[](size_type i) const + { + return data_[i]; + } + + StorageT& storage() + { + return storage_; + } + + const StorageT& storage() const + { + return storage_; + } + + private: + StorageT storage_; + ScalarT* data_{nullptr}; + size_type size_{0}; + bool viewing_{false}; + }; + /** * @brief Model base class for all system constituents */ @@ -65,24 +149,26 @@ namespace GridKit msa = max_steps_; } + // Evaluator requires std::vector references. These return owned storage; when an element is + // bound into a SystemModel, live state is accessed through y_/yp_/f_ instead. std::vector& y() override { - return y_; + return y_.storage(); } const std::vector& y() const override { - return y_; + return y_.storage(); } std::vector& yp() override { - return yp_; + return yp_.storage(); } const std::vector& yp() const override { - return yp_; + return yp_.storage(); } std::vector& tag() override @@ -97,12 +183,31 @@ namespace GridKit std::vector& getResidual() override { - return f_; + return f_.storage(); } const std::vector& getResidual() const override { - return f_; + return f_.storage(); + } + + ScalarT* yData() + { + return y_.data(); + } + + const ScalarT* yData() const + { + return y_.data(); + } + + /// Bind this element's y/yp/f to its slice of the SystemModel's global vectors, so the + /// element reads and writes the global state in place (no scatter/gather). + void bindGlobalState(ScalarT* y, ScalarT* yp, ScalarT* f, std::size_t n) + { + y_.view(n == 0 ? nullptr : y, n); + yp_.view(n == 0 ? nullptr : yp, n); + f_.view(n == 0 ? nullptr : f, n); } MatrixT& getJacobian() override @@ -148,6 +253,13 @@ namespace GridKit } protected: + void allocateState(std::size_t size) + { + f_.resize(size); + y_.resize(size); + yp_.resize(size); + } + IdxT size_{0}; IdxT nnz_{0}; /// Global (system-level) variable indices @@ -155,10 +267,10 @@ namespace GridKit /// Global (system-level) residual indices std::vector residual_indices_; - std::vector y_; - std::vector yp_; + StateVector y_; + StateVector yp_; + StateVector f_; std::vector tag_; - std::vector f_; std::vector g_; MatrixT J_; diff --git a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp index 1344a9b18..d62cad276 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp @@ -61,7 +61,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp index e774074ea..7a08fe557 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp @@ -86,9 +86,7 @@ namespace GridKit // std::cout << "Allocate Load..." << std::endl; auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp index 4cc1d0b47..386d74ce5 100644 --- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp @@ -56,7 +56,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp index 6896662d1..93c38a779 100644 --- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp +++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp @@ -102,9 +102,7 @@ namespace GridKit // std::cout << "Allocate Load..." << std::endl; auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index e428c601a..f13545907 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -157,9 +157,7 @@ namespace GridKit int Ieeest::allocate() { auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp index 18c4cbb36..e60819276 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp @@ -97,7 +97,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index b7ccac6d5..eaecd44e8 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -327,9 +327,7 @@ namespace GridKit { // Resize component model data auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp index a70625e1c..c3c091354 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp @@ -97,7 +97,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp index 12227e69e..19e33eb9f 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp @@ -196,9 +196,7 @@ namespace GridKit { // Resize component model data auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp index c38a88931..05757e3dc 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp @@ -62,7 +62,7 @@ namespace GridKit (bus_->getVariableIndices()).data(), y_.data(), yp_.data(), - (bus_->y()).data(), + bus_->yData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index e3836912f..bd77d269b 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -180,9 +180,7 @@ namespace GridKit { // Resize component model data auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); + this->allocateState(size); tag_.resize(size); variable_indices_.resize(size); residual_indices_.resize(size); diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 25478738c..aeba3c42d 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -441,9 +442,9 @@ namespace GridKit /** * @brief Allocate buses, components, and system objects. * - * This method first allocates bus objects, then component objects, - * and computes system size (number of unknowns). Once the size is - * computed, system global objects are allocated. + * This method computes system size from bus/component sizes, allocates + * global objects, binds each element to its global slice, then allocates + * bus and component objects. * * @post size_ >= 1 * @@ -451,35 +452,16 @@ namespace GridKit int allocate() override { size_ = 0; - - // Allocate all buses for (const auto& bus : buses_) { - bus->allocate(); - for (IdxT j = 0; j < bus->size(); ++j) - { - bus->setVariableIndex(j, size_ + j); - bus->setResidualIndex(j, size_ + j); - } size_ += bus->size(); } - - // Allocate all components for (const auto& component : components_) { - component->allocate(); - for (IdxT j = 0; j < component->size(); ++j) - { - component->setVariableIndex(j, size_ + j); - component->setResidualIndex(j, size_ + j); - } size_ += component->size(); } - // Allocate global vectors - y_.resize(size_); - yp_.resize(size_); - f_.resize(size_); + this->allocateState(static_cast(size_)); tag_.resize(size_); variable_indices_.resize(size_); residual_indices_.resize(size_); @@ -491,6 +473,49 @@ namespace GridKit this->setResidualIndex(j, j); } + auto bind = [this](auto* element, IdxT off, IdxT n) + { + const auto o = static_cast(off); + element->bindGlobalState(y_.data() + o, yp_.data() + o, f_.data() + o, static_cast(n)); + }; + auto checkAllocatedSize = [](auto* element, IdxT expected) + { + if (element->size() != expected) + { + Log::error() << "Element size changed during allocation. Expected " << expected + << ", got " << element->size() << ".\n"; + throw std::runtime_error("SystemModel allocation failed"); + } + }; + + IdxT off = 0; + for (const auto& bus : buses_) + { + const IdxT n = bus->size(); + bind(bus, off, n); + bus->allocate(); + checkAllocatedSize(bus, n); + for (IdxT j = 0; j < n; ++j) + { + bus->setVariableIndex(j, off + j); + bus->setResidualIndex(j, off + j); + } + off += n; + } + for (const auto& component : components_) + { + const IdxT n = component->size(); + bind(component, off, n); + component->allocate(); + checkAllocatedSize(component, n); + for (IdxT j = 0; j < n; ++j) + { + component->setVariableIndex(j, off + j); + component->setResidualIndex(j, off + j); + } + off += n; + } + // Verify component configuration int errorCount = this->verify(); if (errorCount > 0) @@ -577,42 +602,24 @@ namespace GridKit * Also, generators may write to control devices (e.g. governors, * exciters, etc.) during the initialization. * - * @todo Implement writting to system vectors in a thread-safe way. + * @todo Make initialization writes to system vectors thread-safe. * - * @note Currently assuming each component stores variables contiguously in memory and - * that these are simply concateneted in the global system. + * @note Each bus and component is bound to a contiguous slice of the global system state. */ int initialize() override { + // Buses initialize first: components may read bus values during their own + // initialization. Writes land directly in the global vector via the bound spans. for (const auto& bus : buses_) { bus->initialize(); } - for (const auto& bus : buses_) - { - for (IdxT j = 0; j < bus->size(); ++j) - { - y_[bus->getVariableIndex(j)] = bus->y()[j]; - yp_[bus->getVariableIndex(j)] = bus->yp()[j]; - } - } - - // Initialize components for (const auto& component : components_) { component->initialize(); } - for (const auto& component : components_) - { - for (IdxT j = 0; j < component->size(); ++j) - { - y_[component->getVariableIndex(j)] = component->y()[j]; - yp_[component->getVariableIndex(j)] = component->yp()[j]; - } - } - return 0; } @@ -684,24 +691,19 @@ namespace GridKit /** * @brief Compute system residual vector * - * First, update bus and component variables from the system solution - * vector. Next, evaluate residuals in buses and components, and - * then copy values to the global residual vector. + * Buses and components are bound to the system state during allocation. + * Residual evaluation therefore reads the current system values and + * writes residuals directly into the global residual vector. * * @warning Residuals must be computed for buses, before component * residuals are computed. Buses own residuals for currents * Ir and Ii, but the contributions to these residuals come * from components. Buses assign their residual values, while components - * add to those values by in-place adition. This is why (for now) bus + * add to those values by in-place addition. This is why (for now) bus * residuals need to be computed first. - * - * @todo Here, components write to local values, which are then copied - * to global system vectors. Make components write to the system - * vectors directly. */ int evaluateResidual() override { - updateVariables(); for (const auto& bus : buses_) { @@ -713,23 +715,6 @@ namespace GridKit component->evaluateResidual(); } - // Update residual vector - for (const auto& bus : buses_) - { - for (IdxT j = 0; j < bus->size(); ++j) - { - f_[bus->getResidualIndex(j)] = bus->getResidual()[j]; - } - } - - for (const auto& component : components_) - { - for (IdxT j = 0; j < component->size(); ++j) - { - f_[component->getResidualIndex(j)] = component->getResidual()[j]; - } - } - return 0; } @@ -767,7 +752,7 @@ namespace GridKit IdxT nnz_dup = 0; for (const auto& component : components_) { - auto component_jacobian = component->getJacobian(); + auto& component_jacobian = component->getJacobian(); std::tuple&, std::vector&, std::vector&> component_jacobian_entries = component_jacobian.getEntries(false); const auto [rows, columns, values] = component_jacobian_entries; @@ -778,7 +763,7 @@ namespace GridKit } for (const auto& bus : buses_) { - auto bus_jacobian = bus->getJacobian(); + auto& bus_jacobian = bus->getJacobian(); std::tuple&, std::vector&, std::vector&> bus_jacobian_entries = bus_jacobian.getEntries(false); const auto [rows, columns, values] = bus_jacobian_entries; @@ -796,7 +781,7 @@ namespace GridKit IdxT counter = 0; for (const auto& component : components_) { - auto component_jacobian = component->getJacobian(); + auto& component_jacobian = component->getJacobian(); std::tuple&, std::vector&, std::vector&> component_jacobian_entries = component_jacobian.getEntries(false); const auto [rows, columns, values] = component_jacobian_entries; @@ -810,7 +795,7 @@ namespace GridKit } for (const auto& bus : buses_) { - auto bus_jacobian = bus->getJacobian(); + auto& bus_jacobian = bus->getJacobian(); std::tuple&, std::vector&, std::vector&> bus_jacobian_entries = bus_jacobian.getEntries(false); const auto [rows, columns, values] = bus_jacobian_entries; @@ -865,7 +850,7 @@ namespace GridKit IdxT counter = 0; for (const auto& component : components_) { - auto component_jacobian = component->getJacobian(); + auto& component_jacobian = component->getJacobian(); std::tuple&, std::vector&, std::vector&> component_jacobian_entries = component_jacobian.getEntries(false); const auto [rows, columns, values] = component_jacobian_entries; @@ -877,7 +862,7 @@ namespace GridKit } for (const auto& bus : buses_) { - auto bus_jacobian = bus->getJacobian(); + auto& bus_jacobian = bus->getJacobian(); std::tuple&, std::vector&, std::vector&> bus_jacobian_entries = bus_jacobian.getEntries(false); const auto [rows, columns, values] = bus_jacobian_entries; @@ -910,29 +895,6 @@ namespace GridKit monitor_.print(); } - /** - * @brief Update variables in buses and components - */ - void updateVariables() - { - for (const auto& bus : buses_) - { - for (IdxT j = 0; j < bus->size(); ++j) - { - bus->y()[j] = y_[bus->getVariableIndex(j)]; - bus->yp()[j] = yp_[bus->getVariableIndex(j)]; - } - } - for (const auto& component : components_) - { - for (IdxT j = 0; j < component->size(); ++j) - { - component->y()[j] = y_[component->getVariableIndex(j)]; - component->yp()[j] = yp_[component->getVariableIndex(j)]; - } - } - } - /** * @brief Update time * @@ -945,8 +907,6 @@ namespace GridKit { component->updateTime(t, a); } - - updateVariables(); } /** diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 28d8f3299..ca85877f3 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -9,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -151,6 +153,54 @@ namespace GridKit return success.report(__func__); } + TestOutcome elementStateAliasesSystemState() + { + TestStatus success = true; + + PhasorDynamics::SystemModel system; + + PhasorDynamics::Bus bus1(1.0, 2.0); + PhasorDynamics::Bus bus2(3.0, 4.0); + system.addBus(&bus1); + system.addBus(&bus2); + + PhasorDynamics::Branch branch(&bus1, &bus2, 0.0, 1.0, 0.0, 0.0); + system.addComponent(&branch); + + PhasorDynamics::Load load(&bus2, 1.0, 1.0); + system.addComponent(&load); + + system.allocate(); + system.initialize(); + + const auto bus2_vr_index = static_cast(bus2.getVariableIndex(0)); + const auto bus2_vi_index = static_cast(bus2.getVariableIndex(1)); + const auto load_ir_index = static_cast(load.getVariableIndex(0)); + const auto load_ii_index = static_cast(load.getVariableIndex(1)); + + system.y()[bus2_vr_index] = 5.0; + system.y()[bus2_vi_index] = 6.0; + success *= isEqual(bus2.Vr(), static_cast(5.0)); + success *= isEqual(bus2.Vi(), static_cast(6.0)); + + bus2.Vr() = 7.0; + bus2.Vi() = 8.0; + success *= isEqual(system.y()[bus2_vr_index], static_cast(7.0)); + success *= isEqual(system.y()[bus2_vi_index], static_cast(8.0)); + + system.y()[load_ir_index] = 9.0; + system.y()[load_ii_index] = 10.0; + success *= isEqual(load.yData()[0], static_cast(9.0)); + success *= isEqual(load.yData()[1], static_cast(10.0)); + + load.yData()[0] = 11.0; + load.yData()[1] = 12.0; + success *= isEqual(system.y()[load_ir_index], static_cast(11.0)); + success *= isEqual(system.y()[load_ii_index], static_cast(12.0)); + + return success.report(__func__); + } + /** * @brief Test for exception when signals are incorrectly configured */ diff --git a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp index a8b432be2..b337028aa 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp @@ -10,6 +10,7 @@ int main() result += test.constructor(); result += test.composer(); + result += test.elementStateAliasesSystemState(); #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian(); #endif From c1bbacd00f49e79821599e3e5ad041dc8c699da1 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 28 May 2026 17:54:44 -0500 Subject: [PATCH 2/3] Move to LinearAlgebra before further changes --- GridKit/LinearAlgebra/CMakeLists.txt | 1 + GridKit/LinearAlgebra/Vector.hpp | 110 +++++++++++++++++ GridKit/Model/PhasorDynamics/GridElement.hpp | 109 +++-------------- tests/UnitTests/LinearAlgebra/CMakeLists.txt | 1 + .../LinearAlgebra/Vector/CMakeLists.txt | 8 ++ .../LinearAlgebra/Vector/VectorTests.hpp | 115 ++++++++++++++++++ .../LinearAlgebra/Vector/runVectorTests.cpp | 25 ++++ .../UnitTests/PhasorDynamics/SystemTests.hpp | 8 ++ 8 files changed, 285 insertions(+), 92 deletions(-) create mode 100644 GridKit/LinearAlgebra/Vector.hpp create mode 100644 tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt create mode 100644 tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp create mode 100644 tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp diff --git a/GridKit/LinearAlgebra/CMakeLists.txt b/GridKit/LinearAlgebra/CMakeLists.txt index 3d324d485..c8a6609d7 100644 --- a/GridKit/LinearAlgebra/CMakeLists.txt +++ b/GridKit/LinearAlgebra/CMakeLists.txt @@ -5,5 +5,6 @@ add_subdirectory(DenseMatrix) install( FILES ${CMAKE_CURRENT_SOURCE_DIR}/MemoryUtils.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/Vector.hpp DESTINATION include/GridKit/LinearAlgebra) diff --git a/GridKit/LinearAlgebra/Vector.hpp b/GridKit/LinearAlgebra/Vector.hpp new file mode 100644 index 000000000..d428dc7c3 --- /dev/null +++ b/GridKit/LinearAlgebra/Vector.hpp @@ -0,0 +1,110 @@ +#pragma once + +#include +#include +#include +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + /** + * @brief Owns vector storage or aliases externally-owned contiguous data. + */ + template + class Vector + { + public: + using StorageT = std::vector; + using size_type = typename StorageT::size_type; + + Vector() = default; + + Vector(const Vector&) = delete; + Vector& operator=(const Vector&) = delete; + Vector(Vector&&) = delete; + Vector& operator=(Vector&&) = delete; + + void resize(size_type size) + { + if (!owns_data_) + { + if (size != size_) + { + throw std::runtime_error("Cannot resize externally-owned vector data"); + } + return; + } + + storage_.resize(size); + data_ = storage_.empty() ? nullptr : storage_.data(); + size_ = storage_.size(); + } + + void setDataPointer(ScalarT* data, size_type size) + { + if (data == nullptr && size != 0) + { + throw std::runtime_error("Cannot set nonzero vector data to nullptr"); + } + + data_ = data; + size_ = size; + owns_data_ = false; + storage_.clear(); + } + + ScalarT* data() + { + return data_; + } + + const ScalarT* data() const + { + return data_; + } + + size_type size() const + { + return size_; + } + + ScalarT& operator[](size_type i) + { + assert(i < size_); + return data_[i]; + } + + const ScalarT& operator[](size_type i) const + { + assert(i < size_); + return data_[i]; + } + + StorageT& storage() + { + if (!owns_data_) + { + throw std::runtime_error("Vector storage is externally owned"); + } + return storage_; + } + + const StorageT& storage() const + { + if (!owns_data_) + { + throw std::runtime_error("Vector storage is externally owned"); + } + return storage_; + } + + private: + StorageT storage_; + ScalarT* data_{nullptr}; + size_type size_{0}; + bool owns_data_{true}; + }; + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/GridElement.hpp b/GridKit/Model/PhasorDynamics/GridElement.hpp index 652fe01e7..2e8724d7d 100644 --- a/GridKit/Model/PhasorDynamics/GridElement.hpp +++ b/GridKit/Model/PhasorDynamics/GridElement.hpp @@ -1,10 +1,11 @@ #pragma once -#include +#include #include #include #include +#include #include #include #include @@ -15,89 +16,6 @@ namespace GridKit { using Log = ::GridKit::Utilities::Logger; - /// Owns state storage for standalone elements, or aliases a SystemModel-owned global slice. - /// The owned storage remains the std::vector returned through Evaluator; the data accessors - /// are the live state used by component implementations. - template - class StateVector - { - public: - using StorageT = std::vector; - using size_type = typename StorageT::size_type; - - StateVector() = default; - - // Non-copyable/movable: data_ aliases either storage_ or a global slice, and a GridElement - // (which owns these) is never copied or moved -- it is held by raw pointer in SystemModel. - StateVector(const StateVector&) = delete; - StateVector& operator=(const StateVector&) = delete; - StateVector(StateVector&&) = delete; - StateVector& operator=(StateVector&&) = delete; - - template - void resize(SizeT size) - { - const auto sz = static_cast(size); - if (viewing_) - { - assert(sz == size_); - return; - } - - storage_.resize(sz); - data_ = storage_.data(); - size_ = storage_.size(); - } - - void view(ScalarT* data, size_type size) - { - data_ = data; - size_ = size; - viewing_ = true; - } - - ScalarT* data() - { - return data_; - } - - const ScalarT* data() const - { - return data_; - } - - size_type size() const - { - return size_; - } - - ScalarT& operator[](size_type i) - { - return data_[i]; - } - - const ScalarT& operator[](size_type i) const - { - return data_[i]; - } - - StorageT& storage() - { - return storage_; - } - - const StorageT& storage() const - { - return storage_; - } - - private: - StorageT storage_; - ScalarT* data_{nullptr}; - size_type size_{0}; - bool viewing_{false}; - }; - /** * @brief Model base class for all system constituents */ @@ -107,6 +25,7 @@ namespace GridKit public: using RealT = typename Model::Evaluator::RealT; using MatrixT = typename Model::Evaluator::MatrixT; + using VectorT = LinearAlgebra::Vector; GridElement() : size_{0} @@ -149,8 +68,6 @@ namespace GridKit msa = max_steps_; } - // Evaluator requires std::vector references. These return owned storage; when an element is - // bound into a SystemModel, live state is accessed through y_/yp_/f_ instead. std::vector& y() override { return y_.storage(); @@ -205,9 +122,17 @@ namespace GridKit /// element reads and writes the global state in place (no scatter/gather). void bindGlobalState(ScalarT* y, ScalarT* yp, ScalarT* f, std::size_t n) { - y_.view(n == 0 ? nullptr : y, n); - yp_.view(n == 0 ? nullptr : yp, n); - f_.view(n == 0 ? nullptr : f, n); + if (n == 0) + { + y_.setDataPointer(nullptr, 0); + yp_.setDataPointer(nullptr, 0); + f_.setDataPointer(nullptr, 0); + return; + } + + y_.setDataPointer(y, n); + yp_.setDataPointer(yp, n); + f_.setDataPointer(f, n); } MatrixT& getJacobian() override @@ -267,9 +192,9 @@ namespace GridKit /// Global (system-level) residual indices std::vector residual_indices_; - StateVector y_; - StateVector yp_; - StateVector f_; + VectorT y_; + VectorT yp_; + VectorT f_; std::vector tag_; std::vector g_; diff --git a/tests/UnitTests/LinearAlgebra/CMakeLists.txt b/tests/UnitTests/LinearAlgebra/CMakeLists.txt index 58e28ddcb..801a7b3b7 100644 --- a/tests/UnitTests/LinearAlgebra/CMakeLists.txt +++ b/tests/UnitTests/LinearAlgebra/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(SparseMatrix) +add_subdirectory(Vector) diff --git a/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt b/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt new file mode 100644 index 000000000..3fa1b6980 --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt @@ -0,0 +1,8 @@ +add_executable(test_linear_algebra_vector runVectorTests.cpp) +target_link_libraries(test_linear_algebra_vector + GridKit::testing) + +add_test(NAME VectorTest COMMAND $) + +install(TARGETS test_linear_algebra_vector + RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp new file mode 100644 index 000000000..0885e972e --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp @@ -0,0 +1,115 @@ +#pragma once + +#include + +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + using namespace LinearAlgebra; + + template + class VectorTests + { + public: + TestOutcome storage() + { + TestStatus success = true; + + Vector v; + + v.resize(3); + success *= v.size() == 3; + success *= v.data() == v.storage().data(); + + v[0] = static_cast(1.0); + v[1] = static_cast(2.0); + v[2] = static_cast(3.0); + + success *= isEqual(v.storage()[0], static_cast(1.0)); + success *= isEqual(v.storage()[1], static_cast(2.0)); + success *= isEqual(v.storage()[2], static_cast(3.0)); + + return success.report(__func__); + } + + TestOutcome externalView() + { + TestStatus success = true; + + ScalarT data[3] = {static_cast(1.0), + static_cast(2.0), + static_cast(3.0)}; + Vector v; + + v.setDataPointer(data, 3); + success *= v.size() == 3; + success *= v.data() == data; + + v[1] = static_cast(4.0); + success *= isEqual(data[1], static_cast(4.0)); + + data[2] = static_cast(5.0); + success *= isEqual(v[2], static_cast(5.0)); + + success *= throws( + [&]() + { v.storage(); }); + + return success.report(__func__); + } + + TestOutcome externalResize() + { + TestStatus success = true; + + ScalarT data[3] = {static_cast(1.0), + static_cast(2.0), + static_cast(3.0)}; + Vector v; + + v.setDataPointer(data, 3); + v.resize(3); + success *= v.data() == data; + success *= v.size() == 3; + + success *= throws( + [&]() + { v.resize(4); }); + success *= v.data() == data; + success *= v.size() == 3; + + return success.report(__func__); + } + + TestOutcome dependencyTrackingVariable() + { + TestStatus success = true; + + using Variable = DependencyTracking::Variable; + + Vector v; + + v.resize(2); + v[0] = Variable(1.0); + v[1] = Variable(2.0); + + success *= isEqual(v[0].getValue(), 1.0); + success *= isEqual(v[1].getValue(), 2.0); + + Variable external[2] = {Variable(3.0), Variable(4.0)}; + v.setDataPointer(external, 2); + v[0] = Variable(5.0); + + success *= isEqual(external[0].getValue(), 5.0); + success *= isEqual(v[1].getValue(), 4.0); + + return success.report(__func__); + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp b/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp new file mode 100644 index 000000000..edc679a4b --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp @@ -0,0 +1,25 @@ +#include "VectorTests.hpp" + +using namespace GridKit; +using namespace LinearAlgebra; +using namespace Testing; + +template +void runTests(TestingResults& result) +{ + VectorTests test; + + result += test.storage(); + result += test.externalView(); + result += test.externalResize(); + result += test.dependencyTrackingVariable(); +} + +int main(int, char**) +{ + TestingResults result; + + runTests(result); + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index ca85877f3..84e2db9f6 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -198,6 +199,13 @@ namespace GridKit success *= isEqual(system.y()[load_ir_index], static_cast(11.0)); success *= isEqual(system.y()[load_ii_index], static_cast(12.0)); + success *= throws( + [&]() + { bus2.y(); }); + success *= throws( + [&]() + { load.getResidual(); }); + return success.report(__func__); } From 52ec005b82cb30eca3035fc8b030da02e2962a72 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 28 May 2026 18:13:49 -0500 Subject: [PATCH 3/3] add comment about future work --- GridKit/LinearAlgebra/Vector.hpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/GridKit/LinearAlgebra/Vector.hpp b/GridKit/LinearAlgebra/Vector.hpp index d428dc7c3..b139fc623 100644 --- a/GridKit/LinearAlgebra/Vector.hpp +++ b/GridKit/LinearAlgebra/Vector.hpp @@ -11,6 +11,15 @@ namespace GridKit { /** * @brief Owns vector storage or aliases externally-owned contiguous data. + * + * Templated on ScalarT to hold AD types, so it cannot reuse Re::Solve's + * Vector directly. storage() exists only to satisfy + * the std::vector&-based Model::Evaluator interface. + * + * Future work: change Evaluator's state accessors (y/yp/getResidual/tag) to + * return Vector&, which removes storage() and lets SUNDIALS use + * data()/size() directly. Finalize this public API before it becomes an + * interface type. */ template class Vector