diff --git a/GridKit/Model/PartitionEvaluator.hpp b/GridKit/Model/PartitionEvaluator.hpp new file mode 100644 index 000000000..1ce80b513 --- /dev/null +++ b/GridKit/Model/PartitionEvaluator.hpp @@ -0,0 +1,282 @@ +#pragma once + +#include + +#include "Evaluator.hpp" + +namespace GridKit +{ + namespace Model + { + + template + class PartitionEvaluator : public Evaluator + { + public: + using Base = Evaluator; + using RealT = typename Base::RealT; + using MatrixT = typename Base::MatrixT; + using CsrMatrixT = typename Base::CsrMatrixT; + + using ForcingFunc = std::function; + + private: + std::vector y_; + std::vector yp_; + std::vector yB_; + std::vector ypB_; + std::vector param_; + std::vector param_up_; + std::vector param_lo_; + std::vector residual_; + std::vector integrand_; + std::vector adj_integrand_; + std::vector adj_residual_; + std::vector tag_; + + MatrixT jacobian_; + + std::vector coupling_; + std::vector forcing_; + + public: + int allocate() override + { + return 0; + } + + int initialize() override + { + return 0; + } + + int tagDifferentiable() override + { + return 0; + } + + int evaluateResidual() override + { + return 0; + } + + int evaluateJacobian() override + { + return 0; + } + + int evaluateIntegrand() override + { + return 0; + } + + int initializeAdjoint() override + { + return 0; + } + + int evaluateAdjointResidual() override + { + return 0; + } + + int evaluateAdjointIntegrand() override + { + return 0; + } + + IdxT size() override + { + return 0; + } + + IdxT nnz() override + { + return 0; + } + + bool hasJacobian() override + { + return false; + } + + IdxT sizeQuadrature() override + { + return 0; + } + + IdxT sizeParams() override + { + return 0; + } + + void updateTime(RealT, RealT) override + { + } + + void setTolerances(RealT&, RealT&) const override + { + } + + void setMaxSteps(IdxT&) const override + { + } + + 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& yB() override + { + return yB_; + } + + const std::vector& yB() const override + { + return yB_; + } + + std::vector& ypB() override + { + return ypB_; + } + + const std::vector& ypB() const override + { + return ypB_; + } + + std::vector& param() override + { + return param_; + } + + const std::vector& param() const override + { + return param_; + } + + std::vector& param_up() override + { + return param_up_; + } + + const std::vector& param_up() const override + { + return param_up_; + } + + std::vector& param_lo() override + { + return param_lo_; + } + + const std::vector& param_lo() const override + { + return param_lo_; + } + + std::vector& getResidual() override + { + return residual_; + } + + const std::vector& getResidual() const override + { + return residual_; + } + + MatrixT& getJacobian() override + { + return jacobian_; + } + + const MatrixT& getJacobian() const override + { + return jacobian_; + } + + std::vector& getIntegrand() override + { + return integrand_; + } + + const std::vector& getIntegrand() const override + { + return integrand_; + } + + std::vector& getAdjointResidual() override + { + return adj_residual_; + } + + const std::vector& getAdjointResidual() const override + { + return adj_residual_; + } + + std::vector& getAdjointIntegrand() override + { + return adj_integrand_; + } + + const std::vector& getAdjointIntegrand() const override + { + return adj_integrand_; + } + + void addForcing(const ForcingFunc& f) + { + forcing_.push_back(f); + } + + IdxT getCouplingSize() const + { + return static_cast(coupling_.size()); + } + + IdxT getStateSize() const + { + return static_cast(y_.size()); + } + + std::vector& getCoupling() + { + return coupling_; + } + + const std::vector& getCoupling() const + { + return coupling_; + } + }; + + } // namespace Model +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/CMakeLists.txt b/GridKit/Model/PowerElectronics/CMakeLists.txt index 98dce1148..efe1c7f68 100644 --- a/GridKit/Model/PowerElectronics/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/CMakeLists.txt @@ -21,6 +21,7 @@ add_subdirectory(TransmissionLine) add_subdirectory(MicrogridLoad) add_subdirectory(MicrogridLine) add_subdirectory(MicrogridBusDQ) +add_subdirectory(PartitionInterface) install( FILES CircuitComponent.hpp @@ -28,4 +29,5 @@ install( CircuitGraph.hpp SystemModelPowerElectronics.hpp NodeBase.hpp + SubsystemModel.hpp DESTINATION include/GridKit/Model/PowerElectronics) diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index 026db4df8..c7143d446 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -26,6 +26,128 @@ namespace GridKit CircuitComponent() = default; + CircuitComponent(const CircuitComponent& other) + : n_extern_(other.n_extern_), + n_intern_(other.n_intern_), + extern_indices_(other.extern_indices_), + size_(other.size_), + nnz_(other.nnz_), + size_quad_(other.size_quad_), + size_opt_(other.size_opt_), + current_jac_size_(other.current_jac_size_), + + // These pointers refer to storage supplied by a parent system. + // The copied component must be connected to its own storage later. + y_int_(nullptr), + yp_int_(nullptr), + f_int_(nullptr), + + tag_(other.tag_), + time_(other.time_), + alpha_(other.alpha_), + max_steps_(other.max_steps_), + idc_(other.idc_), + allocated_(other.allocated_) + { + /* + * VectorT disables its normal copy constructor and copy-assignment + * operator. Use its provided copyFromExternal() operation to perform + * an independent copy of the vector data. + */ + auto copyVector = [](VectorT& destination, const VectorT& source) + { + const IdxT source_size = source.getSize(); + + if (source_size == 0) + { + return; + } + + destination.resize(source_size); + destination.copyFromExternal(source); + }; + + /* + * Deep-copy the local-to-global connection mapping. + */ + if (other.connection_nodes_) + { + connection_nodes_ = + std::make_unique(static_cast(size_)); + + for (IdxT i = 0; i < size_; ++i) + { + connection_nodes_[static_cast(i)] = + other.connection_nodes_[static_cast(i)]; + } + } + + /* + * Deep-copy the COO Jacobian row indices. + */ + if (other.jacobian_coo_rows_) + { + jacobian_coo_rows_ = + std::make_unique(static_cast(nnz_)); + + for (IdxT i = 0; i < nnz_; ++i) + { + jacobian_coo_rows_[static_cast(i)] = + other.jacobian_coo_rows_[static_cast(i)]; + } + } + + /* + * Deep-copy the COO Jacobian column indices. + */ + if (other.jacobian_coo_cols_) + { + jacobian_coo_cols_ = + std::make_unique(static_cast(nnz_)); + + for (IdxT i = 0; i < nnz_; ++i) + { + jacobian_coo_cols_[static_cast(i)] = + other.jacobian_coo_cols_[static_cast(i)]; + } + } + + /* + * Deep-copy the COO Jacobian values. + */ + if (other.jacobian_coo_values_) + { + jacobian_coo_values_ = + std::make_unique(static_cast(nnz_)); + + for (IdxT i = 0; i < nnz_; ++i) + { + jacobian_coo_values_[static_cast(i)] = + other.jacobian_coo_values_[static_cast(i)]; + } + } + + // State, state derivative, residual, and absolute tolerance. + copyVector(y_, other.y_); + copyVector(yp_, other.yp_); + copyVector(f_, other.f_); + copyVector(abs_tol_, other.abs_tol_); + + // Quadrature vector. + copyVector(g_, other.g_); + + // Adjoint vectors. + copyVector(yB_, other.yB_); + copyVector(ypB_, other.ypB_); + copyVector(fB_, other.fB_); + copyVector(gB_, other.gB_); + + // Parameter vectors. + copyVector(param_, other.param_); + copyVector(param_up_, other.param_up_); + copyVector(param_lo_, other.param_lo_); + } + /** * @note Cannot be marked final, since it is overriden to recurse in the system model. */ @@ -50,7 +172,7 @@ namespace GridKit return this->n_intern_; } - std::set getExternIndices() + std::set getExternIndices() { return this->extern_indices_; } diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp new file mode 100644 index 000000000..578493931 --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp @@ -0,0 +1,271 @@ + +#include "BusPartitionInterface.hpp" + +#include +#include +#include +#include + +namespace GridKit +{ + + /*! + * @brief Construct a bus-level partition interface. + * + * This interface exposes only the two bus variables associated with a + * component that has been separated from the rest of the system during + * partitioning. Internally, it maintains a copy of the original component + * and acts as a proxy between the co-simulation algorithm and the component + * model. + * + * The bus variables become the internal variables of this interface, + * while all remaining component variables are treated as external data + * supplied by neighboring partitions. + * + * @param component Circuit component associated with this interface. + * @param id Unique component identifier. + */ + template + BusPartitionInterface::BusPartitionInterface(node_type& bus, component_type& component, IdxT id) + : component_(component), + bus_(bus) + { + size_ = component_.size(); + n_intern_ = 0; + n_extern_ = static_cast(component_.size()); + idc_ = id; + + for (IdxT i = 0; i < size_; i++) + { + extern_indices_.insert(i); + } + + const IdxT* cooRow = component_.jacobianCooRows(); + const IdxT* cooCol = component_.jacobianCooCols(); + + const IdxT bus_i = bus.getNodeConnection(0); + const IdxT bus_j = bus.getNodeConnection(1); + + nnz_ = 0; + for (IdxT k = 0; k < component_.nnz(); ++k) + { + const IdxT row_node = component_.getNodeConnection(cooRow[k]); + const IdxT col_node = component_.getNodeConnection(cooCol[k]); + + const bool row_is_bus = (row_node == bus_i || row_node == bus_j); + const bool col_is_bus = (col_node == bus_i || col_node == bus_j); + + if (row_is_bus && col_is_bus) + { + ++nnz_; + jac_map_.push_back(k); + } + } + } + + template + BusPartitionInterface::~BusPartitionInterface() + { + } + + /*! + * @brief allocate method + */ + template + int BusPartitionInterface::allocate() + { + + CircuitComponent::allocate(); + + std::fill_n(f_.getData(), size_, 0); + + bool port_i_set = false; + bool port_j_set = false; + + for (size_t i = 0; i < static_cast(size_); i++) + { + if (bus_.getNodeConnection(0) == component_.getNodeConnection(static_cast(i))) + { + bus_port_i_ = i; + port_i_set = true; + } + else if (bus_.getNodeConnection(1) == component_.getNodeConnection(static_cast(i))) + { + bus_port_j_ = i; + port_j_set = true; + } + IdxT global_idx = component_.getNodeConnection(static_cast(i)); + this->setExternalConnectionNodes(static_cast(i), global_idx); + } + + if (!port_i_set || !port_j_set) + { + std::cerr << "ERROR: Invalid partition interface detected. " + << "Bus(ID=" << bus_.busID() + << "), Component(ID=" << component_.getIDcomponent() + << "). Please verify connection-node mappings and internal/external index assignments." + << std::endl; + return 1; + } + + const auto n = component_.getInternalSize(); + + y_ptr = std::make_unique(n); + yp_ptr = std::make_unique(n); + f_ptr = std::make_unique(n); + + component_.setInternalPointer(y_ptr.get()); + component_.setInternalDerivativePointer(yp_ptr.get()); + component_.setInternalResidualPointer(f_ptr.get()); + + return 0; + } + + template + int BusPartitionInterface::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * Initialization of the grid model + */ + template + int BusPartitionInterface::initialize() + { + return 0; + } + + /* + * \brief Identify differential variables + */ + template + int BusPartitionInterface::tagDifferentiable() + { + return 0; + } + + /** + * @brief Eval Internal Residual + */ + template + int BusPartitionInterface::evaluateInternalResidual() + { + return 0; + } + + template + int BusPartitionInterface::evaluateExternalResidual() + { + size_t internal = 0; + size_t external = 0; + + ScalarT* y = y_.getData(); + ScalarT* yp = yp_.getData(); + ScalarT* f = f_.getData(); + + ScalarT* component_y = component_.y().getData(); + ScalarT* component_yp = component_.yp().getData(); + + const auto& extern_indices = component_.getExternIndices(); + + for (size_t i = 0; i < static_cast(component_.size()); ++i) + { + if (extern_indices.contains(static_cast(i))) + { + component_y[external] = y[i]; + component_yp[external] = yp[i]; + ++external; + } + else + { + y_ptr[internal] = y[i]; + yp_ptr[internal] = yp[i]; + ++internal; + } + } + + component_.y().setDataUpdated(); + component_.y().setDataUpdated(); + + component_.evaluateResidual(); + + const auto* residual = component_.getResidual().getData(); + + // TODO: This assumes that external variables are ordered after all internal + // variables in the indexing. To make this more robust, we need to get rid of this assumption + // (although true for all components that we have in GridKit). + f[bus_port_i_] = residual[bus_port_i_]; + f[bus_port_j_] = residual[bus_port_j_]; + + f_.setDataUpdated(); + + return 0; + } + + /** + * @brief Generate Jacobian for + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ + template + int BusPartitionInterface::evaluateJacobian() + { + this->zeroJacMatrix(); + + component_.evaluateJacobian(); + + const IdxT* cooRows = component_.jacobianCooRows(); + const IdxT* cooCols = component_.jacobianCooCols(); + const RealT* cooVals = component_.jacobianCooValues(); + + std::vector r = {}; + std::vector c = {}; + std::vector v = {}; + + for (const auto& index : jac_map_) + { + r.push_back(cooRows[index]); + c.push_back(cooCols[index]); + v.push_back(cooVals[index]); + } + + this->setJacValues(r, c, v); + + return 0; + } + + template + int BusPartitionInterface::evaluateIntegrand() + { + return 0; + } + + template + int BusPartitionInterface::initializeAdjoint() + { + return 0; + } + + template + int BusPartitionInterface::evaluateAdjointResidual() + { + return 0; + } + + template + int BusPartitionInterface::evaluateAdjointIntegrand() + { + return 0; + } + + // Available template instantiations + template class BusPartitionInterface; + template class BusPartitionInterface; + template class BusPartitionInterface; + template class BusPartitionInterface; + +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp new file mode 100644 index 000000000..93c764b8f --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp @@ -0,0 +1,79 @@ + + +#pragma once + +#include +#include + +#include "GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp" + +namespace GridKit +{ + /*! + * @brief Declaration of a passive BusPartitionInterface class. + * + */ + template + class BusPartitionInterface : public PartitionInterface + { + + using component_type = CircuitComponent; + using node_type = typename PowerElectronics::NodeBase; + using RealT = typename CircuitComponent::RealT; + + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_; + using CircuitComponent::yp_int_; + using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_; + using CircuitComponent::f_int_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + using PartitionInterface::bus_port_i_; + using PartitionInterface::bus_port_j_; + + public: + BusPartitionInterface(node_type& bus, component_type& component, IdxT id); + virtual ~BusPartitionInterface(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int setAbsoluteTolerance(RealT); + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + component_type& component_; + node_type& bus_; + + std::unique_ptr y_ptr; + std::unique_ptr yp_ptr; + std::unique_ptr f_ptr; + + std::vector jac_map_; + }; +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt new file mode 100644 index 000000000..b339bfa2f --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt @@ -0,0 +1,4 @@ +gridkit_add_library( + power_elec_partition_interfaces + SOURCES BusPartitionInterface.cpp + HEADERS BusPartitionInterface.hpp PartitionInterface.hpp) diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp b/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp new file mode 100644 index 000000000..f639b0448 --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp @@ -0,0 +1,37 @@ + + +#pragma once + +#include +#include +#include + +namespace GridKit +{ + /*! + * @brief Base class for partition interface components. + * + * PartitionInterface provides the infrastructure needed by components + * that exchange information across partition boundaries during + * co-simulation. It stores the values of external variables received + * from neighboring partitions together with the corresponding global + * indices in the original unpartitioned system. + * + * Derived classes are responsible for implementing the specific + * interface behavior associated with buses, circuit components, or + * other partition boundary entities. + */ + template + class PartitionInterface : public CircuitComponent + { + public: + PartitionInterface() = default; + + ~PartitionInterface() = default; + + protected: + size_t bus_port_i_; + size_t bus_port_j_; + }; + +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/SubsystemModel.hpp b/GridKit/Model/PowerElectronics/SubsystemModel.hpp new file mode 100644 index 000000000..0c0cbcab5 --- /dev/null +++ b/GridKit/Model/PowerElectronics/SubsystemModel.hpp @@ -0,0 +1,657 @@ + + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + + template + class SubsystemModel : public PowerElectronicsModel + { + public: + using ForcingData = std::tuple, std::vector>; + using TimeFunction = std::function; + + protected: + using SystemModel = PowerElectronicsModel; + using RealT = typename CircuitComponent::RealT; + using CsrMatrixT = typename CircuitComponent::CsrMatrixT; + using component_type = CircuitComponent; + using node_type = PowerElectronics::NodeBase; + + using SystemModel::abs_tol_; + using SystemModel::allocated_; + using SystemModel::allocateVectors; + using SystemModel::alpha_; + using SystemModel::f_; + using SystemModel::f_int_; + using SystemModel::n_extern_; + using SystemModel::n_intern_; + using SystemModel::nnz_; + using SystemModel::size_; + using SystemModel::tag_; + using SystemModel::time_; + using SystemModel::y_; + using SystemModel::y_int_; + using SystemModel::yp_; + using SystemModel::yp_int_; + + using SystemModel::components_; + using SystemModel::csr_jac_; + using SystemModel::jac_call_count_; + using SystemModel::map_to_csr_; + using SystemModel::neg1_; + using SystemModel::nodes_; + using SystemModel::use_jac_; + + public: + /** + * @brief Default constructor for the system model + * + * @post System model parameters set as default + */ + + SubsystemModel(bool use_jac = false) + : SystemModel(use_jac) + { + } + + ~SubsystemModel() override + { + // SubsystemModel does not own components_/nodes_ — they belong to (and are + // deleted by) the parent system model. + components_.clear(); + nodes_.clear(); + } + + int allocate() override + { + + hold(); + + n_intern_ = internal_map_.size(); + n_extern_ = external_map_.size(); + size_ = n_intern_; + + CircuitComponent::allocate(); + + // Allocation always rebuilds the system Jacobian and its COO-to-CSR map. + delete csr_jac_; + csr_jac_ = nullptr; + + delete[] map_to_csr_; + map_to_csr_ = nullptr; + + tag_.resize(size_); + + // Allocate subsystem vectors. + y_ext_.resize(n_extern_); + yp_ext_.resize(n_extern_); + f_ext_.resize(n_extern_); + + external_data_indices_.resize(n_extern_); + + // Store the mapping from local subsystem indices back to their global system indices + for (const auto [global_idx, local_idx] : internal_map_) + { + this->setExternalConnectionNodes(local_idx, global_idx); + } + + // Store the global indices of all external coupling variables. + for (const auto [global_idx, local_idx] : external_map_) + { + external_data_indices_[local_idx - n_intern_] = global_idx; + } + + size_t component_internal_idx = 0; + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + auto* f = f_.getData(); + for (component_type* comp : components_) + { + + comp->setInternalPointer(&y[component_internal_idx]); + comp->setInternalDerivativePointer(&yp[component_internal_idx]); + comp->setInternalResidualPointer(&f[component_internal_idx]); + + component_internal_idx += comp->getInternalSize(); + } + + // Evaluate component Jacobians to get sparsity + distributeVectors(); + for (component_type* component : components_) + { + component->evaluateJacobian(); + } + + auto isValidEntry = [this](IdxT row, IdxT col) + { + if (row == neg1_ || col == neg1_) + { + return false; + } + + const bool row_is_internal = row < this->getInternalSize(); + const bool col_is_internal = col < this->getInternalSize(); + + return (row_is_internal && col_is_internal); + }; + + IdxT nnz_dup = 0; + + for (const component_type* component : components_) + { + const IdxT* r = component->jacobianCooRows(); + const IdxT* c = component->jacobianCooCols(); + const IdxT nnz = component->nnz(); + + for (IdxT i = 0; i < nnz; ++i) + { + const IdxT row = component->getNodeConnection(r[i]); + const IdxT col = component->getNodeConnection(c[i]); + + if (isValidEntry(row, col)) + { + ++nnz_dup; + } + } + } + + // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix) + IdxT* rows_dup = new IdxT[nnz_dup]; + IdxT* cols_dup = new IdxT[nnz_dup]; + RealT* vals_dup = new RealT[nnz_dup]; + + IdxT counter = 0; + + for (const component_type* component : components_) + { + const IdxT* r = component->jacobianCooRows(); + const IdxT* c = component->jacobianCooCols(); + const RealT* v = component->jacobianCooValues(); + const IdxT nnz = component->nnz(); + + for (IdxT i = 0; i < nnz; ++i) + { + const IdxT row = component->getNodeConnection(r[i]); + const IdxT col = component->getNodeConnection(c[i]); + + if (!isValidEntry(row, col)) + { + continue; + } + + rows_dup[counter] = row; + cols_dup[counter] = col; + vals_dup[counter] = v[i]; + + ++counter; + } + } + + // Build the system COO Jacobian + LinearAlgebra::CooMatrix jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); + + // Populate CSR data with sort and deduplicate + IdxT* row_ptrs = jac.getCsrRowData(); + + // Deduplicated nnz + nnz_ = jac.getNnz(); + + // Allocate cols/vals with deduplicated nnz + IdxT* cols = new IdxT[nnz_]; + RealT* vals = new RealT[nnz_]; + + std::copy(jac.getColData(), jac.getColData() + nnz_, cols); + std::copy(jac.getValues(), jac.getValues() + nnz_, vals); + + // Create the CSR Jacobian + csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); + + const IdxT* map_to_sorted = jac.getMapToSorted(); + const IdxT* map_to_dedup = jac.getMapToDeduplicated(); + + // Build a mappping from original COO index to CSR index + map_to_csr_ = new IdxT[nnz_dup]; + for (IdxT i = 0; i < nnz_dup; ++i) + { + map_to_csr_[map_to_sorted[i]] = map_to_dedup[i]; + } + + allocated_ = true; + return 0; + } + + /** + * @brief Distribute y and y' to each component based of node connection graph + * + * @post Each component has y and y' set + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int distributeVectors() override + { + + if (forcing_function_) + { + const auto [y_forcing, yp_forcing] = (*forcing_function_)(time_); + + assert(y_forcing.size() == y_ext_.size()); + assert(yp_forcing.size() == yp_ext_.size()); + + std::copy(y_forcing.begin(), y_forcing.end(), y_ext_.begin()); + std::copy(yp_forcing.begin(), yp_forcing.end(), yp_ext_.begin()); + } + + const auto* y_system = y_.getData(); + const auto* yp_system = yp_.getData(); + + for (component_type* component : components_) + { + auto* y = component->y().getData(); + auto* yp = component->yp().getData(); + const auto& externals = component->getExternIndices(); + + for (size_t j : externals) + { + IdxT local_idx = component->getNodeConnection(j); + + if (local_idx == neg1_) + { + y[j] = 0.0; + yp[j] = 0.0; + } + else if (local_idx < this->getInternalSize()) + { + y[j] = y_system[local_idx]; + yp[j] = yp_system[local_idx]; + } + else + { + y[j] = y_ext_[local_idx - this->getInternalSize()]; + yp[j] = yp_ext_[local_idx - this->getInternalSize()]; + } + } + component->y().setDataUpdated(); + component->yp().setDataUpdated(); + } + return 0; + } + + /** + * @brief Add a component to the subsystem. + * + * Rejected while connections are in the local-indexed state, since the + * component's stored connection indices would otherwise be interpreted + * inconsistently with the rest of the subsystem. Call release() first. + * + * @param[in] component Component to add. + */ + void addComponent(component_type* component) + { + if (is_comp_in_local_state) + { + std::cerr << "SubsystemModel::addComponent: cannot add component while " + "in local-indexed state. Call release() first.\n"; + return; + } + + SystemModel::addComponent(component); + } + + /** + * @brief Add a node to the subsystem. + * + * Rejected while connections are in the local-indexed state, since the + * node's stored connection indices would otherwise be interpreted + * inconsistently with the rest of the subsystem. Call release() first. + * + * @param[in] node Node to add. + */ + void addNode(node_type* node) + { + if (is_comp_in_local_state) + { + std::cerr << "SubsystemModel::addNode: cannot add node while in " + "local-indexed state. Call release() first.\n"; + return; + } + + SystemModel::addNode(node); + } + + /** + * @brief Acquire the subsystem's global-to-local index mappings. + * + * Builds the internal/external index maps via buildIndexMappings(), + * then converts component/node connections from global to local + * indices via mapGlobalToLocal(). + * + * Call this after modifying subsystem topology (adding or removing + * components/nodes) and before allocate(), to (re)establish local + * indexing. Pairs with release(), which undoes this. + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int hold() + { + buildIndexMappings(); + + if (int err_code = mapGlobalToLocal()) + { + return err_code; + } + + return 0; + } + + /** + * @brief Release the subsystem's global-to-local index mappings. + * + * If components/nodes currently hold local indices, they are first + * converted back to their original global indices via + * mapLocalToGlobal(). The internal/external index maps built by + * buildIndexMappings() are then cleared, and the model is marked as + * not allocated. + * + * + * Call this before modifying subsystem topology (adding or removing + * components/nodes), so stale local indices aren't left dangling. + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int release() + { + if (int err_code = mapLocalToGlobal()) + { + return err_code; + } + + internal_map_.clear(); + external_map_.clear(); + + allocated_ = false; + + return 0; + } + + const std::vector& getExternalDataIndices() + { + return external_data_indices_; + } + + std::vector& getExternalDataY() + { + return y_ext_; + } + + std::vector& getExternalDataYP() + { + return yp_ext_; + } + + std::vector& getExternalDataF() + { + return f_ext_; + } + + void setTimeFunction(TimeFunction function) + { + forcing_function_ = std::move(function); + } + + const std::unordered_map& getInternalMap() const + { + return internal_map_; + } + + const std::unordered_map& getExternalMap() const + { + return external_map_; + } + + private: + /** + * @brief Replace local subsystem connection indices with their original global indices + * for every component and node. + * + * Inverse of mapGlobalToLocal(). Internal connections (index < internal + * size) are mapped back through this subsystem's own node-connection + * table; external connections are mapped back through + * @c external_data_indices_. Entries equal to @c neg1_ represent no + * connection and are left unchanged. + * + * @return Always returns 0. + */ + int mapLocalToGlobal() + { + if (!is_comp_in_local_state) + { + return 0; + } + + for (auto* component : components_) + { + + for (IdxT i = 0; i < component->size(); i++) + { + IdxT index = component->getNodeConnection(i); + + if (index == neg1_) + { + continue; + } + else if (index < this->getInternalSize()) + { + component->setExternalConnectionNodes(i, this->getNodeConnection(index)); + } + else + { + component->setExternalConnectionNodes(i, external_data_indices_[index - this->getInternalSize()]); + } + } + } + + for (auto* node : nodes_) + { + + for (IdxT i = 0; i < node->size(); i++) + { + IdxT index = node->getNodeConnection(i); + + if (index == neg1_) + { + continue; + } + else if (index < this->getInternalSize()) + { + node->setExternalConnectionNodes(i, this->getNodeConnection(index)); + } + else + { + node->setExternalConnectionNodes(i, external_data_indices_[index - this->getInternalSize()]); + } + } + } + + is_comp_in_local_state = false; + + return 0; + } + + /** + * @brief Replace global connection indices with local subsystem indices for every component + * and node. + * + * After the internal and external maps have been constructed, each component + * and node still stores the original global connection indices. This function + * traverses every connection and replaces each global index with its + * corresponding local subsystem index. + * + * Internal connections are mapped using @c internal_map_, while connections + * that reference variables outside the subsystem are mapped using + * @c external_map_. Entries equal to @c neg1_ represent no connection and + * are left unchanged. + * + * @return Always returns 0. + */ + + int mapGlobalToLocal() + { + + if (is_comp_in_local_state) + { + return 0; + } + + for (auto* component : components_) + { + + for (IdxT i = 0; i < component->size(); i++) + { + IdxT index = component->getNodeConnection(i); + + if (index == neg1_) + { + continue; + } + else if (internal_map_.count(index) > 0) + { + component->setExternalConnectionNodes(i, internal_map_.at(index)); + } + else + { + component->setExternalConnectionNodes(i, external_map_.at(index)); + } + } + } + + for (auto* node : nodes_) + { + + for (IdxT i = 0; i < node->size(); i++) + { + IdxT index = node->getNodeConnection(i); + + if (index == neg1_) + { + continue; + } + else if (internal_map_.count(index) > 0) + { + node->setExternalConnectionNodes(i, internal_map_.at(index)); + } + else + { + node->setExternalConnectionNodes(i, external_map_.at(index)); + } + } + } + + is_comp_in_local_state = true; + + return 0; + } + + /** + * @brief Construct the mappings from global variable indices to local + * subsystem indices. + * + * The subsystem stores its variables in a contiguous local ordering. This + * function assigns each global variable index a corresponding local index and + * separates them into internal and external variables. + * + * Internal variables are those owned by this subsystem, while external + * variables are owned by neighboring subsystems but are required for residual + * and Jacobian evaluation. + */ + void buildIndexMappings() + { + + if (is_comp_in_local_state) + { + return; + } + + internal_map_.clear(); + external_map_.clear(); + + size_t component_internal_idx = 0; + // First pass: Map global component indices to local subsystem indices. + for (component_type* comp : components_) + { + for (IdxT i = 0; i < comp->size(); i++) + { + IdxT index = comp->getNodeConnection(i); + + auto extern_indices = comp->getExternIndices(); + + if (index != neg1_ && !extern_indices.contains(i)) + { + internal_map_[index] = component_internal_idx++; + } + } + } + + for (node_type* node : nodes_) + { + + for (IdxT i = 0; i < node->size(); i++) + { + IdxT index = node->getNodeConnection(i); + + if (index != neg1_) + { + internal_map_[index] = component_internal_idx++; + } + } + } + + for (component_type* comp : components_) + { + auto extern_indices = comp->getExternIndices(); + for (IdxT j = 0; j < comp->size(); j++) + { + if (!extern_indices.contains(j)) + { + continue; + } + IdxT index = comp->getNodeConnection(j); + + if (internal_map_.count(index) < 1 && external_map_.count(index) < 1 && index != neg1_) + { + external_map_[index] = component_internal_idx++; + } + } + } + } + + std::unordered_map internal_map_; + std::unordered_map external_map_; + std::vector external_data_indices_; + + std::vector y_ext_; + std::vector yp_ext_; + std::vector f_ext_; + + std::optional forcing_function_; + + bool is_comp_in_local_state{false}; + + }; // class SubsystemModel + +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 1906b760c..23859d440 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -23,6 +23,7 @@ namespace GridKit using component_type = CircuitComponent; using node_type = PowerElectronics::NodeBase; + protected: using CircuitComponent::size_; using CircuitComponent::n_intern_; using CircuitComponent::n_extern_; @@ -41,17 +42,6 @@ namespace GridKit using CircuitComponent::allocateVectors; public: - /** - * @brief Default constructor for the system model - * - * @post System model parameters set as default - */ - PowerElectronicsModel() - { - // By default don't use the jacobian - use_jac_ = false; - } - /** * @brief Constructor for the system model * @@ -119,7 +109,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int allocate() final + int allocate() override { size_t component_internal_size = 0; for (component_type* comp : components_) @@ -285,7 +275,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int initialize() final + int initialize() { // Initialize components for (const auto& component : components_) @@ -306,7 +296,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int distributeVectors() + virtual int distributeVectors() { const auto* y_system = y_.getData(); const auto* yp_system = yp_.getData(); @@ -336,7 +326,7 @@ namespace GridKit return 0; } - int tagDifferentiable() final + int tagDifferentiable() { return 0; } @@ -364,7 +354,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int evaluateInternalResidual() final + int evaluateInternalResidual() override { auto* f = f_.getData(); @@ -395,7 +385,7 @@ namespace GridKit for (size_t j : externals) { //@todo should do a different grounding check - if (component->getNodeConnection(j) != neg1_) + if (component->getNodeConnection(j) != neg1_ && component->getNodeConnection(j) < this->size()) { f[component->getNodeConnection(j)] += residual[j]; } @@ -423,7 +413,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int evaluateJacobian() final + int evaluateJacobian() override { distributeVectors(); @@ -447,11 +437,18 @@ namespace GridKit for (IdxT i = 0; i < nnz; ++i) { - if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) + const IdxT row = component->getNodeConnection(r[i]); + const IdxT col = component->getNodeConnection(c[i]); + + const bool is_internal_entry = row != neg1_ && col != neg1_ && row < n_intern_ && col < n_intern_; + + if (!is_internal_entry) { - vals[map_to_csr_[counter]] += v[i]; - ++counter; + continue; } + + vals[map_to_csr_[counter]] += v[i]; + ++counter; } } @@ -531,7 +528,7 @@ namespace GridKit allocated_ = false; } - private: + protected: static constexpr IdxT neg1_ = INVALID_INDEX; std::vector components_; diff --git a/GridKit/Solver/Dynamic/MultiStageCosimulation.cpp b/GridKit/Solver/Dynamic/MultiStageCosimulation.cpp new file mode 100644 index 000000000..7b1803bca --- /dev/null +++ b/GridKit/Solver/Dynamic/MultiStageCosimulation.cpp @@ -0,0 +1,192 @@ + +#include "MultiStageCosimulation.hpp"; + +#include + +#include +#include +#include +#include + +#include "GridKit/Solver/Dynamic/Interpolation.hpp" +#include "GridKit/Solver/Dynamic/Rosenbrock.hpp" + +namespace Integrator +{ + + template + MultiStageCosimulation::MultiStageCosimulation( + std::vector partitions, + std::size_t num_stages, + MemorySpace memspace) + : num_partitions_(partitions.size()), + num_stages_(num_stages), + memspace_(memspace), + partitions_(std::move(partitions)) + { + if (num_partitions_ == 0) + { + throw std::invalid_argument("MultiStageCosimulation requires at least one partition."); + } + + if (num_stages_ == 0) + { + throw std::invalid_argument("MultiStageCosimulation requires at least one stage."); + } + + for (auto* partition : partitions_) + { + if (partition == nullptr) + { + throw std::invalid_argument("Partition pointer cannot be null."); + } + } + } + + // TODO: Complete implementation + template + int MultiStageCosimulation::allocate() + { + num_partitions_ = partitions_.size(); + + x0_local_.resize(num_partitions_); + coupling_mat_.resize(num_partitions_); + linear_workspaces_.resize(num_partitions_); + lin_solver_.resize(num_partitions_); + solvers_.resize(num_partitions_); + + for (size_t p = 0; p < num_partitions_; ++p) + { + auto* part = partitions_[p]; + assert(part != nullptr); + + const IdxT num_states = part->getStateSize(); + const IdxT num_coupling = part->getCouplingSize(); + + x0_local_[p] = std::make_unique(num_states); + x0_local_[p]->allocate(memspace_); + + coupling_mat_[p] = std::make_unique(num_coupling); + coupling_mat_[p]->allocate(memspace_); + + linear_workspaces_[p] = std::make_unique(); + + lin_solver_[p] = std::make_unique( + linear_workspaces_[p].get(), + "klu", + "klu", + "klu"); + + lin_solver_[p]->initialize(); + + solvers_[p] = std::make_unique( + Rosenbrock::Tableau::rodas5p(), + part, + lin_solver_[p].get(), + &vector_handler_); + + solvers_[p]->allocate(); + } + + return 0; + } + + template + int MultiStageCosimulation::distributeLocal(const State& global_y, PartitionEvaluator* partition) + { + auto& internal = partition->getInternalIndecies(); + size_t n = partition->size(); + + matrix_handler_.matvec(); + + return 0; + } + + template + int MultiStageCosimulation::distributeCoupling(const State& global_y, int stage) + { + + auto* global_data = global_y.getData(memspace_); + + for (size_t p = 0; p < num_partitions_; ++p) + { + auto* part = partitions_[p]; + auto& external = part->getExternalIndices(); // TODO: add this function to partition evaluate + size_t n = part->size(); + + auto* local_data = coupling_mat_[p].getData(stage, memspace_); + + for (size_t i = 0; i < n; ++i) + { + local_data[i] = global_data[external[i]]; + } + } + + return 0; + } + + // TODO: Complete implementation + template + int MultiStageCosimulation::timestep(const State& global_y) + { + return 0; + } + + // TODO: Complete implementation + template + int MultiStageCosimulation::partitionStageSolve(int stage, double t0, double dt) + { + + std::vector out_times; + + for (size_t iter = 0; iter < max_iterations_; iter++) + { + for (size_t part = 0; part < num_partitions_; part++) + { + + auto* partitionEvaluator = partitions_[part]; + auto* partition_solver = solvers_[part]; + + if (stage == 0 && iter == 0) + { + distributeLocal(y_, partitionEvaluator); + } + + auto& coupling_stage_values = coupling_mat_[part]; + + GridKit::LagrangeInterpolant interp(coupling_stage_values, t0, dt); + partitionEvaluator->addForcing(interp); + + for (size_t i = 0; i < num_partitions_; i++) + { + partition_solution_[i]->setToZero(memspace_); + } + + int counter = 1; + + auto sol_out = [&](double tt, ReSolve::vector::Vector yp) -> void + { + auto* internal_csr = partition_internal_csr_[part].get(); + auto* internal_sol = partition_solution_[counter].get(); + + double alpha = 1.0; + double beta = 1.0; + + matrix_handler_.matvec(internal_csr, &yp, internal_sol, &alpha, &beta, memspace_); + + counter++; + }; + + partition_solver->integrate(out_times, stepController_, params_, sol_out); + } + + if (stage == 0) + { + break; + } + } + + return 0; + } + +} // namespace Integrator diff --git a/GridKit/Solver/Dynamic/MultiStageCosimulation.hpp b/GridKit/Solver/Dynamic/MultiStageCosimulation.hpp new file mode 100644 index 000000000..ab6aabf08 --- /dev/null +++ b/GridKit/Solver/Dynamic/MultiStageCosimulation.hpp @@ -0,0 +1,179 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Integrator +{ + + template + class MultiStageCosimulation + { + using State = ReSolve::vector::Vector; + using RealT = typename GridKit::ScalarTraits::RealT; + using Rosenbrock = Integrator::Rosenbrock; + using PartitionEvaluator = GridKit::Model::PartitionEvaluator; + using MemorySpace = ReSolve::memory::MemorySpace; + + public: + MultiStageCosimulation(std::vector partitions, + size_t num_stages = 2, + size_t max_iterations = 2, + MemorySpace memspace = MemorySpace::HOST); + + int initializeSimulation(RealT t0); + + int runSimulation(RealT tn, IdxT nout, std::optional> step_callback = {}); + + int partitionStageSolve(int, double, double); + + int timestepper(const std::vector&, + std::optional>); + + int distributeLocal(const State&); + + int distributeCoupling(const State&, int); + + int allocate(); + + private: + /** + * @brief + * + */ + std::vector> linear_workspaces_; + + /** + * @brief initial time + */ + RealT t0_; + + /** + * @brief Maximum number of iterations per stage + */ + size_t max_iterations_; + + /** + * + * @brief Number of partitions in the problem + */ + size_t num_partitions_; + + /** + * + * @brief Number of stages used by the method + */ + size_t num_stages_; + + /** + * + * @brief Number of stages used by the method + */ + size_t num_states_; + + /** + * + * @brief The global state of the system + */ + State y_; + + Rosenbrock::Parameters params_; + + Integrator::AdaptiveStep stepController_; + + /** + * + * @brief Memory space used + */ + ReSolve::memory::MemorySpace memspace_; + + /** + * + * @brief Handles vector operations + * + */ + ReSolve::VectorHandler vector_handler_; + + /** + * + * @brief Handles matrix linear algebra operations + * + */ + ReSolve::MatrixHandler matrix_handler_; + + /** + * + * @brief linear solver for the sub-integrator + * + */ + std::vector> lin_solver_; + + /** + * + * @brief Contains vectors of initial conditions for every step + * + */ + std::vector> x0_local_; + + /** + * + * @brief Vectors of coupling matrices for each partition + * + */ + std::vector>> coupling_mat_; + + /** + * + * @brief partition solution multi-dimensional vectors + * + */ + std::vector> partition_solution_; + + /** + * + * @brief partition solution multi-dimensional vectors + * + */ + std::vector partition_external_csr_; + + /** + * + * @brief partition solution multi-dimensional vectors + * + */ + std::vector> partition_internal_csr_; + + /** + * + *@brief Vector of partitions + */ + std::vector partitions_; + + /** + * + *@brief Vector of solvers + */ + std::vector> solvers_; + }; +} // namespace Integrator diff --git a/examples/PowerElectronics/CMakeLists.txt b/examples/PowerElectronics/CMakeLists.txt index 2bf311e22..bdf85eb07 100644 --- a/examples/PowerElectronics/CMakeLists.txt +++ b/examples/PowerElectronics/CMakeLists.txt @@ -10,5 +10,6 @@ if(TARGET SUNDIALS::idas) add_subdirectory(RLCircuit) add_subdirectory(Microgrid) add_subdirectory(ScaleMicrogrid) + add_subdirectory(Partition) endif() endif() diff --git a/examples/PowerElectronics/Partition/CMakeLists.txt b/examples/PowerElectronics/Partition/CMakeLists.txt new file mode 100644 index 000000000..e428f5908 --- /dev/null +++ b/examples/PowerElectronics/Partition/CMakeLists.txt @@ -0,0 +1,40 @@ +find_package(OpenMP) + +add_executable(partition Partition.cpp) + +target_link_libraries( + partition + PRIVATE GridKit::solvers_dyn GridKit::power_elec_partition_interfaces) + +add_executable(partitionMicrogrid PartitionMicrogrid.cpp) + +target_link_libraries( + partitionMicrogrid + PRIVATE GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::solvers_dyn + GridKit::power_elec_microbusdq + GridKit::power_elec_partition_interfaces) + +add_test(NAME PartitionMicrogrid COMMAND partitionMicrogrid) + +install(TARGETS partition partitionMicrogrid RUNTIME DESTINATION bin) + +if(OpenMP_CXX_FOUND) + add_executable(PartitionScaleMicrogrid PartitionScaleMicrogrid.cpp) + + target_link_libraries( + PartitionScaleMicrogrid + PRIVATE GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::solvers_dyn + GridKit::power_elec_microbusdq + GridKit::power_elec_partition_interfaces + OpenMP::OpenMP_CXX) + + add_test(NAME PartitionScaleMicrogrid COMMAND PartitionScaleMicrogrid) + + install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin) +endif() diff --git a/examples/PowerElectronics/Partition/HiresBus.hpp b/examples/PowerElectronics/Partition/HiresBus.hpp new file mode 100644 index 000000000..d37e396d3 --- /dev/null +++ b/examples/PowerElectronics/Partition/HiresBus.hpp @@ -0,0 +1,152 @@ +#pragma once + +#include +#include + +namespace GridKit +{ + /*! + * @brief Declaration of a Hires Bus Component. + * + */ + template + class HiresBus : public CircuitComponent + { + + using RealT = typename CircuitComponent::RealT; + using NodeT = typename PowerElectronics::NodeBase; + + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_; + using CircuitComponent::yp_int_; + using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_; + using CircuitComponent::f_int_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + HiresBus(NodeT* bus, IdxT id) + : node_ref_(bus) + { + size_ = 2; + n_intern_ = 0; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 2; + } + + ~HiresBus() + { + } + + int allocate() + { + CircuitComponent::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1)); + + return 0; + } + + int initialize() + { + return 0; + } + + int tagDifferentiable() + { + return 0; + } + + int evaluateInternalResidual() + { + return 0; + } + + int evaluateExternalResidual() + { + auto* y = y_.getData(); + auto* yp = yp_.getData(); + auto* f = f_.getData(); + + f[0] = -yp[0] - y[0]; + f[1] = -yp[1] - y[1]; + + f_.setDataUpdated(); + + return 0; + } + + int evaluateJacobian() + { + this->zeroJacMatrix(); + + std::vector row = {0, 1}; + std::vector col = {0, 1}; + std::vector val = {-alpha_ - 1.0, -alpha_ - 1.0}; + + this->setJacValues(row, col, val); + + return 0; + } + + int evaluateIntegrand() + { + return 0; + } + + int initializeAdjoint() + { + return 0; + } + + int evaluateAdjointResidual() + { + return 0; + } + + int evaluateAdjointIntegrand() + { + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + private: + NodeT* node_ref_; + }; +} // namespace GridKit diff --git a/examples/PowerElectronics/Partition/HiresComponent1.hpp b/examples/PowerElectronics/Partition/HiresComponent1.hpp new file mode 100644 index 000000000..503bf963d --- /dev/null +++ b/examples/PowerElectronics/Partition/HiresComponent1.hpp @@ -0,0 +1,169 @@ + + +#pragma once + +#include +#include + +namespace GridKit +{ + /*! + * @brief Declaration of a Hires Component 1 class. + * + */ + template + class HiresComponent1 : public CircuitComponent + { + using RealT = typename CircuitComponent::RealT; + using NodeT = typename PowerElectronics::NodeBase; + + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_; + using CircuitComponent::yp_int_; + using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_; + using CircuitComponent::f_int_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + HiresComponent1(NodeT* bus, IdxT id) + : node_ref_(bus) + { + size_ = 5; + n_intern_ = 3; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 12; + } + + ~HiresComponent1() + { + } + + int allocate() + { + CircuitComponent::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1)); + + return 0; + } + + int initialize() + { + return 0; + } + + int tagDifferentiable() + { + return 0; + } + + int evaluateInternalResidual() + { + auto* y = y_.getData(); + + // Internals + f_int_[0] = -yp_int_[0] - 1.71 * y_int_[0] + 0.43 * y_int_[1] + 8.32 * y_int_[2] + 0.0007; + f_int_[1] = -yp_int_[1] + 1.71 * y_int_[0] - 8.75 * y_int_[1]; + f_int_[2] = -yp_int_[2] - 10.03 * y_int_[2] + 0.43 * y[0] + 0.035 * y[1]; + + return 0; + } + + int evaluateExternalResidual() + { + auto* y = y_.getData(); + auto* f = f_.getData(); + + // outputs + f[0] = 8.32 * y_int_[1] + 1.71 * y_int_[2] - 0.1 * y[0]; + f[1] = -0.7 * y[1]; + + f_.setDataUpdated(); + + return 0; + } + + int evaluateJacobian() + { + + this->zeroJacMatrix(); + + // Internal Jacobian Entries + std::vector row = {2, 2, 2, 3, 3, 4, 4, 4}; + std::vector col = {2, 3, 4, 2, 3, 4, 0, 1}; + std::vector val = {-1.71 - alpha_, 0.43, 8.32, 1.71, -8.75 - alpha_, -10.03 - alpha_, 0.43, 0.035}; + + this->setJacValues(row, col, val); + + // External Jacobian Entries + row = {0, 0, 0, 1}; + col = {3, 4, 0, 1}; + val = {8.32, 1.71, -0.1, -0.7}; + + this->setJacValues(row, col, val); + + return 0; + } + + int evaluateIntegrand() + { + return 0; + } + + int initializeAdjoint() + { + return 0; + } + + int evaluateAdjointResidual() + { + return 0; + } + + int evaluateAdjointIntegrand() + { + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + private: + NodeT* node_ref_; + }; +} // namespace GridKit diff --git a/examples/PowerElectronics/Partition/HiresComponent3.hpp b/examples/PowerElectronics/Partition/HiresComponent3.hpp new file mode 100644 index 000000000..5e8ddf4cd --- /dev/null +++ b/examples/PowerElectronics/Partition/HiresComponent3.hpp @@ -0,0 +1,180 @@ +#pragma once + +#include +#include + +namespace GridKit +{ + /*! + * @brief Declaration of a Hires Component 3 class. + * + */ + template + class HiresComponent3 : public CircuitComponent + { + using RealT = typename CircuitComponent::RealT; + using NodeT = typename PowerElectronics::NodeBase; + + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_; + using CircuitComponent::yp_int_; + using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_; + using CircuitComponent::f_int_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + HiresComponent3(NodeT* bus, IdxT id) + : node_ref_(bus) + { + size_ = 5; + n_intern_ = 3; + n_extern_ = 2; + extern_indices_ = {0, 1}; + idc_ = id; + nnz_ = 15; + } + + ~HiresComponent3() + { + } + + int allocate() + { + CircuitComponent::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1)); + + return 0; + } + + int initialize() + { + return 0; + } + + int tagDifferentiable() + { + return 0; + } + + int evaluateInternalResidual() + { + auto* y = y_.getData(); + + // Internals + f_int_[0] = -yp_int_[0] - 280 * y_int_[0] * y_int_[2] + 0.69 * y[0] + 1.71 * y[1] - 0.43 * y_int_[0] + 0.69 * y_int_[1]; + f_int_[1] = -yp_int_[1] + 280 * y_int_[0] * y_int_[2] - 1.81 * y_int_[1]; + f_int_[2] = -yp_int_[2] - 280 * y_int_[0] * y_int_[2] + 1.81 * y_int_[1]; + + return 0; + } + + int evaluateExternalResidual() + { + auto* y = y_.getData(); + auto* f = f_.getData(); + + // Externals + f[0] = -0.02 * y[0]; + f[1] = -0.045 * y[1] + 0.43 * y_int_[0] + 0.43 * y_int_[1]; + + f_.setDataUpdated(); + + return 0; + } + + int evaluateJacobian() + { + this->zeroJacMatrix(); + + // Internal Jacobian Entries [row 1] + std::vector row = {2, 2, 2, 2, 2}; + std::vector col = {2, 3, 4, 0, 1}; + std::vector val = {-280 * y_int_[2] - 0.43 - alpha_, 0.69, -280 * y_int_[0], 0.69, 1.71}; + + this->setJacValues(row, col, val); + + // Internal Jacobian Entries [row 2] + row = {3, 3, 3}; + col = {2, 3, 4}; + val = {280 * y_int_[2], -1.81 - alpha_, 280 * y_int_[0]}; + + this->setJacValues(row, col, val); + + // Internal Jacobian Entries [row 3] + row = {4, 4, 4}; + col = {2, 3, 4}; + val = {-280 * y_int_[2], 1.81, -280 * y_int_[0] - alpha_}; + + this->setJacValues(row, col, val); + + // External Jacobian Entries + row = {0, 1, 1, 1}; + col = {0, 2, 3, 1}; + val = {-0.02, 0.43, 0.43, -0.045}; + + this->setJacValues(row, col, val); + + return 0; + } + + int evaluateIntegrand() + { + return 0; + } + + int initializeAdjoint() + { + return 0; + } + + int evaluateAdjointResidual() + { + return 0; + } + + int evaluateAdjointIntegrand() + { + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * @param rel_tol The relative tolerance which can be used to pick the + * absolute tolerance. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + * @return int 0 if successful, non-zero otherwise. + * + * This represents a "noise" level close to zero for which pure relative + * error cannot be used. + */ + int setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + private: + NodeT* node_ref_; + }; +} // namespace GridKit diff --git a/examples/PowerElectronics/Partition/Partition.cpp b/examples/PowerElectronics/Partition/Partition.cpp new file mode 100644 index 000000000..05feee692 --- /dev/null +++ b/examples/PowerElectronics/Partition/Partition.cpp @@ -0,0 +1,208 @@ +// #include + +#define _USE_MATH_DEFINES +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "HiresBus.hpp" +#include "HiresComponent1.hpp" +#include "HiresComponent3.hpp" + +std::vector hires(const std::vector& y, const std::vector& yp); + +int main(int /* argc */, char const** /* argv */) +{ + + using Bus = GridKit::PowerElectronics::MicrogridBus; + Bus bus; + + GridKit::SubsystemModel* partition1 = new GridKit::SubsystemModel(); + GridKit::SubsystemModel* partition2 = new GridKit::SubsystemModel(); + + GridKit::HiresComponent1* comp1 = new GridKit::HiresComponent1(&bus, 1); + GridKit::HiresBus* bus1 = new GridKit::HiresBus(&bus, 2); + GridKit::HiresComponent3* comp3 = new GridKit::HiresComponent3(&bus, 3); + + bus.allocate(); + + bus.setExternalConnectionNodes(0, 3); + bus.setExternalConnectionNodes(1, 4); + + comp1->allocate(); + comp3->allocate(); + bus1->allocate(); + + comp1->setExternalConnectionNodes(2, 0); + comp1->setExternalConnectionNodes(3, 1); + comp1->setExternalConnectionNodes(4, 2); + + bus1->setExternalConnectionNodes(0, 3); + bus1->setExternalConnectionNodes(1, 4); + + comp3->setExternalConnectionNodes(2, 5); + comp3->setExternalConnectionNodes(3, 6); + comp3->setExternalConnectionNodes(4, 7); + + GridKit::HiresComponent3 comp3copy(*comp3); + GridKit::BusPartitionInterface* busInterface = new GridKit::BusPartitionInterface(bus, comp3copy, 4); + + busInterface->allocate(); + partition1->addComponent(comp1); + partition1->addComponent(bus1); + partition1->addComponent(busInterface); + partition1->addNode(&bus); + + partition2->addComponent(comp3); + + auto subsystems = {partition1, partition2}; + + std::vector y = {1, 2, 3, 4, 5, 6, 7, 8}; + std::vector yp = {1, 2, 3, 4, 5, 6, 7, 8}; + + // Distribute externals to partition 1 + + for (auto* partition : subsystems) + { + partition->allocate(); + partition->updateTime(2, 5); + + for (size_t i = 0; i < partition->getExternSize(); i++) + { + partition->getExternalDataY()[i] = y[partition->getExternalDataIndices()[i]]; + partition->getExternalDataYP()[i] = yp[partition->getExternalDataIndices()[i]]; + } + + auto* partition_y = partition->y().getData(); + auto* partition_yp = partition->yp().getData(); + + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + partition_y[i] = y[partition->getNodeConnection(i)]; + partition_yp[i] = yp[partition->getNodeConnection(i)]; + } + + partition->y().setDataUpdated(); + partition->yp().setDataUpdated(); + + partition->evaluateResidual(); + partition->evaluateJacobian(); + } + + auto printHeader = [](const std::string& title) + { + std::cout << "\n------------- " << title << " -----------\n"; + + std::cout << std::left << std::setw(10) << "Res Values" + << " ---------- " + << std::right << std::setw(10) << "Comp. Index" + << '\n'; + }; + + auto printRow = [](double residual, size_t index) + { + std::cout << std::left << std::setw(10) + << std::defaultfloat << std::setprecision(5) + << residual + << " ---------- " + << std::right << std::setw(10) + << index + << '\n'; + }; + + auto printPartitionResiduals = [&](auto* partition, const std::string& title) + { + printHeader(title); + + auto* residual = partition->getResidual().getData(); + + for (size_t i = 0; i < partition->getInternalSize(); ++i) + { + const auto comp_index = partition->getNodeConnection(i); + printRow(residual[i], comp_index); + } + }; + + printPartitionResiduals(partition1, "Partition 1"); + printPartitionResiduals(partition2, "Partition 2"); + + auto ref = hires(y, yp); + + printHeader("Reference Output"); + for (size_t i = 0; i < ref.size(); ++i) + { + printRow(ref[i], i); + } + + delete comp1; + delete bus1; + delete busInterface; + delete comp3; + + delete partition1; + delete partition2; + + return 0; +} + +std::vector hires(const std::vector& y, const std::vector& yp) +{ + std::vector f(8); + + f[0] = -yp[0] - 1.71 * y[0] + 0.43 * y[1] + 8.32 * y[2] + 0.0007; + f[1] = -yp[1] + 1.71 * y[0] - 8.75 * y[1]; + f[2] = -yp[2] - 10.03 * y[2] + 0.43 * y[3] + 0.035 * y[4]; + f[3] = -yp[3] + 8.32 * y[1] + 1.71 * y[2] - 1.12 * y[3]; + f[4] = -yp[4] - 1.745 * y[4] + 0.43 * y[5] + 0.43 * y[6]; + f[5] = -yp[5] - 280 * y[5] * y[7] + 0.69 * y[3] + 1.71 * y[4] - 0.43 * y[5] + 0.69 * y[6]; + f[6] = -yp[6] + 280 * y[5] * y[7] - 1.81 * y[6]; + f[7] = -yp[7] - 280 * y[5] * y[7] + 1.81 * y[6]; + + return f; +} + +// template +// int jac_hires(const std::vector& y, [[maybe_unused]] const std::vector& yp, GridKit::LinearAlgebra::DenseMatrix& jac, RealT alpha) +// { + +// jac.setValue(0, 0, -alpha - 1.71); +// jac.setValue(0, 1, 0.43); +// jac.setValue(0, 2, 8.32); + +// jac.setValue(1, 0, 1.71); +// jac.setValue(1, 1, -alpha - 8.75); + +// jac.setValue(2, 2, -alpha - 10.03); +// jac.setValue(2, 3, 0.43); +// jac.setValue(2, 4, 0.035); + +// jac.setValue(3, 1, 8.32); +// jac.setValue(3, 2, 1.71); +// jac.setValue(3, 3, -alpha - 1.12); + +// jac.setValue(4, 4, -alpha - 1.745); +// jac.setValue(4, 5, 0.43); +// jac.setValue(4, 6, 0.43); + +// jac.setValue(5, 3, 0.69); +// jac.setValue(5, 4, 1.71); +// jac.setValue(5, 5, -alpha - 280.0 * y[7] - 0.43); +// jac.setValue(5, 6, 0.69); +// jac.setValue(5, 7, -280.0 * y[5]); + +// jac.setValue(6, 5, 280.0 * y[7]); +// jac.setValue(6, 6, -alpha - 1.81); +// jac.setValue(6, 7, 280.0 * y[5]); + +// jac.setValue(7, 5, -280.0 * y[7]); +// jac.setValue(7, 6, 1.81); +// jac.setValue(7, 7, -alpha - 280.0 * y[5]); + +// return 0; +// } diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp new file mode 100644 index 000000000..c9a984950 --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -0,0 +1,332 @@ + +#include +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "jac_test_helper.hpp" + +int main() +{ + /// @todo Needs to be modified. Some components are small relative to others thus + /// there error is high (or could be matlab vector issue) + + bool use_jac = true; + + // Create model + auto* sysmodel = new GridKit::PowerElectronicsModel(use_jac); + + // Modeled after the problem in the paper + double RN = 1.0e4; + + // DG Params + GridKit::DistributedGeneratorParameters parms1; + parms1.wb_ = 2.0 * M_PI * 50.0; + parms1.wc_ = 31.41; + parms1.mp_ = 9.4e-5; + parms1.Vn_ = 380.0; + parms1.nq_ = 1.3e-3; + parms1.F_ = 0.75; + parms1.Kiv_ = 420.0; + parms1.Kpv_ = 0.1; + parms1.Kic_ = 2.0e4; + parms1.Kpc_ = 15.0; + parms1.Cf_ = 5.0e-5; + parms1.rLf_ = 0.1; + parms1.Lf_ = 1.35e-3; + parms1.rLc_ = 0.03; + parms1.Lc_ = 0.35e-3; + + GridKit::DistributedGeneratorParameters parms2; + // Parameters from MATLAB Microgrid code for first DG + parms2.wb_ = 2.0 * M_PI * 50.0; + parms2.wc_ = 31.41; + parms2.mp_ = 12.5e-5; + parms2.Vn_ = 380.0; + parms2.nq_ = 1.5e-3; + parms2.F_ = 0.75; + parms2.Kiv_ = 390.0; + parms2.Kpv_ = 0.05; + parms2.Kic_ = 16.0e3; + parms2.Kpc_ = 10.5; + parms2.Cf_ = 50.0e-6; + parms2.rLf_ = 0.1; + parms2.Lf_ = 1.35e-3; + parms2.rLc_ = 0.03; + parms2.Lc_ = 0.35e-3; + + // Line params + double rline1 = 0.23; + double Lline1 = 0.1 / (2.0 * M_PI * 50.0); + + double rline2 = 0.35; + double Lline2 = 0.58 / (2.0 * M_PI * 50.0); + + double rline3 = 0.23; + double Lline3 = 0.1 / (2.0 * M_PI * 50.0); + + // load parms + double rload1 = 3.0; + double Lload1 = 2.0 / (2.0 * M_PI * 50.0); + + double rload2 = 2.0; + double Lload2 = 1.0 / (2.0 * M_PI * 50.0); + + using SignalNode = GridKit::PowerElectronics::SignalNode; + using Subsystem = GridKit::SubsystemModel; + + SignalNode dg_signal; + + sysmodel->addNode(&dg_signal); + + using Bus = GridKit::PowerElectronics::MicrogridBus; + Bus bus1; + Bus bus2; + Bus bus3; + Bus bus4; + + sysmodel->addNode(&bus1); + sysmodel->addNode(&bus2); + sysmodel->addNode(&bus3); + sysmodel->addNode(&bus4); + + // dg 1 + GridKit::DistributedGenerator* dg1 = new GridKit::DistributedGenerator( + 0, parms1, true, &dg_signal, &bus1); + sysmodel->addComponent(dg1); + + // dg 2 + GridKit::DistributedGenerator* dg2 = new GridKit::DistributedGenerator( + 1, parms1, false, &dg_signal, &bus2); + sysmodel->addComponent(dg2); + + // dg 3 + GridKit::DistributedGenerator* dg3 = new GridKit::DistributedGenerator( + 2, parms2, false, &dg_signal, &bus3); + sysmodel->addComponent(dg3); + + // dg 4 + GridKit::DistributedGenerator* dg4 = new GridKit::DistributedGenerator( + 3, parms2, false, &dg_signal, &bus4); + sysmodel->addComponent(dg4); + + // Lines + + // line 1 + GridKit::MicrogridLine* l1 = new GridKit::MicrogridLine( + 4, rline1, Lline1, &dg_signal, &bus1, &bus2); + sysmodel->addComponent(l1); + + // line 2 + GridKit::MicrogridLine* l2 = new GridKit::MicrogridLine( + 5, rline2, Lline2, &dg_signal, &bus2, &bus3); + sysmodel->addComponent(l2); + + // line 3 + GridKit::MicrogridLine* l3 = new GridKit::MicrogridLine( + 6, rline3, Lline3, &dg_signal, &bus3, &bus4); + sysmodel->addComponent(l3); + + // loads + + // load 1 + GridKit::MicrogridLoad* load1 = new GridKit::MicrogridLoad(7, rload1, Lload1, &dg_signal, &bus1); + sysmodel->addComponent(load1); + + // load 2 + GridKit::MicrogridLoad* load2 = new GridKit::MicrogridLoad(8, rload2, Lload2, &dg_signal, &bus3); + sysmodel->addComponent(load2); + + // Virtual PQ Buses + GridKit::MicrogridBusDQ* bus_para_1 = new GridKit::MicrogridBusDQ(9, RN, &bus1); + sysmodel->addComponent(bus_para_1); + + GridKit::MicrogridBusDQ* bus_para_2 = new GridKit::MicrogridBusDQ(10, RN, &bus2); + sysmodel->addComponent(bus_para_2); + + GridKit::MicrogridBusDQ* bus_para_3 = new GridKit::MicrogridBusDQ(11, RN, &bus3); + sysmodel->addComponent(bus_para_3); + + GridKit::MicrogridBusDQ* bus_para_4 = new GridKit::MicrogridBusDQ(12, RN, &bus4); + sysmodel->addComponent(bus_para_4); + + sysmodel->allocate(); + + Subsystem* partition1 = new Subsystem(); + Subsystem* partition2 = new Subsystem(); + + GridKit::MicrogridLine l2copy(*l2); + GridKit::BusPartitionInterface* busInterface1 = new GridKit::BusPartitionInterface(bus2, l2copy, 14); + + busInterface1->allocate(); + + partition1->addNode(&dg_signal); + partition1->addComponent(dg1); + partition1->addComponent(dg2); + partition1->addComponent(l1); + partition1->addComponent(load1); + partition1->addComponent(busInterface1); + partition1->addComponent(bus_para_1); + partition1->addComponent(bus_para_2); + partition1->addNode(&bus1); + partition1->addNode(&bus2); + + partition2->addComponent(dg3); + partition2->addComponent(dg4); + partition2->addComponent(l2); + partition2->addComponent(l3); + partition2->addComponent(load2); + partition2->addComponent(bus_para_4); + partition2->addComponent(bus_para_3); + partition2->addNode(&bus3); + partition2->addNode(&bus4); + + std::vector y; + std::vector yp; + + for (size_t i = 0; i < sysmodel->size(); i++) + { + y.push_back(static_cast(i + 1)); + yp.push_back(static_cast(i + 1)); + } + + auto* system_y = sysmodel->y().getData(); + auto* system_yp = sysmodel->yp().getData(); + + for (size_t i = 0; i < sysmodel->size(); i++) + { + system_y[i] = y[i]; + system_yp[i] = yp[i]; + } + + sysmodel->y().setDataUpdated(); + sysmodel->yp().setDataUpdated(); + + sysmodel->updateTime(2, 5); + sysmodel->evaluateResidual(); + sysmodel->evaluateJacobian(); + + auto full_jac = sysmodel->getCsrJacobian(); + + auto* f_sysmodel = sysmodel->getResidual().getData(); + + std::vector*> partitions = {partition1, partition2}; + + // Distribute externals to partition 1 + for (auto* partition : partitions) + { + partition->allocate(); + + // test hold and release methods + partition->release(); + partition->hold(); + + partition->updateTime(2, 5); + + // Test setTimeFunction function + auto forcing = [&]([[maybe_unused]] double t) -> Subsystem::ForcingData + { + auto ext_size = partition->getExternSize(); + + std::vector y_ext(ext_size); + std::vector yp_ext(ext_size); + + for (size_t i = 0; i < partition->getExternSize(); i++) + { + y_ext[i] = y[partition->getExternalDataIndices()[i]]; + yp_ext[i] = yp[partition->getExternalDataIndices()[i]]; + } + + return {y_ext, yp_ext}; + }; + + partition->setTimeFunction(forcing); + + auto* partition_y = partition->y().getData(); + auto* partition_yp = partition->yp().getData(); + + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + partition_y[i] = y[partition->getNodeConnection(i)]; + partition_yp[i] = yp[partition->getNodeConnection(i)]; + } + + partition->y().setDataUpdated(); + partition->yp().setDataUpdated(); + + partition->evaluateResidual(); + partition->evaluateJacobian(); + } + + auto partition1_jac = partition1->getCsrJacobian(); + auto partition2_jac = partition2->getCsrJacobian(); + + bool jac_matched_1 = GridKit::Testing::verifySubsystemJacobian(*full_jac, *partition1_jac, *partition1); + bool jac_matched_2 = GridKit::Testing::verifySubsystemJacobian(*full_jac, *partition2_jac, *partition2); + + std::cout << "Jacobian Matched: " << jac_matched_1 * jac_matched_2 << std::endl; + + if (!jac_matched_1 || !jac_matched_2) + { + std::cout << "ERROR: At least one subsystem Jacobian is incorrect!" << std::endl; + return 1; + } + + std::vector f(sysmodel->size(), 0.0); + std::vector error(sysmodel->size(), 1.0); + + // Get internal residuals from partition 1 + for (auto* partition : partitions) + { + const auto* residual = partition->getResidual().getData(); + + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + + f[partition->getNodeConnection(i)] = residual[i]; + } + } + + double max_error = 0; + for (size_t i = 0; i < sysmodel->size(); i++) + { + error[i] = (f_sysmodel[i] - f[i]) / f_sysmodel[i]; + + if (max_error < std::abs(error[i])) + { + max_error = std::abs(error[i]); + } + } + + std::cout << "\nMax Error of Reference and Partition Evaluation: " << max_error << std::endl; + + if (max_error > std::numeric_limits::epsilon()) + { + std::cout << "ERROR: Max Error too high!" << std::endl; + return 1; + } + + delete sysmodel; + delete partition1; + delete partition2; + delete busInterface1; + + return 0; +} diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp new file mode 100644 index 000000000..13b10a4ad --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -0,0 +1,536 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "jac_test_helper.hpp" + +/****************************************************************************** + * Partitioned Residual Evaluation + * + * Construct subsystem models by partitioning the original microgrid into + * independent groups of components. Each subsystem is allocated and receives + * the appropriate subset of the global state vectors (y and yp), consisting + * of: + * + * - External variables (coupling variables from neighboring partitions) + * - Internal variables (states owned by the partition) + * + * After distributing the state information, each partition independently + * evaluates its residual. The local residuals are then scattered back into + * the global residual vector to reconstruct the monolithic residual. + * + * Finally, the reconstructed residual is compared against the reference + * monolithic evaluation to verify the correctness of the partitioning + * implementation. + ******************************************************************************/ + +using index_type = size_t; +using real_type = double; + +struct RunResult +{ + bool success = true; + index_type num_partitions; + real_type partition_time; // seconds + real_type monolithic_time; // seconds + real_type speedup; // monolithic / partition + real_type max_error; + std::string subsys_jac; +}; + +RunResult printMicrogridSystems(index_type N_size, index_type num_partitions); + +/** + * @brief Run Scale Microgrid for a specific N given by the user. + * + * @param[in] argc should be 1 if no N given, 2 if N given + * @param[in] argv can contain the N value (argv[1]) + * @return int + */ +int main() +{ + index_type N_size = 5000; + + std::cout << std::format("{:<16}{:>16}{:>18}{:>12}{:>14}{:>16}\n", + "num_partitions", + "partition_time", + "monolithic_time", + "speedup", + "error", + "Jacobians"); + + std::cout << std::string(93, '-') << "\n"; + + for (index_type p : {10, 48, 100, 500, 1000, 2000, 5000}) + { + assert(p <= 2 * N_size); + + RunResult r = printMicrogridSystems(N_size, p); + + if (!r.success) + { + return 1; + } + + std::cout << std::format("{:<16d}{:>14.4f} s{:>16.4f} s{:>11.2f}x{:>14.3e} {:>16s}\n", + r.num_partitions, + r.partition_time, + r.monolithic_time, + r.speedup, + r.max_error, + r.subsys_jac); + } + + return 0; +} + +/** + * @brief Tests network of distributed generators. + * + * @param[in] N_size - The number of DG line load cobinations to generate for scale + * @return int returns 0 if successful, >0 otherwise + */ +RunResult printMicrogridSystems(index_type N_size, index_type num_partitions) +{ + using namespace GridKit; + + const index_type num_ibrs = 2 * N_size; + + assert(num_partitions <= num_ibrs); + + bool use_jac = true; + + // Create circuit model + auto* sys_model_control = new PowerElectronicsModel(use_jac); + + // Ensure minimum size requirement + assert(N_size >= 1); + + // Modeled after the problem in the paper + // Every Bus has the same virtual resistance. This is due to numerical stability as mentioned in the paper. + real_type RN = 1.0e4; + + // DG Params Vector + // All DGs have the same set of parameters except for the first two. + GridKit::DistributedGeneratorParameters DG_parms1; + DG_parms1.wb_ = 2.0 * M_PI * 50.0; + DG_parms1.wc_ = 31.41; + DG_parms1.mp_ = 9.4e-5; + DG_parms1.Vn_ = 380.0; + DG_parms1.nq_ = 1.3e-3; + DG_parms1.F_ = 0.75; + DG_parms1.Kiv_ = 420.0; + DG_parms1.Kpv_ = 0.1; + DG_parms1.Kic_ = 2.0e4; + DG_parms1.Kpc_ = 15.0; + DG_parms1.Cf_ = 5.0e-5; + DG_parms1.rLf_ = 0.1; + DG_parms1.Lf_ = 1.35e-3; + DG_parms1.rLc_ = 0.03; + DG_parms1.Lc_ = 0.35e-3; + + GridKit::DistributedGeneratorParameters DG_parms2; + DG_parms2.wb_ = 2.0 * M_PI * 50.0; + DG_parms2.wc_ = 31.41; + DG_parms2.mp_ = 12.5e-5; + DG_parms2.Vn_ = 380.0; + DG_parms2.nq_ = 1.5e-3; + DG_parms2.F_ = 0.75; + DG_parms2.Kiv_ = 390.0; + DG_parms2.Kpv_ = 0.05; + DG_parms2.Kic_ = 16.0e3; + DG_parms2.Kpc_ = 10.5; + DG_parms2.Cf_ = 50.0e-6; + DG_parms2.rLf_ = 0.1; + DG_parms2.Lf_ = 1.35e-3; + DG_parms2.rLc_ = 0.03; + DG_parms2.Lc_ = 0.35e-3; + + std::vector> DGParams_list(2 * N_size, DG_parms2); + + // First two generators use parameters 1 + if (DGParams_list.size() >= 1) + { + DGParams_list[0] = DG_parms1; + } + if (DGParams_list.size() >= 2) + { + DGParams_list[1] = DG_parms1; + } + + // line vector params + // Every odd line has the same parameters and every even line has the same parameters + real_type rline1 = 0.23; + real_type Lline1 = 0.1 / (2.0 * M_PI * 50.0); + real_type rline2 = 0.35; + real_type Lline2 = 0.58 / (2.0 * M_PI * 50.0); + std::vector rline_list(2 * N_size - 1, 0.0); + std::vector Lline_list(2 * N_size - 1, 0.0); + for (index_type i = 0; i < rline_list.size(); i++) + { + rline_list[i] = (i % 2) ? rline2 : rline1; + Lline_list[i] = (i % 2) ? Lline2 : Lline1; + } + + // load parms + // Only the first load has the same paramaters. + real_type rload1 = 3.0; + real_type Lload1 = 2.0 / (2.0 * M_PI * 50.0); + real_type rload2 = 2.0; + real_type Lload2 = 1.0 / (2.0 * M_PI * 50.0); + + std::vector rload_list(N_size, rload2); + std::vector Lload_list(N_size, Lload2); + if (rload_list.size() >= 1) + { + rload_list[0] = rload1; + Lload_list[0] = Lload1; + } + + using SignalNode = GridKit::PowerElectronics::SignalNode; + using Bus = GridKit::PowerElectronics::MicrogridBus; + using BusDQ = MicrogridBusDQ; + using Generator = DistributedGenerator; + using Line = MicrogridLine; + using Load = MicrogridLoad; + using PartitionInterface = BusPartitionInterface; + + std::vector buses(num_ibrs, nullptr); + std::vector busesDQ(num_ibrs, nullptr); + std::vector generators(num_ibrs, nullptr); + std::vector lines(num_ibrs, nullptr); + std::vector loads(num_ibrs, nullptr); + std::vector partitionInterface; + std::vector linesCopies; + + SignalNode dg_signal; + // sys_model_control->addNode(&dg_signal); + + for (size_t i = 0; i < 2 * N_size; i++) + { + buses[i] = new Bus(); + // sys_model_control->addNode(buses[i]); + } + + // Create the reference DG + auto* dg_ref = new DistributedGenerator(0, + DGParams_list[0], + true, + &dg_signal, + buses[0]); + // sys_model_control->addComponent(dg_ref); + + generators[0] = dg_ref; + + // Keep track of models and index location + index_type model_id = 1; + // Add all other DGs + for (index_type i = 1; i < 2 * N_size; i++) + { + // current DG to add + auto* dg = new DistributedGenerator(model_id++, + DGParams_list[i], + false, + &dg_signal, + buses[i]); + + generators[i] = dg; + // sys_model_control->addComponent(dg); + } + + // Load all the Line components + for (index_type i = 0; i < 2 * N_size - 1; i++) + { + // line + auto* line_model = new MicrogridLine(model_id++, + rline_list[i], + Lline_list[i], + &dg_signal, + buses[i], + buses[i + 1]); + lines[i + 1] = line_model; + // sys_model_control->addComponent(line_model); + } + + // Load all the Load components + for (index_type i = 0; i < N_size; i++) + { + auto* load_model = new MicrogridLoad(model_id++, + rload_list[i], + Lload_list[i], + &dg_signal, + buses[2 * i]); + loads[2 * i] = load_model; + // sys_model_control->addComponent(load_model); + } + + // Add all the microgrid Virtual DQ Buses + for (index_type i = 0; i < 2 * N_size; i++) + { + auto* virDQbus_model = new MicrogridBusDQ(model_id++, RN, buses[i]); + + busesDQ[i] = virDQbus_model; + // sys_model_control->addComponent(virDQbus_model); + } + + sys_model_control->addNode(&dg_signal); + for (index_type i = 0; i < num_ibrs; ++i) + { + sys_model_control->addComponent(generators[i]); + sys_model_control->addComponent(busesDQ[i]); + + if (loads[i] != nullptr) + { + sys_model_control->addComponent(loads[i]); + } + + if (lines[i] != nullptr) + { + sys_model_control->addComponent(lines[i]); + } + sys_model_control->addNode(buses[i]); + } + + sys_model_control->allocate(); + + std::vector y; + std::vector yp; + + for (size_t i = 0; i < sys_model_control->size(); i++) + { + y.push_back(static_cast(i + 1)); + yp.push_back(static_cast(i + 1)); + } + + auto start_time = std::chrono::high_resolution_clock::now(); + + auto* sys_model_y = sys_model_control->y().getData(); + auto* sys_model_yp = sys_model_control->yp().getData(); + + for (size_t i = 0; i < sys_model_control->size(); i++) + { + sys_model_y[i] = y[i]; + sys_model_yp[i] = yp[i]; + } + + sys_model_control->y().setDataUpdated(); + sys_model_control->yp().setDataUpdated(); + + sys_model_control->updateTime(2, 5); + sys_model_control->evaluateResidual(); + + auto end_time = std::chrono::high_resolution_clock::now(); + auto monolithic_time = std::chrono::duration(end_time - start_time); + + sys_model_control->evaluateJacobian(); + + auto* f_sysmodel_control = sys_model_control->getResidual().getData(); + auto* full_jac = sys_model_control->getCsrJacobian(); + + //--------------------------------------------------------------- + // Partition the monolithic network into independent subsystem + // models. Each subsystem owns a contiguous block of generators, + // buses, lines, and loads. Whenever a partition boundary is + // encountered, a BusPartitionInterface is inserted to expose the + // neighboring partition's boundary variables while maintaining + // the same electrical connectivity as the original network. + //--------------------------------------------------------------- + + std::vector*> subsystems(num_partitions); + + // Partition the system + index_type q = (num_ibrs) / num_partitions; + index_type r = (num_ibrs) % num_partitions; + index_type index = 0; + for (index_type j = 0; j < num_partitions; j++) + { + auto* partition = new SubsystemModel(); + + // add Reference rotor to the first partition + if (j == 0) + { + partition->addNode(&dg_signal); + } + + index_type part_size = q + (j < r ? 1 : 0); + index_type end = std::min(index + part_size, num_ibrs); + + // Add all components belonging to this partition + for (; index < end; ++index) + { + partition->addComponent(generators[index]); + partition->addComponent(busesDQ[index]); + + if (loads[index] != nullptr) + { + partition->addComponent(loads[index]); + } + + if (lines[index] != nullptr) + { + partition->addComponent(lines[index]); + } + + partition->addNode(buses[index]); + } + // Add the partition interface to the left partition + if (index < num_ibrs) + { + auto* linecopy = new GridKit::MicrogridLine(*lines[index]); + auto* busInterface = new GridKit::BusPartitionInterface(*buses[index - 1], *linecopy, model_id++); + + busInterface->allocate(); + partitionInterface.push_back(busInterface); + linesCopies.push_back(linecopy); + + partition->addComponent(busInterface); + } + + subsystems[j] = partition; + } + + //--------------------------------------------------------------- + // Allocate each subsystem and evaluate the partitioned residual. + // + // The global state vectors (y and yp) are scattered into the + // corresponding internal and external vectors of every subsystem. + // After evaluating the residuals independently, the local + // residuals are gathered back into the global residual vector for + // comparison with the monolithic evaluation. + //--------------------------------------------------------------- + + std::vector f(sys_model_control->size(), 1.0); + std::vector error(sys_model_control->size(), 1.0); + + for (auto* partition : subsystems) + { + partition->allocate(); + partition->updateTime(2, 5); + } + + start_time = std::chrono::high_resolution_clock::now(); + +// Distribute externals to partition 1 +#pragma omp parallel for schedule(guided) + for (auto* partition : subsystems) + { + // Distribute external variables + for (size_t i = 0; i < partition->getExternSize(); i++) + { + partition->getExternalDataY()[i] = y[partition->getExternalDataIndices()[i]]; + partition->getExternalDataYP()[i] = yp[partition->getExternalDataIndices()[i]]; + } + + auto* partition_y = partition->y().getData(); + auto* partition_yp = partition->yp().getData(); + + // Distribute internal variables + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + partition_y[i] = y[partition->getNodeConnection(i)]; + partition_yp[i] = yp[partition->getNodeConnection(i)]; + } + + partition->y().setDataUpdated(); + partition->yp().setDataUpdated(); + + // Evaluate residual of this partition + partition->evaluateResidual(); + + auto* residual = partition->getResidual().getData(); + // Reconstructs the monolithic residual from the partition residuals + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + f[partition->getNodeConnection(i)] = residual[i]; + } + } + + end_time = std::chrono::high_resolution_clock::now(); + auto partition_time = std::chrono::duration(end_time - start_time); + + //--------------------------------------------------------------- + // Compare the reconstructed partition residual against the + // reference monolithic residual. The per-entry error is computed + // and the maximum error is reported to verify that the partitioned + // formulation reproduces the monolithic model within numerical + // roundoff. + //--------------------------------------------------------------- + + bool matched = true; + + for (auto* partition : subsystems) + { + partition->evaluateJacobian(); + + matched &= GridKit::Testing::verifySubsystemJacobian( + *full_jac, + *partition->getCsrJacobian(), + *partition); + } + + real_type max_error = 0; + for (size_t i = 0; i < sys_model_control->size(); i++) + { + error[i] = abs(f[i] - f_sysmodel_control[i]) / (f_sysmodel_control[i] + 1); + + if (max_error < error[i]) + { + max_error = error[i]; + } + } + + RunResult result; + result.num_partitions = num_partitions; + result.partition_time = partition_time.count(); + result.monolithic_time = monolithic_time.count(); + result.speedup = monolithic_time.count() / partition_time.count(); + result.max_error = max_error; + result.subsys_jac = matched ? "Correct" : "Wrong"; + + if (!matched) + { + std::cout << "ERROR: At least one subsystem Jacobian is incorrect!" << std::endl; + result.success = false; + } + + if (max_error > std::numeric_limits::epsilon()) + { + std::cout << "ERROR: Max Error too high!" << std::endl; + result.success = false; + } + + for (auto* linescpy : linesCopies) + { + delete linescpy; + } + for (auto* partition : subsystems) + { + delete partition; + } + for (auto* part_interface : partitionInterface) + { + delete part_interface; + } + delete sys_model_control; + + return result; +} diff --git a/examples/PowerElectronics/Partition/jac_test_helper.hpp b/examples/PowerElectronics/Partition/jac_test_helper.hpp new file mode 100644 index 000000000..8851bc2a3 --- /dev/null +++ b/examples/PowerElectronics/Partition/jac_test_helper.hpp @@ -0,0 +1,250 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + + template + bool verifySubsystemJacobian( + GridKit::LinearAlgebra::CsrMatrix& full_jac, + GridKit::LinearAlgebra::CsrMatrix& sub_jac, + GridKit::SubsystemModel& subsystem, + RealT tolerance = static_cast(1e-12)) + { + constexpr auto host = memory::HOST; + // Full-system CSR data. + const IdxT* full_rows = full_jac.getRowData(host); + const IdxT* full_cols = full_jac.getColData(host); + const RealT* full_vals = full_jac.getValues(host); + + // Subsystem CSR data. + const IdxT* sub_rows = sub_jac.getRowData(host); + const IdxT* sub_cols = sub_jac.getColData(host); + const RealT* sub_vals = sub_jac.getValues(host); + + const IdxT num_sub_rows = sub_jac.getNumRows(); + const IdxT num_sub_cols = sub_jac.getNumColumns(); + + /* + * The subsystem internal map stores: + * + * global index -> local index + */ + const auto& global_to_local = subsystem.getInternalMap(); + + if (global_to_local.size() != num_sub_rows) + { + std::cout << "Internal map size mismatch: map has " + << global_to_local.size() + << " entries, but subsystem Jacobian has " + << num_sub_rows + << " rows\n"; + + return false; + } + + if (num_sub_rows != num_sub_cols) + { + std::cout << "Subsystem Jacobian is not square: " + << num_sub_rows + << " rows and " + << num_sub_cols + << " columns\n"; + + return false; + } + + /* + * Build the reverse mapping: + * + * local index -> global index + * + * This is needed because the subsystem Jacobian stores local row + * and column indices. + */ + std::vector local_to_global(num_sub_rows); + std::vector local_index_found(num_sub_rows, false); + + for (const auto& [global_index, local_index] : global_to_local) + { + + if (global_index >= full_jac.getNumRows()) + { + std::cout << "Invalid global index " + << global_index + << " in subsystem internal map\n"; + + return false; + } + + if (local_index >= num_sub_rows) + { + std::cout << "Invalid local index " + << local_index + << " mapped from global index " + << global_index << '\n'; + + return false; + } + + if (local_index_found[local_index]) + { + std::cout << "Duplicate local index " + << local_index + << " in subsystem internal map\n"; + + return false; + } + + local_to_global[local_index] = global_index; + local_index_found[local_index] = true; + } + + // Verify that every local index has a corresponding global index. + for (IdxT local_index = 0; local_index < num_sub_rows; ++local_index) + { + if (!local_index_found[local_index]) + { + std::cout << "No global index maps to local index " + << local_index << '\n'; + + return false; + } + } + + bool matches = true; + + /* + * Compare each subsystem row against the corresponding full-system row. + */ + for (IdxT local_row = 0; local_row < num_sub_rows; ++local_row) + { + const IdxT global_row = local_to_global[local_row]; + + const IdxT sub_begin = sub_rows[local_row]; + const IdxT sub_end = sub_rows[local_row + 1]; + const IdxT full_begin = full_rows[global_row]; + const IdxT full_end = full_rows[global_row + 1]; + + /* + * Store this subsystem row using global column indices. + * + * The subsystem CSR row is sorted by local column index. After + * converting local columns to global columns, the global columns + * may no longer be sorted. Therefore, we use a lookup map instead + * of comparing both CSR rows sequentially. + */ + std::unordered_map sub_row_entries; + + for (IdxT sub_index = sub_begin; sub_index < sub_end; ++sub_index) + { + const IdxT local_column = sub_cols[sub_index]; + + if (local_column >= num_sub_cols) + { + std::cout << "Invalid subsystem column index " + << local_column + << " in local row " + << local_row << '\n'; + + matches = false; + continue; + } + + const IdxT global_column = local_to_global[local_column]; + + // Check for duplicate columns in this subsystem row. + if (sub_row_entries.find(global_column) != sub_row_entries.end()) + { + std::cout << "Duplicate subsystem entry at (" + << global_row << ", " + << global_column << ")\n"; + + matches = false; + continue; + } + + sub_row_entries[global_column] = sub_vals[sub_index]; + } + + /* + * Compare every full-system entry whose column belongs to this + * subsystem. + */ + for (IdxT full_index = full_begin; + full_index < full_end; + ++full_index) + { + const IdxT global_column = full_cols[full_index]; + + // Ignore columns owned by another subsystem. + if (global_to_local.find(global_column) == global_to_local.end()) + { + continue; + } + + auto sub_entry = sub_row_entries.find(global_column); + + if (sub_entry == sub_row_entries.end()) + { + std::cout << "Entry exists only in full Jacobian at (" + << global_row << ", " + << global_column << ")\n"; + + matches = false; + continue; + } + + const RealT full_value = full_vals[full_index]; + const RealT sub_value = sub_entry->second; + const RealT difference = std::abs(full_value - sub_value); + + if (difference > tolerance) + { + std::cout << "Jacobian value mismatch at (" + << global_row << ", " + << global_column << "): " + << "full = " << full_value + << ", subsystem = " << sub_value + << ", difference = " << difference << '\n'; + + matches = false; + } + + /* + * Remove the matched subsystem entry. + * + * After the full row has been checked, any remaining entries + * exist only in the subsystem Jacobian. + */ + sub_row_entries.erase(sub_entry); + } + + // Report subsystem entries that were not found in the full Jacobian. + for (const auto& entry : sub_row_entries) + { + const IdxT global_column = entry.first; + + std::cout << "Entry exists only in subsystem Jacobian at (" + << global_row << ", " + << global_column << ")\n"; + + matches = false; + } + } + + return matches; + } + + } // namespace Testing +} // namespace GridKit