From 3ea9b97822729d678a2557812f0df8a70c10528c Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sun, 21 Jun 2026 11:36:04 -0500 Subject: [PATCH 1/2] Migrate Evaluator state to View; auto-bind in allocate --- GridKit/LinearAlgebra/MemoryUtils.hpp | 10 +- GridKit/LinearAlgebra/Vector/Vector.cpp | 65 ++++++- GridKit/LinearAlgebra/Vector/Vector.hpp | 42 +++- GridKit/Model/Evaluator.hpp | 174 ++++++++++++++--- .../PhasorDynamics/Branch/BranchImpl.hpp | 15 ++ GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp | 27 +-- .../PhasorDynamics/Bus/BusInfiniteImpl.hpp | 17 +- GridKit/Model/PhasorDynamics/BusBase.hpp | 83 ++++---- .../PhasorDynamics/BusFault/BusFaultImpl.hpp | 29 +-- .../BusToSignalAdapterImpl.hpp | 15 ++ GridKit/Model/PhasorDynamics/CMakeLists.txt | 2 + GridKit/Model/PhasorDynamics/Component.hpp | 83 ++++---- .../Exciter/IEEET1/Ieeet1Impl.hpp | 31 +-- .../Exciter/SEXS-PTI/SexsPtiImpl.hpp | 29 +-- .../Governor/Tgov1/Tgov1Impl.hpp | 30 +-- .../PhasorDynamics/Load/LoadZ/LoadZImpl.hpp | 32 ++-- .../Load/LoadZIP/LoadZIPImpl.hpp | 32 ++-- .../Stabilizer/IEEEST/IeeestImpl.hpp | 29 +-- .../GENROUwS/GenrouImpl.hpp | 31 +-- .../GENSALwS/GensalImpl.hpp | 31 +-- .../GenClassical/GenClassicalImpl.hpp | 31 +-- GridKit/Model/PhasorDynamics/SystemModel.hpp | 9 +- .../Model/PhasorDynamics/SystemModelImpl.hpp | 181 ++++++------------ GridKit/Model/PowerElectronics/CMakeLists.txt | 3 + .../PowerElectronics/Capacitor/CMakeLists.txt | 3 +- .../PowerElectronics/CircuitComponent.hpp | 95 ++++----- .../Model/PowerElectronics/CircuitNode.hpp | 98 +++++----- .../DistributedGenerator/CMakeLists.txt | 3 +- .../DistributedGenerator.cpp | 2 +- .../InductionMotor/CMakeLists.txt | 3 +- .../PowerElectronics/Inductor/CMakeLists.txt | 3 +- .../PowerElectronics/Inductor/Inductor.cpp | 2 +- .../LinearTransformer/CMakeLists.txt | 3 +- .../MicrogridBusDQ/CMakeLists.txt | 3 +- .../MicrogridLine/CMakeLists.txt | 3 +- .../MicrogridLine/MicrogridLine.cpp | 2 +- .../MicrogridLoad/CMakeLists.txt | 3 +- .../MicrogridLoad/MicrogridLoad.cpp | 2 +- GridKit/Model/PowerElectronics/NodeBase.hpp | 100 +++++----- .../PowerElectronics/Resistor/CMakeLists.txt | 3 +- .../PowerElectronics/Resistor/Resistor.cpp | 2 +- .../SynchronousMachine/CMakeLists.txt | 3 +- .../SystemModelPowerElectronics.hpp | 65 +++++-- .../TransmissionLine/CMakeLists.txt | 3 +- .../VoltageSource/CMakeLists.txt | 3 +- .../VoltageSource/VoltageSource.cpp | 2 +- GridKit/Model/PowerFlow/Branch/CMakeLists.txt | 3 +- GridKit/Model/PowerFlow/Bus/BusPQ.cpp | 24 ++- GridKit/Model/PowerFlow/Bus/BusPV.cpp | 24 ++- GridKit/Model/PowerFlow/Bus/CMakeLists.txt | 3 +- .../Model/PowerFlow/Generator/CMakeLists.txt | 3 +- .../Model/PowerFlow/Generator2/CMakeLists.txt | 3 +- .../Model/PowerFlow/Generator2/Generator2.cpp | 4 +- .../Model/PowerFlow/Generator4/CMakeLists.txt | 3 +- .../Model/PowerFlow/Generator4/Generator4.cpp | 5 +- .../Generator4Governor/CMakeLists.txt | 3 +- .../Generator4Governor/Generator4Governor.cpp | 5 +- .../PowerFlow/Generator4Param/CMakeLists.txt | 3 +- .../Generator4Param/Generator4Param.cpp | 5 +- GridKit/Model/PowerFlow/Load/CMakeLists.txt | 3 +- .../Model/PowerFlow/MiniGrid/CMakeLists.txt | 3 +- .../Model/PowerFlow/ModelEvaluatorImpl.hpp | 133 ++++++++----- GridKit/Model/PowerFlow/SystemModel.hpp | 49 +++-- .../Model/PowerFlow/SystemModelPowerFlow.hpp | 7 +- GridKit/Solver/Dynamic/Ida.cpp | 49 ++--- GridKit/Solver/Dynamic/Ida.hpp | 8 +- GridKit/Solver/SteadyState/Kinsol.cpp | 11 +- GridKit/Solver/SteadyState/Kinsol.hpp | 8 +- .../TenGen/Classical/TenGenClassical.cpp | 2 +- .../Small/TenGen/Genrou/TenGenGenrou.cpp | 2 +- .../Tiny/ThreeBus/Basic/ThreeBusBasic.cpp | 2 +- .../ThreeBus/Classical/ThreeBusClassical.cpp | 2 +- .../Classical/ThreeBusClassicalJson.cpp | 2 +- .../ThreeBus/ZipLoad/ThreeBusZipLoadJson.cpp | 2 +- .../Tiny/TwoBus/Basic/TwoBusBasic.cpp | 2 +- .../Tiny/TwoBus/Basic/TwoBusBasicJson.cpp | 2 +- .../Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp | 2 +- .../Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp | 2 +- .../Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp | 2 +- .../Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp | 2 +- .../DistributedGeneratorTest/DGTest.cpp | 7 +- .../PowerElectronics/Microgrid/Microgrid.cpp | 4 +- .../PowerElectronics/RLCircuit/RLCircuit.cpp | 7 +- .../ScaleMicrogrid/ScaleMicrogrid.cpp | 4 +- .../LinearAlgebra/Vector/VectorTests.hpp | 46 ++++- .../LinearAlgebra/Vector/runVectorTests.cpp | 1 + .../PhasorDynamics/BusFaultTests.hpp | 12 +- .../PhasorDynamics/ExciterIeeet1Tests.hpp | 6 +- .../PhasorDynamics/ExciterSexsPtiTests.hpp | 6 +- .../PhasorDynamics/GenClassicalTests.hpp | 10 +- .../UnitTests/PhasorDynamics/GenrouTests.hpp | 14 +- .../UnitTests/PhasorDynamics/GensalTests.hpp | 20 +- .../PhasorDynamics/GovernorTgov1Tests.hpp | 14 +- .../UnitTests/PhasorDynamics/LoadZIPTests.hpp | 12 +- tests/UnitTests/PhasorDynamics/LoadZTests.hpp | 2 +- .../PhasorDynamics/StabilizerIeeestTests.hpp | 6 +- .../UnitTests/PhasorDynamics/SystemTests.hpp | 66 ++++++- .../PhasorDynamics/runSystemTests.cpp | 1 + tests/UnitTests/Solver/Dynamic/CMakeLists.txt | 5 +- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 143 ++++++-------- 100 files changed, 1306 insertions(+), 900 deletions(-) diff --git a/GridKit/LinearAlgebra/MemoryUtils.hpp b/GridKit/LinearAlgebra/MemoryUtils.hpp index 9738b890a..bcf12fe53 100644 --- a/GridKit/LinearAlgebra/MemoryUtils.hpp +++ b/GridKit/LinearAlgebra/MemoryUtils.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include namespace GridKit @@ -206,8 +207,7 @@ namespace GridKit template int allocateArrayOnHost(T** v, I n) { - std::size_t arraysize = static_cast(n) * sizeof(T); - *v = new T[arraysize]; + *v = new T[static_cast(n)](); return *v == nullptr ? 1 : 0; } @@ -222,16 +222,14 @@ namespace GridKit template int copyArrayHostToHost(T* dst, const T* src, I n) { - std::size_t arraysize = static_cast(n) * sizeof(T); - memcpy(dst, src, arraysize); + std::copy_n(src, static_cast(n), dst); return 0; } template int setZeroArrayOnHost(T* v, I n) { - std::size_t arraysize = static_cast(n) * sizeof(T); - memset(v, 0, arraysize); + std::fill_n(v, static_cast(n), T{}); return 0; } diff --git a/GridKit/LinearAlgebra/Vector/Vector.cpp b/GridKit/LinearAlgebra/Vector/Vector.cpp index abb77228e..2f5cbb0c2 100644 --- a/GridKit/LinearAlgebra/Vector/Vector.cpp +++ b/GridKit/LinearAlgebra/Vector/Vector.cpp @@ -1,6 +1,7 @@ #include -#include +#include +#include #include #include @@ -39,8 +40,8 @@ namespace GridKit : n_capacity_(n), k_(k), n_size_(n), - gpu_updated_(new bool[k]), - cpu_updated_(new bool[k]) + gpu_updated_(new bool[static_cast(k)]), + cpu_updated_(new bool[static_cast(k)]) { setHostUpdated(false); setDeviceUpdated(false); @@ -52,13 +53,64 @@ namespace GridKit */ template Vector::~Vector() + { + release(); + } + + template + Vector::Vector(Vector&& other) noexcept + { + *this = std::move(other); + } + + template + Vector& Vector::operator=(Vector&& other) noexcept + { + if (this != &other) + { + release(); + + n_capacity_ = other.n_capacity_; + k_ = other.k_; + n_size_ = other.n_size_; + d_data_ = other.d_data_; + h_data_ = other.h_data_; + gpu_updated_ = other.gpu_updated_; + cpu_updated_ = other.cpu_updated_; + owns_gpu_data_ = other.owns_gpu_data_; + owns_cpu_data_ = other.owns_cpu_data_; + mem_ = std::move(other.mem_); + + other.n_capacity_ = 0; + other.k_ = 0; + other.n_size_ = 0; + other.d_data_ = nullptr; + other.h_data_ = nullptr; + other.gpu_updated_ = nullptr; + other.cpu_updated_ = nullptr; + } + + return *this; + } + + template + void Vector::release() { if (owns_cpu_data_ && h_data_) + { mem_.deleteOnHost(h_data_); + } if (owns_gpu_data_ && d_data_) + { mem_.deleteOnDevice(d_data_); + } delete[] gpu_updated_; delete[] cpu_updated_; + + h_data_ = nullptr; + d_data_ = nullptr; + gpu_updated_ = nullptr; + cpu_updated_ = nullptr; } /** @@ -999,9 +1051,12 @@ namespace GridKit std::fill(gpu_updated_, gpu_updated_ + k_, is_updated); } - // template class Vector; + template class Vector; template class Vector; - // template class Vector; + template class Vector; + template class Vector; + template class Vector; + template class Vector; } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/Vector.hpp b/GridKit/LinearAlgebra/Vector/Vector.hpp index 2690fb0d9..017d849e0 100644 --- a/GridKit/LinearAlgebra/Vector/Vector.hpp +++ b/GridKit/LinearAlgebra/Vector/Vector.hpp @@ -1,4 +1,6 @@ #pragma once +#include +#include #include #include @@ -36,10 +38,10 @@ namespace GridKit Vector(IdxT n, IdxT k); ~Vector(); - Vector(const Vector&) = delete; - Vector(Vector&&) = delete; + Vector(const Vector&) = delete; + Vector(Vector&& other) noexcept; Vector& operator=(const Vector&) = delete; - Vector& operator=(Vector&&) = delete; + Vector& operator=(Vector&& other) noexcept; int copyFromExternal(const ScalarT* source, memory::MemorySpace memspaceIn = memory::HOST, @@ -53,6 +55,38 @@ namespace GridKit const ScalarT* getData(memory::MemorySpace memspace = memory::HOST) const; const ScalarT* getData(IdxT i, memory::MemorySpace memspace = memory::HOST) const; + std::size_t size() const + { + return static_cast(getSize()); + } + + bool empty() const + { + return size() == 0; + } + + ScalarT* data(memory::MemorySpace memspace = memory::HOST) + { + return getData(memspace); + } + + const ScalarT* data(memory::MemorySpace memspace = memory::HOST) const + { + return getData(memspace); + } + + ScalarT& operator[](std::size_t i) + { + assert(i < size()); + return data()[i]; + } + + const ScalarT& operator[](std::size_t i) const + { + assert(i < size()); + return data()[i]; + } + IdxT getCapacity() const; IdxT getSize() const; IdxT getNumVectors() const; @@ -77,6 +111,7 @@ namespace GridKit memory::MemorySpace memspaceDst = memory::HOST); private: + void release(); void setHostUpdated(bool is_updated); void setDeviceUpdated(bool is_updated); @@ -93,5 +128,6 @@ namespace GridKit MemoryHandler mem_; ///< Device memory manager object }; + } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index 2600f52af..e26240a78 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -1,10 +1,11 @@ #pragma once -#include +#include #include #include #include +#include #include #include @@ -12,6 +13,8 @@ namespace GridKit { namespace Model { + namespace memory = GridKit::LinearAlgebra::memory; + /*! * @brief Abstract class describing a model. * @@ -25,6 +28,7 @@ namespace GridKit using RealT = typename GridKit::ScalarTraits::RealT; using CsrMatrixT = GridKit::LinearAlgebra::CsrMatrix; using CooMatrixT = GridKit::LinearAlgebra::CooMatrix; + using VectorT = GridKit::LinearAlgebra::Vector; Evaluator() { @@ -122,55 +126,167 @@ namespace GridKit /** * @brief Get the absolute tolerance for each variable in the model * - * @return a reference to the absolute tolerance vector. + * @return a view of the absolute tolerance vector. * * @pre `setAbsoluteTolerance` must have been called first. */ - virtual std::vector& absoluteTolerance() = 0; + virtual VectorT& absoluteTolerance() + { + return abs_tol_; + } + /** * @brief Get the absolute tolerance for each variable in the model * - * @return a const reference to the absolute tolerance vector. + * @return a const view of the absolute tolerance vector. * * @pre `setAbsoluteTolerance` must have been called first. */ - virtual const std::vector& absoluteTolerance() const = 0; + virtual const VectorT& absoluteTolerance() const + { + return abs_tol_; + } - virtual std::vector& y() = 0; - virtual const std::vector& y() const = 0; + virtual VectorT& y() + { + return y_; + } - virtual std::vector& yp() = 0; - virtual const std::vector& yp() const = 0; + virtual const VectorT& y() const + { + return y_; + } - virtual std::vector& tag() = 0; - virtual const std::vector& tag() const = 0; + virtual VectorT& yp() + { + return yp_; + } - virtual std::vector& yB() = 0; - virtual const std::vector& yB() const = 0; + virtual const VectorT& yp() const + { + return yp_; + } - virtual std::vector& ypB() = 0; - virtual const std::vector& ypB() const = 0; + /** + * @brief Get the differential/algebraic tag for each variable. + * + * Entries are scalar values for direct transfer to solver vectors. + */ + virtual VectorT& tag() + { + return tag_; + } - virtual std::vector& param() = 0; - virtual const std::vector& param() const = 0; + /** + * @brief Get the differential/algebraic tag for each variable. + */ + virtual const VectorT& tag() const + { + return tag_; + } + + virtual VectorT& yB() = 0; + virtual const VectorT& yB() const = 0; - virtual std::vector& param_up() = 0; - virtual const std::vector& param_up() const = 0; + virtual VectorT& ypB() = 0; + virtual const VectorT& ypB() const = 0; - virtual std::vector& param_lo() = 0; - virtual const std::vector& param_lo() const = 0; + virtual VectorT& param() = 0; + virtual const VectorT& param() const = 0; - virtual std::vector& getResidual() = 0; - virtual const std::vector& getResidual() const = 0; + virtual VectorT& param_up() = 0; + virtual const VectorT& param_up() const = 0; - virtual std::vector& getIntegrand() = 0; - virtual const std::vector& getIntegrand() const = 0; + virtual VectorT& param_lo() = 0; + virtual const VectorT& param_lo() const = 0; + + virtual VectorT& getResidual() + { + return f_; + } - virtual std::vector& getAdjointResidual() = 0; - virtual const std::vector& getAdjointResidual() const = 0; + virtual const VectorT& getResidual() const + { + return f_; + } + + virtual VectorT& getIntegrand() = 0; + virtual const VectorT& getIntegrand() const = 0; + + virtual VectorT& getAdjointResidual() = 0; + virtual const VectorT& getAdjointResidual() const = 0; + + virtual VectorT& getAdjointIntegrand() = 0; + virtual const VectorT& getAdjointIntegrand() const = 0; + + /** + * @brief Bind this evaluator's state and residual vectors. + */ + virtual void bind(VectorT& y, VectorT& yp, VectorT& f, VectorT& tag, VectorT& abs_tol, IdxT offset) + { + const IdxT n = size(); + assert(static_cast(offset) + static_cast(n) <= y.size()); + assert(static_cast(offset) + static_cast(n) <= yp.size()); + assert(static_cast(offset) + static_cast(n) <= f.size()); + assert(static_cast(offset) + static_cast(n) <= tag.size()); + assert(static_cast(offset) + static_cast(n) <= abs_tol.size()); + + y_ = VectorT(n); + y_.setData(y.data() + offset); + + yp_ = VectorT(n); + yp_.setData(yp.data() + offset); + + f_ = VectorT(n); + f_.setData(f.data() + offset); + + tag_ = VectorT(n); + tag_.setData(tag.data() + offset); + + abs_tol_ = VectorT(n); + abs_tol_.setData(abs_tol.data() + offset); + + offset_ = offset; + allocated_ = true; + } + + protected: + /** + * @brief Allocate this evaluator's state and residual vectors. + */ + void allocateVectors(IdxT n) + { + y_ = VectorT(n); + y_.allocate(memory::HOST); + y_.setDataUpdated(memory::HOST); + + yp_ = VectorT(n); + yp_.allocate(memory::HOST); + yp_.setDataUpdated(memory::HOST); + + f_ = VectorT(n); + f_.allocate(memory::HOST); + f_.setDataUpdated(memory::HOST); + + tag_ = VectorT(n); + tag_.allocate(memory::HOST); + tag_.setDataUpdated(memory::HOST); + + abs_tol_ = VectorT(n); + abs_tol_.allocate(memory::HOST); + abs_tol_.setDataUpdated(memory::HOST); + + offset_ = 0; + allocated_ = true; + } - virtual std::vector& getAdjointIntegrand() = 0; - virtual const std::vector& getAdjointIntegrand() const = 0; + VectorT y_; + VectorT yp_; + VectorT f_; + VectorT tag_; + VectorT abs_tol_; + IdxT offset_{0}; + bool allocated_{false}; }; } // namespace Model diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp index 3beccd4b4..a3680b5c7 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp @@ -119,6 +119,21 @@ namespace GridKit template int Branch::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } + auto size = static_cast(size_); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(this->tag_.size() == size); + assert(this->abs_tol_.size() == size); + + variable_indices_.resize(size); + residual_indices_.resize(size); + wb_.resize(2); h_.resize(2); diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 122102c4b..1a7b5c1af 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -101,28 +101,29 @@ namespace GridKit } /*! - * @brief allocate method resizes local solution and residual vectors. + * @brief Allocate bus storage and index maps. */ template int Bus::allocate() { - // Temporary while we use std::vector in the code + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } size_t size = static_cast(size_); - // Resize component model data - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); - - // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; } return 0; @@ -164,7 +165,7 @@ namespace GridKit template int Bus::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp index e5b0cfa7f..f3bbfcf45 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp @@ -99,11 +99,26 @@ namespace GridKit } /*! - * @brief allocate method resizes local solution and residual vectors. + * @brief Allocate bus storage and index maps. */ template int BusInfinite::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } + auto size = static_cast(size_); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(this->tag_.size() == size); + assert(this->abs_tol_.size() == size); + + variable_indices_.resize(size); + residual_indices_.resize(size); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 3ecfecec2..bcee3ae8e 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -1,5 +1,7 @@ #pragma once +#include +#include #include #include #include @@ -34,6 +36,7 @@ namespace GridKit using CooMatrixT = typename Model::Evaluator::CooMatrixT; using BusTypeT = typename BusData::BusType; using MonitorT = Model::VariableMonitor; + using VectorT = typename Model::Evaluator::VectorT; BusBase() = default; @@ -54,52 +57,32 @@ namespace GridKit return nnz_; } - std::vector& y() override - { - return y_; - } - - const std::vector& y() const override - { - return y_; - } - - std::vector& yp() override - { - return yp_; - } - - const std::vector& yp() const override - { - return yp_; - } - - std::vector& tag() override + VectorT& tag() override { return tag_; } - const std::vector& tag() const override + const VectorT& tag() const override { return tag_; } - std::vector& absoluteTolerance() override + VectorT& absoluteTolerance() override { return abs_tol_; } - const std::vector& absoluteTolerance() const override + const VectorT& absoluteTolerance() const override { return abs_tol_; } - std::vector& getResidual() override + VectorT& getResidual() override { return f_; } - const std::vector& getResidual() const override + const VectorT& getResidual() const override { return f_; } @@ -136,6 +119,17 @@ namespace GridKit return residual_indices_; } + protected: + using Model::Evaluator::y_; + using Model::Evaluator::yp_; + using Model::Evaluator::f_; + using Model::Evaluator::tag_; + using Model::Evaluator::abs_tol_; + using Model::Evaluator::offset_; + using Model::Evaluator::allocated_; + using Model::Evaluator::allocateVectors; + + public: /// Pure virtual function, returns bus type (DEFAULT or SLACK). virtual BusTypeT BusType() const { @@ -188,11 +182,6 @@ namespace GridKit /// Global (system-level) residual indices std::vector residual_indices_; - std::vector y_; - std::vector yp_; - std::vector tag_; - std::vector abs_tol_; - std::vector f_; std::vector g_; CsrMatrixT* csr_jac_{nullptr}; @@ -231,52 +220,52 @@ namespace GridKit throw NotImplementedError(__func__); } - [[noreturn]] std::vector& yB() override + [[noreturn]] VectorT& yB() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& yB() const override + [[noreturn]] const VectorT& yB() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& ypB() override + [[noreturn]] VectorT& ypB() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& ypB() const override + [[noreturn]] const VectorT& ypB() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& param() override + [[noreturn]] VectorT& param() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& param() const override + [[noreturn]] const VectorT& param() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& param_up() override + [[noreturn]] VectorT& param_up() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& param_up() const override + [[noreturn]] const VectorT& param_up() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& param_lo() override + [[noreturn]] VectorT& param_lo() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& param_lo() const override + [[noreturn]] const VectorT& param_lo() const override { throw NotImplementedError(__func__); } @@ -301,32 +290,32 @@ namespace GridKit throw NotImplementedError(__func__); } - [[noreturn]] std::vector& getIntegrand() override + [[noreturn]] VectorT& getIntegrand() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& getIntegrand() const override + [[noreturn]] const VectorT& getIntegrand() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& getAdjointResidual() override + [[noreturn]] VectorT& getAdjointResidual() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& getAdjointResidual() const override + [[noreturn]] const VectorT& getAdjointResidual() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& getAdjointIntegrand() override + [[noreturn]] VectorT& getAdjointIntegrand() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& getAdjointIntegrand() const override + [[noreturn]] const VectorT& getAdjointIntegrand() const override { throw NotImplementedError(__func__); } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 72b45850a..0a898224f 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -110,26 +110,31 @@ namespace GridKit template int BusFault::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } // std::cout << "Allocate BusFault..." << std::endl; + auto size = static_cast(size_); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(this->f_.size() == size); + assert(tag_.size() == size); + assert(this->abs_tol_.size() == size); - auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); - abs_tol_.resize(size); - tag_.resize(size); - variable_indices_.resize(size); - residual_indices_.resize(size); + this->variable_indices_.resize(size); + this->residual_indices_.resize(size); // Resize coupling data wb_.resize(2); h_.resize(2); - // Default variable and residual index mapping to local index + // Default variable and residual index mapping to global indices for (IdxT j = 0; j < size_; ++j) { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); + this->setVariableIndex(j, this->offset_ + j); + this->setResidualIndex(j, this->offset_ + j); } return 0; @@ -190,7 +195,7 @@ namespace GridKit template int BusFault::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterImpl.hpp b/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterImpl.hpp index 4ae676ccc..5158ba412 100644 --- a/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterImpl.hpp @@ -52,6 +52,21 @@ namespace GridKit template int BusToSignalAdapter::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } + auto size = static_cast(size_); + + assert(this->y_.size() == size); + assert(this->yp_.size() == size); + assert(this->f_.size() == size); + assert(this->tag_.size() == size); + assert(this->abs_tol_.size() == size); + + this->variable_indices_.resize(size); + this->residual_indices_.resize(size); + static constexpr auto VREAL = BusToSignalAdapterInternalVariables::VREAL; static constexpr auto VIMAG = BusToSignalAdapterInternalVariables::VIMAG; diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index baf76ab90..0782fd2fc 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -20,6 +20,8 @@ gridkit_add_library( INTERFACE GridKit::utilities_logger INTERFACE + GridKit::dense_vector + INTERFACE GridKit::sparse_matrix) target_link_libraries( diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 20dbd0403..c44a791db 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -1,5 +1,7 @@ #pragma once +#include +#include #include #include @@ -26,6 +28,7 @@ namespace GridKit using RealT = typename Model::Evaluator::RealT; using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; using CooMatrixT = typename Model::Evaluator::CooMatrixT; + using VectorT = typename Model::Evaluator::VectorT; Component() = default; @@ -72,52 +75,32 @@ namespace GridKit return nnz_; } - std::vector& y() override - { - return y_; - } - - const std::vector& y() const override - { - return y_; - } - - std::vector& yp() override - { - return yp_; - } - - const std::vector& yp() const override - { - return yp_; - } - - std::vector& tag() override + VectorT& tag() override { return tag_; } - const std::vector& tag() const override + const VectorT& tag() const override { return tag_; } - std::vector& absoluteTolerance() override + VectorT& absoluteTolerance() override { return abs_tol_; } - const std::vector& absoluteTolerance() const override + const VectorT& absoluteTolerance() const override { return abs_tol_; } - std::vector& getResidual() override + VectorT& getResidual() override { return f_; } - const std::vector& getResidual() const override + const VectorT& getResidual() const override { return f_; } @@ -164,6 +147,17 @@ namespace GridKit return coo_jac_; } + protected: + using Model::Evaluator::y_; + using Model::Evaluator::yp_; + using Model::Evaluator::f_; + using Model::Evaluator::tag_; + using Model::Evaluator::abs_tol_; + using Model::Evaluator::offset_; + using Model::Evaluator::allocated_; + using Model::Evaluator::allocateVectors; + + public: /// @todo Remove this method. It should be part of DynamicSolver class. bool hasJacobian() override { @@ -252,11 +246,6 @@ namespace GridKit /// Global (system-level) residual indices std::vector residual_indices_; - std::vector y_; - std::vector yp_; - std::vector tag_; - std::vector abs_tol_; - std::vector f_; std::vector g_; IdxT* J_rows_buffer_{nullptr}; @@ -305,52 +294,52 @@ namespace GridKit throw NotImplementedError(__func__); } - [[noreturn]] std::vector& yB() override + [[noreturn]] VectorT& yB() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& yB() const override + [[noreturn]] const VectorT& yB() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& ypB() override + [[noreturn]] VectorT& ypB() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& ypB() const override + [[noreturn]] const VectorT& ypB() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& param() override + [[noreturn]] VectorT& param() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& param() const override + [[noreturn]] const VectorT& param() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& param_up() override + [[noreturn]] VectorT& param_up() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& param_up() const override + [[noreturn]] const VectorT& param_up() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& param_lo() override + [[noreturn]] VectorT& param_lo() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& param_lo() const override + [[noreturn]] const VectorT& param_lo() const override { throw NotImplementedError(__func__); } @@ -375,32 +364,32 @@ namespace GridKit throw NotImplementedError(__func__); } - [[noreturn]] std::vector& getIntegrand() override + [[noreturn]] VectorT& getIntegrand() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& getIntegrand() const override + [[noreturn]] const VectorT& getIntegrand() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& getAdjointResidual() override + [[noreturn]] VectorT& getAdjointResidual() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& getAdjointResidual() const override + [[noreturn]] const VectorT& getAdjointResidual() const override { throw NotImplementedError(__func__); } - [[noreturn]] std::vector& getAdjointIntegrand() override + [[noreturn]] VectorT& getAdjointIntegrand() override { throw NotImplementedError(__func__); } - [[noreturn]] const std::vector& getAdjointIntegrand() const override + [[noreturn]] const VectorT& getAdjointIntegrand() const override { throw NotImplementedError(__func__); } diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 87b824dd1..4856b13f4 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -109,15 +109,25 @@ namespace GridKit template int Ieeet1::allocate() { - // Resize component model data + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } // Resize bus data wb_.resize(2); @@ -130,13 +140,6 @@ namespace GridKit ws_[1] = 0.0; ws_indices_[1] = INVALID_INDEX; - // Default variable and residual index mapping to local index - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - // Set output signals if (signals_.template isAssigned()) { @@ -290,7 +293,7 @@ namespace GridKit template int Ieeet1::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index cbbff6399..a3e7d1c2b 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -57,14 +57,25 @@ namespace GridKit template int SexsPti::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } wb_.resize(2); @@ -73,12 +84,6 @@ namespace GridKit ws_[0] = 0.0; ws_indices_[0] = INVALID_INDEX; - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - if (signals_.template isAssigned()) { signals_.template getSignalNode()->set( @@ -196,7 +201,7 @@ namespace GridKit template int SexsPti::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 2d9280fc9..166cd4c2d 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -164,15 +164,26 @@ namespace GridKit template int Tgov1::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } // Allocate local component data auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } // Resize signal variable data ws_.resize(1); @@ -180,13 +191,6 @@ namespace GridKit ws_[0] = 0.0; ws_indices_[0] = INVALID_INDEX; - // Default variable and residual index mapping to local index - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - // Set output signals if (signals_.template isAssigned()) { @@ -279,7 +283,7 @@ namespace GridKit template int Tgov1::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp index bdb001d40..7fcd7eaef 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp @@ -82,26 +82,32 @@ namespace GridKit template int LoadZ::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } + // std::cout << "Allocate Load..." << std::endl; + auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); - abs_tol_.resize(size); - tag_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } // Resize coupling data wb_.resize(2); h_.resize(2); - // Default variable and residual index mapping to local index - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - return 0; } @@ -153,7 +159,7 @@ namespace GridKit template int LoadZ::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp index a41b2df89..7daf01168 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp @@ -102,26 +102,32 @@ namespace GridKit template int LoadZIP::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } + // std::cout << "Allocate Load..." << std::endl; + auto size = static_cast(size_); // avoid compiler warnings - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } // Resize coupling data wb_.resize(2); h_.resize(2); - // Default variable and residual index mapping to local index - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - return 0; } @@ -174,7 +180,7 @@ namespace GridKit template int LoadZIP::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index 97372cf06..c8ba75531 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -156,26 +156,31 @@ namespace GridKit template int Ieeest::allocate() { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } ws_.resize(1); ws_indices_.resize(1); ws_[0] = 0.0; ws_indices_[0] = INVALID_INDEX; - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - if (signals_.template isAssigned()) { signals_.template getSignalNode()->set( @@ -259,7 +264,7 @@ namespace GridKit template int Ieeest::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index cbad7dd0a..75d0b797d 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -328,15 +328,25 @@ namespace GridKit template int Genrou::allocate() { - // Resize component model data + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); - abs_tol_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } // Resize bus data wb_.resize(2); @@ -348,13 +358,6 @@ namespace GridKit ws_indices_[0] = INVALID_INDEX; ws_indices_[1] = INVALID_INDEX; - // Default variable and residual index mapping to local index - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - // Set output signals if (signals_.template isAssigned()) { @@ -543,7 +546,7 @@ namespace GridKit template int Genrou::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp index 793f4af39..1bbe32d6f 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp @@ -197,15 +197,25 @@ namespace GridKit template int Gensal::allocate() { - // Resize component model data + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } // Resize bus data wb_.resize(2); @@ -217,13 +227,6 @@ namespace GridKit ws_indices_[0] = INVALID_INDEX; ws_indices_[1] = INVALID_INDEX; - // Default variable and residual index mapping to local index - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - // Set output signals if (signals_.template isAssigned()) { @@ -366,7 +369,7 @@ namespace GridKit template int Gensal::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index bd5816a8a..2b63d5a9a 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -178,27 +178,30 @@ namespace GridKit template int GenClassical::allocate() { - // Resize component model data + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + variable_indices_.resize(size); residual_indices_.resize(size); + for (IdxT j = 0; j < size_; ++j) + { + variable_indices_[static_cast(j)] = this->offset_ + j; + residual_indices_[static_cast(j)] = this->offset_ + j; + } // Resize coupling data wb_.resize(2); h_.resize(2); - // Default variable and residual index mapping to local index - for (IdxT j = 0; j < size_; ++j) - { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); - } - return 0; } @@ -264,7 +267,7 @@ namespace GridKit template int GenClassical::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 2f0524535..38d5ed1f6 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -55,6 +55,9 @@ namespace GridKit using Component::residual_indices_; using Component::csr_jac_; using Component::map_to_csr_; + using Component::offset_; + using Component::allocated_; + using Component::allocateVectors; public: using ScalarT = scalar_type; @@ -92,7 +95,11 @@ namespace GridKit int evaluateResidual() override; int evaluateJacobian() override; - void updateVariables(); + CsrMatrixT* getCsrJacobian() const override + { + return csr_jac_; + } + void updateTime(RealT t, RealT a) override; void addBus(BusT* bus); diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 287b328dc..716ed49ed 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -1,4 +1,5 @@ #include +#include #include #include @@ -424,11 +425,10 @@ namespace GridKit } /** - * @brief Allocate buses, components, and system objects. + * @brief Allocate buses, components, and system metadata. * - * 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 the system size, allocates the system vectors, + * binds child models to contiguous slices, and allocates child internals. * * @post size_ >= 1 * @@ -437,45 +437,62 @@ namespace GridKit int SystemModel::allocate() { 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_); - tag_.resize(size_); - abs_tol_.resize(size_); - variable_indices_.resize(size_); - residual_indices_.resize(size_); + auto size = static_cast(size_); + + if (!allocated_) + { + allocateVectors(size_); + } + + assert(y_.size() == size); + assert(yp_.size() == size); + assert(f_.size() == size); + assert(tag_.size() == size); + assert(abs_tol_.size() == size); + + IdxT child_offset = 0; + for (const auto& bus : buses_) + { + bus->bind(y_, yp_, f_, tag_, abs_tol_, child_offset); + child_offset += bus->size(); + } + + for (const auto& component : components_) + { + component->bind(y_, yp_, f_, tag_, abs_tol_, child_offset); + child_offset += component->size(); + } + + assert(child_offset == size_); - // Default variable and residual index mapping to local index + variable_indices_.resize(size); + residual_indices_.resize(size); for (IdxT j = 0; j < size_; ++j) { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); + variable_indices_[static_cast(j)] = offset_ + j; + residual_indices_[static_cast(j)] = offset_ + j; + } + + assert(variable_indices_.size() == static_cast(size_)); + assert(residual_indices_.size() == static_cast(size_)); + + for (const auto& bus : buses_) + { + bus->allocate(); + } + + for (const auto& component : components_) + { + component->allocate(); } // Verify component configuration @@ -572,43 +589,25 @@ 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. */ template int SystemModel::initialize() { + // Buses initialize first: components may read bus values during their own + // initialization. Writes land directly in the global vector storage. 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; } @@ -675,19 +674,11 @@ namespace GridKit for (const auto& bus : buses_) { bus->tagDifferentiable(); - for (IdxT j = 0; j < bus->size(); ++j) - { - tag_[bus->getVariableIndex(j)] = bus->tag()[j]; - } } for (const auto& component : components_) { component->tagDifferentiable(); - for (IdxT j = 0; j < component->size(); ++j) - { - tag_[component->getVariableIndex(j)] = component->tag()[j]; - } } return 0; @@ -706,25 +697,14 @@ namespace GridKit template int SystemModel::setAbsoluteTolerance(RealT rel_tol) { - IdxT offset = 0; for (const auto& bus : buses_) { bus->setAbsoluteTolerance(rel_tol); - for (IdxT j = 0; j < bus->size(); ++j) - { - abs_tol_[offset + j] = bus->absoluteTolerance()[j]; - } - offset += bus->size(); } for (const auto& component : components_) { component->setAbsoluteTolerance(rel_tol); - for (IdxT j = 0; j < component->size(); ++j) - { - abs_tol_[offset + j] = component->absoluteTolerance()[j]; - } - offset += component->size(); } return 0; @@ -733,26 +713,20 @@ 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. */ template int SystemModel::evaluateResidual() { - updateVariables(); - for (const auto& bus : buses_) { bus->evaluateResidual(); @@ -763,23 +737,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; } @@ -978,30 +935,6 @@ namespace GridKit return 0; } - /** - * @brief Update variables in buses and components - */ - template - void SystemModel::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 * @@ -1015,8 +948,6 @@ namespace GridKit { component->updateTime(t, a); } - - updateVariables(); } /** diff --git a/GridKit/Model/PowerElectronics/CMakeLists.txt b/GridKit/Model/PowerElectronics/CMakeLists.txt index fee2e752f..98dce1148 100644 --- a/GridKit/Model/PowerElectronics/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/CMakeLists.txt @@ -2,6 +2,9 @@ add_library(power_electronics_circuit_node INTERFACE) target_include_directories( power_electronics_circuit_node INTERFACE $ $) +target_link_libraries( + power_electronics_circuit_node + INTERFACE GridKit::dense_vector) add_library(GridKit::power_electronics_circuit_node ALIAS power_electronics_circuit_node) diff --git a/GridKit/Model/PowerElectronics/Capacitor/CMakeLists.txt b/GridKit/Model/PowerElectronics/Capacitor/CMakeLists.txt index 49aad14b6..4a7bd323c 100644 --- a/GridKit/Model/PowerElectronics/Capacitor/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/Capacitor/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_capacitor SOURCES Capacitor.cpp - HEADERS Capacitor.hpp) + HEADERS Capacitor.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index 995135f08..15e8ca419 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -17,9 +17,19 @@ namespace GridKit template class CircuitComponent : public Model::Evaluator { + protected: + using Model::Evaluator::y_; + using Model::Evaluator::yp_; + using Model::Evaluator::f_; + using Model::Evaluator::tag_; + using Model::Evaluator::abs_tol_; + using Model::Evaluator::allocated_; + using Model::Evaluator::allocateVectors; + public: using RealT = typename Model::Evaluator::RealT; using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; + using VectorT = typename Model::Evaluator::VectorT; CircuitComponent() = default; @@ -96,10 +106,10 @@ namespace GridKit connection_nodes_ = std::make_unique(static_cast(size_)); - y_.resize(static_cast(size_)); - yp_.resize(static_cast(size_)); - abs_tol_.resize(static_cast(size_)); - f_.resize(static_cast(size_)); + if (!allocated_) + { + allocateVectors(size_); + } return 0; } @@ -259,132 +269,132 @@ namespace GridKit return size_opt_; } - std::vector& y() final + VectorT& y() final { return y_; } - const std::vector& y() const final + const VectorT& y() const final { return y_; } - std::vector& yp() final + VectorT& yp() final { return yp_; } - const std::vector& yp() const final + const VectorT& yp() const final { return yp_; } - std::vector& tag() final + VectorT& tag() final { return tag_; } - const std::vector& tag() const final + const VectorT& tag() const final { return tag_; } - std::vector& absoluteTolerance() final + VectorT& absoluteTolerance() final { return abs_tol_; } - const std::vector& absoluteTolerance() const final + const VectorT& absoluteTolerance() const final { return abs_tol_; } - std::vector& yB() final + VectorT& yB() final { return yB_; } - const std::vector& yB() const final + const VectorT& yB() const final { return yB_; } - std::vector& ypB() final + VectorT& ypB() final { return ypB_; } - const std::vector& ypB() const final + const VectorT& ypB() const final { return ypB_; } - std::vector& param() final + VectorT& param() final { return param_; } - const std::vector& param() const final + const VectorT& param() const final { return param_; } - std::vector& param_up() final + VectorT& param_up() final { return param_up_; } - const std::vector& param_up() const final + const VectorT& param_up() const final { return param_up_; } - std::vector& param_lo() final + VectorT& param_lo() final { return param_lo_; } - const std::vector& param_lo() const final + const VectorT& param_lo() const final { return param_lo_; } - std::vector& getResidual() final + VectorT& getResidual() final { return f_; } - const std::vector& getResidual() const final + const VectorT& getResidual() const final { return f_; } - std::vector& getIntegrand() final + VectorT& getIntegrand() final { return g_; } - const std::vector& getIntegrand() const final + const VectorT& getIntegrand() const final { return g_; } - std::vector& getAdjointResidual() final + VectorT& getAdjointResidual() final { return fB_; } - const std::vector& getAdjointResidual() const final + const VectorT& getAdjointResidual() const final { return fB_; } - std::vector& getAdjointIntegrand() final + VectorT& getAdjointIntegrand() final { return gB_; } - const std::vector& getAdjointIntegrand() const final + const VectorT& getAdjointIntegrand() const final { return gB_; } @@ -423,21 +433,16 @@ namespace GridKit /// @brief A pointer to the internal residuals of this component ScalarT* f_int_; - std::vector y_; - std::vector yp_; - std::vector tag_; - std::vector abs_tol_; - std::vector f_; - std::vector g_; - - std::vector yB_; - std::vector ypB_; - std::vector fB_; - std::vector gB_; - - std::vector param_; - std::vector param_up_; - std::vector param_lo_; + VectorT g_; + + VectorT yB_; + VectorT ypB_; + VectorT fB_; + VectorT gB_; + + VectorT param_; + VectorT param_up_; + VectorT param_lo_; RealT time_; RealT alpha_; diff --git a/GridKit/Model/PowerElectronics/CircuitNode.hpp b/GridKit/Model/PowerElectronics/CircuitNode.hpp index 07b14b9df..16d54ce12 100644 --- a/GridKit/Model/PowerElectronics/CircuitNode.hpp +++ b/GridKit/Model/PowerElectronics/CircuitNode.hpp @@ -15,7 +15,16 @@ namespace GridKit template class CircuitNode : public Model::Evaluator { - using RealT = typename Model::Evaluator::RealT; + using Model::Evaluator::y_; + using Model::Evaluator::yp_; + using Model::Evaluator::f_; + using Model::Evaluator::tag_; + using Model::Evaluator::abs_tol_; + using Model::Evaluator::allocated_; + using Model::Evaluator::allocateVectors; + + using RealT = typename Model::Evaluator::RealT; + using VectorT = typename Model::Evaluator::VectorT; public: CircuitNode() @@ -67,13 +76,10 @@ namespace GridKit // Allocate storage for a single-node voltage and KCL residual int allocate() { - size_t size = static_cast(size_); - - y_.resize(size); - yp_.resize(size); - f_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + if (!allocated_) + { + allocateVectors(size_); + } variable_indices_[0] = 0; residual_indices_[0] = 0; @@ -116,7 +122,7 @@ namespace GridKit */ int setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } @@ -178,21 +184,15 @@ namespace GridKit std::map variable_indices_; std::map residual_indices_; - std::vector y_; - std::vector yp_; - std::vector tag_; - std::vector abs_tol_; - std::vector f_; - - std::vector g_{}; - std::vector param_{}; - std::vector param_up_{}; - std::vector param_lo_{}; + VectorT g_{}; + VectorT param_{}; + VectorT param_up_{}; + VectorT param_lo_{}; - std::vector yB_{}; - std::vector ypB_{}; - std::vector fB_{}; - std::vector gB_{}; + VectorT yB_{}; + VectorT ypB_{}; + VectorT fB_{}; + VectorT gB_{}; RealT time_{0}; RealT alpha_{0}; @@ -225,132 +225,132 @@ namespace GridKit // No time to update in node models } - std::vector& y() final + VectorT& y() final { return y_; } - const std::vector& y() const final + const VectorT& y() const final { return y_; } - std::vector& yp() final + VectorT& yp() final { return yp_; } - const std::vector& yp() const final + const VectorT& yp() const final { return yp_; } - std::vector& tag() final + VectorT& tag() final { return tag_; } - const std::vector& tag() const final + const VectorT& tag() const final { return tag_; } - std::vector& absoluteTolerance() final + VectorT& absoluteTolerance() final { return abs_tol_; } - const std::vector& absoluteTolerance() const final + const VectorT& absoluteTolerance() const final { return abs_tol_; } - std::vector& yB() final + VectorT& yB() final { return yB_; } - const std::vector& yB() const final + const VectorT& yB() const final { return yB_; } - std::vector& ypB() final + VectorT& ypB() final { return ypB_; } - const std::vector& ypB() const final + const VectorT& ypB() const final { return ypB_; } - std::vector& param() final + VectorT& param() final { return param_; } - const std::vector& param() const final + const VectorT& param() const final { return param_; } - std::vector& param_up() final + VectorT& param_up() final { return param_up_; } - const std::vector& param_up() const final + const VectorT& param_up() const final { return param_up_; } - std::vector& param_lo() final + VectorT& param_lo() final { return param_lo_; } - const std::vector& param_lo() const final + const VectorT& param_lo() const final { return param_lo_; } - std::vector& getResidual() final + VectorT& getResidual() final { return f_; } - const std::vector& getResidual() const final + const VectorT& getResidual() const final { return f_; } - std::vector& getIntegrand() final + VectorT& getIntegrand() final { return g_; } - const std::vector& getIntegrand() const final + const VectorT& getIntegrand() const final { return g_; } - std::vector& getAdjointResidual() final + VectorT& getAdjointResidual() final { return fB_; } - const std::vector& getAdjointResidual() const final + const VectorT& getAdjointResidual() const final { return fB_; } - std::vector& getAdjointIntegrand() final + VectorT& getAdjointIntegrand() final { return gB_; } - const std::vector& getAdjointIntegrand() const final + const VectorT& getAdjointIntegrand() const final { return gB_; } diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt b/GridKit/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt index 8443f1c62..f42fe6a29 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_disgen SOURCES DistributedGenerator.cpp - HEADERS DistributedGenerator.hpp) + HEADERS DistributedGenerator.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index 797022531..94108bebd 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -92,7 +92,7 @@ namespace GridKit template int DistributedGenerator::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerElectronics/InductionMotor/CMakeLists.txt b/GridKit/Model/PowerElectronics/InductionMotor/CMakeLists.txt index 46e2f645c..df4e46fa4 100644 --- a/GridKit/Model/PowerElectronics/InductionMotor/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/InductionMotor/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_inductionmotor SOURCES InductionMotor.cpp - HEADERS InductionMotor.hpp) + HEADERS InductionMotor.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/Inductor/CMakeLists.txt b/GridKit/Model/PowerElectronics/Inductor/CMakeLists.txt index 8bb22f467..b00478d2d 100644 --- a/GridKit/Model/PowerElectronics/Inductor/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/Inductor/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_inductor SOURCES Inductor.cpp - HEADERS Inductor.hpp) + HEADERS Inductor.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp index 5a8527c55..9125d73f1 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp @@ -66,7 +66,7 @@ namespace GridKit template int Inductor::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerElectronics/LinearTransformer/CMakeLists.txt b/GridKit/Model/PowerElectronics/LinearTransformer/CMakeLists.txt index 0c2c9b7f7..8a2bcb6dd 100644 --- a/GridKit/Model/PowerElectronics/LinearTransformer/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/LinearTransformer/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_lineartrasnformer SOURCES LinearTransformer.cpp - HEADERS LinearTransformer.hpp) + HEADERS LinearTransformer.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt b/GridKit/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt index 2cb3c90e4..7d3626bd0 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_microbusdq SOURCES MicrogridBusDQ.cpp - HEADERS MicrogridBusDQ.hpp) + HEADERS MicrogridBusDQ.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/CMakeLists.txt b/GridKit/Model/PowerElectronics/MicrogridLine/CMakeLists.txt index d5d0b6bee..e4ce2b4e8 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/MicrogridLine/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_microline SOURCES MicrogridLine.cpp - HEADERS MicrogridLine.hpp) + HEADERS MicrogridLine.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index c0896889a..35ca1bebd 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -80,7 +80,7 @@ namespace GridKit template int MicrogridLine::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt b/GridKit/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt index bf45b8ede..83946d13c 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_microload SOURCES MicrogridLoad.cpp - HEADERS MicrogridLoad.hpp) + HEADERS MicrogridLoad.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index d3e9a143d..86808d49f 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -76,7 +76,7 @@ namespace GridKit template int MicrogridLoad::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerElectronics/NodeBase.hpp b/GridKit/Model/PowerElectronics/NodeBase.hpp index 3147f6817..d791328d8 100644 --- a/GridKit/Model/PowerElectronics/NodeBase.hpp +++ b/GridKit/Model/PowerElectronics/NodeBase.hpp @@ -11,8 +11,18 @@ namespace GridKit template class NodeBase : public Model::Evaluator { + protected: + using Model::Evaluator::y_; + using Model::Evaluator::yp_; + using Model::Evaluator::f_; + using Model::Evaluator::tag_; + using Model::Evaluator::abs_tol_; + using Model::Evaluator::allocated_; + using Model::Evaluator::allocateVectors; + public: - using RealT = typename Model::Evaluator::RealT; + using RealT = typename Model::Evaluator::RealT; + using VectorT = typename Model::Evaluator::VectorT; NodeBase(size_t n_intern, size_t n_extern) : n_intern_(n_intern), n_extern_(n_extern) @@ -48,42 +58,42 @@ namespace GridKit { } - std::vector& y() final + VectorT& y() final { return y_; } - const std::vector& y() const final + const VectorT& y() const final { return y_; } - std::vector& yp() final + VectorT& yp() final { return yp_; } - const std::vector& yp() const final + const VectorT& yp() const final { return yp_; } - std::vector& tag() final + VectorT& tag() final { return tag_; } - const std::vector& tag() const final + const VectorT& tag() const final { return tag_; } - std::vector& absoluteTolerance() final + VectorT& absoluteTolerance() final { return abs_tol_; } - const std::vector& absoluteTolerance() const final + const VectorT& absoluteTolerance() const final { return abs_tol_; } @@ -127,15 +137,13 @@ namespace GridKit int allocate() override { - // Temporary while we use std::vector in the code size_t size = static_cast(n_intern_ + n_extern_); - // Resize component model data - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + if (!allocated_) + { + allocateVectors(static_cast(size)); + } + variable_indices_.resize(size); residual_indices_.resize(size); @@ -163,7 +171,7 @@ namespace GridKit */ int setAbsoluteTolerance(RealT rel_tol) final { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } @@ -193,12 +201,6 @@ namespace GridKit std::vector variable_indices_; ///< Global (system-level) variable indices std::vector residual_indices_; ///< Global (system-level) residual indices - std::vector y_; - std::vector yp_; - std::vector tag_; - std::vector abs_tol_; - std::vector f_; - IdxT* J_rows_buffer_{nullptr}; IdxT* J_cols_buffer_{nullptr}; RealT* J_vals_buffer_{nullptr}; @@ -207,15 +209,15 @@ namespace GridKit // Adjoint sensitivity members // - std::vector g_{}; - std::vector yB_{}; - std::vector ypB_{}; - std::vector fB_{}; - std::vector gB_{}; + VectorT g_{}; + VectorT yB_{}; + VectorT ypB_{}; + VectorT fB_{}; + VectorT gB_{}; - std::vector param_{}; - std::vector param_up_{}; - std::vector param_lo_{}; + VectorT param_{}; + VectorT param_up_{}; + VectorT param_lo_{}; std::unique_ptr connection_nodes_; @@ -232,61 +234,61 @@ namespace GridKit return 0; } - std::vector& yB() final + VectorT& yB() final { throw "ERROR: Method not implemented!\n"; return yB_; } - const std::vector& yB() const final + const VectorT& yB() const final { throw "ERROR: Method not implemented!\n"; return yB_; } - std::vector& ypB() final + VectorT& ypB() final { throw "ERROR: Method not implemented!\n"; return ypB_; } - const std::vector& ypB() const final + const VectorT& ypB() const final { throw "ERROR: Method not implemented!\n"; return ypB_; } - std::vector& param() final + VectorT& param() final { throw "ERROR: Method not implemented!\n"; return param_; } - const std::vector& param() const final + const VectorT& param() const final { throw "ERROR: Method not implemented!\n"; return param_; } - std::vector& param_up() final + VectorT& param_up() final { throw "ERROR: Method not implemented!\n"; return param_up_; } - const std::vector& param_up() const final + const VectorT& param_up() const final { throw "ERROR: Method not implemented!\n"; return param_up_; } - std::vector& param_lo() final + VectorT& param_lo() final { throw "ERROR: Method not implemented!\n"; return param_lo_; } - const std::vector& param_lo() const final + const VectorT& param_lo() const final { throw "ERROR: Method not implemented!\n"; return param_lo_; @@ -316,47 +318,47 @@ namespace GridKit return 1; } - std::vector& getResidual() final + VectorT& getResidual() final { return f_; } - const std::vector& getResidual() const final + const VectorT& getResidual() const final { return f_; } - std::vector& getIntegrand() final + VectorT& getIntegrand() final { throw "ERROR: Method not implemented!\n"; return g_; } - const std::vector& getIntegrand() const final + const VectorT& getIntegrand() const final { throw "ERROR: Method not implemented!\n"; return g_; } - std::vector& getAdjointResidual() final + VectorT& getAdjointResidual() final { throw "ERROR: Method not implemented!\n"; return fB_; } - const std::vector& getAdjointResidual() const final + const VectorT& getAdjointResidual() const final { throw "ERROR: Method not implemented!\n"; return fB_; } - std::vector& getAdjointIntegrand() final + VectorT& getAdjointIntegrand() final { throw "ERROR: Method not implemented!\n"; return gB_; } - const std::vector& getAdjointIntegrand() const final + const VectorT& getAdjointIntegrand() const final { throw "ERROR: Method not implemented!\n"; return gB_; diff --git a/GridKit/Model/PowerElectronics/Resistor/CMakeLists.txt b/GridKit/Model/PowerElectronics/Resistor/CMakeLists.txt index aa3c25a80..3bed9d678 100644 --- a/GridKit/Model/PowerElectronics/Resistor/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/Resistor/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_resistor SOURCES Resistor.cpp - HEADERS Resistor.hpp) + HEADERS Resistor.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp index 87d86cab4..979c08834 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp @@ -67,7 +67,7 @@ namespace GridKit template int Resistor::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt b/GridKit/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt index 9a7b81528..4e685097a 100644 --- a/GridKit/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_synmachine SOURCES SynchronousMachine.cpp - HEADERS SynchronousMachine.hpp) + HEADERS SynchronousMachine.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 62cbfd96f..7c6bc9ef4 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include @@ -16,6 +17,48 @@ namespace GridKit { + /** + * Writes a vector to a file in Matrix Market format + * + * @param vec The vector to write + * @param filename The name of the output file + * @param header Additional header information/comments + * @return true if the write was successful, false otherwise + */ + template + void writeVectorToMatrixMarket(const VectorT& vec, const std::string& filename, const std::string& header) + { + std::ofstream outFile(filename); + + if (!outFile.is_open()) + { + std::cerr << "Error: Could not open file " << filename << " for writing." << std::endl; + return; + } + + // Uncomment to write Matrix Market header + // outFile << "%%MatrixMarket vector array real general" << std::endl; + + // Write additional header information as comments + if (!header.empty()) + { + outFile << "% " << header << std::endl; + } + + // Write the vector size + outFile << vec.size() << std::endl; + + // Write the vector elements + outFile << std::scientific << std::setprecision(16); + for (std::size_t i = 0; i < vec.size(); ++i) + { + outFile << vec.data()[i] << std::endl; + } + + outFile.close(); + return; + } + template class PowerElectronicsModel : public CircuitComponent { @@ -38,6 +81,8 @@ namespace GridKit using CircuitComponent::f_int_; using CircuitComponent::tag_; using CircuitComponent::abs_tol_; + using CircuitComponent::allocated_; + using CircuitComponent::allocateVectors; public: /** @@ -136,12 +181,10 @@ namespace GridKit n_extern_ = 0; size_ = n_intern_ + n_extern_; - // Allocate global vectors - y_.resize(size_); - yp_.resize(size_); - f_.resize(size_); - tag_.resize(size_); - abs_tol_.resize(size_); + if (!allocated_) + { + allocateVectors(static_cast(size_)); + } { // Start node internal indexing after all component internals for proper KLU ordering size_t node_internal_idx = component_internal_size; @@ -291,8 +334,8 @@ namespace GridKit { for (component_type* component : components_) { - std::vector& y = component->y(); - std::vector& yp = component->yp(); + auto& y = component->y(); + auto& yp = component->yp(); const std::set& externals = component->getExternIndices(); for (size_t j : externals) @@ -331,7 +374,7 @@ namespace GridKit */ int setAbsoluteTolerance(RealT rel_tol) final { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } @@ -363,8 +406,8 @@ namespace GridKit if (int err_code = component->evaluateExternalResidual()) return err_code; - const std::vector& residual = component->getResidual(); - const std::set& externals = component->getExternIndices(); + auto& residual = component->getResidual(); + const std::set& externals = component->getExternIndices(); for (size_t j : externals) { diff --git a/GridKit/Model/PowerElectronics/TransmissionLine/CMakeLists.txt b/GridKit/Model/PowerElectronics/TransmissionLine/CMakeLists.txt index d24e46a27..5a924a3b4 100644 --- a/GridKit/Model/PowerElectronics/TransmissionLine/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/TransmissionLine/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_tranline SOURCES TransmissionLine.cpp - HEADERS TransmissionLine.hpp) + HEADERS TransmissionLine.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/VoltageSource/CMakeLists.txt b/GridKit/Model/PowerElectronics/VoltageSource/CMakeLists.txt index cb1e2ed81..07be55ffd 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/VoltageSource/CMakeLists.txt @@ -1,4 +1,5 @@ gridkit_add_library( power_elec_voltagesource SOURCES VoltageSource.cpp - HEADERS VoltageSource.hpp) + HEADERS VoltageSource.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp index f2967c7bd..664e9ad65 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp @@ -67,7 +67,7 @@ namespace GridKit template int VoltageSource::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerFlow/Branch/CMakeLists.txt b/GridKit/Model/PowerFlow/Branch/CMakeLists.txt index ed0e7baf7..b1a982071 100644 --- a/GridKit/Model/PowerFlow/Branch/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Branch/CMakeLists.txt @@ -6,4 +6,5 @@ gridkit_add_library( branch SOURCES Branch.cpp - HEADERS Branch.hpp) + HEADERS Branch.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/Bus/BusPQ.cpp b/GridKit/Model/PowerFlow/Bus/BusPQ.cpp index 4cc8a450c..317ba3f60 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPQ.cpp +++ b/GridKit/Model/PowerFlow/Bus/BusPQ.cpp @@ -69,16 +69,22 @@ namespace GridKit template int BusPQ::allocate() { + using VectorT = typename ModelEvaluatorImpl::VectorT; + // std::cout << "Allocate PQ bus ..." << std::endl; - f_.resize(static_cast(size_)); - y_.resize(static_cast(size_)); - yp_.resize(static_cast(size_)); - tag_.resize(static_cast(size_)); - abs_tol_.resize(static_cast(size_)); + this->allocateVectors(size_); + + fB_ = VectorT(size_); + fB_.allocate(GridKit::LinearAlgebra::memory::HOST); + fB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + yB_ = VectorT(size_); + yB_.allocate(GridKit::LinearAlgebra::memory::HOST); + yB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); - fB_.resize(static_cast(size_)); - yB_.resize(static_cast(size_)); - ypB_.resize(static_cast(size_)); + ypB_ = VectorT(size_); + ypB_.allocate(GridKit::LinearAlgebra::memory::HOST); + ypB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); return 0; } @@ -106,7 +112,7 @@ namespace GridKit template int BusPQ::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerFlow/Bus/BusPV.cpp b/GridKit/Model/PowerFlow/Bus/BusPV.cpp index bde6dd49b..73041582b 100644 --- a/GridKit/Model/PowerFlow/Bus/BusPV.cpp +++ b/GridKit/Model/PowerFlow/Bus/BusPV.cpp @@ -67,16 +67,22 @@ namespace GridKit template int BusPV::allocate() { + using VectorT = typename ModelEvaluatorImpl::VectorT; + // std::cout << "Allocate PV bus ..." << std::endl; - f_.resize(static_cast(size_)); - y_.resize(static_cast(size_)); - yp_.resize(static_cast(size_)); - tag_.resize(static_cast(size_)); - abs_tol_.resize(static_cast(size_)); + this->allocateVectors(size_); + + fB_ = VectorT(size_); + fB_.allocate(GridKit::LinearAlgebra::memory::HOST); + fB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + yB_ = VectorT(size_); + yB_.allocate(GridKit::LinearAlgebra::memory::HOST); + yB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); - fB_.resize(static_cast(size_)); - yB_.resize(static_cast(size_)); - ypB_.resize(static_cast(size_)); + ypB_ = VectorT(size_); + ypB_.allocate(GridKit::LinearAlgebra::memory::HOST); + ypB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); return 0; } @@ -103,7 +109,7 @@ namespace GridKit template int BusPV::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerFlow/Bus/CMakeLists.txt b/GridKit/Model/PowerFlow/Bus/CMakeLists.txt index f06511289..0844657ce 100644 --- a/GridKit/Model/PowerFlow/Bus/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Bus/CMakeLists.txt @@ -12,4 +12,5 @@ gridkit_add_library( BusFactory.hpp BusPQ.hpp BusPV.hpp - BusSlack.hpp) + BusSlack.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/Generator/CMakeLists.txt b/GridKit/Model/PowerFlow/Generator/CMakeLists.txt index 1136651f2..d010f9e89 100644 --- a/GridKit/Model/PowerFlow/Generator/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Generator/CMakeLists.txt @@ -10,4 +10,5 @@ gridkit_add_library( GeneratorFactory.hpp GeneratorPQ.hpp GeneratorPV.hpp - GeneratorSlack.hpp) + GeneratorSlack.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/Generator2/CMakeLists.txt b/GridKit/Model/PowerFlow/Generator2/CMakeLists.txt index 0c64cc1a7..35e15178d 100644 --- a/GridKit/Model/PowerFlow/Generator2/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Generator2/CMakeLists.txt @@ -6,4 +6,5 @@ gridkit_add_library( generator2 SOURCES Generator2.cpp - HEADERS Generator2.hpp) + HEADERS Generator2.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/Generator2/Generator2.cpp b/GridKit/Model/PowerFlow/Generator2/Generator2.cpp index eb8bd440e..cdeaeb411 100644 --- a/GridKit/Model/PowerFlow/Generator2/Generator2.cpp +++ b/GridKit/Model/PowerFlow/Generator2/Generator2.cpp @@ -49,8 +49,6 @@ namespace GridKit template int Generator2::allocate() { - tag_.resize(static_cast(size_)); - abs_tol_.resize(static_cast(size_)); return 0; } @@ -65,7 +63,7 @@ namespace GridKit template int Generator2::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerFlow/Generator4/CMakeLists.txt b/GridKit/Model/PowerFlow/Generator4/CMakeLists.txt index ac6252e2c..3abacf840 100644 --- a/GridKit/Model/PowerFlow/Generator4/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Generator4/CMakeLists.txt @@ -8,4 +8,5 @@ gridkit_add_library( generator4 SOURCES Generator4.cpp - HEADERS Generator4.hpp) + HEADERS Generator4.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/Generator4/Generator4.cpp b/GridKit/Model/PowerFlow/Generator4/Generator4.cpp index 296a56867..1f6c6e604 100644 --- a/GridKit/Model/PowerFlow/Generator4/Generator4.cpp +++ b/GridKit/Model/PowerFlow/Generator4/Generator4.cpp @@ -57,9 +57,6 @@ namespace GridKit int Generator4::allocate() { // std::cout << "Allocate Generator4..." << std::endl; - tag_.resize(static_cast(size_)); - abs_tol_.resize(static_cast(size_)); - return 0; } @@ -165,7 +162,7 @@ namespace GridKit template int Generator4::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerFlow/Generator4Governor/CMakeLists.txt b/GridKit/Model/PowerFlow/Generator4Governor/CMakeLists.txt index 41ec72b07..c70ffbe93 100644 --- a/GridKit/Model/PowerFlow/Generator4Governor/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Generator4Governor/CMakeLists.txt @@ -8,4 +8,5 @@ gridkit_add_library( generator4governor SOURCES Generator4Governor.cpp - HEADERS Generator4Governor.hpp) + HEADERS Generator4Governor.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp index 6caaa9674..f91f0cecd 100644 --- a/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp +++ b/GridKit/Model/PowerFlow/Generator4Governor/Generator4Governor.cpp @@ -67,9 +67,6 @@ namespace GridKit int Generator4Governor::allocate() { // std::cout << "Allocate Gen2..." << std::endl; - tag_.resize(static_cast(size_)); - abs_tol_.resize(static_cast(size_)); - return 0; } @@ -172,7 +169,7 @@ namespace GridKit template int Generator4Governor::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerFlow/Generator4Param/CMakeLists.txt b/GridKit/Model/PowerFlow/Generator4Param/CMakeLists.txt index 48f7f5b39..33fd2c59a 100644 --- a/GridKit/Model/PowerFlow/Generator4Param/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Generator4Param/CMakeLists.txt @@ -8,4 +8,5 @@ gridkit_add_library( generator4param SOURCES Generator4Param.cpp - HEADERS Generator4Param.hpp) + HEADERS Generator4Param.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp index e08e33f73..aa17505eb 100644 --- a/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp +++ b/GridKit/Model/PowerFlow/Generator4Param/Generator4Param.cpp @@ -54,9 +54,6 @@ namespace GridKit int Generator4Param::allocate() { // std::cout << "Allocate Generator4Param..." << std::endl; - tag_.resize(static_cast(size_)); - abs_tol_.resize(static_cast(size_)); - return 0; } @@ -154,7 +151,7 @@ namespace GridKit template int Generator4Param::setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Model/PowerFlow/Load/CMakeLists.txt b/GridKit/Model/PowerFlow/Load/CMakeLists.txt index c00bfc959..57e4eb5b6 100644 --- a/GridKit/Model/PowerFlow/Load/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/Load/CMakeLists.txt @@ -6,4 +6,5 @@ gridkit_add_library( load SOURCES Load.cpp - HEADERS Load.hpp) + HEADERS Load.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/MiniGrid/CMakeLists.txt b/GridKit/Model/PowerFlow/MiniGrid/CMakeLists.txt index e814b5852..d76d79ceb 100644 --- a/GridKit/Model/PowerFlow/MiniGrid/CMakeLists.txt +++ b/GridKit/Model/PowerFlow/MiniGrid/CMakeLists.txt @@ -6,4 +6,5 @@ gridkit_add_library( minigrid SOURCES MiniGrid.cpp - HEADERS MiniGrid.hpp) + HEADERS MiniGrid.hpp + LINK_LIBRARIES GridKit::dense_vector) diff --git a/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp b/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp index 7f071904a..a3805d95c 100644 --- a/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp +++ b/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp @@ -15,8 +15,17 @@ namespace GridKit template class ModelEvaluatorImpl : public Model::Evaluator { + protected: + using Model::Evaluator::y_; + using Model::Evaluator::yp_; + using Model::Evaluator::f_; + using Model::Evaluator::tag_; + using Model::Evaluator::abs_tol_; + using Model::Evaluator::allocateVectors; + public: - using RealT = typename Model::Evaluator::RealT; + using RealT = typename Model::Evaluator::RealT; + using VectorT = typename Model::Evaluator::VectorT; ModelEvaluatorImpl() : size_(0), @@ -27,20 +36,43 @@ namespace GridKit ModelEvaluatorImpl(IdxT size, IdxT size_quad, IdxT size_opt) : size_(size), + nnz_(0), size_quad_(size_quad), size_opt_(size_opt), - y_(static_cast(size_)), - yp_(static_cast(size_)), - f_(static_cast(size_)), - g_(static_cast(size_quad_)), - yB_(static_cast(size_)), - ypB_(static_cast(size_)), - fB_(static_cast(size_)), - gB_(static_cast(size_opt_)), - param_(static_cast(size_opt_)), - param_up_(static_cast(size_opt_)), - param_lo_(static_cast(size_opt_)) + g_(size_quad_), + yB_(size_), + ypB_(size_), + fB_(size_), + gB_(size_opt_), + param_(size_opt_), + param_up_(size_opt_), + param_lo_(size_opt_) { + allocateVectors(size_); + + g_.allocate(GridKit::LinearAlgebra::memory::HOST); + g_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + yB_.allocate(GridKit::LinearAlgebra::memory::HOST); + yB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + ypB_.allocate(GridKit::LinearAlgebra::memory::HOST); + ypB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + fB_.allocate(GridKit::LinearAlgebra::memory::HOST); + fB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + gB_.allocate(GridKit::LinearAlgebra::memory::HOST); + gB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + param_.allocate(GridKit::LinearAlgebra::memory::HOST); + param_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + param_up_.allocate(GridKit::LinearAlgebra::memory::HOST); + param_up_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + param_lo_.allocate(GridKit::LinearAlgebra::memory::HOST); + param_lo_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); } virtual IdxT size() @@ -80,132 +112,132 @@ namespace GridKit msa = max_steps_; } - std::vector& y() + VectorT& y() { return y_; } - const std::vector& y() const + const VectorT& y() const { return y_; } - std::vector& yp() + VectorT& yp() { return yp_; } - const std::vector& yp() const + const VectorT& yp() const { return yp_; } - std::vector& tag() + VectorT& tag() { return tag_; } - const std::vector& tag() const + const VectorT& tag() const { return tag_; } - std::vector& absoluteTolerance() + VectorT& absoluteTolerance() { return abs_tol_; } - const std::vector& absoluteTolerance() const + const VectorT& absoluteTolerance() const { return abs_tol_; } - std::vector& yB() + VectorT& yB() { return yB_; } - const std::vector& yB() const + const VectorT& yB() const { return yB_; } - std::vector& ypB() + VectorT& ypB() { return ypB_; } - const std::vector& ypB() const + const VectorT& ypB() const { return ypB_; } - std::vector& param() + VectorT& param() { return param_; } - const std::vector& param() const + const VectorT& param() const { return param_; } - std::vector& param_up() + VectorT& param_up() { return param_up_; } - const std::vector& param_up() const + const VectorT& param_up() const { return param_up_; } - std::vector& param_lo() + VectorT& param_lo() { return param_lo_; } - const std::vector& param_lo() const + const VectorT& param_lo() const { return param_lo_; } - std::vector& getResidual() + VectorT& getResidual() { return f_; } - const std::vector& getResidual() const + const VectorT& getResidual() const { return f_; } - std::vector& getIntegrand() + VectorT& getIntegrand() { return g_; } - const std::vector& getIntegrand() const + const VectorT& getIntegrand() const { return g_; } - std::vector& getAdjointResidual() + VectorT& getAdjointResidual() { return fB_; } - const std::vector& getAdjointResidual() const + const VectorT& getAdjointResidual() const { return fB_; } - std::vector& getAdjointIntegrand() + VectorT& getAdjointIntegrand() { return gB_; } - const std::vector& getAdjointIntegrand() const + const VectorT& getAdjointIntegrand() const { return gB_; } @@ -222,21 +254,16 @@ namespace GridKit IdxT size_quad_; IdxT size_opt_; - std::vector y_; - std::vector yp_; - std::vector tag_; - std::vector abs_tol_; - std::vector f_; - std::vector g_; - - std::vector yB_; - std::vector ypB_; - std::vector fB_; - std::vector gB_; - - std::vector param_; - std::vector param_up_; - std::vector param_lo_; + VectorT g_; + + VectorT yB_; + VectorT ypB_; + VectorT fB_; + VectorT gB_; + + VectorT param_; + VectorT param_up_; + VectorT param_lo_; RealT time_; RealT alpha_; diff --git a/GridKit/Model/PowerFlow/SystemModel.hpp b/GridKit/Model/PowerFlow/SystemModel.hpp index e174b4346..e7fa41563 100644 --- a/GridKit/Model/PowerFlow/SystemModel.hpp +++ b/GridKit/Model/PowerFlow/SystemModel.hpp @@ -27,6 +27,7 @@ namespace GridKit using bus_type = Model::Evaluator; using component_type = Model::Evaluator; using RealT = typename ModelEvaluatorImpl::RealT; + using VectorT = typename ModelEvaluatorImpl::VectorT; using ModelEvaluatorImpl::size_; using ModelEvaluatorImpl::size_quad_; @@ -101,21 +102,39 @@ namespace GridKit } // Allocate global vectors - y_.resize(size_); - yp_.resize(size_); - yB_.resize(size_); - ypB_.resize(size_); - f_.resize(size_); - fB_.resize(size_); - tag_.resize(size_); - abs_tol_.resize(size_); - - g_.resize(size_quad_); - gB_.resize(size_quad_ * size_opt_); - - param_.resize(size_opt_); - param_lo_.resize(size_opt_); - param_up_.resize(size_opt_); + this->allocateVectors(size_); + + yB_ = VectorT(size_); + yB_.allocate(GridKit::LinearAlgebra::memory::HOST); + yB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + ypB_ = VectorT(size_); + ypB_.allocate(GridKit::LinearAlgebra::memory::HOST); + ypB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + fB_ = VectorT(size_); + fB_.allocate(GridKit::LinearAlgebra::memory::HOST); + fB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + g_ = VectorT(size_quad_); + g_.allocate(GridKit::LinearAlgebra::memory::HOST); + g_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + gB_ = VectorT(size_quad_ * size_opt_); + gB_.allocate(GridKit::LinearAlgebra::memory::HOST); + gB_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + param_ = VectorT(size_opt_); + param_.allocate(GridKit::LinearAlgebra::memory::HOST); + param_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + param_lo_ = VectorT(size_opt_); + param_lo_.allocate(GridKit::LinearAlgebra::memory::HOST); + param_lo_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); + + param_up_ = VectorT(size_opt_); + param_up_.allocate(GridKit::LinearAlgebra::memory::HOST); + param_up_.setDataUpdated(GridKit::LinearAlgebra::memory::HOST); assert(size_quad_ == 1 or size_quad_ == 0); diff --git a/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp b/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp index 44b8a75fb..6ee079ed6 100644 --- a/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp +++ b/GridKit/Model/PowerFlow/SystemModelPowerFlow.hpp @@ -148,10 +148,7 @@ namespace GridKit } // Allocate global vectors - y_.resize(size_); - f_.resize(size_); - tag_.resize(size_); - abs_tol_.resize(size_); + this->allocateVectors(size_); return 0; } @@ -233,7 +230,7 @@ namespace GridKit */ int setAbsoluteTolerance(RealT rel_tol) { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index 9007bbf44..462783251 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -728,8 +728,7 @@ namespace AnalysisManager model->updateTime(tres, 0.0); model->evaluateResidual(); - const std::vector& f = model->getResidual(); - copyVec(f, rr); + copyVec(model->getResidual(), rr); return 0; } @@ -794,8 +793,7 @@ namespace AnalysisManager model->updateTime(tt, 0.0); model->evaluateIntegrand(); - const std::vector& g = model->getIntegrand(); - copyVec(g, rhsQ); + copyVec(model->getIntegrand(), rhsQ); return 0; } @@ -818,8 +816,7 @@ namespace AnalysisManager model->updateTime(tt, 0.0); model->evaluateAdjointResidual(); - const std::vector& fB = model->getAdjointResidual(); - copyVec(fB, rrB); + copyVec(model->getAdjointResidual(), rrB); return 0; } @@ -842,73 +839,51 @@ namespace AnalysisManager model->updateTime(tt, 0.0); model->evaluateAdjointIntegrand(); - const std::vector& gB = model->getAdjointIntegrand(); - copyVec(gB, rhsQB); + copyVec(model->getAdjointIntegrand(), rhsQB); return 0; } /** - * @brief Copy SUNDIALS N_Vector to std::vector + * @brief Copy SUNDIALS N_Vector to Vector * * @tparam ScalarT * @tparam IdxT */ template - void Ida::copyVec(const N_Vector x, std::vector& y) + void Ida::copyVec(const N_Vector x, VectorT& y) { const auto xsize = static_cast(N_VGetLength(x)); if (xsize != y.size()) { - std::cerr << "\nN_Vector size (" << xsize << ") does not match std::vector size (" + std::cerr << "\nN_Vector size (" << xsize << ") does not match vector size (" << y.size() << ").\n\n"; throw SundialsException(); } const ScalarT* xdata = N_VGetArrayPointer(x); - std::copy_n(xdata, y.size(), y.begin()); + std::copy_n(xdata, y.size(), y.data()); } /** - * @brief Copy std::vector to SUNDIALS N_Vector + * @brief Copy Vector to SUNDIALS N_Vector * * @tparam ScalarT * @tparam IdxT */ template - void Ida::copyVec(const std::vector& x, N_Vector y) + void Ida::copyVec(const VectorT& x, N_Vector y) { const auto ysize = static_cast(N_VGetLength(y)); if (x.size() != ysize) { - std::cerr << "\nstd::vector size (" << x.size() << ") does not match N_Vector size (" + std::cerr << "\nvector size (" << x.size() << ") does not match N_Vector size (" << ysize << ").\n\n"; throw SundialsException(); } ScalarT* ydata = N_VGetArrayPointer(y); - std::copy(x.cbegin(), x.cend(), ydata); - } - - /** - * @brief Copy std::vector to SUNDIALS N_Vector - * - * @tparam ScalarT - * @tparam IdxT - */ - template - void Ida::copyVec(const std::vector& x, N_Vector y) - { - const auto ysize = static_cast(N_VGetLength(y)); - if (x.size() != ysize) - { - std::cerr << "\nstd::vector size (" << x.size() << ") does not match N_Vector size (" - << ysize << ").\n\n"; - throw SundialsException(); - } - - ScalarT* ydata = N_VGetArrayPointer(y); - std::copy(x.cbegin(), x.cend(), ydata); + std::copy_n(x.data(), x.size(), ydata); } /** diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 1cdcd1269..fe421577f 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -45,7 +45,8 @@ namespace AnalysisManager { using DynamicSolver::model_; - using RealT = typename GridKit::ScalarTraits::RealT; + using RealT = typename GridKit::ScalarTraits::RealT; + using VectorT = typename GridKit::Model::Evaluator::VectorT; public: Ida(GridKit::Model::Evaluator* model); @@ -235,9 +236,8 @@ namespace AnalysisManager private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); - static void copyVec(const N_Vector x, std::vector& y); - static void copyVec(const std::vector& x, N_Vector y); - static void copyVec(const std::vector& x, N_Vector y); + static void copyVec(const N_Vector x, VectorT& y); + static void copyVec(const VectorT& x, N_Vector y); // int check_flag(void *flagvalue, const char *funcname, int opt); static void checkAllocation(void* v, const char* functionName); diff --git a/GridKit/Solver/SteadyState/Kinsol.cpp b/GridKit/Solver/SteadyState/Kinsol.cpp index aa04fe4e6..e7c77483b 100644 --- a/GridKit/Solver/SteadyState/Kinsol.cpp +++ b/GridKit/Solver/SteadyState/Kinsol.cpp @@ -139,24 +139,23 @@ namespace AnalysisManager copyVec(yy, model->y()); model->evaluateResidual(); - const std::vector& f = model->getResidual(); - copyVec(f, rr); + copyVec(model->getResidual(), rr); return 0; } template - void Kinsol::copyVec(const N_Vector x, std::vector& y) + void Kinsol::copyVec(const N_Vector x, VectorT& y) { const ScalarT* xdata = N_VGetArrayPointer(x); - std::copy_n(xdata, y.size(), y.begin()); + std::copy_n(xdata, y.size(), y.data()); } template - void Kinsol::copyVec(const std::vector& x, N_Vector y) + void Kinsol::copyVec(const VectorT& x, N_Vector y) { ScalarT* ydata = N_VGetArrayPointer(y); - std::copy_n(x.cbegin(), x.size(), ydata); + std::copy_n(x.data(), x.size(), ydata); } template diff --git a/GridKit/Solver/SteadyState/Kinsol.hpp b/GridKit/Solver/SteadyState/Kinsol.hpp index 89fc6c720..b867d936e 100644 --- a/GridKit/Solver/SteadyState/Kinsol.hpp +++ b/GridKit/Solver/SteadyState/Kinsol.hpp @@ -29,7 +29,8 @@ namespace AnalysisManager { using SteadyStateSolver::model_; - using RealT = typename GridKit::ScalarTraits::RealT; + using RealT = typename GridKit::ScalarTraits::RealT; + using VectorT = typename GridKit::Model::Evaluator::VectorT; static_assert(std::is_same_v, "RealT must be the same type as sunrealtype"); @@ -140,9 +141,8 @@ namespace AnalysisManager private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); - static void copyVec(const N_Vector x, std::vector& y); - static void copyVec(const std::vector& x, N_Vector y); - static void copyVec(const std::vector& x, N_Vector y); + static void copyVec(const N_Vector x, VectorT& y); + static void copyVec(const VectorT& x, N_Vector y); // int check_flag(void *flagvalue, const char *funcname, int opt); static inline void checkAllocation(void* v, const char* functionName); diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp index 4f687a5bc..65e3308db 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp @@ -123,7 +123,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& yval = sys.y(); + auto& yval = sys.y(); // Output time out << t << ","; diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp index e75edf180..823e67116 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp @@ -118,7 +118,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& yval = sys.y(); + auto& yval = sys.y(); // Output time out << t << ","; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp index 0204549d9..57abcad77 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp @@ -210,7 +210,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); // Bus 1 -> +2 // Bus 2 -> +2 diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp index 7e1bbc550..f4c2faf8f 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp @@ -187,7 +187,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); output.push_back(OutputData{t, 1 + static_cast(y_val[7]), diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp index 4892d1b1b..ba61a51e3 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp @@ -134,7 +134,7 @@ int main(int argc, const char* argv[]) auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); output.push_back(OutputData{t, 1 + static_cast(y_val[7]), diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/ZipLoad/ThreeBusZipLoadJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/ZipLoad/ThreeBusZipLoadJson.cpp index 3caf9d8bf..baa1c4ee4 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/ZipLoad/ThreeBusZipLoadJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/ZipLoad/ThreeBusZipLoadJson.cpp @@ -132,7 +132,7 @@ int main(int argc, const char* argv[]) auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); output.push_back(OutputData{t, 1 + static_cast(y_val[7]), diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp index 00d3fb551..6c5716097 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp @@ -133,7 +133,7 @@ int main() // push it into output, which is updated outside the callback. auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); output.push_back(OutputData{t, static_cast(y_val[0]), diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp index 1efb9ba3f..aca5e0661 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp @@ -113,7 +113,7 @@ int main(int argc, const char* argv[]) // push it into output, which is updated outside the callback. auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); output.push_back(OutputData{t, static_cast(y_val[0]), diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp index 432c05cce..b830ea07a 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp @@ -182,7 +182,7 @@ int main() // push it into output, which is updated outside the callback. auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); // Note Omega of gen is at state index 5! (Each added signal shifted by 1) // Bus -> 2 States diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp index 692e9431a..9bdf2225c 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp @@ -113,7 +113,7 @@ int main(int argc, const char* argv[]) // push it into output, which is updated outside the callback. auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); // Note Omega of gen is at state index 5! (Each added signal shifted by 1) // Bus -> 2 States diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp index 54e05ef5d..f7ef82e90 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp @@ -156,7 +156,7 @@ int main() // push it into output, which is updated outside the callback. auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); output.push_back(OutputData{t, static_cast(y_val[0]), diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp index 55d2ebd14..66aaf0a89 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp @@ -117,7 +117,7 @@ int main(int argc, const char* argv[]) // push it into output, which is updated outside the callback. auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + auto& y_val = sys.y(); output.push_back(OutputData{t, static_cast(y_val[0]), diff --git a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp index 15709c95f..e8032da66 100644 --- a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp +++ b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp @@ -2,6 +2,7 @@ #define _USE_MATH_DEFINES #include +#include #include #include #include @@ -73,9 +74,9 @@ int main(int /* argc */, char const** /* argv */) bus.allocate(); dg.allocate(); - dg.y() = t2; - dg.yp() = t1; - dg.getResidual() = res; + std::copy(t2.begin(), t2.end(), dg.y().data()); + std::copy(t1.begin(), t1.end(), dg.yp().data()); + std::copy(res.begin(), res.end(), dg.getResidual().data()); dg.setInternalPointer(&t2[dg.getExternSize()]); dg.setInternalDerivativePointer(&t1[dg.getExternSize()]); dg.setInternalResidualPointer(&dg.getResidual()[dg.getExternSize()]); diff --git a/examples/PowerElectronics/Microgrid/Microgrid.cpp b/examples/PowerElectronics/Microgrid/Microgrid.cpp index 035a8df0f..2fba1cc5d 100644 --- a/examples/PowerElectronics/Microgrid/Microgrid.cpp +++ b/examples/PowerElectronics/Microgrid/Microgrid.cpp @@ -195,7 +195,7 @@ int main(int /* argc */, char const** /* argv */) // Optional debuging output if (debug_output) { - std::vector& fres = sysmodel->getResidual(); + auto& fres = sysmodel->getResidual(); std::cout << "Verify initial resisdual is zero: {\n"; for (size_t i = 0; i < fres.size(); i++) { @@ -229,7 +229,7 @@ int main(int /* argc */, char const** /* argv */) idas->runSimulation(t_final); - std::vector& yfinial = sysmodel->y(); + auto& yfinial = sysmodel->y(); // Optional debugging output if (debug_output) diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index f6ada5c01..720f2bb39 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -79,9 +79,10 @@ int main(int /* argc */, char const** /* argv */) sysmodel.evaluateResidual(); std::cout << "Verify initial resisdual is zero: {"; - for (double i : sysmodel.getResidual()) + auto& residual = sysmodel.getResidual(); + for (std::size_t i = 0; i < residual.size(); ++i) { - std::cout << i << ", "; + std::cout << residual.data()[i] << ", "; } std::cout << "}\n"; @@ -104,7 +105,7 @@ int main(int /* argc */, char const** /* argv */) idas.runSimulation(t_final); - std::vector& yfinial = sysmodel.y(); + auto& yfinial = sysmodel.y(); std::cout << "Final vector y\n"; for (size_t i = 0; i < yfinial.size(); i++) diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp index a8f0aaae2..d7dca615e 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -258,7 +258,7 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) sys_model->initialize(); sys_model->evaluateResidual(); - std::vector& fres = sys_model->getResidual(); + auto& fres = sys_model->getResidual(); if (debug_output) { std::cout << "Verify initial resisdual is zero: {\n"; @@ -290,7 +290,7 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) idas->runSimulation(t_final); - std::vector& yfinal = sys_model->y(); + auto& yfinal = sys_model->y(); if (debug_output) { diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp index 8cb36bcd7..83bcc507d 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp @@ -108,6 +108,46 @@ namespace GridKit return status.report(__func__); } + /** + * @brief Test non-owning Vector storage over host vector data. + * + * @param[in] N Number of elements in the vector. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome externalData(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + x.allocate(memory::HOST); + x.setToConst(ScalarT{}, memory::HOST); + + Vector alias(N); + alias.setData(x.data()); + + status *= (alias.size() == static_cast(N)); + status *= (!alias.empty()); + status *= (alias.data() == x.data()); + status *= (alias.size() == static_cast(N)); + + for (IdxT i = 0; i < N; ++i) + { + alias[static_cast(i)] = static_cast(i); + } + + const auto& const_alias = alias; + status *= (const_alias.size() == alias.size()); + status *= (const_alias.data() == alias.data()); + for (IdxT i = 0; i < N; ++i) + { + status *= isEqual(x[static_cast(i)], static_cast(i)); + status *= isEqual(const_alias[static_cast(i)], static_cast(i)); + } + + return status.report(__func__); + } + /** * @brief Test setting data in a vector from array. * @@ -191,7 +231,7 @@ namespace GridKit z.allocate(memory::HOST); z.copyFromExternal(y, memspace_, memory::HOST); - const ScalarT* z_data = z.getData(memory::HOST); + const ScalarT* z_data = z.data(); if (z_data == nullptr) { @@ -381,11 +421,11 @@ namespace GridKit for (IdxT i = 0; i < x.getSize(); ++i) { // std::cout << x->getData("cpu")[i] << "\n"; - if (!isEqual(x.getData(memory::HOST)[i], answer)) + if (!isEqual(x.data()[i], answer)) { std::cout << std::setprecision(16); success = false; - std::cout << "Solution vector element x[" << i << "] = " << x.getData(memory::HOST)[i] + std::cout << "Solution vector element x[" << i << "] = " << x.data()[i] << ", expected: " << answer << "\n"; break; } diff --git a/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp b/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp index 53bacdb87..b9454a67f 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp +++ b/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp @@ -21,6 +21,7 @@ int main(int, char**) // result += test.copyToExternal(50); result += test.resize(100, 50); + result += test.externalData(50); result += test.setToConst(50); result += test.setToConst(50); diff --git a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp index 3e5ea1083..27eb4871e 100644 --- a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp @@ -64,7 +64,7 @@ namespace GridKit fault.allocate(); fault.initialize(); fault.evaluateResidual(); - std::vector res = fault.getResidual(); + auto& res = fault.getResidual(); for (size_t i = 0; i < res.size(); ++i) { @@ -135,7 +135,10 @@ namespace GridKit bus.evaluateResidual(); fault.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_y = fault.getResidual(); + auto& residual_y_view = fault.getResidual(); + std::vector residual_y( + residual_y_view.data(), + residual_y_view.data() + residual_y_view.size()); // Get d/dy' bus.initialize(); @@ -149,7 +152,10 @@ namespace GridKit bus.evaluateResidual(); fault.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_yp = fault.getResidual(); + auto& residual_yp_view = fault.getResidual(); + std::vector residual_yp( + residual_yp_view.data(), + residual_yp_view.data() + residual_yp_view.size()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp index fa440b3b4..d922ceea0 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp @@ -128,7 +128,8 @@ namespace GridKit bus.evaluateResidual(); exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_y = exciter.getResidual(); + auto& residual_y_view = exciter.getResidual(); + std::vector residual_y(residual_y_view.data(), residual_y_view.data() + residual_y_view.size()); // Get d/dy' bus.initialize(); @@ -142,7 +143,8 @@ namespace GridKit bus.evaluateResidual(); exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_yp = exciter.getResidual(); + auto& residual_yp_view = exciter.getResidual(); + std::vector residual_yp(residual_yp_view.data(), residual_yp_view.data() + residual_yp_view.size()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp index d6e2d7369..b74ccefe4 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp @@ -345,7 +345,8 @@ namespace GridKit bus.evaluateResidual(); exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_y = exciter.getResidual(); + auto& residual_y_view = exciter.getResidual(); + std::vector residual_y(residual_y_view.data(), residual_y_view.data() + residual_y_view.size()); // Get d/dy' bus.initialize(); @@ -359,7 +360,8 @@ namespace GridKit bus.evaluateResidual(); exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_yp = exciter.getResidual(); + auto& residual_yp_view = exciter.getResidual(); + std::vector residual_yp(residual_yp_view.data(), residual_yp_view.data() + residual_yp_view.size()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index 73b9cd3dc..77773e5c7 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -130,7 +130,7 @@ namespace GridKit gen.yp()[4] = 0; gen.evaluateResidual(); - std::vector& residual = gen.getResidual(); + auto& residual = gen.getResidual(); for (size_t i = 0; i < res_answer.size(); ++i) { @@ -277,7 +277,7 @@ namespace GridKit gen.allocate(); gen.initialize(); gen.evaluateResidual(); - std::vector res = gen.getResidual(); + auto& res = gen.getResidual(); for (size_t i = 0; i < res.size(); ++i) { @@ -351,7 +351,8 @@ namespace GridKit bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_y = gen.getResidual(); + auto& residual_y_view = gen.getResidual(); + std::vector residual_y(residual_y_view.data(), residual_y_view.data() + residual_y_view.size()); // Get d/dy' bus.initialize(); @@ -365,7 +366,8 @@ namespace GridKit bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_yp = gen.getResidual(); + auto& residual_yp_view = gen.getResidual(); + std::vector residual_yp(residual_yp_view.data(), residual_yp_view.data() + residual_yp_view.size()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 84a888b7c..2dfe9d176 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -127,10 +127,10 @@ namespace GridKit // Require results to be within machine precision auto tol = 10 * std::numeric_limits::epsilon(); - const std::vector& f = gen.getResidual(); - for (const auto& f_val : f) + const auto& f = gen.getResidual(); + for (std::size_t i = 0; i < f.size(); ++i) { - if (!isEqual(f_val, 0.0, tol)) + if (!isEqual(f.data()[i], 0.0, tol)) success = false; } @@ -280,7 +280,7 @@ namespace GridKit gen.yp()[5] = 1; // Edp_dot gen.evaluateResidual(); - std::vector& residual = gen.getResidual(); + auto& residual = gen.getResidual(); for (size_t i = 0; i < res_answer.size(); ++i) { @@ -375,7 +375,8 @@ namespace GridKit bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_y = gen.getResidual(); + auto& residual_y_view = gen.getResidual(); + std::vector residual_y(residual_y_view.data(), residual_y_view.data() + residual_y_view.size()); // Get d/dy' bus.initialize(); @@ -389,7 +390,8 @@ namespace GridKit bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual_yp = gen.getResidual(); + auto& residual_yp_view = gen.getResidual(); + std::vector residual_yp(residual_yp_view.data(), residual_yp_view.data() + residual_yp_view.size()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/GensalTests.hpp b/tests/UnitTests/PhasorDynamics/GensalTests.hpp index a4b3982d8..ecf2756d8 100644 --- a/tests/UnitTests/PhasorDynamics/GensalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GensalTests.hpp @@ -99,10 +99,10 @@ namespace GridKit gen.initialize(); gen.evaluateResidual(); - const std::vector& f = gen.getResidual(); - for (const auto& f_val : f) + const auto& f = gen.getResidual(); + for (std::size_t i = 0; i < f.size(); ++i) { - if (!isEqual(f_val, 0.0, tol_)) + if (!isEqual(f.data()[i], 0.0, tol_)) { success = false; break; @@ -139,10 +139,10 @@ namespace GridKit gen.initialize(); gen.evaluateResidual(); - const std::vector& f = gen.getResidual(); - for (const auto& f_val : f) + const auto& f = gen.getResidual(); + for (std::size_t i = 0; i < f.size(); ++i) { - if (!isEqual(f_val, 0.0, tol_)) + if (!isEqual(f.data()[i], 0.0, tol_)) { success = false; break; @@ -307,7 +307,7 @@ namespace GridKit gen.yp()[4] = 0.9; // psiqpp_dot gen.evaluateResidual(); - std::vector& residual = gen.getResidual(); + auto& residual = gen.getResidual(); for (size_t i = 0; i < res_answer.size(); ++i) { @@ -370,7 +370,8 @@ namespace GridKit bus.evaluateResidual(); gen.evaluateResidual(); - std::vector residual_y = gen.getResidual(); + auto& residual_y_view = gen.getResidual(); + std::vector residual_y(residual_y_view.data(), residual_y_view.data() + residual_y_view.size()); bus.initialize(); gen.initialize(); @@ -382,7 +383,8 @@ namespace GridKit bus.evaluateResidual(); gen.evaluateResidual(); - std::vector residual_yp = gen.getResidual(); + auto& residual_yp_view = gen.getResidual(); + std::vector residual_yp(residual_yp_view.data(), residual_yp_view.data() + residual_yp_view.size()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 054c542e9..ee0fcb050 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -128,7 +128,7 @@ namespace GridKit gov.yp()[1] = 2.0; // pv_dot gov.evaluateResidual(); - std::vector& residual = gov.getResidual(); + auto& residual = gov.getResidual(); for (size_t i = 0; i < res_answer.size(); ++i) { @@ -213,10 +213,10 @@ namespace GridKit // Require results to be within machine precision auto tol = 10 * std::numeric_limits::epsilon(); - const std::vector& f = gov.getResidual(); - for (const auto& f_val : f) + const auto& f = gov.getResidual(); + for (std::size_t i = 0; i < f.size(); ++i) { - if (!isEqual(f_val, 0.0, tol)) + if (!isEqual(f.data()[i], 0.0, tol)) success = false; } @@ -317,7 +317,8 @@ namespace GridKit gen.evaluateResidual(); gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking // the dependencies - std::vector residual_y = gov.getResidual(); + auto& residual_y_view = gov.getResidual(); + std::vector residual_y(residual_y_view.data(), residual_y_view.data() + residual_y_view.size()); // Get d/dy' bus.initialize(); @@ -333,7 +334,8 @@ namespace GridKit gen.evaluateResidual(); gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking // the dependencies - std::vector residual_yp = gov.getResidual(); + auto& residual_yp_view = gov.getResidual(); + std::vector residual_yp(residual_yp_view.data(), residual_yp_view.data() + residual_yp_view.size()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp b/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp index bd3f2f688..c7e1640a3 100644 --- a/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp +++ b/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp @@ -191,8 +191,10 @@ namespace GridKit } bus.evaluateResidual(); - load.evaluateResidual(); - auto residual_y = load.getResidual(); + 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.data(), residual_y_view.data() + residual_y_view.size()); bus.initialize(); load.initialize(); @@ -203,8 +205,10 @@ namespace GridKit } bus.evaluateResidual(); - load.evaluateResidual(); - auto residual_yp = load.getResidual(); + 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.data(), residual_yp_view.data() + residual_yp_view.size()); std::vector dependencies(residual_y.size()); for (IdxT i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/LoadZTests.hpp b/tests/UnitTests/PhasorDynamics/LoadZTests.hpp index 5df7cd27f..8065b5f58 100644 --- a/tests/UnitTests/PhasorDynamics/LoadZTests.hpp +++ b/tests/UnitTests/PhasorDynamics/LoadZTests.hpp @@ -117,7 +117,7 @@ namespace GridKit load.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residuals = load.getResidual(); + auto& residuals = load.getResidual(); std::vector ref = analyticalJacobian(R, X); /// Compare dependencies computed automatically to the ones computed analytically diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 15b64876f..077fd4122 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -223,7 +223,8 @@ namespace GridKit setStatePointDep(stab); stab.evaluateResidual(); - std::vector residual_y = stab.getResidual(); + auto& residual_y_view = stab.getResidual(); + std::vector residual_y(residual_y_view.data(), residual_y_view.data() + residual_y_view.size()); // --- d/dy': tag derivatives as independent --- stab.initialize(); @@ -236,7 +237,8 @@ namespace GridKit setStatePointDep(stab); stab.evaluateResidual(); - std::vector residual_yp = stab.getResidual(); + auto& residual_yp_view = stab.getResidual(); + std::vector residual_yp(residual_yp_view.data(), residual_yp_view.data() + residual_yp_view.size()); // Print dependencies for debugging for (size_t i = 0; i < residual_y.size(); ++i) diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 94bd68412..5fc117493 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -1,8 +1,10 @@ #pragma once +#include #include #include #include +#include #include #include @@ -11,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -153,6 +156,66 @@ 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::LoadZ 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.y().data()[0], static_cast(9.0)); + success *= isEqual(load.y().data()[1], static_cast(10.0)); + + load.y().data()[0] = 11.0; + load.y().data()[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)); + + bus2.tag()[0] = 1.0; + success *= isEqual(system.tag()[bus2_vr_index], static_cast(1.0)); + + system.tag()[bus2_vr_index] = 0.0; + success *= isEqual(bus2.tag()[0], static_cast(0.0)); + + load.absoluteTolerance()[0] = 0.123; + success *= isEqual(system.absoluteTolerance()[load_ir_index], static_cast(0.123)); + + system.absoluteTolerance()[load_ii_index] = 0.456; + success *= isEqual(load.absoluteTolerance()[1], static_cast(0.456)); + + return success.report(__func__); + } + /** * @brief Test for exception when signals are incorrectly configured */ @@ -168,7 +231,6 @@ namespace GridKit status *= throws( [&]() { sys.allocate(); }); - return status.report(__func__); } @@ -239,7 +301,7 @@ namespace GridKit // Evaluate and get the system residuals system.evaluateResidual(); - std::vector residual = system.getResidual(); + auto& residual = system.getResidual(); // Print the dependencies for (size_t i = 0; i < residual.size(); ++i) 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 diff --git a/tests/UnitTests/Solver/Dynamic/CMakeLists.txt b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt index b6ed121e8..fb2de0fcf 100644 --- a/tests/UnitTests/Solver/Dynamic/CMakeLists.txt +++ b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt @@ -2,7 +2,10 @@ add_executable(test_ida runIdaTests.cpp) target_link_libraries( - test_ida GridKit::solvers_dyn GridKit::testing) + test_ida + GridKit::solvers_dyn + GridKit::dense_vector + GridKit::testing) # add_test(NAME IDATest COMMAND $) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index 3962d76db..1956556a1 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -14,8 +14,18 @@ namespace GridKit template class NullEvaluator : public Model::Evaluator { + protected: + using Model::Evaluator::y_; + using Model::Evaluator::yp_; + using Model::Evaluator::f_; + using Model::Evaluator::tag_; + using Model::Evaluator::abs_tol_; + using Model::Evaluator::allocated_; + using Model::Evaluator::allocateVectors; + public: - using RealT = typename Model::Evaluator::RealT; + using RealT = typename Model::Evaluator::RealT; + using VectorT = typename Model::Evaluator::VectorT; NullEvaluator() { @@ -23,19 +33,25 @@ namespace GridKit int allocate() override { + if (!allocated_) + { + allocateVectors(size()); + } return 0; } int initialize() override { - y_ = {0}; - yp_ = {0}; - - tag_ = {false}; - abs_tol_ = {0}; - - f_ = {0}; - g_ = {0}; + if (!allocated_) + { + allocate(); + } + + y_[0] = 0.0; + yp_[0] = 0.0; + tag_[0] = 0.0; + abs_tol_[0] = 0.0; + f_[0] = 0.0; return 0; } @@ -71,13 +87,13 @@ namespace GridKit int setAbsoluteTolerance(RealT rel_tol) override { - std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + std::fill(abs_tol_.data(), abs_tol_.data() + abs_tol_.size(), rel_tol); return 0; } int evaluateResidual() override { - f_ = y_; + f_[0] = y_[0]; return 0; } @@ -110,137 +126,97 @@ namespace GridKit { } - std::vector& y() override - { - return y_; - } - - const std::vector& y() const override - { - return y_; - } - - std::vector& yp() override - { - return yp_; - } - - const std::vector& yp() const override - { - return yp_; - } - - std::vector& tag() override - { - return tag_; - } - - const std::vector& tag() const override - { - return tag_; - } - - std::vector& absoluteTolerance() override - { - return abs_tol_; - } - - const std::vector& absoluteTolerance() const override - { - return abs_tol_; - } - - std::vector& yB() override + VectorT& yB() override { return yB_; } - const std::vector& yB() const override + const VectorT& yB() const override { return yB_; } - std::vector& ypB() override + VectorT& ypB() override { return ypB_; } - const std::vector& ypB() const override + const VectorT& ypB() const override { return ypB_; } - std::vector& param() override + VectorT& param() override { return param_; } - const std::vector& param() const override + const VectorT& param() const override { return param_; } - std::vector& param_up() override + VectorT& param_up() override { return param_up_; } - const std::vector& param_up() const override + const VectorT& param_up() const override { return param_up_; } - std::vector& param_lo() override + VectorT& param_lo() override { return param_lo_; } - const std::vector& param_lo() const override + const VectorT& param_lo() const override { return param_lo_; } - std::vector& getResidual() override + VectorT& getResidual() override { return f_; } - const std::vector& getResidual() const override + const VectorT& getResidual() const override { return f_; } GridKit::LinearAlgebra::CsrMatrix* getCsrJacobian() const override { - return csr_jac_; + return nullptr; } - std::vector& getIntegrand() override + VectorT& getIntegrand() override { return g_; } - const std::vector& getIntegrand() const override + const VectorT& getIntegrand() const override { return g_; } - std::vector& getAdjointResidual() override + VectorT& getAdjointResidual() override { return fB_; } - const std::vector& getAdjointResidual() const override + const VectorT& getAdjointResidual() const override { return fB_; } - std::vector& getAdjointIntegrand() override + VectorT& getAdjointIntegrand() override { return gB_; } - const std::vector& getAdjointIntegrand() const override + const VectorT& getAdjointIntegrand() const override { return gB_; } @@ -251,23 +227,16 @@ namespace GridKit } protected: - std::vector y_; - std::vector yp_; - std::vector tag_; - std::vector abs_tol_; - std::vector f_; - std::vector g_; - - std::vector yB_; - std::vector ypB_; - std::vector fB_; - std::vector gB_; - - GridKit::LinearAlgebra::CsrMatrix* csr_jac_; - - std::vector param_; - std::vector param_up_; - std::vector param_lo_; + VectorT g_; + + VectorT yB_; + VectorT ypB_; + VectorT fB_; + VectorT gB_; + + VectorT param_; + VectorT param_up_; + VectorT param_lo_; }; template From 3efa5b9ee699db2c1c4e4c9052f2383ed732cd66 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Fri, 26 Jun 2026 19:16:51 +0000 Subject: [PATCH 2/2] Apply pre-commit fixes --- GridKit/Model/PhasorDynamics/BusBase.hpp | 2 +- GridKit/Model/PhasorDynamics/Component.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index bcee3ae8e..5795765fc 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -36,7 +36,7 @@ namespace GridKit using CooMatrixT = typename Model::Evaluator::CooMatrixT; using BusTypeT = typename BusData::BusType; using MonitorT = Model::VariableMonitor; - using VectorT = typename Model::Evaluator::VectorT; + using VectorT = typename Model::Evaluator::VectorT; BusBase() = default; diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index c44a791db..c0716f2cd 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -28,7 +28,7 @@ namespace GridKit using RealT = typename Model::Evaluator::RealT; using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; using CooMatrixT = typename Model::Evaluator::CooMatrixT; - using VectorT = typename Model::Evaluator::VectorT; + using VectorT = typename Model::Evaluator::VectorT; Component() = default;