diff --git a/CHANGELOG.md b/CHANGELOG.md index 8766c4b80..f30c515db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -55,6 +55,7 @@ - Added 200 Bus Synthetic Illinois Case - Added node objects to `PowerElectronics` module & updated all examples to make use of them. - Separated internal and external residuals of `PowerElectronics` models. +- Remove data copying between system and components in `PowerElectronics` models. ## v0.1 diff --git a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp index eb5f4fe68..f21ebaa6a 100644 --- a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp +++ b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp @@ -59,7 +59,7 @@ namespace GridKit template int Capacitor::evaluateInternalResidual() { - f_[2] = -C_ * yp_[2] + y_[0] - y_[1] - y_[2]; + f_int_[0] = -C_ * yp_int_[0] + y_[0] - y_[1] - y_int_[0]; return 0; } @@ -67,9 +67,9 @@ namespace GridKit int Capacitor::evaluateExternalResidual() { // input - f_[0] = C_ * yp_[2]; + f_[0] = C_ * yp_int_[0]; // output - f_[1] = -C_ * yp_[2]; + f_[1] = -C_ * yp_int_[0]; return 0; } diff --git a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp index fafdbf79a..a45b55535 100644 --- a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp +++ b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp @@ -26,9 +26,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index b787b48d8..a99257454 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -166,6 +166,21 @@ namespace GridKit */ virtual int evaluateExternalResidual() = 0; + void setInternalPointer(const ScalarT* internals) + { + y_int_ = internals; + } + + void setInternalDerivativePointer(const ScalarT* internals_p) + { + yp_int_ = internals_p; + } + + void setInternalResidualPointer(ScalarT* internal_res) + { + f_int_ = internal_res; + } + protected: /** * @brief Reset the Jacobian so it can be constructed. Helper method for \ref setJacValues(). @@ -412,6 +427,13 @@ namespace GridKit /// The number of non-zero elements currently inserted into the Jacobian. See \ref setJacValues() size_t current_jac_size_{0}; + /// @brief A pointer to the internal variables of this component. + const ScalarT* y_int_; + /// @brief A pointer to the internal derivatives of this component. + const ScalarT* yp_int_; + /// @brief A pointer to the internal residuals of this component + ScalarT* f_int_; + std::vector y_; std::vector yp_; std::vector tag_; diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index e77ebb5c7..1458c3db2 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -45,20 +45,13 @@ namespace GridKit assert(node_bus_->size() == 2); // internals [Pi, Qi, phi_di, phi_qi, gamma_di, gamma_qi, il_di, il_qi, vo_di, vo_qi, io_di, io_qi, \delta_i] // externals [\omega_ref, vba_out, vbb_out] - size_ = 16; + size_ = refframe_ ? 15 : 16; n_intern_ = refframe_ ? 12 : 13; - n_extern_ = refframe_ ? 4 : 3; + n_extern_ = 3; idc_ = id; - nnz_ = refframe_ ? 77 : 78; + nnz_ = refframe_ ? 73 : 78; - if (refframe_) - { - extern_indices_ = {0, 1, 2, 15}; - } - else - { - extern_indices_ = {0, 1, 2}; - } + extern_indices_ = {0, 1, 2}; } template @@ -91,48 +84,49 @@ namespace GridKit template int DistributedGenerator::evaluateInternalResidual() { - ScalarT omega = wb_ - mp_ * y_[3]; + ScalarT omega = wb_ - mp_ * y_int_[0]; + ScalarT delta = refframe_ ? ScalarT(0.0) : y_int_[12]; // Take incoming voltages to current rotator reference frame - ScalarT vbd_in = std::cos(y_[15]) * y_[1] + std::sin(y_[15]) * y_[2]; - ScalarT vbq_in = -std::sin(y_[15]) * y_[1] + std::cos(y_[15]) * y_[2]; + ScalarT vbd_in = std::cos(delta) * y_[1] + std::sin(delta) * y_[2]; + ScalarT vbq_in = -std::sin(delta) * y_[1] + std::cos(delta) * y_[2]; // ### Internal Componenets ## // P and Q equations - f_[3] = -yp_[3] + wc_ * (y_[11] * y_[13] + y_[12] * y_[14] - y_[3]); - f_[4] = -yp_[4] + wc_ * (-y_[11] * y_[14] + y_[12] * y_[13] - y_[4]); + f_int_[0] = -yp_int_[0] + wc_ * (y_int_[8] * y_int_[10] + y_int_[9] * y_int_[11] - y_int_[0]); + f_int_[1] = -yp_int_[1] + wc_ * (-y_int_[8] * y_int_[11] + y_int_[9] * y_int_[10] - y_int_[1]); // Voltage control - ScalarT vod_star = Vn_ - nq_ * y_[4]; + ScalarT vod_star = Vn_ - nq_ * y_int_[1]; ScalarT voq_star = static_cast(0.0); - f_[5] = -yp_[5] + vod_star - y_[11]; - f_[6] = -yp_[6] + voq_star - y_[12]; + f_int_[2] = -yp_int_[2] + vod_star - y_int_[8]; + f_int_[3] = -yp_int_[3] + voq_star - y_int_[9]; - ScalarT ild_star = F_ * y_[13] - wb_ * Cf_ * y_[12] + Kpv_ * (vod_star - y_[11]) + Kiv_ * y_[5]; - ScalarT ilq_star = F_ * y_[14] + wb_ * Cf_ * y_[11] + Kpv_ * (voq_star - y_[12]) + Kiv_ * y_[6]; + ScalarT ild_star = F_ * y_int_[10] - wb_ * Cf_ * y_int_[9] + Kpv_ * (vod_star - y_int_[8]) + Kiv_ * y_int_[2]; + ScalarT ilq_star = F_ * y_int_[11] + wb_ * Cf_ * y_int_[8] + Kpv_ * (voq_star - y_int_[9]) + Kiv_ * y_int_[3]; // Current control - f_[7] = -yp_[7] + ild_star - y_[9]; - f_[8] = -yp_[8] + ilq_star - y_[10]; + f_int_[4] = -yp_int_[4] + ild_star - y_int_[6]; + f_int_[5] = -yp_int_[5] + ilq_star - y_int_[7]; - ScalarT vid_star = -wb_ * Lf_ * y_[10] + Kpc_ * (ild_star - y_[9]) + Kic_ * y_[7]; - ScalarT viq_star = wb_ * Lf_ * y_[9] + Kpc_ * (ilq_star - y_[10]) + Kic_ * y_[8]; + ScalarT vid_star = -wb_ * Lf_ * y_int_[7] + Kpc_ * (ild_star - y_int_[6]) + Kic_ * y_int_[4]; + ScalarT viq_star = wb_ * Lf_ * y_int_[6] + Kpc_ * (ilq_star - y_int_[7]) + Kic_ * y_int_[5]; // Output LC Filter - f_[9] = -yp_[9] - (rLf_ / Lf_) * y_[9] + omega * y_[10] + (vid_star - y_[11]) / Lf_; - f_[10] = -yp_[10] - (rLf_ / Lf_) * y_[10] - omega * y_[9] + (viq_star - y_[12]) / Lf_; + f_int_[6] = -yp_int_[6] - (rLf_ / Lf_) * y_int_[6] + omega * y_int_[7] + (vid_star - y_int_[8]) / Lf_; + f_int_[7] = -yp_int_[7] - (rLf_ / Lf_) * y_int_[7] - omega * y_int_[6] + (viq_star - y_int_[9]) / Lf_; - f_[11] = -yp_[11] + omega * y_[12] + (y_[9] - y_[13]) / Cf_; - f_[12] = -yp_[12] - omega * y_[11] + (y_[10] - y_[14]) / Cf_; + f_int_[8] = -yp_int_[8] + omega * y_int_[9] + (y_int_[6] - y_int_[10]) / Cf_; + f_int_[9] = -yp_int_[9] - omega * y_int_[8] + (y_int_[7] - y_int_[11]) / Cf_; // Output Connector - f_[13] = -yp_[13] - (rLc_ / Lc_) * y_[13] + omega * y_[14] + (y_[11] - vbd_in) / Lc_; - f_[14] = -yp_[14] - (rLc_ / Lc_) * y_[14] - omega * y_[13] + (y_[12] - vbq_in) / Lc_; + f_int_[10] = -yp_int_[10] - (rLc_ / Lc_) * y_int_[10] + omega * y_int_[11] + (y_int_[8] - vbd_in) / Lc_; + f_int_[11] = -yp_int_[11] - (rLc_ / Lc_) * y_int_[11] - omega * y_int_[10] + (y_int_[9] - vbq_in) / Lc_; // Rotor difference angle if (!refframe_) - f_[15] = -yp_[15] + omega - y_[0]; + f_int_[12] = -yp_int_[12] + omega - y_[0]; return 0; } @@ -140,7 +134,8 @@ namespace GridKit template int DistributedGenerator::evaluateExternalResidual() { - ScalarT omega = wb_ - mp_ * y_[3]; + ScalarT omega = wb_ - mp_ * y_int_[0]; + ScalarT delta = refframe_ ? ScalarT(0.0) : y_int_[12]; // ref common ref motor angle if (refframe_) { @@ -153,8 +148,8 @@ namespace GridKit // output // current transformed to common frame - f_[1] = std::cos(y_[15]) * y_[13] - std::sin(y_[15]) * y_[14]; - f_[2] = std::sin(y_[15]) * y_[13] + std::cos(y_[15]) * y_[14]; + f_[1] = std::cos(delta) * y_int_[10] - std::sin(delta) * y_int_[11]; + f_[2] = std::sin(delta) * y_int_[10] + std::cos(delta) * y_int_[11]; return 0; } @@ -194,31 +189,39 @@ namespace GridKit std::vector rtemp{}; std::vector valtemp{}; + RealT delta = refframe_ ? RealT(0.0) : static_cast(y_int_[12]); + // Create dF/dy // r = 1 - ctemp = {13, 14, 15}; + ctemp = {13, 14}; + if (!refframe_) + ctemp.push_back(15); rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(1); valtemp = { - cos(static_cast(y_[15])), - -sin(static_cast(y_[15])), - -sin(static_cast(y_[15])) * static_cast(y_[13]) - cos(static_cast(y_[15])) * static_cast(y_[14]), + cos(delta), + -sin(delta), }; + if (!refframe_) + valtemp.push_back(-sin(delta) * static_cast(y_int_[10]) - cos(delta) * static_cast(y_int_[11])); this->setJacValues(rtemp, ctemp, valtemp); // r = 2 - ctemp = {13, 14, 15}; + ctemp = {13, 14}; + if (!refframe_) + ctemp.push_back(15); rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(2); valtemp = { - sin(static_cast(y_[15])), - cos(static_cast(y_[15])), - cos(static_cast(y_[15])) * static_cast(y_[13]) - sin(static_cast(y_[15])) * static_cast(y_[14]), + sin(delta), + cos(delta), }; + if (!refframe_) + valtemp.push_back(cos(delta) * static_cast(y_int_[10]) - sin(delta) * static_cast(y_int_[11])); this->setJacValues(rtemp, ctemp, valtemp); // r = 0 @@ -238,10 +241,10 @@ namespace GridKit for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(3); valtemp = {-wc_ - alpha_, - wc_ * static_cast(y_[13]), - wc_ * static_cast(y_[14]), - wc_ * static_cast(y_[11]), - wc_ * static_cast(y_[12])}; + wc_ * static_cast(y_int_[10]), + wc_ * static_cast(y_int_[11]), + wc_ * static_cast(y_int_[8]), + wc_ * static_cast(y_int_[9])}; this->setJacValues(rtemp, ctemp, valtemp); // r = 4 @@ -250,10 +253,10 @@ namespace GridKit for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(4); valtemp = {-wc_ - alpha_, - -wc_ * static_cast(y_[14]), - wc_ * static_cast(y_[13]), - wc_ * static_cast(y_[12]), - -wc_ * static_cast(y_[11])}; + -wc_ * static_cast(y_int_[11]), + wc_ * static_cast(y_int_[10]), + wc_ * static_cast(y_int_[9]), + -wc_ * static_cast(y_int_[8])}; this->setJacValues(rtemp, ctemp, valtemp); // r = 5 @@ -293,12 +296,12 @@ namespace GridKit rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(9); - valtemp = {-mp_ * static_cast(y_[10]), + valtemp = {-mp_ * static_cast(y_int_[7]), -(Kpc_ * Kpv_ * nq_) / Lf_, (Kpc_ * Kiv_) / Lf_, Kic_ / Lf_, -(Kpc_ + rLf_) / Lf_ - alpha_, - -mp_ * static_cast(y_[3]), + -mp_ * static_cast(y_int_[0]), -(Kpc_ * Kpv_ + 1.0) / Lf_, -(Cf_ * Kpc_ * wb_) / Lf_, (F_ * Kpc_) / Lf_}; @@ -309,10 +312,10 @@ namespace GridKit rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(10); - valtemp = {mp_ * static_cast(y_[9]), + valtemp = {mp_ * static_cast(y_int_[6]), (Kiv_ * Kpc_) / Lf_, Kic_ / Lf_, - mp_ * static_cast(y_[3]), + mp_ * static_cast(y_int_[0]), -(Kpc_ + rLf_) / Lf_ - alpha_, (Cf_ * Kpc_ * wb_) / Lf_, -(Kpc_ * Kpv_ + 1.0) / Lf_, @@ -324,10 +327,10 @@ namespace GridKit rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(11); - valtemp = {-mp_ * static_cast(y_[12]), + valtemp = {-mp_ * static_cast(y_int_[9]), 1.0 / Cf_, -alpha_, - wb_ - mp_ * static_cast(y_[3]), + wb_ - mp_ * static_cast(y_int_[0]), -1.0 / Cf_}; this->setJacValues(rtemp, ctemp, valtemp); @@ -336,39 +339,45 @@ namespace GridKit rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(12); - valtemp = {mp_ * static_cast(y_[11]), 1.0 / Cf_, -wb_ + mp_ * static_cast(y_[3]), -alpha_, -1.0 / Cf_}; + valtemp = {mp_ * static_cast(y_int_[8]), 1.0 / Cf_, -wb_ + mp_ * static_cast(y_int_[0]), -alpha_, -1.0 / Cf_}; this->setJacValues(rtemp, ctemp, valtemp); // r = 13 - ctemp = {1, 2, 3, 11, 13, 14, 15}; + ctemp = {1, 2, 3, 11, 13, 14}; + if (!refframe_) + ctemp.push_back(15); rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(13); valtemp = { - (1.0 / Lc_) * -cos(static_cast(y_[15])), - (1.0 / Lc_) * -sin(static_cast(y_[15])), - -mp_ * static_cast(y_[14]), + (1.0 / Lc_) * -cos(delta), + (1.0 / Lc_) * -sin(delta), + -mp_ * static_cast(y_int_[11]), 1.0 / Lc_, -rLc_ / Lc_ - alpha_, - wb_ - mp_ * static_cast(y_[3]), - (1.0 / Lc_) * (sin(static_cast(y_[15])) * static_cast(y_[1]) - cos(static_cast(y_[15])) * static_cast(y_[2])), + wb_ - mp_ * static_cast(y_int_[0]), }; + if (!refframe_) + valtemp.push_back((1.0 / Lc_) * (sin(delta) * static_cast(y_[1]) - cos(delta) * static_cast(y_[2]))); this->setJacValues(rtemp, ctemp, valtemp); // r = 14 - ctemp = {1, 2, 3, 12, 13, 14, 15}; + ctemp = {1, 2, 3, 12, 13, 14}; + if (!refframe_) + ctemp.push_back(15); rtemp.clear(); for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(14); valtemp = { - (1.0 / Lc_) * sin(static_cast(y_[15])), - (1.0 / Lc_) * -cos(static_cast(y_[15])), - mp_ * static_cast(y_[13]), + (1.0 / Lc_) * sin(delta), + (1.0 / Lc_) * -cos(delta), + mp_ * static_cast(y_int_[10]), 1.0 / Lc_, - -wb_ + mp_ * static_cast(y_[3]), + -wb_ + mp_ * static_cast(y_int_[0]), -rLc_ / Lc_ - alpha_, - (1.0 / Lc_) * (cos(static_cast(y_[15])) * static_cast(y_[1]) + sin(static_cast(y_[15])) * static_cast(y_[2])), }; + if (!refframe_) + valtemp.push_back((1.0 / Lc_) * (cos(delta) * static_cast(y_[1]) + sin(delta) * static_cast(y_[2]))); this->setJacValues(rtemp, ctemp, valtemp); if (!refframe_) @@ -394,11 +403,6 @@ namespace GridKit this->setExternalConnectionNodes(1, node_bus_->getNodeConnection(0)); this->setExternalConnectionNodes(2, node_bus_->getNodeConnection(1)); - if (refframe_) - { - this->setExternalConnectionNodes(15, INVALID_INDEX); - } - return 0; } diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp index 3b6aaba3f..3632bde33 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -48,9 +48,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp index 20f0cc5d3..337c4fe7c 100644 --- a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp +++ b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp @@ -70,21 +70,21 @@ namespace GridKit template int InductionMotor::evaluateInternalResidual() { - 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 / std::sqrt(3.0)) * (-y_[1] + y_[2]) - Rs_ * y_[6] - (Lls_ + Lms_) * yp_[6] - Lms_ * yp_[5]; - f_[7] = (y_[0] + y_[1] + y_[2]) / 3.0 - Rs_ * y_[7] - Lls_ * yp_[7]; - f_[8] = Rr_ * y_[8] + (Llr_ + Lms_) * yp_[8] + Lms_ * yp_[5] - (P_ / 2.0) * y_[3] * ((Llr_ + Lms_) * y_[9] + Lms_ * y_[6]); - f_[9] = Rr_ * y_[9] + (Llr_ + Lms_) * yp_[9] + Lms_ * yp_[6] + (P_ / 2.0) * y_[3] * ((Llr_ + Lms_) * y_[8] + Lms_ * y_[5]); + f_int_[0] = (1.0 / 3.0) * (2.0 * y_[0] - y_[1] - y_[2]) - Rs_ * y_int_[0] - (Lls_ + Lms_) * yp_int_[0] - Lms_ * yp_int_[1]; + f_int_[1] = (1.0 / std::sqrt(3.0)) * (-y_[1] + y_[2]) - Rs_ * y_int_[1] - (Lls_ + Lms_) * yp_int_[1] - Lms_ * yp_int_[0]; + f_int_[2] = (y_[0] + y_[1] + y_[2]) / 3.0 - Rs_ * y_int_[2] - Lls_ * yp_int_[7]; + f_int_[3] = Rr_ * y_int_[3] + (Llr_ + Lms_) * yp_int_[8] + Lms_ * yp_int_[0] - (P_ / 2.0) * y_[3] * ((Llr_ + Lms_) * y_int_[4] + Lms_ * y_int_[1]); + f_int_[4] = Rr_ * y_int_[4] + (Llr_ + Lms_) * yp_int_[9] + Lms_ * yp_int_[1] + (P_ / 2.0) * y_[3] * ((Llr_ + Lms_) * y_int_[3] + Lms_ * y_int_[0]); return 0; } template int InductionMotor::evaluateExternalResidual() { - f_[0] = y_[5] + y_[7]; - f_[1] = (-1.0 / 2.0) * y_[5] - (std::sqrt(3.0) / 2.0) * y_[6] + y_[7]; - f_[2] = (-1.0 / 2.0) * y_[5] + (std::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_[0] = y_int_[0] + y_int_[2]; + f_[1] = (-1.0 / 2.0) * y_int_[0] - (std::sqrt(3.0) / 2.0) * y_int_[1] + y_int_[2]; + f_[2] = (-1.0 / 2.0) * y_int_[0] + (std::sqrt(3.0) / 2.0) * y_int_[1] + y_int_[2]; + f_[3] = RJ_ * yp_[3] - (3.0 / 4.0) * P_ * Lms_ * (y_[5] * y_int_[4] - y_int_[1] * y_int_[3]); f_[4] = yp_[4] - y_[3]; return 0; } diff --git a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp index e4f0c2f6a..89ea72419 100644 --- a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp +++ b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp @@ -26,9 +26,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp index e7fcdfb22..00248ad29 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp @@ -58,7 +58,7 @@ namespace GridKit template int Inductor::evaluateInternalResidual() { - f_[2] = -L_ * yp_[2] + y_[1] - y_[0]; + f_int_[0] = -L_ * yp_int_[0] + y_[1] - y_[0]; return 0; } @@ -66,9 +66,9 @@ namespace GridKit int Inductor::evaluateExternalResidual() { // input - f_[0] = -y_[2]; + f_[0] = -y_int_[0]; // output - f_[1] = y_[2]; + f_[1] = y_int_[0]; return 0; } diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp index 8ccd72612..ba02ff7a9 100644 --- a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp +++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp @@ -28,9 +28,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp index 3a10fc7da..615b43578 100644 --- a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp +++ b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp @@ -71,16 +71,16 @@ namespace GridKit template int LinearTransformer::evaluateInternalResidual() { - f_[2] = y_[0] - R0_ * y_[2] - L0_ * yp_[2] - M_ * yp_[3]; - f_[3] = y_[1] - R1_ * y_[3] - M_ * yp_[2] - L1_ * yp_[3]; + f_int_[0] = y_[0] - R0_ * y_int_[0] - L0_ * yp_int_[0] - M_ * yp_int_[1]; + f_int_[1] = y_[1] - R1_ * y_int_[1] - M_ * yp_int_[0] - L1_ * yp_int_[1]; return 0; } template int LinearTransformer::evaluateExternalResidual() { - f_[0] = y_[2]; - f_[1] = y_[3]; + f_[0] = y_int_[0]; + f_[1] = y_int_[1]; return 0; } diff --git a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp index a83bb8879..6e30b185c 100644 --- a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp +++ b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp @@ -26,9 +26,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index 2c7ff4d18..e44c131c4 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -28,9 +28,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index 712c3af17..b5f6652dc 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -72,8 +72,8 @@ namespace GridKit template int MicrogridLine::evaluateInternalResidual() { - f_[5] = -yp_[5] - (R_ / L_) * y_[5] + y_[0] * y_[6] + (y_[1] - y_[3]) / L_; - f_[6] = -yp_[6] - (R_ / L_) * y_[6] - y_[0] * y_[5] + (y_[2] - y_[4]) / L_; + f_int_[0] = -yp_int_[0] - (R_ / L_) * y_int_[0] + y_[0] * y_int_[1] + (y_[1] - y_[3]) / L_; + f_int_[1] = -yp_int_[1] - (R_ / L_) * y_int_[1] - y_[0] * y_int_[0] + (y_[2] - y_[4]) / L_; return 0; } @@ -85,12 +85,12 @@ namespace GridKit f_[0] = 0.0; // Port 1 - f_[1] = -y_[5]; - f_[2] = -y_[6]; + f_[1] = -y_int_[0]; + f_[2] = -y_int_[1]; // Port 2 - f_[3] = y_[5]; - f_[4] = y_[6]; + f_[3] = y_int_[0]; + f_[4] = y_int_[1]; return 0; } @@ -117,12 +117,12 @@ namespace GridKit std::vector rcord(ccord.size(), 5); std::vector vals{}; - vals = {static_cast(y_[6]), (1.0 / L_), -(1.0 / L_), -(R_ / L_) - alpha_, static_cast(y_[0])}; + vals = {static_cast(y_int_[1]), (1.0 / L_), -(1.0 / L_), -(R_ / L_) - alpha_, static_cast(y_[0])}; this->setJacValues(rcord, ccord, vals); std::vector ccor2{0, 2, 4, 5, 6}; std::fill(rcord.begin(), rcord.end(), 6); - vals = {-static_cast(y_[5]), (1.0 / L_), -(1.0 / L_), -static_cast(y_[0]), -(R_ / L_) - alpha_}; + vals = {-static_cast(y_int_[0]), (1.0 / L_), -(1.0 / L_), -static_cast(y_[0]), -(R_ / L_) - alpha_}; this->setJacValues(rcord, ccor2, vals); return 0; diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp index 79fa5334c..c36ce9632 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -28,9 +28,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index 7ba54d06a..caf203c44 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -67,8 +67,8 @@ namespace GridKit template int MicrogridLoad::evaluateInternalResidual() { - f_[3] = -yp_[3] - (R_ / L_) * y_[3] + y_[0] * y_[4] + y_[1] / L_; - f_[4] = -yp_[4] - (R_ / L_) * y_[4] - y_[0] * y_[3] + y_[2] / L_; + f_int_[0] = -yp_int_[0] - (R_ / L_) * y_int_[0] + y_[0] * y_int_[1] + y_[1] / L_; + f_int_[1] = -yp_int_[1] - (R_ / L_) * y_int_[1] - y_[0] * y_int_[0] + y_[2] / L_; return 0; } @@ -82,8 +82,8 @@ namespace GridKit // only input for loads // input - f_[1] = -y_[3]; - f_[2] = -y_[4]; + f_[1] = -y_int_[0]; + f_[2] = -y_int_[1]; return 0; } @@ -110,12 +110,12 @@ namespace GridKit std::vector rcord(ccord.size(), 3); std::vector vals{}; - vals = {static_cast(y_[4]), (1.0 / L_), -(R_ / L_) - alpha_, static_cast(y_[0])}; + vals = {static_cast(y_int_[1]), (1.0 / L_), -(R_ / L_) - alpha_, static_cast(y_[0])}; this->setJacValues(rcord, ccord, vals); std::vector ccor2{0, 2, 3, 4}; std::fill(rcord.begin(), rcord.end(), 4); - vals = {-static_cast(y_[3]), (1.0 / L_), -static_cast(y_[0]), -(R_ / L_) - alpha_}; + vals = {-static_cast(y_int_[0]), (1.0 / L_), -static_cast(y_[0]), -(R_ / L_) - alpha_}; this->setJacValues(rcord, ccor2, vals); return 0; diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 63a5df20c..99687d3da 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -28,9 +28,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/NodeBase.hpp b/GridKit/Model/PowerElectronics/NodeBase.hpp index 51f40f975..9e9748529 100644 --- a/GridKit/Model/PowerElectronics/NodeBase.hpp +++ b/GridKit/Model/PowerElectronics/NodeBase.hpp @@ -103,6 +103,7 @@ namespace GridKit int setBusID(IdxT bus_id) { bus_id_ = bus_id; + return 0; } IdxT busID() const diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp index 542e3d786..b5ad84e80 100644 --- a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp +++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp @@ -28,9 +28,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp index 3f22e37bf..6231421e5 100644 --- a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp +++ b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp @@ -94,18 +94,18 @@ namespace GridKit ScalarT llkq1 = static_cast(std::get<0>(Llkq_)); [[maybe_unused]] ScalarT llkq2 = static_cast(std::get<1>(Llkq_)); - ScalarT cos1 = std::cos((P_ / 2.0) * y_[5]); - ScalarT sin1 = std::sin((P_ / 2.0) * y_[5]); - ScalarT cos23m = std::cos((P_ / 2.0) * y_[5] - (2.0 / 3.0) * M_PI); - ScalarT sin23m = std::sin((P_ / 2.0) * y_[5] - (2.0 / 3.0) * M_PI); - ScalarT cos23p = std::cos((P_ / 2.0) * y_[5] + (2.0 / 3.0) * M_PI); - ScalarT sin23p = std::sin((P_ / 2.0) * y_[5] + (2.0 / 3.0) * M_PI); - - 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]); - f_[7] = (-1.0 / 3.0) * (y_[0] + y_[1] + y_[2]) + Rs_ * y_[8] + Lls_ * yp_[8]; - f_[8] = rkq1 * y_[9] + (llkq1 + Lmq_) * yp_[9] + Lmq_ * yp_[6] + Lmq_ * yp_[10]; - f_[9] = rkq1 * y_[9] + (llkq1 + Lmq_) * yp_[9] + Lmq_ * yp_[6] + Lmq_ * yp_[10]; + ScalarT cos1 = std::cos((P_ / 2.0) * y_int_[0]); + ScalarT sin1 = std::sin((P_ / 2.0) * y_int_[0]); + ScalarT cos23m = std::cos((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * M_PI); + ScalarT sin23m = std::sin((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * M_PI); + ScalarT cos23p = std::cos((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * M_PI); + ScalarT sin23p = std::sin((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * M_PI); + + f_int_[0] = (-2.0 / 3.0) * (y_[0] * cos1 + y_[1] * cos23m + y_[2] * cos23p) + Rs_ * y_int_[1] + (Lls_ + Lmq_) * yp_int_[1] + Lmq_ * yp_int_[4] + Lmq_ * yp_int_[5] + y_[4] * (P_ / 2.0) * ((Lls_ + Lmd_) * y_int_[2] + Lmd_ * y_int_[6] + Lmd_ * y_int_[7]); + f_int_[1] = (-2.0 / 3.0) * (y_[0] * sin1 - y_[1] * sin23m - y_[2] * sin23p) + Rs_ * y_int_[2] + (Lls_ + Lmd_) * yp_int_[2] + Lmd_ * yp_int_[6] + Lmd_ * yp_int_[7] - y_[4] * (P_ / 2.0) * ((Lls_ + Lmq_) * y_int_[1] + Lmq_ * y_int_[4] + Lmq_ * y_int_[5]); + f_int_[2] = (-1.0 / 3.0) * (y_[0] + y_[1] + y_[2]) + Rs_ * y_int_[3] + Lls_ * yp_int_[3]; + f_int_[3] = rkq1 * y_int_[4] + (llkq1 + Lmq_) * yp_int_[4] + Lmq_ * yp_int_[1] + Lmq_ * yp_int_[5]; + f_int_[4] = rkq1 * y_int_[4] + (llkq1 + Lmq_) * yp_int_[4] + Lmq_ * yp_int_[1] + Lmq_ * yp_int_[5]; return 0; } @@ -115,18 +115,18 @@ namespace GridKit [[maybe_unused]] ScalarT rkq2 = static_cast(std::get<1>(Rkq_)); [[maybe_unused]] ScalarT llkq2 = static_cast(std::get<1>(Llkq_)); - ScalarT cos1 = std::cos((P_ / 2.0) * y_[5]); - ScalarT sin1 = std::sin((P_ / 2.0) * y_[5]); - ScalarT cos23m = std::cos((P_ / 2.0) * y_[5] - (2.0 / 3.0) * M_PI); - ScalarT sin23m = std::sin((P_ / 2.0) * y_[5] - (2.0 / 3.0) * M_PI); - ScalarT cos23p = std::cos((P_ / 2.0) * y_[5] + (2.0 / 3.0) * M_PI); - ScalarT sin23p = std::sin((P_ / 2.0) * y_[5] + (2.0 / 3.0) * M_PI); - - 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_[4] = yp_[5] - y_[4]; + ScalarT cos1 = std::cos((P_ / 2.0) * y_int_[0]); + ScalarT sin1 = std::sin((P_ / 2.0) * y_int_[0]); + ScalarT cos23m = std::cos((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * M_PI); + ScalarT sin23m = std::sin((P_ / 2.0) * y_int_[0] - (2.0 / 3.0) * M_PI); + ScalarT cos23p = std::cos((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * M_PI); + ScalarT sin23p = std::sin((P_ / 2.0) * y_int_[0] + (2.0 / 3.0) * M_PI); + + f_[0] = y_int_[1] * cos1 + y_int_[2] * sin1 + y_int_[3]; + f_[1] = y_int_[1] * cos23m + y_int_[2] * sin23m + y_int_[3]; + f_[2] = y_int_[1] * cos23p + y_int_[2] * sin23p + y_int_[3]; + f_[3] = RJ_ * yp_[4] - (3.0 / 4.0) * P_ * (Lmd_ * y_int_[1] * (y_int_[2] + y_int_[6] + y_int_[7]) - Lmq_ * y_int_[2] * (y_int_[1] + y_int_[4] + y_[0])); + f_[4] = yp_int_[0] - y_[4]; return 0; } diff --git a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp index c24eed38a..027efcf21 100644 --- a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp +++ b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp @@ -28,9 +28,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index a10607b34..751ab8441 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -73,8 +73,11 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::rel_tol_; using CircuitComponent::abs_tol_; @@ -165,6 +168,9 @@ namespace GridKit * @post System model vectors allocated with size s * @post CSR Jacobian sparsity pattern is computed * @post COO->CSR mapping is computed + * @post Every component's \ref CircuitComponent::y_int_, \ref CircuitComponent::yp_int_, and \ref CircuitComponent::f_int_ pointers + * are set to their appropriate offsets in the system vector, allowing them to directly access their internal variables, derivatives, + * and residuals. * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ @@ -206,11 +212,18 @@ namespace GridKit } { + // The offset for each component's internal variables in the system vector. + // They start at 0, and are stacked of each other. size_t component_internal_idx = 0; for (component_type* comp : components_) { comp->allocate(); + // Update component internal pointers to their correct offsets + comp->setInternalPointer(&y_[component_internal_idx]); + comp->setInternalDerivativePointer(&yp_[component_internal_idx]); + comp->setInternalResidualPointer(&f_[component_internal_idx]); + const auto& external_indices = comp->getExternIndices(); for (size_t i = 0; i < comp->size(); i++) { @@ -330,13 +343,13 @@ namespace GridKit */ int distributeVectors() { - for (const auto& component : components_) + for (component_type* component : components_) { - IdxT size = component->size(); - std::vector& y = component->y(); - std::vector& yp = component->yp(); + std::vector& y = component->y(); + std::vector& yp = component->yp(); + const std::set& externals = component->getExternIndices(); - for (IdxT j = 0; j < size; ++j) + for (size_t j : externals) { if (component->getNodeConnection(j) != neg1_) { @@ -381,14 +394,15 @@ namespace GridKit return err_code; } - for (const auto& component : components_) + for (component_type* component : components_) { if (int err_code = component->evaluateExternalResidual()) return err_code; - IdxT size = component->size(); - const std::vector& residual = component->getResidual(); - for (IdxT j = 0; j < size; ++j) + const std::vector& residual = component->getResidual(); + const std::set& externals = component->getExternIndices(); + + for (size_t j : externals) { //@todo should do a different grounding check if (component->getNodeConnection(j) != neg1_) diff --git a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp index f73c06042..9f518ddb8 100644 --- a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp +++ b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp @@ -92,12 +92,12 @@ namespace GridKit // Internal variables // row 1 - f_[8] = YReMat_ * (V1re - V2re) - (YImMatDi_ * V1im + YImMatOff_ * V2im) - y_[8]; - f_[9] = YReMat_ * (V1im - V2im) + (YImMatDi_ * V1re + YImMatOff_ * V2re) - y_[9]; + f_int_[0] = YReMat_ * (V1re - V2re) - (YImMatDi_ * V1im + YImMatOff_ * V2im) - y_int_[0]; + f_int_[1] = YReMat_ * (V1im - V2im) + (YImMatDi_ * V1re + YImMatOff_ * V2re) - y_int_[1]; // row2 - f_[10] = YReMat_ * (V2re - V1re) - (YImMatOff_ * V1im + YImMatDi_ * V2im) - y_[10]; - f_[11] = YReMat_ * (V2im - V1im) + (YImMatOff_ * V1re + YImMatDi_ * V2re) - y_[11]; + f_int_[2] = YReMat_ * (V2re - V1re) - (YImMatOff_ * V1im + YImMatDi_ * V2im) - y_int_[2]; + f_int_[3] = YReMat_ * (V2im - V1im) + (YImMatOff_ * V1re + YImMatDi_ * V2re) - y_int_[3]; return 0; } @@ -106,17 +106,17 @@ namespace GridKit int TransmissionLine::evaluateExternalResidual() { // input - f_[0] = y_[8]; - f_[1] = y_[9]; + f_[0] = y_int_[0]; + f_[1] = y_int_[1]; - f_[2] = y_[10]; - f_[3] = y_[11]; + f_[2] = y_int_[2]; + f_[3] = y_int_[3]; // ouput - f_[4] = -y_[8]; - f_[5] = -y_[9]; + f_[4] = -y_int_[0]; + f_[5] = -y_int_[1]; - f_[6] = -y_[10]; - f_[7] = -y_[11]; + f_[6] = -y_int_[2]; + f_[7] = -y_int_[3]; return 0; } diff --git a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp index e557b5d7c..d2441a1c7 100644 --- a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp +++ b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp @@ -30,9 +30,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp index 58b3173ad..1f742bade 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp @@ -59,7 +59,7 @@ namespace GridKit int VoltageSource::evaluateInternalResidual() { // internal - f_[2] = y_[1] - y_[0] - V_; + f_int_[0] = y_[1] - y_[0] - V_; return 0; } @@ -67,9 +67,9 @@ namespace GridKit int VoltageSource::evaluateExternalResidual() { // input - f_[0] = -y_[2]; + f_[0] = -y_int_[0]; // ouput - f_[1] = y_[2]; + f_[1] = y_int_[0]; return 0; } diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp index 234513a73..aa79ddd14 100644 --- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp +++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp @@ -28,9 +28,12 @@ namespace GridKit using CircuitComponent::time_; using CircuitComponent::alpha_; using CircuitComponent::y_; + using CircuitComponent::y_int_; using CircuitComponent::yp_; + using CircuitComponent::yp_int_; using CircuitComponent::tag_; using CircuitComponent::f_; + using CircuitComponent::f_int_; using CircuitComponent::g_; using CircuitComponent::yB_; using CircuitComponent::ypB_; diff --git a/examples/Enzyme/PowerElectronics/main.cpp b/examples/Enzyme/PowerElectronics/main.cpp index bb5f2d245..246cb92b6 100644 --- a/examples/Enzyme/PowerElectronics/main.cpp +++ b/examples/Enzyme/PowerElectronics/main.cpp @@ -109,22 +109,18 @@ void evaluateResidual(std::vector y_, std::vector yp_, std::vect // Function that computes the Jacobian via automatic differentiation template -void EnzymeModelJacobian(T* model, DenseMatrix& jac) +void EnzymeModelJacobian(T* model, DenseMatrix& jac, const std::vector& y, const std::vector& res) { size_t N = model->size(); - std::vector y(N); std::vector yp(N); std::vector v(N); - std::vector res(N); std::vector d_res(N); for (size_t idy = 0; idy < N; ++idy) { // Elementary vector for Jacobian-vector product for (size_t idx = 0; idx < N; ++idx) { - y[idx] = (model->y())[idx]; - res[idx] = (model->getResidual())[idx]; - v[idx] = 0.0; + v[idx] = 0.0; } v[idy] = 1.0; @@ -180,19 +176,30 @@ int main() dg->initialize(); dg->updateTime(0.0, 0.0); + std::vector y(dg->size()); + std::vector yp(dg->size()); + std::vector res(dg->size()); + + dg->setInternalPointer(&y[dg->getExternSize()]); + dg->setInternalDerivativePointer(&yp[dg->getExternSize()]); + dg->setInternalResidualPointer(&res[dg->getExternSize()]); + // Residual evaluation and reference Jacobian dg->evaluateResidual(); dg->evaluateJacobian(); - std::vector y = dg->y(); - std::vector yp = dg->yp(); - std::vector res = dg->getResidual(); + + for (size_t i = 0; i < dg->getExternSize(); i++) + { + y[i] = dg->y()[i]; + res[i] = dg->getResidual()[i]; + } DenseMatrix jac_ref_dense(dg->size(), dg->size()); jac_ref_dense.setValues(dg->nnz(), dg->jacobianCooRows(), dg->jacobianCooCols(), dg->jacobianCooValues()); // Enzyme Jacobian DenseMatrix jac_autodiff(dg->size(), dg->size()); - EnzymeModelJacobian(dg, jac_autodiff); + EnzymeModelJacobian(dg, jac_autodiff, y, res); // Check int fail = 0; diff --git a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp index fb927c0be..15709c95f 100644 --- a/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp +++ b/examples/PowerElectronics/DistributedGeneratorTest/DGTest.cpp @@ -46,7 +46,7 @@ int main(int /* argc */, char const** /* argv */) using Bus = GridKit::PowerElectronics::MicrogridBus; Bus bus; - GridKit::DistributedGenerator dg(0, parms, true, &dg_signal, &bus); + GridKit::DistributedGenerator dg(0, parms, false, &dg_signal, &bus); std::vector t1(16, 0.0); std::vector t2{ @@ -67,18 +67,23 @@ int main(int /* argc */, char const** /* argv */) 1.5, 0.3, }; + std::vector res(16, 0.0); dg_signal.allocate(); bus.allocate(); dg.allocate(); - dg.y() = t2; - dg.yp() = t1; + dg.y() = t2; + dg.yp() = t1; + dg.getResidual() = res; + dg.setInternalPointer(&t2[dg.getExternSize()]); + dg.setInternalDerivativePointer(&t1[dg.getExternSize()]); + dg.setInternalResidualPointer(&dg.getResidual()[dg.getExternSize()]); dg.evaluateResidual(); // Generated from matlab code with same parameters and inputs - std::vector true_vec{3.141592277589793e+02, + std::vector true_vec{0, 8.941907747838389e-01, 1.846733023014284e+00, 1.014543000000000e+02, diff --git a/third-party/nlohmann-json b/third-party/nlohmann-json index d70e46bc6..a0e9fb1e6 160000 --- a/third-party/nlohmann-json +++ b/third-party/nlohmann-json @@ -1 +1 @@ -Subproject commit d70e46bc65a21a8eb1ef0e3e9ba707b7b9b78378 +Subproject commit a0e9fb1e638cfbb5b8b556b7c51eaa81977bad48