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/ExamplesHelper/CMakeLists.txt b/examples/PowerElectronics/ExamplesHelper/CMakeLists.txt index d80aa95e8..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) - -add_library(GridKit::power_elec_microgrid_network ALIAS power_elec_microgrid_network) +add_library( + power_elec_microgrid_network + INTERFACE + MicrogridNetwork.hpp + SystemAssembler.hpp) target_link_libraries( power_elec_microgrid_network @@ -8,3 +10,16 @@ target_link_libraries( GridKit::power_elec_microbusdq GridKit::power_elec_microline GridKit::power_elec_microload) + +add_library( + power_elec_partition_examples_helper + INTERFACE + PartitionUtilities.hpp + JacTestHelper.hpp) + +target_link_libraries( + power_elec_partition_examples_helper + INTERFACE GridKit::power_elec_partition_interfaces + GridKit::testing + GridKit::sparse_matrix + power_elec_microgrid_network) diff --git a/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp new file mode 100644 index 000000000..38f6e486b --- /dev/null +++ b/examples/PowerElectronics/ExamplesHelper/JacTestHelper.hpp @@ -0,0 +1,261 @@ +// JacTestHelper.hpp +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +/** + * @brief Verify that the subsystem Jacobian is an exact subset of the + * full-system Jacobian to ensure accuracy of subsystem Jacobians. + * + * Each subsystem Jacobian entry should match the corresponding entry in the + * full-system Jacobian after converting local subsystem indices back to global + * indices. The check verifies that: + * + * - every subsystem entry exists in the full-system Jacobian, + * - every full-system entry whose row and column belong to the subsystem + * exists in the subsystem Jacobian, and + * - matching entries agree within the specified tolerance. + * + * @note When the components in a subsystem are evaluated in a different order + * than in the monolithic reference system, a larger tolerance may be + * required to account for floating-point roundoff introduced by the + * different summation order. + * + * @param full_jac Full-system Jacobian. + * @param sub_jac Subsystem Jacobian. + * @param subsystem Subsystem associated with sub_jac. + * @param tolerance Maximum allowed difference between matching entries. + * + * @return true if the subsystem Jacobian is an exact subset of the + * full-system Jacobian, false otherwise. + */ + +template +bool verifySubsystemJacobian( + GridKit::LinearAlgebra::CsrMatrix& full_jac, + GridKit::LinearAlgebra::CsrMatrix& sub_jac, + GridKit::SubsystemModel& subsystem, + ScalarT tolerance = std::numeric_limits::epsilon()) +{ + constexpr auto host = GridKit::memory::HOST; + // Full-system CSR data. + const IdxT* full_rows = full_jac.getRowData(host); + const IdxT* full_cols = full_jac.getColData(host); + const ScalarT* full_vals = full_jac.getValues(host); + + // Subsystem CSR data. + const IdxT* sub_rows = sub_jac.getRowData(host); + const IdxT* sub_cols = sub_jac.getColData(host); + const ScalarT* sub_vals = sub_jac.getValues(host); + + const IdxT num_sub_rows = sub_jac.getNumRows(); + const IdxT num_sub_cols = sub_jac.getNumColumns(); + + // The subsystem internal map stores global-to-local indices. + const auto& global_to_local = subsystem.getInternalMap(); + + // The internal map should contain one entry for every subsystem row. + if (global_to_local.size() != num_sub_rows) + { + std::cout << "Internal map size mismatch: map has " + << global_to_local.size() + << " entries, but subsystem Jacobian has " + << num_sub_rows + << " rows\n"; + + return false; + } + + // The subsystem Jacobian is expected to be square. + if (num_sub_rows != num_sub_cols) + { + std::cout << "Subsystem Jacobian is not square: " + << num_sub_rows + << " rows and " + << num_sub_cols + << " columns\n"; + + return false; + } + + // Build the reverse map from local subsystem indices to global indices. + // This is needed to compare subsystem rows and columns with the + // corresponding entries in the full-system Jacobian. + std::vector local_to_global(num_sub_rows); + std::vector local_index_found(num_sub_rows, false); + + for (const auto& [global_index, local_index] : global_to_local) + { + + // Make sure the global index is valid for the full-system Jacobian. + if (global_index >= full_jac.getNumRows()) + { + std::cout << "Invalid global index " + << global_index + << " in subsystem internal map\n"; + + return false; + } + + // Make sure the local index is valid for the subsystem Jacobian. + if (local_index >= num_sub_rows) + { + std::cout << "Invalid local index " + << local_index + << " mapped from global index " + << global_index << '\n'; + + return false; + } + + // Each local index should correspond to only one global index. + if (local_index_found[local_index]) + { + std::cout << "Duplicate local index " + << local_index + << " in subsystem internal map\n"; + + return false; + } + + local_to_global[local_index] = global_index; + local_index_found[local_index] = true; + } + + // Make sure every subsystem local index has a global index. + for (IdxT local_index = 0; local_index < num_sub_rows; ++local_index) + { + if (!local_index_found[local_index]) + { + std::cout << "No global index maps to local index " + << local_index << '\n'; + + return false; + } + } + + bool matches = true; + + // Compare each subsystem row with the corresponding full-system row. + for (IdxT local_row = 0; local_row < num_sub_rows; ++local_row) + { + const IdxT global_row = local_to_global[local_row]; + + const IdxT sub_begin = sub_rows[local_row]; + const IdxT sub_end = sub_rows[local_row + 1]; + const IdxT full_begin = full_rows[global_row]; + const IdxT full_end = full_rows[global_row + 1]; + + /* + * Store the current subsystem row using global column indices. + * + * The subsystem Jacobian uses local column indices, while the full-system + * Jacobian uses global column indices. Convert each local column to its + * corresponding global column so the two rows can be compared directly. + */ + std::unordered_map sub_row_entries; + + for (IdxT sub_index = sub_begin; sub_index < sub_end; ++sub_index) + { + const IdxT local_column = sub_cols[sub_index]; + + // Fail if the subsystem column index is outside the valid local range. + if (local_column >= num_sub_cols) + { + std::cout << "Invalid subsystem column index " + << local_column + << " in local row " + << local_row << '\n'; + + matches = false; + continue; + } + + const IdxT global_column = local_to_global[local_column]; + + // Fail if the subsystem row contains more than one entry for the same column. + if (sub_row_entries.find(global_column) != sub_row_entries.end()) + { + std::cout << "Duplicate subsystem entry at (" + << global_row << ", " + << global_column << ")\n"; + + matches = false; + continue; + } + + sub_row_entries[global_column] = sub_vals[sub_index]; + } + + // Compare full-system entries whose columns belong to the subsystem. + for (IdxT full_index = full_begin; full_index < full_end; ++full_index) + { + const IdxT global_column = full_cols[full_index]; + + // No need to check if column belongs outside the subsystem. + if (global_to_local.find(global_column) == global_to_local.end()) + { + continue; + } + + auto sub_entry = sub_row_entries.find(global_column); + + // Then it must exist in the subsystem Jacobian; fail otherwise. + if (sub_entry == sub_row_entries.end()) + { + std::cout << "Entry exists only in full Jacobian at (" + << global_row << ", " + << global_column << ")\n"; + + matches = false; + continue; + } + + const ScalarT full_value = full_vals[full_index]; + const ScalarT sub_value = sub_entry->second; + const ScalarT difference = std::abs(full_value - sub_value); + + // Different component evaluation orders can change the order of floating-point + // summation and introduce small roundoff differences. Use an appropriate + // tolerance when comparing against the monolithic reference. + if (!GridKit::Testing::isEqual(sub_value, full_value, tolerance)) + { + std::cout << "Jacobian value mismatch at (" + << global_row << ", " + << global_column << "): " + << "full = " << full_value + << ", subsystem = " << sub_value + << ", difference = " << difference << '\n'; + + matches = false; + } + + // Remove matched entry + sub_row_entries.erase(sub_entry); + } + + // Any entries remaining must be missing from the + // full-system Jacobian, so this subsystem row contains incorrect entries, fail! + for (const auto& entry : sub_row_entries) + { + const IdxT global_column = entry.first; + + std::cout << "Entry exists only in subsystem Jacobian at (" + << global_row << ", " + << global_column << ")\n"; + + matches = false; + } + } + + return matches; +} diff --git a/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp new file mode 100644 index 000000000..6a7511afb --- /dev/null +++ b/examples/PowerElectronics/ExamplesHelper/PartitionUtilities.hpp @@ -0,0 +1,256 @@ +// PartitionUtilities.hpp +#pragma once + +#include + +#include +#include +#include + +#include +#include + +#include + +/** + * @brief Partition a scaled microgrid network into subsystem models. + * + * Divides the network into contiguous groups of IBRs and creates one + * subsystem for each group. The IBRs are distributed as evenly as possible, + * with any remainder assigned to the first partitions. + * + * The reference signal node is added to the first subsystem. A partition + * interface is added at the right boundary of every subsystem except the + * last to represent its connection to the neighboring subsystem. + * + * @param[in,out] network Scaled microgrid network to partition. + * @param[out] subsystems Vector populated with the created subsystem models. + * @param[in] num_partitions Number of subsystems to create. + * + * @pre The network has been constructed and contains @c 2*N_size IBRs. + * @pre @p num_partitions is greater than zero and does not exceed the + * number of IBRs in the network. + * @pre The generators, buses, virtual DQ buses, loads, and lines referenced + * by the network have valid lifetimes for use by the created subsystems. + * + * @post @p subsystems contains exactly @p num_partitions subsystem models. + * @post Every IBR in the network belongs to exactly one subsystem. + * @post Every subsystem except the last has a partition interface at its + * right boundary. + * + * @note The subsystem models are dynamically allocated. The caller is + * responsible for releasing and deleting them. + */ +template +void partitionNetwork( + ScaleMicrogridNetwork& network, + std::vector*>& subsystems, + size_t num_partitions) +{ + const size_t num_ibrs = 2 * network.N_size; + + assert(num_partitions <= num_ibrs); + + subsystems.resize(num_partitions); + + IdxT q = num_ibrs / num_partitions; + IdxT r = num_ibrs % num_partitions; + + IdxT partition_begin = 0; + + for (IdxT j = 0; j < num_partitions; j++) + { + auto* partition = new GridKit::SubsystemModel(); + + // Add the reference signal node to the first partition. + if (j == 0) + { + partition->addNode(&network.dg_signal); + } + + // Divide the IBRs as evenly as possible among the partitions. + // Each partition receives q IBRs, and the first r partitions receive + // one additional IBR to account for any remainder. The partition spans + // [partition_begin, partition_end). + const IdxT partition_size = q + ((j < r) ? 1 : 0); + const IdxT partition_end = partition_begin + partition_size; + + // Add all components belonging to this partition. + for (IdxT i = partition_begin; i < partition_end; ++i) + { + partition->addComponent(network.generators[i]); + partition->addComponent(network.busesDQ[i]); + + if (network.loads[i] != nullptr) + { + partition->addComponent(network.loads[i]); + } + + // Each line is owned by the partition containing its right bus. + // line[i - 1] connects bus[i - 1] to bus[i]. + if (i > 0) + { + partition->addComponent(network.lines[i - 1]); + } + + partition->addNode(&network.buses[i]); + } + + // If another partition exists to the right, the line crossing this + // partition's right boundary is owned by that next partition. Add an + // interface here to the left partition. + if (partition_end < num_ibrs) + { + + const IdxT boundary_line_index = partition_end - 1; + + auto* bus_interface = new GridKit::BusPartitionInterface( + &network.buses[boundary_line_index], + network.lines[boundary_line_index], + network.model_id_next++); + + bus_interface->allocate(); + partition->addInterface(bus_interface); + } + + subsystems[j] = partition; + partition_begin = partition_end; + } +} + +/** + * @brief Evaluate subsystem residuals and reconstruct the global residual. + * + * Distributes the global state and state-derivative vectors to each + * subsystem using its internal and external index mappings. Each subsystem + * residual is then evaluated in parallel, and its internal residual entries + * are gathered into the global residual vector. + * + * @param[in] subsystems Subsystem models to evaluate. + * @param[in] y Global state vector. + * @param[in] yp Global state-derivative vector. + * @param[out] f Global residual vector reconstructed from the subsystem + * residuals. + * @param[in] time Current simulation time. + * @param[in] alpha Jacobian scaling parameter associated with the current + * time-integration evaluation. + * + * @pre All subsystem models have been allocated. + * @pre @p y and @p yp contain all global entries + * @pre @p f is large enough to contain every global residual entry referenced + * by the subsystems. + * @pre The subsystems have disjoint internal index sets so that parallel + * writes to @p f do not overlap. + * + * @post Each subsystem residual has been evaluated at the supplied + * @p time and @p alpha. + * @post @p f contains the reconstructed residual for all internal variables + * represented by the subsystems. + */ +template +void evaluatePartitionResiduals( + const std::vector*>& subsystems, + const std::vector& y, + const std::vector& yp, + std::vector& f, + ScalarT time, + ScalarT alpha) +{ +#ifdef _OPENMP +#pragma omp parallel for schedule(guided) +#endif + for (auto* partition : subsystems) + { + partition->updateTime(time, alpha); + + auto* external_y = partition->getExternalDataY().getData(); + auto* external_yp = partition->getExternalDataYP().getData(); + + // Supply external variable values required by this partition from neighboring subsystems. + for (size_t i = 0; i < partition->getExternSize(); ++i) + { + const auto global_index = partition->getExternalDataIndices()[i]; + + external_y[i] = y[global_index]; + external_yp[i] = yp[global_index]; + } + + partition->getExternalDataY().setDataUpdated(); + partition->getExternalDataYP().setDataUpdated(); + + auto* partition_y = partition->y().getData(); + auto* partition_yp = partition->yp().getData(); + + // Supply external variable values required by this partition. + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + partition_y[i] = y[partition->getNodeConnection(i)]; + partition_yp[i] = yp[partition->getNodeConnection(i)]; + } + + partition->y().setDataUpdated(); + partition->yp().setDataUpdated(); + + // Evaluate this partition's residuals + partition->evaluateResidual(); + + auto* residual = partition->getResidual().getData(); + + // Gather the residuals from the partition into the full monolithic vector + for (size_t i = 0; i < partition->getInternalSize(); i++) + { + f[partition->getNodeConnection(i)] = residual[i]; + } + } +} + +/** + * @brief Assemble the scaled microgrid from left to right. + * + * Adds the network components to @p sys_model in physical left-to-right + * order. For each bus location, the generator, virtual DQ bus, load, line, + * and bus associated with that location are added before moving to the next + * location in the network. + * + * This ordering is useful when comparing the monolithic system with + * partition-based evaluations that traverse the network from left to right. + * + * @param[in] network Constructed scaled microgrid network. + * @param[in,out] sys_model Power electronics model to populate. + * + * @pre @c network.N_size is greater than zero. + * @pre @p network has been constructed by buildScaleMicrogridNetwork(). + * @pre All component and node pointers stored in @p network are valid. + * + * @post All nodes and components in @p network have been added to + * @p sys_model in left-to-right network order. + */ +template +void assembleSystemLeftToRight( + ScaleMicrogridNetwork& network, + GridKit::PowerElectronicsModel& sys_model) +{ + const size_t num_ibrs = 2 * network.N_size; + + assert(network.N_size > 0); + + sys_model.addNode(&network.dg_signal); + + for (IdxT i = 0; i < num_ibrs; ++i) + { + sys_model.addComponent(network.generators[i]); + sys_model.addComponent(network.busesDQ[i]); + + if (network.loads[i] != nullptr) + { + sys_model.addComponent(network.loads[i]); + } + + if (i > 0) + { + sys_model.addComponent(network.lines[i - 1]); + } + + sys_model.addNode(&network.buses[i]); + } +} diff --git a/examples/PowerElectronics/Partition/CMakeLists.txt b/examples/PowerElectronics/Partition/CMakeLists.txt new file mode 100644 index 000000000..6d52b0249 --- /dev/null +++ b/examples/PowerElectronics/Partition/CMakeLists.txt @@ -0,0 +1,30 @@ +add_executable(PartitionMicrogrid PartitionMicrogrid.cpp) +add_executable(PartitionScaleMicrogrid PartitionScaleMicrogrid.cpp) + +target_link_libraries( + PartitionMicrogrid + PRIVATE GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::solvers_dyn + GridKit::power_elec_microbusdq + power_elec_partition_examples_helper) + +target_link_libraries( + PartitionScaleMicrogrid + PRIVATE GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::solvers_dyn + GridKit::power_elec_microbusdq + power_elec_partition_examples_helper) + +if(GRIDKIT_ENABLE_OPENMP) + target_link_libraries(PartitionScaleMicrogrid PRIVATE OpenMP::OpenMP_CXX) +endif() + +add_test(NAME PartitionMicrogrid COMMAND PartitionMicrogrid) +add_test(NAME PartitionScaleMicrogrid COMMAND PartitionScaleMicrogrid) + +install(TARGETS PartitionMicrogrid RUNTIME DESTINATION bin) +install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp new file mode 100644 index 000000000..90bb5270a --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp @@ -0,0 +1,293 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +using Component = GridKit::CircuitComponent; +using Node = GridKit::PowerElectronics::NodeBase; +using Subsystem = GridKit::SubsystemModel; +using System = GridKit::PowerElectronicsModel; + +std::vector getComponentConnections(const std::vector& components); +std::vector getNodeConnections(const std::vector& nodes); + +/* + * Verify partitioned residual and Jacobian evaluation against the + * monolithic microgrid model. + * + * The microgrid is manually divided into two subsystems. A partition interface + * is introduced at the boundary between bus 1 and the line connecting buses 1 + * and 2. Each subsystem is evaluated independently, and its residuals are gathered + * and compared with the monolithic reference. + */ +int main() +{ + + constexpr size_t N_size = 2; + constexpr double time = 0.1; + constexpr double alpha = 0.1; + + bool use_jac = true; + + // --------------------------------------------------------------------------- + // build the grid network and assemble the system model + // --------------------------------------------------------------------------- + ScaleMicrogridNetwork network(N_size); + + auto* system = new System(use_jac); + + assembleSystemLeftToRight(network, *system); + system->allocate(); + + std::vector y(system->size()); + std::vector yp(system->size()); + + for (size_t i = 0; i < system->size(); ++i) + { + y[i] = static_cast(i + 1); + yp[i] = static_cast(i + 1); + } + + auto* system_y = system->y().getData(); + auto* system_yp = system->yp().getData(); + + for (size_t i = 0; i < system->size(); ++i) + { + system_y[i] = y[i]; + system_yp[i] = yp[i]; + } + + system->y().setDataUpdated(); + system->yp().setDataUpdated(); + + system->updateTime(time, alpha); + system->evaluateResidual(); + system->evaluateJacobian(); + + auto* system_jacobian = system->getCsrJacobian(); + auto* system_residual = system->getResidual().getData(); + + //------------------------------------------------------------------------------ + // Gather all global indices belonging partition 1 and 2 to test release() later + //------------------------------------------------------------------------------ + std::vector components = { + network.generators[0], + network.generators[1], + network.lines[0], + network.loads[0], + network.busesDQ[0], + network.busesDQ[1], + network.generators[2], + network.generators[3], + network.lines[1], + network.lines[2], + network.loads[2], + network.busesDQ[2], + network.busesDQ[3]}; + + std::vector nodes = { + &network.dg_signal, + &network.buses[0], + &network.buses[1], + &network.buses[2], + &network.buses[3]}; + + const auto original_component_connections = getComponentConnections(components); + const auto original_node_connections = getNodeConnections(nodes); + + // -------------------------------------------------------------------------------- + // Create 2 Partitions and a Partition interfaces to split to network + // -------------------------------------------------------------------------------- + + auto* bus_interface = new GridKit::BusPartitionInterface( + &network.buses[1], + network.lines[1], + 14); + + if (int err = bus_interface->allocate()) + { + return err; + } + + auto* partition1 = new Subsystem(); + auto* partition2 = new Subsystem(); + + // -------------------------------------------------------------------------------- + // Manually add components, nodes and a bus partition interface to Partition 1 + // -------------------------------------------------------------------------------- + + partition1->addNode(&network.dg_signal); + partition1->addComponent(network.generators[0]); + partition1->addComponent(network.busesDQ[0]); + partition1->addComponent(network.loads[0]); + partition1->addNode(&network.buses[0]); + partition1->addComponent(network.lines[0]); + partition1->addComponent(network.generators[1]); + partition1->addComponent(network.busesDQ[1]); + partition1->addInterface(bus_interface); + partition1->addNode(&network.buses[1]); + + // --------------------------------------------------------------------------- + // Manually add components and nodes to Partition 2 + // --------------------------------------------------------------------------- + + partition2->addComponent(network.generators[2]); + partition2->addComponent(network.busesDQ[2]); + partition2->addComponent(network.loads[2]); + partition2->addComponent(network.lines[1]); + partition2->addNode(&network.buses[2]); + partition2->addComponent(network.generators[3]); + partition2->addComponent(network.busesDQ[3]); + partition2->addComponent(network.lines[2]); + partition2->addNode(&network.buses[3]); + + std::vector partitions = {partition1, partition2}; + + for (auto* partition : partitions) + { + partition->allocate(); + } + + std::vector partition_residual(system->size(), 0.0); + + evaluatePartitionResiduals(partitions, y, yp, partition_residual, time, alpha); + + // --------------------------------------------------------------------------- + // Verify the subsystem Jacobians + // --------------------------------------------------------------------------- + + bool jacobians_match = true; + + for (auto* partition : partitions) + { + partition->evaluateJacobian(); + + auto* partition_jacobian = partition->getCsrJacobian(); + + jacobians_match = jacobians_match + && verifySubsystemJacobian(*system_jacobian, + *partition_jacobian, + *partition, + 1e-13); + } + + if (!jacobians_match) + { + std::cout << "ERROR: At least one subsystem Jacobian is incorrect!\n"; + return 1; + } + + // --------------------------------------------------------------------------- + // Gather and verify the subsystem residuals against monolithic residuals + // --------------------------------------------------------------------------- + + double max_error = 0.0; + + for (size_t i = 0; i < system->size(); ++i) + { + const double error = std::abs(system_residual[i] - partition_residual[i]) / std::abs(system_residual[i] + 1); + + max_error = std::max(max_error, error); + } + + const double machine_epsilon = + std::numeric_limits::epsilon(); + + const bool residuals_match = max_error <= machine_epsilon; + + std::cout << "\nPartition Microgrid Validation\n"; + std::cout << "------------------------------\n"; + + std::cout << std::left + << std::setw(32) << "Maximum residual error:" + << std::setprecision(16) + << max_error << '\n'; + + std::cout << std::left + << std::setw(32) << "Machine epsilon:" + << machine_epsilon << '\n'; + + std::cout << std::left + << std::setw(32) << "Residuals matched:" + << (residuals_match ? "True" : "False") + << '\n'; + + // --------------------------------------------------------------------------- + // Verify subsystem release() from SubsystemModel + // --------------------------------------------------------------------------- + + for (auto* partition : partitions) + { + partition->release(); + } + + const bool components_restored = getComponentConnections(components) == original_component_connections; + const bool nodes_restored = getNodeConnections(nodes) == original_node_connections; + + if (!components_restored || !nodes_restored) + { + std::cout << "ERROR: Subsystem release did not restore " + "the original global connection indices!\n"; + + return 1; + } + + // --------------------------------------------------------------------------- + // Clean up + // --------------------------------------------------------------------------- + delete system; + + for (auto* partition : partitions) + { + delete partition; + } + + return residuals_match ? 0 : 1; +} + +/** + * @brief Collect the connection indices of a set of components. + */ +std::vector getComponentConnections(const std::vector& components) +{ + std::vector connections; + + for (const auto* component : components) + { + for (size_t i = 0; i < component->size(); ++i) + { + connections.push_back(component->getNodeConnection(i)); + } + } + + return connections; +} + +/** + * @brief Collect the connection indices of a set of nodes. + */ +std::vector getNodeConnections(const std::vector& nodes) +{ + std::vector connections; + + for (auto* node : nodes) + + for (size_t i = 0; i < node->size(); ++i) + { + connections.push_back(node->getNodeConnection(i).idx_); + } + + return connections; +} diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp new file mode 100644 index 000000000..b16c1701a --- /dev/null +++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp @@ -0,0 +1,449 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +using index_type = size_t; +using real_type = double; + +using SignalNode = GridKit::PowerElectronics::SignalNode; +using Bus = GridKit::PowerElectronics::MicrogridBus; +using BusDQ = GridKit::MicrogridBusDQ; +using DGGenerator = GridKit::DistributedGenerator; +using Line = GridKit::MicrogridLine; +using Load = GridKit::MicrogridLoad; + +/* + * Output data for each parallel function evaluation run + */ +struct RunResult +{ + bool success = true; + index_type num_partitions; + real_type partition_eval_time; // seconds + real_type monolithic_eval_time; // seconds + real_type speedup; // monolithic / partition + real_type max_error; + std::string jacobian_status; +}; + +/* + * Reference data from the monolithic system evaluation. + * + * The monolithic system is evaluated once and its states, residual, and + * evaluation time are stored here. Each partition configuration reuses this + * data for validation and performance comparison. + */ +struct MonolithicReference +{ + std::vector y; + std::vector yp; + std::vector residual; + real_type monolithic_eval_time; +}; + +/** + * @brief Evaluate the monolithic system and store reference data. + * + * Initializes the global state and state-derivative vectors with deterministic + * values, evaluates and times the monolithic residual, evaluates the full + * Jacobian, and stores the resulting data for later comparison with + * partitioned evaluations. + * + * The residual evaluation time is measured independently from the Jacobian + * evaluation so that the stored timing represents only the cost of the + * monolithic residual evaluation. + * + * @param[in,out] system Allocated monolithic power electronics system. + * + * @return MonolithicReference containing the state vector, state-derivative + * vector, residual, and residual evaluation time. + * + * @pre @p system is non-null. + * @pre @p system has already been assembled and allocated. + * + * @post The state and state-derivative vectors of @p system contain the + * deterministic reference values. + * @post The monolithic residual and Jacobian have been evaluated. + * @post The returned reference contains an independent copy of the + * monolithic residual. + */ +MonolithicReference evaluateMonolithicSystem(GridKit::PowerElectronicsModel* system) +{ + MonolithicReference reference; + + // --------------------------------------------------------------------------- + // Initialize the reference state and state-derivative vectors + // --------------------------------------------------------------------------- + reference.y.resize(system->size()); + reference.yp.resize(system->size()); + + for (index_type i = 0; i < system->size(); ++i) + { + reference.y[i] = static_cast(i + 1); + reference.yp[i] = static_cast(i + 1); + } + + // Copy the reference values into the vectors owned by the system model. + auto* system_y = system->y().getData(); + auto* system_yp = system->yp().getData(); + + for (index_type i = 0; i < system->size(); ++i) + { + system_y[i] = reference.y[i]; + system_yp[i] = reference.yp[i]; + } + + system->y().setDataUpdated(); + system->yp().setDataUpdated(); + + // --------------------------------------------------------------------------- + // Evaluate and time the monolithic residual + // --------------------------------------------------------------------------- + + auto start_time = std::chrono::high_resolution_clock::now(); + + system->updateTime(0.1, 0.1); + system->evaluateResidual(); + + auto end_time = std::chrono::high_resolution_clock::now(); + + reference.monolithic_eval_time = std::chrono::duration(end_time - start_time).count(); + + // evaluate the monolithic Jacobian + system->evaluateJacobian(); + + // Store a copy of the residual return reference data + auto* residual = system->getResidual().getData(); + + reference.residual.assign(residual, residual + system->size()); + + return reference; +} + +/** + * @brief Evaluate one partition configuration and validate it against the + * monolithic reference system. + * + * Creates the requested subsystem decomposition, allocates each subsystem, + * evaluates the partitioned residual, verifies each subsystem Jacobian against + * the monolithic Jacobian, and compares the reconstructed global residual with + * the stored monolithic residual. + * + * The partitioned residual evaluation is timed independently and compared with + * the previously measured monolithic residual evaluation time to compute the + * resulting speedup. + * + * @param[in,out] network Scaled microgrid network to partition. + * @param[in] system Allocated monolithic system used as the reference. + * @param[in] reference Stored monolithic state, derivative, residual, and + * timing information. + * @param[in] num_partitions Number of subsystem partitions to create. + * + * @return Performance and validation results for the requested partition count. + * + * @pre @p system is non-null and has already been evaluated. + * @pre @p network has been constructed using buildScaleMicrogridNetwork(). + * @pre @p num_partitions is greater than zero and does not exceed the number + * of IBRs in the network. + * @pre @p reference contains state and derivative vectors consistent with the + * size of @p system. + * + * @post All temporary subsystem models created by this function are released + * and deleted before returning. + * @post The returned result reports the partition timing, speedup, residual + * error, and subsystem Jacobian validation status. + */ +template +RunResult evaluatePartitioning( + ScaleMicrogridNetwork& network, + GridKit::PowerElectronicsModel* system, + const MonolithicReference& reference, + index_type num_partitions) +{ + + // --------------------------------------------------------------------------- + // Create and allocate subsystem partitions + // --------------------------------------------------------------------------- + + std::vector*> subsystems; + + partitionNetwork(network, subsystems, num_partitions); + + for (auto* partition : subsystems) + { + partition->allocate(); + } + + // Global residual reconstructed from the subsystem residuals. + std::vector f(system->size(), 1.0); + // Elementwise error between the monolithic and reconstructed residuals. + std::vector error(system->size(), 1.0); + + // --------------------------------------------------------------------------- + // Evaluate and time the partitioned residual + // --------------------------------------------------------------------------- + auto start_time = std::chrono::high_resolution_clock::now(); + + evaluatePartitionResiduals(subsystems, reference.y, reference.yp, f, 0.1, 0.1); + + auto end_time = std::chrono::high_resolution_clock::now(); + + auto partition_eval_time = std::chrono::duration(end_time - start_time); + + // --------------------------------------------------------------------------- + // Verify subsystem Jacobians against the monolithic Jacobian + // --------------------------------------------------------------------------- + + auto* system_jacobian = system->getCsrJacobian(); + + bool jacobian_match = true; + real_type Jac_tol = 1e-13; + + for (auto* partition : subsystems) + { + partition->evaluateJacobian(); + + jacobian_match = jacobian_match + && verifySubsystemJacobian(*system_jacobian, + *partition->getCsrJacobian(), + *partition, + Jac_tol); + } + + // --------------------------------------------------------------------------- + // Compare the reconstructed and monolithic residuals + // --------------------------------------------------------------------------- + real_type max_error = 0.0; + + for (index_type i = 0; i < system->size(); ++i) + { + error[i] = std::abs(f[i] - reference.residual[i]) / (std::abs(reference.residual[i]) + 1.0); + + if (max_error < error[i]) + { + max_error = error[i]; + } + } + + // --------------------------------------------------------------------------- + // Store performance and validation results + // --------------------------------------------------------------------------- + RunResult result; + + result.num_partitions = num_partitions; + result.partition_eval_time = partition_eval_time.count(); + result.monolithic_eval_time = reference.monolithic_eval_time; + result.speedup = reference.monolithic_eval_time / partition_eval_time.count(); + result.max_error = max_error; + result.jacobian_status = jacobian_match ? "Correct" : "Wrong"; + + // --------------------------------------------------------------------------- + // Check validation results + // --------------------------------------------------------------------------- + if (!jacobian_match) + { + std::cout << "ERROR: At least one subsystem Jacobian is incorrect!" + << std::endl; + + result.success = false; + } + + if (max_error > std::numeric_limits::epsilon()) + { + std::cout << "ERROR: Max Error too high!: " << max_error << std::endl; + result.success = false; + } + + // --------------------------------------------------------------------------- + // Release and destroy temporary subsystem models + // --------------------------------------------------------------------------- + + for (auto* partition : subsystems) + { + partition->release(); + delete partition; + } + + return result; +} + +/** + * @brief Print the header for the partition evaluation results table. + */ +void printResultsHeader() +{ + std::cout << std::left + << std::setw(16) << "num_partitions" + << std::right + << std::setw(16) << "partition_time" + << std::setw(18) << "monolithic_time" + << std::setw(12) << "speedup" + << std::setw(14) << "error" + << std::setw(16) << "Jacobians" + << '\n'; + + std::cout << std::string(92, '-') << '\n'; +} + +/** + * @brief Print one row of partition evaluation results. + * + * @param[in] result Results from the partitioned system evaluation. + */ +void printResult(const RunResult& result) +{ + std::cout << std::left + << std::setw(16) << result.num_partitions + << std::right + << std::setw(14) << std::fixed << std::setprecision(4) + << result.partition_eval_time << " s" + << std::setw(16) << result.monolithic_eval_time << " s" + << std::setw(11) << std::setprecision(2) + << result.speedup << "x" + << std::setw(14) << std::scientific << std::setprecision(3) + << result.max_error << ' ' + << std::setw(16) << result.jacobian_status + << '\n'; + + std::cout << std::defaultfloat; +} + +/** + * @brief Benchmark and validate partitioned residual evaluation of a large + * scaled microgrid. It also confirms the correctness of the subsystem + * Jacobian. + * + * Builds a scaled microgrid once, constructs a monolithic reference system, + * and compares several subsystem decompositions of the same physical network. + * + * For each requested partition count, the benchmark measures the partitioned + * residual evaluation time, computes its speedup relative to the monolithic + * evaluation, verifies the reconstructed residual, and checks every subsystem + * Jacobian against the monolithic Jacobian. + * + * The physical network and monolithic system are constructed only once. + * Individual subsystem models are created and destroyed separately for each + * partition count. + * + * @return 0 if all residual and Jacobian validation checks pass; otherwise 1. + */ +int main(int argc, char const* argv[]) +{ + index_type N_size = 5000; + + std::vector num_partitions_list = {500}; + + /* + * If command-line arguments are provided, the first argument specifies + * N_size and all remaining arguments specify partition counts. + * + * Example: + * ./PartitionScaleMicrogrid 5000 10 48 100 500 + */ + if (argc > 1) + { + try + { + N_size = static_cast(std::stoull(argv[1])); + + if (N_size < 1) + { + std::cerr << "ERROR: N_size must be at least 1.\n"; + return 1; + } + + // When N_size is supplied explicitly, at least one partition count + // must also be supplied. + if (argc < 3) + { + std::cerr << "ERROR: At least one partition count must be provided " + << "when N_size is specified.\n"; + return 1; + } + + num_partitions_list.clear(); + + for (int i = 2; i < argc; ++i) + { + index_type num_partitions = static_cast(std::stoull(argv[i])); + + if (num_partitions < 1) + { + std::cerr << "ERROR: Number of partitions must be at least 1.\n"; + return 1; + } + + if (num_partitions > 2 * N_size) + { + std::cerr << "ERROR: Number of partitions (" + << num_partitions + << ") cannot exceed the number of IBRs (" + << 2 * N_size + << ").\n"; + return 1; + } + + num_partitions_list.push_back(num_partitions); + } + } + catch (const std::exception& e) + { + std::cerr + << "ERROR: Invalid command-line argument: " + << e.what() + << "\n"; + + return 1; + } + } + + const bool use_jac = true; + + // Build the physical network once. + ScaleMicrogridNetwork network(N_size); + + // Build, assemble and allocate the monolithic system once. + auto* system = new GridKit::PowerElectronicsModel(use_jac); + + assembleSystemLeftToRight(network, *system); + system->allocate(); + + // Evaluate the monolithic reference once. + const MonolithicReference reference = evaluateMonolithicSystem(system); + + printResultsHeader(); + + // Only the partitioned system is rebuilt and evaluated for each partition count. + for (const index_type p : num_partitions_list) + { + assert(p <= 2 * N_size); + + // Partition the network and perform the parallel function evaluation. + const RunResult r = evaluatePartitioning(network, system, reference, p); + + if (!r.success) + { + delete system; + return 1; + } + + // Output the partition evaluation results. + printResult(r); + } + + delete system; + + return 0; +}