diff --git a/CircuitGraph.hpp b/CircuitGraph.hpp index 070d838..c0dde1a 100644 --- a/CircuitGraph.hpp +++ b/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/ComponentLib/PowerElectronics/Capacitor/Capacitor.cpp b/ComponentLib/PowerElectronics/Capacitor/Capacitor.cpp index 779b8a6..6e02cd2 100644 --- a/ComponentLib/PowerElectronics/Capacitor/Capacitor.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/Capacitor/Capacitor.hpp b/ComponentLib/PowerElectronics/Capacitor/Capacitor.hpp index 9a54b66..3b218a3 100644 --- a/ComponentLib/PowerElectronics/Capacitor/Capacitor.hpp +++ b/ComponentLib/PowerElectronics/Capacitor/Capacitor.hpp @@ -36,7 +36,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/ComponentLib/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index 00d174a..2cbb39d 100644 --- a/ComponentLib/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/ComponentLib/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp index 69f8e41..2277b2d 100644 --- a/ComponentLib/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/ComponentLib/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -56,7 +56,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/InductionMotor/InductionMotor.cpp b/ComponentLib/PowerElectronics/InductionMotor/InductionMotor.cpp index fa0da7f..1382ae5 100644 --- a/ComponentLib/PowerElectronics/InductionMotor/InductionMotor.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/InductionMotor/InductionMotor.hpp b/ComponentLib/PowerElectronics/InductionMotor/InductionMotor.hpp index 5072a02..94de912 100644 --- a/ComponentLib/PowerElectronics/InductionMotor/InductionMotor.hpp +++ b/ComponentLib/PowerElectronics/InductionMotor/InductionMotor.hpp @@ -36,7 +36,7 @@ namespace ModelLib 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 ModelLib ScalarT Llr_; ScalarT Rr_; ScalarT Lms_; - ScalarT RJ_; + ScalarT Rjac_; ScalarT P_; }; } diff --git a/ComponentLib/PowerElectronics/Inductor/Inductor.cpp b/ComponentLib/PowerElectronics/Inductor/Inductor.cpp index 2bb4dd1..b4b31da 100644 --- a/ComponentLib/PowerElectronics/Inductor/Inductor.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/Inductor/Inductor.hpp b/ComponentLib/PowerElectronics/Inductor/Inductor.hpp index e69e2aa..ba2a56f 100644 --- a/ComponentLib/PowerElectronics/Inductor/Inductor.hpp +++ b/ComponentLib/PowerElectronics/Inductor/Inductor.hpp @@ -36,7 +36,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/LinearTransformer/LinearTransformer.hpp b/ComponentLib/PowerElectronics/LinearTransformer/LinearTransformer.hpp index 487f2de..e531035 100644 --- a/ComponentLib/PowerElectronics/LinearTransformer/LinearTransformer.hpp +++ b/ComponentLib/PowerElectronics/LinearTransformer/LinearTransformer.hpp @@ -36,7 +36,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/ComponentLib/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp index f486e17..b1590c8 100644 --- a/ComponentLib/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/ComponentLib/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index 63d33ef..021e81d 100644 --- a/ComponentLib/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/ComponentLib/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -36,7 +36,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/ComponentLib/PowerElectronics/MicrogridLine/MicrogridLine.cpp index 3198af2..ea20d49 100644 --- a/ComponentLib/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/ComponentLib/PowerElectronics/MicrogridLine/MicrogridLine.hpp index ff25bd9..d4b17b8 100644 --- a/ComponentLib/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/ComponentLib/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -35,7 +35,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/ComponentLib/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index ddfa9ab..fd674c9 100644 --- a/ComponentLib/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/ComponentLib/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 1ecd18f..1d3a00a 100644 --- a/ComponentLib/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/ComponentLib/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -36,7 +36,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/Resistor/Resistor.cpp b/ComponentLib/PowerElectronics/Resistor/Resistor.cpp index caf0335..8c7453a 100644 --- a/ComponentLib/PowerElectronics/Resistor/Resistor.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/Resistor/Resistor.hpp b/ComponentLib/PowerElectronics/Resistor/Resistor.hpp index 64995cb..0c27bf9 100644 --- a/ComponentLib/PowerElectronics/Resistor/Resistor.hpp +++ b/ComponentLib/PowerElectronics/Resistor/Resistor.hpp @@ -36,7 +36,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/ComponentLib/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index e3f65be..a1e8957 100644 --- a/ComponentLib/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp b/ComponentLib/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp index 53cc60b..ad1730e 100644 --- a/ComponentLib/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp +++ b/ComponentLib/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp @@ -37,7 +37,7 @@ namespace ModelLib 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 ModelLib std::tuple Rkq_; ScalarT Rfd_; ScalarT Rkd_; - ScalarT RJ_; + ScalarT Rjac_; ScalarT P_; ScalarT mub_; }; diff --git a/ComponentLib/PowerElectronics/TransmissionLine/TransmissionLine.cpp b/ComponentLib/PowerElectronics/TransmissionLine/TransmissionLine.cpp index 010f38e..c937d61 100644 --- a/ComponentLib/PowerElectronics/TransmissionLine/TransmissionLine.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/TransmissionLine/TransmissionLine.hpp b/ComponentLib/PowerElectronics/TransmissionLine/TransmissionLine.hpp index c510d4c..3a42b62 100644 --- a/ComponentLib/PowerElectronics/TransmissionLine/TransmissionLine.hpp +++ b/ComponentLib/PowerElectronics/TransmissionLine/TransmissionLine.hpp @@ -40,7 +40,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerElectronics/VoltageSource/VoltageSource.cpp b/ComponentLib/PowerElectronics/VoltageSource/VoltageSource.cpp index 2e40415..0e5bfb1 100644 --- a/ComponentLib/PowerElectronics/VoltageSource/VoltageSource.cpp +++ b/ComponentLib/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/ComponentLib/PowerElectronics/VoltageSource/VoltageSource.hpp b/ComponentLib/PowerElectronics/VoltageSource/VoltageSource.hpp index e464feb..b06f2a2 100644 --- a/ComponentLib/PowerElectronics/VoltageSource/VoltageSource.hpp +++ b/ComponentLib/PowerElectronics/VoltageSource/VoltageSource.hpp @@ -36,7 +36,7 @@ namespace ModelLib using CircuitComponent::ypB_; using CircuitComponent::fB_; using CircuitComponent::gB_; - using CircuitComponent::J_; + using CircuitComponent::jac_; using CircuitComponent::param_; using CircuitComponent::idc_; diff --git a/ComponentLib/PowerFlow/Bus/BaseBus.hpp b/ComponentLib/PowerFlow/Bus/BaseBus.hpp index 4e8cbd5..9ef7780 100644 --- a/ComponentLib/PowerFlow/Bus/BaseBus.hpp +++ b/ComponentLib/PowerFlow/Bus/BaseBus.hpp @@ -84,8 +84,8 @@ namespace ModelLib 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/ComponentLib/PowerFlow/Bus/BusSlack.hpp b/ComponentLib/PowerFlow/Bus/BusSlack.hpp index 19c27d9..ddbd495 100644 --- a/ComponentLib/PowerFlow/Bus/BusSlack.hpp +++ b/ComponentLib/PowerFlow/Bus/BusSlack.hpp @@ -82,8 +82,8 @@ namespace ModelLib 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/ComponentLib/PowerFlow/MiniGrid/MiniGrid.cpp b/ComponentLib/PowerFlow/MiniGrid/MiniGrid.cpp index 1774b2c..83e2ca2 100644 --- a/ComponentLib/PowerFlow/MiniGrid/MiniGrid.cpp +++ b/ComponentLib/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/ComponentLib/PowerFlow/MiniGrid/MiniGrid.hpp b/ComponentLib/PowerFlow/MiniGrid/MiniGrid.hpp index 1ec9e43..01db2ee 100644 --- a/ComponentLib/PowerFlow/MiniGrid/MiniGrid.hpp +++ b/ComponentLib/PowerFlow/MiniGrid/MiniGrid.hpp @@ -75,8 +75,8 @@ namespace ModelLib 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/Examples/Grid3Bus/README.md b/Examples/Grid3Bus/README.md index 39d2774..d547b21 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 0483b2b..6ea533d 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 c563dbd..34930a5 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/ModelEvaluator.hpp b/ModelEvaluator.hpp index ad56fbd..4da3b8a 100644 --- a/ModelEvaluator.hpp +++ b/ModelEvaluator.hpp @@ -106,7 +106,7 @@ namespace ModelLib 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/ModelEvaluatorImpl.hpp b/ModelEvaluatorImpl.hpp index 0f8969f..afb8023 100644 --- a/ModelEvaluatorImpl.hpp +++ b/ModelEvaluatorImpl.hpp @@ -94,7 +94,7 @@ namespace ModelLib 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 ModelLib // 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 ModelLib COO_Matrix& getJacobian() { - return J_; + return jac_; } const COO_Matrix& getJacobian() const { - return J_; + return jac_; } std::vector& getIntegrand() @@ -299,7 +299,7 @@ namespace ModelLib std::vector fB_; std::vector gB_; - COO_Matrix J_; + COO_Matrix jac_; std::vector param_; std::vector param_up_; @@ -308,8 +308,8 @@ namespace ModelLib 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/PowerElectronicsModel.hpp b/PowerElectronicsModel.hpp index 539d497..57521f1 100644 --- a/PowerElectronicsModel.hpp +++ b/PowerElectronicsModel.hpp @@ -34,9 +34,9 @@ namespace ModelLib // 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 ModelLib 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 ModelLib /** * @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 ModelLib /** * @brief Destructor for the system model + * + * @pre System components are allocated + * + * @post System components are deallocated + * */ virtual ~PowerElectronicsModel() { @@ -94,8 +108,8 @@ namespace ModelLib * @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 ModelLib * @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 ModelLib /** * @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 ModelLib /** * @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 ModelLib /** * @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 ModelLib /** * @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 ModelLib } // 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/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index 2c7df1d..5c8dc2a 100644 --- a/Solver/Dynamic/Ida.cpp +++ b/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/SparseMatrix/COO_Matrix.hpp b/SparseMatrix/COO_Matrix.hpp index 646f4f2..62bcd66 100644 --- a/SparseMatrix/COO_Matrix.hpp +++ b/SparseMatrix/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/SystemModel.hpp b/SystemModel.hpp index 14b7e4e..6d438de 100644 --- a/SystemModel.hpp +++ b/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; } diff --git a/SystemSteadyStateModel.hpp b/SystemSteadyStateModel.hpp index ca1af16..b102b6c 100644 --- a/SystemSteadyStateModel.hpp +++ b/SystemSteadyStateModel.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)