From 90f8ae5f126754d5d4532a9311cab0dbdf76a7a5 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 6 Nov 2025 11:57:11 -0500 Subject: [PATCH 01/29] Remove short-circuit from COO::axpy Previously, axpy would do no work if a == 0, but this can change sparsity patterns. Typically, a > 0 during simulations anyway, so this only affects when axpy's are done e.g. before simulation, when memory is unset. --- GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp index 80f80ecb..1cd2f6e9 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp @@ -292,11 +292,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, COO_Matrix& a) { - if (alpha == 0) - { - return; - } - if (!this->sorted_) { this->sortSparse(); @@ -372,9 +367,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, std::vector r, std::vector c, std::vector v) { - if (alpha == 0) - return; - if (!this->sorted_) { this->sortSparse(); From b4c18d47c786c4026a2c024df5817968688edf80 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 6 Nov 2025 13:24:11 -0500 Subject: [PATCH 02/29] Add sparsity pattern computation Also added a unit test which compares sparsity pattern generated by DependencyTracking::Variable and new sparsity pattern on microgrid problem. --- .../SystemModelPowerElectronics.hpp | 189 ++++++++- tests/UnitTests/CMakeLists.txt | 1 + .../UnitTests/PowerElectronics/CMakeLists.txt | 11 + .../PowerElectronics/SystemTests.hpp | 373 ++++++++++++++++++ .../PowerElectronics/runSystemTests.cpp | 14 + 5 files changed, 574 insertions(+), 14 deletions(-) create mode 100644 tests/UnitTests/PowerElectronics/CMakeLists.txt create mode 100644 tests/UnitTests/PowerElectronics/SystemTests.hpp create mode 100644 tests/UnitTests/PowerElectronics/runSystemTests.cpp diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 736f9d00..00e1aac8 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -81,11 +81,13 @@ namespace GridKit PowerElectronicsModel() { // Set system model parameters as default - rel_tol_ = 1e-4; - abs_tol_ = 1e-4; - this->max_steps_ = 2000; + rel_tol_ = 1e-4; + abs_tol_ = 1e-4; + this->max_steps_ = 2000; // By default don't use the jacobian - use_jac_ = false; + use_jac_ = false; + jac_sparsity_row_indices = nullptr; + jac_sparsity_col_indices = nullptr; } /** @@ -104,11 +106,13 @@ namespace GridKit IdxT max_steps = 2000) { // Set system model tolerances from input - rel_tol_ = rel_tol; - abs_tol_ = abs_tol; - this->max_steps_ = max_steps; + rel_tol_ = rel_tol; + abs_tol_ = abs_tol; + this->max_steps_ = max_steps; // Can choose if to use jacobian - use_jac_ = use_jac; + use_jac_ = use_jac; + jac_sparsity_row_indices = nullptr; + jac_sparsity_col_indices = nullptr; } /** @@ -123,18 +127,140 @@ namespace GridKit { for (auto comp : this->components_) delete comp; + + if (jac_sparsity_row_indices != nullptr) + { + delete[] jac_sparsity_row_indices; + } + if (jac_sparsity_col_indices != nullptr) + { + delete[] jac_sparsity_col_indices; + } } /** - * @brief allocator default - * - * @todo this should throw an exception as no allocation without a graph is allowed. - * Or needs to be removed from the base class - * - * @return int + * @brief Allocates and constructs sparsity patterns for system jacobian. + * To do this, sparsity patterns of components are needed, so + * each component's jacobian is evaluated, and the sparsity pattern + * of that component's jacobian is expected to match the sparsity pattern + * for the rest of the simulation. + * @pre `size_` must be correctly initialized */ int allocate() { + // Helper struct for identifying a particular component's local system variable + struct ComponentContribution + { + // The index of the component in `components_`. + size_t comp_idx_; + // The local system variable index + size_t local_row_idx_; + }; + + // Helper struct for saving CSR sparsities of a particular component, since + // right now they only report COO. Can be removed once components have CSR matrices. + struct CsrSparsity + { + std::vector row_indices_; + std::vector col_indices_; + }; + + // A reverse mapping from system variables -> component variables + std::vector> reverse_map(size_); + std::vector component_sparsities; + component_sparsities.reserve(components_.size()); + + // Loop over all components, evaluate their jacobians, save their sparsity information, + // and construct the reverse variable mapping. + for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) + { + component_type* component = components_[comp_idx]; + component->evaluateJacobian(); + + auto [row_indices, col_indices, _] = component->getJacobian().setDataToCSR(); + + component_sparsities.push_back({ + .row_indices_ = std::move(row_indices), + .col_indices_ = std::move(col_indices), + }); + + for (IdxT local_row = 0; local_row < component->size(); local_row++) + { + IdxT global_row = component->getNodeConnection(local_row); + + // Not all local variables map to a global variable + if (global_row != neg1_) + { + reverse_map[global_row].push_back({.comp_idx_ = comp_idx, .local_row_idx_ = local_row}); + } + } + } + + // Allocate the final sparsity pattern info + IdxT* global_row_indices = new IdxT[size_ + 1]; + std::vector global_col_indices; + global_col_indices.reserve(size_); + + // Construct sparsity pattern row-by-row. In the future, can be batched + // into contiguous blocks of rows which can be calculated in parallel, + // then joined sequentially + for (size_t row = 0; row < size_; row++) + { + global_row_indices[row] = global_col_indices.size(); + + // Collect columns from each component which has a row which contributes to this row + for (ComponentContribution contrib : reverse_map[row]) + { + component_type* component = components_[contrib.comp_idx_]; + CsrSparsity& comp_sparsity = component_sparsities[contrib.comp_idx_]; + + IdxT row_start = comp_sparsity.row_indices_[contrib.local_row_idx_]; + IdxT row_end = comp_sparsity.row_indices_[contrib.local_row_idx_ + 1]; + for (size_t local_elem_idx = row_start; local_elem_idx < row_end; local_elem_idx++) + { + IdxT local_col = comp_sparsity.col_indices_[local_elem_idx]; + IdxT global_col = component->getNodeConnection(local_col); + + // Not all columns map to columns in the system jacobian + if (global_col != neg1_) + { + global_col_indices.push_back(global_col); + } + } + } + + // Sort the row by column indices. Since the mapping from local indices to global indices isn't monotonically increasing, + // this is necessary. + auto start = std::next(global_col_indices.begin(), global_row_indices[row]); + std::sort(start, global_col_indices.end()); + + // If there were multiple components which contributed columns, then we definitely need + // to de-duplicate the columns. + if (reverse_map[row].size() > 1) + { + global_col_indices.erase(std::unique(start, global_col_indices.end()), global_col_indices.end()); + } + } + // Final row index (beyond the end) keeps track of nnz, which is also the size of the column indices. + global_row_indices[size_] = global_col_indices.size(); + + // Delete old sparsity buffers, if they exist + if (jac_sparsity_row_indices != nullptr) + { + delete[] jac_sparsity_row_indices; + } + if (jac_sparsity_col_indices != nullptr) + { + delete[] jac_sparsity_col_indices; + } + + // Allocate new sparsity buffers + jac_sparsity_row_indices = global_row_indices; + jac_sparsity_col_indices = new IdxT[global_col_indices.size()]; + + // Copy column indices + std::copy(global_col_indices.begin(), global_col_indices.end(), jac_sparsity_col_indices); + return 1; } @@ -184,6 +310,8 @@ namespace GridKit yp_.resize(size_); f_.resize(size_); + allocate(); + return 0; } @@ -401,9 +529,42 @@ namespace GridKit components_.push_back(component); } + /** + * @brief Returns a pointer to the buffer of Jacobian CSR row indices. + * @todo Can be removed once `getJacobian()` returns a `CsrMatrix` + */ + IdxT* getJacRowIndices() + { + return jac_sparsity_row_indices; + } + + /** + * @brief Returns a pointer to the buffer of Jacobian CSR column indices. + * @todo Can be removed once `getJacobian()` returns a `CsrMatrix` + */ + IdxT* getJacColIndices() + { + return jac_sparsity_col_indices; + } + private: static constexpr IdxT neg1_ = static_cast(-1); + /** + * @brief Keeps track of the CSR row indices of the system jacobian. + * `nullptr` is used to indicate an un-allocated buffer. + * @todo Unneeded once the jacobian is in CSR format. This can be stored + * in the jacobian itself. + */ + IdxT* jac_sparsity_row_indices; + /** + * @brief Keeps track of the CSR col indices of the system jacobian. + * `nullptr` is used to indicate an un-allocated buffer. + * @todo Unneeded once the jacobian is in CSR format. This can be stored + * in the jacobian itself. + */ + IdxT* jac_sparsity_col_indices; + std::vector components_; int jac_call_count_{0}; diff --git a/tests/UnitTests/CMakeLists.txt b/tests/UnitTests/CMakeLists.txt index c002660b..3c7e0892 100644 --- a/tests/UnitTests/CMakeLists.txt +++ b/tests/UnitTests/CMakeLists.txt @@ -2,5 +2,6 @@ add_subdirectory(AutomaticDifferentiation) add_subdirectory(LinearAlgebra) add_subdirectory(PhasorDynamics) +add_subdirectory(PowerElectronics) add_subdirectory(Solver) add_subdirectory(Utilities) diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt new file mode 100644 index 00000000..b8ea63c1 --- /dev/null +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -0,0 +1,11 @@ +add_executable(test_power_elec_system runSystemTests.cpp) +target_link_libraries(test_power_elec_system GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::solvers_dyn + GridKit::power_elec_microbusdq) + +add_test(NAME PowerElectronicsSystemTest COMMAND test_power_elec_system) + +install(TARGETS test_power_elec_system + RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/PowerElectronics/SystemTests.hpp b/tests/UnitTests/PowerElectronics/SystemTests.hpp new file mode 100644 index 00000000..618346d3 --- /dev/null +++ b/tests/UnitTests/PowerElectronics/SystemTests.hpp @@ -0,0 +1,373 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + + template + class SystemTests + { + public: + SystemTests() = default; + ~SystemTests() = default; + + TestOutcome jacobianSparsity() + { + TestStatus success = true; + + double abs_tol = 1.0e-8; + double rel_tol = 1.0e-8; + size_t max_step_number = 3000; + bool use_jac = true; + + // Create model + auto* sysmodel = new GridKit::PowerElectronicsModel(rel_tol, abs_tol, use_jac, max_step_number); + + // Modeled after the problem in the paper + double RN = 1.0e4; + + // DG Params + GridKit::DistributedGeneratorParameters parms1; + parms1.wb_ = 2.0 * M_PI * 50.0; + parms1.wc_ = 31.41; + parms1.mp_ = 9.4e-5; + parms1.Vn_ = 380.0; + parms1.nq_ = 1.3e-3; + parms1.F_ = 0.75; + parms1.Kiv_ = 420.0; + parms1.Kpv_ = 0.1; + parms1.Kic_ = 2.0e4; + parms1.Kpc_ = 15.0; + parms1.Cf_ = 5.0e-5; + parms1.rLf_ = 0.1; + parms1.Lf_ = 1.35e-3; + parms1.rLc_ = 0.03; + parms1.Lc_ = 0.35e-3; + + GridKit::DistributedGeneratorParameters parms2; + // Parameters from MATLAB Microgrid code for first DG + parms2.wb_ = 2.0 * M_PI * 50.0; + parms2.wc_ = 31.41; + parms2.mp_ = 12.5e-5; + parms2.Vn_ = 380.0; + parms2.nq_ = 1.5e-3; + parms2.F_ = 0.75; + parms2.Kiv_ = 390.0; + parms2.Kpv_ = 0.05; + parms2.Kic_ = 16.0e3; + parms2.Kpc_ = 10.5; + parms2.Cf_ = 50.0e-6; + parms2.rLf_ = 0.1; + parms2.Lf_ = 1.35e-3; + parms2.rLc_ = 0.03; + parms2.Lc_ = 0.35e-3; + + // Line params + double rline1 = 0.23; + double Lline1 = 0.1 / (2.0 * M_PI * 50.0); + + double rline2 = 0.35; + double Lline2 = 0.58 / (2.0 * M_PI * 50.0); + + double rline3 = 0.23; + double Lline3 = 0.1 / (2.0 * M_PI * 50.0); + + // load parms + double rload1 = 3.0; + double Lload1 = 2.0 / (2.0 * M_PI * 50.0); + + double rload2 = 2.0; + double Lload2 = 1.0 / (2.0 * M_PI * 50.0); + + // indexing sets + size_t Nsize = 2; + // DGs + - refframe Lines + Loads + size_t vec_size_internals = 13 * (2 * Nsize) - 1 + (2 + 4 * (Nsize - 1)) + 2 * Nsize; + // \omegaref + BusDQ + size_t vec_size_externals = 1 + 2 * (2 * Nsize); + size_t dqbus1 = vec_size_internals + 1; + size_t dqbus2 = vec_size_internals + 3; + size_t dqbus3 = vec_size_internals + 5; + size_t dqbus4 = vec_size_internals + 7; + + size_t vec_size_total = vec_size_internals + vec_size_externals; + + size_t indexv = 0; + + // dg 1 + GridKit::DistributedGenerator* dg1 = new GridKit::DistributedGenerator(0, parms1, true); + // ref motor + dg1->setExternalConnectionNodes(0, vec_size_internals); + // outputs + dg1->setExternalConnectionNodes(1, dqbus1); + dg1->setExternalConnectionNodes(2, dqbus1 + 1); + //"grounding" of the difference + dg1->setExternalConnectionNodes(3, static_cast(-1)); + // internal connections + for (size_t i = 0; i < 12; i++) + { + dg1->setExternalConnectionNodes(4 + i, indexv + i); + } + indexv += 12; + sysmodel->addComponent(dg1); + + // dg 2 + GridKit::DistributedGenerator* dg2 = new GridKit::DistributedGenerator(1, parms1, false); + // ref motor + dg2->setExternalConnectionNodes(0, vec_size_internals); + // outputs + dg2->setExternalConnectionNodes(1, dqbus2); + dg2->setExternalConnectionNodes(2, dqbus2 + 1); + // internal connections + for (size_t i = 0; i < 13; i++) + { + dg2->setExternalConnectionNodes(3 + i, indexv + i); + } + indexv += 13; + sysmodel->addComponent(dg2); + + // dg 3 + GridKit::DistributedGenerator* dg3 = new GridKit::DistributedGenerator(2, parms2, false); + // ref motor + dg3->setExternalConnectionNodes(0, vec_size_internals); + // outputs + dg3->setExternalConnectionNodes(1, dqbus3); + dg3->setExternalConnectionNodes(2, dqbus3 + 1); + // internal connections + for (size_t i = 0; i < 13; i++) + { + dg3->setExternalConnectionNodes(3 + i, indexv + i); + } + indexv += 13; + sysmodel->addComponent(dg3); + + // dg 4 + GridKit::DistributedGenerator* dg4 = new GridKit::DistributedGenerator(3, parms2, false); + // ref motor + dg4->setExternalConnectionNodes(0, vec_size_internals); + // outputs + dg4->setExternalConnectionNodes(1, dqbus4); + dg4->setExternalConnectionNodes(2, dqbus4 + 1); + + // internal connections + for (size_t i = 0; i < 13; i++) + { + dg4->setExternalConnectionNodes(3 + i, indexv + i); + } + indexv += 13; + sysmodel->addComponent(dg4); + + // Lines + + // line 1 + GridKit::MicrogridLine* l1 = new GridKit::MicrogridLine(4, rline1, Lline1); + // ref motor + l1->setExternalConnectionNodes(0, vec_size_internals); + // input connections + l1->setExternalConnectionNodes(1, dqbus1); + l1->setExternalConnectionNodes(2, dqbus1 + 1); + // output connections + l1->setExternalConnectionNodes(3, dqbus2); + l1->setExternalConnectionNodes(4, dqbus2 + 1); + // internal connections + for (size_t i = 0; i < 2; i++) + { + l1->setExternalConnectionNodes(5 + i, indexv + i); + } + indexv += 2; + sysmodel->addComponent(l1); + + // line 2 + GridKit::MicrogridLine* l2 = new GridKit::MicrogridLine(5, rline2, Lline2); + // ref motor + l2->setExternalConnectionNodes(0, vec_size_internals); + // input connections + l2->setExternalConnectionNodes(1, dqbus2); + l2->setExternalConnectionNodes(2, dqbus2 + 1); + // output connections + l2->setExternalConnectionNodes(3, dqbus3); + l2->setExternalConnectionNodes(4, dqbus3 + 1); + // internal connections + for (size_t i = 0; i < 2; i++) + { + l2->setExternalConnectionNodes(5 + i, indexv + i); + } + indexv += 2; + sysmodel->addComponent(l2); + + // line 3 + GridKit::MicrogridLine* l3 = new GridKit::MicrogridLine(6, rline3, Lline3); + // ref motor + l3->setExternalConnectionNodes(0, vec_size_internals); + // input connections + l3->setExternalConnectionNodes(1, dqbus3); + l3->setExternalConnectionNodes(2, dqbus3 + 1); + // output connections + l3->setExternalConnectionNodes(3, dqbus4); + l3->setExternalConnectionNodes(4, dqbus4 + 1); + // internal connections + for (size_t i = 0; i < 2; i++) + { + l3->setExternalConnectionNodes(5 + i, indexv + i); + } + indexv += 2; + sysmodel->addComponent(l3); + + // loads + + // load 1 + GridKit::MicrogridLoad* load1 = new GridKit::MicrogridLoad(7, rload1, Lload1); + // ref motor + load1->setExternalConnectionNodes(0, vec_size_internals); + // input connections + load1->setExternalConnectionNodes(1, dqbus1); + load1->setExternalConnectionNodes(2, dqbus1 + 1); + // internal connections + for (size_t i = 0; i < 2; i++) + { + load1->setExternalConnectionNodes(3 + i, indexv + i); + } + indexv += 2; + sysmodel->addComponent(load1); + + // load 2 + GridKit::MicrogridLoad* load2 = new GridKit::MicrogridLoad(8, rload2, Lload2); + // ref motor + load2->setExternalConnectionNodes(0, vec_size_internals); + // input connections + load2->setExternalConnectionNodes(1, dqbus3); + load2->setExternalConnectionNodes(2, dqbus3 + 1); + // internal connections + for (size_t i = 0; i < 2; i++) + { + load2->setExternalConnectionNodes(3 + i, indexv + i); + } + indexv += 2; + sysmodel->addComponent(load2); + + // Virtual PQ Buses + GridKit::MicrogridBusDQ* bus1 = new GridKit::MicrogridBusDQ(9, RN); + + bus1->setExternalConnectionNodes(0, dqbus1); + bus1->setExternalConnectionNodes(1, dqbus1 + 1); + sysmodel->addComponent(bus1); + + GridKit::MicrogridBusDQ* bus2 = new GridKit::MicrogridBusDQ(10, RN); + + bus2->setExternalConnectionNodes(0, dqbus2); + bus2->setExternalConnectionNodes(1, dqbus2 + 1); + sysmodel->addComponent(bus2); + + GridKit::MicrogridBusDQ* bus3 = new GridKit::MicrogridBusDQ(11, RN); + + bus3->setExternalConnectionNodes(0, dqbus3); + bus3->setExternalConnectionNodes(1, dqbus3 + 1); + sysmodel->addComponent(bus3); + + GridKit::MicrogridBusDQ* bus4 = new GridKit::MicrogridBusDQ(12, RN); + + bus4->setExternalConnectionNodes(0, dqbus4); + bus4->setExternalConnectionNodes(1, dqbus4 + 1); + sysmodel->addComponent(bus4); + + sysmodel->allocate(vec_size_total); + + // Create Intial points for states + for (size_t i = 0; i < vec_size_total; i++) + { + sysmodel->y()[i] = 0.0; + sysmodel->yp()[i] = 0.0; + sysmodel->y()[i].setVariableNumber(i); + } + + // Create Intial derivatives specifics generated in MATLAB + // DGs 1 + sysmodel->yp()[2] = parms1.Vn_; + sysmodel->yp()[4] = parms1.Kpv_ * parms1.Vn_; + sysmodel->yp()[6] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_; + sysmodel->yp()[12 + 3] = parms1.Vn_; + sysmodel->yp()[12 + 5] = parms1.Kpv_ * parms1.Vn_; + sysmodel->yp()[12 + 7] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_; + for (size_t i = 2; i < 4; i++) + { + sysmodel->yp()[13 * i - 1 + 3] = parms2.Vn_; + sysmodel->yp()[13 * i - 1 + 5] = parms2.Kpv_ * parms2.Vn_; + sysmodel->yp()[13 * i - 1 + 7] = (parms2.Kpc_ * parms2.Kpv_ * parms2.Vn_) / parms2.Lf_; + } + + // since the intial P_com = 0 + sysmodel->y()[vec_size_internals] = parms1.wb_; + sysmodel->y()[vec_size_internals].setVariableNumber(vec_size_internals); + + for (size_t i = 0; i < vec_size_total; i++) + { + sysmodel->yp()[i].setVariableNumber(i); + } + + sysmodel->initialize(); + sysmodel->evaluateResidual(); + + std::vector& residuals = sysmodel->getResidual(); + + size_t* row_indices = sysmodel->getJacRowIndices(); + size_t* col_indices = sysmodel->getJacColIndices(); + + for (size_t row = 0; row < residuals.size(); row++) + { + bool row_correct = true; + + DependencyTracking::Variable& row_residual = residuals[row]; + auto& dependency_map = row_residual.getDependencies(); + auto dependencies = std::vector(dependency_map.begin(), dependency_map.end()); + + // Verify both rows are the same size + row_correct &= (row_indices[row + 1] - row_indices[row]) == dependencies.size(); + + if (row_correct) + { + for (size_t col_idx = 0; col_idx < dependencies.size(); col_idx++) + { + row_correct &= col_indices[col_idx + row_indices[row]] == std::get<0>(dependencies[col_idx]); + } + } + + if (!row_correct) + { + std::cerr << "Row " << row << " is incorrect:\n"; + std::cerr << " CSR Row: {"; + for (size_t col_idx = row_indices[row]; col_idx < row_indices[row + 1]; col_idx++) + { + std::cerr << std::format("{:3}{}", col_indices[col_idx], col_idx == row_indices[row + 1] - 1 ? "" : ", "); + } + std::cerr << "}\n" + << " Dep Tracker Row: {"; + for (size_t col_idx = 0; col_idx < dependencies.size(); col_idx++) + { + std::cerr << std::format("{:3}{}", std::get<0>(dependencies[col_idx]), col_idx == dependencies.size() - 1 ? "" : ", "); + } + std::cerr << "}\n\n"; + } + + success *= row_correct; + } + + delete sysmodel; + + return success.report(__func__); + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PowerElectronics/runSystemTests.cpp b/tests/UnitTests/PowerElectronics/runSystemTests.cpp new file mode 100644 index 00000000..d924e1b8 --- /dev/null +++ b/tests/UnitTests/PowerElectronics/runSystemTests.cpp @@ -0,0 +1,14 @@ +#include "SystemTests.hpp" + +int main() +{ + using namespace GridKit; + using namespace GridKit::Testing; + + GridKit::Testing::TestingResults result; + GridKit::Testing::SystemTests test; + + result += test.jacobianSparsity(); + + return result.summary(); +} From 6edf7d034f9af519daad3a861737c5425687ce65 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 11 Nov 2025 14:17:44 -0500 Subject: [PATCH 03/29] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2999c732..60eba3d9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,6 +37,7 @@ - Added support for DependencyTracking::Variable in PowerElectronics models. - Updated Jacobian value storage from `ScalarT` to `RealT`. - Added a header file defining constants to be used throughout the code. +- Added Jacobian sparsity pattern computation into `PowerElectronics`. ## v0.1 From 75bc32e0417e62e1a94a2e2dda4caf52106911c0 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 11 Nov 2025 14:25:28 -0500 Subject: [PATCH 04/29] Remove unneeded dependency --- tests/UnitTests/PowerElectronics/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt index b8ea63c1..556de655 100644 --- a/tests/UnitTests/PowerElectronics/CMakeLists.txt +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -2,7 +2,6 @@ add_executable(test_power_elec_system runSystemTests.cpp) target_link_libraries(test_power_elec_system GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload - GridKit::solvers_dyn GridKit::power_elec_microbusdq) add_test(NAME PowerElectronicsSystemTest COMMAND test_power_elec_system) From 669d9898396175e10d8d3af46cbdbb4ecb390aa1 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 11 Nov 2025 16:06:14 -0500 Subject: [PATCH 05/29] Remove unnecessary this-> --- .../SystemModelPowerElectronics.hpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 00e1aac8..745604dc 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -71,6 +71,7 @@ namespace GridKit using CircuitComponent::jac_; using CircuitComponent::rel_tol_; using CircuitComponent::abs_tol_; + using CircuitComponent::max_steps_; public: /** @@ -83,7 +84,7 @@ namespace GridKit // Set system model parameters as default rel_tol_ = 1e-4; abs_tol_ = 1e-4; - this->max_steps_ = 2000; + max_steps_ = 2000; // By default don't use the jacobian use_jac_ = false; jac_sparsity_row_indices = nullptr; @@ -108,7 +109,7 @@ namespace GridKit // Set system model tolerances from input rel_tol_ = rel_tol; abs_tol_ = abs_tol; - this->max_steps_ = max_steps; + max_steps_ = max_steps; // Can choose if to use jacobian use_jac_ = use_jac; jac_sparsity_row_indices = nullptr; @@ -125,7 +126,7 @@ namespace GridKit */ virtual ~PowerElectronicsModel() { - for (auto comp : this->components_) + for (auto comp : components_) delete comp; if (jac_sparsity_row_indices != nullptr) @@ -273,7 +274,7 @@ namespace GridKit */ bool hasJacobian() { - if (!this->use_jac_) + if (!use_jac_) return false; for (const auto& component : components_) @@ -327,7 +328,7 @@ namespace GridKit { component->initialize(); } - this->distributeVectors(); + distributeVectors(); return 0; } @@ -372,12 +373,12 @@ namespace GridKit */ int evaluateResidual() { - for (IdxT i = 0; i < this->f_.size(); i++) + for (IdxT i = 0; i < f_.size(); i++) { f_[i] = 0.0; } - this->distributeVectors(); + distributeVectors(); // Update system residual vector From 759f22f5c03e7a56f760bf383cd047039de5d04e Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 11 Nov 2025 16:21:38 -0500 Subject: [PATCH 06/29] Remove std::vector for component_sparsities --- .../PowerElectronics/SystemModelPowerElectronics.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 745604dc..b3b0b07d 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -168,8 +168,7 @@ namespace GridKit // A reverse mapping from system variables -> component variables std::vector> reverse_map(size_); - std::vector component_sparsities; - component_sparsities.reserve(components_.size()); + CsrSparsity* component_sparsities = new CsrSparsity[components_.size()]; // Loop over all components, evaluate their jacobians, save their sparsity information, // and construct the reverse variable mapping. @@ -180,10 +179,10 @@ namespace GridKit auto [row_indices, col_indices, _] = component->getJacobian().setDataToCSR(); - component_sparsities.push_back({ + component_sparsities[comp_idx] = CsrSparsity{ .row_indices_ = std::move(row_indices), .col_indices_ = std::move(col_indices), - }); + }; for (IdxT local_row = 0; local_row < component->size(); local_row++) { @@ -262,6 +261,8 @@ namespace GridKit // Copy column indices std::copy(global_col_indices.begin(), global_col_indices.end(), jac_sparsity_col_indices); + delete[] component_sparsities; + return 1; } From c0c6c379b9cdb313c819d421dfda59182c72fc0b Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 11:54:22 -0500 Subject: [PATCH 07/29] Update member names for guidelines --- .../SystemModelPowerElectronics.hpp | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index b3b0b07d..ec221442 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -82,13 +82,13 @@ namespace GridKit PowerElectronicsModel() { // Set system model parameters as default - rel_tol_ = 1e-4; - abs_tol_ = 1e-4; - max_steps_ = 2000; + rel_tol_ = 1e-4; + abs_tol_ = 1e-4; + max_steps_ = 2000; // By default don't use the jacobian - use_jac_ = false; - jac_sparsity_row_indices = nullptr; - jac_sparsity_col_indices = nullptr; + use_jac_ = false; + jac_sparsity_row_indices_ = nullptr; + jac_sparsity_col_indices_ = nullptr; } /** @@ -107,13 +107,13 @@ namespace GridKit IdxT max_steps = 2000) { // Set system model tolerances from input - rel_tol_ = rel_tol; - abs_tol_ = abs_tol; - max_steps_ = max_steps; + rel_tol_ = rel_tol; + abs_tol_ = abs_tol; + max_steps_ = max_steps; // Can choose if to use jacobian - use_jac_ = use_jac; - jac_sparsity_row_indices = nullptr; - jac_sparsity_col_indices = nullptr; + use_jac_ = use_jac; + jac_sparsity_row_indices_ = nullptr; + jac_sparsity_col_indices_ = nullptr; } /** @@ -129,13 +129,13 @@ namespace GridKit for (auto comp : components_) delete comp; - if (jac_sparsity_row_indices != nullptr) + if (jac_sparsity_row_indices_ != nullptr) { - delete[] jac_sparsity_row_indices; + delete[] jac_sparsity_row_indices_; } - if (jac_sparsity_col_indices != nullptr) + if (jac_sparsity_col_indices_ != nullptr) { - delete[] jac_sparsity_col_indices; + delete[] jac_sparsity_col_indices_; } } @@ -245,21 +245,21 @@ namespace GridKit global_row_indices[size_] = global_col_indices.size(); // Delete old sparsity buffers, if they exist - if (jac_sparsity_row_indices != nullptr) + if (jac_sparsity_row_indices_ != nullptr) { - delete[] jac_sparsity_row_indices; + delete[] jac_sparsity_row_indices_; } - if (jac_sparsity_col_indices != nullptr) + if (jac_sparsity_col_indices_ != nullptr) { - delete[] jac_sparsity_col_indices; + delete[] jac_sparsity_col_indices_; } // Allocate new sparsity buffers - jac_sparsity_row_indices = global_row_indices; - jac_sparsity_col_indices = new IdxT[global_col_indices.size()]; + jac_sparsity_row_indices_ = global_row_indices; + jac_sparsity_col_indices_ = new IdxT[global_col_indices.size()]; // Copy column indices - std::copy(global_col_indices.begin(), global_col_indices.end(), jac_sparsity_col_indices); + std::copy(global_col_indices.begin(), global_col_indices.end(), jac_sparsity_col_indices_); delete[] component_sparsities; @@ -537,7 +537,7 @@ namespace GridKit */ IdxT* getJacRowIndices() { - return jac_sparsity_row_indices; + return jac_sparsity_row_indices_; } /** @@ -546,7 +546,7 @@ namespace GridKit */ IdxT* getJacColIndices() { - return jac_sparsity_col_indices; + return jac_sparsity_col_indices_; } private: @@ -558,14 +558,14 @@ namespace GridKit * @todo Unneeded once the jacobian is in CSR format. This can be stored * in the jacobian itself. */ - IdxT* jac_sparsity_row_indices; + IdxT* jac_sparsity_row_indices_; /** * @brief Keeps track of the CSR col indices of the system jacobian. * `nullptr` is used to indicate an un-allocated buffer. * @todo Unneeded once the jacobian is in CSR format. This can be stored * in the jacobian itself. */ - IdxT* jac_sparsity_col_indices; + IdxT* jac_sparsity_col_indices_; std::vector components_; From 963c12756f3bc946aecfefffa3f0936dabb6187a Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 11:54:51 -0500 Subject: [PATCH 08/29] Change reverse_map buckets to linked list --- .../SystemModelPowerElectronics.hpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index ec221442..ef61d975 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -3,6 +3,7 @@ #pragma once #include +#include #include #include #include @@ -167,8 +168,8 @@ namespace GridKit }; // A reverse mapping from system variables -> component variables - std::vector> reverse_map(size_); - CsrSparsity* component_sparsities = new CsrSparsity[components_.size()]; + auto reverse_map = new std::forward_list[size_]; + auto component_sparsities = new CsrSparsity[components_.size()]; // Loop over all components, evaluate their jacobians, save their sparsity information, // and construct the reverse variable mapping. @@ -191,7 +192,7 @@ namespace GridKit // Not all local variables map to a global variable if (global_row != neg1_) { - reverse_map[global_row].push_back({.comp_idx_ = comp_idx, .local_row_idx_ = local_row}); + reverse_map[global_row].push_front({.comp_idx_ = comp_idx, .local_row_idx_ = local_row}); } } } @@ -208,12 +209,17 @@ namespace GridKit { global_row_indices[row] = global_col_indices.size(); + // How many components contribute to this row + size_t num_contrib = 0; + // Collect columns from each component which has a row which contributes to this row for (ComponentContribution contrib : reverse_map[row]) { component_type* component = components_[contrib.comp_idx_]; CsrSparsity& comp_sparsity = component_sparsities[contrib.comp_idx_]; + num_contrib++; + IdxT row_start = comp_sparsity.row_indices_[contrib.local_row_idx_]; IdxT row_end = comp_sparsity.row_indices_[contrib.local_row_idx_ + 1]; for (size_t local_elem_idx = row_start; local_elem_idx < row_end; local_elem_idx++) @@ -236,7 +242,7 @@ namespace GridKit // If there were multiple components which contributed columns, then we definitely need // to de-duplicate the columns. - if (reverse_map[row].size() > 1) + if (num_contrib > 1) { global_col_indices.erase(std::unique(start, global_col_indices.end()), global_col_indices.end()); } @@ -262,6 +268,7 @@ namespace GridKit std::copy(global_col_indices.begin(), global_col_indices.end(), jac_sparsity_col_indices_); delete[] component_sparsities; + delete[] reverse_map; return 1; } From 2cd92485f7a669c873112c2405899bea97ad5d56 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 11:58:51 -0500 Subject: [PATCH 09/29] Change sparsity test to use scalable microgrid --- .../PowerElectronics/SystemTests.hpp | 376 +++++++----------- .../PowerElectronics/runSystemTests.cpp | 3 +- 2 files changed, 154 insertions(+), 225 deletions(-) diff --git a/tests/UnitTests/PowerElectronics/SystemTests.hpp b/tests/UnitTests/PowerElectronics/SystemTests.hpp index 618346d3..1ea52759 100644 --- a/tests/UnitTests/PowerElectronics/SystemTests.hpp +++ b/tests/UnitTests/PowerElectronics/SystemTests.hpp @@ -24,7 +24,7 @@ namespace GridKit SystemTests() = default; ~SystemTests() = default; - TestOutcome jacobianSparsity() + TestOutcome jacobianSparsity(size_t Nsize) { TestStatus success = true; @@ -39,249 +39,183 @@ namespace GridKit // Modeled after the problem in the paper double RN = 1.0e4; - // DG Params - GridKit::DistributedGeneratorParameters parms1; - parms1.wb_ = 2.0 * M_PI * 50.0; - parms1.wc_ = 31.41; - parms1.mp_ = 9.4e-5; - parms1.Vn_ = 380.0; - parms1.nq_ = 1.3e-3; - parms1.F_ = 0.75; - parms1.Kiv_ = 420.0; - parms1.Kpv_ = 0.1; - parms1.Kic_ = 2.0e4; - parms1.Kpc_ = 15.0; - parms1.Cf_ = 5.0e-5; - parms1.rLf_ = 0.1; - parms1.Lf_ = 1.35e-3; - parms1.rLc_ = 0.03; - parms1.Lc_ = 0.35e-3; - - GridKit::DistributedGeneratorParameters parms2; - // Parameters from MATLAB Microgrid code for first DG - parms2.wb_ = 2.0 * M_PI * 50.0; - parms2.wc_ = 31.41; - parms2.mp_ = 12.5e-5; - parms2.Vn_ = 380.0; - parms2.nq_ = 1.5e-3; - parms2.F_ = 0.75; - parms2.Kiv_ = 390.0; - parms2.Kpv_ = 0.05; - parms2.Kic_ = 16.0e3; - parms2.Kpc_ = 10.5; - parms2.Cf_ = 50.0e-6; - parms2.rLf_ = 0.1; - parms2.Lf_ = 1.35e-3; - parms2.rLc_ = 0.03; - parms2.Lc_ = 0.35e-3; - - // Line params - double rline1 = 0.23; - double Lline1 = 0.1 / (2.0 * M_PI * 50.0); - - double rline2 = 0.35; - double Lline2 = 0.58 / (2.0 * M_PI * 50.0); - - double rline3 = 0.23; - double Lline3 = 0.1 / (2.0 * M_PI * 50.0); + // DG Params Vector + // All DGs have the same set of parameters except for the first two. + GridKit::DistributedGeneratorParameters DG_parms1; + DG_parms1.wb_ = 2.0 * M_PI * 50.0; + DG_parms1.wc_ = 31.41; + DG_parms1.mp_ = 9.4e-5; + DG_parms1.Vn_ = 380.0; + DG_parms1.nq_ = 1.3e-3; + DG_parms1.F_ = 0.75; + DG_parms1.Kiv_ = 420.0; + DG_parms1.Kpv_ = 0.1; + DG_parms1.Kic_ = 2.0e4; + DG_parms1.Kpc_ = 15.0; + DG_parms1.Cf_ = 5.0e-5; + DG_parms1.rLf_ = 0.1; + DG_parms1.Lf_ = 1.35e-3; + DG_parms1.rLc_ = 0.03; + DG_parms1.Lc_ = 0.35e-3; + + GridKit::DistributedGeneratorParameters DG_parms2; + DG_parms2.wb_ = 2.0 * M_PI * 50.0; + DG_parms2.wc_ = 31.41; + DG_parms2.mp_ = 12.5e-5; + DG_parms2.Vn_ = 380.0; + DG_parms2.nq_ = 1.5e-3; + DG_parms2.F_ = 0.75; + DG_parms2.Kiv_ = 390.0; + DG_parms2.Kpv_ = 0.05; + DG_parms2.Kic_ = 16.0e3; + DG_parms2.Kpc_ = 10.5; + DG_parms2.Cf_ = 50.0e-6; + DG_parms2.rLf_ = 0.1; + DG_parms2.Lf_ = 1.35e-3; + DG_parms2.rLc_ = 0.03; + DG_parms2.Lc_ = 0.35e-3; + + std::vector> DGParams_list(2 * Nsize, DG_parms2); + + DGParams_list[0] = DG_parms1; + DGParams_list[1] = DG_parms1; + + // line vector params + // Every odd line has the same parameters and every even line has the same parameters + double rline1 = 0.23; + double Lline1 = 0.1 / (2.0 * M_PI * 50.0); + double rline2 = 0.35; + double Lline2 = 0.58 / (2.0 * M_PI * 50.0); + std::vector rline_list(2 * Nsize - 1, 0.0); + std::vector Lline_list(2 * Nsize - 1, 0.0); + for (size_t i = 0; i < rline_list.size(); i++) + { + rline_list[i] = (i % 2) ? rline2 : rline1; + Lline_list[i] = (i % 2) ? Lline2 : Lline1; + } // load parms + // Only the first load has the same paramaters. double rload1 = 3.0; double Lload1 = 2.0 / (2.0 * M_PI * 50.0); - double rload2 = 2.0; double Lload2 = 1.0 / (2.0 * M_PI * 50.0); - // indexing sets - size_t Nsize = 2; + std::vector rload_list(Nsize, rload2); + std::vector Lload_list(Nsize, Lload2); + rload_list[0] = rload1; + Lload_list[0] = Lload1; + // DGs + - refframe Lines + Loads size_t vec_size_internals = 13 * (2 * Nsize) - 1 + (2 + 4 * (Nsize - 1)) + 2 * Nsize; // \omegaref + BusDQ size_t vec_size_externals = 1 + 2 * (2 * Nsize); - size_t dqbus1 = vec_size_internals + 1; - size_t dqbus2 = vec_size_internals + 3; - size_t dqbus3 = vec_size_internals + 5; - size_t dqbus4 = vec_size_internals + 7; - - size_t vec_size_total = vec_size_internals + vec_size_externals; - - size_t indexv = 0; - // dg 1 - GridKit::DistributedGenerator* dg1 = new GridKit::DistributedGenerator(0, parms1, true); - // ref motor - dg1->setExternalConnectionNodes(0, vec_size_internals); - // outputs - dg1->setExternalConnectionNodes(1, dqbus1); - dg1->setExternalConnectionNodes(2, dqbus1 + 1); - //"grounding" of the difference - dg1->setExternalConnectionNodes(3, static_cast(-1)); - // internal connections - for (size_t i = 0; i < 12; i++) + std::vector vdqbus_index(2 * Nsize, 0); + vdqbus_index[0] = vec_size_internals + 1; + for (size_t i = 1; i < vdqbus_index.size(); i++) { - dg1->setExternalConnectionNodes(4 + i, indexv + i); + vdqbus_index[i] = vdqbus_index[i - 1] + 2; } - indexv += 12; - sysmodel->addComponent(dg1); - // dg 2 - GridKit::DistributedGenerator* dg2 = new GridKit::DistributedGenerator(1, parms1, false); - // ref motor - dg2->setExternalConnectionNodes(0, vec_size_internals); - // outputs - dg2->setExternalConnectionNodes(1, dqbus2); - dg2->setExternalConnectionNodes(2, dqbus2 + 1); - // internal connections - for (size_t i = 0; i < 13; i++) - { - dg2->setExternalConnectionNodes(3 + i, indexv + i); - } - indexv += 13; - sysmodel->addComponent(dg2); + // Total size of the vector setup + size_t vec_size_total = vec_size_internals + vec_size_externals; - // dg 3 - GridKit::DistributedGenerator* dg3 = new GridKit::DistributedGenerator(2, parms2, false); + // Create the reference DG + auto* dg_ref = new DistributedGenerator(0, + DGParams_list[0], + true); // ref motor - dg3->setExternalConnectionNodes(0, vec_size_internals); + dg_ref->setExternalConnectionNodes(0, vec_size_internals); // outputs - dg3->setExternalConnectionNodes(1, dqbus3); - dg3->setExternalConnectionNodes(2, dqbus3 + 1); + dg_ref->setExternalConnectionNodes(1, vdqbus_index[0]); + dg_ref->setExternalConnectionNodes(2, vdqbus_index[0] + 1); + //"grounding" of the difference + dg_ref->setExternalConnectionNodes(3, static_cast(-1)); // internal connections - for (size_t i = 0; i < 13; i++) + for (size_t i = 0; i < 12; i++) { - dg3->setExternalConnectionNodes(3 + i, indexv + i); + dg_ref->setExternalConnectionNodes(4 + i, i); } - indexv += 13; - sysmodel->addComponent(dg3); - - // dg 4 - GridKit::DistributedGenerator* dg4 = new GridKit::DistributedGenerator(3, parms2, false); - // ref motor - dg4->setExternalConnectionNodes(0, vec_size_internals); - // outputs - dg4->setExternalConnectionNodes(1, dqbus4); - dg4->setExternalConnectionNodes(2, dqbus4 + 1); + sysmodel->addComponent(dg_ref); - // internal connections - for (size_t i = 0; i < 13; i++) + // Keep track of models and index location + size_t indexv = 12; + size_t model_id = 1; + // Add all other DGs + for (size_t i = 1; i < 2 * Nsize; i++) { - dg4->setExternalConnectionNodes(3 + i, indexv + i); - } - indexv += 13; - sysmodel->addComponent(dg4); - - // Lines - - // line 1 - GridKit::MicrogridLine* l1 = new GridKit::MicrogridLine(4, rline1, Lline1); - // ref motor - l1->setExternalConnectionNodes(0, vec_size_internals); - // input connections - l1->setExternalConnectionNodes(1, dqbus1); - l1->setExternalConnectionNodes(2, dqbus1 + 1); - // output connections - l1->setExternalConnectionNodes(3, dqbus2); - l1->setExternalConnectionNodes(4, dqbus2 + 1); - // internal connections - for (size_t i = 0; i < 2; i++) - { - l1->setExternalConnectionNodes(5 + i, indexv + i); + // current DG to add + auto* dg = new DistributedGenerator(model_id++, + DGParams_list[i], + false); + // ref motor + dg->setExternalConnectionNodes(0, vec_size_internals); + // outputs + dg->setExternalConnectionNodes(1, vdqbus_index[i]); + dg->setExternalConnectionNodes(2, vdqbus_index[i] + 1); + // internal connections + for (size_t j = 0; j < 13; j++) + { + dg->setExternalConnectionNodes(3 + j, indexv + j); + } + indexv += 13; + sysmodel->addComponent(dg); } - indexv += 2; - sysmodel->addComponent(l1); - // line 2 - GridKit::MicrogridLine* l2 = new GridKit::MicrogridLine(5, rline2, Lline2); - // ref motor - l2->setExternalConnectionNodes(0, vec_size_internals); - // input connections - l2->setExternalConnectionNodes(1, dqbus2); - l2->setExternalConnectionNodes(2, dqbus2 + 1); - // output connections - l2->setExternalConnectionNodes(3, dqbus3); - l2->setExternalConnectionNodes(4, dqbus3 + 1); - // internal connections - for (size_t i = 0; i < 2; i++) + // Load all the Line compoenents + for (size_t i = 0; i < 2 * Nsize - 1; i++) { - l2->setExternalConnectionNodes(5 + i, indexv + i); + // line + auto* line_model = new MicrogridLine(model_id++, + rline_list[i], + Lline_list[i]); + // ref motor + line_model->setExternalConnectionNodes(0, vec_size_internals); + // input connections + line_model->setExternalConnectionNodes(1, vdqbus_index[i]); + line_model->setExternalConnectionNodes(2, vdqbus_index[i] + 1); + // output connections + line_model->setExternalConnectionNodes(3, vdqbus_index[i + 1]); + line_model->setExternalConnectionNodes(4, vdqbus_index[i + 1] + 1); + // internal connections + for (size_t j = 0; j < 2; j++) + { + line_model->setExternalConnectionNodes(5 + j, indexv + j); + } + indexv += 2; + sysmodel->addComponent(line_model); } - indexv += 2; - sysmodel->addComponent(l2); - // line 3 - GridKit::MicrogridLine* l3 = new GridKit::MicrogridLine(6, rline3, Lline3); - // ref motor - l3->setExternalConnectionNodes(0, vec_size_internals); - // input connections - l3->setExternalConnectionNodes(1, dqbus3); - l3->setExternalConnectionNodes(2, dqbus3 + 1); - // output connections - l3->setExternalConnectionNodes(3, dqbus4); - l3->setExternalConnectionNodes(4, dqbus4 + 1); - // internal connections - for (size_t i = 0; i < 2; i++) + // Load all the Load components + for (size_t i = 0; i < Nsize; i++) { - l3->setExternalConnectionNodes(5 + i, indexv + i); + auto* load_model = new MicrogridLoad(model_id++, + rload_list[i], + Lload_list[i]); + // ref motor + load_model->setExternalConnectionNodes(0, vec_size_internals); + // input connections + load_model->setExternalConnectionNodes(1, vdqbus_index[2 * i]); + load_model->setExternalConnectionNodes(2, vdqbus_index[2 * i] + 1); + // internal connections + for (size_t j = 0; j < 2; j++) + { + load_model->setExternalConnectionNodes(3 + j, indexv + j); + } + indexv += 2; + sysmodel->addComponent(load_model); } - indexv += 2; - sysmodel->addComponent(l3); - - // loads - // load 1 - GridKit::MicrogridLoad* load1 = new GridKit::MicrogridLoad(7, rload1, Lload1); - // ref motor - load1->setExternalConnectionNodes(0, vec_size_internals); - // input connections - load1->setExternalConnectionNodes(1, dqbus1); - load1->setExternalConnectionNodes(2, dqbus1 + 1); - // internal connections - for (size_t i = 0; i < 2; i++) + // Add all the microgrid Virtual DQ Buses + for (size_t i = 0; i < 2 * Nsize; i++) { - load1->setExternalConnectionNodes(3 + i, indexv + i); - } - indexv += 2; - sysmodel->addComponent(load1); + auto* virDQbus_model = new MicrogridBusDQ(model_id++, RN); - // load 2 - GridKit::MicrogridLoad* load2 = new GridKit::MicrogridLoad(8, rload2, Lload2); - // ref motor - load2->setExternalConnectionNodes(0, vec_size_internals); - // input connections - load2->setExternalConnectionNodes(1, dqbus3); - load2->setExternalConnectionNodes(2, dqbus3 + 1); - // internal connections - for (size_t i = 0; i < 2; i++) - { - load2->setExternalConnectionNodes(3 + i, indexv + i); + virDQbus_model->setExternalConnectionNodes(0, vdqbus_index[i]); + virDQbus_model->setExternalConnectionNodes(1, vdqbus_index[i] + 1); + sysmodel->addComponent(virDQbus_model); } - indexv += 2; - sysmodel->addComponent(load2); - - // Virtual PQ Buses - GridKit::MicrogridBusDQ* bus1 = new GridKit::MicrogridBusDQ(9, RN); - - bus1->setExternalConnectionNodes(0, dqbus1); - bus1->setExternalConnectionNodes(1, dqbus1 + 1); - sysmodel->addComponent(bus1); - - GridKit::MicrogridBusDQ* bus2 = new GridKit::MicrogridBusDQ(10, RN); - - bus2->setExternalConnectionNodes(0, dqbus2); - bus2->setExternalConnectionNodes(1, dqbus2 + 1); - sysmodel->addComponent(bus2); - - GridKit::MicrogridBusDQ* bus3 = new GridKit::MicrogridBusDQ(11, RN); - - bus3->setExternalConnectionNodes(0, dqbus3); - bus3->setExternalConnectionNodes(1, dqbus3 + 1); - sysmodel->addComponent(bus3); - - GridKit::MicrogridBusDQ* bus4 = new GridKit::MicrogridBusDQ(12, RN); - - bus4->setExternalConnectionNodes(0, dqbus4); - bus4->setExternalConnectionNodes(1, dqbus4 + 1); - sysmodel->addComponent(bus4); sysmodel->allocate(vec_size_total); @@ -294,22 +228,15 @@ namespace GridKit } // Create Intial derivatives specifics generated in MATLAB - // DGs 1 - sysmodel->yp()[2] = parms1.Vn_; - sysmodel->yp()[4] = parms1.Kpv_ * parms1.Vn_; - sysmodel->yp()[6] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_; - sysmodel->yp()[12 + 3] = parms1.Vn_; - sysmodel->yp()[12 + 5] = parms1.Kpv_ * parms1.Vn_; - sysmodel->yp()[12 + 7] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_; - for (size_t i = 2; i < 4; i++) + for (size_t i = 0; i < 2 * Nsize; i++) { - sysmodel->yp()[13 * i - 1 + 3] = parms2.Vn_; - sysmodel->yp()[13 * i - 1 + 5] = parms2.Kpv_ * parms2.Vn_; - sysmodel->yp()[13 * i - 1 + 7] = (parms2.Kpc_ * parms2.Kpv_ * parms2.Vn_) / parms2.Lf_; + sysmodel->yp()[13 * i - 1 + 3] = DGParams_list[i].Vn_; + sysmodel->yp()[13 * i - 1 + 5] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_; + sysmodel->yp()[13 * i - 1 + 7] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_; } - // since the intial P_com = 0 - sysmodel->y()[vec_size_internals] = parms1.wb_; + // since the intial P_com = 0, the set the intial vector to the reference frame + sysmodel->y()[vec_size_internals] = DG_parms1.wb_; sysmodel->y()[vec_size_internals].setVariableNumber(vec_size_internals); for (size_t i = 0; i < vec_size_total; i++) @@ -366,7 +293,8 @@ namespace GridKit delete sysmodel; - return success.report(__func__); + std::string test_name = __func__ + std::format(" (Nsize = {})", Nsize); + return success.report(test_name.c_str()); } }; } // namespace Testing diff --git a/tests/UnitTests/PowerElectronics/runSystemTests.cpp b/tests/UnitTests/PowerElectronics/runSystemTests.cpp index d924e1b8..de40ecfd 100644 --- a/tests/UnitTests/PowerElectronics/runSystemTests.cpp +++ b/tests/UnitTests/PowerElectronics/runSystemTests.cpp @@ -8,7 +8,8 @@ int main() GridKit::Testing::TestingResults result; GridKit::Testing::SystemTests test; - result += test.jacobianSparsity(); + result += test.jacobianSparsity(2); + result += test.jacobianSparsity(128); return result.summary(); } From 7c793594909d34c3c0f811a952245b7926fa1c74 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 13:19:43 -0500 Subject: [PATCH 10/29] Calculate sparsity directly into CsrMatrix Due to the way that assignment operators are currently implemented in CsrMatrix, there is a memory issue. --- .../SystemModelPowerElectronics.hpp | 82 +++++-------------- .../UnitTests/PowerElectronics/CMakeLists.txt | 3 +- .../PowerElectronics/SystemTests.hpp | 4 +- 3 files changed, 25 insertions(+), 64 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index ef61d975..1106f2c4 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -8,6 +8,7 @@ #include #include +#include #include #include #include @@ -83,13 +84,11 @@ namespace GridKit PowerElectronicsModel() { // Set system model parameters as default - rel_tol_ = 1e-4; - abs_tol_ = 1e-4; - max_steps_ = 2000; + rel_tol_ = 1e-4; + abs_tol_ = 1e-4; + max_steps_ = 2000; // By default don't use the jacobian - use_jac_ = false; - jac_sparsity_row_indices_ = nullptr; - jac_sparsity_col_indices_ = nullptr; + use_jac_ = false; } /** @@ -108,13 +107,11 @@ namespace GridKit IdxT max_steps = 2000) { // Set system model tolerances from input - rel_tol_ = rel_tol; - abs_tol_ = abs_tol; - max_steps_ = max_steps; + rel_tol_ = rel_tol; + abs_tol_ = abs_tol; + max_steps_ = max_steps; // Can choose if to use jacobian - use_jac_ = use_jac; - jac_sparsity_row_indices_ = nullptr; - jac_sparsity_col_indices_ = nullptr; + use_jac_ = use_jac; } /** @@ -129,15 +126,6 @@ namespace GridKit { for (auto comp : components_) delete comp; - - if (jac_sparsity_row_indices_ != nullptr) - { - delete[] jac_sparsity_row_indices_; - } - if (jac_sparsity_col_indices_ != nullptr) - { - delete[] jac_sparsity_col_indices_; - } } /** @@ -248,24 +236,18 @@ namespace GridKit } } // Final row index (beyond the end) keeps track of nnz, which is also the size of the column indices. - global_row_indices[size_] = global_col_indices.size(); - - // Delete old sparsity buffers, if they exist - if (jac_sparsity_row_indices_ != nullptr) - { - delete[] jac_sparsity_row_indices_; - } - if (jac_sparsity_col_indices_ != nullptr) - { - delete[] jac_sparsity_col_indices_; - } + global_row_indices[size_] = static_cast(global_col_indices.size()); // Allocate new sparsity buffers - jac_sparsity_row_indices_ = global_row_indices; - jac_sparsity_col_indices_ = new IdxT[global_col_indices.size()]; + IdxT nnz = static_cast(global_col_indices.size()); + auto jac_row_indices = global_row_indices; + auto jac_col_indices = new IdxT[nnz]; + auto jac_values = new RealT[nnz]; // Copy column indices - std::copy(global_col_indices.begin(), global_col_indices.end(), jac_sparsity_col_indices_); + std::copy(global_col_indices.begin(), global_col_indices.end(), jac_col_indices); + + csr_jac_ = LinearAlgebra::CsrMatrix(size_, size_, nnz, &jac_row_indices, &jac_col_indices, &jac_values); delete[] component_sparsities; delete[] reverse_map; @@ -539,40 +521,18 @@ namespace GridKit } /** - * @brief Returns a pointer to the buffer of Jacobian CSR row indices. + * @brief Returns a reference to the CSR Jacobian. * @todo Can be removed once `getJacobian()` returns a `CsrMatrix` */ - IdxT* getJacRowIndices() + LinearAlgebra::CsrMatrix& getCsrJac() { - return jac_sparsity_row_indices_; - } - - /** - * @brief Returns a pointer to the buffer of Jacobian CSR column indices. - * @todo Can be removed once `getJacobian()` returns a `CsrMatrix` - */ - IdxT* getJacColIndices() - { - return jac_sparsity_col_indices_; + return csr_jac_; } private: static constexpr IdxT neg1_ = static_cast(-1); - /** - * @brief Keeps track of the CSR row indices of the system jacobian. - * `nullptr` is used to indicate an un-allocated buffer. - * @todo Unneeded once the jacobian is in CSR format. This can be stored - * in the jacobian itself. - */ - IdxT* jac_sparsity_row_indices_; - /** - * @brief Keeps track of the CSR col indices of the system jacobian. - * `nullptr` is used to indicate an un-allocated buffer. - * @todo Unneeded once the jacobian is in CSR format. This can be stored - * in the jacobian itself. - */ - IdxT* jac_sparsity_col_indices_; + LinearAlgebra::CsrMatrix csr_jac_; std::vector components_; diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt index 556de655..24e1905c 100644 --- a/tests/UnitTests/PowerElectronics/CMakeLists.txt +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -2,7 +2,8 @@ add_executable(test_power_elec_system runSystemTests.cpp) target_link_libraries(test_power_elec_system GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload - GridKit::power_elec_microbusdq) + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) add_test(NAME PowerElectronicsSystemTest COMMAND test_power_elec_system) diff --git a/tests/UnitTests/PowerElectronics/SystemTests.hpp b/tests/UnitTests/PowerElectronics/SystemTests.hpp index 1ea52759..e8a7d628 100644 --- a/tests/UnitTests/PowerElectronics/SystemTests.hpp +++ b/tests/UnitTests/PowerElectronics/SystemTests.hpp @@ -249,8 +249,8 @@ namespace GridKit std::vector& residuals = sysmodel->getResidual(); - size_t* row_indices = sysmodel->getJacRowIndices(); - size_t* col_indices = sysmodel->getJacColIndices(); + size_t* row_indices = sysmodel->getCsrJac().getRowData(); + size_t* col_indices = sysmodel->getCsrJac().getColData(); for (size_t row = 0; row < residuals.size(); row++) { From 56d2f6b68531d078c15b086afb4ca58a85e7fe0a Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 14:01:12 -0500 Subject: [PATCH 11/29] Add levvel of indirection to csr_jac_ --- .../SystemModelPowerElectronics.hpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 1106f2c4..9226bc2c 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -89,6 +89,7 @@ namespace GridKit max_steps_ = 2000; // By default don't use the jacobian use_jac_ = false; + csr_jac_ = nullptr; } /** @@ -112,6 +113,7 @@ namespace GridKit max_steps_ = max_steps; // Can choose if to use jacobian use_jac_ = use_jac; + csr_jac_ = nullptr; } /** @@ -126,6 +128,9 @@ namespace GridKit { for (auto comp : components_) delete comp; + + if (!csr_jac_) + delete csr_jac_; } /** @@ -238,6 +243,9 @@ namespace GridKit // Final row index (beyond the end) keeps track of nnz, which is also the size of the column indices. global_row_indices[size_] = static_cast(global_col_indices.size()); + if (!csr_jac_) + delete csr_jac_; + // Allocate new sparsity buffers IdxT nnz = static_cast(global_col_indices.size()); auto jac_row_indices = global_row_indices; @@ -247,7 +255,7 @@ namespace GridKit // Copy column indices std::copy(global_col_indices.begin(), global_col_indices.end(), jac_col_indices); - csr_jac_ = LinearAlgebra::CsrMatrix(size_, size_, nnz, &jac_row_indices, &jac_col_indices, &jac_values); + csr_jac_ = new LinearAlgebra::CsrMatrix(size_, size_, nnz, &jac_row_indices, &jac_col_indices, &jac_values); delete[] component_sparsities; delete[] reverse_map; @@ -526,13 +534,13 @@ namespace GridKit */ LinearAlgebra::CsrMatrix& getCsrJac() { - return csr_jac_; + return *csr_jac_; } private: static constexpr IdxT neg1_ = static_cast(-1); - LinearAlgebra::CsrMatrix csr_jac_; + LinearAlgebra::CsrMatrix* csr_jac_; std::vector components_; From 0266df9b7c930a81a3a7f0ca078089b4c1cd88b5 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 14:53:36 -0500 Subject: [PATCH 12/29] Add resize() to CsrMatrix --- .../LinearAlgebra/SparseMatrix/CsrMatrix.cpp | 19 +++++++++++++++++++ .../LinearAlgebra/SparseMatrix/CsrMatrix.hpp | 2 ++ .../SystemModelPowerElectronics.hpp | 15 ++++----------- 3 files changed, 25 insertions(+), 11 deletions(-) diff --git a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp index a1296c57..7a3331b5 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp @@ -818,6 +818,25 @@ namespace GridKit } } + /** + * @brief Changes the size of the matrix. De-allocates the matrix, if needed. + * Does not re-allocate any memory. + * @param n The new number of rows + * @param m The new number of columns + * @return 0 if succesful + */ + template + int CsrMatrix::resize(IdxT n, IdxT m) + { + destroyMatrixData(memory::HOST); + destroyMatrixData(memory::DEVICE); + + n_ = n; + m_ = m; + + return 0; + } + template class CsrMatrix; template class CsrMatrix; } // namespace LinearAlgebra diff --git a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp index 975db353..cd7e18ff 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp @@ -76,6 +76,8 @@ namespace GridKit virtual int setValuesPointer(RealT* new_vals, memory::MemorySpace memspace); + int resize(IdxT n, IdxT m); + private: IdxT n_{0}; ///< number of rows IdxT m_{0}; ///< number of columns diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 9226bc2c..a51db5db 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -89,7 +89,6 @@ namespace GridKit max_steps_ = 2000; // By default don't use the jacobian use_jac_ = false; - csr_jac_ = nullptr; } /** @@ -113,7 +112,6 @@ namespace GridKit max_steps_ = max_steps; // Can choose if to use jacobian use_jac_ = use_jac; - csr_jac_ = nullptr; } /** @@ -128,9 +126,6 @@ namespace GridKit { for (auto comp : components_) delete comp; - - if (!csr_jac_) - delete csr_jac_; } /** @@ -243,9 +238,6 @@ namespace GridKit // Final row index (beyond the end) keeps track of nnz, which is also the size of the column indices. global_row_indices[size_] = static_cast(global_col_indices.size()); - if (!csr_jac_) - delete csr_jac_; - // Allocate new sparsity buffers IdxT nnz = static_cast(global_col_indices.size()); auto jac_row_indices = global_row_indices; @@ -255,7 +247,8 @@ namespace GridKit // Copy column indices std::copy(global_col_indices.begin(), global_col_indices.end(), jac_col_indices); - csr_jac_ = new LinearAlgebra::CsrMatrix(size_, size_, nnz, &jac_row_indices, &jac_col_indices, &jac_values); + csr_jac_.resize(size_, size_); + csr_jac_.setDataPointers(jac_row_indices, jac_col_indices, jac_values, LinearAlgebra::memory::HOST); delete[] component_sparsities; delete[] reverse_map; @@ -534,13 +527,13 @@ namespace GridKit */ LinearAlgebra::CsrMatrix& getCsrJac() { - return *csr_jac_; + return csr_jac_; } private: static constexpr IdxT neg1_ = static_cast(-1); - LinearAlgebra::CsrMatrix* csr_jac_; + LinearAlgebra::CsrMatrix csr_jac_; std::vector components_; From 6a999c0b1481992c69e54080eebc3718b0bdecd0 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 15:57:09 -0500 Subject: [PATCH 13/29] Add CsrMatrix dependency to power electronics examples --- .../DistributedGeneratorTest/CMakeLists.txt | 3 ++- examples/PowerElectronics/Microgrid/CMakeLists.txt | 3 ++- examples/PowerElectronics/RLCircuit/CMakeLists.txt | 3 ++- examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt | 6 ++++-- 4 files changed, 10 insertions(+), 5 deletions(-) diff --git a/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt b/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt index 225dd73f..d17aa846 100644 --- a/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt +++ b/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt @@ -5,7 +5,8 @@ add_executable(dgtest DGTest.cpp) target_link_libraries(dgtest GridKit::power_elec_disgen GridKit::power_elec_microline - GridKit::power_elec_microload) + GridKit::power_elec_microload + GridKit::CsrMatrix) add_test(NAME DistributedGeneratorTest COMMAND $) install(TARGETS dgtest RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/Microgrid/CMakeLists.txt b/examples/PowerElectronics/Microgrid/CMakeLists.txt index 1c67f74f..ef82d2e4 100644 --- a/examples/PowerElectronics/Microgrid/CMakeLists.txt +++ b/examples/PowerElectronics/Microgrid/CMakeLists.txt @@ -7,7 +7,8 @@ target_link_libraries(microgrid GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload GridKit::solvers_dyn - GridKit::power_elec_microbusdq) + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) add_test(NAME Microgrid COMMAND $) install(TARGETS microgrid RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/RLCircuit/CMakeLists.txt b/examples/PowerElectronics/RLCircuit/CMakeLists.txt index 42ffce87..52af279f 100644 --- a/examples/PowerElectronics/RLCircuit/CMakeLists.txt +++ b/examples/PowerElectronics/RLCircuit/CMakeLists.txt @@ -7,7 +7,8 @@ target_link_libraries(rlcircuit GridKit::power_elec_capacitor GridKit::power_elec_inductor GridKit::power_elec_resistor GridKit::power_elec_voltagesource - GridKit::solvers_dyn) + GridKit::solvers_dyn + GridKit::CsrMatrix) add_test(NAME RLCircuit COMMAND $) install(TARGETS rlcircuit RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt b/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt index aeea5f01..c711d0d6 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt +++ b/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt @@ -8,13 +8,15 @@ target_link_libraries(scalemicrogrid GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload GridKit::solvers_dyn - GridKit::power_elec_microbusdq) + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) target_link_libraries(scalemicrogridarbitrary GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload GridKit::solvers_dyn - GridKit::power_elec_microbusdq) + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) add_test(NAME ScaleMicrogrid COMMAND $) install(TARGETS scalemicrogrid RUNTIME DESTINATION bin) From f6fcd97789f324bcb0338672023d97d87c91f39d Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 16:28:37 -0500 Subject: [PATCH 14/29] Allocate csr matrix properly --- .../SystemModelPowerElectronics.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index a51db5db..515ed638 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -239,19 +239,19 @@ namespace GridKit global_row_indices[size_] = static_cast(global_col_indices.size()); // Allocate new sparsity buffers - IdxT nnz = static_cast(global_col_indices.size()); - auto jac_row_indices = global_row_indices; - auto jac_col_indices = new IdxT[nnz]; - auto jac_values = new RealT[nnz]; - - // Copy column indices - std::copy(global_col_indices.begin(), global_col_indices.end(), jac_col_indices); + IdxT nnz = static_cast(global_col_indices.size()); csr_jac_.resize(size_, size_); - csr_jac_.setDataPointers(jac_row_indices, jac_col_indices, jac_values, LinearAlgebra::memory::HOST); + csr_jac_.setNnz(nnz); + csr_jac_.allocateMatrixData(LinearAlgebra::memory::HOST); + + // Copy column indices + std::copy(global_col_indices.begin(), global_col_indices.end(), csr_jac_.getColData()); + std::copy(global_row_indices, global_row_indices + size_ + 1, csr_jac_.getRowData()); delete[] component_sparsities; delete[] reverse_map; + delete[] global_row_indices; return 1; } From 2522abc32acecc6ef0369958150c6dfe2482303f Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 13 Nov 2025 16:30:19 -0500 Subject: [PATCH 15/29] [skip ci] Update documentation --- GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp index 7a3331b5..a6291f2b 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp @@ -819,7 +819,8 @@ namespace GridKit } /** - * @brief Changes the size of the matrix. De-allocates the matrix, if needed. + * @brief Changes the size of the matrix. De-allocates the matrix, if needed, since the + * old indices corresponded to a differently-size matrix and no longer make sense. * Does not re-allocate any memory. * @param n The new number of rows * @param m The new number of columns From 22edecbb93a206e3abaaef362141b8d54f9b08f8 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Wed, 19 Nov 2025 13:04:05 -0500 Subject: [PATCH 16/29] Ensure system vector is mapped correctly --- .../SystemModelPowerElectronics.hpp | 43 ++++++++++++++++++- 1 file changed, 42 insertions(+), 1 deletion(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 515ed638..feabcde3 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -64,6 +64,8 @@ namespace GridKit using component_type = CircuitComponent; using CircuitComponent::size_; + using CircuitComponent::n_intern_; + using CircuitComponent::n_extern_; using CircuitComponent::nnz_; using CircuitComponent::time_; using CircuitComponent::alpha_; @@ -291,10 +293,49 @@ namespace GridKit int allocate(IdxT s) { // Allocate all components - size_ = s; + size_ = s; + n_intern_ = 0; for (const auto& component : components_) { component->allocate(); + + // Count up the amount of internal variables which get mapped to system variables + auto extern_indices = component->getExternIndices(); + for (IdxT comp_var_idx = 0; comp_var_idx < component->size(); comp_var_idx++) + { + IdxT sys_var_idx = component->getNodeConnection(comp_var_idx); + if (!extern_indices.contains(comp_var_idx) && sys_var_idx != neg1_) + { + n_intern_++; + } + } + } + n_extern_ = size_ - n_intern_; + + // Ensure that all component locals are mapped to system locals + // and all component externals are mapped to system externals. + // System locals are stored first, in 0..n_intern_ and externals + // are stored after, in n_intern_.. + for (auto component : components_) + { + auto extern_indices = component->getExternIndices(); + for (IdxT comp_var_idx = 0; comp_var_idx < component->size(); comp_var_idx++) + { + IdxT sys_var_idx = component->getNodeConnection(comp_var_idx); + + // Ignore local variables which aren't mapped to sytem variables + if (sys_var_idx == neg1_) + continue; + + if (extern_indices.contains(comp_var_idx)) + { + assert(sys_var_idx >= n_intern_); + } + else + { + assert(sys_var_idx < n_intern_); + } + } } // Allocate global vectors From 756d9a953419730e567171f62206fe4613f62173 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Wed, 19 Nov 2025 13:07:37 -0500 Subject: [PATCH 17/29] Ensure local variable grouping in system vector --- .../SystemModelPowerElectronics.hpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index feabcde3..20383ab5 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -319,6 +319,12 @@ namespace GridKit for (auto component : components_) { auto extern_indices = component->getExternIndices(); + + // Whether or not we've seen a local variable yet + bool has_seen_local = false; + // The system index of the last local we've seen + IdxT last_local_sys; + for (IdxT comp_var_idx = 0; comp_var_idx < component->size(); comp_var_idx++) { IdxT sys_var_idx = component->getNodeConnection(comp_var_idx); @@ -334,6 +340,16 @@ namespace GridKit else { assert(sys_var_idx < n_intern_); + + // Make sure that all of the locals for a particular component are grouped + // together in the sytem vector, and have increasing indices. + if (has_seen_local) + { + assert(sys_var_idx == last_local_sys + 1); + } + + has_seen_local = true; + last_local_sys = sys_var_idx; } } } From a68ec9d1d4cd87b4937f4c4a3105b3c4fc7db83c Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Wed, 19 Nov 2025 13:30:57 -0500 Subject: [PATCH 18/29] Ensure components are sorted in system vector --- .../PowerElectronics/SystemModelPowerElectronics.hpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 20383ab5..2d0e21d2 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -316,8 +316,9 @@ namespace GridKit // and all component externals are mapped to system externals. // System locals are stored first, in 0..n_intern_ and externals // are stored after, in n_intern_.. - for (auto component : components_) + for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) { + auto component = components_[comp_idx]; auto extern_indices = component->getExternIndices(); // Whether or not we've seen a local variable yet @@ -347,6 +348,11 @@ namespace GridKit { assert(sys_var_idx == last_local_sys + 1); } + else if (comp_idx > 0) + { + // Ensure increasing indices in components - so no need to sort components + assert(sys_var_idx > last_local_sys); + } has_seen_local = true; last_local_sys = sys_var_idx; From c031309e2452378785976f4687bd9b95de60c366 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Wed, 19 Nov 2025 13:32:12 -0500 Subject: [PATCH 19/29] Comment --- GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 2d0e21d2..f800c90f 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -316,6 +316,9 @@ namespace GridKit // and all component externals are mapped to system externals. // System locals are stored first, in 0..n_intern_ and externals // are stored after, in n_intern_.. + // Additionally, ensure that components locals are grouped in the system vectors - no other + // variables are between locals from a single component. As well, ensure that components are + // sorted by these groupings, so the first component is the first block and so on. for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) { auto component = components_[comp_idx]; From 1dbdc34e9f0ee3addd26951b0f7bf432b8f87ffa Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 20 Nov 2025 11:01:29 -0500 Subject: [PATCH 20/29] Construct internal and external row sparsities separately Also removed vector for global_col_indices, instead allocating once with an upper bound of the combined nnz of the components. This shouldn't be *too* much of an over-estimation. --- .../SystemModelPowerElectronics.hpp | 118 ++++++++++++------ 1 file changed, 80 insertions(+), 38 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index f800c90f..9c40cdd1 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -157,9 +157,10 @@ namespace GridKit std::vector col_indices_; }; - // A reverse mapping from system variables -> component variables - auto reverse_map = new std::forward_list[size_]; - auto component_sparsities = new CsrSparsity[components_.size()]; + // A reverse mapping from external system variables -> component variables + auto reverse_extern_map = new std::forward_list[n_extern_]; + auto component_sparsities = new CsrSparsity[components_.size()]; + size_t component_nnz = 0; // Loop over all components, evaluate their jacobians, save their sparsity information, // and construct the reverse variable mapping. @@ -170,46 +171,90 @@ namespace GridKit auto [row_indices, col_indices, _] = component->getJacobian().setDataToCSR(); - component_sparsities[comp_idx] = CsrSparsity{ - .row_indices_ = std::move(row_indices), - .col_indices_ = std::move(col_indices), + auto& comp_sparsity = component_sparsities[comp_idx]; + comp_sparsity = CsrSparsity{ + .row_indices_ = std::move(row_indices), + .col_indices_ = std::move(col_indices), }; - for (IdxT local_row = 0; local_row < component->size(); local_row++) + for (IdxT local_external_row : component->getExternIndices()) { - IdxT global_row = component->getNodeConnection(local_row); + IdxT global_row = component->getNodeConnection(local_external_row); // Not all local variables map to a global variable if (global_row != neg1_) { - reverse_map[global_row].push_front({.comp_idx_ = comp_idx, .local_row_idx_ = local_row}); + // global_row >= n_intern_ (see allocate(s)) + reverse_extern_map[global_row - n_intern_].push_front({.comp_idx_ = comp_idx, .local_row_idx_ = local_external_row}); } } + + component_nnz += comp_sparsity.row_indices_.back(); } // Allocate the final sparsity pattern info - IdxT* global_row_indices = new IdxT[size_ + 1]; - std::vector global_col_indices; - global_col_indices.reserve(size_); - - // Construct sparsity pattern row-by-row. In the future, can be batched - // into contiguous blocks of rows which can be calculated in parallel, - // then joined sequentially - for (size_t row = 0; row < size_; row++) + IdxT* global_row_indices = new IdxT[size_ + 1]; + IdxT* global_col_indices = new IdxT[component_nnz]; // Use component_nnz as an upper bound on nnz + global_row_indices[0] = 0; + + // Construct Jacobian sparsity pattern + + // Start with internal rows, which must be before external rows, are grouped by component, + // strictly increasing wrt component internals, and must be sorted by component (see allocate(s)) + size_t curr_internal_row = 0; + for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) { - global_row_indices[row] = global_col_indices.size(); + auto component = components_[comp_idx]; + auto comp_externals = component->getExternIndices(); + auto& comp_sparsity = component_sparsities[comp_idx]; + + for (IdxT local_row = 0; local_row < component->size(); local_row++) + { + // Skip external variables and variables which aren't system variables + if (comp_externals.contains(local_row) || component->getNodeConnection(local_row) == neg1_) + continue; + + global_row_indices[curr_internal_row + 1] = global_row_indices[curr_internal_row]; + + assert(component->getNodeConnection(local_row) == curr_internal_row); + + IdxT row_start = comp_sparsity.row_indices_[local_row]; + IdxT row_end = comp_sparsity.row_indices_[local_row + 1]; + for (IdxT local_col_idx = row_start; local_col_idx < row_end; local_col_idx++) + { + IdxT global_col_idx = component->getNodeConnection(comp_sparsity.col_indices_[local_col_idx]); - // How many components contribute to this row - size_t num_contrib = 0; + if (global_col_idx == neg1_) + continue; + + global_col_indices[global_row_indices[curr_internal_row + 1]] = global_col_idx; + global_row_indices[curr_internal_row + 1]++; + } + + // Need to sort internal rows because even though the mapping from component internal to system internal + // is monotonic, the mapping from component external to system external may not be, and internal rows + // may contain external columns. + auto global_row_start = global_col_indices + global_row_indices[curr_internal_row]; + auto global_row_end = global_col_indices + global_row_indices[curr_internal_row + 1]; + std::sort(global_row_start, global_row_end); + + curr_internal_row++; + } + } + + assert(curr_internal_row == n_intern_); + + // Then move on to external rows, which come after internal rows + for (size_t row = n_intern_; row < size_; row++) + { + global_row_indices[row + 1] = global_row_indices[row]; // Collect columns from each component which has a row which contributes to this row - for (ComponentContribution contrib : reverse_map[row]) + for (ComponentContribution contrib : reverse_extern_map[row - n_intern_]) { component_type* component = components_[contrib.comp_idx_]; CsrSparsity& comp_sparsity = component_sparsities[contrib.comp_idx_]; - num_contrib++; - IdxT row_start = comp_sparsity.row_indices_[contrib.local_row_idx_]; IdxT row_end = comp_sparsity.row_indices_[contrib.local_row_idx_ + 1]; for (size_t local_elem_idx = row_start; local_elem_idx < row_end; local_elem_idx++) @@ -220,40 +265,37 @@ namespace GridKit // Not all columns map to columns in the system jacobian if (global_col != neg1_) { - global_col_indices.push_back(global_col); + global_col_indices[global_row_indices[row + 1]] = global_col; + global_row_indices[row + 1]++; } } } // Sort the row by column indices. Since the mapping from local indices to global indices isn't monotonically increasing, // this is necessary. - auto start = std::next(global_col_indices.begin(), global_row_indices[row]); - std::sort(start, global_col_indices.end()); + auto start = std::next(global_col_indices, global_row_indices[row]); + auto end = std::next(global_col_indices, global_row_indices[row + 1]); + std::sort(start, end); - // If there were multiple components which contributed columns, then we definitely need - // to de-duplicate the columns. - if (num_contrib > 1) - { - global_col_indices.erase(std::unique(start, global_col_indices.end()), global_col_indices.end()); - } + // De-duplicate the columns + auto new_end = std::unique(start, end); + global_row_indices[row + 1] = global_row_indices[row] + std::distance(start, new_end); } - // Final row index (beyond the end) keeps track of nnz, which is also the size of the column indices. - global_row_indices[size_] = static_cast(global_col_indices.size()); - // Allocate new sparsity buffers - IdxT nnz = static_cast(global_col_indices.size()); + IdxT nnz = global_row_indices[size_]; csr_jac_.resize(size_, size_); csr_jac_.setNnz(nnz); csr_jac_.allocateMatrixData(LinearAlgebra::memory::HOST); // Copy column indices - std::copy(global_col_indices.begin(), global_col_indices.end(), csr_jac_.getColData()); + std::copy(global_col_indices, global_col_indices + nnz, csr_jac_.getColData()); std::copy(global_row_indices, global_row_indices + size_ + 1, csr_jac_.getRowData()); delete[] component_sparsities; - delete[] reverse_map; + delete[] reverse_extern_map; delete[] global_row_indices; + delete[] global_col_indices; return 1; } From ba0794ee51d98436538b004bc1258dd653482322 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 20 Nov 2025 11:35:57 -0500 Subject: [PATCH 21/29] Removed some magic numbers --- .../DistributedGenerator.cpp | 6 +++--- .../DistributedGenerator.hpp | 3 +++ .../MicrogridBusDQ/MicrogridBusDQ.cpp | 6 +++--- .../MicrogridBusDQ/MicrogridBusDQ.hpp | 3 +++ .../MicrogridLine/MicrogridLine.cpp | 6 +++--- .../MicrogridLine/MicrogridLine.hpp | 3 +++ .../MicrogridLoad/MicrogridLoad.cpp | 6 +++--- .../MicrogridLoad/MicrogridLoad.hpp | 3 +++ .../PowerElectronics/SystemTests.hpp | 19 +++++++++++++++---- 9 files changed, 39 insertions(+), 16 deletions(-) diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index 5e1eee4d..88ed090e 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -39,9 +39,9 @@ namespace GridKit { // internals [\delta_i, Pi, Qi, phi_di, phi_qi, gamma_di, gamma_qi, il_di, il_qi, vo_di, vo_qi, io_di, io_qi] // externals [\omega_ref, vba_out, vbb_out] - size_ = 16; - n_intern_ = 13; - n_extern_ = 3; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1, 2}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp index d2c930ac..8b54c9aa 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -64,6 +64,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 13; + static constexpr size_t NUM_EXTERNALS = 3; + DistributedGenerator(IdxT id, DistributedGeneratorParameters parm, bool reference_frame); diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp index fd5dbe8f..2f159448 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -24,9 +24,9 @@ namespace GridKit : RN_(RN) { // externals [vbus_d, vbus_q] - size_ = 2; - n_intern_ = 0; - n_extern_ = 2; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index 8c0d8dad..646e56fd 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -43,6 +43,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 0; + static constexpr size_t NUM_EXTERNALS = 2; + MicrogridBusDQ(IdxT id, RealT RN); virtual ~MicrogridBusDQ(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index 44f9073d..b5a8f6b8 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -28,9 +28,9 @@ namespace GridKit { // internals [id, iq] // externals [\omegaref, vbd_in, vbq_in, vbd_out, vbq_out] - size_ = 7; - n_intern_ = 2; - n_extern_ = 5; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1, 2, 3, 4}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp index 64992c5c..1f646251 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -44,6 +44,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 2; + static constexpr size_t NUM_EXTERNALS = 5; + MicrogridLine(IdxT id, RealT R, RealT L); virtual ~MicrogridLine(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index e16b70b6..2ca50f90 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -26,9 +26,9 @@ namespace GridKit { // internals [id, iq] // externals [\omegaref, vbd_out, vbq_out] - size_ = 5; - n_intern_ = 2; - n_extern_ = 3; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1, 2}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 0a8f34bb..b8e0432c 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -44,6 +44,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 2; + static constexpr size_t NUM_EXTERNALS = 3; + MicrogridLoad(IdxT id, RealT R, RealT L); virtual ~MicrogridLoad(); diff --git a/tests/UnitTests/PowerElectronics/SystemTests.hpp b/tests/UnitTests/PowerElectronics/SystemTests.hpp index e8a7d628..9cfef285 100644 --- a/tests/UnitTests/PowerElectronics/SystemTests.hpp +++ b/tests/UnitTests/PowerElectronics/SystemTests.hpp @@ -106,10 +106,21 @@ namespace GridKit rload_list[0] = rload1; Lload_list[0] = Lload1; - // DGs + - refframe Lines + Loads - size_t vec_size_internals = 13 * (2 * Nsize) - 1 + (2 + 4 * (Nsize - 1)) + 2 * Nsize; - // \omegaref + BusDQ - size_t vec_size_externals = 1 + 2 * (2 * Nsize); + size_t num_dgs = 2 * Nsize; + size_t num_lines = 2 * Nsize - 1; + size_t num_loads = Nsize; + size_t num_buses = 2 * Nsize; + + size_t vec_size_internals = + // - refframe + num_dgs * DistributedGenerator::NUM_INTERNALS - 1 + + num_lines * MicrogridLine::NUM_INTERNALS + + num_loads * MicrogridLoad::NUM_INTERNALS + + num_buses * MicrogridBusDQ::NUM_INTERNALS; + + size_t vec_size_externals = + // + refframe + num_buses * MicrogridBusDQ::NUM_EXTERNALS + 1; std::vector vdqbus_index(2 * Nsize, 0); vdqbus_index[0] = vec_size_internals + 1; From a371239256a269249a730903eaf01046daceea9b Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 20 Nov 2025 11:44:47 -0500 Subject: [PATCH 22/29] Fix RLCircuit example's system vector ordering --- .../PowerElectronics/RLCircuit/RLCircuit.cpp | 48 ++++++++++++------- 1 file changed, 30 insertions(+), 18 deletions(-) diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index dce13668..f3bac178 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -31,15 +31,27 @@ int main(int /* argc */, char const** /* argv */) double linit = 1.0; double vinit = 1.0; + // Current system map + // 0 -> Resistor/VSource external + // 1 -> Inductor/Resistor external + // 2 -> Inductor internal + // 3 -> VSource internal + + // New system map + // 0 -> Inductor internal (old 2) + // 1 -> VSource internal (old 3) + // 2 -> Resistor/VSource external (old 0) + // 3 -> Inductor/Resistor external (old 1) + // inductor GridKit::Inductor* induct = new GridKit::Inductor(idoff, linit); // Form index to node uid realations // input - induct->setExternalConnectionNodes(0, 1); + induct->setExternalConnectionNodes(0, 3); // output induct->setExternalConnectionNodes(1, static_cast(-1)); // internal - induct->setExternalConnectionNodes(2, 2); + induct->setExternalConnectionNodes(2, 0); // add component sysmodel.addComponent(induct); @@ -48,9 +60,9 @@ int main(int /* argc */, char const** /* argv */) GridKit::Resistor* resis = new GridKit::Resistor(idoff, rinit); // Form index to node uid realations // input - resis->setExternalConnectionNodes(0, 0); + resis->setExternalConnectionNodes(0, 2); // output - resis->setExternalConnectionNodes(1, 1); + resis->setExternalConnectionNodes(1, 3); // add sysmodel.addComponent(resis); @@ -61,9 +73,9 @@ int main(int /* argc */, char const** /* argv */) // input vsource->setExternalConnectionNodes(0, static_cast(-1)); // output - vsource->setExternalConnectionNodes(1, 0); + vsource->setExternalConnectionNodes(1, 2); // internal - vsource->setExternalConnectionNodes(2, 3); + vsource->setExternalConnectionNodes(2, 1); sysmodel.addComponent(vsource); @@ -74,15 +86,15 @@ int main(int /* argc */, char const** /* argv */) // Grounding for IDA. If no grounding then circuit is \mu > 1 // v_0 (grounded) // Create Intial points - sysmodel.y()[0] = vinit; // v_1 - sysmodel.y()[1] = vinit; // v_2 - sysmodel.y()[2] = 0.0; // i_L - sysmodel.y()[3] = 0.0; // i_s + sysmodel.y()[0] = 0.0; // i_L + sysmodel.y()[1] = 0.0; // i_s + sysmodel.y()[2] = vinit; // v_1 + sysmodel.y()[3] = vinit; // v_2 - sysmodel.yp()[0] = 0.0; // v'_1 - sysmodel.yp()[1] = 0.0; // v'_2 - sysmodel.yp()[2] = -vinit / linit; // i'_s - sysmodel.yp()[3] = -vinit / linit; // i'_L + sysmodel.yp()[0] = -vinit / linit; // i'_s + sysmodel.yp()[1] = -vinit / linit; // i'_L + sysmodel.yp()[2] = 0.0; // v'_1 + sysmodel.yp()[3] = 0.0; // v'_2 sysmodel.initialize(); sysmodel.evaluateResidual(); @@ -123,10 +135,10 @@ int main(int /* argc */, char const** /* argv */) std::vector yexact(4); // analytical solution to the circuit - yexact[0] = vinit; - yexact[2] = (vinit / rinit) * (exp(-(rinit / linit) * t_final) - 1.0); - yexact[3] = yexact[2]; - yexact[1] = vinit + rinit * yexact[2]; + yexact[2] = vinit; + yexact[0] = (vinit / rinit) * (exp(-(rinit / linit) * t_final) - 1.0); + yexact[1] = yexact[0]; + yexact[3] = vinit + rinit * yexact[0]; std::cout << "Element-wise Relative error at t=" << t_final << "\n"; for (size_t i = 0; i < yfinial.size(); i++) From e2b8ba228e277321a6e94a626b6851ccedb80874 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Thu, 20 Nov 2025 12:01:32 -0500 Subject: [PATCH 23/29] [skip ci] Remove extra comment --- .../PowerElectronics/RLCircuit/RLCircuit.cpp | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index f3bac178..733aabbe 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -31,17 +31,11 @@ int main(int /* argc */, char const** /* argv */) double linit = 1.0; double vinit = 1.0; - // Current system map - // 0 -> Resistor/VSource external - // 1 -> Inductor/Resistor external - // 2 -> Inductor internal - // 3 -> VSource internal - - // New system map - // 0 -> Inductor internal (old 2) - // 1 -> VSource internal (old 3) - // 2 -> Resistor/VSource external (old 0) - // 3 -> Inductor/Resistor external (old 1) + // System vector map + // 0 -> Inductor internal + // 1 -> VSource internal + // 2 -> Resistor/VSource external + // 3 -> Inductor/Resistor external // inductor GridKit::Inductor* induct = new GridKit::Inductor(idoff, linit); From e87b59be3635fe843e9134ce2e4cb411d61e522d Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Fri, 12 Dec 2025 12:50:37 -0500 Subject: [PATCH 24/29] PowerElectronics test 13-->dg_ref.NUM_INTERNALS. --- tests/UnitTests/PowerElectronics/SystemTests.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/UnitTests/PowerElectronics/SystemTests.hpp b/tests/UnitTests/PowerElectronics/SystemTests.hpp index 9cfef285..44780a46 100644 --- a/tests/UnitTests/PowerElectronics/SystemTests.hpp +++ b/tests/UnitTests/PowerElectronics/SystemTests.hpp @@ -144,14 +144,14 @@ namespace GridKit //"grounding" of the difference dg_ref->setExternalConnectionNodes(3, static_cast(-1)); // internal connections - for (size_t i = 0; i < 12; i++) + for (size_t i = 0; i < dg_ref->NUM_INTERNALS - 1; i++) { dg_ref->setExternalConnectionNodes(4 + i, i); } sysmodel->addComponent(dg_ref); // Keep track of models and index location - size_t indexv = 12; + size_t indexv = dg_ref->NUM_INTERNALS - 1; size_t model_id = 1; // Add all other DGs for (size_t i = 1; i < 2 * Nsize; i++) @@ -166,11 +166,11 @@ namespace GridKit dg->setExternalConnectionNodes(1, vdqbus_index[i]); dg->setExternalConnectionNodes(2, vdqbus_index[i] + 1); // internal connections - for (size_t j = 0; j < 13; j++) + for (size_t j = 0; j < dg_ref->NUM_INTERNALS; j++) { dg->setExternalConnectionNodes(3 + j, indexv + j); } - indexv += 13; + indexv += dg_ref->NUM_INTERNALS; sysmodel->addComponent(dg); } @@ -241,9 +241,9 @@ namespace GridKit // Create Intial derivatives specifics generated in MATLAB for (size_t i = 0; i < 2 * Nsize; i++) { - sysmodel->yp()[13 * i - 1 + 3] = DGParams_list[i].Vn_; - sysmodel->yp()[13 * i - 1 + 5] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_; - sysmodel->yp()[13 * i - 1 + 7] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_; + sysmodel->yp()[dg_ref->NUM_INTERNALS * i - 1 + 3] = DGParams_list[i].Vn_; + sysmodel->yp()[dg_ref->NUM_INTERNALS * i - 1 + 5] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_; + sysmodel->yp()[dg_ref->NUM_INTERNALS * i - 1 + 7] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_; } // since the intial P_com = 0, the set the intial vector to the reference frame From 7975429f3edbdf34961b8042d7d9ef8f4cd1490e Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Fri, 12 Dec 2025 14:13:15 -0500 Subject: [PATCH 25/29] Workaround warnings. --- .../Model/PowerElectronics/SystemModelPowerElectronics.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 9c40cdd1..3449a14c 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -273,13 +273,13 @@ namespace GridKit // Sort the row by column indices. Since the mapping from local indices to global indices isn't monotonically increasing, // this is necessary. - auto start = std::next(global_col_indices, global_row_indices[row]); - auto end = std::next(global_col_indices, global_row_indices[row + 1]); + auto start = std::next(global_col_indices, static_cast(global_row_indices[row])); + auto end = std::next(global_col_indices, static_cast(global_row_indices[row + 1])); std::sort(start, end); // De-duplicate the columns auto new_end = std::unique(start, end); - global_row_indices[row + 1] = global_row_indices[row] + std::distance(start, new_end); + global_row_indices[row + 1] = global_row_indices[row] + static_cast(std::distance(start, new_end)); } // Allocate new sparsity buffers IdxT nnz = global_row_indices[size_]; From 295b5721da0667bc721efcb6cdbd3150940a14eb Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 22 Dec 2025 13:46:29 -0500 Subject: [PATCH 26/29] Remove std::next --- .../Model/PowerElectronics/SystemModelPowerElectronics.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 3449a14c..097355ae 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -273,8 +273,8 @@ namespace GridKit // Sort the row by column indices. Since the mapping from local indices to global indices isn't monotonically increasing, // this is necessary. - auto start = std::next(global_col_indices, static_cast(global_row_indices[row])); - auto end = std::next(global_col_indices, static_cast(global_row_indices[row + 1])); + auto start = global_col_indices + global_row_indices[row]; + auto end = global_col_indices + global_row_indices[row + 1]; std::sort(start, end); // De-duplicate the columns From 00360c05603426cfd8b5e131ba124bca0947bd00 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 22 Dec 2025 13:47:56 -0500 Subject: [PATCH 27/29] Remove component CSR sparsities from alloc sparsity calculations --- .../SystemModelPowerElectronics.hpp | 110 +++++++++--------- 1 file changed, 58 insertions(+), 52 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 097355ae..1e2377ed 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -2,6 +2,7 @@ #pragma once +#include #include #include #include @@ -149,18 +150,9 @@ namespace GridKit size_t local_row_idx_; }; - // Helper struct for saving CSR sparsities of a particular component, since - // right now they only report COO. Can be removed once components have CSR matrices. - struct CsrSparsity - { - std::vector row_indices_; - std::vector col_indices_; - }; - // A reverse mapping from external system variables -> component variables - auto reverse_extern_map = new std::forward_list[n_extern_]; - auto component_sparsities = new CsrSparsity[components_.size()]; - size_t component_nnz = 0; + auto reverse_extern_map = new std::forward_list[n_extern_]; + size_t component_nnz = 0; // Loop over all components, evaluate their jacobians, save their sparsity information, // and construct the reverse variable mapping. @@ -169,13 +161,7 @@ namespace GridKit component_type* component = components_[comp_idx]; component->evaluateJacobian(); - auto [row_indices, col_indices, _] = component->getJacobian().setDataToCSR(); - - auto& comp_sparsity = component_sparsities[comp_idx]; - comp_sparsity = CsrSparsity{ - .row_indices_ = std::move(row_indices), - .col_indices_ = std::move(col_indices), - }; + auto& comp_jacobian = component->getJacobian(); for (IdxT local_external_row : component->getExternIndices()) { @@ -189,7 +175,7 @@ namespace GridKit } } - component_nnz += comp_sparsity.row_indices_.back(); + component_nnz += comp_jacobian.nnz(); } // Allocate the final sparsity pattern info @@ -204,44 +190,56 @@ namespace GridKit size_t curr_internal_row = 0; for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) { - auto component = components_[comp_idx]; - auto comp_externals = component->getExternIndices(); - auto& comp_sparsity = component_sparsities[comp_idx]; + component_type* component = components_[comp_idx]; + auto comp_externals = component->getExternIndices(); + const auto& [rows, columns, values] = component->getJacobian().getEntries(); - for (IdxT local_row = 0; local_row < component->size(); local_row++) - { - // Skip external variables and variables which aren't system variables - if (comp_externals.contains(local_row) || component->getNodeConnection(local_row) == neg1_) - continue; + if (rows.empty()) + continue; - global_row_indices[curr_internal_row + 1] = global_row_indices[curr_internal_row]; + // Loop through all elements of the component jacobian, adding them to the system jacobian if needed + for (size_t elem_idx = 0; elem_idx < rows.size(); elem_idx++) + { + IdxT local_row = rows[elem_idx]; + IdxT global_row = component->getNodeConnection(local_row); - assert(component->getNodeConnection(local_row) == curr_internal_row); + // Skip variables which aren't system variables + if (comp_externals.contains(local_row) || global_row == neg1_) + continue; - IdxT row_start = comp_sparsity.row_indices_[local_row]; - IdxT row_end = comp_sparsity.row_indices_[local_row + 1]; - for (IdxT local_col_idx = row_start; local_col_idx < row_end; local_col_idx++) + if (global_row > curr_internal_row) { - IdxT global_col_idx = component->getNodeConnection(comp_sparsity.col_indices_[local_col_idx]); + curr_internal_row++; + global_row_indices[curr_internal_row + 1] = global_row_indices[curr_internal_row]; + } - if (global_col_idx == neg1_) - continue; + assert(global_row == curr_internal_row); - global_col_indices[global_row_indices[curr_internal_row + 1]] = global_col_idx; - global_row_indices[curr_internal_row + 1]++; - } + IdxT local_col_idx = columns[elem_idx]; + IdxT global_col_idx = component->getNodeConnection(local_col_idx); - // Need to sort internal rows because even though the mapping from component internal to system internal - // is monotonic, the mapping from component external to system external may not be, and internal rows - // may contain external columns. - auto global_row_start = global_col_indices + global_row_indices[curr_internal_row]; - auto global_row_end = global_col_indices + global_row_indices[curr_internal_row + 1]; - std::sort(global_row_start, global_row_end); + if (global_col_idx == neg1_) + continue; - curr_internal_row++; + global_col_indices[global_row_indices[curr_internal_row + 1]] = global_col_idx; + global_row_indices[curr_internal_row + 1]++; } } + // One last row after + curr_internal_row++; + global_row_indices[curr_internal_row + 1] = global_row_indices[curr_internal_row]; + + // Need to sort internal rows because even though the mapping from component internal to system internal + // is monotonic, the mapping from component external to system external may not be, and internal rows + // may contain external columns. + for (size_t row_idx = 0; row_idx < curr_internal_row; row_idx++) + { + auto global_row_start = global_col_indices + global_row_indices[row_idx]; + auto global_row_end = global_col_indices + global_row_indices[row_idx + 1]; + std::sort(global_row_start, global_row_end); + } + assert(curr_internal_row == n_intern_); // Then move on to external rows, which come after internal rows @@ -252,14 +250,23 @@ namespace GridKit // Collect columns from each component which has a row which contributes to this row for (ComponentContribution contrib : reverse_extern_map[row - n_intern_]) { - component_type* component = components_[contrib.comp_idx_]; - CsrSparsity& comp_sparsity = component_sparsities[contrib.comp_idx_]; + component_type* component = components_[contrib.comp_idx_]; + const auto& [rows, columns, values] = component->getJacobian().getEntries(); + + auto row_start = std::ranges::lower_bound(rows, contrib.local_row_idx_); + + // It can happen where external contributions are only constants, and do not appear in the jacobian. + // If that is the case, we won't be able to find local_row_idx_ and must skip this contribution + if (row_start == rows.end() || *row_start != contrib.local_row_idx_) + continue; + + auto row_end = std::upper_bound(row_start, rows.end(), contrib.local_row_idx_); - IdxT row_start = comp_sparsity.row_indices_[contrib.local_row_idx_]; - IdxT row_end = comp_sparsity.row_indices_[contrib.local_row_idx_ + 1]; - for (size_t local_elem_idx = row_start; local_elem_idx < row_end; local_elem_idx++) + for (size_t local_elem_idx = std::distance(rows.begin(), row_start); + local_elem_idx < static_cast(std::distance(rows.begin(), row_end)); + local_elem_idx++) { - IdxT local_col = comp_sparsity.col_indices_[local_elem_idx]; + IdxT local_col = columns[local_elem_idx]; IdxT global_col = component->getNodeConnection(local_col); // Not all columns map to columns in the system jacobian @@ -292,7 +299,6 @@ namespace GridKit std::copy(global_col_indices, global_col_indices + nnz, csr_jac_.getColData()); std::copy(global_row_indices, global_row_indices + size_ + 1, csr_jac_.getRowData()); - delete[] component_sparsities; delete[] reverse_extern_map; delete[] global_row_indices; delete[] global_col_indices; From 9d33cfe07a80ddeb9bf3a12aa7f57d407a0148eb Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 Jan 2026 15:44:03 -0500 Subject: [PATCH 28/29] Remove auto keyword usage --- .../SystemModelPowerElectronics.hpp | 41 +++++++++++-------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 1e2377ed..443c0426 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -63,6 +63,7 @@ namespace GridKit { using RealT = typename CircuitComponent::RealT; using component_type = CircuitComponent; + using MatrixT = CircuitComponent::MatrixT; using CircuitComponent::size_; using CircuitComponent::n_intern_; @@ -151,8 +152,8 @@ namespace GridKit }; // A reverse mapping from external system variables -> component variables - auto reverse_extern_map = new std::forward_list[n_extern_]; - size_t component_nnz = 0; + std::forward_list* reverse_extern_map = new std::forward_list[n_extern_]; + size_t component_nnz = 0; // Loop over all components, evaluate their jacobians, save their sparsity information, // and construct the reverse variable mapping. @@ -161,7 +162,7 @@ namespace GridKit component_type* component = components_[comp_idx]; component->evaluateJacobian(); - auto& comp_jacobian = component->getJacobian(); + MatrixT& comp_jacobian = component->getJacobian(); for (IdxT local_external_row : component->getExternIndices()) { @@ -190,9 +191,11 @@ namespace GridKit size_t curr_internal_row = 0; for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) { - component_type* component = components_[comp_idx]; - auto comp_externals = component->getExternIndices(); - const auto& [rows, columns, values] = component->getJacobian().getEntries(); + component_type* component = components_[comp_idx]; + std::set comp_externals = component->getExternIndices(); + std::tuple&, std::vector&, std::vector&> entries = component->getJacobian().getEntries(); + const std::vector& rows = std::get<0>(entries); + const std::vector& columns = std::get<1>(entries); if (rows.empty()) continue; @@ -235,8 +238,8 @@ namespace GridKit // may contain external columns. for (size_t row_idx = 0; row_idx < curr_internal_row; row_idx++) { - auto global_row_start = global_col_indices + global_row_indices[row_idx]; - auto global_row_end = global_col_indices + global_row_indices[row_idx + 1]; + IdxT* global_row_start = global_col_indices + global_row_indices[row_idx]; + IdxT* global_row_end = global_col_indices + global_row_indices[row_idx + 1]; std::sort(global_row_start, global_row_end); } @@ -250,17 +253,19 @@ namespace GridKit // Collect columns from each component which has a row which contributes to this row for (ComponentContribution contrib : reverse_extern_map[row - n_intern_]) { - component_type* component = components_[contrib.comp_idx_]; - const auto& [rows, columns, values] = component->getJacobian().getEntries(); + component_type* component = components_[contrib.comp_idx_]; + std::tuple&, std::vector&, std::vector&> entries = component->getJacobian().getEntries(); + const std::vector& rows = std::get<0>(entries); + const std::vector& columns = std::get<1>(entries); - auto row_start = std::ranges::lower_bound(rows, contrib.local_row_idx_); + typename std::vector::const_iterator row_start = std::ranges::lower_bound(rows, contrib.local_row_idx_); // It can happen where external contributions are only constants, and do not appear in the jacobian. // If that is the case, we won't be able to find local_row_idx_ and must skip this contribution if (row_start == rows.end() || *row_start != contrib.local_row_idx_) continue; - auto row_end = std::upper_bound(row_start, rows.end(), contrib.local_row_idx_); + typename std::vector::const_iterator row_end = std::upper_bound(row_start, rows.end(), contrib.local_row_idx_); for (size_t local_elem_idx = std::distance(rows.begin(), row_start); local_elem_idx < static_cast(std::distance(rows.begin(), row_end)); @@ -280,12 +285,12 @@ namespace GridKit // Sort the row by column indices. Since the mapping from local indices to global indices isn't monotonically increasing, // this is necessary. - auto start = global_col_indices + global_row_indices[row]; - auto end = global_col_indices + global_row_indices[row + 1]; + IdxT* start = global_col_indices + global_row_indices[row]; + IdxT* end = global_col_indices + global_row_indices[row + 1]; std::sort(start, end); // De-duplicate the columns - auto new_end = std::unique(start, end); + IdxT* new_end = std::unique(start, end); global_row_indices[row + 1] = global_row_indices[row] + static_cast(std::distance(start, new_end)); } // Allocate new sparsity buffers @@ -348,7 +353,7 @@ namespace GridKit component->allocate(); // Count up the amount of internal variables which get mapped to system variables - auto extern_indices = component->getExternIndices(); + std::set extern_indices = component->getExternIndices(); for (IdxT comp_var_idx = 0; comp_var_idx < component->size(); comp_var_idx++) { IdxT sys_var_idx = component->getNodeConnection(comp_var_idx); @@ -369,8 +374,8 @@ namespace GridKit // sorted by these groupings, so the first component is the first block and so on. for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) { - auto component = components_[comp_idx]; - auto extern_indices = component->getExternIndices(); + component_type* component = components_[comp_idx]; + std::set extern_indices = component->getExternIndices(); // Whether or not we've seen a local variable yet bool has_seen_local = false; From bcb1357a2394dd9a1115b15df9dc8ece3be26152 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 Jan 2026 16:25:13 -0500 Subject: [PATCH 29/29] Put csr_jac_ on the heap Allows for use of "stealing" constructors by constructing a new csr matrix object whenever `allocate()` is called. Known Issue: A warning is erroneously generated because of some methods incorrectly labeled as virtual in `CsrMatrix`. --- .../LinearAlgebra/SparseMatrix/CsrMatrix.cpp | 20 ---------------- .../SystemModelPowerElectronics.hpp | 24 ++++++++++--------- 2 files changed, 13 insertions(+), 31 deletions(-) diff --git a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp index a6291f2b..a1296c57 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp @@ -818,26 +818,6 @@ namespace GridKit } } - /** - * @brief Changes the size of the matrix. De-allocates the matrix, if needed, since the - * old indices corresponded to a differently-size matrix and no longer make sense. - * Does not re-allocate any memory. - * @param n The new number of rows - * @param m The new number of columns - * @return 0 if succesful - */ - template - int CsrMatrix::resize(IdxT n, IdxT m) - { - destroyMatrixData(memory::HOST); - destroyMatrixData(memory::DEVICE); - - n_ = n; - m_ = m; - - return 0; - } - template class CsrMatrix; template class CsrMatrix; } // namespace LinearAlgebra diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 443c0426..98e92243 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -9,6 +9,7 @@ #include #include +#include #include #include #include @@ -116,6 +117,7 @@ namespace GridKit max_steps_ = max_steps; // Can choose if to use jacobian use_jac_ = use_jac; + csr_jac_ = nullptr; } /** @@ -130,6 +132,9 @@ namespace GridKit { for (auto comp : components_) delete comp; + + if (csr_jac_ != nullptr) + delete csr_jac_; } /** @@ -179,7 +184,7 @@ namespace GridKit component_nnz += comp_jacobian.nnz(); } - // Allocate the final sparsity pattern info + // Allocate the final sparsity pattern info. Not deleted, because ownership is stolen by the jacobian IdxT* global_row_indices = new IdxT[size_ + 1]; IdxT* global_col_indices = new IdxT[component_nnz]; // Use component_nnz as an upper bound on nnz global_row_indices[0] = 0; @@ -296,17 +301,14 @@ namespace GridKit // Allocate new sparsity buffers IdxT nnz = global_row_indices[size_]; - csr_jac_.resize(size_, size_); - csr_jac_.setNnz(nnz); - csr_jac_.allocateMatrixData(LinearAlgebra::memory::HOST); + if (csr_jac_ != nullptr) + delete csr_jac_; - // Copy column indices - std::copy(global_col_indices, global_col_indices + nnz, csr_jac_.getColData()); - std::copy(global_row_indices, global_row_indices + size_ + 1, csr_jac_.getRowData()); + RealT* value_buffer = new RealT[nnz]; + csr_jac_ = new LinearAlgebra::CsrMatrix( + size_, size_, nnz, &global_row_indices, &global_col_indices, &value_buffer, LinearAlgebra::memory::HOST, LinearAlgebra::memory::HOST); delete[] reverse_extern_map; - delete[] global_row_indices; - delete[] global_col_indices; return 1; } @@ -646,13 +648,13 @@ namespace GridKit */ LinearAlgebra::CsrMatrix& getCsrJac() { - return csr_jac_; + return *csr_jac_; } private: static constexpr IdxT neg1_ = static_cast(-1); - LinearAlgebra::CsrMatrix csr_jac_; + LinearAlgebra::CsrMatrix* csr_jac_; std::vector components_;