From cbb3294514cd5cdbaf1d649eb238cf209e5421a7 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Thu, 16 Jan 2025 09:53:13 -0500 Subject: [PATCH 01/18] fixed some style issues and variable name selection --- examples/Microgrid/Microgrid.cpp | 4 +- examples/RLCircuit/RLCircuit.cpp | 4 +- examples/ScaleMicrogrid/ScaleMicrogrid.cpp | 4 +- src/LinearAlgebra/COO_Matrix.hpp | 129 +++++++++--------- .../PowerElectronics/CircuitComponent.hpp | 2 +- .../SystemModelPowerElectronics.hpp | 24 ++-- 6 files changed, 84 insertions(+), 83 deletions(-) diff --git a/examples/Microgrid/Microgrid.cpp b/examples/Microgrid/Microgrid.cpp index 901433b1..85b8f15c 100644 --- a/examples/Microgrid/Microgrid.cpp +++ b/examples/Microgrid/Microgrid.cpp @@ -24,10 +24,10 @@ int main(int argc, char const *argv[]) double abstol = 1.0e-8; double reltol = 1.0e-8; size_t max_step_amount = 3000; - bool usejac = true; + bool use_jac = true; //Create model - GridKit::PowerElectronicsModel* sysmodel = new GridKit::PowerElectronicsModel(reltol, abstol, usejac, max_step_amount); + ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(reltol, abstol, use_jac, max_step_amount); //Modeled after the problem in the paper double RN = 1.0e4; diff --git a/examples/RLCircuit/RLCircuit.cpp b/examples/RLCircuit/RLCircuit.cpp index 46afe872..11092959 100644 --- a/examples/RLCircuit/RLCircuit.cpp +++ b/examples/RLCircuit/RLCircuit.cpp @@ -20,11 +20,11 @@ int main(int argc, char const *argv[]) { double abstol = 1.0e-8; double reltol = 1.0e-8; - bool usejac = true; + bool use_jac = true; //TODO:setup as named parameters //Create circuit model - GridKit::PowerElectronicsModel* sysmodel = new GridKit::PowerElectronicsModel(reltol, abstol, usejac); + ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(reltol, abstol, use_jac); size_t idoff = 0; diff --git a/examples/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/ScaleMicrogrid/ScaleMicrogrid.cpp index a145c0be..93fc176c 100644 --- a/examples/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -66,7 +66,7 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) { using namespace GridKit; - bool usejac = true; + bool use_jac = true; real_type t_init = 0.0; real_type t_final = 1.0; @@ -77,7 +77,7 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) // Create circuit model auto* sysmodel = new PowerElectronicsModel(reltol, abstol, - usejac, + use_jac, SCALE_MICROGRID_MAX_STEPS); const std::vector* true_vec = &answer_key_N8; diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 9036085a..646f4f20 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -21,8 +21,8 @@ class COO_Matrix { private: std::vector values_; - std::vector row_indexes_; - std::vector column_indexes_; + std::vector row_indices_; + std::vector column_indices_; Intdx rows_size_; Intdx columns_size_; bool sorted_; @@ -39,8 +39,8 @@ class COO_Matrix // --- Functions which call sort --- std::tuple, std::vector> getRowCopy(Intdx r); std::tuple&, std::vector&, std::vector&> getEntries(); - std::tuple, std::vector, std::vector> getEntrieCopies(); - std::tuple, std::vector, std::vector> getEntrieCopiesSubMatrix(std::vector submap); + std::tuple, std::vector, std::vector> getEntryCopies(); + std::tuple, std::vector, std::vector> getEntryCopiesSubMatrix(std::vector submap); std::tuple, std::vector, std::vector> getDataToCSR(); std::vector getCSRRowData(); @@ -53,7 +53,7 @@ class COO_Matrix ScalarT frobnorm(); // --- Permutation Operations --- - //No sorting is actually done. Only done when nesscary + //No sorting is actually done. Only done when necessary void permutation(std::vector row_perm, std::vector col_perm); void permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n); @@ -107,13 +107,13 @@ inline std::tuple, std::vector> COO_Matrixvalues_.size() && this->row_indexes_[rsize] == r); + } while (rsize < this->values_.size() && this->row_indices_[rsize] == r); - return {{this->column_indexes_.begin() + rowindex, this->column_indexes_.begin() + rsize},{this->values_.begin() + rowindex, this->values_.begin() + rsize}}; + return {{this->column_indices_.begin() + rowindex, this->column_indices_.begin() + rsize},{this->values_.begin() + rowindex, this->values_.begin() + rsize}}; } /** - * @brief Get all Entries pointers. Will sort before returnings + * @brief Get all entry pointers. Will sort before returning * * @tparam ScalarT * @tparam Intdx @@ -126,28 +126,28 @@ inline std::tuple&, std::vector&, std::vector { this->sortSparse(); } - return {this->row_indexes_, this->column_indexes_, this->values_}; + return {this->row_indices_, this->column_indices_, this->values_}; } /** - * @brief Get copies of the data. Sorted_ before returning + * @brief Sorts the data if it's not already sorted * * @tparam ScalarT * @tparam Intdx * @return std::tuple, std::vector, std::vector> */ template -inline std::tuple, std::vector, std::vector> COO_Matrix::getEntrieCopies() +inline std::tuple, std::vector, std::vector> COO_Matrix::getEntryCopies() { if (!this->sorted_) { this->sortSparse(); } - return {this->row_indexes_, this->column_indexes_, this->values_}; + return {this->row_indices_, this->column_indices_, this->values_}; } /** - * @brief Returns the data into CSR Format + * @brief Returns the data in CSR Format * * @tparam ScalarT * @tparam Intdx @@ -162,13 +162,13 @@ inline std::tuple, std::vector, std::vector> for (Intdx i = 0; i < static_cast(rowsizevec.size() - 1); i++) { rowsizevec[i + 1] = rowsizevec[i]; - while (counter < static_cast(this->row_indexes_.size()) && i == this->row_indexes_[counter]) + while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) { rowsizevec[i+1]++; counter++; } } - return {rowsizevec, this->column_indexes_, this->values_}; + return {rowsizevec, this->column_indices_, this->values_}; } /** @@ -190,7 +190,7 @@ inline std::vector COO_Matrix::getCSRRowData() for (Intdx i = 0; i < static_cast(rowsizevec.size() - 1); i++) { rowsizevec[i + 1] = rowsizevec[i]; - while (counter < static_cast(this->row_indexes_.size()) && i == this->row_indexes_[counter]) + while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) { rowsizevec[i+1]++; counter++; @@ -218,13 +218,13 @@ inline void COO_Matrix::setValues(std::vector r, std::vec //Duplicated with axpy. Could replace with function depdent on lambda expression Intdx aiter = 0; //iterate for all current values_ in matrix - for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values_ when they are not in current matrix - while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indexes_[i] || (r[aiter] == this->row_indexes_[i] && c[aiter] < this->column_indexes_[i]))) + while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indices_[i] || (r[aiter] == this->row_indices_[i] && c[aiter] < this->column_indices_[i]))) { - this->row_indexes_.push_back(r[aiter]); - this->column_indexes_.push_back(c[aiter]); + this->row_indices_.push_back(r[aiter]); + this->column_indices_.push_back(c[aiter]); this->values_.push_back(v[aiter]); this->checkIncreaseSize(r[aiter], c[aiter]); aiter++; @@ -235,7 +235,7 @@ inline void COO_Matrix::setValues(std::vector r, std::vec } - if (r[aiter] == this->row_indexes_[i] && c[aiter] == this->column_indexes_[i]) + if (r[aiter] == this->row_indices_[i] && c[aiter] == this->column_indices_[i]) { this->values_[i] = v[aiter]; aiter++; @@ -244,8 +244,8 @@ inline void COO_Matrix::setValues(std::vector r, std::vec //push back rest that was not found sorted for (Intdx i = aiter; i < static_cast(r.size()); i++) { - this->row_indexes_.push_back(r[i]); - this->column_indexes_.push_back(c[i]); + this->row_indices_.push_back(r[i]); + this->column_indices_.push_back(c[i]); this->values_.push_back(v[i]); this->checkIncreaseSize(r[i], c[i]); @@ -285,19 +285,19 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrixrows_size_ = this->rows_size_ > m ? this->rows_size_ : m; this->columns_size_ = this->columns_size_ > n ? this->columns_size_ : n; Intdx aiter = 0; //iterate for all current values_ in matrix - for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values_ when they are not in current matrix - while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indexes_[i] || (r[aiter] == this->row_indexes_[i] && c[aiter] < this->column_indexes_[i]))) + while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indices_[i] || (r[aiter] == this->row_indices_[i] && c[aiter] < this->column_indices_[i]))) { - this->row_indexes_.push_back(r[aiter]); - this->column_indexes_.push_back(c[aiter]); + this->row_indices_.push_back(r[aiter]); + this->column_indices_.push_back(c[aiter]); this->values_.push_back(alpha * val[aiter]); this->checkIncreaseSize(r[aiter], c[aiter]); @@ -309,7 +309,7 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrixrow_indexes_[i] && c[aiter] == this->column_indexes_[i]) + if (r[aiter] == this->row_indices_[i] && c[aiter] == this->column_indices_[i]) { this->values_[i] += alpha * val[aiter]; aiter++; @@ -318,8 +318,8 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix(r.size()); i++) { - this->row_indexes_.push_back(r[i]); - this->column_indexes_.push_back(c[i]); + this->row_indices_.push_back(r[i]); + this->column_indices_.push_back(c[i]); this->values_.push_back(alpha * val[i]); this->checkIncreaseSize(r[i], c[i]); @@ -353,13 +353,13 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r Intdx aiter = 0; //iterate for all current values_ in matrix - for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values_ when they are not in current matrix - while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indexes_[i] || (r[aiter] == this->row_indexes_[i] && c[aiter] < this->column_indexes_[i]))) + while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indices_[i] || (r[aiter] == this->row_indices_[i] && c[aiter] < this->column_indices_[i]))) { - this->row_indexes_.push_back(r[aiter]); - this->column_indexes_.push_back(c[aiter]); + this->row_indices_.push_back(r[aiter]); + this->column_indices_.push_back(c[aiter]); this->values_.push_back(alpha * v[aiter]); this->checkIncreaseSize(r[aiter], c[aiter]); @@ -371,7 +371,7 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r } - if (r[aiter] == this->row_indexes_[i] && c[aiter] == this->column_indexes_[i]) + if (r[aiter] == this->row_indices_[i] && c[aiter] == this->column_indices_[i]) { this->values_[i] += alpha * v[aiter]; aiter++; @@ -380,8 +380,8 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r //push back rest that was not found sorted_ for (Intdx i = aiter; i < static_cast(r.size()); i++) { - this->row_indexes_.push_back(r[i]); - this->column_indexes_.push_back(c[i]); + this->row_indices_.push_back(r[i]); + this->column_indices_.push_back(c[i]); this->values_.push_back(alpha * v[i]); this->checkIncreaseSize(r[i], c[i]); @@ -440,16 +440,16 @@ inline void COO_Matrix::permutation(std::vector row_perm, for (int i = 0; i < this->values_.size(); i++) { - this->row_indexes_[i] = row_perm[this->row_indexes_[i]]; - this->column_indexes_[i] = col_perm[this->column_indexes_[i]]; + this->row_indices_[i] = row_perm[this->row_indices_[i]]; + this->column_indices_[i] = col_perm[this->column_indices_[i]]; } this->sorted_ = false; //cycle sorting maybe useful since permutations are already known } /** - * @brief Permutates the matrix and can change its size efficently - * if size is shrinking and value is to be removed the negative one + * @brief Permutes the matrix and can change its size efficently + * if size is shrinking and value is to be removed, it's set to -1 * * @tparam ScalarT * @tparam Intdx @@ -469,21 +469,21 @@ inline void COO_Matrix::permutationSizeMap(std::vector ro for (int i = 0; i < this->values_.size(); i++) { - if (row_perm[this->row_indexes_[i]] == -1 || col_perm[this->column_indexes_[i]] == -1) + if (row_perm[this->row_indices_[i]] == -1 || col_perm[this->column_indices_[i]] == -1) { this->values_[i] = 0; } else { - this->row_indexes_[i] = row_perm[this->row_indexes_[i]]; - this->column_indexes_[i] = col_perm[this->column_indexes_[i]]; + this->row_indices_[i] = row_perm[this->row_indices_[i]]; + this->column_indices_[i] = col_perm[this->column_indices_[i]]; } } this->sorted_ = false; } /** - * @brief Turn matrix into the zero matrix. Does not actual delete memory + * @brief Turn matrix into the zero matrix. Does not actually delete memory * * @tparam ScalarT * @tparam Intdx @@ -492,8 +492,8 @@ template inline void COO_Matrix::zeroMatrix() { //resize doesn't effect capacity if smaller - this->column_indexes_.resize(0); - this->row_indexes_.resize(0); + this->column_indices_.resize(0); + this->row_indices_.resize(0); this->values_.resize(0); this->sorted_ = true; } @@ -505,8 +505,8 @@ inline void COO_Matrix::identityMatrix(Intdx n) this->zeroMatrix(); for (Intdx i = 0; i < n; i++) { - this->column_indexes_[i] = i; - this->row_indexes_[i] = i; + this->column_indices_[i] = i; + this->row_indices_[i] = i; this->values_[i] = 1.0; } this->sorted_ = true; @@ -521,7 +521,7 @@ inline void COO_Matrix::identityMatrix(Intdx n) template inline void COO_Matrix::sortSparse() { - this->sortSparseCOO(this->row_indexes_, this->column_indexes_, this->values_); + this->sortSparseCOO(this->row_indices_, this->column_indices_, this->values_); this->sorted_ = true; } @@ -544,7 +544,7 @@ inline std::tuple COO_Matrix::getDimensions() } /** - * @brief Print matrix that is sorted_ + * @brief Print matrix in sorted order * * @tparam ScalarT * @tparam Intdx @@ -561,17 +561,17 @@ inline void COO_Matrix::printMatrix() std::cout << "(x , y, value)\n"; for (size_t i = 0; i < this->values_.size(); i++) { - std::cout << "(" << this->row_indexes_[i] - << ", " << this->column_indexes_[i] + std::cout << "(" << this->row_indices_[i] + << ", " << this->column_indices_[i] << ", " << this->values_[i] << ")\n"; } std::cout << std::flush; } /** - * @brief Find the lowest row cordinate from set of provided cordinates + * @brief Find the lowest row cordinate from set of provided coordinates * - * Assumes rows and columns are sorted_ + * Assumes rows and columns are sorted * @tparam ScalarT * @tparam Intdx * @param r @@ -679,9 +679,10 @@ inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) } /** - * @brief Sort a disoreded set of values_. Assume nothing on order. + * @brief Sort a disordered set of values_. Assume nothing on order. * - * @todo simple setup. Should add stable sorting since list are pre-sorted_ + * @todo simple setup. Should add stable sorting since list are pre-sorted_ + * SR: are we assuming something on the order, or not? * * @tparam ScalarT * @tparam Intdx @@ -730,8 +731,8 @@ template inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n) { this->values_ = v; - this->row_indexes_ = r; - this->column_indexes_ = c; + this->row_indices_ = r; + this->column_indices_ = c; this->rows_size_ = m; this->columns_size_ = n; this->sorted_ = false; @@ -743,8 +744,8 @@ inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) this->rows_size_ = m; this->columns_size_ = n; this->values_ = std::vector(); - this->row_indexes_ = std::vector(); - this->column_indexes_ = std::vector(); + this->row_indices_ = std::vector(); + this->column_indices_ = std::vector(); this->sorted_ = false; } @@ -754,8 +755,8 @@ inline COO_Matrix::COO_Matrix() this->rows_size_ = 0; this->columns_size_ = 0; this->values_ = std::vector(); - this->row_indexes_ = std::vector(); - this->column_indexes_ = std::vector(); + this->row_indices_ = std::vector(); + this->column_indices_ = std::vector(); this->sorted_ = false; } diff --git a/src/Model/PowerElectronics/CircuitComponent.hpp b/src/Model/PowerElectronics/CircuitComponent.hpp index 290054c8..fd88925c 100644 --- a/src/Model/PowerElectronics/CircuitComponent.hpp +++ b/src/Model/PowerElectronics/CircuitComponent.hpp @@ -46,7 +46,7 @@ namespace GridKit } /** - * @brief Create the mappings from local to global indexes + * @brief Create the mappings from local to global indices * * @param local_index * @param global_index diff --git a/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 50128b93..5549a84c 100644 --- a/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -51,21 +51,21 @@ namespace GridKit rtol_ = 1e-4; atol_ = 1e-4; this->max_steps_ = 2000; - // By default not use jacobian - usejac_ = false; + // By default don't use the jacobian + use_jac_ = false; } /** * @brief Constructor for the system model */ - PowerElectronicsModel(double rt = 1e-4, double at = 1e-4, bool ju = false, IdxT msa = 2000) : ModelEvaluatorImpl(0, 0, 0) + PowerElectronicsModel(double rtol = 1e-4, double atol = 1e-4, bool use_jac = false, IdxT max_steps = 2000) : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances from input - rtol_ = rt; - atol_ = at; - this->max_steps_ = msa; - // Can choose if to use jacobain - usejac_ = ju; + rtol_ = rtol; + atol_ = atol; + this->max_steps_ = max_steps; + // Can choose if to use jacobian + use_jac_ = use_jac; } /** @@ -91,7 +91,7 @@ namespace GridKit } /** - * @brief Will check if each component has jacobian avalible. If one doesn't then jacobain is false. + * @brief Will check if each component has jacobian avalible. If one doesn't have it, return false. * * * @return true @@ -99,7 +99,7 @@ namespace GridKit */ bool hasJacobian() { - if (!this->usejac_) + if (!this->use_jac_) return false; for (const auto &component : components_) @@ -234,7 +234,7 @@ namespace GridKit { component->evaluateJacobian(); - // get references to local jacobain + // get references to local jacobian std::tuple&, std::vector&, std::vector&> tpm = component->getJacobian().getEntries(); const auto& [r, c, v] = tpm; @@ -325,7 +325,7 @@ namespace GridKit static constexpr IdxT neg1_ = static_cast(-1); std::vector components_; - bool usejac_; + bool use_jac_; }; // class PowerElectronicsModel From 8aed1cb423e4f7ea2163cfdde0bd70e4d6cbaa20 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 17 Jan 2025 16:41:07 -0500 Subject: [PATCH 02/18] further code cleanup and commenting --- examples/Grid3Bus/README.md | 2 +- examples/Microgrid/Microgrid.cpp | 2 +- examples/SparseTest/SparseTest.cpp | 2 +- src/CircuitGraph.hpp | 5 +- src/LinearAlgebra/COO_Matrix.hpp | 304 ++++++++++++------ .../PowerElectronics/Capacitor/Capacitor.cpp | 6 +- .../PowerElectronics/Capacitor/Capacitor.hpp | 2 +- .../DistributedGenerator.cpp | 36 +-- .../DistributedGenerator.hpp | 2 +- .../InductionMotor/InductionMotor.cpp | 4 +- .../InductionMotor/InductionMotor.hpp | 4 +- .../PowerElectronics/Inductor/Inductor.cpp | 6 +- .../PowerElectronics/Inductor/Inductor.hpp | 2 +- .../LinearTransformer/LinearTransformer.hpp | 2 +- .../MicrogridBusDQ/MicrogridBusDQ.cpp | 4 +- .../MicrogridBusDQ/MicrogridBusDQ.hpp | 2 +- .../MicrogridLine/MicrogridLine.cpp | 10 +- .../MicrogridLine/MicrogridLine.hpp | 2 +- .../MicrogridLoad/MicrogridLoad.cpp | 10 +- .../MicrogridLoad/MicrogridLoad.hpp | 2 +- .../PowerElectronics/Resistor/Resistor.cpp | 2 +- .../PowerElectronics/Resistor/Resistor.hpp | 2 +- .../SynchronousMachine/SynchronousMachine.cpp | 4 +- .../SynchronousMachine/SynchronousMachine.hpp | 4 +- .../SystemModelPowerElectronics.hpp | 58 ++-- .../TransmissionLine/TransmissionLine.cpp | 10 +- .../TransmissionLine/TransmissionLine.hpp | 2 +- .../VoltageSource/VoltageSource.cpp | 2 +- .../VoltageSource/VoltageSource.hpp | 2 +- src/Model/PowerFlow/Bus/BaseBus.hpp | 4 +- src/Model/PowerFlow/Bus/BusSlack.hpp | 4 +- src/Model/PowerFlow/MiniGrid/MiniGrid.cpp | 4 +- src/Model/PowerFlow/MiniGrid/MiniGrid.hpp | 4 +- src/Model/PowerFlow/SystemModelPowerFlow.hpp | 12 +- src/ModelEvaluator.hpp | 2 +- src/ModelEvaluatorImpl.hpp | 18 +- src/Solver/Dynamic/Ida.cpp | 26 +- src/SystemModel.hpp | 8 +- 38 files changed, 352 insertions(+), 225 deletions(-) diff --git a/examples/Grid3Bus/README.md b/examples/Grid3Bus/README.md index b628ad42..b420c373 100644 --- a/examples/Grid3Bus/README.md +++ b/examples/Grid3Bus/README.md @@ -57,7 +57,7 @@ with variables $`\theta_2, |V_2|`$ and $`\theta_3`$. Nonlinear solver can approximate Jacobian numerically, however this is computationaly expensive and scales poorly with the size of the problem. Typically, one needs to provide Jacobian in addition to residual to the nonlinear solver. For nonlinear problem defined by (vector) function $`\mathbf{f}(\mathbf{x})=0`$, Jacobian matrix is defined as ```math -J_{i,j}=\frac{\partial f_i}{\partial x_j}, ~~~ i,j=1,\ldots,N +jac_{i,j}=\frac{\partial f_i}{\partial x_j}, ~~~ i,j=1,\ldots,N ``` where $`N`$ is the size of vectors $`\mathbf{f}`$ and $`\mathbf{x}`$. Jacobian for our model is evaluated as ```math diff --git a/examples/Microgrid/Microgrid.cpp b/examples/Microgrid/Microgrid.cpp index 85b8f15c..589b3787 100644 --- a/examples/Microgrid/Microgrid.cpp +++ b/examples/Microgrid/Microgrid.cpp @@ -336,7 +336,7 @@ int main(int argc, char const *argv[]) sysmodel->updateTime(0.0, 1.0e-8); sysmodel->evaluateJacobian(); std::cout << "Intial Jacobian with alpha:\n"; - // std::cout << sysmodel->getJacobian().frobnorm() << "\n"; + // std::cout << sysmodel->getJacobian().frobNorm() << "\n"; //Create numerical integrator and configure it for the generator model diff --git a/examples/SparseTest/SparseTest.cpp b/examples/SparseTest/SparseTest.cpp index c7c07195..32c075c4 100644 --- a/examples/SparseTest/SparseTest.cpp +++ b/examples/SparseTest/SparseTest.cpp @@ -50,7 +50,7 @@ int main(int argc, char const *argv[]) std::vector r; std::vector c; std::vector v; - std::tie(r,c,v) = A.getDataToCSR(); + std::tie(r,c,v) = A.setDataToCSR(); for (size_t i = 0; i < r.size() - 1; i++) { diff --git a/src/CircuitGraph.hpp b/src/CircuitGraph.hpp index 070d838f..c0dde1a2 100644 --- a/src/CircuitGraph.hpp +++ b/src/CircuitGraph.hpp @@ -95,12 +95,13 @@ size_t CircuitGraph::amountHyperEdges() } /** - * @brief Printing + * @brief Print the bipartite graph * * @todo need to add verbose printing for connections display * * @tparam IdxT - * @param verbose + * @param[in] verbose if true will print connections, + * otherwise just the number of nodes and edges */ template diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 646f4f20..62bcd662 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -12,7 +12,7 @@ /** * @brief Quick class to provide sparse matrices of COO type. Simplifies data movement * - * @todo add functionality to keep track of multiple sorted_ list. Faster adding of new entries and will have a threshold to sort completely. + * @todo add functionality to keep track of multiple sorted lists. Faster adding of new entries and will have a threshold to sort completely. * * m x n sparse matrix */ @@ -42,7 +42,7 @@ class COO_Matrix std::tuple, std::vector, std::vector> getEntryCopies(); std::tuple, std::vector, std::vector> getEntryCopiesSubMatrix(std::vector submap); - std::tuple, std::vector, std::vector> getDataToCSR(); + std::tuple, std::vector, std::vector> setDataToCSR(); std::vector getCSRRowData(); // BLAS. Will sort before running @@ -50,10 +50,10 @@ class COO_Matrix void axpy(ScalarT alpha, COO_Matrix& a); void axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v); void scal(ScalarT alpha); - ScalarT frobnorm(); + ScalarT frobNorm(); // --- Permutation Operations --- - //No sorting is actually done. Only done when necessary + //Sorting is only done if not already sorted. void permutation(std::vector row_perm, std::vector col_perm); void permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n); @@ -71,7 +71,7 @@ class COO_Matrix void printMatrix(); - static void sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &vals); + static void sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &values); private: Intdx indexStartRow(const std::vector &rows, Intdx r); @@ -81,11 +81,11 @@ class COO_Matrix }; /** - * @brief Get copy of row values_ + * @brief Get copy of row index * * @tparam ScalarT * @tparam Intdx - * @param r + * @param[in] r row index * @return std::tuple, std::vector> */ template @@ -95,21 +95,21 @@ inline std::tuple, std::vector> COO_MatrixsortSparse(); } - Intdx rowindex = this->indexStartRow(r); + Intdx row_index = this->indexStartRow(r); - if (rowindex == -1) + if (row_index == -1) { return {std::vector(),std::vector()}; } - Intdx rsize = rowindex; + Intdx rsize = row_index; do { rsize++; } while (rsize < this->values_.size() && this->row_indices_[rsize] == r); - return {{this->column_indices_.begin() + rowindex, this->column_indices_.begin() + rsize},{this->values_.begin() + rowindex, this->values_.begin() + rsize}}; + return {{this->column_indices_.begin() + row_index, this->column_indices_.begin() + rsize},{this->values_.begin() + row_index, this->values_.begin() + rsize}}; } /** @@ -154,21 +154,21 @@ inline std::tuple, std::vector, std::vector> * @return std::tuple, std::vector, std::vector> */ template -inline std::tuple, std::vector, std::vector> COO_Matrix::getDataToCSR() +inline std::tuple, std::vector, std::vector> COO_Matrix::setDataToCSR() { if (!this->isSorted()) this->sortSparse(); - std::vector rowsizevec(this->rows_size_ + 1, 0); + std::vector row_size_vec(this->rows_size_ + 1, 0); Intdx counter = 0; - for (Intdx i = 0; i < static_cast(rowsizevec.size() - 1); i++) + for (Intdx i = 0; i < static_cast(row_size_vec.size() - 1); i++) { - rowsizevec[i + 1] = rowsizevec[i]; + row_size_vec[i + 1] = row_size_vec[i]; while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) { - rowsizevec[i+1]++; + row_size_vec[i+1]++; counter++; } } - return {rowsizevec, this->column_indices_, this->values_}; + return {row_size_vec, this->column_indices_, this->values_}; } /** @@ -185,28 +185,33 @@ template inline std::vector COO_Matrix::getCSRRowData() { if (!this->isSorted()) this->sortSparse(); - std::vector rowsizevec(this->rows_size_ + 1, 0); + std::vector row_size_vec(this->rows_size_ + 1, 0); Intdx counter = 0; - for (Intdx i = 0; i < static_cast(rowsizevec.size() - 1); i++) + for (Intdx i = 0; i < static_cast(row_size_vec.size() - 1); i++) { - rowsizevec[i + 1] = rowsizevec[i]; + row_size_vec[i + 1] = row_size_vec[i]; while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) { - rowsizevec[i+1]++; + row_size_vec[i+1]++; counter++; } } - return rowsizevec; + return row_size_vec; } /** - * @brief Given set of vector data it will set the values_ into the matrix + * @brief Set coordinates and values of the matrix. Will sort before returning * * @tparam ScalarT * @tparam Intdx - * @param r - * @param c - * @param v + * @param[in] r row indices of the matrix + * @param[in] c column indices of the matrix + * @param[in] v values of the matrix + * + * @pre r.size() == c.size() == v.size() + * @pre r,c,v represent an array in COO format + * + * @post Coordinates and values are set in the matrix. */ template inline void COO_Matrix::setValues(std::vector r, std::vector c, std::vector v) @@ -216,33 +221,33 @@ inline void COO_Matrix::setValues(std::vector r, std::vec //Duplicated with axpy. Could replace with function depdent on lambda expression - Intdx aiter = 0; + Intdx a_iter = 0; //iterate for all current values_ in matrix for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values_ when they are not in current matrix - while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indices_[i] || (r[aiter] == this->row_indices_[i] && c[aiter] < this->column_indices_[i]))) + while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) { - this->row_indices_.push_back(r[aiter]); - this->column_indices_.push_back(c[aiter]); - this->values_.push_back(v[aiter]); - this->checkIncreaseSize(r[aiter], c[aiter]); - aiter++; + this->row_indices_.push_back(r[a_iter]); + this->column_indices_.push_back(c[a_iter]); + this->values_.push_back(v[a_iter]); + this->checkIncreaseSize(r[a_iter], c[a_iter]); + a_iter++; } - if (aiter >= static_cast(r.size())) + if (a_iter >= static_cast(r.size())) { break; } - if (r[aiter] == this->row_indices_[i] && c[aiter] == this->column_indices_[i]) + if (r[a_iter] == this->row_indices_[i] && c[a_iter] == this->column_indices_[i]) { - this->values_[i] = v[aiter]; - aiter++; + this->values_[i] = v[a_iter]; + a_iter++; } } //push back rest that was not found sorted - for (Intdx i = aiter; i < static_cast(r.size()); i++) + for (Intdx i = a_iter; i < static_cast(r.size()); i++) { this->row_indices_.push_back(r[i]); this->column_indices_.push_back(c[i]); @@ -256,12 +261,14 @@ inline void COO_Matrix::setValues(std::vector r, std::vec } /** - * @brief BLAS axpy operation on another COO matrix. Will sort both matrices before acting + * @brief Implements axpy this += alpha * a. Will sort before running * * @tparam ScalarT * @tparam Intdx - * @param alpha - * @param a + * @param[in] alpha matrix to be added + * @param[in] a scalar to multiply by + * + * @post this = this + alpha * a */ template inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix& a) @@ -289,34 +296,34 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrixrows_size_ = this->rows_size_ > m ? this->rows_size_ : m; this->columns_size_ = this->columns_size_ > n ? this->columns_size_ : n; - Intdx aiter = 0; - //iterate for all current values_ in matrix + Intdx a_iter = 0; + //iterate for all current values in matrix for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) { - //pushback values_ when they are not in current matrix - while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indices_[i] || (r[aiter] == this->row_indices_[i] && c[aiter] < this->column_indices_[i]))) + //pushback values when they are not in current matrix + while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) { - this->row_indices_.push_back(r[aiter]); - this->column_indices_.push_back(c[aiter]); - this->values_.push_back(alpha * val[aiter]); + this->row_indices_.push_back(r[a_iter]); + this->column_indices_.push_back(c[a_iter]); + this->values_.push_back(alpha * val[a_iter]); - this->checkIncreaseSize(r[aiter], c[aiter]); - aiter++; + this->checkIncreaseSize(r[a_iter], c[a_iter]); + a_iter++; } - if (aiter >= static_cast(r.size())) + if (a_iter >= static_cast(r.size())) { break; } - if (r[aiter] == this->row_indices_[i] && c[aiter] == this->column_indices_[i]) + if (r[a_iter] == this->row_indices_[i] && c[a_iter] == this->column_indices_[i]) { - this->values_[i] += alpha * val[aiter]; - aiter++; + this->values_[i] += alpha * val[a_iter]; + a_iter++; } } //push back rest that was not found sorted_ - for (Intdx i = aiter; i < static_cast(r.size()); i++) + for (Intdx i = a_iter; i < static_cast(r.size()); i++) { this->row_indices_.push_back(r[i]); this->column_indices_.push_back(c[i]); @@ -329,14 +336,19 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix inline void COO_Matrix::axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v) @@ -351,34 +363,34 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r //sort input this->sortSparseCOO(r, c, v); - Intdx aiter = 0; + Intdx a_iter = 0; //iterate for all current values_ in matrix for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values_ when they are not in current matrix - while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indices_[i] || (r[aiter] == this->row_indices_[i] && c[aiter] < this->column_indices_[i]))) + while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) { - this->row_indices_.push_back(r[aiter]); - this->column_indices_.push_back(c[aiter]); - this->values_.push_back(alpha * v[aiter]); + this->row_indices_.push_back(r[a_iter]); + this->column_indices_.push_back(c[a_iter]); + this->values_.push_back(alpha * v[a_iter]); - this->checkIncreaseSize(r[aiter], c[aiter]); - aiter++; + this->checkIncreaseSize(r[a_iter], c[a_iter]); + a_iter++; } - if (aiter >= static_cast(r.size())) + if (a_iter >= static_cast(r.size())) { break; } - if (r[aiter] == this->row_indices_[i] && c[aiter] == this->column_indices_[i]) + if (r[a_iter] == this->row_indices_[i] && c[a_iter] == this->column_indices_[i]) { - this->values_[i] += alpha * v[aiter]; - aiter++; + this->values_[i] += alpha * v[a_iter]; + a_iter++; } } //push back rest that was not found sorted_ - for (Intdx i = aiter; i < static_cast(r.size()); i++) + for (Intdx i = a_iter; i < static_cast(r.size()); i++) { this->row_indices_.push_back(r[i]); this->column_indices_.push_back(c[i]); @@ -391,11 +403,11 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r } /** - * @brief Scale all values_ + * @brief Scale all values by alpha * * @tparam ScalarT * @tparam Intdx - * @param alpha + * @param[in] alpha scalar to scale by */ template inline void COO_Matrix::scal(ScalarT alpha) @@ -407,14 +419,14 @@ inline void COO_Matrix::scal(ScalarT alpha) } /** - * @brief Frobenius Norm of the Matrix + * @brief Calculates the Frobenius Norm of the matrix * * @tparam ScalarT * @tparam Intdx - * @return ScalarT + * @return ScalarT - Frobenius Norm of the matrix */ template -inline ScalarT COO_Matrix::frobnorm() +inline ScalarT COO_Matrix::frobNorm() { ScalarT totsum = 0.0; for (auto i = this->values_.begin(); i < this->values_.end(); i++) @@ -429,8 +441,12 @@ inline ScalarT COO_Matrix::frobnorm() * * @tparam ScalarT * @tparam Intdx - * @param row_perm - * @param col_perm + * @param[in] row_perm + * @param[out] col_perm + * + * @pre row_perm.size() == this->rows_size_ = col_perm.size() == this->columns_size_ + * + * @post this = this(row_perm, col_perm) */ template inline void COO_Matrix::permutation(std::vector row_perm, std::vector col_perm) @@ -449,14 +465,19 @@ inline void COO_Matrix::permutation(std::vector row_perm, /** * @brief Permutes the matrix and can change its size efficently - * if size is shrinking and value is to be removed, it's set to -1 * * @tparam ScalarT * @tparam Intdx - * @param row_perm size of m - * @param col_perm size of n - * @param m - * @param n + * @param[in] row_perm row permutation + * @param[in] col_perm column permutation + * @param[in] m number of rows + * @param[in] n number of columns + * + * @pre row_perm.size() == this->rows_size_ + * @pre col_perm.size() == this->columns_size_ + * @pre indices are set to -1 if they are to be removed + * + * @post this = this(row_perm, col_perm) and removed indices have corresponding values set to 0 */ template inline void COO_Matrix::permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n) @@ -485,8 +506,11 @@ inline void COO_Matrix::permutationSizeMap(std::vector ro /** * @brief Turn matrix into the zero matrix. Does not actually delete memory * + * @SR - If it's not deleteing memory, what's the point? + * * @tparam ScalarT * @tparam Intdx + * */ template inline void COO_Matrix::zeroMatrix() @@ -498,6 +522,19 @@ inline void COO_Matrix::zeroMatrix() this->sorted_ = true; } +/** + * @brief Turn matrix into the identity matrix + * + * @tparam ScalarT + * @tparam Intdx + * + * @param[in] n size of the identity matrix + * + * @post this = I_n + * + * @SR - it might be better to explicitly zero out the matrix and require to be so in preconditions + */ + template inline void COO_Matrix::identityMatrix(Intdx n) { @@ -525,12 +562,29 @@ inline void COO_Matrix::sortSparse() this->sorted_ = true; } +/** + * @brief Check if the matrix is sorted + * + * @tparam ScalarT + * @tparam Intdx + * + * @param[out] bool - true if sorted, false otherwise + */ + template inline bool COO_Matrix::isSorted() { return this->sorted_; } +/** + * @brief Get the number of non-zero elements in the matrix + * + * @tparam ScalarT + * @tparam Intdx + * + * @param[out] Intdx - number of non-zero elements in the matrix + */ template inline Intdx COO_Matrix::nnz() { @@ -574,8 +628,11 @@ inline void COO_Matrix::printMatrix() * Assumes rows and columns are sorted * @tparam ScalarT * @tparam Intdx - * @param r - * @return Intdx + * + * @param[in] rows - row indices + * @param[in] r - row index + * + * @return Intdx - index of lowest row */ template inline Intdx COO_Matrix::indexStartRow(const std::vector &rows, Intdx r) @@ -617,11 +674,11 @@ inline Intdx COO_Matrix::indexStartRow(const std::vector * * @tparam ScalarT * @tparam Intdx - * @param rows - * @param columns - * @param ri - * @param ci - * @return Intdx + * @param[in] rows - row indices + * @param[in] columns - column indices + * @param[in] ri - row index + * @param[in] ci - column index + * @return Intdx - returns the index of the coordinate */ template inline Intdx COO_Matrix::sparseCordBinarySearch(const std::vector &rows, const std::vector &columns, Intdx ri, Intdx ci) @@ -660,6 +717,16 @@ inline Intdx COO_Matrix::sparseCordBinarySearch(const std::vecto return m; } +/** + * @brief Check if the size of the matrix needs to be increased + * + * @tparam ScalarT + * @tparam Intdx + * @param[in] r row index + * @param[in] c column index + * @return true if size was increased + */ + template inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) { @@ -679,19 +746,19 @@ inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) } /** - * @brief Sort a disordered set of values_. Assume nothing on order. + * @brief Sort a disordered set of values. Assume nothing on order. * - * @todo simple setup. Should add stable sorting since list are pre-sorted_ + * @todo simple setup. Should add stable sorting since lists are pre-sorted_ * SR: are we assuming something on the order, or not? * * @tparam ScalarT * @tparam Intdx * @param rows * @param columns - * @param values_ + * @param values */ template -inline void COO_Matrix::sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &vals) +inline void COO_Matrix::sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &values) { //index based sort code @@ -718,7 +785,7 @@ inline void COO_Matrix::sortSparseCOO(std::vector &rows, { std::swap(rows[ordervec[i]], rows[ordervec[ordervec[i]]]); std::swap(columns[ordervec[i]], columns[ordervec[ordervec[i]]]); - std::swap(vals[ordervec[i]], vals[ordervec[ordervec[i]]]); + std::swap(values[ordervec[i]], values[ordervec[ordervec[i]]]); //swap orderings std::swap(ordervec[i], ordervec[ordervec[i]]); @@ -727,6 +794,24 @@ inline void COO_Matrix::sortSparseCOO(std::vector &rows, } } +/** + * @brief Constructor for COO Matrix with given cooridnates and values + * + * @tparam ScalarT + * @tparam Intdx + * + * @param[in] r row indices + * @param[in] c column indices + * @param[in] v values + * @param[in] m number of rows + * @param[in] n number of columns + * + * @pre r.size() == c.size() == v.size() + * @pre r,c,v represent an array in COO format + * + * @post COO_Matrix is created with given coordinates and values + */ + template inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n) { @@ -735,9 +820,21 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vector< this->column_indices_ = c; this->rows_size_ = m; this->columns_size_ = n; - this->sorted_ = false; + this->sorted_ = false; - //SR: Why is this assumed false? } +/** + * @brief Constructor for empty COO Matrix of a given size + * + * @tparam ScalarT + * @tparam Intdx + * + * @param[in] m number of rows + * @param[in] n number of columns + * + * @post empty COO Matrix is created with given size + */ + template inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) { @@ -746,9 +843,18 @@ inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = false; + this->sorted_ = false; //SR: I think this would be true, since it's empty. } +/** + * @brief Constructor for empty COO Matrix of size 0 + * + * @tparam ScalarT + * @tparam Intdx + * + * @post empty COO Matrix of size 0 is created + */ + template inline COO_Matrix::COO_Matrix() { @@ -757,7 +863,7 @@ inline COO_Matrix::COO_Matrix() this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = false; + this->sorted_ = false; //SR: I think this would be true, since it's empty. } template diff --git a/src/Model/PowerElectronics/Capacitor/Capacitor.cpp b/src/Model/PowerElectronics/Capacitor/Capacitor.cpp index 6a884446..ca09e389 100644 --- a/src/Model/PowerElectronics/Capacitor/Capacitor.cpp +++ b/src/Model/PowerElectronics/Capacitor/Capacitor.cpp @@ -90,12 +90,12 @@ int Capacitor::evaluateResidual() template int Capacitor::evaluateJacobian() { - J_.zeroMatrix(); + jac_.zeroMatrix(); //Create dF/dy std::vector rcord{2,2,2}; std::vector ccord{0,1,2}; std::vector vals{1.0, -1.0, -1.0}; - J_.setValues(rcord, ccord, vals); + jac_.setValues(rcord, ccord, vals); //Create dF/dy' std::vector rcordder{0,1,2}; @@ -104,7 +104,7 @@ int Capacitor::evaluateJacobian() COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,3,3); //Perform dF/dy + \alpha dF/dy' - J_.axpy(alpha_, Jacder); + jac_.axpy(alpha_, Jacder); return 0; } diff --git a/src/Model/PowerElectronics/Capacitor/Capacitor.hpp b/src/Model/PowerElectronics/Capacitor/Capacitor.hpp index 145ad807..ad2f9bdb 100644 --- a/src/Model/PowerElectronics/Capacitor/Capacitor.hpp +++ b/src/Model/PowerElectronics/Capacitor/Capacitor.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index 4ab6ddd3..68cd9ae2 100644 --- a/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -180,7 +180,7 @@ int DistributedGenerator::evaluateResidual() template int DistributedGenerator::evaluateJacobian() { - J_.zeroMatrix(); + jac_.zeroMatrix(); //Create dF/dy' std::vector rcordder(13); std::vector valsder(13, -1.0); @@ -203,7 +203,7 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(1); valtemp = { - sin(y_[3]) * y_[14] - cos(y_[3]) * y_[15], cos(y_[3]),-sin(y_[3])}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 2 @@ -211,7 +211,7 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(2); valtemp = { cos(y_[3]) * y_[14] - sin(y_[3]) * y_[15], sin(y_[3]),cos(y_[3])}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 3 @@ -219,7 +219,7 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(3); valtemp = {-1.0, -mp_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 0 if (refframe_) @@ -228,7 +228,7 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(0); valtemp = {-1.0, -mp_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); } @@ -237,63 +237,63 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(4); valtemp = {-wc_, wc_*y_[14], wc_*y_[15], wc_*y_[12], wc_*y_[13]}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 5 ctemp = {5, 12, 13, 14, 15}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(5); valtemp = {-wc_, -wc_*y_[15], wc_*y_[14], wc_*y_[13], -wc_*y_[12]}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 6 ctemp = {5, 12}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(6); valtemp = {-nq_, -1.0}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 7 ctemp = {13}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(7); valtemp = {-1.0}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 8 ctemp = {5,6,10,12,13,14}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(8); valtemp = {-Kpv_*nq_, Kiv_, -1.0, -Kpv_, -Cf_*wb_, F_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 9 ctemp = {7, 11, 12, 13, 15}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(9); valtemp = {Kiv_, -1.0, Cf_*wb_,-Kpv_,F_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 10 ctemp = {4, 5, 6, 8, 10, 11, 12, 13, 14}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(10); valtemp = {-mp_ * y_[11], -(Kpc_ * Kpv_ * nq_) / Lf_, (Kpc_ * Kiv_) / Lf_, Kic_ / Lf_, -(Kpc_ + rLf_) / Lf_, -mp_ * y_[4], -(Kpc_ * Kpv_ + 1.0) / Lf_, -(Cf_ * Kpc_ * wb_) / Lf_, (F_ * Kpc_) / Lf_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 11 ctemp = {4, 7, 9, 10, 11, 12, 13, 15}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(11); valtemp = {mp_ * y_[10], (Kiv_ * Kpc_) / Lf_, Kic_ / Lf_, mp_ * y_[4], -(Kpc_ + rLf_) / Lf_, (Cf_ * Kpc_ * wb_) / Lf_, -(Kpc_ * Kpv_ + 1.0) / Lf_, (F_ * Kpc_) / Lf_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 12 ctemp = {4, 10, 13, 14}; rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(12); valtemp = {-mp_ * y_[13], 1.0 / Cf_, wb_ - mp_ * y_[4], -1.0 / Cf_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 13 @@ -301,7 +301,7 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(13); valtemp = {mp_ * y_[12], 1.0 / Cf_, -wb_ + mp_ * y_[4], -1.0 / Cf_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 14 @@ -309,7 +309,7 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(14); valtemp = {(1.0/Lc_) * -cos(y_[3]) , (1.0/Lc_) * -sin(y_[3]) , (1.0/Lc_) * (sin(y_[3]) * y_[1] - cos(y_[3]) * y_[2]), -mp_ * y_[15], 1.0 / Lc_, -rLc_ / Lc_, wb_ - mp_ * y_[4]}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //r = 15 @@ -317,12 +317,12 @@ int DistributedGenerator::evaluateJacobian() rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(15); valtemp = {(1.0/Lc_) * sin(y_[3]) , (1.0/Lc_) * -cos(y_[3]), (1.0/Lc_) * (cos(y_[3]) * y_[1] + sin(y_[3]) * y_[2]), mp_ * y_[14], 1.0 / Lc_, -wb_ + mp_ * y_[4], -rLc_ / Lc_}; - J_.setValues(rtemp, ctemp, valtemp); + jac_.setValues(rtemp, ctemp, valtemp); //Perform dF/dy + \alpha dF/dy' - J_.axpy(alpha_, Jacder); + jac_.axpy(alpha_, Jacder); return 0; } diff --git a/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp index b2a96cde..b5cd28ae 100644 --- a/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/src/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -56,7 +56,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp b/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp index 4b66e7e9..1ec8f0ca 100644 --- a/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp +++ b/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp @@ -26,7 +26,7 @@ InductionMotor::InductionMotor(IdxT id, ScalarT Lls, ScalarT Rs, Llr_(Llr), Rr_(Rr), Lms_(Lms), - RJ_(RJ), + Rjac_(RJ), P_(P) { size_ = 10; @@ -82,7 +82,7 @@ int InductionMotor::evaluateResidual() f_[0] = y_[5] + y_[7]; f_[1] = (-1.0/2.0) * y_[5] - (sqrt(3.0)/2.0)*y_[6] + y_[7]; f_[2] = (-1.0/2.0) * y_[5] + (sqrt(3.0)/2.0)*y_[6] + y_[7]; - f_[3] = RJ_ * yp_[3] - (3.0/4.0)*P_*Lms_ * (y_[5]*y_[9] - y_[6]*y_[8]); + f_[3] = Rjac_ * yp_[3] - (3.0/4.0)*P_*Lms_ * (y_[5]*y_[9] - y_[6]*y_[8]); f_[4] = yp_[4] - y_[3]; f_[5] = (1.0/3.0)*(2.0* y_[0] - y_[1] - y_[2]) - Rs_*y_[5] - (Lls_ + Lms_) * yp_[5] - Lms_ * yp_[6]; f_[6] = (1.0/sqrt(3.0))*(-y_[1] + y_[2]) - Rs_*y_[6] - (Lls_ + Lms_) * yp_[6] - Lms_ * yp_[5]; diff --git a/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp b/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp index 5a3002d5..76fd1513 100644 --- a/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp +++ b/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; @@ -66,7 +66,7 @@ namespace GridKit ScalarT Llr_; ScalarT Rr_; ScalarT Lms_; - ScalarT RJ_; + ScalarT Rjac_; ScalarT P_; }; } diff --git a/src/Model/PowerElectronics/Inductor/Inductor.cpp b/src/Model/PowerElectronics/Inductor/Inductor.cpp index b157af0e..6da117fc 100644 --- a/src/Model/PowerElectronics/Inductor/Inductor.cpp +++ b/src/Model/PowerElectronics/Inductor/Inductor.cpp @@ -88,13 +88,13 @@ int Inductor::evaluateResidual() template int Inductor::evaluateJacobian() { - J_.zeroMatrix(); + jac_.zeroMatrix(); //Create dF/dy std::vector rcord{0,1,2,2}; std::vector ccord{2,2,0,1}; std::vector vals{-1.0, 1.0, -1.0, 1.0}; - J_.setValues(rcord, ccord, vals); + jac_.setValues(rcord, ccord, vals); //Create dF/dy' std::vector rcordder{2}; @@ -103,7 +103,7 @@ int Inductor::evaluateJacobian() COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,3,3); //Perform dF/dy + \alpha dF/dy' - J_.axpy(alpha_, Jacder); + jac_.axpy(alpha_, Jacder); return 0; } diff --git a/src/Model/PowerElectronics/Inductor/Inductor.hpp b/src/Model/PowerElectronics/Inductor/Inductor.hpp index ec196fb6..8f2a0485 100644 --- a/src/Model/PowerElectronics/Inductor/Inductor.hpp +++ b/src/Model/PowerElectronics/Inductor/Inductor.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp b/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp index f11f5e31..c2e18248 100644 --- a/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp +++ b/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp index 063b9e0e..75c41be9 100644 --- a/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp +++ b/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -93,13 +93,13 @@ int MicrogridBusDQ::evaluateResidual() template int MicrogridBusDQ::evaluateJacobian() { - J_.zeroMatrix(); + jac_.zeroMatrix(); //Create dF/dy std::vector rtemp{0,1}; std::vector ctemp{0,1}; std::vector vals{-1.0 / RN_,-1.0 / RN_}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); return 0; } diff --git a/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index f8dd201d..880cc494 100644 --- a/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/src/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index 21deeb80..222b8871 100644 --- a/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -105,26 +105,26 @@ int MicrogridLine::evaluateResidual() template int MicrogridLine::evaluateJacobian() { - J_.zeroMatrix(); + jac_.zeroMatrix(); //Create dF/dy std::vector rtemp{1,2,3,4}; std::vector ctemp{5,6,5,6}; std::vector vals{-1.0,-1.0,1.0,1.0}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); std::vector ccord{0, 1, 3, 5, 6}; std::vector rcord(ccord.size(),5); vals = {y_[6], (1.0 / L_) , -(1.0 / L_), - (R_ / L_) , y_[0]}; - J_.setValues(rcord, ccord, vals); + jac_.setValues(rcord, ccord, vals); std::vector ccor2{0, 2, 4, 5, 6}; std::fill(rcord.begin(), rcord.end(), 6); vals = {-y_[5], (1.0 / L_) , -(1.0 / L_), -y_[0], - (R_ / L_)}; - J_.setValues(rcord, ccor2, vals); + jac_.setValues(rcord, ccor2, vals); //Create -dF/dy' @@ -134,7 +134,7 @@ int MicrogridLine::evaluateJacobian() COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,7,7); //Perform dF/dy + \alpha dF/dy' - J_.axpy(alpha_, Jacder); + jac_.axpy(alpha_, Jacder); return 0; diff --git a/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp index cf81f888..01affa64 100644 --- a/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/src/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -35,7 +35,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index 36508ee5..9b4b8643 100644 --- a/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -101,26 +101,26 @@ int MicrogridLoad::evaluateResidual() template int MicrogridLoad::evaluateJacobian() { - J_.zeroMatrix(); + jac_.zeroMatrix(); //Create dF/dy std::vector rtemp{1,2}; std::vector ctemp{3,4}; std::vector vals{-1.0,-1.0}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); std::vector ccord{0, 1, 3, 4}; std::vector rcord(ccord.size(),3); vals = {y_[4], (1.0 / L_) , - (R_ / L_) , y_[0]}; - J_.setValues(rcord, ccord, vals); + jac_.setValues(rcord, ccord, vals); std::vector ccor2{0, 2, 3, 4}; std::fill(rcord.begin(), rcord.end(), 4); vals = {-y_[3], (1.0 / L_) , -y_[0], - (R_ / L_)}; - J_.setValues(rcord, ccor2, vals); + jac_.setValues(rcord, ccor2, vals); //Create -dF/dy' @@ -130,7 +130,7 @@ int MicrogridLoad::evaluateJacobian() COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,5,5); //Perform dF/dy + \alpha dF/dy' - J_.axpy(alpha_, Jacder); + jac_.axpy(alpha_, Jacder); return 0; } diff --git a/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 9fc5befb..e292bbf1 100644 --- a/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/src/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/Resistor/Resistor.cpp b/src/Model/PowerElectronics/Resistor/Resistor.cpp index 564c1219..843d161d 100644 --- a/src/Model/PowerElectronics/Resistor/Resistor.cpp +++ b/src/Model/PowerElectronics/Resistor/Resistor.cpp @@ -84,7 +84,7 @@ int Resistor::evaluateJacobian() std::vector rcord{0,0,1,1}; std::vector ccord{0,1,0,1}; std::vector vals{1.0 / R_, -1.0 / R_, -1.0 / R_, 1.0 / R_}; - J_.setValues(rcord, ccord, vals); + jac_.setValues(rcord, ccord, vals); return 0; } diff --git a/src/Model/PowerElectronics/Resistor/Resistor.hpp b/src/Model/PowerElectronics/Resistor/Resistor.hpp index 885373f8..822f5b48 100644 --- a/src/Model/PowerElectronics/Resistor/Resistor.hpp +++ b/src/Model/PowerElectronics/Resistor/Resistor.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index cc0e0d84..762268bc 100644 --- a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp @@ -31,7 +31,7 @@ SynchronousMachine::SynchronousMachine(IdxT id, ScalarT Lls, std: Rkq_(Rkq), Rfd_(Rfd), Rkd_(Rkd), - RJ_(RJ), + Rjac_(RJ), P_(P), mub_(mub) { @@ -100,7 +100,7 @@ int SynchronousMachine::evaluateResidual() f_[0] = y_[6]*cos1 + y_[7]*sin1 + y_[8]; f_[1] = y_[6]*cos23m + y_[7]*sin23m + y_[8]; f_[2] = y_[6]*cos23p + y_[7]*sin23p + y_[8]; - f_[3] = RJ_ * yp_[4] - (3.0/4.0)*P_*(Lmd_ *y_[6]* (y_[7] + y_[11] + y_[12]) - Lmq_*y_[7]*(y_[6] + y_[9] + y_[0])); + f_[3] = Rjac_ * yp_[4] - (3.0/4.0)*P_*(Lmd_ *y_[6]* (y_[7] + y_[11] + y_[12]) - Lmq_*y_[7]*(y_[6] + y_[9] + y_[0])); f_[4] = yp_[5] - y_[4]; f_[5] = (-2.0/3.0)*(y_[0]*cos1 +y_[1]*cos23m + y_[2]*cos23p) + Rs_*y_[6] + (Lls_ + Lmq_)*yp_[6] + Lmq_*yp_[9] + Lmq_*yp_[10] + y_[4]*(P_/2.0)*((Lls_ + Lmd_)*y_[7] + Lmd_*y_[11] + Lmd_*y_[12]); f_[6] = (-2.0/3.0)*(y_[0]*sin1 -y_[1]*sin23m - y_[2]*sin23p) + Rs_*y_[7] + (Lls_ + Lmd_)*yp_[7] + Lmd_*yp_[11] + Lmd_*yp_[12] - y_[4]*(P_/2.0)*((Lls_ + Lmq_)*y_[6] + Lmq_*y_[9] + Lmq_*y_[10]); diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp index 190ebad6..5bf88b6c 100644 --- a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp +++ b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp @@ -37,7 +37,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; @@ -72,7 +72,7 @@ namespace GridKit std::tuple Rkq_; ScalarT Rfd_; ScalarT Rkd_; - ScalarT RJ_; + ScalarT Rjac_; ScalarT P_; ScalarT mub_; }; diff --git a/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 5549a84c..c1c1a220 100644 --- a/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -34,9 +34,9 @@ namespace GridKit // using ModelEvaluatorImpl::fB_; // using ModelEvaluatorImpl::g_; // using ModelEvaluatorImpl::gB_; - using ModelEvaluatorImpl::J_; - using ModelEvaluatorImpl::rtol_; - using ModelEvaluatorImpl::atol_; + using ModelEvaluatorImpl::jac_; + using ModelEvaluatorImpl::rel_tol_; + using ModelEvaluatorImpl::abs_tol_; // using ModelEvaluatorImpl::param_; // using ModelEvaluatorImpl::param_up_; // using ModelEvaluatorImpl::param_lo_; @@ -44,12 +44,14 @@ namespace GridKit public: /** * @brief Default constructor for the system model + * + * @post System model parameters set as default */ PowerElectronicsModel() : ModelEvaluatorImpl(0, 0, 0) { - // Set system model tolerances as default - rtol_ = 1e-4; - atol_ = 1e-4; + // Set system model parameters as default + rel_tol_ = 1e-4; + abs_tol_ = 1e-4; this->max_steps_ = 2000; // By default don't use the jacobian use_jac_ = false; @@ -57,12 +59,19 @@ namespace GridKit /** * @brief Constructor for the system model + * + * @param[in] rel_tol Relative tolerance for the system model + * @param[in] abs_tol Absolute tolerance for the system model + * @param[in] use_jac Boolean to choose if to use jacobian + * @param[in] max_steps Maximum number of steps for the system model + * + * @post System model parameters set as input */ - PowerElectronicsModel(double rtol = 1e-4, double atol = 1e-4, bool use_jac = false, IdxT max_steps = 2000) : ModelEvaluatorImpl(0, 0, 0) + PowerElectronicsModel(double rel_tol = 1e-4, double abs_tol = 1e-4, bool use_jac = false, IdxT max_steps = 2000) : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances from input - rtol_ = rtol; - atol_ = atol; + rel_tol_ = rel_tol; + abs_tol_ = abs_tol; this->max_steps_ = max_steps; // Can choose if to use jacobian use_jac_ = use_jac; @@ -70,6 +79,11 @@ namespace GridKit /** * @brief Destructor for the system model + * + * @pre System components are allocated + * + * @post System components are deallocated + * */ virtual ~PowerElectronicsModel() { @@ -94,8 +108,8 @@ namespace GridKit * @brief Will check if each component has jacobian avalible. If one doesn't have it, return false. * * - * @return true - * @return false + * @return true if all components have jacobian + * @return false otherwise */ bool hasJacobian() { @@ -116,8 +130,11 @@ namespace GridKit * @brief Allocate the vector data with size amount * @todo Add capability to go through component model connection to get the size of the actual vector * - * @param s - * @return int + * @param[in] s size of the vector + * + * @post System model vectors allocated with size s + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ int allocate(IdxT s) { @@ -140,7 +157,7 @@ namespace GridKit /** * @brief Set intial y and y' of each component * - * @return int + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ int initialize() { @@ -158,7 +175,9 @@ namespace GridKit /** * @brief Distribute y and y' to each component based of node connection graph * - * @return int + * @post Each component has y and y' set + * + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ int distributeVectors() { @@ -189,7 +208,7 @@ namespace GridKit /** * @brief Evaluate Residuals at each component then collect them * - * @return int + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ int evaluateResidual() { @@ -222,11 +241,11 @@ namespace GridKit /** * @brief Creates the Sparse COO Jacobian representing \alpha dF/dy' + dF/dy * - * @return int + * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ int evaluateJacobian() { - J_.zeroMatrix(); + jac_.zeroMatrix(); distributeVectors(); // Evaluate component jacs @@ -253,7 +272,8 @@ namespace GridKit } // AXPY to Global Jacobian - J_.axpy(1.0, rgr, cgr, vgr); + // elementwise jac_(rgr, cgr) += vgr + jac_.axpy(1.0, rgr, cgr, vgr); } return 0; diff --git a/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp b/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp index 4dfc3320..b836217a 100644 --- a/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp +++ b/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp @@ -141,29 +141,29 @@ int TransmissionLine::evaluateJacobian() std::vector rtemp{0,1,2,3,4,5,6,7,8,9,10,11}; std::vector ctemp{8,9,10,11,8,9,10,11,8,9,10,11}; std::vector vals{1.0,1.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); std::vector ccord{0,1,2,3,4,5,6,7}; std::vector rcord(ccord.size(),8); vals = {YReMat_, -YImMatDi_ ,-YReMat_, -YImMatOff_,-YReMat_, YImMatDi_ ,YReMat_, YImMatOff_}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); std::fill(rcord.begin(), rcord.end(), 9); vals = {YImMatDi_ ,YReMat_, YImMatOff_, -YReMat_,-YImMatDi_ ,-YReMat_, -YImMatOff_, YReMat_}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); std::fill(rcord.begin(), rcord.end(), 10); vals = {-YReMat_, -YImMatDi_ ,YReMat_, -YImMatOff_,YReMat_, YImMatDi_ ,-YReMat_, YImMatOff_}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); std::fill(rcord.begin(), rcord.end(), 11); vals = {YImMatDi_ ,-YReMat_, YImMatOff_, YReMat_,-YImMatDi_ ,YReMat_, -YImMatOff_, -YReMat_}; - J_.setValues(rtemp, ctemp, vals); + jac_.setValues(rtemp, ctemp, vals); return 0; } diff --git a/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp b/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp index 4da149e9..0ab0b340 100644 --- a/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp +++ b/src/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp @@ -40,7 +40,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerElectronics/VoltageSource/VoltageSource.cpp b/src/Model/PowerElectronics/VoltageSource/VoltageSource.cpp index ec9d5b19..cddb777b 100644 --- a/src/Model/PowerElectronics/VoltageSource/VoltageSource.cpp +++ b/src/Model/PowerElectronics/VoltageSource/VoltageSource.cpp @@ -83,7 +83,7 @@ int VoltageSource::evaluateJacobian() std::vector rcord{0,1,2,2}; std::vector ccord{2,2,0,1}; std::vector vals{-1.0, 1.0, -1.0, 1.0}; - J_.setValues(rcord, ccord, vals); + jac_.setValues(rcord, ccord, vals); return 0; } diff --git a/src/Model/PowerElectronics/VoltageSource/VoltageSource.hpp b/src/Model/PowerElectronics/VoltageSource/VoltageSource.hpp index 8b15381f..591865fa 100644 --- a/src/Model/PowerElectronics/VoltageSource/VoltageSource.hpp +++ b/src/Model/PowerElectronics/VoltageSource/VoltageSource.hpp @@ -36,7 +36,7 @@ namespace GridKit using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/src/Model/PowerFlow/Bus/BaseBus.hpp b/src/Model/PowerFlow/Bus/BaseBus.hpp index 87247e81..b3b2ed5c 100644 --- a/src/Model/PowerFlow/Bus/BaseBus.hpp +++ b/src/Model/PowerFlow/Bus/BaseBus.hpp @@ -84,8 +84,8 @@ namespace GridKit using ModelEvaluatorImpl::nnz_; using ModelEvaluatorImpl::time_; using ModelEvaluatorImpl::alpha_; - using ModelEvaluatorImpl::rtol_; - using ModelEvaluatorImpl::atol_; + using ModelEvaluatorImpl::rel_tol_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::yp_; using ModelEvaluatorImpl::tag_; diff --git a/src/Model/PowerFlow/Bus/BusSlack.hpp b/src/Model/PowerFlow/Bus/BusSlack.hpp index 61e56a17..4a88594e 100644 --- a/src/Model/PowerFlow/Bus/BusSlack.hpp +++ b/src/Model/PowerFlow/Bus/BusSlack.hpp @@ -82,8 +82,8 @@ namespace GridKit using BaseBus::yp_; using BaseBus::f_; using BaseBus::g_; - using BaseBus::atol_; - using BaseBus::rtol_; + using BaseBus::abs_tol_; + using BaseBus::rel_tol_; public: using real_type = typename ModelEvaluatorImpl::real_type; diff --git a/src/Model/PowerFlow/MiniGrid/MiniGrid.cpp b/src/Model/PowerFlow/MiniGrid/MiniGrid.cpp index f50ebb0b..462abbce 100644 --- a/src/Model/PowerFlow/MiniGrid/MiniGrid.cpp +++ b/src/Model/PowerFlow/MiniGrid/MiniGrid.cpp @@ -87,8 +87,8 @@ MiniGrid::MiniGrid() B23_( 12.0) { //std::cout << "Create a load model with " << size_ << " variables ...\n"; - rtol_ = 1e-5; - atol_ = 1e-5; + rel_tol_ = 1e-5; + abs_tol_ = 1e-5; } template diff --git a/src/Model/PowerFlow/MiniGrid/MiniGrid.hpp b/src/Model/PowerFlow/MiniGrid/MiniGrid.hpp index ea81e416..5cd33bd6 100644 --- a/src/Model/PowerFlow/MiniGrid/MiniGrid.hpp +++ b/src/Model/PowerFlow/MiniGrid/MiniGrid.hpp @@ -75,8 +75,8 @@ namespace GridKit using ModelEvaluatorImpl::time_; using ModelEvaluatorImpl::y_; using ModelEvaluatorImpl::f_; - using ModelEvaluatorImpl::rtol_; - using ModelEvaluatorImpl::atol_; + using ModelEvaluatorImpl::rel_tol_; + using ModelEvaluatorImpl::abs_tol_; typedef typename ModelEvaluatorImpl::real_type real_type; diff --git a/src/Model/PowerFlow/SystemModelPowerFlow.hpp b/src/Model/PowerFlow/SystemModelPowerFlow.hpp index 281668f6..1a7bd214 100644 --- a/src/Model/PowerFlow/SystemModelPowerFlow.hpp +++ b/src/Model/PowerFlow/SystemModelPowerFlow.hpp @@ -110,8 +110,8 @@ class SystemSteadyStateModel : public ModelEvaluatorImpl // using ModelEvaluatorImpl::fB_; // using ModelEvaluatorImpl::g_; // using ModelEvaluatorImpl::gB_; - using ModelEvaluatorImpl::rtol_; - using ModelEvaluatorImpl::atol_; + using ModelEvaluatorImpl::rel_tol_; + using ModelEvaluatorImpl::abs_tol_; // using ModelEvaluatorImpl::param_; // using ModelEvaluatorImpl::param_up_; // using ModelEvaluatorImpl::param_lo_; @@ -123,8 +123,8 @@ class SystemSteadyStateModel : public ModelEvaluatorImpl SystemSteadyStateModel() : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances - rtol_ = 1e-5; - atol_ = 1e-5; + rel_tol_ = 1e-5; + abs_tol_ = 1e-5; } /** @@ -134,8 +134,8 @@ class SystemSteadyStateModel : public ModelEvaluatorImpl */ SystemSteadyStateModel(GridKit::PowerSystemData::SystemModelData mp) : ModelEvaluatorImpl(0,0,0) { - rtol_ = 1e-5; - atol_ = 1e-5; + rel_tol_ = 1e-5; + abs_tol_ = 1e-5; //add buses for(auto busdata : mp.bus) diff --git a/src/ModelEvaluator.hpp b/src/ModelEvaluator.hpp index 28ed6531..dde5130b 100644 --- a/src/ModelEvaluator.hpp +++ b/src/ModelEvaluator.hpp @@ -106,7 +106,7 @@ namespace GridKit virtual IdxT size_quad() = 0; virtual IdxT size_opt() = 0; virtual void updateTime(real_type t, real_type a) = 0; - virtual void setTolerances(real_type& rtol, real_type& atol) const = 0; + virtual void setTolerances(real_type& rel_tol, real_type& abs_tol) const = 0; virtual void setMaxSteps(IdxT& msa) const = 0; virtual std::vector& y() = 0; diff --git a/src/ModelEvaluatorImpl.hpp b/src/ModelEvaluatorImpl.hpp index 87bc6a40..e3b15863 100644 --- a/src/ModelEvaluatorImpl.hpp +++ b/src/ModelEvaluatorImpl.hpp @@ -94,7 +94,7 @@ namespace GridKit ypB_(size_), fB_(size_), gB_(size_opt_), - J_(COO_Matrix()), + jac_(COO_Matrix()), param_(size_opt_), param_up_(size_opt_), param_lo_(size_opt_) @@ -133,10 +133,10 @@ namespace GridKit // std::cout << "updateTime: t = " << time_ << "\n"; // } - virtual void setTolerances(real_type& rtol, real_type& atol) const + virtual void setTolerances(real_type& rel_tol, real_type& abs_tol) const { - rtol = rtol_; - atol = atol_; + rel_tol = rel_tol_; + abs_tol = abs_tol_; } virtual void setMaxSteps(IdxT& msa) const @@ -236,12 +236,12 @@ namespace GridKit COO_Matrix& getJacobian() { - return J_; + return jac_; } const COO_Matrix& getJacobian() const { - return J_; + return jac_; } std::vector& getIntegrand() @@ -299,7 +299,7 @@ namespace GridKit std::vector fB_; std::vector gB_; - COO_Matrix J_; + COO_Matrix jac_; std::vector param_; std::vector param_up_; @@ -308,8 +308,8 @@ namespace GridKit real_type time_; real_type alpha_; - real_type rtol_; - real_type atol_; + real_type rel_tol_; + real_type abs_tol_; IdxT max_steps_; diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index bc38376e..0c81ada1 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -149,11 +149,11 @@ namespace Sundials checkOutput(retval, "IDASetUserData"); // Set tolerances - sunrealtype rtol; - sunrealtype atol; + sunrealtype rel_tol; + sunrealtype abs_tol; - model_->setTolerances(rtol, atol); ///< \todo Function name should be "getTolerances"! - retval = IDASStolerances(solver_, rtol, atol); + model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"! + retval = IDASStolerances(solver_, rel_tol, abs_tol); checkOutput(retval, "IDASStolerances"); IdxT msa; @@ -324,11 +324,11 @@ namespace Sundials checkOutput(retval, "IDAQuadInit"); // Set tolerances and error control for quadratures - real_type rtol, atol; - model_->setTolerances(rtol, atol); + real_type rel_tol, abs_tol; + model_->setTolerances(rel_tol, abs_tol); // Set tolerances for quadrature stricter than for integration - retval = IDAQuadSStolerances(solver_, rtol*0.1, atol*0.1); + retval = IDAQuadSStolerances(solver_, rel_tol*0.1, abs_tol*0.1); checkOutput(retval, "IDAQuadSStolerances"); // Include quadrature in eror checking @@ -434,8 +434,8 @@ namespace Sundials int Ida::initializeBackwardSimulation(real_type tf) { int retval = 0; - sunrealtype rtol; - sunrealtype atol; + sunrealtype rel_tol; + sunrealtype abs_tol; model_->initializeAdjoint(); @@ -450,8 +450,8 @@ namespace Sundials retval = IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_); checkOutput(retval, "IDAInitB"); - model_->setTolerances(rtol, atol); - retval = IDASStolerancesB(solver_, backwardID_, rtol, atol); + model_->setTolerances(rel_tol, abs_tol); + retval = IDASStolerancesB(solver_, backwardID_, rel_tol, abs_tol); checkOutput(retval, "IDASStolerancesB"); retval = IDASetUserDataB(solver_, backwardID_, model_); @@ -476,8 +476,8 @@ namespace Sundials retval = IDAQuadInitB(solver_, backwardID_, this->adjointIntegrand, qB_); checkOutput(retval, "IDAQuadInitB"); - //retval = IDAQuadSStolerancesB(solver_, backwardID_, rtol*1.1, atol*1.1); - retval = IDAQuadSStolerancesB(solver_, backwardID_, rtol*0.1, atol*0.1); + //retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol*1.1, abs_tol*1.1); + retval = IDAQuadSStolerancesB(solver_, backwardID_, rel_tol*0.1, abs_tol*0.1); checkOutput(retval, "IDAQuadSStolerancesB"); // Include quadratures in error control diff --git a/src/SystemModel.hpp b/src/SystemModel.hpp index 4dee38e0..9efeef12 100644 --- a/src/SystemModel.hpp +++ b/src/SystemModel.hpp @@ -102,8 +102,8 @@ class SystemModel : public ModelEvaluatorImpl using ModelEvaluatorImpl::fB_; using ModelEvaluatorImpl::g_; using ModelEvaluatorImpl::gB_; - using ModelEvaluatorImpl::rtol_; - using ModelEvaluatorImpl::atol_; + using ModelEvaluatorImpl::rel_tol_; + using ModelEvaluatorImpl::abs_tol_; using ModelEvaluatorImpl::param_; using ModelEvaluatorImpl::param_up_; using ModelEvaluatorImpl::param_lo_; @@ -115,8 +115,8 @@ class SystemModel : public ModelEvaluatorImpl SystemModel() : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances - rtol_ = 1e-7; - atol_ = 1e-9; + rel_tol_ = 1e-7; + abs_tol_ = 1e-9; this->max_steps_=2000; } From 82f37525887ff3e370dd08c4fd93acca54a4e75f Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Wed, 22 Jan 2025 10:00:25 -0500 Subject: [PATCH 03/18] installed enzyme and ran checks --- src/LinearAlgebra/COO_Matrix.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 62bcd662..9cc65254 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -820,7 +820,7 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vector< this->column_indices_ = c; this->rows_size_ = m; this->columns_size_ = n; - this->sorted_ = false; - //SR: Why is this assumed false? + this->sorted_ = false; //SR: Why is this assumed false? } /** From f81acaafaf6c0d258785a478c0ae6177d5a9a353 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Wed, 22 Jan 2025 14:40:38 -0500 Subject: [PATCH 04/18] addressing comments --- examples/Microgrid/Microgrid.cpp | 2 -- examples/RLCircuit/RLCircuit.cpp | 6 +++--- examples/ScaleMicrogrid/ScaleMicrogrid.cpp | 8 ++++---- examples/ScaleMicrogrid/SolutionKeys.hpp | 2 +- src/LinearAlgebra/COO_Matrix.hpp | 10 +++++++--- .../PowerElectronics/SystemModelPowerElectronics.hpp | 5 ++++- 6 files changed, 19 insertions(+), 14 deletions(-) diff --git a/examples/Microgrid/Microgrid.cpp b/examples/Microgrid/Microgrid.cpp index 589b3787..326fd5e1 100644 --- a/examples/Microgrid/Microgrid.cpp +++ b/examples/Microgrid/Microgrid.cpp @@ -336,8 +336,6 @@ int main(int argc, char const *argv[]) sysmodel->updateTime(0.0, 1.0e-8); sysmodel->evaluateJacobian(); std::cout << "Intial Jacobian with alpha:\n"; - // std::cout << sysmodel->getJacobian().frobNorm() << "\n"; - //Create numerical integrator and configure it for the generator model AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(sysmodel); diff --git a/examples/RLCircuit/RLCircuit.cpp b/examples/RLCircuit/RLCircuit.cpp index 11092959..bfe54e82 100644 --- a/examples/RLCircuit/RLCircuit.cpp +++ b/examples/RLCircuit/RLCircuit.cpp @@ -18,13 +18,13 @@ int main(int argc, char const *argv[]) { - double abstol = 1.0e-8; - double reltol = 1.0e-8; + double abs_tol = 1.0e-8; + double rel_tol = 1.0e-8; bool use_jac = true; //TODO:setup as named parameters //Create circuit model - ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(reltol, abstol, use_jac); + auto* sysmodel = new ModelLib::PowerElectronicsModel(rel_tol, abs_tol, use_jac); size_t idoff = 0; diff --git a/examples/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/ScaleMicrogrid/ScaleMicrogrid.cpp index 93fc176c..496eb8c0 100644 --- a/examples/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -71,12 +71,12 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) real_type t_init = 0.0; real_type t_final = 1.0; - real_type reltol = SCALE_MICROGRID_REL_TOL; - real_type abstol = SCALE_MICROGRID_ABS_TOL; + real_type rel_tol = SCALE_MICROGRID_REL_TOL; + real_type abs_tol = SCALE_MICROGRID_ABS_TOL; // Create circuit model - auto* sysmodel = new PowerElectronicsModel(reltol, - abstol, + auto* sysmodel = new PowerElectronicsModel(rel_tol, + abs_tol, use_jac, SCALE_MICROGRID_MAX_STEPS); diff --git a/examples/ScaleMicrogrid/SolutionKeys.hpp b/examples/ScaleMicrogrid/SolutionKeys.hpp index 1d0ca8e6..6a74dc0d 100644 --- a/examples/ScaleMicrogrid/SolutionKeys.hpp +++ b/examples/ScaleMicrogrid/SolutionKeys.hpp @@ -5,7 +5,7 @@ * @brief Answer keys for Scaled Microgrid test with Nsize = 2, 4, 8 * * Data generated with Matlab ode23tb solver with tolerances set to - * abstol = 1e-12 and reltol = 1e-12 for the ODE derivation of the model. + * abs_tol = 1e-12 and rel_tol = 1e-12 for the ODE derivation of the model. * No index reduction was preformed to get to the ODE model. * * @note This file is only to be included in ScaleMicrogrid.cpp. It has no diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 9cc65254..5286cde5 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -200,7 +200,9 @@ inline std::vector COO_Matrix::getCSRRowData() } /** - * @brief Set coordinates and values of the matrix. Will sort before returning + * @brief Set coordinates and values of the matrix. + * + * Matrix entries will be sorted in row-major order before the method returns. * * @tparam ScalarT * @tparam Intdx @@ -746,10 +748,12 @@ inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) } /** - * @brief Sort a disordered set of values. Assume nothing on order. + * @brief Sorts unordered COO matrix + * + * Matrix entries can appear in arbitrary order and will be sorted in + * row-major order before the method returns. * * @todo simple setup. Should add stable sorting since lists are pre-sorted_ - * SR: are we assuming something on the order, or not? * * @tparam ScalarT * @tparam Intdx diff --git a/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp index c1c1a220..29fa037b 100644 --- a/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -67,7 +67,10 @@ namespace GridKit * * @post System model parameters set as input */ - PowerElectronicsModel(double rel_tol = 1e-4, double abs_tol = 1e-4, bool use_jac = false, IdxT max_steps = 2000) : ModelEvaluatorImpl(0, 0, 0) + PowerElectronicsModel(double rel_tol = 1e-4, + double abs_tol = 1e-4, + bool use_jac = false, + IdxT max_steps = 2000) : ModelEvaluatorImpl(0, 0, 0) { // Set system model tolerances from input rel_tol_ = rel_tol; From ff90b921bd3dee1e91e94bdb7512eecc8b43d505 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Wed, 22 Jan 2025 14:42:39 -0500 Subject: [PATCH 05/18] addressing comments --- src/LinearAlgebra/COO_Matrix.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 5286cde5..43cdbf1c 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -824,7 +824,7 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vector< this->column_indices_ = c; this->rows_size_ = m; this->columns_size_ = n; - this->sorted_ = false; //SR: Why is this assumed false? + this->sorted_ = true; //Any new entry can assume that the matrix is sorted when being added. } /** @@ -847,7 +847,7 @@ inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = false; //SR: I think this would be true, since it's empty. + this->sorted_ = true; //Any new entry can assume that the matrix is sorted when being added. } /** @@ -867,7 +867,7 @@ inline COO_Matrix::COO_Matrix() this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = false; //SR: I think this would be true, since it's empty. + this->sorted_ = true; //Any new entry can assume that the matrix is sorted when being added. } template From ccaa6c48e021a8f78ac82978e647cac25d2adfef Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Wed, 22 Jan 2025 15:54:36 -0500 Subject: [PATCH 06/18] fixed sparse matrix --- src/LinearAlgebra/COO_Matrix.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 43cdbf1c..2e15d1e6 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -824,7 +824,7 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vector< this->column_indices_ = c; this->rows_size_ = m; this->columns_size_ = n; - this->sorted_ = true; //Any new entry can assume that the matrix is sorted when being added. + this->sorted_ = false; //Any new entry can assume that the matrix is sorted when being added. } /** From 5ab9bbc2146ab8115a08e570c5b7451057402ab2 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 09:31:08 -0500 Subject: [PATCH 07/18] addressing comments --- examples/Microgrid/Microgrid.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/Microgrid/Microgrid.cpp b/examples/Microgrid/Microgrid.cpp index 326fd5e1..e6254720 100644 --- a/examples/Microgrid/Microgrid.cpp +++ b/examples/Microgrid/Microgrid.cpp @@ -21,13 +21,13 @@ int main(int argc, char const *argv[]) { ///@todo Needs to be modified. Some components are small relative to others thus there error is high (or could be matlab vector issue) - double abstol = 1.0e-8; - double reltol = 1.0e-8; + double abs_tol = 1.0e-8; + double rel_tol = 1.0e-8; size_t max_step_amount = 3000; bool use_jac = true; //Create model - ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(reltol, abstol, use_jac, max_step_amount); + ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(rel_tol, abs_tol, use_jac, max_step_amount); //Modeled after the problem in the paper double RN = 1.0e4; From 5d7fc9623e2d23e131f0e164216f0e6c5ae77608 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 10:03:57 -0500 Subject: [PATCH 08/18] addressing comments --- src/LinearAlgebra/COO_Matrix.hpp | 16 +++++++------ .../SynchronousMachine/SynchronousMachine.cpp | 23 +++++++++++++++---- 2 files changed, 28 insertions(+), 11 deletions(-) diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 2e15d1e6..9839c42c 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -1,4 +1,4 @@ - +#pragma once #include #include @@ -508,8 +508,6 @@ inline void COO_Matrix::permutationSizeMap(std::vector ro /** * @brief Turn matrix into the zero matrix. Does not actually delete memory * - * @SR - If it's not deleteing memory, what's the point? - * * @tparam ScalarT * @tparam Intdx * @@ -534,7 +532,7 @@ inline void COO_Matrix::zeroMatrix() * * @post this = I_n * - * @SR - it might be better to explicitly zero out the matrix and require to be so in preconditions + * @todo - it might be better to explicitly zero out the matrix and require to be so in preconditions */ template @@ -752,6 +750,10 @@ inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) * * Matrix entries can appear in arbitrary order and will be sorted in * row-major order before the method returns. + * Duplicate entries are not allowed and should be pre-summed. + * + * @pre rows, columns, and values are of the same size and represent a COO matrix with no duplicates + * @post Matrix entries are sorted in row-major order * * @todo simple setup. Should add stable sorting since lists are pre-sorted_ * @@ -824,7 +826,7 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vector< this->column_indices_ = c; this->rows_size_ = m; this->columns_size_ = n; - this->sorted_ = false; //Any new entry can assume that the matrix is sorted when being added. + this->sorted_ = false; // An empty matrix is a sorted matrix. } /** @@ -847,7 +849,7 @@ inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = true; //Any new entry can assume that the matrix is sorted when being added. + this->sorted_ = true; // An empty matrix is a sorted matrix. } /** @@ -867,7 +869,7 @@ inline COO_Matrix::COO_Matrix() this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = true; //Any new entry can assume that the matrix is sorted when being added. + this->sorted_ = true; // An empty matrix is a sorted matrix. } template diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index 762268bc..feb1952c 100644 --- a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp @@ -1,6 +1,3 @@ - - - #include #include #include @@ -15,8 +12,26 @@ namespace GridKit { * @brief Constructor for a constant SynchronousMachine model * * Calls default ModelEvaluatorImpl constructor. - * @todo This models equations are not finish + * @todo This model's equations are not finished * @todo needs to be tested for correctness + * + * @tparam ScalarT - floating point type for the model + * @tparam IdxT - integer index type for the model + * + * @param[in] id - unique identifier for the component + * @param[in] Lls - stator leakage inductance + * @param[in] Llkq - tuple of rotor leakage inductances + * @param[in] Llfd - field leakage inductance + * @param[in] Llkd - damper leakage inductance + * @param[in] Lmq - quadrature axis magnetizing inductance + * @param[in] Lmd - direct axis magnetizing inductance + * @param[in] Rs - stator resistance + * @param[in] Rkq - tuple of rotor resistances + * @param[in] Rfd - field resistance + * @param[in] Rkd - damper resistance + * @param[in] RJ - rotor moment of inertia + * @param[in] P - number of pole pairs + * @param[in] mub - magnetic saturation coefficient */ template From b259843f8672972d7fc198396c636435d5771003 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 10:08:26 -0500 Subject: [PATCH 09/18] fixed index type --- src/LinearAlgebra/COO_Matrix.hpp | 310 +++++++++++++++---------------- 1 file changed, 155 insertions(+), 155 deletions(-) diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 9839c42c..95f2ce03 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -16,20 +16,20 @@ * * m x n sparse matrix */ -template +template class COO_Matrix { private: std::vector values_; - std::vector row_indices_; - std::vector column_indices_; - Intdx rows_size_; - Intdx columns_size_; + std::vector row_indices_; + std::vector column_indices_; + IdxT rows_size_; + IdxT columns_size_; bool sorted_; public: //Constructors - COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n); - COO_Matrix(Intdx m, Intdx n); + COO_Matrix(std::vector r, std::vector c, std::vector v, IdxT m, IdxT n); + COO_Matrix(IdxT m, IdxT n); COO_Matrix(); ~COO_Matrix(); @@ -37,46 +37,46 @@ class COO_Matrix //Operations // --- Functions which call sort --- - std::tuple, std::vector> getRowCopy(Intdx r); - std::tuple&, std::vector&, std::vector&> getEntries(); - std::tuple, std::vector, std::vector> getEntryCopies(); - std::tuple, std::vector, std::vector> getEntryCopiesSubMatrix(std::vector submap); + std::tuple, std::vector> getRowCopy(IdxT r); + std::tuple&, std::vector&, std::vector&> getEntries(); + std::tuple, std::vector, std::vector> getEntryCopies(); + std::tuple, std::vector, std::vector> getEntryCopiesSubMatrix(std::vector submap); - std::tuple, std::vector, std::vector> setDataToCSR(); - std::vector getCSRRowData(); + std::tuple, std::vector, std::vector> setDataToCSR(); + std::vector getCSRRowData(); // BLAS. Will sort before running - void setValues(std::vector r, std::vector c, std::vector v); - void axpy(ScalarT alpha, COO_Matrix& a); - void axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v); + void setValues(std::vector r, std::vector c, std::vector v); + void axpy(ScalarT alpha, COO_Matrix& a); + void axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v); void scal(ScalarT alpha); ScalarT frobNorm(); // --- Permutation Operations --- //Sorting is only done if not already sorted. - void permutation(std::vector row_perm, std::vector col_perm); - void permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n); + void permutation(std::vector row_perm, std::vector col_perm); + void permutationSizeMap(std::vector row_perm, std::vector col_perm, IdxT m, IdxT n); void zeroMatrix(); - void identityMatrix(Intdx n); + void identityMatrix(IdxT n); //Resort values_ void sortSparse(); bool isSorted(); - Intdx nnz(); + IdxT nnz(); - std::tuple getDimensions(); + std::tuple getDimensions(); void printMatrix(); - static void sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &values); + static void sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &values); private: - Intdx indexStartRow(const std::vector &rows, Intdx r); - Intdx sparseCordBinarySearch(const std::vector &rows, const std::vector &columns, Intdx ri, Intdx ci); - bool checkIncreaseSize(Intdx r, Intdx c); + IdxT indexStartRow(const std::vector &rows, IdxT r); + IdxT sparseCordBinarySearch(const std::vector &rows, const std::vector &columns, IdxT ri, IdxT ci); + bool checkIncreaseSize(IdxT r, IdxT c); }; @@ -84,26 +84,26 @@ class COO_Matrix * @brief Get copy of row index * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] r row index - * @return std::tuple, std::vector> + * @return std::tuple, std::vector> */ -template -inline std::tuple, std::vector> COO_Matrix::getRowCopy(Intdx r) +template +inline std::tuple, std::vector> COO_Matrix::getRowCopy(IdxT r) { if (!this->sorted_) { this->sortSparse(); } - Intdx row_index = this->indexStartRow(r); + IdxT row_index = this->indexStartRow(r); if (row_index == -1) { - return {std::vector(),std::vector()}; + return {std::vector(),std::vector()}; } - Intdx rsize = row_index; + IdxT rsize = row_index; do { rsize++; @@ -116,11 +116,11 @@ inline std::tuple, std::vector> COO_Matrix, std::vector, std::vector> + * @tparam IdxT + * @return std::tuple, std::vector, std::vector> */ -template -inline std::tuple&, std::vector&, std::vector&> COO_Matrix::getEntries() +template +inline std::tuple&, std::vector&, std::vector&> COO_Matrix::getEntries() { if (!this->sorted_) { @@ -133,11 +133,11 @@ inline std::tuple&, std::vector&, std::vector * @brief Sorts the data if it's not already sorted * * @tparam ScalarT - * @tparam Intdx - * @return std::tuple, std::vector, std::vector> + * @tparam IdxT + * @return std::tuple, std::vector, std::vector> */ -template -inline std::tuple, std::vector, std::vector> COO_Matrix::getEntryCopies() +template +inline std::tuple, std::vector, std::vector> COO_Matrix::getEntryCopies() { if (!this->sorted_) { @@ -150,19 +150,19 @@ inline std::tuple, std::vector, std::vector> * @brief Returns the data in CSR Format * * @tparam ScalarT - * @tparam Intdx - * @return std::tuple, std::vector, std::vector> + * @tparam IdxT + * @return std::tuple, std::vector, std::vector> */ -template -inline std::tuple, std::vector, std::vector> COO_Matrix::setDataToCSR() +template +inline std::tuple, std::vector, std::vector> COO_Matrix::setDataToCSR() { if (!this->isSorted()) this->sortSparse(); - std::vector row_size_vec(this->rows_size_ + 1, 0); - Intdx counter = 0; - for (Intdx i = 0; i < static_cast(row_size_vec.size() - 1); i++) + std::vector row_size_vec(this->rows_size_ + 1, 0); + IdxT counter = 0; + for (IdxT i = 0; i < static_cast(row_size_vec.size() - 1); i++) { row_size_vec[i + 1] = row_size_vec[i]; - while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) + while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) { row_size_vec[i+1]++; counter++; @@ -178,19 +178,19 @@ inline std::tuple, std::vector, std::vector> * * * @tparam ScalarT - * @tparam Intdx - * @return std::vector + * @tparam IdxT + * @return std::vector */ -template -inline std::vector COO_Matrix::getCSRRowData() +template +inline std::vector COO_Matrix::getCSRRowData() { if (!this->isSorted()) this->sortSparse(); - std::vector row_size_vec(this->rows_size_ + 1, 0); - Intdx counter = 0; - for (Intdx i = 0; i < static_cast(row_size_vec.size() - 1); i++) + std::vector row_size_vec(this->rows_size_ + 1, 0); + IdxT counter = 0; + for (IdxT i = 0; i < static_cast(row_size_vec.size() - 1); i++) { row_size_vec[i + 1] = row_size_vec[i]; - while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) + while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) { row_size_vec[i+1]++; counter++; @@ -205,7 +205,7 @@ inline std::vector COO_Matrix::getCSRRowData() * Matrix entries will be sorted in row-major order before the method returns. * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] r row indices of the matrix * @param[in] c column indices of the matrix * @param[in] v values of the matrix @@ -215,20 +215,20 @@ inline std::vector COO_Matrix::getCSRRowData() * * @post Coordinates and values are set in the matrix. */ -template -inline void COO_Matrix::setValues(std::vector r, std::vector c, std::vector v) +template +inline void COO_Matrix::setValues(std::vector r, std::vector c, std::vector v) { //sort input this->sortSparseCOO(r, c, v); //Duplicated with axpy. Could replace with function depdent on lambda expression - Intdx a_iter = 0; + IdxT a_iter = 0; //iterate for all current values_ in matrix - for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) + for (IdxT i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values_ when they are not in current matrix - while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) + while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) { this->row_indices_.push_back(r[a_iter]); this->column_indices_.push_back(c[a_iter]); @@ -236,7 +236,7 @@ inline void COO_Matrix::setValues(std::vector r, std::vec this->checkIncreaseSize(r[a_iter], c[a_iter]); a_iter++; } - if (a_iter >= static_cast(r.size())) + if (a_iter >= static_cast(r.size())) { break; } @@ -249,7 +249,7 @@ inline void COO_Matrix::setValues(std::vector r, std::vec } } //push back rest that was not found sorted - for (Intdx i = a_iter; i < static_cast(r.size()); i++) + for (IdxT i = a_iter; i < static_cast(r.size()); i++) { this->row_indices_.push_back(r[i]); this->column_indices_.push_back(c[i]); @@ -266,14 +266,14 @@ inline void COO_Matrix::setValues(std::vector r, std::vec * @brief Implements axpy this += alpha * a. Will sort before running * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] alpha matrix to be added * @param[in] a scalar to multiply by * * @post this = this + alpha * a */ -template -inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix& a) +template +inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix& a) { if (alpha == 0) { @@ -288,9 +288,9 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix&, std::vector&, std::vector&> tpm = a.getEntries(); + IdxT m = 0; + IdxT n = 0; + std::tuple&, std::vector&, std::vector&> tpm = a.getEntries(); const auto& [r, c, val] = tpm; std::tie(m,n) = a.getDimensions(); @@ -298,12 +298,12 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrixrows_size_ = this->rows_size_ > m ? this->rows_size_ : m; this->columns_size_ = this->columns_size_ > n ? this->columns_size_ : n; - Intdx a_iter = 0; + IdxT a_iter = 0; //iterate for all current values in matrix - for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) + for (IdxT i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values when they are not in current matrix - while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) + while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) { this->row_indices_.push_back(r[a_iter]); this->column_indices_.push_back(c[a_iter]); @@ -312,7 +312,7 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_MatrixcheckIncreaseSize(r[a_iter], c[a_iter]); a_iter++; } - if (a_iter >= static_cast(r.size())) + if (a_iter >= static_cast(r.size())) { break; } @@ -325,7 +325,7 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix(r.size()); i++) + for (IdxT i = a_iter; i < static_cast(r.size()); i++) { this->row_indices_.push_back(r[i]); this->column_indices_.push_back(c[i]); @@ -341,7 +341,7 @@ inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix::axpy(ScalarT alpha, COO_Matrix -inline void COO_Matrix::axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v) +template +inline void COO_Matrix::axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v) { if (alpha == 0) return; @@ -365,12 +365,12 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r //sort input this->sortSparseCOO(r, c, v); - Intdx a_iter = 0; + IdxT a_iter = 0; //iterate for all current values_ in matrix - for (Intdx i = 0; i < static_cast(this->row_indices_.size()); i++) + for (IdxT i = 0; i < static_cast(this->row_indices_.size()); i++) { //pushback values_ when they are not in current matrix - while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) + while(a_iter < static_cast(r.size()) && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) { this->row_indices_.push_back(r[a_iter]); this->column_indices_.push_back(c[a_iter]); @@ -379,7 +379,7 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r this->checkIncreaseSize(r[a_iter], c[a_iter]); a_iter++; } - if (a_iter >= static_cast(r.size())) + if (a_iter >= static_cast(r.size())) { break; } @@ -392,7 +392,7 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r } } //push back rest that was not found sorted_ - for (Intdx i = a_iter; i < static_cast(r.size()); i++) + for (IdxT i = a_iter; i < static_cast(r.size()); i++) { this->row_indices_.push_back(r[i]); this->column_indices_.push_back(c[i]); @@ -408,11 +408,11 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r * @brief Scale all values by alpha * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] alpha scalar to scale by */ -template -inline void COO_Matrix::scal(ScalarT alpha) +template +inline void COO_Matrix::scal(ScalarT alpha) { for (auto i = this->values_.begin(); i < this->values_.end(); i++) { @@ -424,11 +424,11 @@ inline void COO_Matrix::scal(ScalarT alpha) * @brief Calculates the Frobenius Norm of the matrix * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @return ScalarT - Frobenius Norm of the matrix */ -template -inline ScalarT COO_Matrix::frobNorm() +template +inline ScalarT COO_Matrix::frobNorm() { ScalarT totsum = 0.0; for (auto i = this->values_.begin(); i < this->values_.end(); i++) @@ -442,7 +442,7 @@ inline ScalarT COO_Matrix::frobNorm() * @brief Permutate the matrix to a different one. Only changes the coordinates * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] row_perm * @param[out] col_perm * @@ -450,8 +450,8 @@ inline ScalarT COO_Matrix::frobNorm() * * @post this = this(row_perm, col_perm) */ -template -inline void COO_Matrix::permutation(std::vector row_perm, std::vector col_perm) +template +inline void COO_Matrix::permutation(std::vector row_perm, std::vector col_perm) { assert(row_perm.size() = this->rows_size_); assert(col_perm.size() = this->columns_size_); @@ -469,7 +469,7 @@ inline void COO_Matrix::permutation(std::vector row_perm, * @brief Permutes the matrix and can change its size efficently * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] row_perm row permutation * @param[in] col_perm column permutation * @param[in] m number of rows @@ -481,8 +481,8 @@ inline void COO_Matrix::permutation(std::vector row_perm, * * @post this = this(row_perm, col_perm) and removed indices have corresponding values set to 0 */ -template -inline void COO_Matrix::permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n) +template +inline void COO_Matrix::permutationSizeMap(std::vector row_perm, std::vector col_perm, IdxT m, IdxT n) { assert(row_perm.size() == this->rows_size_); assert(col_perm.size() == this->columns_size_); @@ -509,11 +509,11 @@ inline void COO_Matrix::permutationSizeMap(std::vector ro * @brief Turn matrix into the zero matrix. Does not actually delete memory * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * */ -template -inline void COO_Matrix::zeroMatrix() +template +inline void COO_Matrix::zeroMatrix() { //resize doesn't effect capacity if smaller this->column_indices_.resize(0); @@ -526,7 +526,7 @@ inline void COO_Matrix::zeroMatrix() * @brief Turn matrix into the identity matrix * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * * @param[in] n size of the identity matrix * @@ -535,12 +535,12 @@ inline void COO_Matrix::zeroMatrix() * @todo - it might be better to explicitly zero out the matrix and require to be so in preconditions */ -template -inline void COO_Matrix::identityMatrix(Intdx n) +template +inline void COO_Matrix::identityMatrix(IdxT n) { //Reset Matrix this->zeroMatrix(); - for (Intdx i = 0; i < n; i++) + for (IdxT i = 0; i < n; i++) { this->column_indices_[i] = i; this->row_indices_[i] = i; @@ -553,10 +553,10 @@ inline void COO_Matrix::identityMatrix(Intdx n) * @brief Restructure the sparse matrix for faster accesses and modifications * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT */ -template -inline void COO_Matrix::sortSparse() +template +inline void COO_Matrix::sortSparse() { this->sortSparseCOO(this->row_indices_, this->column_indices_, this->values_); this->sorted_ = true; @@ -566,13 +566,13 @@ inline void COO_Matrix::sortSparse() * @brief Check if the matrix is sorted * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * * @param[out] bool - true if sorted, false otherwise */ -template -inline bool COO_Matrix::isSorted() +template +inline bool COO_Matrix::isSorted() { return this->sorted_; } @@ -581,30 +581,30 @@ inline bool COO_Matrix::isSorted() * @brief Get the number of non-zero elements in the matrix * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * - * @param[out] Intdx - number of non-zero elements in the matrix + * @param[out] IdxT - number of non-zero elements in the matrix */ -template -inline Intdx COO_Matrix::nnz() +template +inline IdxT COO_Matrix::nnz() { - return static_cast(this->values_.size); + return static_cast(this->values_.size); } -template -inline std::tuple COO_Matrix::getDimensions() +template +inline std::tuple COO_Matrix::getDimensions() { - return std::tuple(this->rows_size_, this->columns_size_); + return std::tuple(this->rows_size_, this->columns_size_); } /** * @brief Print matrix in sorted order * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT */ -template -inline void COO_Matrix::printMatrix() +template +inline void COO_Matrix::printMatrix() { if (this->sorted_ == false) { @@ -627,21 +627,21 @@ inline void COO_Matrix::printMatrix() * * Assumes rows and columns are sorted * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * * @param[in] rows - row indices * @param[in] r - row index * - * @return Intdx - index of lowest row + * @return IdxT - index of lowest row */ -template -inline Intdx COO_Matrix::indexStartRow(const std::vector &rows, Intdx r) +template +inline IdxT COO_Matrix::indexStartRow(const std::vector &rows, IdxT r) { //Specialized Binary Search for Lowest Row - Intdx i1 = 0; - Intdx i2 = rows->size()-1; - Intdx m_smallest = -1; - Intdx m = -1; + IdxT i1 = 0; + IdxT i2 = rows->size()-1; + IdxT m_smallest = -1; + IdxT m = -1; while (i1 <= i2) { m = (i2 + i1) / 2; @@ -673,21 +673,21 @@ inline Intdx COO_Matrix::indexStartRow(const std::vector * @brief Basic binary search * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] rows - row indices * @param[in] columns - column indices * @param[in] ri - row index * @param[in] ci - column index - * @return Intdx - returns the index of the coordinate + * @return IdxT - returns the index of the coordinate */ -template -inline Intdx COO_Matrix::sparseCordBinarySearch(const std::vector &rows, const std::vector &columns, Intdx ri, Intdx ci) +template +inline IdxT COO_Matrix::sparseCordBinarySearch(const std::vector &rows, const std::vector &columns, IdxT ri, IdxT ci) { assert(rows.size() == columns.size()); //basic binary search - Intdx i1 = 0; - Intdx i2 = rows.size()-1; - Intdx m = 0; + IdxT i1 = 0; + IdxT i2 = rows.size()-1; + IdxT m = 0; while (i1 <= i2) { m = (i2 + i1) / 2; @@ -721,14 +721,14 @@ inline Intdx COO_Matrix::sparseCordBinarySearch(const std::vecto * @brief Check if the size of the matrix needs to be increased * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param[in] r row index * @param[in] c column index * @return true if size was increased */ -template -inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) +template +inline bool COO_Matrix::checkIncreaseSize(IdxT r, IdxT c) { bool changed = false; if (r + 1 > this->rows_size_) @@ -758,13 +758,13 @@ inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) * @todo simple setup. Should add stable sorting since lists are pre-sorted_ * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param rows * @param columns * @param values */ -template -inline void COO_Matrix::sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &values) +template +inline void COO_Matrix::sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &values) { //index based sort code @@ -804,7 +804,7 @@ inline void COO_Matrix::sortSparseCOO(std::vector &rows, * @brief Constructor for COO Matrix with given cooridnates and values * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * * @param[in] r row indices * @param[in] c column indices @@ -818,8 +818,8 @@ inline void COO_Matrix::sortSparseCOO(std::vector &rows, * @post COO_Matrix is created with given coordinates and values */ -template -inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n) +template +inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, IdxT m, IdxT n) { this->values_ = v; this->row_indices_ = r; @@ -833,7 +833,7 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vector< * @brief Constructor for empty COO Matrix of a given size * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * * @param[in] m number of rows * @param[in] n number of columns @@ -841,14 +841,14 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vector< * @post empty COO Matrix is created with given size */ -template -inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) +template +inline COO_Matrix::COO_Matrix(IdxT m, IdxT n) { this->rows_size_ = m; this->columns_size_ = n; this->values_ = std::vector(); - this->row_indices_ = std::vector(); - this->column_indices_ = std::vector(); + this->row_indices_ = std::vector(); + this->column_indices_ = std::vector(); this->sorted_ = true; // An empty matrix is a sorted matrix. } @@ -856,24 +856,24 @@ inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) * @brief Constructor for empty COO Matrix of size 0 * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * * @post empty COO Matrix of size 0 is created */ -template -inline COO_Matrix::COO_Matrix() +template +inline COO_Matrix::COO_Matrix() { this->rows_size_ = 0; this->columns_size_ = 0; this->values_ = std::vector(); - this->row_indices_ = std::vector(); - this->column_indices_ = std::vector(); + this->row_indices_ = std::vector(); + this->column_indices_ = std::vector(); this->sorted_ = true; // An empty matrix is a sorted matrix. } -template -COO_Matrix::~COO_Matrix() +template +COO_Matrix::~COO_Matrix() { } From c9bd81cf47ee0c18005bfa4f485f4c311222fd8a Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 10:17:33 -0500 Subject: [PATCH 10/18] fixed sorting inconsistency --- src/LinearAlgebra/COO_Matrix.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 95f2ce03..6d611d86 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -826,7 +826,7 @@ inline COO_Matrix::COO_Matrix(std::vector r, std::vectorcolumn_indices_ = c; this->rows_size_ = m; this->columns_size_ = n; - this->sorted_ = false; // An empty matrix is a sorted matrix. + this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. } /** @@ -849,7 +849,7 @@ inline COO_Matrix::COO_Matrix(IdxT m, IdxT n) this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = true; // An empty matrix is a sorted matrix. + this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. } /** @@ -869,7 +869,7 @@ inline COO_Matrix::COO_Matrix() this->values_ = std::vector(); this->row_indices_ = std::vector(); this->column_indices_ = std::vector(); - this->sorted_ = true; // An empty matrix is a sorted matrix. + this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. } template From 171cd648eace3d4a0f72e7f4df40e2a07e83439c Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 13:12:49 -0500 Subject: [PATCH 11/18] updated synchronous machine names --- .../SynchronousMachine/SynchronousMachine.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index feb1952c..0ae8ad4e 100644 --- a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp @@ -20,17 +20,17 @@ namespace GridKit { * * @param[in] id - unique identifier for the component * @param[in] Lls - stator leakage inductance - * @param[in] Llkq - tuple of rotor leakage inductances - * @param[in] Llfd - field leakage inductance - * @param[in] Llkd - damper leakage inductance - * @param[in] Lmq - quadrature axis magnetizing inductance - * @param[in] Lmd - direct axis magnetizing inductance + * @param[in] Llkq - tuple of damper leakage reactances + * @param[in] Llfd - field leakage reactance + * @param[in] Llkd - damper leakage reactance + * @param[in] Lmq - quadrature axis magnetizing reactance + * @param[in] Lmd - direct axis magnetizing reactance * @param[in] Rs - stator resistance - * @param[in] Rkq - tuple of rotor resistances + * @param[in] Rkq - tuple of damper resistances * @param[in] Rfd - field resistance * @param[in] Rkd - damper resistance * @param[in] RJ - rotor moment of inertia - * @param[in] P - number of pole pairs + * @param[in] P - number of poles * @param[in] mub - magnetic saturation coefficient */ From 49fee8fd3df5d46299c03c92a8b4323f19b81be9 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 13:13:40 -0500 Subject: [PATCH 12/18] updated synchronous machine names --- .../PowerElectronics/SynchronousMachine/SynchronousMachine.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index 0ae8ad4e..960113b6 100644 --- a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp @@ -31,7 +31,7 @@ namespace GridKit { * @param[in] Rkd - damper resistance * @param[in] RJ - rotor moment of inertia * @param[in] P - number of poles - * @param[in] mub - magnetic saturation coefficient + * @param[in] mub - rated frequency */ template From bf69b7061050ba0d4019c8706c073b8b5f47371a Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 13:22:53 -0500 Subject: [PATCH 13/18] fixed notation and updated documentation for InductionMotor and LinearTransformer --- examples/Grid3Bus/README.md | 6 ++--- .../InductionMotor/InductionMotor.cpp | 10 +++++++-- .../InductionMotor/InductionMotor.hpp | 2 +- .../LinearTransformer/LinearTransformer.cpp | 22 ++++++++++++++----- .../LinearTransformer/LinearTransformer.hpp | 6 ++--- .../SynchronousMachine/SynchronousMachine.cpp | 4 ++-- .../SynchronousMachine/SynchronousMachine.hpp | 2 +- 7 files changed, 34 insertions(+), 18 deletions(-) diff --git a/examples/Grid3Bus/README.md b/examples/Grid3Bus/README.md index b420c373..329a66ab 100644 --- a/examples/Grid3Bus/README.md +++ b/examples/Grid3Bus/README.md @@ -22,13 +22,13 @@ Problem variables are voltage magnitudes and phases; they are stored in bus obje **Bus 1**: Slack bus, does not store variables no residuals. Voltage and phase are set to $`V_1 \equiv 1`$p.u. and $`\theta_1 \equiv 0`$, respectively. -**Bus 2**: PQ bus, stores variables $`V_2, \theta_2`$ and residuals $`P_2, Q_2`$. Load $`P_{L2} = 2.5`$p.u., $`Q_{L2} = -j0.8`$p.u. is attached to it. From the equations for [branch](../../src/Model/PowerFlow/Branch/README.md) and [load](../../src/Model/PowerFlow/Load/README.md) components, we assemble Bus 2 residuals as: +**Bus 2**: PQ bus, stores variables $`V_2, \theta_2`$ and residuals $`P_2, Q_2`$. Load $`P_{L1} = 2.5`$p.u., $`Q_{L1} = -j0.8`$p.u. is attached to it. From the equations for [branch](../../src/Model/PowerFlow/Branch/README.md) and [load](../../src/Model/PowerFlow/Load/README.md) components, we assemble Bus 2 residuals as: ```math \begin{array}{rcll} -P_2 & = &-P_{L2} &~~~\mathrm{(load ~2)} \\ +P_2 & = &-P_{L1} &~~~\mathrm{(load ~2)} \\ &&+ b_{12}|V_1||V_2|\sin(\theta_2-\theta_1) &~~~\mathrm{(branch ~12)} \\ &&+ b_{23}|V_2||V_3|\sin(\theta_2-\theta_3) &~~~\mathrm{(branch ~23)} \\ -Q_2 & = & -Q_{L2} &~~~\mathrm{(load ~2)} \\ +Q_2 & = & -Q_{L1} &~~~\mathrm{(load ~2)} \\ &&+ b_{12}|V_2|^2 - b_{12}|V_1||V_2|\cos(\theta_2-\theta_1) &~~~\mathrm{(branch ~12)} \\ &&+ b_{23}|V_2|^2 - b_{23}|V_2||V_3|\cos(\theta_2-\theta_3) &~~~\mathrm{(branch ~23)} \end{array} diff --git a/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp b/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp index 1ec8f0ca..dd0493a1 100644 --- a/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp +++ b/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp @@ -17,6 +17,12 @@ namespace GridKit { * Calls default ModelEvaluatorImpl constructor. * @todo create a test case utilizing the component. * @todo create a unit test to check correctness of component + * + * @tparam ScalarT - floating point type for the model + * @tparam IdxT - integer index type for the model + * + * @param[in] id - unique identifier for the component + * @param[in] Lls - stator leakage inductance */ template @@ -26,7 +32,7 @@ InductionMotor::InductionMotor(IdxT id, ScalarT Lls, ScalarT Rs, Llr_(Llr), Rr_(Rr), Lms_(Lms), - Rjac_(RJ), + RJ_(RJ), P_(P) { size_ = 10; @@ -82,7 +88,7 @@ int InductionMotor::evaluateResidual() f_[0] = y_[5] + y_[7]; f_[1] = (-1.0/2.0) * y_[5] - (sqrt(3.0)/2.0)*y_[6] + y_[7]; f_[2] = (-1.0/2.0) * y_[5] + (sqrt(3.0)/2.0)*y_[6] + y_[7]; - f_[3] = Rjac_ * yp_[3] - (3.0/4.0)*P_*Lms_ * (y_[5]*y_[9] - y_[6]*y_[8]); + f_[3] = RJ_ * yp_[3] - (3.0/4.0)*P_*Lms_ * (y_[5]*y_[9] - y_[6]*y_[8]); f_[4] = yp_[4] - y_[3]; f_[5] = (1.0/3.0)*(2.0* y_[0] - y_[1] - y_[2]) - Rs_*y_[5] - (Lls_ + Lms_) * yp_[5] - Lms_ * yp_[6]; f_[6] = (1.0/sqrt(3.0))*(-y_[1] + y_[2]) - Rs_*y_[6] - (Lls_ + Lms_) * yp_[6] - Lms_ * yp_[5]; diff --git a/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp b/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp index 76fd1513..97ef633e 100644 --- a/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp +++ b/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp @@ -66,7 +66,7 @@ namespace GridKit ScalarT Llr_; ScalarT Rr_; ScalarT Lms_; - ScalarT Rjac_; + ScalarT RJ_; ScalarT P_; }; } diff --git a/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp b/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp index 84a7c5a7..9e42482e 100644 --- a/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp +++ b/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp @@ -14,14 +14,24 @@ namespace GridKit { * Calls default ModelEvaluatorImpl constructor. * @todo Not tested in any model yet. Should be * @todo Has not been tested for correctness + * + * @tparam ScalarT - floating point type for the model + * @tparam IdxT - integer index type for the model + * + * @param[in] id - unique identifier for the component + * @param[in] L0 - inductance 0 + * @param[in] L1 - inductance 1 + * @param[in] R0 - resistance 0 + * @param[in] R1 - resistance 1 + * @param[in] M - mutual inductance */ template -LinearTransformer::LinearTransformer(IdxT id, ScalarT L1, ScalarT L2, ScalarT R1, ScalarT R2, ScalarT M) - : L1_(L1), - L2_(L2), +LinearTransformer::LinearTransformer(IdxT id, ScalarT L0, ScalarT L1, ScalarT R0, ScalarT R1, ScalarT M) + : L0_(L0), + L1_(L1), + R0_(R0), R1_(R1), - R2_(R2), M_(M) { size_ = 4; @@ -75,8 +85,8 @@ int LinearTransformer::evaluateResidual() { f_[0] = y_[2]; f_[1] = y_[3]; - f_[2] = y_[0] - R1_ * y_[2] - L1_ * yp_[2] - M_ * yp_[3]; - f_[2] = y_[1] - R2_ * y_[3] - M_ * yp_[2] - L2_ * yp_[3]; + f_[2] = y_[0] - R0_ * y_[2] - L0_ * yp_[2] - M_ * yp_[3]; + f_[2] = y_[1] - R1_ * y_[3] - M_ * yp_[2] - L1_ * yp_[3]; return 0; } diff --git a/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp b/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp index c2e18248..b09373d0 100644 --- a/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp +++ b/src/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp @@ -45,7 +45,7 @@ namespace GridKit using CircuitComponent::n_intern_; public: - LinearTransformer(IdxT id, ScalarT L1, ScalarT L2, ScalarT R1, ScalarT R2, ScalarT M); + LinearTransformer(IdxT id, ScalarT L0, ScalarT L1, ScalarT R0, ScalarT R1, ScalarT M); virtual ~LinearTransformer(); int allocate(); @@ -61,10 +61,10 @@ namespace GridKit int evaluateAdjointIntegrand(); private: + ScalarT L0_; ScalarT L1_; - ScalarT L2_; + ScalarT R0_; ScalarT R1_; - ScalarT R2_; ScalarT M_; }; } diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index 960113b6..6a3e6b52 100644 --- a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp @@ -46,7 +46,7 @@ SynchronousMachine::SynchronousMachine(IdxT id, ScalarT Lls, std: Rkq_(Rkq), Rfd_(Rfd), Rkd_(Rkd), - Rjac_(RJ), + RJ_(RJ), P_(P), mub_(mub) { @@ -115,7 +115,7 @@ int SynchronousMachine::evaluateResidual() f_[0] = y_[6]*cos1 + y_[7]*sin1 + y_[8]; f_[1] = y_[6]*cos23m + y_[7]*sin23m + y_[8]; f_[2] = y_[6]*cos23p + y_[7]*sin23p + y_[8]; - f_[3] = Rjac_ * yp_[4] - (3.0/4.0)*P_*(Lmd_ *y_[6]* (y_[7] + y_[11] + y_[12]) - Lmq_*y_[7]*(y_[6] + y_[9] + y_[0])); + f_[3] = RJ_ * yp_[4] - (3.0/4.0)*P_*(Lmd_ *y_[6]* (y_[7] + y_[11] + y_[12]) - Lmq_*y_[7]*(y_[6] + y_[9] + y_[0])); f_[4] = yp_[5] - y_[4]; f_[5] = (-2.0/3.0)*(y_[0]*cos1 +y_[1]*cos23m + y_[2]*cos23p) + Rs_*y_[6] + (Lls_ + Lmq_)*yp_[6] + Lmq_*yp_[9] + Lmq_*yp_[10] + y_[4]*(P_/2.0)*((Lls_ + Lmd_)*y_[7] + Lmd_*y_[11] + Lmd_*y_[12]); f_[6] = (-2.0/3.0)*(y_[0]*sin1 -y_[1]*sin23m - y_[2]*sin23p) + Rs_*y_[7] + (Lls_ + Lmd_)*yp_[7] + Lmd_*yp_[11] + Lmd_*yp_[12] - y_[4]*(P_/2.0)*((Lls_ + Lmq_)*y_[6] + Lmq_*y_[9] + Lmq_*y_[10]); diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp index 5bf88b6c..ec1d7486 100644 --- a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp +++ b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp @@ -72,7 +72,7 @@ namespace GridKit std::tuple Rkq_; ScalarT Rfd_; ScalarT Rkd_; - ScalarT Rjac_; + ScalarT RJ_; ScalarT P_; ScalarT mub_; }; From 5ade72f5fbee9d4cfc870e5c848b43f5efa66c19 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 13:25:07 -0500 Subject: [PATCH 14/18] fixed guard for hpp file --- src/LinearAlgebra/COO_Matrix.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/COO_Matrix.hpp index 6d611d86..b84e239e 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -1,4 +1,5 @@ -#pragma once +#ifndef COO_MATRIX_HPP +#define COO_MATRIX_HPP #include #include @@ -877,3 +878,5 @@ COO_Matrix::~COO_Matrix() { } + +#endif \ No newline at end of file From 01b87ab2fa743fa815ae2408bb3f10252ad7bc27 Mon Sep 17 00:00:00 2001 From: Shaked Regev Date: Fri, 14 Feb 2025 15:13:10 -0500 Subject: [PATCH 15/18] straggling typo --- examples/Microgrid/Microgrid.cpp | 2 +- examples/RLCircuit/RLCircuit.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/Microgrid/Microgrid.cpp b/examples/Microgrid/Microgrid.cpp index e6254720..75c70fdd 100644 --- a/examples/Microgrid/Microgrid.cpp +++ b/examples/Microgrid/Microgrid.cpp @@ -27,7 +27,7 @@ int main(int argc, char const *argv[]) bool use_jac = true; //Create model - ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(rel_tol, abs_tol, use_jac, max_step_amount); + GridKit::PowerElectronicsModel* sysmodel = new GridKit::PowerElectronicsModel(rel_tol, abs_tol, use_jac, max_step_amount); //Modeled after the problem in the paper double RN = 1.0e4; diff --git a/examples/RLCircuit/RLCircuit.cpp b/examples/RLCircuit/RLCircuit.cpp index bfe54e82..b6c2c7da 100644 --- a/examples/RLCircuit/RLCircuit.cpp +++ b/examples/RLCircuit/RLCircuit.cpp @@ -24,7 +24,7 @@ int main(int argc, char const *argv[]) //TODO:setup as named parameters //Create circuit model - auto* sysmodel = new ModelLib::PowerElectronicsModel(rel_tol, abs_tol, use_jac); + auto* sysmodel = new GridKit::PowerElectronicsModel(rel_tol, abs_tol, use_jac); size_t idoff = 0; From 3e9b2dd6f327f4c95b7d6262161080ab3361c799 Mon Sep 17 00:00:00 2001 From: pelesh Date: Mon, 17 Feb 2025 17:41:30 -0500 Subject: [PATCH 16/18] [skip ci] Update comment in InductionMotor.cpp --- src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp b/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp index dd0493a1..b5791014 100644 --- a/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp +++ b/src/Model/PowerElectronics/InductionMotor/InductionMotor.cpp @@ -18,7 +18,7 @@ namespace GridKit { * @todo create a test case utilizing the component. * @todo create a unit test to check correctness of component * - * @tparam ScalarT - floating point type for the model + * @tparam ScalarT - data type for scalar variables in the model * @tparam IdxT - integer index type for the model * * @param[in] id - unique identifier for the component From 7a0081f393a435dd638da43f0307c2c49eba3550 Mon Sep 17 00:00:00 2001 From: pelesh Date: Mon, 17 Feb 2025 17:42:15 -0500 Subject: [PATCH 17/18] [skip ci]Update equation in Grid3Bus/README.md --- examples/Grid3Bus/README.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/Grid3Bus/README.md b/examples/Grid3Bus/README.md index 329a66ab..19e23325 100644 --- a/examples/Grid3Bus/README.md +++ b/examples/Grid3Bus/README.md @@ -57,8 +57,7 @@ with variables $`\theta_2, |V_2|`$ and $`\theta_3`$. Nonlinear solver can approximate Jacobian numerically, however this is computationaly expensive and scales poorly with the size of the problem. Typically, one needs to provide Jacobian in addition to residual to the nonlinear solver. For nonlinear problem defined by (vector) function $`\mathbf{f}(\mathbf{x})=0`$, Jacobian matrix is defined as ```math -jac_{i,j}=\frac{\partial f_i}{\partial x_j}, ~~~ i,j=1,\ldots,N -``` +J_{i,j}=\frac{\partial f_i}{\partial x_j}, ~~~ i,j=1,\ldots,N where $`N`$ is the size of vectors $`\mathbf{f}`$ and $`\mathbf{x}`$. Jacobian for our model is evaluated as ```math \mathbf{J} = From 2af2f9e620073291a05d44040b54885d7caf6c07 Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Mon, 17 Feb 2025 17:49:54 -0500 Subject: [PATCH 18/18] [skip ci] Fix latex error. --- examples/Grid3Bus/README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/Grid3Bus/README.md b/examples/Grid3Bus/README.md index 19e23325..31bd6f9b 100644 --- a/examples/Grid3Bus/README.md +++ b/examples/Grid3Bus/README.md @@ -58,6 +58,7 @@ with variables $`\theta_2, |V_2|`$ and $`\theta_3`$. Nonlinear solver can approximate Jacobian numerically, however this is computationaly expensive and scales poorly with the size of the problem. Typically, one needs to provide Jacobian in addition to residual to the nonlinear solver. For nonlinear problem defined by (vector) function $`\mathbf{f}(\mathbf{x})=0`$, Jacobian matrix is defined as ```math J_{i,j}=\frac{\partial f_i}{\partial x_j}, ~~~ i,j=1,\ldots,N +``` where $`N`$ is the size of vectors $`\mathbf{f}`$ and $`\mathbf{x}`$. Jacobian for our model is evaluated as ```math \mathbf{J} =