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