diff --git a/examples/Grid3Bus/README.md b/examples/Grid3Bus/README.md index b628ad42..31bd6f9b 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/examples/Microgrid/Microgrid.cpp b/examples/Microgrid/Microgrid.cpp index 901433b1..75c70fdd 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 usejac = true; + bool use_jac = true; //Create model - GridKit::PowerElectronicsModel* sysmodel = new GridKit::PowerElectronicsModel(reltol, abstol, usejac, 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; @@ -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 46afe872..b6c2c7da 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; - bool usejac = true; + double abs_tol = 1.0e-8; + double rel_tol = 1.0e-8; + bool use_jac = true; //TODO:setup as named parameters //Create circuit model - GridKit::PowerElectronicsModel* sysmodel = new GridKit::PowerElectronicsModel(reltol, abstol, usejac); + auto* sysmodel = new GridKit::PowerElectronicsModel(rel_tol, abs_tol, use_jac); size_t idoff = 0; diff --git a/examples/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/ScaleMicrogrid/ScaleMicrogrid.cpp index a145c0be..496eb8c0 100644 --- a/examples/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -66,18 +66,18 @@ 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; - 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, - usejac, + auto* sysmodel = new PowerElectronicsModel(rel_tol, + abs_tol, + use_jac, SCALE_MICROGRID_MAX_STEPS); const std::vector* true_vec = &answer_key_N8; 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/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 9036085a..b84e239e 100644 --- a/src/LinearAlgebra/COO_Matrix.hpp +++ b/src/LinearAlgebra/COO_Matrix.hpp @@ -1,4 +1,5 @@ - +#ifndef COO_MATRIX_HPP +#define COO_MATRIX_HPP #include #include @@ -12,24 +13,24 @@ /** * @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 */ -template +template class COO_Matrix { private: std::vector values_; - std::vector row_indexes_; - std::vector column_indexes_; - 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,138 +38,138 @@ 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> getEntrieCopies(); - std::tuple, std::vector, std::vector> getEntrieCopiesSubMatrix(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> getDataToCSR(); - 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(); + ScalarT frobNorm(); // --- Permutation Operations --- - //No sorting is actually done. Only done when nesscary - void permutation(std::vector row_perm, std::vector col_perm); - void permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n); + //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, 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 &vals); + 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); }; /** - * @brief Get copy of row values_ + * @brief Get copy of row index * * @tparam ScalarT - * @tparam Intdx - * @param r - * @return std::tuple, std::vector> + * @tparam IdxT + * @param[in] r row index + * @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 rowindex = this->indexStartRow(r); + IdxT row_index = this->indexStartRow(r); - if (rowindex == -1) + if (row_index == -1) { - return {std::vector(),std::vector()}; + return {std::vector(),std::vector()}; } - Intdx rsize = rowindex; + IdxT rsize = row_index; do { rsize++; - } while (rsize < this->values_.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() + row_index, this->column_indices_.begin() + rsize},{this->values_.begin() + row_index, 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 - * @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::getEntries() +template +inline std::tuple&, std::vector&, std::vector&> COO_Matrix::getEntries() { if (!this->sorted_) { 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> + * @tparam IdxT + * @return std::tuple, std::vector, std::vector> */ -template -inline std::tuple, std::vector, std::vector> COO_Matrix::getEntrieCopies() +template +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 - * @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::getDataToCSR() +template +inline std::tuple, std::vector, std::vector> COO_Matrix::setDataToCSR() { if (!this->isSorted()) this->sortSparse(); - std::vector rowsizevec(this->rows_size_ + 1, 0); - Intdx counter = 0; - for (Intdx i = 0; i < static_cast(rowsizevec.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++) { - rowsizevec[i + 1] = rowsizevec[i]; - while (counter < static_cast(this->row_indexes_.size()) && i == this->row_indexes_[counter]) + 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_indexes_, this->values_}; + return {row_size_vec, this->column_indices_, this->values_}; } /** @@ -178,74 +179,81 @@ 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 rowsizevec(this->rows_size_ + 1, 0); - Intdx counter = 0; - for (Intdx i = 0; i < static_cast(rowsizevec.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++) { - rowsizevec[i + 1] = rowsizevec[i]; - while (counter < static_cast(this->row_indexes_.size()) && i == this->row_indexes_[counter]) + 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. + * + * Matrix entries will be sorted in row-major order before the method returns. * * @tparam ScalarT - * @tparam Intdx - * @param r - * @param c - * @param v + * @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 + * + * @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) +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 aiter = 0; + IdxT a_iter = 0; //iterate for all current values_ in matrix - for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + for (IdxT 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(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_indexes_.push_back(r[aiter]); - this->column_indexes_.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_indexes_[i] && c[aiter] == this->column_indexes_[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 (IdxT i = a_iter; 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]); @@ -256,15 +264,17 @@ 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 + * @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) { @@ -279,47 +289,47 @@ 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(); - //Increase size as nesscary + //Increase size as necessary this->rows_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++) + IdxT a_iter = 0; + //iterate for all current values in matrix + for (IdxT 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]))) + //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_indexes_.push_back(r[aiter]); - this->column_indexes_.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_indexes_[i] && c[aiter] == this->column_indexes_[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 (IdxT i = a_iter; 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 * val[i]); this->checkIncreaseSize(r[i], c[i]); @@ -329,17 +339,22 @@ 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) +template +inline void COO_Matrix::axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v) { if (alpha == 0) return; @@ -351,37 +366,37 @@ inline void COO_Matrix::axpy(ScalarT alpha, std::vector r //sort input this->sortSparseCOO(r, c, v); - Intdx aiter = 0; + IdxT a_iter = 0; //iterate for all current values_ in matrix - for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + for (IdxT 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(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_indexes_.push_back(r[aiter]); - this->column_indexes_.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_indexes_[i] && c[aiter] == this->column_indexes_[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 (IdxT i = a_iter; 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]); @@ -391,14 +406,14 @@ 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 + * @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++) { @@ -407,14 +422,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 + * @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++) @@ -428,38 +443,47 @@ inline ScalarT COO_Matrix::frobnorm() * @brief Permutate the matrix to a different one. Only changes the coordinates * * @tparam ScalarT - * @tparam Intdx - * @param row_perm - * @param col_perm + * @tparam IdxT + * @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) +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_); 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 * * @tparam ScalarT - * @tparam Intdx - * @param row_perm size of m - * @param col_perm size of n - * @param m - * @param n + * @tparam IdxT + * @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) +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_); @@ -469,44 +493,58 @@ 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 + * @tparam IdxT + * */ -template -inline void COO_Matrix::zeroMatrix() +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; } -template -inline void COO_Matrix::identityMatrix(Intdx n) +/** + * @brief Turn matrix into the identity matrix + * + * @tparam ScalarT + * @tparam IdxT + * + * @param[in] n size of the identity matrix + * + * @post this = I_n + * + * @todo - it might be better to explicitly zero out the matrix and require to be so in preconditions + */ + +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_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; @@ -516,41 +554,58 @@ 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_indexes_, this->column_indexes_, this->values_); + this->sortSparseCOO(this->row_indices_, this->column_indices_, this->values_); this->sorted_ = true; } -template -inline bool COO_Matrix::isSorted() +/** + * @brief Check if the matrix is sorted + * + * @tparam ScalarT + * @tparam IdxT + * + * @param[out] bool - true if sorted, false otherwise + */ + +template +inline bool COO_Matrix::isSorted() { return this->sorted_; } -template -inline Intdx COO_Matrix::nnz() +/** + * @brief Get the number of non-zero elements in the matrix + * + * @tparam ScalarT + * @tparam IdxT + * + * @param[out] IdxT - number of non-zero elements in the matrix + */ +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 that is sorted_ + * @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) { @@ -561,30 +616,33 @@ 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 - * @return Intdx + * @tparam IdxT + * + * @param[in] rows - row indices + * @param[in] r - row index + * + * @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; @@ -616,21 +674,21 @@ inline Intdx COO_Matrix::indexStartRow(const std::vector * @brief Basic binary search * * @tparam ScalarT - * @tparam Intdx - * @param rows - * @param columns - * @param ri - * @param ci - * @return Intdx + * @tparam IdxT + * @param[in] rows - row indices + * @param[in] columns - column indices + * @param[in] ri - row index + * @param[in] ci - column index + * @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; @@ -660,8 +718,18 @@ inline Intdx COO_Matrix::sparseCordBinarySearch(const std::vecto return m; } -template -inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) +/** + * @brief Check if the size of the matrix needs to be increased + * + * @tparam ScalarT + * @tparam IdxT + * @param[in] r row index + * @param[in] c column index + * @return true if size was increased + */ + +template +inline bool COO_Matrix::checkIncreaseSize(IdxT r, IdxT c) { bool changed = false; if (r + 1 > this->rows_size_) @@ -679,18 +747,25 @@ inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) } /** - * @brief Sort a disoreded 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. + * Duplicate entries are not allowed and should be pre-summed. * - * @todo simple setup. Should add stable sorting since list are pre-sorted_ + * @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_ * * @tparam ScalarT - * @tparam Intdx + * @tparam IdxT * @param rows * @param columns - * @param values_ + * @param values */ -template -inline void COO_Matrix::sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &vals) +template +inline void COO_Matrix::sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &values) { //index based sort code @@ -717,7 +792,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]]); @@ -726,41 +801,82 @@ inline void COO_Matrix::sortSparseCOO(std::vector &rows, } } -template -inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n) +/** + * @brief Constructor for COO Matrix with given cooridnates and values + * + * @tparam ScalarT + * @tparam IdxT + * + * @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, IdxT m, IdxT 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; + this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. } -template -inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) +/** + * @brief Constructor for empty COO Matrix of a given size + * + * @tparam ScalarT + * @tparam IdxT + * + * @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(IdxT m, IdxT n) { this->rows_size_ = m; this->columns_size_ = n; this->values_ = std::vector(); - this->row_indexes_ = std::vector(); - this->column_indexes_ = std::vector(); - this->sorted_ = false; + this->row_indices_ = std::vector(); + this->column_indices_ = std::vector(); + this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. } -template -inline COO_Matrix::COO_Matrix() +/** + * @brief Constructor for empty COO Matrix of size 0 + * + * @tparam ScalarT + * @tparam IdxT + * + * @post empty COO Matrix of size 0 is created + */ + +template +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->sorted_ = false; + this->row_indices_ = std::vector(); + this->column_indices_ = std::vector(); + this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. } -template -COO_Matrix::~COO_Matrix() +template +COO_Matrix::~COO_Matrix() { } + +#endif \ No newline at end of file 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/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/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..b5791014 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 - data type for scalar variables in 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 diff --git a/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp b/src/Model/PowerElectronics/InductionMotor/InductionMotor.hpp index 5a3002d5..97ef633e 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_; 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.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 f11f5e31..b09373d0 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_; @@ -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/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..6a3e6b52 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 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 damper resistances + * @param[in] Rfd - field resistance + * @param[in] Rkd - damper resistance + * @param[in] RJ - rotor moment of inertia + * @param[in] P - number of poles + * @param[in] mub - rated frequency */ template diff --git a/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp b/src/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp index 190ebad6..ec1d7486 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_; diff --git a/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/src/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 50128b93..29fa037b 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,32 +44,49 @@ 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 not use jacobian - usejac_ = false; + // By default don't use the jacobian + use_jac_ = false; } /** * @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 rt = 1e-4, double at = 1e-4, bool ju = false, IdxT msa = 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_ = rt; - atol_ = at; - this->max_steps_ = msa; - // Can choose if to use jacobain - usejac_ = ju; + rel_tol_ = rel_tol; + abs_tol_ = abs_tol; + this->max_steps_ = max_steps; + // Can choose if to use jacobian + use_jac_ = use_jac; } /** * @brief Destructor for the system model + * + * @pre System components are allocated + * + * @post System components are deallocated + * */ virtual ~PowerElectronicsModel() { @@ -91,15 +108,15 @@ 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 - * @return false + * @return true if all components have jacobian + * @return false otherwise */ bool hasJacobian() { - if (!this->usejac_) + if (!this->use_jac_) return false; for (const auto &component : components_) @@ -116,8 +133,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 +160,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 +178,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 +211,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 +244,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 @@ -234,7 +256,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; @@ -253,7 +275,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; @@ -325,7 +348,7 @@ namespace GridKit static constexpr IdxT neg1_ = static_cast(-1); std::vector components_; - bool usejac_; + bool use_jac_; }; // class PowerElectronicsModel 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; }