From 2990d162b1056110ea77a285bca188d359576a13 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 07:36:54 +0000 Subject: [PATCH 1/9] Add BusPartitionInterface to define partition boundaries --- .../BusPartitionInterface.cpp | 365 ++++++++++++++++++ .../BusPartitionInterface.hpp | 164 ++++++++ .../PartitionInterface/CMakeLists.txt | 9 + .../PartitionInterface/PartitionInterface.hpp | 29 ++ 4 files changed, 567 insertions(+) create mode 100644 GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp create mode 100644 GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp create mode 100644 GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt create mode 100644 GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp new file mode 100644 index 000000000..d4619f3f1 --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp @@ -0,0 +1,365 @@ + +#include "BusPartitionInterface.hpp" + +#include +#include +#include + +#include +#include + +namespace GridKit +{ + + /** + * @brief Construct a partition interface between a bus and a circuit component. + * + * The interface wraps a copy of a component located across a partition + * boundary and evaluates the contributions that the component makes to the + * connected bus. All interface variables are treated as external variables + * and have the same size as the wrapped component. + * + * + * @param bus Bus associated with the partition interface. + * @param component Copy of the component across the partition boundary. + * @param id Unique identifier for the interface. + */ + template + BusPartitionInterface::BusPartitionInterface(node_type* bus, component_type* component, IdxT id) + : component_(component->clone()), + bus_(bus) + { + size_ = component_->size(); + n_intern_ = 0; + n_extern_ = static_cast(component_->size()); + idc_ = id; + + // All variables of the bus interface are external to the interface. + for (IdxT i = 0; i < size_; i++) + { + extern_indices_.insert(i); + } + + // Map each global bus connection index to its position in the bus. + std::unordered_map bus_connections; + + for (size_t i = 0; i < static_cast(bus_->size()); ++i) + { + bus_connections[bus_->getNodeConnection(i).idx_] = i; + } + + // The interface Jacobian contains only rows corresponding to variables + // owned by the bus. + const IdxT* coo_rows = component_->jacobianCooRows(); + + nnz_ = 0; + + for (IdxT k = 0; k < component_->nnz(); ++k) + { + const IdxT row_node = component_->getNodeConnection(static_cast(coo_rows[k])); + + if (bus_connections.contains(row_node)) + { + ++nnz_; + jac_map_.push_back(k); // Keep track of the entries so they can be easily extracted + } + } + } + + template + BusPartitionInterface::~BusPartitionInterface() + { + delete component_; + } + + /** + * @brief Allocate storage and initialize the interface mappings. + * + * Identifies the external variables of the wrapped component that are + * connected to the bus and builds the mappings used to transfer residual + * contributions between the component and the interface. Private storage is + * also allocated for the internal variables of the wrapped component. + * + * @return 0 on success and a nonzero value if the component does not contain + * all variables required by the bus interface. + */ + template + int BusPartitionInterface::allocate() + { + + CircuitComponent::allocate(); + + // Build a lookup of the global indices belonging to the bus. + std::unordered_map bus_connections; + + for (size_t i = 0; i < static_cast(bus_->size()); ++i) + { + bus_connections[bus_->getNodeConnection(i).idx_] = i; + } + + const auto& external_indices = component_->getExternIndices(); + + is_external_.assign(static_cast(size_), false); + + bus_input_ports_.clear(); + bus_output_ports_.clear(); + + size_t external_index = 0; + + // Identify which component variables are external and which of those external + // variables are connected to the bus. bus_input_ports_ stores the corresponding + // variable position in this interface, while bus_output_ports_ stores its position + // in the wrapped component's external residual vector. + for (size_t i = 0; i < static_cast(size_); ++i) + { + const IdxT connection_index = component_->getNodeConnection(i); + + this->setConnectionNodes(i, connection_index); + + if (!external_indices.contains(static_cast(i))) + { + continue; + } + + is_external_[i] = true; + + if (bus_connections.contains(connection_index)) + { + bus_input_ports_.push_back(i); + bus_output_ports_.push_back(external_index); + } + + ++external_index; + } + + // A valid bus interface must contain every variable belonging to the bus. + // If fewer bus variables are found, the wrapped component does not provide + // the complete coupling required by this interface. + const size_t bus_size = static_cast(bus_->size()); + + if (bus_input_ports_.size() != bus_size) + { + GridKit::Utilities::Logger::error() << "ERROR: Invalid partition interface detected. " + << "Bus(ID=" << bus_->busID() + << "), Component(ID=" << component_->getIDcomponent() + << "). Expected " << bus_size + << " bus connections, but found " + << bus_input_ports_.size() << "." + << std::endl; + + return 1; + } + + // The wrapped component is evaluated independently by the interface. + // Allocate storage for its internal variables and residuals and redirect + // its internal pointers to this private storage. + const size_t internal_size = static_cast(component_->getInternalSize()); + + component_y_int_ = std::make_unique(internal_size); + component_yp_int_ = std::make_unique(internal_size); + component_f_int_ = std::make_unique(internal_size); + + component_->setInternalPointer(component_y_int_.get()); + component_->setInternalDerivativePointer(component_yp_int_.get()); + component_->setInternalResidualPointer(component_f_int_.get()); + + component_f_ext_ = std::make_unique(component_->getExternSize()); + + return 0; + } + + template + int BusPartitionInterface::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Initialize the partition interface. + * + * @return 0 on success. + */ + 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; + } + + /** + * @brief Evaluate the wrapped component's contributions to the bus residual. + * + * The wrapped component is evaluated using state information supplied through + * the interface. Residual contributions associated with bus variables are + * extracted from the component's external residual and accumulated into the + * corresponding interface residual entries. + * + * @return 0 on success, or the error code returned by the wrapped component. + */ + template + int BusPartitionInterface::evaluateExternalResidual() + { + std::fill_n(component_f_ext_.get(), component_->getExternSize(), ScalarT{}); + + updateComponentPointers(); + + if (int err_code = component_->evaluateExternalResidual()) + { + return err_code; + } + + // Only contributions associated with the bus are accumulated + // into the interface residual below. + for (size_t i = 0; i < bus_input_ports_.size(); ++i) + { + *f_ext_[bus_input_ports_[i]] += component_f_ext_[bus_output_ports_[i]]; + } + + return 0; + } + + /** + * @brief Evaluate the Jacobian contributions associated with the bus. + * + * Evaluates the Jacobian of the wrapped component and extracts the entries + * whose residual rows correspond to variables belonging to the connected bus. + * The selected entries are then used to assemble the interface Jacobian. + * + * @return 0 on success. + */ + template + int BusPartitionInterface::evaluateJacobian() + { + + this->zeroJacMatrix(); + + updateComponentPointers(); + + component_->evaluateJacobian(); + + const IdxT* cooRows = component_->jacobianCooRows(); + const IdxT* cooCols = component_->jacobianCooCols(); + const RealT* cooVals = component_->jacobianCooValues(); + + std::vector rows; + std::vector cols; + std::vector vals; + + rows.reserve(jac_map_.size()); + cols.reserve(jac_map_.size()); + vals.reserve(jac_map_.size()); + + // Extract only the Jacobian entries whose residual rows belong to the bus. + for (const IdxT index : jac_map_) + { + rows.push_back(cooRows[index]); + cols.push_back(cooCols[index]); + vals.push_back(cooVals[index]); + } + + this->setJacValues(rows, cols, vals); + + return 0; + } + + /** + * @brief Update the wrapped component with state and output residual pointers. + * + * Reconstructs the state required to evaluate the wrapped component using + * data supplied through the interface. External component variables are + * connected directly to the interface data, while internal component + * variables are copied into the interface's private storage, which the wrapped + * component already keeps track of for its internal state. + * + * @pre The interface must be allocated and its external state pointers must + * reference valid state and derivative data. + * + * @post The wrapped component is configured with the state, derivative, and + * residual pointers required for evaluation. + * + * @param residual Storage for the component's external residual contributions, + * or nullptr when residual output is not required. + * + * @return 0 on success. + */ + template + int BusPartitionInterface::updateComponentPointers() + { + size_t internal_index = 0; + size_t external_index = 0; + + // Reconstruct the state expected by the wrapped component from the + // interface's external data. Route the external residual to component_f_ext_. + for (size_t i = 0; i < static_cast(component_->size()); ++i) + { + if (is_external_[i]) + { + ExternalConnection connection{ + .y_ = y_ext_[i], + .yp_ = yp_ext_[i], + .f_ = &component_f_ext_[external_index], + .idx_ = component_->getNodeConnection(i)}; + + component_->setExternalConnectionNodes(i, connection); + ++external_index; + } + else + { + component_y_int_[internal_index] = *y_ext_[i]; + component_yp_int_[internal_index] = *yp_ext_[i]; + ++internal_index; + } + } + + 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..aa832095b --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp @@ -0,0 +1,164 @@ + + +#pragma once + +#include + +#include +#include +#include + +namespace GridKit +{ + /** + * @brief Partition interface that provides a wrapped component's current + * injection contributions to a bus. + * + * A BusPartitionInterface is used when a bus and one of its connected + * components belong to different partitions. The interface represents the + * component within the partition containing the bus and provides the current + * injection contributions that the wrapped component makes to that bus. + * + * The interface has the same number of variables as the wrapped component, + * with a one-to-one correspondence between interface variables and wrapped + * component variables. However, all variables are external from the + * interface's point of view, regardless of whether the corresponding variable + * is internal or external to the wrapped component: + * + * @code + * Interface variables: [ext_, ext_, ext_, ext_, ext_, ext_, ext_] + * | | | | | | | + * Wrapped component: [int_, int_, ext_, int_, ext_, ext_, ext_] + * @endcode + * + * This allows the partition containing the wrapped component to provide all + * state information required to reconstruct and evaluate the component on the + * bus side of the partition boundary. Variables that are internal to the + * wrapped component are stored in private interface storage during evaluation, + * while variables that are external to the wrapped component are connected + * directly to the corresponding interface data. + * + * After evaluating the wrapped component, the interface extracts the residual + * and Jacobian contributions associated with the bus and exposes them to the + * partition containing the bus. + * + * * @note The interface residual contributions are zero for all variables except + * those connected to the bus. Only variables connected to the bus expose + * the wrapped component's contribution to the partition. + * + * @note This interface must be added to the partition containing the bus + * to ensure that the system is properly partitioned. + * + * @tparam ScalarT Scalar type used by the model. + * @tparam IdxT Index type used for variables and connections. + */ + 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_ext_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_ext_; + using CircuitComponent::yp_int_; + using CircuitComponent::tag_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_ext_; + 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 CircuitComponent::connection_nodes_; + + public: + BusPartitionInterface(node_type* bus, component_type* component, IdxT id); + virtual ~BusPartitionInterface(); + + int allocate() final; + int initialize() final; + int tagDifferentiable() final; + int setAbsoluteTolerance(RealT) final; + int evaluateInternalResidual() final; + int evaluateExternalResidual() final; + int evaluateJacobian() final; + int evaluateIntegrand() final; + + int initializeAdjoint() final; + int evaluateAdjointResidual() final; + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand() final; + + private: + int updateComponentPointers(); + + /** + * @brief Circuit component participating in the partition split. + */ + component_type* component_; + + /** + * @brief Bus associated with the partition boundary. + */ + node_type* bus_; + + /** + * @brief Storage for the wrapped component's internal state variables. + */ + std::unique_ptr component_y_int_; + + /** + * @brief Storage for the wrapped component's internal state derivatives. + */ + std::unique_ptr component_yp_int_; + + /** + * @brief Storage for the wrapped component's internal residual values. + */ + std::unique_ptr component_f_int_; + + /** + * @brief Storage for the wrapped component's external residual values. + */ + std::unique_ptr component_f_ext_; + + /** + * @brief Wrapped component ports that receive variables from the bus. + */ + std::vector bus_input_ports_; + + /** + * @brief Wrapped component ports whose residual contributions are exposed + * to the bus. + */ + std::vector bus_output_ports_; + + /** + * @brief Mapping from interface Jacobian entries to wrapped component + * Jacobian entries. + */ + std::vector jac_map_; + + /** + * @brief Cached lookup indicating whether each wrapped component variable is external. + * + * Avoids repeatedly searching `extern_indices_` when determining whether + * a wrapped component variable is internal or external. + */ + std::vector is_external_; + }; +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt new file mode 100644 index 000000000..8d833c6dd --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt @@ -0,0 +1,9 @@ +gridkit_add_library( + power_elec_partition_interfaces + SOURCES BusPartitionInterface.cpp + HEADERS BusPartitionInterface.hpp PartitionInterface.hpp + LINK_LIBRARIES + PUBLIC + GridKit::dense_vector + PUBLIC + GridKit::utilities_logger) diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp b/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp new file mode 100644 index 000000000..e76930e64 --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp @@ -0,0 +1,29 @@ + + +#pragma once + +#include +#include +#include + +namespace GridKit +{ + /*! + * @brief Base class for partition interface components. + * + * A partition interface is simply a component that is designed to aid in + * marking partition boundary to facilitate partion evaluation and + * Jacobian assembly. All partition interface that will be implemented + * in the future should inherit from this class and all common features + * should be factored here. + */ + template + class PartitionInterface : public CircuitComponent + { + public: + PartitionInterface() = default; + + ~PartitionInterface() = default; + }; + +} // namespace GridKit From 9c567c0681c2642b1c4e43103b265b7d4600afcc Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 07:44:41 +0000 Subject: [PATCH 2/9] Add SubsystemModel for partition evaluation --- GridKit/Model/PowerElectronics/CMakeLists.txt | 2 + .../Model/PowerElectronics/SubsystemModel.hpp | 844 ++++++++++++++++++ .../SystemModelPowerElectronics.hpp | 38 +- 3 files changed, 863 insertions(+), 21 deletions(-) create mode 100644 GridKit/Model/PowerElectronics/SubsystemModel.hpp diff --git a/GridKit/Model/PowerElectronics/CMakeLists.txt b/GridKit/Model/PowerElectronics/CMakeLists.txt index c985d69d6..368f89515 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 @@ -29,4 +30,5 @@ install( SystemModelPowerElectronics.hpp NodeBase.hpp ExternalConnection.hpp + SubsystemModel.hpp DESTINATION include/GridKit/Model/PowerElectronics) diff --git a/GridKit/Model/PowerElectronics/SubsystemModel.hpp b/GridKit/Model/PowerElectronics/SubsystemModel.hpp new file mode 100644 index 000000000..9ee1eaee1 --- /dev/null +++ b/GridKit/Model/PowerElectronics/SubsystemModel.hpp @@ -0,0 +1,844 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + + /** + * @brief Represents a subset of a PowerElectronicsModel that can be evaluated + * independently. + * + * A SubsystemModel contains a collection of existing GridKit components and + * nodes taken from a larger system. Variables owned by those components and + * nodes become internal variables of the subsystem. Variables needed by those + * components but owned outside the subsystem become external coupling + * variables. + * + * Components normally store connection indices in the global system indexing. + * During subsystem allocation, these indices are temporarily replaced with a + * contiguous local subsystem indexing so that the subsystem can be evaluated + * like an independent PowerElectronicsModel. + * + * External coupling values must be supplied before residual or Jacobian + * evaluation, either directly through the external-data vectors or through a + * forcing function. + * + * @tparam ScalarT Scalar type used by the model. + * @tparam IdxT Index type used for variable and connection indices. + */ + template + class SubsystemModel : public PowerElectronicsModel + { + public: + struct ForcingData + { + std::vector y; + std::vector yp; + }; + + 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 interface_type = PartitionInterface; + + using SystemModel::abs_tol_; + using SystemModel::allocated_; + using SystemModel::allocateVectors; + using SystemModel::alpha_; + using SystemModel::connection_nodes_; + 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_int_; + 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() + : SystemModel(false) + { + } + + /** + * @brief Default constructor for the system model + * + * @post System model parameters set as default + */ + + SubsystemModel(bool use_jac) + : SystemModel(use_jac) + { + } + + ~SubsystemModel() override + { + for (auto* comp : interfaces_) + { + delete comp; + } + // SubsystemModel does not own components_/nodes_ — they belong to (and are + // deleted by) the parent system model. + components_.clear(); + nodes_.clear(); + } + + /** + * @brief Allocate the subsystem for independent evaluation. + * + * Converts the selected components and nodes from global system indexing to + * local subsystem indexing, allocates storage for the subsystem's internal + * variables, and creates storage for coupling variables that lie outside the + * subsystem. + * + * Component pointers are then redirected to either the subsystem's internal + * vectors or its external coupling-data vectors. Finally, the subsystem + * Jacobian structure is assembled using only entries whose row and column + * belong to internal subsystem variables. + * + * @pre Components and nodes added to the subsystem must already be allocated + * by the parent system. + * + * @post Component and node connection indices use local subsystem indexing, + * and the subsystem is ready for residual and Jacobian evaluation. + * + * @return 0 on success, otherwise an error code returned during allocation. + */ + int allocate() override + { + if (int err_code = buildIndexMappings()) + { + return err_code; + } + + if (int err_code = mapGlobalToLocal()) + { + return err_code; + } + + n_intern_ = internal_map_.size(); + n_extern_ = external_map_.size(); + size_ = n_intern_; + + // Allocate subsystem vectors. + y_ext_data_.resize(n_extern_); + yp_ext_data_.resize(n_extern_); + f_ext_data_.resize(n_extern_); + + connection_nodes_ = std::make_unique(size_); + if (!allocated_) + { + allocateVectors(static_cast(size_), true); + abs_tol_.setToZero(memory::HOST); + } + + tag_.resize(size_); + + 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->setConnectionNodes(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; + } + + { // Start node internal indexing after all component internals for proper KLU ordering + size_t node_internal_idx; + for (node_type* node : nodes_) + { + for (size_t i = 0; i < node->getInternalSize(); i++) + { + node_internal_idx = node->getNodeConnection(i).idx_; + + ExternalConnection node_connection{ + .y_ = y_int_ + node_internal_idx, + .yp_ = yp_int_ + node_internal_idx, + .f_ = f_int_ + node_internal_idx, + .idx_ = static_cast(node_internal_idx)}; + + node->setExternalConnectionNodes(i, node_connection); + } + } + } + + { + // The offset for each component's internal variables in the system vector. + // They start at 0, and are stacked on top of each other. + size_t component_internal_idx = 0; + for (component_type* comp : components_) + { + // Update component internal pointers to their correct offsets + comp->setInternalPointer(&y_int_[component_internal_idx]); + comp->setInternalDerivativePointer(&yp_int_[component_internal_idx]); + comp->setInternalResidualPointer(&f_int_[component_internal_idx]); + + component_internal_idx += comp->getInternalSize(); + + const auto& external_indices = comp->getExternIndices(); + + for (IdxT local_index : external_indices) + { + const IdxT connection_index = comp->getNodeConnection(local_index); + + // A variable can be external from the component's point of view while still + // being owned by this subsystem. In that case, connect the component directly + // to the corresponding entry in the subsystem's internal vectors. + if (connection_index < n_intern_) + { + ExternalConnection connection{ + .y_ = y_int_ + connection_index, + .yp_ = yp_int_ + connection_index, + .f_ = f_int_ + connection_index, + .idx_ = connection_index}; + + comp->setExternalConnectionNodes(local_index, connection); + + continue; + } + + // Otherwise the variable is owned outside this subsystem. Connect the + // component to the subsystem's external coupling-data storage instead. + const IdxT external_offset = connection_index - static_cast(n_intern_); + + ExternalConnection connection{ + .y_ = &y_ext_data_[external_offset], + .yp_ = &yp_ext_data_[external_offset], + .f_ = &f_ext_data_[external_offset], + .idx_ = connection_index}; + + comp->setExternalConnectionNodes(local_index, connection); + } + } + } + + // 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; + + // Evaluate component Jacobians to get sparsity + for (component_type* component : components_) + { + component->evaluateJacobian(); + } + + // Check whether a Jacobian entry belongs to the subsystem Jacobian. + // Only entries whose row and column are both internal subsystem + // variables are retained. + 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 Update the subsystem external state and derivative data. + * + * If a forcing function is provided, evaluate it at the current subsystem + * time and copy the returned coupling values into the external state vectors. + * + * @post y_ext_data_ and yp_ext_data_ contain the external coupling values + * returned by the forcing function, if one is set. + * + * @throws std::runtime_error If the forcing function returns vectors whose + * sizes do not match the subsystem external-data vectors. + * + * @return 0 on success. + */ + int distributeExternalVectors() + { + + if (forcing_function_) + { + const auto forcing = (*forcing_function_)(time_); + + if (forcing.y.size() != y_ext_data_.size() || forcing.yp.size() != yp_ext_data_.size()) + { + throw std::runtime_error( + "SubsystemModel::distributeExternalVectors: forcing function " + "returned vectors with incorrect sizes."); + } + + std::copy(forcing.y.begin(), forcing.y.end(), y_ext_data_.begin()); + std::copy(forcing.yp.begin(), forcing.yp.end(), yp_ext_data_.begin()); + } + + return 0; + } + + /** + * @brief Evaluate Residuals at each component then collect them + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int evaluateInternalResidual() override + { + if (int err_code = distributeExternalVectors()) + { + return err_code; + } + + return SystemModel::evaluateInternalResidual(); + } + + /** + * @brief Creates the system Jacobian representing \f$\alpha dF/dy' + dF/dy\f$ + * + * Updates the CSR Jacobian values using the per-component mappings + * computed during allocate(). + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable + */ + int evaluateJacobian() override + { + if (int err_code = distributeExternalVectors()) + { + return err_code; + } + + return SystemModel::evaluateJacobian(); + } + + /** + * @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 (!component->isAllocated()) + { + throw std::logic_error( + "SubsystemModel::addComponent: cannot add an unallocated component."); + } + + if (connections_are_local_) + { + throw std::logic_error( + "SubsystemModel::addComponent: cannot add component while " + "in local-indexed state. Call release() first."); + } + + 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 (!node->isAllocated()) + { + throw std::logic_error( + "SubsystemModel::addNode: cannot add an unallocated node."); + } + + if (connections_are_local_) + { + throw std::logic_error( + "SubsystemModel::addNode: cannot add node while in " + "local-indexed state. Call release() first."); + } + + SystemModel::addNode(node); + } + + /** + * @brief Add a partition interface to the subsystem. + * + * Adds the partition interface to the subsystem's component list and keeps a + * separate reference to it in the interface list. + * + * @param component Pointer to the interface component to add. + */ + void addInterface(interface_type* component) + { + addComponent(component); + interfaces_.push_back(component); + } + + /** + * @brief Restore the subsystem topology to global system indexing. + * + * Reverses local subsystem connection indices with the original + * global system indices. The local internal/external mappings are + * then cleared so that components or nodes can safely be added or removed. + * + * @pre Component and node connections may use local subsystem indices. + * + * @post Component and node connections use their original global indices, + * the subsystem mappings are cleared, and the subsystem is marked + * unallocated. + * + * @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() const + { + return external_data_indices_; + } + + std::vector& getExternalDataY() + { + return y_ext_data_; + } + + std::vector& getExternalDataYP() + { + return yp_ext_data_; + } + + std::vector& getExternalDataF() + { + return f_ext_data_; + } + + void setForcingFunction(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 Restore local subsystem connection indices to global system indices. + * + * Reverses \ref mapGlobalToLocal(). Internal subsystem indices are translated + * through the subsystem connection-node table, while external subsystem + * indices are translated through external_data_indices_. + * + * The component/node connectivity is unchanged; only the index representation + * is restored. + * + * @pre Component and node connections use local subsystem indices. + * + * @post Component and node connections use their original global indices. + * + * @return 0 on success. + */ + int mapLocalToGlobal() + { + if (!connections_are_local_) + { + return 0; + } + + for (component_type* component : components_) + { + + for (IdxT i = 0; i < component->size(); i++) + { + const IdxT index = component->getNodeConnection(i); + + if (index == neg1_) + { + continue; + } + else if (index < this->getInternalSize()) + { + component->setConnectionNodes(i, this->getNodeConnection(index)); + } + else + { + component->setConnectionNodes(i, external_data_indices_[index - this->getInternalSize()]); + } + } + } + + for (node_type* node : nodes_) + { + + for (IdxT i = 0; i < node->size(); i++) + { + const IdxT index = node->getNodeConnection(i).idx_; + + if (index == neg1_) + { + continue; + } + node->setConnectionNodes(i, this->getNodeConnection(index)); + } + } + + connections_are_local_ = false; + + return 0; + } + + /** + * @brief Replace global connection indices with subsystem-local indices. + * + * Uses the mappings created by \ref buildIndexMappings() to rewrite the connection + * indices stored by every component and node. Variables owned by the subsystem + * use `internal_map_`; variables owned outside the subsystem use `external_map_`. + * + * This changes only the indexing used to identify connections; it does not + * change the physical component/node connectivity. + * + * @pre `internal_map_` and `external_map_` have been constructed from global + * connection indices. + * + * @post All valid component and node connections use local subsystem indices. + * + * @return 0 on success. + */ + int mapGlobalToLocal() + { + if (connections_are_local_) + { + return 0; + } + + for (component_type* component : components_) + { + for (IdxT i = 0; i < component->size(); ++i) + { + const IdxT index = component->getNodeConnection(i); + + if (index == neg1_) + { + continue; + } + + if (internal_map_.contains(index)) + { + component->setConnectionNodes(i, internal_map_.at(index)); + } + else + { + component->setConnectionNodes(i, external_map_.at(index)); + } + } + } + + for (node_type* node : nodes_) + { + for (IdxT i = 0; i < node->size(); ++i) + { + const IdxT index = node->getNodeConnection(i).idx_; + + if (index == neg1_) + { + continue; + } + + node->setConnectionNodes(i, internal_map_.at(index)); + } + } + + connections_are_local_ = true; + + return 0; + } + + /** + * @brief Build the global-to-local variable mappings for the subsystem. + * + * Examines the connection indices of all components and nodes in the + * subsystem and divides the referenced variables into two groups: + * + * - Internal variables are owned by a component or node in this subsystem. + * - External variables are required by a subsystem component but are owned + * outside the subsystem. + * + * Internal variables receive local indices first. External coupling variables + * are then assigned indices immediately after the internal range. This gives + * every variable referenced by the subsystem a unique local index while + * preserving its original global index in the corresponding map. + * + * @pre Component and node connections use global system indices. + * + * @post internal_map_ and external_map_ contain the global-to-local mappings + * needed to convert the subsystem topology to local indexing. + * + * @return 0 on success. + */ + int buildIndexMappings() + { + + if (connections_are_local_) + { + return 0; + } + + internal_map_.clear(); + external_map_.clear(); + + size_t component_internal_idx = 0; + // Pass 1: Add variables owned internally by subsystem components. + for (component_type* comp : components_) + { + const auto& extern_indices = comp->getExternIndices(); + + for (IdxT i = 0; i < comp->size(); i++) + { + const IdxT index = comp->getNodeConnection(i); + + if (index != neg1_ && !extern_indices.contains(i)) + { + internal_map_[index] = component_internal_idx++; + } + } + } + + // Pass 2: Add variables owned by subsystem nodes. + for (node_type* node : nodes_) + { + + for (IdxT i = 0; i < node->size(); i++) + { + const IdxT index = node->getNodeConnection(i).idx_; + + if (index != neg1_) + { + internal_map_[index] = component_internal_idx++; + } + } + } + + // Pass 3: Add component dependencies that are owned outside the subsystem. + size_t component_external_idx = 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; + } + + const IdxT index = comp->getNodeConnection(j); + + if (index != neg1_ && !internal_map_.contains(index) && !external_map_.contains(index)) + { + external_map_[index] = component_external_idx++; + } + } + } + + return 0; + } + + /** + *@brief Maps global system indices to local subsystem indices for internal variables. + */ + std::unordered_map internal_map_; + + /** + * @brief Maps global system indices to local subsystem indices for external variables. + */ + std::unordered_map external_map_; + + /** + * @brief Global system index corresponding to each entry in the external subsystem vectors. + */ + std::vector external_data_indices_; + + /** + * @brief subsystem external State, derivative, and residual vectors. + */ + std::vector y_ext_data_; + + /** + * @brief subsystem external derivative + */ + std::vector yp_ext_data_; + + /** + * @brief subsystem external derivative + */ + std::vector f_ext_data_; + + /** + * @brief Optional forcing function used to provide external subsystem data. + * + * It is continuous function that can be sampled at different times + */ + std::optional forcing_function_; + + /** + * @brief Partition interfaces owned by the subsystem. + * + * These interfaces expose the contributions of components participating + * in the partition split. + */ + std::vector interfaces_; + + /** + * @brief Keeps track of whether the components are in local state or global state + */ + bool connections_are_local_{false}; + + }; // class SubsystemModel + +} // namespace GridKit \ No newline at end of file diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 51c8dcee1..20111ace9 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -40,17 +40,6 @@ namespace GridKit using CircuitComponent::yp_int_; 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 * @@ -116,7 +105,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_) @@ -281,7 +270,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_) @@ -304,7 +293,7 @@ namespace GridKit * stored contiguously in the same order as \ref components_, with node variables at the end. * Each internal variable's tag in the system is set to its tag in the component. */ - int tagDifferentiable() final + int tagDifferentiable() { // Ask all component to tag their differentiables for (size_t i = 0; i < components_.size(); i++) @@ -367,7 +356,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int evaluateInternalResidual() final + int evaluateInternalResidual() override { for (IdxT i = 0; i < size_; i++) { @@ -408,7 +397,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int evaluateJacobian() final + int evaluateJacobian() override { // Zero out values RealT* vals = csr_jac_->getValues(); @@ -430,11 +419,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; } } @@ -514,7 +510,7 @@ namespace GridKit allocated_ = false; } - private: + protected: static constexpr IdxT neg1_ = INVALID_INDEX; std::vector components_; @@ -528,4 +524,4 @@ namespace GridKit }; // class PowerElectronicsModel -} // namespace GridKit +} // namespace GridKit \ No newline at end of file From 5d6c5f5bc2de2944e2260c794c5e4b52b5210185 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 07:46:16 +0000 Subject: [PATCH 3/9] Add HIRES test for SubsystemModel --- .../UnitTests/PowerElectronics/CMakeLists.txt | 12 +- tests/UnitTests/PowerElectronics/README.md | 139 ++++ .../SubsystemModelWithHiresTest.hpp | 775 ++++++++++++++++++ .../runSubsystemModelWithHiresTest.cpp | 12 + 4 files changed, 936 insertions(+), 2 deletions(-) create mode 100644 tests/UnitTests/PowerElectronics/README.md create mode 100644 tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp create mode 100644 tests/UnitTests/PowerElectronics/runSubsystemModelWithHiresTest.cpp diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt index 94cd4fd91..88a0bc7f6 100644 --- a/tests/UnitTests/PowerElectronics/CMakeLists.txt +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -3,6 +3,13 @@ target_link_libraries( test_power_electronics_node PRIVATE GridKit::power_electronics_circuit_node GridKit::testing) +add_executable(test_subsystem_model_with_hires runSubsystemModelWithHiresTest.cpp) +target_link_libraries( + test_subsystem_model_with_hires + PRIVATE GridKit::power_elec_partition_interfaces + GridKit::testing + GridKit::sparse_matrix) + add_executable(test_power_electronics_component_clone runComponentCloneTests.cpp) target_link_libraries( test_power_electronics_component_clone @@ -13,8 +20,9 @@ target_link_libraries( GridKit::testing) add_test(NAME PowerElectronicsNodeTest COMMAND $) -add_test(NAME PowerElectronicsComponentCloneTest - COMMAND $) +add_test(NAME SubsystemModelWithHires COMMAND $) +add_test(NAME PowerElectronicsComponentCloneTest COMMAND $) install(TARGETS test_power_electronics_node RUNTIME DESTINATION bin) +install(TARGETS test_subsystem_model_with_hires RUNTIME DESTINATION bin) install(TARGETS test_power_electronics_component_clone RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/PowerElectronics/README.md b/tests/UnitTests/PowerElectronics/README.md new file mode 100644 index 000000000..be5c30248 --- /dev/null +++ b/tests/UnitTests/PowerElectronics/README.md @@ -0,0 +1,139 @@ +## HIRES Partitioning Test Problem + +The HIRES test problem is a simple ODE system with eight variables. To +demonstrate the partitioning machinery in GridKit, the system is divided into +three components. Equations 1--3 belong to **Component 1**, equations 4--5 +belong to **Component 2**, which is modeled as a bus, and equations 6--8 belong +to **Component 3**. + +> **Note:** HIRES is not a circuit problem. In this example, it is modeled to +> resemble GridKit circuit components so that the existing partitioning +> machinery can be used. The example also provides a simple test problem for +> verifying the order of co-simulation methods. + +The full HIRES system is + +$$ +\begin{aligned} +f_1 &= \frac{dy_1}{dt} +1.71y_1 -0.43y_2 -8.32y_3 -0.0007, \\ +f_2 &= \frac{dy_2}{dt} -1.71y_1 +8.75y_2, \\ +f_3 &= \frac{dy_3}{dt} +10.03y_3 -0.43y_4 -0.035y_5, \\ +f_4 &= \frac{dy_4}{dt} -8.32y_2 -1.71y_3 +1.12y_4, \\ +f_5 &= \frac{dy_5}{dt} +1.745y_5 -0.43y_6 -0.43y_7, \\ +f_6 &= \frac{dy_6}{dt} +280y_6y_8 -0.69y_4 -1.71y_5 + +0.43y_6 -0.69y_7, \\ +f_7 &= \frac{dy_7}{dt} -280y_6y_8 +1.81y_7, \\ +f_8 &= \frac{dy_8}{dt} +280y_6y_8 -1.81y_7. +\end{aligned} +$$ + +For the component representation, equations $f_4$ and $f_5$ are decomposed +to expose the contributions from each component: + +$$ +\begin{aligned} +f_4 +&= \frac{dy_4}{dt} -8.32y_2 -1.71y_3 +1.12y_4 \qquad +&\longrightarrow +\left(\frac{dy_4}{dt}+y_4\right) ++\left(0.1y_4-8.32y_2-1.71y_3\right) ++0.02y_4, +\\[6pt] +f_5 +&= \frac{dy_5}{dt}+1.745y_5-0.43y_6-0.43y_7 \qquad +&\longrightarrow +\left(\frac{dy_5}{dt}+y_5\right) ++0.7y_5 ++\left[0.045y_5-0.43y_6-0.43y_7\right]. +\end{aligned} +$$ + +The decomposition does not change the HIRES system. It separates the terms in +the bus equations (conviniently choosen to be equation 4 and 5) according to the component responsible for each contribution. + +### Component 1 + +Component 1 has three internal equations: + +$$ +\begin{aligned} +f_1 &= \frac{dy_1}{dt} +1.71y_1 -0.43y_2 -8.32y_3 -0.0007, \\ +f_2 &= \frac{dy_2}{dt} -1.71y_1 +8.75y_2, \\ +f_3 &= \frac{dy_3}{dt} +10.03y_3 -0.43y_4 -0.035y_5. +\end{aligned} +$$ + +It also contributes the following terms to the bus equations as its external contribution: + +$$ +\begin{aligned} +f_4^{(1)} &= 0.1y_4 -8.32y_2 -1.71y_3, \\ +f_5^{(1)} &= 0.7y_5. +\end{aligned} +$$ + +### Component 2 (HiresBus) + +Component 2 represents the bus and owns the following contributions to +equations 4 and 5: + +$$ +\begin{aligned} +f_4^{(2)} &= \frac{dy_4}{dt} + y_4, \\ +f_5^{(2)} &= \frac{dy_5}{dt} + y_5. +\end{aligned} +$$ + +The remaining terms in these equations are supplied by the components +connected to the bus. + +### Component 3 + +Component 3 has three internal equations: + +$$ +\begin{aligned} +f_6 &= \frac{dy_6}{dt} -280y_6y_8 +0.69y_4 +1.71y_5 + -0.43y_6 +0.69y_7, \\ +f_7 &= \frac{dy_7}{dt} +280y_6y_8 -1.81y_7, \\ +f_8 &= \frac{dy_8}{dt} -280y_6y_8 +1.81y_7. +\end{aligned} +$$ + +It also contributes the following terms to the bus equations as its external contribution: + +$$ +\begin{aligned} +f_4^{(3)} &= 0.02y_4, \\ +f_5^{(3)} &= 0.045y_5 -0.43y_6 -0.43y_7. +\end{aligned} +$$ + +### Bus Residual Assembly + +The complete bus residuals are obtained by adding the contributions from +Components 1, 2, and 3: + +$$ +f_4 = f_4^{(1)} + f_4^{(2)} + f_4^{(3)}, +$$ + +$$ +f_5 = f_5^{(1)} + f_5^{(2)} + f_5^{(3)}. +$$ + +This reconstruction gives exactly the corresponding equations in the original +HIRES system. + +### Partitioning + +The full HIRES system in component form looks like this: + +```text +Component 1 -------- Component 2 (HiresBus) -------- Component 3 +``` + +The system is divided between **Component 2 (HiresBus)** and **Component 3**, and +a `BusPartitionInterface` is added to the first partition. The resulting +partition residuals are then evaluated independently and are then compared with the full-system residual to verify +that the partitioned evaluation reproduces the original system. diff --git a/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp b/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp new file mode 100644 index 000000000..b0c98b804 --- /dev/null +++ b/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp @@ -0,0 +1,775 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + + /*! + * @brief 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::alpha_; + using CircuitComponent::y_ext_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_ext_; + using CircuitComponent::yp_int_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_ext_; + using CircuitComponent::f_int_; + 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() final + { + CircuitComponent::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1)); + + return 0; + } + + int initialize() final + { + return 0; + } + + int tagDifferentiable() final + { + return 0; + } + + int evaluateInternalResidual() final + { + // 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_ext_[0] + 0.035 * *y_ext_[1]; + + return 0; + } + + int evaluateExternalResidual() + { + // outputs + *f_ext_[0] += 8.32 * y_int_[1] + 1.71 * y_int_[2] - 0.1 * *y_ext_[0]; + *f_ext_[1] += -0.7 * *y_ext_[1]; + + return 0; + } + + int evaluateJacobian() final + { + + 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() final + { + return 0; + } + + int initializeAdjoint() final + { + return 0; + } + + int evaluateAdjointResidual() final + { + return 0; + } + + int evaluateAdjointIntegrand() final + { + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + */ + int setAbsoluteTolerance(RealT rel_tol) final + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + CircuitComponent* clone() const + { + return new HiresComponent1(*this); + } + + private: + NodeT* node_ref_; + }; + + /*! + * @brief Hires Bus Component (Component 2). + * + */ + template + class HiresBus : public CircuitComponent + { + + using RealT = typename CircuitComponent::RealT; + using NodeT = typename PowerElectronics::NodeBase; + + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::alpha_; + using CircuitComponent::y_ext_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_ext_; + using CircuitComponent::yp_int_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_ext_; + using CircuitComponent::f_int_; + 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() final + { + CircuitComponent::allocate(); + + this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0)); + this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1)); + + return 0; + } + + int initialize() final + { + return 0; + } + + int tagDifferentiable() final + { + return 0; + } + + int evaluateInternalResidual() final + { + return 0; + } + + int evaluateExternalResidual() final + { + *f_ext_[0] += -*yp_ext_[0] - *y_ext_[0]; + *f_ext_[1] += -*yp_ext_[1] - *y_ext_[1]; + + return 0; + } + + int evaluateJacobian() final + { + 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() final + { + return 0; + } + + int initializeAdjoint() final + { + return 0; + } + + int evaluateAdjointResidual() final + { + return 0; + } + + int evaluateAdjointIntegrand() final + { + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + */ + int setAbsoluteTolerance(RealT rel_tol) final + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + CircuitComponent* clone() const + { + return new HiresBus(*this); + } + + private: + NodeT* node_ref_; + }; + + /*! + * @brief 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::alpha_; + using CircuitComponent::y_ext_; + using CircuitComponent::y_int_; + using CircuitComponent::yp_ext_; + using CircuitComponent::yp_int_; + using CircuitComponent::abs_tol_; + using CircuitComponent::f_ext_; + using CircuitComponent::f_int_; + 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() + { + + // Internals + f_int_[0] = -yp_int_[0] - 280 * y_int_[0] * y_int_[2] + 0.69 * *y_ext_[0] + 1.71 * *y_ext_[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() + { + // Externals + *f_ext_[0] += -0.02 * *y_ext_[0]; + *f_ext_[1] += -0.045 * *y_ext_[1] + 0.43 * y_int_[0] + 0.43 * y_int_[1]; + + 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 + * + */ + int setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + CircuitComponent* clone() const + { + return new HiresComponent3(*this); + } + + private: + NodeT* node_ref_; + }; + + namespace Testing + { + + /** + * This example assembles the HIRES problem and partitions it to demonstrate + * that the partitioned system can reconstruct the full system. + * + * The HIRES problem consists of eight ODE equations and is divided into three + * components. HiresComponent1 contains equations 1, 2, and 3, + * HiresComponent3 contains equations 6, 7, and 8, and HiresBus contains + * equations 4 and 5. HiresBus plays a role similar to MicrogridBusDQ, where + * other components connected to the bus contribute terms to its equations. + * + * In the equations below, terms enclosed in parentheses ( ) are contributions + * from HiresComponent1, terms enclosed in square brackets [ ] are contributions + * from HiresComponent3, and the remaining terms in equations 4 and 5 belong + * to HiresBus. + *\f[ + * f_1 = dy_1/dt + 1.71y_1 - 0.43y_2 - 8.32y_3 - 0.0007 + * + * f_2 = dy_2/dt - 1.71y_1 + 8.75y_2 + * + * f_3 = dy_3/dt + 10.03y_3 - 0.43y_4 - 0.035y_5 + * + * f_4 = dy_4/dt + y_4 + (0.1y_4 - 8.32y_2 - 1.71y_3) + [0.02y_4] + * + * f_5 = dy_5/dt + y_5 + \left(0.7y_5 \right) + \left[0.045y_5 - 0.43y_6 - 0.43y_7 \right] + * + * f_6 = dy_6/dt - 280y_6y_8 + 0.69y_4 + 1.71y_5 - 0.43y_6 + 0.69y_7 + * + * f_7 = dy_7/dt + 280y_6y_8 - 1.81y_7 + * + * f_8 = dy_8/dt - 280y_6y_8 + 1.81y_7 + *\f] + * + * The assembled system has the following structure: + * + * (HiresComponent1) -------- (HiresBus) -------- (HiresComponent3) + * + * The system is partitioned between HiresBus and HiresComponent3. A bus + * partition interface is introduced in the partition containing HiresBus to + * preserve the contribution of HiresComponent3 to the bus equations. The + * residuals evaluated independently by the partitions are then printed with + * the residual of the full system for eye-ball comparision. + * + * This example will also be useful later for testing the order of accuracy of co-simulation methods. + */ + template + class SubsystemModelWithHiresTests + { + using RealT = typename CircuitComponent::RealT; + using Bus = PowerElectronics::MicrogridBus; + using Subsystem = SubsystemModel; + using System = PowerElectronicsModel; + + public: + /** + * @brief Construct the HIRES reference system and subsystem partitions. + */ + SubsystemModelWithHiresTests() + : system_(new System()), + partition1_(new Subsystem()), + partition2_(new Subsystem()), + comp1_(new HiresComponent1(&bus_, 1)), + bus1_(new HiresBus(&bus_, 2)), + comp3_(new HiresComponent3(&bus_, 3)) + { + + y_ = {1, 2, 3, 4, 5, 6, 7, 8}; + yp_ = {1, 2, 3, 4, 5, 6, 7, 8}; + // --------------------------------------------------------------------- + // Assemble and allocate the monolithic reference system + // --------------------------------------------------------------------- + + system_->addComponent(comp1_); + system_->addComponent(comp3_); + system_->addComponent(bus1_); + system_->addNode(&bus_); + + system_->allocate(); + + auto* system_y = system_->y().getData(); + auto* system_yp = system_->yp().getData(); + + for (size_t i = 0; i < system_->size(); ++i) + { + system_y[i] = y_[i]; + system_yp[i] = yp_[i]; + } + + system_->y().setDataUpdated(); + system_->yp().setDataUpdated(); + + system_->updateTime(1.0, 2.0); + system_->evaluateResidual(); + + // --------------------------------------------------------------------- + // Construct the partition interface + // --------------------------------------------------------------------- + + auto* comp3_copy = new HiresComponent3(*comp3_); + bus_interface_ = new BusPartitionInterface(&bus_, comp3_copy, 4); + + bus_interface_->allocate(); + + // --------------------------------------------------------------------- + // Assemble the subsystem partitions + // --------------------------------------------------------------------- + + partition1_->addComponent(comp1_); + partition1_->addComponent(bus1_); + partition1_->addInterface(bus_interface_); + partition1_->addNode(&bus_); + + partition2_->addComponent(comp3_); + + partition1_->allocate(); + partition2_->allocate(); + partitions_ = {partition1_, partition2_}; + + // Distribute variables to all partitions + distributeVariables(y_, yp_); + } + + ~SubsystemModelWithHiresTests() + { + delete partition1_; + delete partition2_; + delete system_; + } + + void distributeVariables(const std::vector& y, const std::vector& yp) + { + for (auto* partition : partitions_) + { + + // Test the forcing function mechanism. + auto forcing_function = [partition, &y, &yp]([[maybe_unused]] ScalarT t) + { + typename Subsystem::ForcingData forcing_data; + + forcing_data.y.resize(partition->getExternSize()); + forcing_data.yp.resize(partition->getExternSize()); + + for (size_t i = 0; i < partition->getExternSize(); ++i) + { + const auto global_index = partition->getExternalDataIndices()[i]; + + forcing_data.y[i] = y[global_index]; + forcing_data.yp[i] = yp[global_index]; + } + + return forcing_data; + }; + + partition->setForcingFunction(std::move(forcing_function)); + + auto* partition_y = partition->y().getData(); + auto* partition_yp = partition->yp().getData(); + + for (size_t i = 0; i < partition->getInternalSize(); ++i) + { + const auto global_index = partition->getNodeConnection(i); + + partition_y[i] = y[global_index]; + partition_yp[i] = yp[global_index]; + } + + partition->y().setDataUpdated(); + partition->yp().setDataUpdated(); + } + } + + /** + * @brief Verify that the partitioned HIRES residual matches the + * monolithic residual. + */ + TestOutcome residual() + { + TestStatus success = true; + + std::vector partition_residual(system_->size(), 0.0); + + for (auto* partition : partitions_) + { + partition->updateTime(1.0, 2.0); + partition->evaluateResidual(); + + auto* residual = partition->getResidual().getData(); + + for (size_t i = 0; i < partition->getInternalSize(); ++i) + { + partition_residual[partition->getNodeConnection(i)] = residual[i]; + } + + partition->getResidual().setDataUpdated(); + } + + auto* reference_residual = system_->getResidual().getData(); + + RealT max_error = 0.0; + + for (size_t i = 0; i < system_->size(); ++i) + { + double error = std::abs(partition_residual[i] - reference_residual[i]) / std::abs(reference_residual[i] + 1); + max_error = std::max(max_error, error); + } + + std::cout << "max error " << max_error << std::endl; + + success *= max_error <= std::numeric_limits::epsilon(); + + return success.report(__func__); + } + + /** + * @brief Verify that each subsystem Jacobian matches the corresponding + * entries of the monolithic HIRES Jacobian. + */ + TestOutcome jacobian() + { + TestStatus success = true; + + RealT alpha = 2.0; + + /* + * GridKit stores the HIRES variables in component assembly order: + * + * System index: 0 1 2 3 4 5 6 7 + * HIRES index: 0 1 2 5 6 7 3 4 + * Variable: y1 y2 y3 y6 y7 y8 y4 y5 + */ + const std::array sysmodel_to_hires = { + 0, 1, 2, 5, 6, 7, 3, 4}; + + const std::array hires_to_sysmodel = { + 0, 1, 2, 6, 7, 3, 4, 5}; + + const RealT y5 = y_[hires_to_sysmodel[5]]; + const RealT y7 = y_[hires_to_sysmodel[7]]; + + std::array, 8> reference_jac = + {{{-alpha - 1.71, 0.43, 8.32, 0.0, 0.0, 0.0, 0.0, 0.0}, + {1.71, -alpha - 8.75, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, -alpha - 10.03, 0.43, 0.035, 0.0, 0.0, 0.0}, + {0.0, 8.32, 1.71, -alpha - 1.12, 0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0, -alpha - 1.745, 0.43, 0.43, 0.0}, + {0.0, 0.0, 0.0, 0.69, 1.71, -alpha - 280.0 * y7 - 0.43, 0.69, -280.0 * y5}, + {0.0, 0.0, 0.0, 0.0, 0.0, 280.0 * y7, -alpha - 1.81, 280.0 * y5}, + {0.0, 0.0, 0.0, 0.0, 0.0, -280.0 * y7, 1.81, -alpha - 280.0 * y5}}}; + + for (auto* partition : partitions_) + { + partition->updateTime(1.0, alpha); + partition->evaluateJacobian(); + + auto* partition_jac = partition->getCsrJacobian(); + + const auto* row_ptr = partition_jac->getRowData(); + const auto* cols = partition_jac->getColData(); + const auto* vals = partition_jac->getValues(); + + const size_t n = partition->getInternalSize(); + + std::vector> dense_jac(n, std::vector(n, 0.0)); + + // Convert the partition CSR Jacobian to a dense matrix. + for (size_t row = 0; row < n; ++row) + { + for (IdxT k = row_ptr[row]; k < row_ptr[row + 1]; ++k) + { + dense_jac[row][cols[k]] = vals[k]; + } + } + + // Compare with the corresponding entries of the full Jacobian. + for (size_t row = 0; row < n; ++row) + { + const IdxT global_row = partition->getNodeConnection(row); + const IdxT ref_row = sysmodel_to_hires[global_row]; + + for (size_t col = 0; col < n; ++col) + { + const IdxT global_col = partition->getNodeConnection(col); + const IdxT ref_col = sysmodel_to_hires[global_col]; + + success *= std::abs(dense_jac[row][col] - reference_jac[ref_row][ref_col]) <= std::numeric_limits::epsilon(); + } + } + } + + return success.report(__func__); + } + + private: + /// @brief Bus used to define the connection between the HIRES components. + Bus bus_; + + /// @brief Complete monolithic HIRES system. + System* system_; + + /// @brief First subsystem of the partitioned HIRES system. + Subsystem* partition1_; + + /// @brief Second subsystem of the partitioned HIRES system. + Subsystem* partition2_; + + /// @brief First component of the HIRES system. + HiresComponent1* comp1_; + + /// @brief Bus component connecting the two parts of the HIRES system. + HiresBus* bus1_; + + /// @brief Third component of the HIRES system. + HiresComponent3* comp3_; + + /// @brief Partition interface representing the bus connection across the subsystem boundary. + BusPartitionInterface* bus_interface_; + + /// @brief Collection of subsystems comprising the partitioned HIRES system. + std::vector partitions_; + + /// @brief State vector used to initialize and evaluate the HIRES system. + std::vector y_; + + /// @brief State derivative vector used to initialize and evaluate the HIRES system. + std::vector yp_; + }; + + } // namespace Testing + +} // namespace GridKit diff --git a/tests/UnitTests/PowerElectronics/runSubsystemModelWithHiresTest.cpp b/tests/UnitTests/PowerElectronics/runSubsystemModelWithHiresTest.cpp new file mode 100644 index 000000000..2b8ef11b0 --- /dev/null +++ b/tests/UnitTests/PowerElectronics/runSubsystemModelWithHiresTest.cpp @@ -0,0 +1,12 @@ +#include "SubsystemModelWithHiresTest.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + GridKit::Testing::SubsystemModelWithHiresTests test; + + result += test.residual(); + result += test.jacobian(); + + return result.summary(); +} From b8d126f10b798742f628d6f89c600467190a81cb Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 08:00:00 +0000 Subject: [PATCH 4/9] Add protected access to SystemModelPowerElectronics --- .../Model/PowerElectronics/SystemModelPowerElectronics.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 20111ace9..1c07ebf5b 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -22,6 +22,7 @@ namespace GridKit using component_type = CircuitComponent; using node_type = PowerElectronics::NodeBase; + protected: using CircuitComponent::abs_tol_; using CircuitComponent::allocated_; using CircuitComponent::allocateVectors; @@ -270,7 +271,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int initialize() + int initialize() override { // Initialize components for (const auto& component : components_) @@ -293,7 +294,7 @@ namespace GridKit * stored contiguously in the same order as \ref components_, with node variables at the end. * Each internal variable's tag in the system is set to its tag in the component. */ - int tagDifferentiable() + int tagDifferentiable() override { // Ask all component to tag their differentiables for (size_t i = 0; i < components_.size(); i++) From 25ccb26abb153386b3bdf557514448b76347164a Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 01:52:32 +0000 Subject: [PATCH 5/9] Change external data vector to use GridKit vectors --- .../Model/PowerElectronics/SubsystemModel.hpp | 42 ++++++++++++------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SubsystemModel.hpp b/GridKit/Model/PowerElectronics/SubsystemModel.hpp index 9ee1eaee1..479551c40 100644 --- a/GridKit/Model/PowerElectronics/SubsystemModel.hpp +++ b/GridKit/Model/PowerElectronics/SubsystemModel.hpp @@ -1,3 +1,8 @@ +/** + * @file BusPartitionInterface.hpp + * @author Abdourahman Barry (abdourahman@vt.edu) + * + */ #pragma once #include @@ -38,6 +43,8 @@ namespace GridKit * evaluation, either directly through the external-data vectors or through a * forcing function. * + * @todo Find a better name for this class and its base class. + * * @tparam ScalarT Scalar type used by the model. * @tparam IdxT Index type used for variable and connection indices. */ @@ -55,6 +62,7 @@ namespace GridKit protected: using SystemModel = PowerElectronicsModel; + using VectorT = typename SystemModel::VectorT; using RealT = typename CircuitComponent::RealT; using CsrMatrixT = typename CircuitComponent::CsrMatrixT; using component_type = CircuitComponent; @@ -204,6 +212,10 @@ namespace GridKit } { + ScalarT* y_ext_data_ptr = y_ext_data_.getData(); + ScalarT* yp_ext_data_ptr = yp_ext_data_.getData(); + ScalarT* f_ext_data_ptr = f_ext_data_.getData(); + // The offset for each component's internal variables in the system vector. // They start at 0, and are stacked on top of each other. size_t component_internal_idx = 0; @@ -243,9 +255,9 @@ namespace GridKit const IdxT external_offset = connection_index - static_cast(n_intern_); ExternalConnection connection{ - .y_ = &y_ext_data_[external_offset], - .yp_ = &yp_ext_data_[external_offset], - .f_ = &f_ext_data_[external_offset], + .y_ = y_ext_data_ptr + external_offset, + .yp_ = yp_ext_data_ptr + external_offset, + .f_ = f_ext_data_ptr + external_offset, .idx_ = connection_index}; comp->setExternalConnectionNodes(local_index, connection); @@ -388,15 +400,15 @@ namespace GridKit { const auto forcing = (*forcing_function_)(time_); - if (forcing.y.size() != y_ext_data_.size() || forcing.yp.size() != yp_ext_data_.size()) + if (forcing.y.size() != y_ext_data_.getSize() || forcing.yp.size() != yp_ext_data_.getSize()) { throw std::runtime_error( "SubsystemModel::distributeExternalVectors: forcing function " "returned vectors with incorrect sizes."); } - std::copy(forcing.y.begin(), forcing.y.end(), y_ext_data_.begin()); - std::copy(forcing.yp.begin(), forcing.yp.end(), yp_ext_data_.begin()); + std::copy(forcing.y.begin(), forcing.y.end(), y_ext_data_.getData()); + std::copy(forcing.yp.begin(), forcing.yp.end(), yp_ext_data_.getData()); } return 0; @@ -538,17 +550,17 @@ namespace GridKit return external_data_indices_; } - std::vector& getExternalDataY() + VectorT& getExternalDataY() { return y_ext_data_; } - std::vector& getExternalDataYP() + VectorT& getExternalDataYP() { return yp_ext_data_; } - std::vector& getExternalDataF() + VectorT& getExternalDataF() { return f_ext_data_; } @@ -805,19 +817,19 @@ namespace GridKit std::vector external_data_indices_; /** - * @brief subsystem external State, derivative, and residual vectors. + * @brief subsystem external state data. */ - std::vector y_ext_data_; + VectorT y_ext_data_; /** - * @brief subsystem external derivative + * @brief subsystem external state derivative */ - std::vector yp_ext_data_; + VectorT yp_ext_data_; /** - * @brief subsystem external derivative + * @brief subsystem external residual data */ - std::vector f_ext_data_; + VectorT f_ext_data_; /** * @brief Optional forcing function used to provide external subsystem data. From 90776bafcabcf23828196517692d14aed560a987 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 01:53:58 +0000 Subject: [PATCH 6/9] Updated SystemModelPowerElectronics.hpp --- .../PowerElectronics/SystemModelPowerElectronics.hpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 1c07ebf5b..9764010d6 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -14,6 +14,10 @@ namespace GridKit { + /** + * @todo Move all PowerElectronics models into the PowerElectronics namespace + * for consistency. + */ template class PowerElectronicsModel : public CircuitComponent { @@ -425,13 +429,11 @@ namespace GridKit const bool is_internal_entry = row != neg1_ && col != neg1_ && row < n_intern_ && col < n_intern_; - if (!is_internal_entry) + if (is_internal_entry) { - continue; + vals[map_to_csr_[counter]] += v[i]; + ++counter; } - - vals[map_to_csr_[counter]] += v[i]; - ++counter; } } From 08e2b763bffc597b236359156284279dd36a9599 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 01:58:41 +0000 Subject: [PATCH 7/9] Update README.md for Hires problem --- tests/UnitTests/PowerElectronics/README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/UnitTests/PowerElectronics/README.md b/tests/UnitTests/PowerElectronics/README.md index be5c30248..e5b4368b3 100644 --- a/tests/UnitTests/PowerElectronics/README.md +++ b/tests/UnitTests/PowerElectronics/README.md @@ -37,13 +37,13 @@ f_4 &\longrightarrow \left(\frac{dy_4}{dt}+y_4\right) +\left(0.1y_4-8.32y_2-1.71y_3\right) -+0.02y_4, -\\[6pt] ++\left[0.02y_4\right], +\\ \\ f_5 &= \frac{dy_5}{dt}+1.745y_5-0.43y_6-0.43y_7 \qquad &\longrightarrow \left(\frac{dy_5}{dt}+y_5\right) -+0.7y_5 ++\left(0.7y_5\right) +\left[0.045y_5-0.43y_6-0.43y_7\right]. \end{aligned} $$ From c90bae298312e0d7a01fcd069a3f8092a4106731 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 02:01:22 +0000 Subject: [PATCH 8/9] Minor update --- .../PartitionInterface/BusPartitionInterface.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp index aa832095b..1f1109720 100644 --- a/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp +++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp @@ -1,4 +1,8 @@ - +/** + * @file BusPartitionInterface.hpp + * @author Abdourahman Barry (abdourahman@vt.edu) + * + */ #pragma once From 6b976decf5f76eb99615d3e76269f521a382734f Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Thu, 3 Sep 2026 18:57:17 +0000 Subject: [PATCH 9/9] Apply pre-commit fixes --- .../PowerElectronics/PartitionInterface/CMakeLists.txt | 10 +++++----- GridKit/Model/PowerElectronics/SubsystemModel.hpp | 2 +- .../PowerElectronics/SystemModelPowerElectronics.hpp | 2 +- tests/UnitTests/PowerElectronics/CMakeLists.txt | 7 +++---- 4 files changed, 10 insertions(+), 11 deletions(-) diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt index 8d833c6dd..9f013ad67 100644 --- a/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt +++ b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt @@ -2,8 +2,8 @@ gridkit_add_library( power_elec_partition_interfaces SOURCES BusPartitionInterface.cpp HEADERS BusPartitionInterface.hpp PartitionInterface.hpp - LINK_LIBRARIES - PUBLIC - GridKit::dense_vector - PUBLIC - GridKit::utilities_logger) + LINK_LIBRARIES + PUBLIC + GridKit::dense_vector + PUBLIC + GridKit::utilities_logger) diff --git a/GridKit/Model/PowerElectronics/SubsystemModel.hpp b/GridKit/Model/PowerElectronics/SubsystemModel.hpp index 479551c40..9ff9e41fd 100644 --- a/GridKit/Model/PowerElectronics/SubsystemModel.hpp +++ b/GridKit/Model/PowerElectronics/SubsystemModel.hpp @@ -853,4 +853,4 @@ namespace GridKit }; // class SubsystemModel -} // namespace GridKit \ No newline at end of file +} // namespace GridKit diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 9764010d6..9796ee9f1 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -527,4 +527,4 @@ namespace GridKit }; // class PowerElectronicsModel -} // namespace GridKit \ No newline at end of file +} // namespace GridKit diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt index 88a0bc7f6..756278102 100644 --- a/tests/UnitTests/PowerElectronics/CMakeLists.txt +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -6,9 +6,7 @@ target_link_libraries( add_executable(test_subsystem_model_with_hires runSubsystemModelWithHiresTest.cpp) target_link_libraries( test_subsystem_model_with_hires - PRIVATE GridKit::power_elec_partition_interfaces - GridKit::testing - GridKit::sparse_matrix) + PRIVATE GridKit::power_elec_partition_interfaces GridKit::testing GridKit::sparse_matrix) add_executable(test_power_electronics_component_clone runComponentCloneTests.cpp) target_link_libraries( @@ -21,7 +19,8 @@ target_link_libraries( add_test(NAME PowerElectronicsNodeTest COMMAND $) add_test(NAME SubsystemModelWithHires COMMAND $) -add_test(NAME PowerElectronicsComponentCloneTest COMMAND $) +add_test(NAME PowerElectronicsComponentCloneTest + COMMAND $) install(TARGETS test_power_electronics_node RUNTIME DESTINATION bin) install(TARGETS test_subsystem_model_with_hires RUNTIME DESTINATION bin)