From c0eddb4e17a0c0c9f8f7ec1951a8dad4e146e3ba Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 07:36:54 +0000 Subject: [PATCH 01/18] 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 0000000000..d4619f3f16 --- /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 0000000000..aa832095b7 --- /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 0000000000..8d833c6dd6 --- /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 0000000000..e76930e649 --- /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 6b607d4be263d914a14cae432e4133c8b53bdf3a Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 07:44:41 +0000 Subject: [PATCH 02/18] 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 c985d69d6e..368f89515d 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 0000000000..9ee1eaee18 --- /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 51c8dcee18..20111ace90 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 b54e31dac8e81aa58beac7119fd9e605d124333b Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 07:46:16 +0000 Subject: [PATCH 03/18] 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 94cd4fd912..88a0bc7f6a 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 0000000000..be5c30248b --- /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 0000000000..b0c98b8048 --- /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 0000000000..2b8ef11b08 --- /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 81387ddc63fc737794d54b82cb215d28744127e5 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 08:00:00 +0000 Subject: [PATCH 04/18] 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 20111ace90..1c07ebf5b0 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 dc794523d90939cabaf4048ad0b98082b707e3ca Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 01:52:32 +0000 Subject: [PATCH 05/18] 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 9ee1eaee18..479551c404 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 d10f40fcbd12f8cebf363579200b79aafdf7c0e3 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 01:53:58 +0000 Subject: [PATCH 06/18] 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 1c07ebf5b0..9764010d6a 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 c002f61a364248bee4ed5e69d6f4a0f83d71a1e4 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 01:58:41 +0000 Subject: [PATCH 07/18] 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 be5c30248b..e5b4368b3b 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 f74f1f16d3360de210790cb4745551fecc33ad2c Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Fri, 28 Aug 2026 02:01:22 +0000 Subject: [PATCH 08/18] 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 aa832095b7..1f1109720d 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 e23c7babd06d8535974e3c3471a27149da674c3e Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Thu, 3 Sep 2026 18:57:17 +0000 Subject: [PATCH 09/18] 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 8d833c6dd6..9f013ad673 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 479551c404..9ff9e41fd3 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 9764010d6a..9796ee9f16 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 88a0bc7f6a..756278102a 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) From 81873143457211aebe4d4364d087647a50dd6a03 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Tue, 8 Sep 2026 19:24:08 +0000 Subject: [PATCH 10/18] Address leak in test --- .../UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp b/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp index b0c98b8048..71a1aaaefd 100644 --- a/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp +++ b/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp @@ -541,8 +541,7 @@ namespace GridKit // Construct the partition interface // --------------------------------------------------------------------- - auto* comp3_copy = new HiresComponent3(*comp3_); - bus_interface_ = new BusPartitionInterface(&bus_, comp3_copy, 4); + bus_interface_ = new BusPartitionInterface(&bus_, comp3_, 4); bus_interface_->allocate(); From 79287e6993571d719cfa33356ad11282b1f3d5a1 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 08:10:53 +0000 Subject: [PATCH 11/18] Add helper files for partitioned microgrid problems --- .../ExamplesHelper/JacTestHelper.hpp | 281 ++++++++++++++++++ .../ExamplesHelper/PartitionUtilities.hpp | 238 +++++++++++++++ 2 files changed, 519 insertions(+) create mode 100644 examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp create mode 100644 examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp diff --git a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp new file mode 100644 index 0000000000..55889f5658 --- /dev/null +++ b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp @@ -0,0 +1,281 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + /** + * @brief Verify that the subsystem Jacobian is an exact subset of the + * full-system Jacobian to ensure accuracy of subsystem Jacobians. + * + * Each subsystem Jacobian entry should match the corresponding entry in the + * full-system Jacobian after converting local subsystem indices back to global + * indices. The check verifies that: + * + * - every subsystem entry exists in the full-system Jacobian, + * - every full-system entry whose row and column belong to the subsystem + * exists in the subsystem Jacobian, and + * - matching entries agree within the specified tolerance. + * + * @note When the components in a subsystem are evaluated in a different order + * than in the monolithic reference system, a larger tolerance may be + * required to account for floating-point roundoff introduced by the + * different summation order. + * + * @param full_jac Full-system Jacobian. + * @param sub_jac Subsystem Jacobian. + * @param subsystem Subsystem associated with sub_jac. + * @param tolerance Maximum allowed difference between matching entries. + * + * @return true if the subsystem Jacobian is an exact subset of the + * full-system Jacobian, false otherwise. + */ + + template + bool verifySubsystemJacobian( + GridKit::LinearAlgebra::CsrMatrix& full_jac, + GridKit::LinearAlgebra::CsrMatrix& sub_jac, + GridKit::SubsystemModel& subsystem, + std::optional tolerance = std::nullopt) + { + 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-to-local indices. + const auto& global_to_local = subsystem.getInternalMap(); + + // The internal map should contain one entry for every subsystem row. + 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; + } + + // The subsystem Jacobian is expected to be square. + 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 map from local subsystem indices to global indices. + // This is needed to compare subsystem rows and columns with the + // corresponding entries in the full-system Jacobian. + 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) + { + + // Make sure the global index is valid for the full-system Jacobian. + if (global_index >= full_jac.getNumRows()) + { + std::cout << "Invalid global index " + << global_index + << " in subsystem internal map\n"; + + return false; + } + + // Make sure the local index is valid for the subsystem Jacobian. + if (local_index >= num_sub_rows) + { + std::cout << "Invalid local index " + << local_index + << " mapped from global index " + << global_index << '\n'; + + return false; + } + + // Each local index should correspond to only one global index. + 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; + } + + // Make sure every subsystem local index has a 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 with 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 the current subsystem row using global column indices. + * + * The subsystem Jacobian uses local column indices, while the full-system + * Jacobian uses global column indices. Convert each local column to its + * corresponding global column so the two rows can be compared directly. + */ + 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]; + + // Fail if the subsystem column index is outside the valid local range. + 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]; + + // Fail if the subsystem row contains more than one entry for the same column. + 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 full-system entries whose columns belong to the subsystem. + for (IdxT full_index = full_begin; full_index < full_end; ++full_index) + { + const IdxT global_column = full_cols[full_index]; + + // No need to check if column belongs outside the subsystem. + if (global_to_local.find(global_column) == global_to_local.end()) + { + continue; + } + + auto sub_entry = sub_row_entries.find(global_column); + + // Then it must exist in the subsystem Jacobian; fail otherwise. + 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); + + // Different component evaluation orders can change the order of floating-point + // summation and introduce small roundoff differences. Use an appropriate + // tolerance when comparing against the monolithic reference. + // if tolerance is not supplied we simply use default machine precision provided + // by GridKit's Test::isEqual + auto isEqual = [&tolerance](RealT value, RealT reference) + { + if (tolerance) + { + return GridKit::Testing::isEqual(value, reference, *tolerance); + } + + return GridKit::Testing::isEqual(value, reference); + }; + + // Then the values must agree, fail otherwise + if (!isEqual(sub_value, full_value)) + { + std::cout << "Jacobian value mismatch at (" + << global_row << ", " + << global_column << "): " + << "full = " << full_value + << ", subsystem = " << sub_value + << ", difference = " << difference << '\n'; + + matches = false; + } + + // Remove matched entry + sub_row_entries.erase(sub_entry); + } + + // Any entries remaining must be missing from the + // full-system Jacobian, so this subsystem row contains incorrect entries, fail! + 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 diff --git a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp new file mode 100644 index 0000000000..d55119da4c --- /dev/null +++ b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp @@ -0,0 +1,238 @@ +#pragma once + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include "MicrogridNetwork.hpp" + +namespace GridKit +{ + + /** + * @brief Partition a scaled microgrid network into subsystem models. + * + * Divides the network into contiguous groups of IBRs and creates one + * subsystem for each group. The IBRs are distributed as evenly as possible, + * with any remainder assigned to the first partitions. + * + * The reference signal node is added to the first subsystem. A partition + * interface is added at the right boundary of every subsystem except the + * last to represent its connection to the neighboring subsystem. + * + * @param[in,out] network Scaled microgrid network to partition. + * @param[out] subsystems Vector populated with the created subsystem models. + * @param[in] num_partitions Number of subsystems to create. + * + * @pre The network has been constructed and contains @c 2*N_size IBRs. + * @pre @p num_partitions is greater than zero and does not exceed the + * number of IBRs in the network. + * @pre The generators, buses, virtual DQ buses, loads, and lines referenced + * by the network have valid lifetimes for use by the created subsystems. + * + * @post @p subsystems contains exactly @p num_partitions subsystem models. + * @post Every IBR in the network belongs to exactly one subsystem. + * @post Every subsystem except the last has a partition interface at its + * right boundary. + * + * @note The subsystem models are dynamically allocated. The caller is + * responsible for releasing and deleting them. + */ + template + void partitionNetwork( + ScaleMicrogridNetwork& network, + std::vector*>& subsystems, + size_t num_partitions) + { + const size_t num_ibrs = 2 * network.N_size; + + assert(num_partitions <= num_ibrs); + + subsystems.resize(num_partitions); + + IdxT q = num_ibrs / num_partitions; + IdxT r = num_ibrs % num_partitions; + IdxT index = 0; + + for (IdxT j = 0; j < num_partitions; j++) + { + auto* partition = + new GridKit::SubsystemModel(); + + // Add the reference signal node to the first partition. + if (j == 0) + { + partition->addNode(&network.dg_signal); + } + + IdxT part_size = q + (j < r ? 1 : 0); + IdxT end = std::min(index + part_size, num_ibrs); + + // Add all components belonging to this partition. + for (; index < end; ++index) + { + partition->addComponent(network.generators[index]); + partition->addComponent(network.busesDQ[index]); + + if (network.loads[index] != nullptr) + { + partition->addComponent(network.loads[index]); + } + + if (network.lines[index] != nullptr) + { + partition->addComponent(network.lines[index]); + } + + partition->addNode(&network.buses[index]); + } + + // Add the interface at the right boundary of the partition. + if (index < num_ibrs) + { + + auto* busInterface = new GridKit::BusPartitionInterface( + &network.buses[index - 1], + network.lines[index], + network.model_id_next++); + + busInterface->allocate(); + partition->addInterface(busInterface); + } + + subsystems[j] = partition; + } + } + + /** + * @brief Evaluate subsystem residuals and reconstruct the global residual. + * + * Distributes the global state and state-derivative vectors to each + * subsystem using its internal and external index mappings. Each subsystem + * residual is then evaluated in parallel, and its internal residual entries + * are gathered into the global residual vector. + * + * @param[in] subsystems Subsystem models to evaluate. + * @param[in] y Global state vector. + * @param[in] yp Global state-derivative vector. + * @param[out] f Global residual vector reconstructed from the subsystem + * residuals. + * @param[in] time Current simulation time. + * @param[in] alpha Jacobian scaling parameter associated with the current + * time-integration evaluation. + * + * @pre All subsystem models have been allocated. + * @pre @p y and @p yp contain all global entries + * @pre @p f is large enough to contain every global residual entry referenced + * by the subsystems. + * @pre The subsystems have disjoint internal index sets so that parallel + * writes to @p f do not overlap. + * + * @post Each subsystem residual has been evaluated at the supplied + * @p time and @p alpha. + * @post @p f contains the reconstructed residual for all internal variables + * represented by the subsystems. + */ + template + void evaluatePartitionResiduals( + const std::vector*>& subsystems, + const std::vector& y, + const std::vector& yp, + std::vector& f, + ScalarT time, + ScalarT alpha) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(guided) +#endif + for (auto* partition : subsystems) + { + partition->updateTime(time, alpha); + + 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(); + + auto* residual = partition->getResidual().getData(); + + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + f[partition->getNodeConnection(i)] = residual[i]; + } + } + } + + /** + * @brief Assemble the scaled microgrid from left to right. + * + * Adds the network components to @p sys_model in physical left-to-right + * order. For each bus location, the generator, virtual DQ bus, load, line, + * and bus associated with that location are added before moving to the next + * location in the network. + * + * This ordering is useful when comparing the monolithic system with + * partition-based evaluations that traverse the network from left to right. + * + * @param[in] network Constructed scaled microgrid network. + * @param[in,out] sys_model Power electronics model to populate. + * + * @pre @c network.N_size is greater than zero. + * @pre @p network has been constructed by buildScaleMicrogridNetwork(). + * @pre All component and node pointers stored in @p network are valid. + * + * @post All nodes and components in @p network have been added to + * @p sys_model in left-to-right network order. + */ + template + void assembleSystemLeftToRight( + ScaleMicrogridNetwork& network, + GridKit::PowerElectronicsModel& sys_model) + { + const size_t num_ibrs = 2 * network.N_size; + + assert(network.N_size > 0); + + sys_model.addNode(&network.dg_signal); + + for (IdxT i = 0; i < num_ibrs; ++i) + { + sys_model.addComponent(network.generators[i]); + sys_model.addComponent(network.busesDQ[i]); + + if (network.loads[i] != nullptr) + { + sys_model.addComponent(network.loads[i]); + } + + if (network.lines[i] != nullptr) + { + sys_model.addComponent(network.lines[i]); + } + + sys_model.addNode(&network.buses[i]); + } + } +} // namespace GridKit From a92322fdf1c98d780e755a8d6951276f8071fed1 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 08:13:17 +0000 Subject: [PATCH 12/18] Add partitioned microgrid examples --- examples/PowerElectronics/CMakeLists.txt | 1 + .../PowerElectronics/Partition/CMakeLists.txt | 39 ++ .../Partition/PartitionMicrogrid.cpp | 298 ++++++++++++ .../Partition/PartitionScaleMicrogrid.cpp | 428 ++++++++++++++++++ 4 files changed, 766 insertions(+) create mode 100644 examples/PowerElectronics/Partition/CMakeLists.txt create mode 100644 examples/PowerElectronics/Partition/PartitionMicrogrid.cpp create mode 100644 examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp diff --git a/examples/PowerElectronics/CMakeLists.txt b/examples/PowerElectronics/CMakeLists.txt index 90592fece9..22b8b86883 100644 --- a/examples/PowerElectronics/CMakeLists.txt +++ b/examples/PowerElectronics/CMakeLists.txt @@ -11,5 +11,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 0000000000..cd182dcd34 --- /dev/null +++ b/examples/PowerElectronics/Partition/CMakeLists.txt @@ -0,0 +1,39 @@ +find_package(OpenMP) + +add_executable(PartitionMicrogrid PartitionMicrogrid.cpp) +target_include_directories( + PartitionMicrogrid + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..) + +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_executable(PartitionScaleMicrogrid PartitionScaleMicrogrid.cpp) +target_include_directories( + PartitionScaleMicrogrid + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..) + +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) + +if(OpenMP_CXX_FOUND) + target_link_libraries(PartitionScaleMicrogrid PRIVATE OpenMP::OpenMP_CXX) +endif() + +add_test(NAME PartitionMicrogrid COMMAND PartitionMicrogrid) +add_test(NAME PartitionScaleMicrogrid COMMAND PartitionScaleMicrogrid) + +install(TARGETS PartitionMicrogrid RUNTIME DESTINATION bin) +install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp new file mode 100644 index 0000000000..c9a0470d5f --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -0,0 +1,298 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "PowerElectronicsExamplesHelper/JacTestHelper.hpp" +#include "PowerElectronicsExamplesHelper/MicrogridNetwork.hpp" +#include "PowerElectronicsExamplesHelper/PartitionUtilities.hpp" + +using Component = GridKit::CircuitComponent; +using Node = GridKit::PowerElectronics::NodeBase; +using Subsystem = GridKit::SubsystemModel; +using System = GridKit::PowerElectronicsModel; + +std::vector getComponentConnections(const std::vector& components); +std::vector getNodeConnections(const std::vector& nodes); + +/* + * Verify partitioned residual and Jacobian evaluation against the + * monolithic microgrid model. + * + * The microgrid is manually divided into two subsystems. A partition interface + * is introduced at the boundary between bus 1 and the line connecting buses 1 + * and 2. Each subsystem is evaluated independently, and its residuals are gathered + * and compared with the monolithic reference. + */ +int main() +{ + + constexpr size_t num_network_sections = 2; + constexpr double time = 0.1; + constexpr double alpha = 0.1; + + bool use_jac = true; + + // --------------------------------------------------------------------------- + // build the grid network and assemble the system model + // --------------------------------------------------------------------------- + GridKit::ScaleMicrogridNetwork network(num_network_sections); + + GridKit::buildScaleMicrogridNetwork(network); + + auto* system = new System(use_jac); + + GridKit::assembleSystemLeftToRight(network, *system); + system->allocate(); + + std::vector y(system->size()); + std::vector yp(system->size()); + + for (size_t i = 0; i < system->size(); ++i) + { + y[i] = static_cast(i + 1); + yp[i] = static_cast(i + 1); + } + + 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(time, alpha); + system->evaluateResidual(); + system->evaluateJacobian(); + + auto* system_jacobian = system->getCsrJacobian(); + auto* system_residual = system->getResidual().getData(); + + //------------------------------------------------------------------------------ + // Gather all global indices belonging partition 1 and 2 to test release() later + //------------------------------------------------------------------------------ + std::vector components = { + network.generators[0], + network.generators[1], + network.lines[1], + network.loads[0], + network.busesDQ[0], + network.busesDQ[1], + network.generators[2], + network.generators[3], + network.lines[2], + network.lines[3], + network.loads[2], + network.busesDQ[2], + network.busesDQ[3]}; + + std::vector nodes = { + &network.dg_signal, + &network.buses[0], + &network.buses[1], + &network.buses[2], + &network.buses[3]}; + + const auto original_component_connections = getComponentConnections(components); + const auto original_node_connections = getNodeConnections(nodes); + + // -------------------------------------------------------------------------------- + // Create 2 Partitions and a Partition interfaces to split to network + // -------------------------------------------------------------------------------- + + auto* bus_interface = new GridKit::BusPartitionInterface( + &network.buses[1], + network.lines[2], + 14); + + if (int err = bus_interface->allocate()) + { + return err; + } + + auto* partition1 = new Subsystem(); + auto* partition2 = new Subsystem(); + + // -------------------------------------------------------------------------------- + // Manually add components, nodes and a bus partition interface to Partition 1 + // -------------------------------------------------------------------------------- + + partition1->addNode(&network.dg_signal); + partition1->addComponent(network.generators[0]); + partition1->addComponent(network.busesDQ[0]); + partition1->addComponent(network.loads[0]); + partition1->addNode(&network.buses[0]); + partition1->addComponent(network.lines[1]); + partition1->addComponent(network.generators[1]); + partition1->addComponent(network.busesDQ[1]); + partition1->addInterface(bus_interface); + partition1->addNode(&network.buses[1]); + + // --------------------------------------------------------------------------- + // Manually add components and nodes to Partition 2 + // --------------------------------------------------------------------------- + + partition2->addComponent(network.generators[2]); + partition2->addComponent(network.busesDQ[2]); + partition2->addComponent(network.loads[2]); + partition2->addComponent(network.lines[2]); + partition2->addNode(&network.buses[2]); + partition2->addComponent(network.generators[3]); + partition2->addComponent(network.busesDQ[3]); + partition2->addComponent(network.lines[3]); + partition2->addNode(&network.buses[3]); + + std::vector partitions = {partition1, partition2}; + + for (auto* partition : partitions) + { + partition->allocate(); + } + + std::vector partition_residual(system->size(), 0.0); + + GridKit::evaluatePartitionResiduals(partitions, y, yp, partition_residual, time, alpha); + + // --------------------------------------------------------------------------- + // Verify the subsystem Jacobians + // --------------------------------------------------------------------------- + + bool jacobians_match = true; + + for (auto* partition : partitions) + { + partition->evaluateJacobian(); + + auto* partition_jacobian = partition->getCsrJacobian(); + + jacobians_match = jacobians_match && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, *partition_jacobian, *partition); + } + + if (!jacobians_match) + { + std::cout << "ERROR: At least one subsystem Jacobian is incorrect!\n"; + return 1; + } + + // --------------------------------------------------------------------------- + // Gather and verify the subsystem residuals against monolithic residuals + // --------------------------------------------------------------------------- + + double max_error = 0.0; + + for (size_t i = 0; i < system->size(); ++i) + { + const double error = std::abs(system_residual[i] - partition_residual[i]) / std::abs(system_residual[i] + 1); + + max_error = std::max(max_error, error); + } + + const double machine_epsilon = + std::numeric_limits::epsilon(); + + const bool residuals_match = max_error <= machine_epsilon; + + std::cout << "\nPartition Microgrid Validation\n"; + std::cout << "------------------------------\n"; + + std::cout << std::left + << std::setw(32) << "Maximum residual error:" + << std::setprecision(16) + << max_error << '\n'; + + std::cout << std::left + << std::setw(32) << "Machine epsilon:" + << machine_epsilon << '\n'; + + std::cout << std::left + << std::setw(32) << "Residuals matched:" + << (residuals_match ? "True" : "False") + << '\n'; + + // --------------------------------------------------------------------------- + // Verify subsystem release() from SubsystemModel + // --------------------------------------------------------------------------- + + for (auto* partition : partitions) + { + partition->release(); + } + + const bool components_restored = getComponentConnections(components) == original_component_connections; + const bool nodes_restored = getNodeConnections(nodes) == original_node_connections; + + if (!components_restored || !nodes_restored) + { + std::cout << "ERROR: Subsystem release did not restore " + "the original global connection indices!\n"; + + return 1; + } + + // --------------------------------------------------------------------------- + // Clean up + // --------------------------------------------------------------------------- + delete system; + + for (auto* partition : partitions) + { + delete partition; + } + + return residuals_match ? 0 : 1; +} + +/** + * @brief Collect the connection indices of a set of components. + */ +std::vector getComponentConnections(const std::vector& components) +{ + std::vector connections; + + for (const auto* component : components) + { + for (size_t i = 0; i < component->size(); ++i) + { + connections.push_back(component->getNodeConnection(i)); + } + } + + return connections; +} + +/** + * @brief Collect the connection indices of a set of nodes. + */ +std::vector getNodeConnections(const std::vector& nodes) +{ + std::vector connections; + + for (auto* node : nodes) + + for (size_t i = 0; i < node->size(); ++i) + { + connections.push_back(node->getNodeConnection(i).idx_); + } + + return connections; +} \ No newline at end of file diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp new file mode 100644 index 0000000000..73590bbf84 --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -0,0 +1,428 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Common/JacTestHelper.hpp" +#include "Common/MicrogridNetwork.hpp" +#include "Common/PartitionUtilities.hpp" + +using index_type = size_t; +using real_type = double; + +using SignalNode = GridKit::PowerElectronics::SignalNode; +using Bus = GridKit::PowerElectronics::MicrogridBus; +using BusDQ = GridKit::MicrogridBusDQ; +using DGGenerator = GridKit::DistributedGenerator; +using Line = GridKit::MicrogridLine; +using Load = GridKit::MicrogridLoad; + +/* + * Output data for each parallel function evaluation run + */ +struct RunResult +{ + bool success = true; + index_type num_partitions; + real_type partition_eval_time; // seconds + real_type monolithic_eval_time; // seconds + real_type speedup; // monolithic / partition + real_type max_error; + std::string jacobian_status; +}; + +/* + * Reference data from the monolithic system evaluation. + * + * The monolithic system is evaluated once and its states, residual, and + * evaluation time are stored here. Each partition configuration reuses this + * data for validation and performance comparison. + */ +struct MonolithicReference +{ + std::vector y; + std::vector yp; + std::vector residual; + real_type monolithic_eval_time; +}; + +/** + * @brief Evaluate the monolithic system and store reference data. + * + * Initializes the global state and state-derivative vectors with deterministic + * values, evaluates and times the monolithic residual, evaluates the full + * Jacobian, and stores the resulting data for later comparison with + * partitioned evaluations. + * + * The residual evaluation time is measured independently from the Jacobian + * evaluation so that the stored timing represents only the cost of the + * monolithic residual evaluation. + * + * @param[in,out] system Allocated monolithic power electronics system. + * + * @return MonolithicReference containing the state vector, state-derivative + * vector, residual, and residual evaluation time. + * + * @pre @p system is non-null. + * @pre @p system has already been assembled and allocated. + * + * @post The state and state-derivative vectors of @p system contain the + * deterministic reference values. + * @post The monolithic residual and Jacobian have been evaluated. + * @post The returned reference contains an independent copy of the + * monolithic residual. + */ +MonolithicReference evaluateMonolithicSystem(GridKit::PowerElectronicsModel* system) +{ + MonolithicReference reference; + + // --------------------------------------------------------------------------- + // Initialize the reference state and state-derivative vectors + // --------------------------------------------------------------------------- + reference.y.resize(system->size()); + reference.yp.resize(system->size()); + + for (size_t i = 0; i < system->size(); ++i) + { + reference.y[i] = static_cast(i + 1); + reference.yp[i] = static_cast(i + 1); + } + + // Copy the reference values into the vectors owned by the system model. + auto* system_y = system->y().getData(); + auto* system_yp = system->yp().getData(); + + for (size_t i = 0; i < system->size(); ++i) + { + system_y[i] = reference.y[i]; + system_yp[i] = reference.yp[i]; + } + + system->y().setDataUpdated(); + system->yp().setDataUpdated(); + + // --------------------------------------------------------------------------- + // Evaluate and time the monolithic residual + // --------------------------------------------------------------------------- + + auto start_time = std::chrono::high_resolution_clock::now(); + + system->updateTime(0.1, 0.1); + system->evaluateResidual(); + + auto end_time = std::chrono::high_resolution_clock::now(); + + reference.monolithic_eval_time = std::chrono::duration(end_time - start_time).count(); + + // evaluate the monolithic Jacobian + system->evaluateJacobian(); + + // Store a copy of the residual return reference data + auto* residual = system->getResidual().getData(); + + reference.residual.assign(residual, residual + system->size()); + + return reference; +} + +/** + * @brief Evaluate one partition configuration and validate it against the + * monolithic reference system. + * + * Creates the requested subsystem decomposition, allocates each subsystem, + * evaluates the partitioned residual, verifies each subsystem Jacobian against + * the monolithic Jacobian, and compares the reconstructed global residual with + * the stored monolithic residual. + * + * The partitioned residual evaluation is timed independently and compared with + * the previously measured monolithic residual evaluation time to compute the + * resulting speedup. + * + * @param[in,out] network Scaled microgrid network to partition. + * @param[in] system Allocated monolithic system used as the reference. + * @param[in] reference Stored monolithic state, derivative, residual, and + * timing information. + * @param[in] num_partitions Number of subsystem partitions to create. + * + * @return Performance and validation results for the requested partition count. + * + * @pre @p system is non-null and has already been evaluated. + * @pre @p network has been constructed using buildScaleMicrogridNetwork(). + * @pre @p num_partitions is greater than zero and does not exceed the number + * of IBRs in the network. + * @pre @p reference contains state and derivative vectors consistent with the + * size of @p system. + * + * @post All temporary subsystem models created by this function are released + * and deleted before returning. + * @post The returned result reports the partition timing, speedup, residual + * error, and subsystem Jacobian validation status. + */ +RunResult evaluatePartitioning( + GridKit::ScaleMicrogridNetwork& network, + GridKit::PowerElectronicsModel* system, + const MonolithicReference& reference, + index_type num_partitions) +{ + + // --------------------------------------------------------------------------- + // Create and allocate subsystem partitions + // --------------------------------------------------------------------------- + + std::vector*> subsystems; + + GridKit::partitionNetwork(network, subsystems, num_partitions); + + for (auto* partition : subsystems) + { + partition->allocate(); + } + + // Global residual reconstructed from the subsystem residuals. + std::vector f(system->size(), 1.0); + // Elementwise error between the monolithic and reconstructed residuals. + std::vector error(system->size(), 1.0); + + // --------------------------------------------------------------------------- + // Evaluate and time the partitioned residual + // --------------------------------------------------------------------------- + auto start_time = std::chrono::high_resolution_clock::now(); + + GridKit::evaluatePartitionResiduals(subsystems, reference.y, reference.yp, f, 0.1, 0.1); + + auto end_time = std::chrono::high_resolution_clock::now(); + + auto partition_eval_time = std::chrono::duration(end_time - start_time); + + // --------------------------------------------------------------------------- + // Verify subsystem Jacobians against the monolithic Jacobian + // --------------------------------------------------------------------------- + + auto* system_jacobian = system->getCsrJacobian(); + + bool jacobian_match = true; + + for (auto* partition : subsystems) + { + partition->evaluateJacobian(); + + jacobian_match = jacobian_match && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, *partition->getCsrJacobian(), *partition); + } + + // --------------------------------------------------------------------------- + // Compare the reconstructed and monolithic residuals + // --------------------------------------------------------------------------- + real_type max_error = 0.0; + + for (size_t i = 0; i < system->size(); ++i) + { + error[i] = std::abs(f[i] - reference.residual[i]) / (reference.residual[i] + 1.0); + + if (max_error < error[i]) + { + max_error = error[i]; + } + } + + // --------------------------------------------------------------------------- + // Store performance and validation results + // --------------------------------------------------------------------------- + RunResult result; + + result.num_partitions = num_partitions; + result.partition_eval_time = partition_eval_time.count(); + result.monolithic_eval_time = reference.monolithic_eval_time; + result.speedup = reference.monolithic_eval_time / partition_eval_time.count(); + result.max_error = max_error; + result.jacobian_status = jacobian_match ? "Correct" : "Wrong"; + + // --------------------------------------------------------------------------- + // Check validation results + // --------------------------------------------------------------------------- + if (!jacobian_match) + { + 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!: " << max_error << std::endl; + result.success = false; + } + + // --------------------------------------------------------------------------- + // Release and destroy temporary subsystem models + // --------------------------------------------------------------------------- + + for (auto* partition : subsystems) + { + partition->release(); + delete partition; + } + + return result; +} + +/** + * @brief Benchmark and validate partitioned residual evaluation of a large + * scaled microgrid. It also confirms the correctness of the subsystem + * Jacobian. + * + * Builds a scaled microgrid once, constructs a monolithic reference system, + * and compares several subsystem decompositions of the same physical network. + * + * For each requested partition count, the benchmark measures the partitioned + * residual evaluation time, computes its speedup relative to the monolithic + * evaluation, verifies the reconstructed residual, and checks every subsystem + * Jacobian against the monolithic Jacobian. + * + * The physical network and monolithic system are constructed only once. + * Individual subsystem models are created and destroyed separately for each + * partition count. + * + * @return 0 if all residual and Jacobian validation checks pass; otherwise 1. + */ +int main(int argc, char const* argv[]) +{ + index_type N_size = 5000; + + std::vector num_partitions_list = {500}; + + /* + * If command-line arguments are provided, the first argument specifies + * N_size and all remaining arguments specify partition counts. + * + * Example: + * ./PartitionScaleMicrogrid 5000 10 48 100 500 + */ + if (argc > 1) + { + try + { + N_size = static_cast(std::stoull(argv[1])); + + if (N_size < 1) + { + std::cerr << "ERROR: N_size must be at least 1.\n"; + return 1; + } + + // When N_size is supplied explicitly, at least one partition count + // must also be supplied. + if (argc < 3) + { + std::cerr << "ERROR: At least one partition count must be provided " + << "when N_size is specified.\n"; + return 1; + } + + num_partitions_list.clear(); + + for (int i = 2; i < argc; ++i) + { + index_type num_partitions = static_cast(std::stoull(argv[i])); + + if (num_partitions < 1) + { + std::cerr << "ERROR: Number of partitions must be at least 1.\n"; + return 1; + } + + if (num_partitions > 2 * N_size) + { + std::cerr << "ERROR: Number of partitions (" + << num_partitions + << ") cannot exceed the number of IBRs (" + << 2 * N_size + << ").\n"; + return 1; + } + + num_partitions_list.push_back(num_partitions); + } + } + catch (const std::exception& e) + { + std::cerr + << "ERROR: Invalid command-line argument: " + << e.what() + << "\n"; + + return 1; + } + } + + bool use_jac = true; + + // Build the physical network once. + GridKit::ScaleMicrogridNetwork network(N_size); + GridKit::buildScaleMicrogridNetwork(network); + + // Build, assemble and allocate the monolithic system once. + auto* system = new GridKit::PowerElectronicsModel(use_jac); + + GridKit::assembleSystemLeftToRight(network, *system); + system->allocate(); + + // Evaluate the monolithic reference once. + MonolithicReference reference = evaluateMonolithicSystem(system); + + 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"; + + // Only the partitioned system is rebuilt and evaluated for each partition count. + for (index_type p : num_partitions_list) + { + assert(p <= 2 * N_size); + + // Takes in the network, partition it into p partitions, and perform parallel function eval + RunResult r = evaluatePartitioning(network, system, reference, p); + + if (!r.success) + { + delete system; + return 1; + } + + // Output the results from partition evaluation + std::cout << std::format("{:<16d}{:>14.4f} s{:>16.4f} s{:>11.2f}x{:>14.3e} {:>16s}\n", + r.num_partitions, + r.partition_eval_time, + r.monolithic_eval_time, + r.speedup, + r.max_error, + r.jacobian_status); + } + + delete system; + + return 0; +} From d902f6f1635d071a6e67a6ebf5308331a52d3319 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Wed, 26 Aug 2026 09:01:31 +0000 Subject: [PATCH 13/18] Fixed include path in partitioned examples --- .../Partition/PartitionMicrogrid.cpp | 7 ------- .../Partition/PartitionScaleMicrogrid.cpp | 13 +++---------- 2 files changed, 3 insertions(+), 17 deletions(-) diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp index c9a0470d5f..8eaa7d9d97 100644 --- a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -6,15 +6,8 @@ #include #include -#include -#include -#include -#include -#include -#include #include #include -#include #include #include diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp index 73590bbf84..81491b9bc0 100644 --- a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -8,22 +8,15 @@ #include #include -#include -#include -#include -#include -#include -#include #include #include -#include #include #include #include -#include "Common/JacTestHelper.hpp" -#include "Common/MicrogridNetwork.hpp" -#include "Common/PartitionUtilities.hpp" +#include "PowerElectronicsExamplesHelper/JacTestHelper.hpp" +#include "PowerElectronicsExamplesHelper/MicrogridNetwork.hpp" +#include "PowerElectronicsExamplesHelper/PartitionUtilities.hpp" using index_type = size_t; using real_type = double; From ab6b9b543134c68c1d8cea04e85ffc57f6a8ea17 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Mon, 7 Sep 2026 01:06:12 +0000 Subject: [PATCH 14/18] Update partition utility files --- .../ExamplesHelper/JacTestHelper.hpp | 40 +- .../ExamplesHelper/PartitionUtilities.hpp | 382 +++++++++--------- 2 files changed, 214 insertions(+), 208 deletions(-) diff --git a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp index 55889f5658..e4411e2c2c 100644 --- a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp +++ b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include @@ -43,23 +43,23 @@ namespace GridKit * full-system Jacobian, false otherwise. */ - template + template bool verifySubsystemJacobian( - GridKit::LinearAlgebra::CsrMatrix& full_jac, - GridKit::LinearAlgebra::CsrMatrix& sub_jac, - GridKit::SubsystemModel& subsystem, - std::optional tolerance = std::nullopt) + GridKit::LinearAlgebra::CsrMatrix& full_jac, + GridKit::LinearAlgebra::CsrMatrix& sub_jac, + GridKit::SubsystemModel& subsystem, + ScalarT tolerance = std::numeric_limits::epsilon()) { 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); + const ScalarT* 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* sub_rows = sub_jac.getRowData(host); + const IdxT* sub_cols = sub_jac.getColData(host); + const ScalarT* sub_vals = sub_jac.getValues(host); const IdxT num_sub_rows = sub_jac.getNumRows(); const IdxT num_sub_cols = sub_jac.getNumColumns(); @@ -166,7 +166,7 @@ namespace GridKit * Jacobian uses global column indices. Convert each local column to its * corresponding global column so the two rows can be compared directly. */ - std::unordered_map sub_row_entries; + std::unordered_map sub_row_entries; for (IdxT sub_index = sub_begin; sub_index < sub_end; ++sub_index) { @@ -224,27 +224,17 @@ namespace GridKit 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); + const ScalarT full_value = full_vals[full_index]; + const ScalarT sub_value = sub_entry->second; + const ScalarT difference = std::abs(full_value - sub_value); // Different component evaluation orders can change the order of floating-point // summation and introduce small roundoff differences. Use an appropriate // tolerance when comparing against the monolithic reference. // if tolerance is not supplied we simply use default machine precision provided // by GridKit's Test::isEqual - auto isEqual = [&tolerance](RealT value, RealT reference) - { - if (tolerance) - { - return GridKit::Testing::isEqual(value, reference, *tolerance); - } - - return GridKit::Testing::isEqual(value, reference); - }; - // Then the values must agree, fail otherwise - if (!isEqual(sub_value, full_value)) + if (!GridKit::Testing::isEqual(sub_value, full_value, tolerance)) { std::cout << "Jacobian value mismatch at (" << global_row << ", " diff --git a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp index d55119da4c..0ccc66138f 100644 --- a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp +++ b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp @@ -13,226 +13,242 @@ #include "MicrogridNetwork.hpp" -namespace GridKit +/** + * @brief Partition a scaled microgrid network into subsystem models. + * + * Divides the network into contiguous groups of IBRs and creates one + * subsystem for each group. The IBRs are distributed as evenly as possible, + * with any remainder assigned to the first partitions. + * + * The reference signal node is added to the first subsystem. A partition + * interface is added at the right boundary of every subsystem except the + * last to represent its connection to the neighboring subsystem. + * + * @param[in,out] network Scaled microgrid network to partition. + * @param[out] subsystems Vector populated with the created subsystem models. + * @param[in] num_partitions Number of subsystems to create. + * + * @pre The network has been constructed and contains @c 2*N_size IBRs. + * @pre @p num_partitions is greater than zero and does not exceed the + * number of IBRs in the network. + * @pre The generators, buses, virtual DQ buses, loads, and lines referenced + * by the network have valid lifetimes for use by the created subsystems. + * + * @post @p subsystems contains exactly @p num_partitions subsystem models. + * @post Every IBR in the network belongs to exactly one subsystem. + * @post Every subsystem except the last has a partition interface at its + * right boundary. + * + * @note The subsystem models are dynamically allocated. The caller is + * responsible for releasing and deleting them. + */ +template +void partitionNetwork( + ScaleMicrogridNetwork& network, + std::vector*>& subsystems, + size_t num_partitions) { + const size_t num_ibrs = 2 * network.N_size; - /** - * @brief Partition a scaled microgrid network into subsystem models. - * - * Divides the network into contiguous groups of IBRs and creates one - * subsystem for each group. The IBRs are distributed as evenly as possible, - * with any remainder assigned to the first partitions. - * - * The reference signal node is added to the first subsystem. A partition - * interface is added at the right boundary of every subsystem except the - * last to represent its connection to the neighboring subsystem. - * - * @param[in,out] network Scaled microgrid network to partition. - * @param[out] subsystems Vector populated with the created subsystem models. - * @param[in] num_partitions Number of subsystems to create. - * - * @pre The network has been constructed and contains @c 2*N_size IBRs. - * @pre @p num_partitions is greater than zero and does not exceed the - * number of IBRs in the network. - * @pre The generators, buses, virtual DQ buses, loads, and lines referenced - * by the network have valid lifetimes for use by the created subsystems. - * - * @post @p subsystems contains exactly @p num_partitions subsystem models. - * @post Every IBR in the network belongs to exactly one subsystem. - * @post Every subsystem except the last has a partition interface at its - * right boundary. - * - * @note The subsystem models are dynamically allocated. The caller is - * responsible for releasing and deleting them. - */ - template - void partitionNetwork( - ScaleMicrogridNetwork& network, - std::vector*>& subsystems, - size_t num_partitions) - { - const size_t num_ibrs = 2 * network.N_size; + assert(num_partitions <= num_ibrs); + + subsystems.resize(num_partitions); - assert(num_partitions <= num_ibrs); + IdxT q = num_ibrs / num_partitions; + IdxT r = num_ibrs % num_partitions; - subsystems.resize(num_partitions); + IdxT partition_begin = 0; - IdxT q = num_ibrs / num_partitions; - IdxT r = num_ibrs % num_partitions; - IdxT index = 0; + for (IdxT j = 0; j < num_partitions; j++) + { + auto* partition = new GridKit::SubsystemModel(); - for (IdxT j = 0; j < num_partitions; j++) + // Add the reference signal node to the first partition. + if (j == 0) { - auto* partition = - new GridKit::SubsystemModel(); + partition->addNode(&network.dg_signal); + } - // Add the reference signal node to the first partition. - if (j == 0) - { - partition->addNode(&network.dg_signal); - } + // Divide the IBRs as evenly as possible among the partitions. + // Each partition receives q IBRs, and the first r partitions receive + // one additional IBR to account for any remainder. The partition spans + // [partition_begin, partition_end), where partition_end is the first + // IBR index not included in this partition. + const IdxT partition_size = q + ((j < r) ? 1 : 0); + const IdxT partition_end = partition_begin + partition_size; - IdxT part_size = q + (j < r ? 1 : 0); - IdxT end = std::min(index + part_size, num_ibrs); + // Add all components belonging to this partition. + for (IdxT i = partition_begin; i < partition_end; ++i) + { + partition->addComponent(network.generators[i]); + partition->addComponent(network.busesDQ[i]); - // Add all components belonging to this partition. - for (; index < end; ++index) + if (network.loads[i] != nullptr) { - partition->addComponent(network.generators[index]); - partition->addComponent(network.busesDQ[index]); - - if (network.loads[index] != nullptr) - { - partition->addComponent(network.loads[index]); - } - - if (network.lines[index] != nullptr) - { - partition->addComponent(network.lines[index]); - } - - partition->addNode(&network.buses[index]); + partition->addComponent(network.loads[i]); } - // Add the interface at the right boundary of the partition. - if (index < num_ibrs) + // Each line is owned by the partition containing its right bus. + // line[i - 1] connects bus[i - 1] to bus[i]. + if (i > 0) { + partition->addComponent(network.lines[i - 1]); + } - auto* busInterface = new GridKit::BusPartitionInterface( - &network.buses[index - 1], - network.lines[index], - network.model_id_next++); + partition->addNode(&network.buses[i]); + } - busInterface->allocate(); - partition->addInterface(busInterface); - } + // If another partition exists to the right, the line crossing this + // partition's right boundary is owned by that next partition. Add an + // interface here to the left partition. + if (partition_end < num_ibrs) + { - subsystems[j] = partition; + const IdxT boundary_line_index = partition_end - 1; + + auto* bus_interface = new GridKit::BusPartitionInterface( + &network.buses[boundary_line_index], + network.lines[boundary_line_index], + network.model_id_next++); + + bus_interface->allocate(); + partition->addInterface(bus_interface); } - } - /** - * @brief Evaluate subsystem residuals and reconstruct the global residual. - * - * Distributes the global state and state-derivative vectors to each - * subsystem using its internal and external index mappings. Each subsystem - * residual is then evaluated in parallel, and its internal residual entries - * are gathered into the global residual vector. - * - * @param[in] subsystems Subsystem models to evaluate. - * @param[in] y Global state vector. - * @param[in] yp Global state-derivative vector. - * @param[out] f Global residual vector reconstructed from the subsystem - * residuals. - * @param[in] time Current simulation time. - * @param[in] alpha Jacobian scaling parameter associated with the current - * time-integration evaluation. - * - * @pre All subsystem models have been allocated. - * @pre @p y and @p yp contain all global entries - * @pre @p f is large enough to contain every global residual entry referenced - * by the subsystems. - * @pre The subsystems have disjoint internal index sets so that parallel - * writes to @p f do not overlap. - * - * @post Each subsystem residual has been evaluated at the supplied - * @p time and @p alpha. - * @post @p f contains the reconstructed residual for all internal variables - * represented by the subsystems. - */ - template - void evaluatePartitionResiduals( - const std::vector*>& subsystems, - const std::vector& y, - const std::vector& yp, - std::vector& f, - ScalarT time, - ScalarT alpha) - { + subsystems[j] = partition; + partition_begin = partition_end; + } +} + +/** + * @brief Evaluate subsystem residuals and reconstruct the global residual. + * + * Distributes the global state and state-derivative vectors to each + * subsystem using its internal and external index mappings. Each subsystem + * residual is then evaluated in parallel, and its internal residual entries + * are gathered into the global residual vector. + * + * @param[in] subsystems Subsystem models to evaluate. + * @param[in] y Global state vector. + * @param[in] yp Global state-derivative vector. + * @param[out] f Global residual vector reconstructed from the subsystem + * residuals. + * @param[in] time Current simulation time. + * @param[in] alpha Jacobian scaling parameter associated with the current + * time-integration evaluation. + * + * @pre All subsystem models have been allocated. + * @pre @p y and @p yp contain all global entries + * @pre @p f is large enough to contain every global residual entry referenced + * by the subsystems. + * @pre The subsystems have disjoint internal index sets so that parallel + * writes to @p f do not overlap. + * + * @post Each subsystem residual has been evaluated at the supplied + * @p time and @p alpha. + * @post @p f contains the reconstructed residual for all internal variables + * represented by the subsystems. + */ +template +void evaluatePartitionResiduals( + const std::vector*>& subsystems, + const std::vector& y, + const std::vector& yp, + std::vector& f, + ScalarT time, + ScalarT alpha) +{ #ifdef _OPENMP #pragma omp parallel for schedule(guided) #endif - for (auto* partition : subsystems) + for (auto* partition : subsystems) + { + partition->updateTime(time, alpha); + + auto* external_y = partition->getExternalDataY().getData(); + auto* external_yp = partition->getExternalDataYP().getData(); + + for (size_t i = 0; i < partition->getExternSize(); ++i) { - partition->updateTime(time, alpha); + const auto global_index = partition->getExternalDataIndices()[i]; - for (size_t i = 0; i < partition->getExternSize(); i++) - { - partition->getExternalDataY()[i] = y[partition->getExternalDataIndices()[i]]; - partition->getExternalDataYP()[i] = yp[partition->getExternalDataIndices()[i]]; - } + external_y[i] = y[global_index]; + external_yp[i] = yp[global_index]; + } - auto* partition_y = partition->y().getData(); - auto* partition_yp = partition->yp().getData(); + partition->getExternalDataY().setDataUpdated(); + partition->getExternalDataYP().setDataUpdated(); - for (size_t i = 0; i < partition->getInternalSize(); i++) - { - partition_y[i] = y[partition->getNodeConnection(i)]; - partition_yp[i] = yp[partition->getNodeConnection(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->y().setDataUpdated(); + partition->yp().setDataUpdated(); - partition->evaluateResidual(); + partition->evaluateResidual(); - auto* residual = partition->getResidual().getData(); + auto* residual = partition->getResidual().getData(); - for (size_t i = 0; i < partition->getInternalSize(); i++) - { - f[partition->getNodeConnection(i)] = residual[i]; - } + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + f[partition->getNodeConnection(i)] = residual[i]; } } +} + +/** + * @brief Assemble the scaled microgrid from left to right. + * + * Adds the network components to @p sys_model in physical left-to-right + * order. For each bus location, the generator, virtual DQ bus, load, line, + * and bus associated with that location are added before moving to the next + * location in the network. + * + * This ordering is useful when comparing the monolithic system with + * partition-based evaluations that traverse the network from left to right. + * + * @param[in] network Constructed scaled microgrid network. + * @param[in,out] sys_model Power electronics model to populate. + * + * @pre @c network.N_size is greater than zero. + * @pre @p network has been constructed by buildScaleMicrogridNetwork(). + * @pre All component and node pointers stored in @p network are valid. + * + * @post All nodes and components in @p network have been added to + * @p sys_model in left-to-right network order. + */ +template +void assembleSystemLeftToRight( + ScaleMicrogridNetwork& network, + GridKit::PowerElectronicsModel& sys_model) +{ + const size_t num_ibrs = 2 * network.N_size; - /** - * @brief Assemble the scaled microgrid from left to right. - * - * Adds the network components to @p sys_model in physical left-to-right - * order. For each bus location, the generator, virtual DQ bus, load, line, - * and bus associated with that location are added before moving to the next - * location in the network. - * - * This ordering is useful when comparing the monolithic system with - * partition-based evaluations that traverse the network from left to right. - * - * @param[in] network Constructed scaled microgrid network. - * @param[in,out] sys_model Power electronics model to populate. - * - * @pre @c network.N_size is greater than zero. - * @pre @p network has been constructed by buildScaleMicrogridNetwork(). - * @pre All component and node pointers stored in @p network are valid. - * - * @post All nodes and components in @p network have been added to - * @p sys_model in left-to-right network order. - */ - template - void assembleSystemLeftToRight( - ScaleMicrogridNetwork& network, - GridKit::PowerElectronicsModel& sys_model) - { - const size_t num_ibrs = 2 * network.N_size; + assert(network.N_size > 0); - assert(network.N_size > 0); + sys_model.addNode(&network.dg_signal); - sys_model.addNode(&network.dg_signal); + for (IdxT i = 0; i < num_ibrs; ++i) + { + sys_model.addComponent(network.generators[i]); + sys_model.addComponent(network.busesDQ[i]); - for (IdxT i = 0; i < num_ibrs; ++i) + if (network.loads[i] != nullptr) { - sys_model.addComponent(network.generators[i]); - sys_model.addComponent(network.busesDQ[i]); - - if (network.loads[i] != nullptr) - { - sys_model.addComponent(network.loads[i]); - } - - if (network.lines[i] != nullptr) - { - sys_model.addComponent(network.lines[i]); - } + sys_model.addComponent(network.loads[i]); + } - sys_model.addNode(&network.buses[i]); + if (i > 0) + { + sys_model.addComponent(network.lines[i - 1]); } + + sys_model.addNode(&network.buses[i]); } -} // namespace GridKit +} \ No newline at end of file From e356c31621656c8009d391636cbf427ecaf6996b Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Mon, 7 Sep 2026 01:06:53 +0000 Subject: [PATCH 15/18] Update partition microgrid examples --- .../Partition/PartitionMicrogrid.cpp | 33 +++--- .../Partition/PartitionScaleMicrogrid.cpp | 104 ++++++++++++------ 2 files changed, 86 insertions(+), 51 deletions(-) diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp index 8eaa7d9d97..8b28930e3d 100644 --- a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -35,22 +36,20 @@ std::vector getNodeConnections(const std::vector& nodes); int main() { - constexpr size_t num_network_sections = 2; - constexpr double time = 0.1; - constexpr double alpha = 0.1; + constexpr size_t N_size = 2; + constexpr double time = 0.1; + constexpr double alpha = 0.1; bool use_jac = true; // --------------------------------------------------------------------------- // build the grid network and assemble the system model // --------------------------------------------------------------------------- - GridKit::ScaleMicrogridNetwork network(num_network_sections); - - GridKit::buildScaleMicrogridNetwork(network); + ScaleMicrogridNetwork network(N_size); auto* system = new System(use_jac); - GridKit::assembleSystemLeftToRight(network, *system); + assembleSystemLeftToRight(network, *system); system->allocate(); std::vector y(system->size()); @@ -87,14 +86,14 @@ int main() std::vector components = { network.generators[0], network.generators[1], - network.lines[1], + network.lines[0], network.loads[0], network.busesDQ[0], network.busesDQ[1], network.generators[2], network.generators[3], + network.lines[1], network.lines[2], - network.lines[3], network.loads[2], network.busesDQ[2], network.busesDQ[3]}; @@ -115,7 +114,7 @@ int main() auto* bus_interface = new GridKit::BusPartitionInterface( &network.buses[1], - network.lines[2], + network.lines[1], 14); if (int err = bus_interface->allocate()) @@ -135,7 +134,7 @@ int main() partition1->addComponent(network.busesDQ[0]); partition1->addComponent(network.loads[0]); partition1->addNode(&network.buses[0]); - partition1->addComponent(network.lines[1]); + partition1->addComponent(network.lines[0]); partition1->addComponent(network.generators[1]); partition1->addComponent(network.busesDQ[1]); partition1->addInterface(bus_interface); @@ -148,11 +147,11 @@ int main() partition2->addComponent(network.generators[2]); partition2->addComponent(network.busesDQ[2]); partition2->addComponent(network.loads[2]); - partition2->addComponent(network.lines[2]); + partition2->addComponent(network.lines[1]); partition2->addNode(&network.buses[2]); partition2->addComponent(network.generators[3]); partition2->addComponent(network.busesDQ[3]); - partition2->addComponent(network.lines[3]); + partition2->addComponent(network.lines[2]); partition2->addNode(&network.buses[3]); std::vector partitions = {partition1, partition2}; @@ -164,7 +163,7 @@ int main() std::vector partition_residual(system->size(), 0.0); - GridKit::evaluatePartitionResiduals(partitions, y, yp, partition_residual, time, alpha); + evaluatePartitionResiduals(partitions, y, yp, partition_residual, time, alpha); // --------------------------------------------------------------------------- // Verify the subsystem Jacobians @@ -178,7 +177,11 @@ int main() auto* partition_jacobian = partition->getCsrJacobian(); - jacobians_match = jacobians_match && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, *partition_jacobian, *partition); + jacobians_match = jacobians_match + && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, + *partition_jacobian, + *partition, + 1e-13); } if (!jacobians_match) diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp index 81491b9bc0..b61b99acaf 100644 --- a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -4,7 +4,7 @@ #include #include #include -#include +#include #include #include @@ -93,7 +93,7 @@ MonolithicReference evaluateMonolithicSystem(GridKit::PowerElectronicsModelsize()); reference.yp.resize(system->size()); - for (size_t i = 0; i < system->size(); ++i) + for (index_type i = 0; i < system->size(); ++i) { reference.y[i] = static_cast(i + 1); reference.yp[i] = static_cast(i + 1); @@ -103,7 +103,7 @@ MonolithicReference evaluateMonolithicSystem(GridKit::PowerElectronicsModely().getData(); auto* system_yp = system->yp().getData(); - for (size_t i = 0; i < system->size(); ++i) + for (index_type i = 0; i < system->size(); ++i) { system_y[i] = reference.y[i]; system_yp[i] = reference.yp[i]; @@ -169,8 +169,9 @@ MonolithicReference evaluateMonolithicSystem(GridKit::PowerElectronicsModel RunResult evaluatePartitioning( - GridKit::ScaleMicrogridNetwork& network, + ScaleMicrogridNetwork& network, GridKit::PowerElectronicsModel* system, const MonolithicReference& reference, index_type num_partitions) @@ -182,7 +183,7 @@ RunResult evaluatePartitioning( std::vector*> subsystems; - GridKit::partitionNetwork(network, subsystems, num_partitions); + partitionNetwork(network, subsystems, num_partitions); for (auto* partition : subsystems) { @@ -199,7 +200,7 @@ RunResult evaluatePartitioning( // --------------------------------------------------------------------------- auto start_time = std::chrono::high_resolution_clock::now(); - GridKit::evaluatePartitionResiduals(subsystems, reference.y, reference.yp, f, 0.1, 0.1); + evaluatePartitionResiduals(subsystems, reference.y, reference.yp, f, 0.1, 0.1); auto end_time = std::chrono::high_resolution_clock::now(); @@ -211,13 +212,18 @@ RunResult evaluatePartitioning( auto* system_jacobian = system->getCsrJacobian(); - bool jacobian_match = true; + bool jacobian_match = true; + real_type Jac_tol = 1e-13; for (auto* partition : subsystems) { partition->evaluateJacobian(); - jacobian_match = jacobian_match && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, *partition->getCsrJacobian(), *partition); + jacobian_match = jacobian_match + && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, + *partition->getCsrJacobian(), + *partition, + Jac_tol); } // --------------------------------------------------------------------------- @@ -225,9 +231,9 @@ RunResult evaluatePartitioning( // --------------------------------------------------------------------------- real_type max_error = 0.0; - for (size_t i = 0; i < system->size(); ++i) + for (index_type i = 0; i < system->size(); ++i) { - error[i] = std::abs(f[i] - reference.residual[i]) / (reference.residual[i] + 1.0); + error[i] = std::abs(f[i] - reference.residual[i]) / (std::abs(reference.residual[i]) + 1.0); if (max_error < error[i]) { @@ -277,6 +283,47 @@ RunResult evaluatePartitioning( return result; } +/** + * @brief Print the header for the partition evaluation results table. + */ +void printResultsHeader() +{ + std::cout << std::left + << std::setw(16) << "num_partitions" + << std::right + << std::setw(16) << "partition_time" + << std::setw(18) << "monolithic_time" + << std::setw(12) << "speedup" + << std::setw(14) << "error" + << std::setw(16) << "Jacobians" + << '\n'; + + std::cout << std::string(92, '-') << '\n'; +} + +/** + * @brief Print one row of partition evaluation results. + * + * @param[in] result Results from the partitioned system evaluation. + */ +void printResult(const RunResult& result) +{ + std::cout << std::left + << std::setw(16) << result.num_partitions + << std::right + << std::setw(14) << std::fixed << std::setprecision(4) + << result.partition_eval_time << " s" + << std::setw(16) << result.monolithic_eval_time << " s" + << std::setw(11) << std::setprecision(2) + << result.speedup << "x" + << std::setw(14) << std::scientific << std::setprecision(3) + << result.max_error << ' ' + << std::setw(16) << result.jacobian_status + << '\n'; + + std::cout << std::defaultfloat; +} + /** * @brief Benchmark and validate partitioned residual evaluation of a large * scaled microgrid. It also confirms the correctness of the subsystem @@ -366,38 +413,29 @@ int main(int argc, char const* argv[]) } } - bool use_jac = true; + const bool use_jac = true; // Build the physical network once. - GridKit::ScaleMicrogridNetwork network(N_size); - GridKit::buildScaleMicrogridNetwork(network); + ScaleMicrogridNetwork network(N_size); // Build, assemble and allocate the monolithic system once. auto* system = new GridKit::PowerElectronicsModel(use_jac); - GridKit::assembleSystemLeftToRight(network, *system); + assembleSystemLeftToRight(network, *system); system->allocate(); // Evaluate the monolithic reference once. - MonolithicReference reference = evaluateMonolithicSystem(system); - - std::cout << std::format("{:<16}{:>16}{:>18}{:>12}{:>14}{:>16}\n", - "num_partitions", - "partition_time", - "monolithic_time", - "speedup", - "error", - "Jacobians"); + const MonolithicReference reference = evaluateMonolithicSystem(system); - std::cout << std::string(93, '-') << "\n"; + printResultsHeader(); // Only the partitioned system is rebuilt and evaluated for each partition count. - for (index_type p : num_partitions_list) + for (const index_type p : num_partitions_list) { assert(p <= 2 * N_size); - // Takes in the network, partition it into p partitions, and perform parallel function eval - RunResult r = evaluatePartitioning(network, system, reference, p); + // Partition the network and perform the parallel function evaluation. + const RunResult r = evaluatePartitioning(network, system, reference, p); if (!r.success) { @@ -405,17 +443,11 @@ int main(int argc, char const* argv[]) return 1; } - // Output the results from partition evaluation - std::cout << std::format("{:<16d}{:>14.4f} s{:>16.4f} s{:>11.2f}x{:>14.3e} {:>16s}\n", - r.num_partitions, - r.partition_eval_time, - r.monolithic_eval_time, - r.speedup, - r.max_error, - r.jacobian_status); + // Output the partition evaluation results. + printResult(r); } delete system; return 0; -} +} \ No newline at end of file From 4e3806c7057826d60579ab3d1f4e9a0f713cdd74 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Tue, 8 Sep 2026 14:44:02 +0000 Subject: [PATCH 16/18] Updated Partition helper files --- .../ExamplesHelper/CMakeLists.txt | 17 +- .../ExamplesHelper/JacTestHelper.hpp | 446 +++++++++--------- .../ExamplesHelper/PartitionUtilities.hpp | 12 +- 3 files changed, 239 insertions(+), 236 deletions(-) diff --git a/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt b/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt index d80aa95e86..1996050d64 100644 --- a/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt +++ b/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt @@ -1,6 +1,6 @@ -add_library(power_elec_microgrid_network INTERFACE) - -add_library(GridKit::power_elec_microgrid_network ALIAS power_elec_microgrid_network) +add_library(power_elec_microgrid_network + INTERFACE MicrogridNetwork.hpp + SystemAssembler.hpp) target_link_libraries( power_elec_microgrid_network @@ -8,3 +8,14 @@ target_link_libraries( GridKit::power_elec_microbusdq GridKit::power_elec_microline GridKit::power_elec_microload) + +add_library(power_elec_partition_examples_helper + INTERFACE PartitionUtilities.hpp + JacTestHelper.hpp) + +target_link_libraries( + power_elec_partition_examples_helper + INTERFACE GridKit::power_elec_partition_interfaces + GridKit::testing + GridKit::sparse_matrix + power_elec_microgrid_network) diff --git a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp index e4411e2c2c..5b33cd18e2 100644 --- a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp +++ b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp @@ -1,3 +1,4 @@ +// JacTestHelper.hpp #pragma once #include @@ -9,263 +10,252 @@ #include #include #include -#include #include -namespace GridKit +/** + * @brief Verify that the subsystem Jacobian is an exact subset of the + * full-system Jacobian to ensure accuracy of subsystem Jacobians. + * + * Each subsystem Jacobian entry should match the corresponding entry in the + * full-system Jacobian after converting local subsystem indices back to global + * indices. The check verifies that: + * + * - every subsystem entry exists in the full-system Jacobian, + * - every full-system entry whose row and column belong to the subsystem + * exists in the subsystem Jacobian, and + * - matching entries agree within the specified tolerance. + * + * @note When the components in a subsystem are evaluated in a different order + * than in the monolithic reference system, a larger tolerance may be + * required to account for floating-point roundoff introduced by the + * different summation order. + * + * @param full_jac Full-system Jacobian. + * @param sub_jac Subsystem Jacobian. + * @param subsystem Subsystem associated with sub_jac. + * @param tolerance Maximum allowed difference between matching entries. + * + * @return true if the subsystem Jacobian is an exact subset of the + * full-system Jacobian, false otherwise. + */ + +template +bool verifySubsystemJacobian( + GridKit::LinearAlgebra::CsrMatrix& full_jac, + GridKit::LinearAlgebra::CsrMatrix& sub_jac, + GridKit::SubsystemModel& subsystem, + ScalarT tolerance = std::numeric_limits::epsilon()) { - namespace Testing + constexpr auto host = GridKit::memory::HOST; + // Full-system CSR data. + const IdxT* full_rows = full_jac.getRowData(host); + const IdxT* full_cols = full_jac.getColData(host); + const ScalarT* 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 ScalarT* 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-to-local indices. + const auto& global_to_local = subsystem.getInternalMap(); + + // The internal map should contain one entry for every subsystem row. + if (global_to_local.size() != num_sub_rows) { - /** - * @brief Verify that the subsystem Jacobian is an exact subset of the - * full-system Jacobian to ensure accuracy of subsystem Jacobians. - * - * Each subsystem Jacobian entry should match the corresponding entry in the - * full-system Jacobian after converting local subsystem indices back to global - * indices. The check verifies that: - * - * - every subsystem entry exists in the full-system Jacobian, - * - every full-system entry whose row and column belong to the subsystem - * exists in the subsystem Jacobian, and - * - matching entries agree within the specified tolerance. - * - * @note When the components in a subsystem are evaluated in a different order - * than in the monolithic reference system, a larger tolerance may be - * required to account for floating-point roundoff introduced by the - * different summation order. - * - * @param full_jac Full-system Jacobian. - * @param sub_jac Subsystem Jacobian. - * @param subsystem Subsystem associated with sub_jac. - * @param tolerance Maximum allowed difference between matching entries. - * - * @return true if the subsystem Jacobian is an exact subset of the - * full-system Jacobian, false otherwise. - */ + std::cout << "Internal map size mismatch: map has " + << global_to_local.size() + << " entries, but subsystem Jacobian has " + << num_sub_rows + << " rows\n"; + + return false; + } + + // The subsystem Jacobian is expected to be square. + 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 map from local subsystem indices to global indices. + // This is needed to compare subsystem rows and columns with the + // corresponding entries in the full-system Jacobian. + 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) + { + + // Make sure the global index is valid for the full-system Jacobian. + if (global_index >= full_jac.getNumRows()) + { + std::cout << "Invalid global index " + << global_index + << " in subsystem internal map\n"; + + return false; + } - template - bool verifySubsystemJacobian( - GridKit::LinearAlgebra::CsrMatrix& full_jac, - GridKit::LinearAlgebra::CsrMatrix& sub_jac, - GridKit::SubsystemModel& subsystem, - ScalarT tolerance = std::numeric_limits::epsilon()) + // Make sure the local index is valid for the subsystem Jacobian. + if (local_index >= num_sub_rows) { - 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 ScalarT* full_vals = full_jac.getValues(host); + std::cout << "Invalid local index " + << local_index + << " mapped from global index " + << global_index << '\n'; - // Subsystem CSR data. - const IdxT* sub_rows = sub_jac.getRowData(host); - const IdxT* sub_cols = sub_jac.getColData(host); - const ScalarT* sub_vals = sub_jac.getValues(host); + return false; + } + + // Each local index should correspond to only one global index. + if (local_index_found[local_index]) + { + std::cout << "Duplicate local index " + << local_index + << " in subsystem internal map\n"; - const IdxT num_sub_rows = sub_jac.getNumRows(); - const IdxT num_sub_cols = sub_jac.getNumColumns(); + return false; + } - // The subsystem internal map stores global-to-local indices. - const auto& global_to_local = subsystem.getInternalMap(); + local_to_global[local_index] = global_index; + local_index_found[local_index] = true; + } - // The internal map should contain one entry for every subsystem row. - if (global_to_local.size() != num_sub_rows) + // Make sure every subsystem local index has a 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 with 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 the current subsystem row using global column indices. + * + * The subsystem Jacobian uses local column indices, while the full-system + * Jacobian uses global column indices. Convert each local column to its + * corresponding global column so the two rows can be compared directly. + */ + 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]; + + // Fail if the subsystem column index is outside the valid local range. + if (local_column >= num_sub_cols) { - std::cout << "Internal map size mismatch: map has " - << global_to_local.size() - << " entries, but subsystem Jacobian has " - << num_sub_rows - << " rows\n"; + std::cout << "Invalid subsystem column index " + << local_column + << " in local row " + << local_row << '\n'; - return false; + matches = false; + continue; } - // The subsystem Jacobian is expected to be square. - if (num_sub_rows != num_sub_cols) + const IdxT global_column = local_to_global[local_column]; + + // Fail if the subsystem row contains more than one entry for the same column. + if (sub_row_entries.find(global_column) != sub_row_entries.end()) { - std::cout << "Subsystem Jacobian is not square: " - << num_sub_rows - << " rows and " - << num_sub_cols - << " columns\n"; + std::cout << "Duplicate subsystem entry at (" + << global_row << ", " + << global_column << ")\n"; - return false; + matches = false; + continue; } - // Build the reverse map from local subsystem indices to global indices. - // This is needed to compare subsystem rows and columns with the - // corresponding entries in the full-system Jacobian. - std::vector local_to_global(num_sub_rows); - std::vector local_index_found(num_sub_rows, false); + sub_row_entries[global_column] = sub_vals[sub_index]; + } - for (const auto& [global_index, local_index] : global_to_local) - { + // Compare full-system entries whose columns belong to the subsystem. + for (IdxT full_index = full_begin; full_index < full_end; ++full_index) + { + const IdxT global_column = full_cols[full_index]; - // Make sure the global index is valid for the full-system Jacobian. - if (global_index >= full_jac.getNumRows()) - { - std::cout << "Invalid global index " - << global_index - << " in subsystem internal map\n"; - - return false; - } - - // Make sure the local index is valid for the subsystem Jacobian. - if (local_index >= num_sub_rows) - { - std::cout << "Invalid local index " - << local_index - << " mapped from global index " - << global_index << '\n'; - - return false; - } - - // Each local index should correspond to only one global index. - 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; + // No need to check if column belongs outside the subsystem. + if (global_to_local.find(global_column) == global_to_local.end()) + { + continue; } - // Make sure every subsystem local index has a global index. - for (IdxT local_index = 0; local_index < num_sub_rows; ++local_index) + auto sub_entry = sub_row_entries.find(global_column); + + // Then it must exist in the subsystem Jacobian; fail otherwise. + if (sub_entry == sub_row_entries.end()) { - if (!local_index_found[local_index]) - { - std::cout << "No global index maps to local index " - << local_index << '\n'; + std::cout << "Entry exists only in full Jacobian at (" + << global_row << ", " + << global_column << ")\n"; - return false; - } + matches = false; + continue; } - bool matches = true; + const ScalarT full_value = full_vals[full_index]; + const ScalarT sub_value = sub_entry->second; + const ScalarT difference = std::abs(full_value - sub_value); - // Compare each subsystem row with the corresponding full-system row. - for (IdxT local_row = 0; local_row < num_sub_rows; ++local_row) + // Different component evaluation orders can change the order of floating-point + // summation and introduce small roundoff differences. Use an appropriate + // tolerance when comparing against the monolithic reference. + if (!GridKit::Testing::isEqual(sub_value, full_value, tolerance)) { - 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 the current subsystem row using global column indices. - * - * The subsystem Jacobian uses local column indices, while the full-system - * Jacobian uses global column indices. Convert each local column to its - * corresponding global column so the two rows can be compared directly. - */ - 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]; - - // Fail if the subsystem column index is outside the valid local range. - 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]; - - // Fail if the subsystem row contains more than one entry for the same column. - 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 full-system entries whose columns belong to the subsystem. - for (IdxT full_index = full_begin; full_index < full_end; ++full_index) - { - const IdxT global_column = full_cols[full_index]; - - // No need to check if column belongs outside the subsystem. - if (global_to_local.find(global_column) == global_to_local.end()) - { - continue; - } - - auto sub_entry = sub_row_entries.find(global_column); - - // Then it must exist in the subsystem Jacobian; fail otherwise. - 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 ScalarT full_value = full_vals[full_index]; - const ScalarT sub_value = sub_entry->second; - const ScalarT difference = std::abs(full_value - sub_value); - - // Different component evaluation orders can change the order of floating-point - // summation and introduce small roundoff differences. Use an appropriate - // tolerance when comparing against the monolithic reference. - // if tolerance is not supplied we simply use default machine precision provided - // by GridKit's Test::isEqual - // Then the values must agree, fail otherwise - if (!GridKit::Testing::isEqual(sub_value, full_value, tolerance)) - { - std::cout << "Jacobian value mismatch at (" - << global_row << ", " - << global_column << "): " - << "full = " << full_value - << ", subsystem = " << sub_value - << ", difference = " << difference << '\n'; - - matches = false; - } - - // Remove matched entry - sub_row_entries.erase(sub_entry); - } - - // Any entries remaining must be missing from the - // full-system Jacobian, so this subsystem row contains incorrect entries, fail! - 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; - } + std::cout << "Jacobian value mismatch at (" + << global_row << ", " + << global_column << "): " + << "full = " << full_value + << ", subsystem = " << sub_value + << ", difference = " << difference << '\n'; + + matches = false; } - return matches; + // Remove matched entry + sub_row_entries.erase(sub_entry); + } + + // Any entries remaining must be missing from the + // full-system Jacobian, so this subsystem row contains incorrect entries, fail! + 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; } + } - } // namespace Testing -} // namespace GridKit + return matches; +} \ No newline at end of file diff --git a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp index 0ccc66138f..9ef165cd6d 100644 --- a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp +++ b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp @@ -1,17 +1,16 @@ +// PartitionUtilities.hpp #pragma once #include -#include #include #include #include -#include #include #include -#include "MicrogridNetwork.hpp" +#include /** * @brief Partition a scaled microgrid network into subsystem models. @@ -72,8 +71,7 @@ void partitionNetwork( // Divide the IBRs as evenly as possible among the partitions. // Each partition receives q IBRs, and the first r partitions receive // one additional IBR to account for any remainder. The partition spans - // [partition_begin, partition_end), where partition_end is the first - // IBR index not included in this partition. + // [partition_begin, partition_end). const IdxT partition_size = q + ((j < r) ? 1 : 0); const IdxT partition_end = partition_begin + partition_size; @@ -168,6 +166,7 @@ void evaluatePartitionResiduals( auto* external_y = partition->getExternalDataY().getData(); auto* external_yp = partition->getExternalDataYP().getData(); + // Supply external variable values required by this partition from neighboring subsystems. for (size_t i = 0; i < partition->getExternSize(); ++i) { const auto global_index = partition->getExternalDataIndices()[i]; @@ -182,6 +181,7 @@ void evaluatePartitionResiduals( auto* partition_y = partition->y().getData(); auto* partition_yp = partition->yp().getData(); + // Supply external variable values required by this partition. for (size_t i = 0; i < partition->getInternalSize(); i++) { partition_y[i] = y[partition->getNodeConnection(i)]; @@ -191,10 +191,12 @@ void evaluatePartitionResiduals( partition->y().setDataUpdated(); partition->yp().setDataUpdated(); + // Evaluate this partition's residuals partition->evaluateResidual(); auto* residual = partition->getResidual().getData(); + // Gather the residuals from the partition into the full monolithic vector for (size_t i = 0; i < partition->getInternalSize(); i++) { f[partition->getNodeConnection(i)] = residual[i]; From eae6899bbfd2eea58886f7eec5b41e1439886019 Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Tue, 8 Sep 2026 15:20:32 +0000 Subject: [PATCH 17/18] Updated partitioned examples --- .../PowerElectronics/Partition/CMakeLists.txt | 19 +++++-------------- .../Partition/PartitionMicrogrid.cpp | 13 ++++++------- .../Partition/PartitionScaleMicrogrid.cpp | 16 ++++++---------- 3 files changed, 17 insertions(+), 31 deletions(-) diff --git a/examples/PowerElectronics/Partition/CMakeLists.txt b/examples/PowerElectronics/Partition/CMakeLists.txt index cd182dcd34..c2a70273ea 100644 --- a/examples/PowerElectronics/Partition/CMakeLists.txt +++ b/examples/PowerElectronics/Partition/CMakeLists.txt @@ -1,9 +1,5 @@ -find_package(OpenMP) - add_executable(PartitionMicrogrid PartitionMicrogrid.cpp) -target_include_directories( - PartitionMicrogrid - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..) +add_executable(PartitionScaleMicrogrid PartitionScaleMicrogrid.cpp) target_link_libraries( PartitionMicrogrid @@ -12,12 +8,7 @@ target_link_libraries( GridKit::power_elec_microload GridKit::solvers_dyn GridKit::power_elec_microbusdq - GridKit::power_elec_partition_interfaces) - -add_executable(PartitionScaleMicrogrid PartitionScaleMicrogrid.cpp) -target_include_directories( - PartitionScaleMicrogrid - PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..) + power_elec_partition_examples_helper) target_link_libraries( PartitionScaleMicrogrid @@ -26,9 +17,9 @@ target_link_libraries( GridKit::power_elec_microload GridKit::solvers_dyn GridKit::power_elec_microbusdq - GridKit::power_elec_partition_interfaces) + power_elec_partition_examples_helper) -if(OpenMP_CXX_FOUND) +if(GRIDKIT_ENABLE_OPENMP) target_link_libraries(PartitionScaleMicrogrid PRIVATE OpenMP::OpenMP_CXX) endif() @@ -36,4 +27,4 @@ add_test(NAME PartitionMicrogrid COMMAND PartitionMicrogrid) add_test(NAME PartitionScaleMicrogrid COMMAND PartitionScaleMicrogrid) install(TARGETS PartitionMicrogrid RUNTIME DESTINATION bin) -install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin) +install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin) \ No newline at end of file diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp index 8b28930e3d..a619ceff68 100644 --- a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -12,9 +12,8 @@ #include #include -#include "PowerElectronicsExamplesHelper/JacTestHelper.hpp" -#include "PowerElectronicsExamplesHelper/MicrogridNetwork.hpp" -#include "PowerElectronicsExamplesHelper/PartitionUtilities.hpp" +#include +#include using Component = GridKit::CircuitComponent; using Node = GridKit::PowerElectronics::NodeBase; @@ -178,10 +177,10 @@ int main() auto* partition_jacobian = partition->getCsrJacobian(); jacobians_match = jacobians_match - && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, - *partition_jacobian, - *partition, - 1e-13); + && verifySubsystemJacobian(*system_jacobian, + *partition_jacobian, + *partition, + 1e-13); } if (!jacobians_match) diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp index b61b99acaf..08604be54e 100644 --- a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -8,15 +8,11 @@ #include #include -#include #include -#include -#include #include -#include "PowerElectronicsExamplesHelper/JacTestHelper.hpp" -#include "PowerElectronicsExamplesHelper/MicrogridNetwork.hpp" -#include "PowerElectronicsExamplesHelper/PartitionUtilities.hpp" +#include +#include using index_type = size_t; using real_type = double; @@ -220,10 +216,10 @@ RunResult evaluatePartitioning( partition->evaluateJacobian(); jacobian_match = jacobian_match - && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, - *partition->getCsrJacobian(), - *partition, - Jac_tol); + && verifySubsystemJacobian(*system_jacobian, + *partition->getCsrJacobian(), + *partition, + Jac_tol); } // --------------------------------------------------------------------------- From 4770822667a0365e1ccfe0d1559ef333d44f63fb Mon Sep 17 00:00:00 2001 From: abdourahmanbarry Date: Tue, 8 Sep 2026 20:30:42 +0000 Subject: [PATCH 18/18] Apply pre-commit fixes --- .../ExamplesHelper/CMakeLists.txt | 16 ++++++++++------ .../ExamplesHelper/JacTestHelper.hpp | 2 +- .../ExamplesHelper/PartitionUtilities.hpp | 2 +- .../PowerElectronics/Partition/CMakeLists.txt | 2 +- .../Partition/PartitionMicrogrid.cpp | 2 +- .../Partition/PartitionScaleMicrogrid.cpp | 2 +- 6 files changed, 15 insertions(+), 11 deletions(-) diff --git a/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt b/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt index 1996050d64..bb228893f1 100644 --- a/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt +++ b/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt @@ -1,6 +1,8 @@ -add_library(power_elec_microgrid_network - INTERFACE MicrogridNetwork.hpp - SystemAssembler.hpp) +add_library( + power_elec_microgrid_network + INTERFACE + MicrogridNetwork.hpp + SystemAssembler.hpp) target_link_libraries( power_elec_microgrid_network @@ -9,9 +11,11 @@ target_link_libraries( GridKit::power_elec_microline GridKit::power_elec_microload) -add_library(power_elec_partition_examples_helper - INTERFACE PartitionUtilities.hpp - JacTestHelper.hpp) +add_library( + power_elec_partition_examples_helper + INTERFACE + PartitionUtilities.hpp + JacTestHelper.hpp) target_link_libraries( power_elec_partition_examples_helper diff --git a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp index 5b33cd18e2..38f6e486b7 100644 --- a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp +++ b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp @@ -258,4 +258,4 @@ bool verifySubsystemJacobian( } return matches; -} \ No newline at end of file +} diff --git a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp index 9ef165cd6d..6a7511afb3 100644 --- a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp +++ b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp @@ -253,4 +253,4 @@ void assembleSystemLeftToRight( sys_model.addNode(&network.buses[i]); } -} \ No newline at end of file +} diff --git a/examples/PowerElectronics/Partition/CMakeLists.txt b/examples/PowerElectronics/Partition/CMakeLists.txt index c2a70273ea..6d52b02493 100644 --- a/examples/PowerElectronics/Partition/CMakeLists.txt +++ b/examples/PowerElectronics/Partition/CMakeLists.txt @@ -27,4 +27,4 @@ add_test(NAME PartitionMicrogrid COMMAND PartitionMicrogrid) add_test(NAME PartitionScaleMicrogrid COMMAND PartitionScaleMicrogrid) install(TARGETS PartitionMicrogrid RUNTIME DESTINATION bin) -install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin) \ No newline at end of file +install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp index a619ceff68..90bb5270a3 100644 --- a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -290,4 +290,4 @@ std::vector getNodeConnections(const std::vector& nodes) } return connections; -} \ No newline at end of file +} diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp index 08604be54e..b16c1701af 100644 --- a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -446,4 +446,4 @@ int main(int argc, char const* argv[]) delete system; return 0; -} \ No newline at end of file +}