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/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..1f1109720d --- /dev/null +++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp @@ -0,0 +1,168 @@ +/** + * @file BusPartitionInterface.hpp + * @author Abdourahman Barry (abdourahman@vt.edu) + * + */ + +#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..9f013ad673 --- /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 diff --git a/GridKit/Model/PowerElectronics/SubsystemModel.hpp b/GridKit/Model/PowerElectronics/SubsystemModel.hpp new file mode 100644 index 0000000000..9ff9e41fd3 --- /dev/null +++ b/GridKit/Model/PowerElectronics/SubsystemModel.hpp @@ -0,0 +1,856 @@ +/** + * @file BusPartitionInterface.hpp + * @author Abdourahman Barry (abdourahman@vt.edu) + * + */ +#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. + * + * @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. + */ + template + class SubsystemModel : public PowerElectronicsModel + { + public: + struct ForcingData + { + std::vector y; + std::vector yp; + }; + + using TimeFunction = std::function; + + protected: + using SystemModel = PowerElectronicsModel; + using VectorT = typename SystemModel::VectorT; + 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); + } + } + } + + { + 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; + 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_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); + } + } + } + + // 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_.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_.getData()); + std::copy(forcing.yp.begin(), forcing.yp.end(), yp_ext_data_.getData()); + } + + 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_; + } + + VectorT& getExternalDataY() + { + return y_ext_data_; + } + + VectorT& getExternalDataYP() + { + return yp_ext_data_; + } + + VectorT& 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 data. + */ + VectorT y_ext_data_; + + /** + * @brief subsystem external state derivative + */ + VectorT yp_ext_data_; + + /** + * @brief subsystem external residual data + */ + VectorT 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 diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 51c8dcee18..9796ee9f16 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 { @@ -22,6 +26,7 @@ namespace GridKit using component_type = CircuitComponent; using node_type = PowerElectronics::NodeBase; + protected: using CircuitComponent::abs_tol_; using CircuitComponent::allocated_; using CircuitComponent::allocateVectors; @@ -40,17 +45,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 +110,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 +275,7 @@ namespace GridKit * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ - int initialize() final + int initialize() override { // Initialize components for (const auto& component : components_) @@ -304,7 +298,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() override { // Ask all component to tag their differentiables for (size_t i = 0; i < components_.size(); i++) @@ -367,7 +361,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 +402,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,7 +424,12 @@ 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; @@ -514,7 +513,7 @@ namespace GridKit allocated_ = false; } - private: + protected: static constexpr IdxT neg1_ = INVALID_INDEX; std::vector components_; 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/ExamplesHelper/CMakeLists.txt b/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt index d80aa95e86..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) - -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 +10,16 @@ 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 new file mode 100644 index 0000000000..38f6e486b7 --- /dev/null +++ b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp @@ -0,0 +1,261 @@ +// JacTestHelper.hpp +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +/** + * @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()) +{ + 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) + { + 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 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 (!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; + } + } + + return matches; +} diff --git a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp new file mode 100644 index 0000000000..6a7511afb3 --- /dev/null +++ b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp @@ -0,0 +1,256 @@ +// PartitionUtilities.hpp +#pragma once + +#include + +#include +#include +#include + +#include +#include + +#include + +/** + * @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 partition_begin = 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); + } + + // 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). + const IdxT partition_size = q + ((j < r) ? 1 : 0); + const IdxT partition_end = partition_begin + partition_size; + + // 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]); + + if (network.loads[i] != nullptr) + { + partition->addComponent(network.loads[i]); + } + + // 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]); + } + + partition->addNode(&network.buses[i]); + } + + // 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) + { + + 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); + } + + 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) + { + partition->updateTime(time, alpha); + + 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]; + + external_y[i] = y[global_index]; + external_yp[i] = yp[global_index]; + } + + partition->getExternalDataY().setDataUpdated(); + partition->getExternalDataYP().setDataUpdated(); + + 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)]; + partition_yp[i] = yp[partition->getNodeConnection(i)]; + } + + 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]; + } + } +} + +/** + * @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 (i > 0) + { + sys_model.addComponent(network.lines[i - 1]); + } + + sys_model.addNode(&network.buses[i]); + } +} diff --git a/examples/PowerElectronics/Partition/CMakeLists.txt b/examples/PowerElectronics/Partition/CMakeLists.txt new file mode 100644 index 0000000000..6d52b02493 --- /dev/null +++ b/examples/PowerElectronics/Partition/CMakeLists.txt @@ -0,0 +1,30 @@ +add_executable(PartitionMicrogrid PartitionMicrogrid.cpp) +add_executable(PartitionScaleMicrogrid PartitionScaleMicrogrid.cpp) + +target_link_libraries( + PartitionMicrogrid + PRIVATE GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::solvers_dyn + GridKit::power_elec_microbusdq + power_elec_partition_examples_helper) + +target_link_libraries( + PartitionScaleMicrogrid + PRIVATE GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::solvers_dyn + GridKit::power_elec_microbusdq + power_elec_partition_examples_helper) + +if(GRIDKIT_ENABLE_OPENMP) + 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..90bb5270a3 --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -0,0 +1,293 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +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 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 + // --------------------------------------------------------------------------- + ScaleMicrogridNetwork network(N_size); + + auto* system = new System(use_jac); + + 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[0], + network.loads[0], + network.busesDQ[0], + network.busesDQ[1], + network.generators[2], + network.generators[3], + network.lines[1], + network.lines[2], + 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[1], + 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[0]); + 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[1]); + partition2->addNode(&network.buses[2]); + partition2->addComponent(network.generators[3]); + partition2->addComponent(network.busesDQ[3]); + partition2->addComponent(network.lines[2]); + partition2->addNode(&network.buses[3]); + + std::vector partitions = {partition1, partition2}; + + for (auto* partition : partitions) + { + partition->allocate(); + } + + std::vector partition_residual(system->size(), 0.0); + + 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 + && verifySubsystemJacobian(*system_jacobian, + *partition_jacobian, + *partition, + 1e-13); + } + + 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; +} diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp new file mode 100644 index 0000000000..b16c1701af --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -0,0 +1,449 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +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 (index_type 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 (index_type 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. + */ +template +RunResult evaluatePartitioning( + ScaleMicrogridNetwork& network, + GridKit::PowerElectronicsModel* system, + const MonolithicReference& reference, + index_type num_partitions) +{ + + // --------------------------------------------------------------------------- + // Create and allocate subsystem partitions + // --------------------------------------------------------------------------- + + std::vector*> subsystems; + + 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(); + + 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; + real_type Jac_tol = 1e-13; + + for (auto* partition : subsystems) + { + partition->evaluateJacobian(); + + jacobian_match = jacobian_match + && verifySubsystemJacobian(*system_jacobian, + *partition->getCsrJacobian(), + *partition, + Jac_tol); + } + + // --------------------------------------------------------------------------- + // Compare the reconstructed and monolithic residuals + // --------------------------------------------------------------------------- + real_type max_error = 0.0; + + for (index_type i = 0; i < system->size(); ++i) + { + error[i] = std::abs(f[i] - reference.residual[i]) / (std::abs(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 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 + * 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; + } + } + + const bool use_jac = true; + + // Build the physical network once. + ScaleMicrogridNetwork network(N_size); + + // Build, assemble and allocate the monolithic system once. + auto* system = new GridKit::PowerElectronicsModel(use_jac); + + assembleSystemLeftToRight(network, *system); + system->allocate(); + + // Evaluate the monolithic reference once. + const MonolithicReference reference = evaluateMonolithicSystem(system); + + printResultsHeader(); + + // Only the partitioned system is rebuilt and evaluated for each partition count. + for (const index_type p : num_partitions_list) + { + assert(p <= 2 * N_size); + + // Partition the network and perform the parallel function evaluation. + const RunResult r = evaluatePartitioning(network, system, reference, p); + + if (!r.success) + { + delete system; + return 1; + } + + // Output the partition evaluation results. + printResult(r); + } + + delete system; + + return 0; +} diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt index 94cd4fd912..756278102a 100644 --- a/tests/UnitTests/PowerElectronics/CMakeLists.txt +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -3,6 +3,11 @@ 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 +18,10 @@ target_link_libraries( GridKit::testing) add_test(NAME PowerElectronicsNodeTest 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..e5b4368b3b --- /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) ++\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) ++\left(0.7y_5\right) ++\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..71a1aaaefd --- /dev/null +++ b/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp @@ -0,0 +1,774 @@ +#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 + // --------------------------------------------------------------------- + + bus_interface_ = new BusPartitionInterface(&bus_, comp3_, 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(); +}