From 308423841a0d7b21e1ed6fe8caa9d6c104616664 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sat, 8 Aug 2026 17:48:00 -0500 Subject: [PATCH 1/5] Collapse enzyme wb, ws, h --- .../AutomaticDifferentiation/Enzyme/DfDwb.hpp | 176 ------------------ .../AutomaticDifferentiation/Enzyme/DfDws.hpp | 107 ----------- .../AutomaticDifferentiation/Enzyme/DfDy.hpp | 79 +------- .../Enzyme/{DhDwb.hpp => DfDyExt.hpp} | 12 +- .../AutomaticDifferentiation/Enzyme/DfDyp.hpp | 81 +------- .../AutomaticDifferentiation/Enzyme/DhDy.hpp | 103 ---------- .../Enzyme/LowerSparseStorage.hpp | 8 +- .../Enzyme/ModelWrappers.hpp | 159 ++-------------- .../Enzyme/SparseJacobians.hpp | 5 +- .../Model/PhasorDynamics/Branch/Branch.hpp | 31 +-- .../PhasorDynamics/Branch/BranchEnzyme.cpp | 80 ++------ .../PhasorDynamics/Branch/BranchImpl.hpp | 128 ++++++------- .../PhasorDynamics/BusFault/BusFault.hpp | 18 +- .../BusFault/BusFaultEnzyme.cpp | 75 ++++---- .../PhasorDynamics/BusFault/BusFaultImpl.hpp | 113 +++++++---- GridKit/Model/PhasorDynamics/Component.hpp | 38 +++- .../PhasorDynamics/Controller/REECB/Reecb.hpp | 15 +- .../Controller/REECB/ReecbEnzyme.cpp | 45 ++--- .../Controller/REECB/ReecbImpl.hpp | 121 +++++++----- .../PhasorDynamics/Controller/REPCA/Repca.hpp | 12 +- .../Controller/REPCA/RepcaEnzyme.cpp | 42 ++--- .../Controller/REPCA/RepcaImpl.hpp | 167 +++++++++++------ .../PhasorDynamics/Converter/REGCA/Regca.hpp | 19 +- .../Converter/REGCA/RegcaEnzyme.cpp | 130 ++++++------- .../Converter/REGCA/RegcaImpl.hpp | 149 +++++++++------ .../PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp | 22 ++- .../Exciter/ESDC1A/Esdc1aEnzyme.cpp | 100 ++++------ .../Exciter/ESDC1A/Esdc1aImpl.hpp | 134 +++++++------ .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 10 +- .../Exciter/IEEET1/Ieeet1Enzyme.cpp | 101 +++++----- .../Exciter/IEEET1/Ieeet1Impl.hpp | 74 +++++--- .../Exciter/SEXS-PTI/SexsPti.hpp | 11 +- .../Exciter/SEXS-PTI/SexsPtiEnzyme.cpp | 101 +++++----- .../Exciter/SEXS-PTI/SexsPtiImpl.hpp | 73 +++++--- .../Governor/GASTPTI/GastPti.hpp | 10 +- .../Governor/GASTPTI/GastPtiEnzyme.cpp | 84 ++++----- .../Governor/GASTPTI/GastPtiImpl.hpp | 84 +++++---- .../PhasorDynamics/Governor/HYGOV/Hygov.hpp | 10 +- .../Governor/HYGOV/HygovEnzyme.cpp | 82 ++++---- .../Governor/HYGOV/HygovImpl.hpp | 92 +++++---- .../PhasorDynamics/Governor/Tgov1/Tgov1.hpp | 10 +- .../Governor/Tgov1/Tgov1Enzyme.cpp | 85 +++++---- .../Governor/Tgov1/Tgov1Impl.hpp | 64 ++++--- .../Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp | 18 +- .../PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp | 61 +++--- .../PhasorDynamics/Load/LoadZ/LoadZImpl.hpp | 87 ++++++--- .../PhasorDynamics/Load/LoadZIP/LoadZIP.hpp | 18 +- .../Load/LoadZIP/LoadZIPEnzyme.cpp | 62 +++--- .../Load/LoadZIP/LoadZIPImpl.hpp | 91 ++++++--- .../Stabilizer/IEEEST/Ieeest.hpp | 9 +- .../Stabilizer/IEEEST/IeeestEnzyme.cpp | 85 +++++---- .../Stabilizer/IEEEST/IeeestImpl.hpp | 53 ++++-- .../SynchronousMachine/GENROU/Genrou.hpp | 15 +- .../GENROU/GenrouEnzyme.cpp | 128 ++++++------- .../SynchronousMachine/GENROU/GenrouImpl.hpp | 122 +++++++----- .../SynchronousMachine/GENSAL/Gensal.hpp | 15 +- .../GENSAL/GensalEnzyme.cpp | 128 ++++++------- .../SynchronousMachine/GENSAL/GensalImpl.hpp | 122 +++++++----- .../GenClassical/GenClassical.hpp | 15 +- .../GenClassical/GenClassicalEnzyme.cpp | 128 ++++++------- .../GenClassical/GenClassicalImpl.hpp | 123 ++++++++---- GridKit/Model/PhasorDynamics/SystemModel.hpp | 2 + .../Model/PhasorDynamics/SystemModelImpl.hpp | 46 +++-- .../UnitTests/PhasorDynamics/BranchTests.hpp | 4 + .../PhasorDynamics/BusFaultTests.hpp | 47 +++-- .../UnitTests/PhasorDynamics/SystemTests.hpp | 7 + 66 files changed, 2108 insertions(+), 2338 deletions(-) delete mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp delete mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp rename GridKit/AutomaticDifferentiation/Enzyme/{DhDwb.hpp => DfDyExt.hpp} (95%) delete mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp deleted file mode 100644 index fb28f4954..000000000 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp +++ /dev/null @@ -1,176 +0,0 @@ -/** - * @file DfDwb.hpp - * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) - * - */ - -#pragma once - -#include -#include -#include - -namespace GridKit -{ - namespace Enzyme - { - namespace Sparse - { - /** - * @brief Enzyme automatic differentiation Jacobian evaluator : df/dwb - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - */ - template - struct DfDwb - { - using ScalarT = typename ModelT::ScalarT; - using IdxT = typename ModelT::IdxT; - using RealT = typename ModelT::RealT; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Global residual indices - * @param[in] var_indices - Global variable indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] rows - Row indices - * @param[out] cols - Column indices - * @param[out] vals - Values - * @param[out] nnz - Number of nonzeros - */ - static void eval(ModelT* model, - const size_t n_res, - const size_t n_var, - const IdxT* res_indices, - const IdxT* var_indices, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - IdxT* rows, - IdxT* cols, - RealT* vals, - IdxT& nnz) - { - if (n_res > 0 && n_var > 0) - { - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage. @see LowerSparseStorage.hpp - ScalarT* output = __enzyme_todense((void*) ident_load, - (void*) ident_store, - var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, - (void*) sparse_store, - var_i, - 1.0, // value scaling - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic differentiation intrinsic that will be replaced by a derivative - __enzyme_fwddiff((void*) ModelWrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - } - } - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Global residual indices - * @param[in] var_indices - Global variable indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[out] rows - Row indices - * @param[out] cols - Column indices - * @param[out] vals - Values - * @param[out] nnz - Number of nonzeros - */ - static void eval(ModelT* model, - const size_t n_res, - const size_t n_var, - const IdxT* res_indices, - const IdxT* var_indices, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, - IdxT* rows, - IdxT* cols, - RealT* vals, - IdxT& nnz) - { - if (n_res > 0 && n_var > 0) - { - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage. @see LowerSparseStorage.hpp - ScalarT* output = __enzyme_todense((void*) ident_load, - (void*) ident_store, - var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, - (void*) sparse_store, - var_i, - 1.0, // value scaling - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic differentiation intrinsic that will be replaced by a derivative - __enzyme_fwddiff((void*) ModelWrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_const, - ws, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - } - } - }; - } // namespace Sparse - } // namespace Enzyme -} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp deleted file mode 100644 index 5640f7626..000000000 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp +++ /dev/null @@ -1,107 +0,0 @@ -/** - * @file DfDws.hpp - * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) - * - */ - -#pragma once - -#include -#include -#include - -namespace GridKit -{ - namespace Enzyme - { - namespace Sparse - { - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: df/dws - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - */ - template - struct DfDws - { - using ScalarT = typename ModelT::ScalarT; - using IdxT = typename ModelT::IdxT; - using RealT = typename ModelT::RealT; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Global residual indices - * @param[in] var_indices - Global variable indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[out] rows - Row indices - * @param[out] cols - Column indices - * @param[out] vals - Values - * @param[out] nnz - Number of nonzeros - */ - static void eval(ModelT* model, - const size_t n_res, - const size_t n_var, - const IdxT* res_indices, - const IdxT* var_indices, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, - IdxT* rows, - IdxT* cols, - RealT* vals, - IdxT& nnz) - { - if (n_res > 0 && n_var > 0) - { - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage. @see LowerSparseStorage.hpp - ScalarT* output = __enzyme_todense((void*) ident_load, - (void*) ident_store, - var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, - (void*) sparse_store, - var_i, - 1.0, // value scaling - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic differentiation intrinsic that will be replaced by a derivative - __enzyme_fwddiff((void*) ModelWrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_dup, - ws, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - } - } - }; - } // namespace Sparse - } // namespace Enzyme -} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp index 90c26e5e9..8b3fc0d7b 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp @@ -37,7 +37,7 @@ namespace GridKit * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables + * @param[in] y_ext - External variables * @param[out] rows - Row indices * @param[out] cols - Column indices * @param[out] vals - Values @@ -50,7 +50,7 @@ namespace GridKit const IdxT* var_indices, const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, + const ScalarT* y_ext, IdxT* rows, IdxT* cols, RealT* vals, @@ -90,80 +90,7 @@ namespace GridKit enzyme_const, yp, enzyme_const, - wb, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - } - } - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[out] rows - Row indices - * @param[out] cols - Column indices - * @param[out] vals - Values - * @param[out] nnz - Number of nonzeros - */ - static void eval(ModelT* model, - const size_t n_res, - const size_t n_var, - const IdxT* res_indices, - const IdxT* var_indices, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, - IdxT* rows, - IdxT* cols, - RealT* vals, - IdxT& nnz) - { - if (n_res > 0 && n_var > 0) - { - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage. @see LowerSparseStorage.hpp - ScalarT* output = __enzyme_todense((void*) ident_load, - (void*) ident_store, - var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, - (void*) sparse_store, - var_i, - 1.0, // value scaling - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic differentiation intrinsic that will be replaced by a derivative - __enzyme_fwddiff((void*) ModelWrapper::eval, - enzyme_const, - model, - enzyme_dup, - y, - output, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_const, - ws, + y_ext, enzyme_dupnoneed, elementary_v.data(), d_output); diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDyExt.hpp similarity index 95% rename from GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp rename to GridKit/AutomaticDifferentiation/Enzyme/DfDyExt.hpp index b481c80ae..56549b1b0 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDyExt.hpp @@ -1,5 +1,5 @@ /** - * @file DhDwb.hpp + * @file DfDyExt.hpp * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) * */ @@ -17,13 +17,13 @@ namespace GridKit namespace Sparse { /** - * @brief Enzyme automatic differentiation Jacobian evaluator: dh/dwb + * @brief Enzyme automatic differentiation Jacobian evaluator : df/dy_ext * * @tparam ModelT - model type * @tparam MemberFunctions - member function parameter key */ template - struct DhDwb + struct DfDyExt { using ScalarT = typename ModelT::ScalarT; using IdxT = typename ModelT::IdxT; @@ -37,7 +37,7 @@ namespace GridKit * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables + * @param[in] y_ext - External variables * @param[out] rows - Row indices * @param[out] cols - Column indices * @param[out] vals - Values @@ -50,7 +50,7 @@ namespace GridKit const IdxT* var_indices, const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, + const ScalarT* y_ext, IdxT* rows, IdxT* cols, RealT* vals, @@ -89,7 +89,7 @@ namespace GridKit enzyme_const, yp, enzyme_dup, - wb, + y_ext, output, enzyme_dupnoneed, elementary_v.data(), diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDyp.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDyp.hpp index 7da906eef..b2bc21ae5 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDyp.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDyp.hpp @@ -37,7 +37,7 @@ namespace GridKit * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables + * @param[in] y_ext - External variables * @param[in] alpha - Time derivative jacobian coefficient * @param[out] rows - Row indices * @param[out] cols - Column indices @@ -51,7 +51,7 @@ namespace GridKit const IdxT* var_indices, ScalarT* y, ScalarT* yp, - ScalarT* wb, + ScalarT* y_ext, RealT alpha, IdxT* rows, IdxT* cols, @@ -92,82 +92,7 @@ namespace GridKit yp, output, enzyme_const, - wb, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - } - } - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[in] alpha - Time derivative jacobian coefficient - * @param[out] rows - Row indices - * @param[out] cols - Column indices - * @param[out] vals - Values - * @param[out] nnz - Number of nonzeros - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const IdxT* res_indices, - const IdxT* var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - RealT alpha, - IdxT* rows, - IdxT* cols, - RealT* vals, - IdxT& nnz) - { - if (n_res > 0 && n_var > 0) - { - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage. @see LowerSparseStorage.hpp - ScalarT* output = __enzyme_todense((void*) ident_load, - (void*) ident_store, - var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, - (void*) sparse_store, - var_i, - alpha, // value scaling - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic differentiation intrinsic that will be replaced by a derivative - __enzyme_fwddiff((void*) ModelWrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_dup, - yp, - output, - enzyme_const, - wb, - enzyme_const, - ws, + y_ext, enzyme_dupnoneed, elementary_v.data(), d_output); diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp deleted file mode 100644 index 82bdc7873..000000000 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp +++ /dev/null @@ -1,103 +0,0 @@ -/** - * @file DhDy.hpp - * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) - * - */ - -#pragma once - -#include -#include -#include - -namespace GridKit -{ - namespace Enzyme - { - namespace Sparse - { - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: dh/dy - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - */ - template - struct DhDy - { - using ScalarT = typename ModelT::ScalarT; - using IdxT = typename ModelT::IdxT; - using RealT = typename ModelT::RealT; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Global residual indices - * @param[in] var_indices - Global variable indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] rows - Row indices - * @param[out] cols - Column indices - * @param[out] vals - Values - * @param[out] nnz - Number of nonzeros - */ - static void eval(ModelT* model, - const size_t n_res, - const size_t n_var, - const IdxT* res_indices, - const IdxT* var_indices, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - IdxT* rows, - IdxT* cols, - RealT* vals, - IdxT& nnz) - { - if (n_res > 0 && n_var > 0) - { - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage. @see LowerSparseStorage.hpp - ScalarT* output = __enzyme_todense((void*) ident_load, - (void*) ident_store, - var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, - (void*) sparse_store, - var_i, - 1.0, // value scaling - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic differentiation intrinsic that will be replaced by a derivative - __enzyme_fwddiff((void*) ModelWrapper::eval, - enzyme_const, - model, - enzyme_dup, - y, - output, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - } - } - }; - } // namespace Sparse - } // namespace Enzyme -} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp index 6714a28a7..5194a0fab 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp @@ -60,7 +60,7 @@ namespace GridKit { const auto row_mapped = row_indices[static_cast(row)]; const auto col_mapped = col_indices[static_cast(col)]; - if (col_mapped != INVALID_INDEX) + if (row_mapped != INVALID_INDEX && col_mapped != INVALID_INDEX) { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; @@ -99,7 +99,7 @@ namespace GridKit { const auto row_mapped = row_indices[static_cast(row)]; const auto col_mapped = col_indices[static_cast(col)]; - if (col_mapped != INVALID_INDEX) + if (row_mapped != INVALID_INDEX && col_mapped != INVALID_INDEX) { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; @@ -138,7 +138,7 @@ namespace GridKit { const auto row_mapped = row_indices[static_cast(row)]; const auto col_mapped = col_indices[static_cast(col)]; - if (col_mapped != INVALID_INDEX) + if (row_mapped != INVALID_INDEX && col_mapped != INVALID_INDEX) { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; @@ -177,7 +177,7 @@ namespace GridKit { const auto row_mapped = row_indices[static_cast(row)]; const auto col_mapped = col_indices[static_cast(col)]; - if (col_mapped != INVALID_INDEX) + if (row_mapped != INVALID_INDEX && col_mapped != INVALID_INDEX) { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; diff --git a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp index d53a653e6..fe81f6de5 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp @@ -21,12 +21,7 @@ namespace GridKit enum class MemberFunctions { InternalResidual, - InternalResidualWithSignal, - BusResidual, - BusResidual11, //< Special case for branches that are connected to two buses - BusResidual12, //< Special case for branches that are connected to two buses - BusResidual21, //< Special case for branches that are connected to two buses - BusResidual22 //< Special case for branches that are connected to two buses + ExternalResidual }; /** @@ -54,25 +49,25 @@ namespace GridKit * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables + * @param[in] y_ext - External variables * @param[out] f - Internal residual */ static void eval(ModelT* model, const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, + const ScalarT* y_ext, ScalarT* f) { - model->evaluateInternalResidual(y, yp, wb, f); + model->evaluateInternalResidual(y, yp, y_ext, f); } }; /** - * @brief Residual wrapper partial template specialization for InternalResidualWithSignal + * @brief Residual wrapper partial template specialization for ExternalResidual * */ template - struct ModelWrapper + struct ModelWrapper { using ScalarT = typename ModelT::ScalarT; @@ -80,151 +75,19 @@ namespace GridKit * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[out] f - Internal residual - */ - static void eval(ModelT* model, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, - ScalarT* f) - { - model->evaluateInternalResidual(y, yp, wb, ws, f); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual - * - */ - template - struct ModelWrapper - { - using ScalarT = typename ModelT::ScalarT; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) - { - model->evaluateBusResidual(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual11 (branch member function) - * - */ - template - struct ModelWrapper - { - using ScalarT = typename ModelT::ScalarT; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) - { - model->evaluateBusResidual11(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual12 (branch member function) - * - */ - template - struct ModelWrapper - { - using ScalarT = typename ModelT::ScalarT; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) - { - model->evaluateBusResidual12(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual21 (branch member function) - * - */ - template - struct ModelWrapper - { - using ScalarT = typename ModelT::ScalarT; - - /** - * @param[in] model - Pointer to the model to be differentiated - __enzyme_fwddiff((void*) ModelWrapper::eval, - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual + * @param[in] y_ext - External variables + * @param[out] f_ext - External residual */ static void eval(ModelT* model, const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) + const ScalarT* y_ext, + ScalarT* f_ext) { - model->evaluateBusResidual21(y, yp, wb, h); + model->evaluateExternalResidual(y, yp, y_ext, f_ext); } }; - /** - * @brief Residual wrapper partial template specialization for BusResidual22 (branch member function) - * - */ - template - struct ModelWrapper - { - using ScalarT = typename ModelT::ScalarT; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) - { - model->evaluateBusResidual22(y, yp, wb, h); - } - }; } // namespace Sparse } // namespace Enzyme } // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp b/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp index 3ef328159..ddd532b11 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp @@ -6,9 +6,6 @@ #pragma once -#include -#include #include +#include #include -#include -#include diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index b3f3fff8e..bd0c0475c 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -17,6 +17,15 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `Branch` + enum class BranchExternalVariables : size_t + { + VR1, ///< \f$V_{r1}\f$ + VI1, ///< \f$V_{i1}\f$ + VR2, ///< \f$V_{r2}\f$ + VI2, ///< \f$V_{i2}\f$ + }; + /** * @brief Implementation of a line or off-nominal transformer branch between two buses. * @@ -36,8 +45,10 @@ namespace GridKit using Component::yp_; using Component::tag_; using Component::f_; - using Component::wb_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; + using Component::residual_indices_ext_; + using Component::f_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -70,7 +81,9 @@ namespace GridKit virtual int initialize() override final; virtual int tagDifferentiable() override final; virtual int setAbsoluteTolerance(RealT rel_tol) override final; + virtual int evaluateInternalResidual() override final; virtual int evaluateResidual() override final; + virtual int evaluateExternalResidual() override final; virtual int evaluateJacobian() override final; virtual int verify() const override final; @@ -115,6 +128,7 @@ namespace GridKit private: void initializeParameters(const ModelDataT& data); void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); void terminalCurrent1(ScalarT& Ir, ScalarT& Ii); void terminalCurrent2(ScalarT& Ir, ScalarT& Ii); @@ -129,11 +143,6 @@ namespace GridKit ScalarT& Ir, ScalarT& Ii); - static __attribute__((always_inline)) inline void evaluateAdmittanceBlock(const RealT G, - const RealT B, - const ScalarT* wb, - ScalarT* h); - ScalarT& Vr1() { return bus1_->Vr(); @@ -175,13 +184,7 @@ namespace GridKit } public: - __attribute__((always_inline)) inline int evaluateBusResidual11( - const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); - __attribute__((always_inline)) inline int evaluateBusResidual12( - const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); - __attribute__((always_inline)) inline int evaluateBusResidual21( - const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); - __attribute__((always_inline)) inline int evaluateBusResidual22( + __attribute__((always_inline)) inline int evaluateExternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); private: diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp index 9094a60ac..52a940135 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp @@ -23,14 +23,16 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. - auto bus1_size = static_cast(bus1_->size()); - auto bus2_size = static_cast(bus2_->size()); - auto buffer_size = (bus1_size + bus2_size) * (bus1_size + bus2_size); + auto f_ext_size = static_cast(f_ext_.getSize()); + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = f_ext_size * y_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -38,65 +40,19 @@ namespace GridKit nnz_ = 0; - // Bus 1 diagonal Jacobian block owned by the bus - GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11>::eval(this, - static_cast(bus1_->size()), - static_cast((bus1_->y()).getSize()), - (bus1_->getResidualIndices()).data(), - (bus1_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus1_->y().getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - // Bus 2 diagonal Jacobian block owned by the bus - GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22>::eval(this, - static_cast(bus2_->size()), - static_cast((bus2_->y()).getSize()), - (bus2_->getResidualIndices()).data(), - (bus2_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus2_->y().getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - // Off-diagonal Jacobian block (Bus2 variables) owned by the branch - GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12>::eval(this, - static_cast(bus1_->size()), - static_cast((bus2_->y()).getSize()), - (bus1_->getResidualIndices()).data(), - (bus2_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus2_->y().getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - // Off-diagonal Jacobian block (Bus1 variables) owned by the branch - GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21>::eval(this, - static_cast(bus2_->size()), - static_cast((bus1_->y()).getSize()), - (bus2_->getResidualIndices()).data(), - (bus1_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus1_->y().getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_ext_.getSize()), + residual_indices_ext_.data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp index f7b1391f1..68d1000ce 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp @@ -11,6 +11,7 @@ #include +#include #include #include #include @@ -127,8 +128,10 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(4); + f_ext_.setToZero(); + residual_indices_ext_.assign(4, INVALID_INDEX); allocated_ = true; return 0; @@ -197,20 +200,6 @@ namespace GridKit Ii += B * Vr + G * Vi; } - template - __attribute__((always_inline)) inline void Branch::evaluateAdmittanceBlock( - const RealT G, - const RealT B, - const ScalarT* wb, - ScalarT* h) - { - const ScalarT Vr = wb[0]; - const ScalarT Vi = wb[1]; - - h[0] = G * Vr - B * Vi; - h[1] = B * Vr + G * Vi; - } - /** * @brief Compute the absolute tolerance for each variable in the model * @@ -230,88 +219,88 @@ namespace GridKit } /** - * @brief Bus 1 residual contribution from bus 1 variables + * @brief External residual contributions to both terminal buses. * */ template - __attribute__((always_inline)) inline int Branch::evaluateBusResidual11( + __attribute__((always_inline)) inline int Branch::evaluateExternalResidual( [[maybe_unused]] const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) + const ScalarT* y_ext, + ScalarT* f_ext) { - evaluateAdmittanceBlock(g11_, b11_, wb, h); + const ScalarT Vr1 = y_ext[0]; + const ScalarT Vi1 = y_ext[1]; + const ScalarT Vr2 = y_ext[2]; + const ScalarT Vi2 = y_ext[3]; - return 0; - } - - /** - * @brief Bus 1 residual contribution from bus 2 variables - * - */ - template - __attribute__((always_inline)) inline int Branch::evaluateBusResidual12( - [[maybe_unused]] const ScalarT* y, - [[maybe_unused]] const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) - { - evaluateAdmittanceBlock(g12_, b12_, wb, h); + f_ext[0] = (g11_ * Vr1 - b11_ * Vi1) + (g12_ * Vr2 - b12_ * Vi2); + f_ext[1] = (b11_ * Vr1 + g11_ * Vi1) + (b12_ * Vr2 + g12_ * Vi2); + f_ext[2] = (g21_ * Vr1 - b21_ * Vi1) + (g22_ * Vr2 - b22_ * Vi2); + f_ext[3] = (b21_ * Vr1 + g21_ * Vi1) + (b22_ * Vr2 + g22_ * Vi2); return 0; } /** - * @brief Bus 2 residual contribution from bus 1 variables + * @brief The branch owns no internal residual equations. * */ template - __attribute__((always_inline)) int Branch::evaluateBusResidual21( - [[maybe_unused]] const ScalarT* y, - [[maybe_unused]] const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) + int Branch::evaluateInternalResidual() { - evaluateAdmittanceBlock(g21_, b21_, wb, h); - return 0; } /** - * @brief Bus 2 residual contribution from bus 2 variables + * @brief Gather external variables and index maps. * */ template - __attribute__((always_inline)) int Branch::evaluateBusResidual22( - [[maybe_unused]] const ScalarT* y, - [[maybe_unused]] const ScalarT* yp, - const ScalarT* wb, - ScalarT* h) + void Branch::gatherExternalVariables() { - evaluateAdmittanceBlock(g22_, b22_, wb, h); + auto* y_ext = y_ext_.getData(); - return 0; + y_ext[0] = Vr1(); + y_ext[1] = Vi1(); + y_ext[2] = Vr2(); + y_ext[3] = Vi2(); + if (bus1_->size() > 0) + { + variable_indices_ext_[0] = bus1_->getVariableIndex(0); + variable_indices_ext_[1] = bus1_->getVariableIndex(1); + residual_indices_ext_[0] = bus1_->getResidualIndex(0); + residual_indices_ext_[1] = bus1_->getResidualIndex(1); + } + if (bus2_->size() > 0) + { + variable_indices_ext_[2] = bus2_->getVariableIndex(0); + variable_indices_ext_[3] = bus2_->getVariableIndex(1); + residual_indices_ext_[2] = bus2_->getResidualIndex(0); + residual_indices_ext_[3] = bus2_->getResidualIndex(1); + } } /** - * @brief Residual contribution of the branch is computed and pushed to the terminal buses. + * @brief External residual contributions to the two terminal buses. * */ template - int Branch::evaluateResidual() + int Branch::evaluateExternalResidual() { - ScalarT ir1{0.0}; - ScalarT ii1{0.0}; - ScalarT ir2{0.0}; - ScalarT ii2{0.0}; + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); + + gatherExternalVariables(); - terminalCurrent1(ir1, ii1); - terminalCurrent2(ir2, ii2); + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext, f_ext); - Ir1() += ir1; - Ii1() += ii1; - Ir2() += ir2; - Ii2() += ii2; + Ir1() += f_ext[0]; + Ii1() += f_ext[1]; + Ir2() += f_ext[2]; + Ii2() += f_ext[3]; if (bus1_->size() > 0) { @@ -325,6 +314,17 @@ namespace GridKit return 0; } + /** + * @brief Residual contribution of the branch is computed and pushed to the terminal buses. + * + */ + template + int Branch::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } + template void Branch::terminalCurrent1(ScalarT& Ir, ScalarT& Ii) { diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp index 405f78e6c..035ae1ac7 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp @@ -10,6 +10,13 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `BusFault` + enum class BusFaultExternalVariables : size_t + { + VR, ///< \f$V_r\f$ + VI, ///< \f$V_i\f$ + }; + template class BusFault : public Component { @@ -22,8 +29,10 @@ namespace GridKit using Component::yp_; using Component::abs_tol_; using Component::tag_; - using Component::wb_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; + using Component::residual_indices_ext_; + using Component::f_ext_; using Component::f_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; @@ -50,7 +59,9 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; + int evaluateExternalResidual() override final; int evaluateJacobian() override final; int verify() const override final @@ -83,6 +94,7 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; private: + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -106,7 +118,7 @@ namespace GridKit } public: - __attribute__((always_inline)) inline int evaluateBusResidual( + __attribute__((always_inline)) inline int evaluateExternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); __attribute__((always_inline)) inline int evaluateInternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp index 359e79f23..ad0a9d344 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp @@ -23,14 +23,17 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto buffer_size = size * size + 2 * size * bus_size; + auto f_ext_size = static_cast(f_ext_.getSize()); + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = size * size + size * y_ext_size + size * f_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -46,43 +49,51 @@ namespace GridKit (this->getVariableIndices()).data(), y_.getData(), yp_.getData(), - wb_.getData(), + y_ext_.getData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz_); - IdxT nnz_tmp = nnz_; - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus_->y().getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, - static_cast(bus_->size()), - static_cast(y_.getSize()), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + const IdxT internal_external_begin = nnz_; + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + if (!status_) // Value contributions from DfDyExt only when status_ + { + for (IdxT i = internal_external_begin; i < nnz_; ++i) + { + J_vals_buffer_[i] = 0.0; + } + } - if (!status_) // Value contributions from DfDwb and DhDy only when status_ + const IdxT external_residual_begin = nnz_; + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + if (!status_) // External bus-row contributions only when status_ { - for (IdxT i = nnz_tmp; i < nnz_; ++i) + for (IdxT i = external_residual_begin; i < nnz_; ++i) { J_vals_buffer_[i] = 0.0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 39d825a8b..4c0df2de4 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -2,6 +2,7 @@ #include +#include #include #include #include @@ -125,8 +126,10 @@ namespace GridKit residual_indices_.resize(size); // Resize coupling data - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(2); + f_ext_.setToZero(); + residual_indices_ext_.assign(2, INVALID_INDEX); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) @@ -205,20 +208,20 @@ namespace GridKit } /** - * @brief Bus residual + * @brief External residual * */ template - __attribute__((always_inline)) int BusFault::evaluateBusResidual( + __attribute__((always_inline)) int BusFault::evaluateExternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - ScalarT* h) + [[maybe_unused]] const ScalarT* y_ext, + ScalarT* f_ext) { const ScalarT Ir = y[0]; const ScalarT Ii = y[1]; - h[0] = Ir; - h[1] = Ii; + f_ext[0] = Ir; + f_ext[1] = Ii; return 0; } @@ -231,11 +234,11 @@ namespace GridKit __attribute__((always_inline)) int BusFault::evaluateInternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - const ScalarT* wb, + const ScalarT* y_ext, ScalarT* f) { - const ScalarT Vr = wb[0]; - const ScalarT Vi = wb[1]; + const ScalarT Vr = y_ext[0]; + const ScalarT Vi = y_ext[1]; const ScalarT Ir = y[0]; const ScalarT Ii = y[1]; f[0] = Ir + Vr * G_ - Vi * B_; @@ -245,47 +248,89 @@ namespace GridKit } /** - * \brief Residual contribution of the branch is pushed to the - * two terminal buses. + * @brief Gather external variables and index maps. * */ template - int BusFault::evaluateResidual() + void BusFault::gatherExternalVariables() { - auto* wb = wb_.getData(); + auto* y_ext = y_ext_.getData(); + + if (status_) + { + y_ext[0] = Vr(); + y_ext[1] = Vi(); + } + else + { + y_ext[0] = 0.0; + y_ext[1] = 0.0; + } + if (bus_->size() > 0) + { + variable_indices_ext_[0] = bus_->getVariableIndex(0); + variable_indices_ext_[1] = bus_->getVariableIndex(1); + residual_indices_ext_[0] = bus_->getResidualIndex(0); + residual_indices_ext_[1] = bus_->getResidualIndex(1); + } + } + + /** + * \brief Internal residual for the fault model. + * + */ + template + int BusFault::evaluateInternalResidual() + { + gatherExternalVariables(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + auto* f = f_.getData(); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * \brief External residual contributions to the bus while the fault is + * active. + * + */ + template + int BusFault::evaluateExternalResidual() + { + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); if (status_) { - wb[0] = Vr(); - wb[1] = Vi(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); - auto* f = f_.getData(); - auto* h = h_.getData(); - evaluateInternalResidual(y, yp, wb, f); - evaluateBusResidual(y, yp, wb, h); - Ir() += h[0]; - Ii() += h[1]; + evaluateExternalResidual(y, yp, y_ext, f_ext); + Ir() += f_ext[0]; + Ii() += f_ext[1]; if (bus_->size() > 0) { bus_->getResidual().setDataUpdated(); } } - else - { - wb[0] = 0.0; - wb[1] = 0.0; - const auto* y = y_.getData(); - const auto* yp = yp_.getData(); - auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb, f); - } - - f_.setDataUpdated(); return 0; } + /** + * \brief Residual contribution of the fault is pushed to the bus. + * + */ + template + int BusFault::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } + template const Model::VariableMonitorBase* BusFault::getMonitor() const { diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 707902dfc..be975b72a 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include @@ -63,6 +64,17 @@ namespace GridKit virtual int verify() const = 0; + // @todo This is much clearer if we write it this way + virtual int evaluateInternalResidual() + { + return this->evaluateResidual(); + } + + virtual int evaluateExternalResidual() + { + return 0; + } + IdxT size() override final { return size_; @@ -316,6 +328,18 @@ namespace GridKit abs_tol_.resize(n); } + /** + * @brief Allocate this component's external variable vectors. + */ + void allocateExternalVectors(IdxT n) + { + y_ext_.resize(n); + y_ext_.setToZero(); + yp_ext_.resize(n); + yp_ext_.setToZero(); + variable_indices_ext_.assign(static_cast(n), INVALID_INDEX); + } + int constructCoo() { if (coo_jac_ == nullptr) @@ -344,14 +368,19 @@ namespace GridKit IdxT nnz_{0}; /// Global (system-level) variable indices std::vector variable_indices_; + std::vector variable_indices_ext_; /// Global (system-level) residual indices std::vector residual_indices_; + std::vector residual_indices_ext_; VectorT y_; + VectorT y_ext_; VectorT yp_; + VectorT yp_ext_; std::vector tag_; VectorT abs_tol_; VectorT f_; + VectorT f_ext_; bool allocated_{false}; std::vector g_; @@ -378,15 +407,6 @@ namespace GridKit IdxT gridkit_component_id_{0}; - /// Bus interface values seen by the internal residual - VectorT wb_; - /// External residual equations vector - VectorT h_; - /// External signal values seen by the internal residual - VectorT ws_; - /// Global indices of attached external signals - std::vector ws_indices_; - RealT time_; RealT alpha_; diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp index 302361cbe..309972ff8 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp @@ -54,9 +54,11 @@ namespace GridKit IPCMD, ///< \f$I_p^\mathrm{cmd}\f$ Algebraic active-current command output on system base [p.u.] }; - /// External signal variables read or initialized by a `Reecb`. + /// External variables read by a `Reecb`. enum class ReecbExternalVariables : size_t { + VR, ///< \f$V_\mathrm{r}\f$ Terminal-bus real voltage [p.u.] + VI, ///< \f$V_\mathrm{i}\f$ Terminal-bus imaginary voltage [p.u.] PE, ///< \f$P_e\f$ Optional Known active-power feedback input on system base [p.u.] QGEN, ///< \f$Q^\mathrm{gen}\f$ Optional Known reactive-power feedback input on system base [p.u.] QEXT, ///< \f$Q^\mathrm{ext}\f$ Optional Unknown Volt/VAr reference input: system-base reactive power [p.u.], or the terminal-voltage reference [p.u.] when \f$s_Q=1\f$ and \f$s_V=0\f$ @@ -89,10 +91,9 @@ namespace GridKit using Component::va_component_base_; using Component::va_system_base_; using Component::variable_indices_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -121,6 +122,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -134,8 +136,7 @@ namespace GridKit [[gnu::always_inline]] inline int evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f); private: @@ -144,6 +145,8 @@ namespace GridKit static constexpr RealT CURRENT_CIRCLE_KNEE = INITIALIZATION_TOLERANCE / Math::MU; + void gatherExternalVariables(); + struct InitialPoint; struct InitialCurrentLimit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp index 682b18700..f74f0df69 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp @@ -18,12 +18,10 @@ namespace GridKit * @brief Assemble the sparse REECB Jacobian with Enzyme. * * Differentiates the internal residual with respect to state, derivative, - * terminal-bus, and linked signal variables, then constructs the model + * external variables, then constructs the model * COO matrix. * * @pre allocate() has sized the model and Jacobian index maps. - * @pre evaluateResidual() has refreshed the current bus/signal values and - * signal indices. * @pre The containing solver has set the current integration coefficient * and global variable/residual indices. */ @@ -33,12 +31,13 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Reecb...\n"; Log::misc() << "Jacobian evaluation is experimental!\n"; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { const auto size = static_cast(size_); - const auto bus_size = static_cast(bus_->size()); - const auto signal_size = static_cast(ws_.getSize()); - const auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + const auto y_ext_size = static_cast(y_ext_.getSize()); + const auto buffer_size = 2 * size * size + size * y_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; @@ -50,7 +49,7 @@ namespace GridKit nnz_ = 0; - GridKit::Enzyme::Sparse::DfDy::eval( + GridKit::Enzyme::Sparse::DfDy::eval( this, static_cast(f_.getSize()), static_cast(y_.getSize()), @@ -58,14 +57,13 @@ namespace GridKit this->getVariableIndices().data(), y_.getData(), yp_.getData(), - wb_.getData(), - ws_.getData(), + y_ext_.getData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz_); - GridKit::Enzyme::Sparse::DfDyp::eval( + GridKit::Enzyme::Sparse::DfDyp::eval( this, static_cast(f_.getSize()), static_cast(y_.getSize()), @@ -73,39 +71,22 @@ namespace GridKit this->getVariableIndices().data(), y_.getData(), yp_.getData(), - wb_.getData(), - ws_.getData(), + y_ext_.getData(), alpha_, J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz_); - GridKit::Enzyme::Sparse::DfDwb::eval( - this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - this->getResidualIndices().data(), - bus_->getVariableIndices().data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDws::eval( + GridKit::Enzyme::Sparse::DfDyExt::eval( this, static_cast(f_.getSize()), - static_cast(ws_.getSize()), + static_cast(y_ext_.getSize()), this->getResidualIndices().data(), - ws_indices_.data(), + variable_indices_ext_.data(), y_.getData(), yp_.getData(), - wb_.getData(), - ws_.getData(), + y_ext_.getData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp index 0b8fd4615..21604d7b2 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -123,13 +123,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.resize(2); - wb_.setToZero(); - - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); for (IdxT j = 0; j < size_; ++j) { @@ -788,65 +782,104 @@ namespace GridKit } /** - * @brief Evaluate the model residuals + * @brief Gather external variables and index maps. * - * Starts from latched values, refreshes attached signals and their indices, - * refreshes terminal-bus voltage, and evaluates the internal residual. - * REECB contributes no bus residual. + * Unattached signal inputs retain the values latched by initialize(). */ template - int Reecb::evaluateResidual() + void Reecb::gatherExternalVariables() { + auto* y_ext = y_ext_.getData(); + + const auto VR_EXT = static_cast(ReecbExternalVariables::VR); + const auto VI_EXT = static_cast(ReecbExternalVariables::VI); const auto PE = static_cast(ReecbExternalVariables::PE); const auto QGEN = static_cast(ReecbExternalVariables::QGEN); const auto QEXT = static_cast(ReecbExternalVariables::QEXT); const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); const auto PREF = static_cast(ReecbExternalVariables::PREF); - auto* ws = ws_.getData(); + y_ext[VR_EXT] = Vr(); + y_ext[VI_EXT] = Vi(); + variable_indices_ext_[VR_EXT] = INVALID_INDEX; + variable_indices_ext_[VI_EXT] = INVALID_INDEX; + if (bus_->size() > 0) + { + variable_indices_ext_[VR_EXT] = bus_->getVariableIndex(0); + variable_indices_ext_[VI_EXT] = bus_->getVariableIndex(1); + } - ws[PE] = pe_set_; - ws[QGEN] = qgen_set_; - ws[QEXT] = qext_set_; - ws[PFAREF] = pfaref_set_; - ws[PREF] = pref_set_; - std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + y_ext[PE] = pe_set_; + y_ext[QGEN] = qgen_set_; + y_ext[QEXT] = qext_set_; + y_ext[PFAREF] = pfaref_set_; + y_ext[PREF] = pref_set_; + variable_indices_ext_[PE] = INVALID_INDEX; + variable_indices_ext_[QGEN] = INVALID_INDEX; + variable_indices_ext_[QEXT] = INVALID_INDEX; + variable_indices_ext_[PFAREF] = INVALID_INDEX; + variable_indices_ext_[PREF] = INVALID_INDEX; if (auto pe_port = ports_.in.template port()) { - ws[PE] = pe_port.readSignal(); - ws_indices_[PE] = pe_port.signalVariableIndex(); + y_ext[PE] = ports_.in.template port().readSignal(); + variable_indices_ext_[PE] = + ports_.in.template port().signalVariableIndex(); } if (auto qgen_port = ports_.in.template port()) { - ws[QGEN] = qgen_port.readSignal(); - ws_indices_[QGEN] = qgen_port.signalVariableIndex(); + y_ext[QGEN] = ports_.in.template port().readSignal(); + variable_indices_ext_[QGEN] = + ports_.in.template port().signalVariableIndex(); } if (auto qext_port = ports_.in.template port()) { - ws[QEXT] = qext_port.readSignal(); - ws_indices_[QEXT] = qext_port.signalVariableIndex(); + y_ext[QEXT] = ports_.in.template port().readSignal(); + variable_indices_ext_[QEXT] = + ports_.in.template port().signalVariableIndex(); } if (auto pfaref_port = ports_.in.template port()) { - ws[PFAREF] = pfaref_port.readSignal(); - ws_indices_[PFAREF] = pfaref_port.signalVariableIndex(); + y_ext[PFAREF] = + ports_.in.template port().readSignal(); + variable_indices_ext_[PFAREF] = + ports_.in.template port().signalVariableIndex(); } if (auto pref_port = ports_.in.template port()) { - ws[PREF] = pref_port.readSignal(); - ws_indices_[PREF] = pref_port.signalVariableIndex(); + y_ext[PREF] = ports_.in.template port().readSignal(); + variable_indices_ext_[PREF] = + ports_.in.template port().signalVariableIndex(); } + } - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + /** + * @brief Evaluate the internal REECB residual equations. + */ + template + int Reecb::evaluateInternalResidual() + { + gatherExternalVariables(); - evaluateInternalResidual(y_.getData(), yp_.getData(), wb, ws, f_.getData()); + evaluateInternalResidual(y_.getData(), yp_.getData(), y_ext_.getData(), f_.getData()); f_.setDataUpdated(); return 0; } + /** + * @brief Evaluate internal equations and external contributions. + * + * REECB contributes no external residual, so the base implementation + * returns zero after the internal equations are evaluated. + */ + template + int Reecb::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + + /** * @brief Access the optional variable monitor * @@ -866,8 +899,7 @@ namespace GridKit * * @param[in] y Internal variables. * @param[in] yp Internal variable derivatives. - * @param[in] wb Terminal-bus voltage components. - * @param[in] ws External signal values in their documented port units and bases. + * @param[in] y_ext External variable values. * @param[out] f Internal residuals. */ template @@ -875,8 +907,7 @@ namespace GridKit Reecb::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); @@ -902,6 +933,8 @@ namespace GridKit const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + const auto VR_EXT = static_cast(ReecbExternalVariables::VR); + const auto VI_EXT = static_cast(ReecbExternalVariables::VI); const auto PE = static_cast(ReecbExternalVariables::PE); const auto QGEN = static_cast(ReecbExternalVariables::QGEN); const auto QEXT = static_cast(ReecbExternalVariables::QEXT); @@ -938,14 +971,14 @@ namespace GridKit const ScalarT qv_dot = yp[QV]; const ScalarT pord_dot = yp[PORD]; - const ScalarT vr = wb[0]; - const ScalarT vi = wb[1]; + const ScalarT vr = y_ext[VR_EXT]; + const ScalarT vi = y_ext[VI_EXT]; - const ScalarT pe = this->toComponentBase(ws[PE]); - const ScalarT qgen = this->toComponentBase(ws[QGEN]); - const ScalarT extref = ws[QEXT]; - const ScalarT pfaref = ws[PFAREF]; - const ScalarT pref = this->toComponentBase(ws[PREF]); + const ScalarT pe = this->toComponentBase(y_ext[PE]); + const ScalarT qgen = this->toComponentBase(y_ext[QGEN]); + const ScalarT extref = y_ext[QEXT]; + const ScalarT pfaref = y_ext[PFAREF]; + const ScalarT pref = this->toComponentBase(y_ext[PREF]); const ScalarT iqcmd = this->toComponentBase(iqcmd_system); const ScalarT ipcmd = this->toComponentBase(ipcmd_system); diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp index 16f9af8de..fc5d6e3db 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp @@ -57,6 +57,8 @@ namespace GridKit /// External variables of `Repca`. enum class RepcaExternalVariables : size_t { + VR, ///< \f$V_\mathrm{r}\f$ Regulated-bus real voltage component [p.u.] + VI, ///< \f$V_\mathrm{i}\f$ Regulated-bus imaginary voltage component [p.u.] IR, ///< \f$I_\mathrm{r}\f$ Required branch-current real component on system base [p.u.] II, ///< \f$I_\mathrm{i}\f$ Required branch-current imaginary component on system base [p.u.] P, ///< \f$P\f$ Required branch active power on system base [p.u.] @@ -94,10 +96,9 @@ namespace GridKit using Component::va_component_base_; using Component::va_system_base_; using Component::variable_indices_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -125,6 +126,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -138,8 +140,7 @@ namespace GridKit [[gnu::always_inline]] inline int evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f); private: @@ -148,6 +149,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParameters(); bool invertClamp(ScalarT output, RealT lower, RealT upper, ScalarT& input) const; diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp index 1b0242296..b0a2d7c42 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp @@ -22,8 +22,6 @@ namespace GridKit * COO matrix. * * @pre allocate() has sized the model and Jacobian index maps. - * @pre evaluateResidual() has refreshed the current bus/signal values and - * signal indices. * @pre The containing solver has set the current integration coefficient * and global variable/residual indices. */ @@ -33,12 +31,12 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Repca...\n"; Log::misc() << "Jacobian evaluation is experimental!\n"; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { const auto size = static_cast(size_); - const auto bus_size = static_cast(bus_->size()); - const auto signal_size = static_cast(ws_.getSize()); - const auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + const auto buffer_size = 2 * size * size + size * static_cast(y_ext_.getSize()); J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; @@ -50,7 +48,7 @@ namespace GridKit nnz_ = 0; - GridKit::Enzyme::Sparse::DfDy::eval( + GridKit::Enzyme::Sparse::DfDy::eval( this, static_cast(f_.getSize()), static_cast(y_.getSize()), @@ -58,14 +56,13 @@ namespace GridKit this->getVariableIndices().data(), y_.getData(), yp_.getData(), - wb_.getData(), - ws_.getData(), + y_ext_.getData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz_); - GridKit::Enzyme::Sparse::DfDyp::eval( + GridKit::Enzyme::Sparse::DfDyp::eval( this, static_cast(f_.getSize()), static_cast(y_.getSize()), @@ -73,39 +70,22 @@ namespace GridKit this->getVariableIndices().data(), y_.getData(), yp_.getData(), - wb_.getData(), - ws_.getData(), + y_ext_.getData(), alpha_, J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz_); - GridKit::Enzyme::Sparse::DfDwb::eval( - this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - this->getResidualIndices().data(), - bus_->getVariableIndices().data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDws::eval( + GridKit::Enzyme::Sparse::DfDyExt::eval( this, static_cast(f_.getSize()), - static_cast(ws_.getSize()), + static_cast(y_ext_.getSize()), this->getResidualIndices().data(), - ws_indices_.data(), + variable_indices_ext_.data(), y_.getData(), yp_.getData(), - wb_.getData(), - ws_.getData(), + y_ext_.getData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp index dbfd8c6bc..8f607a0e4 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp @@ -83,7 +83,7 @@ namespace GridKit /** * @brief Allocate model vectors and wire assigned command outputs * - * Sizes the state, residual, bus, and signal-interface buffers, initializes + * Sizes the state, residual, and external-variable buffers, initializes * identity index maps, and points assigned `qext` and `pext` nodes at * the internal system-base states that REPCA publishes. Repeated * allocation reuses the existing model vectors and signal links. @@ -104,13 +104,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.resize(2); - wb_.setToZero(); - - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); for (IdxT j = 0; j < size_; ++j) { @@ -586,15 +580,19 @@ namespace GridKit } /** - * @brief Evaluate the REPCA-owned residual rows + * @brief Gather regulated-bus and external-signal values and indices. * - * Refreshes required bus and measurement inputs, starts optional - * references from the values latched by initialize(), then overwrites - * them from attached signals. REPCA contributes no bus residual. + * Required measurements are read from their attached signals. Optional + * inputs start from the values latched by initialize() and are overwritten + * when attached. */ template - int Repca::evaluateResidual() + void Repca::gatherExternalVariables() { + auto* y_ext = y_ext_.getData(); + + const auto VR = static_cast(RepcaExternalVariables::VR); + const auto VI = static_cast(RepcaExternalVariables::VI); const auto IR = static_cast(RepcaExternalVariables::IR); const auto II = static_cast(RepcaExternalVariables::II); const auto P = static_cast(RepcaExternalVariables::P); @@ -605,63 +603,109 @@ namespace GridKit const auto QREF = static_cast(RepcaExternalVariables::QREF); const auto FREQREF = static_cast(RepcaExternalVariables::FREQREF); - auto* ws = ws_.getData(); - - ws[VREF] = vref_set_; - ws[PREF] = pref_set_; - ws[QREF] = qref_set_; - ws[FREQ] = static_cast(ONE); - ws[FREQREF] = freqref_set_; - std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); - - ws[IR] = ports_.in.template port().readSignal(); - ws_indices_[IR] = ports_.in.template port().signalVariableIndex(); - ws[II] = ports_.in.template port().readSignal(); - ws_indices_[II] = ports_.in.template port().signalVariableIndex(); - ws[P] = ports_.in.template port().readSignal(); - ws_indices_[P] = ports_.in.template port().signalVariableIndex(); - ws[Q] = ports_.in.template port().readSignal(); - ws_indices_[Q] = ports_.in.template port().signalVariableIndex(); - if (auto freq_port = ports_.in.template port()) + std::fill(variable_indices_ext_.begin(), + variable_indices_ext_.end(), + INVALID_INDEX); + + y_ext[VR] = Vr(); + y_ext[VI] = Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[VR] = bus_->getVariableIndex(0); + variable_indices_ext_[VI] = bus_->getVariableIndex(1); + } + + y_ext[VREF] = vref_set_; + y_ext[PREF] = pref_set_; + y_ext[QREF] = qref_set_; + y_ext[FREQ] = static_cast(ONE); + y_ext[FREQREF] = freqref_set_; + + y_ext[IR] = + ports_.in.template port().readSignal(); + variable_indices_ext_[IR] = + ports_.in.template port().signalVariableIndex(); + y_ext[II] = + ports_.in.template port().readSignal(); + variable_indices_ext_[II] = + ports_.in.template port().signalVariableIndex(); + y_ext[P] = + ports_.in.template port().readSignal(); + variable_indices_ext_[P] = + ports_.in.template port().signalVariableIndex(); + y_ext[Q] = + ports_.in.template port().readSignal(); + variable_indices_ext_[Q] = + ports_.in.template port().signalVariableIndex(); + if (ports_.in.template port().connected()) { - ws[FREQ] = freq_port.readSignal(); - ws_indices_[FREQ] = freq_port.signalVariableIndex(); + y_ext[FREQ] = + ports_.in.template port().readSignal(); + variable_indices_ext_[FREQ] = + ports_.in.template port().signalVariableIndex(); } if (auto vref_port = ports_.in.template port()) { - ws[VREF] = vref_port.readSignal(); - ws_indices_[VREF] = vref_port.signalVariableIndex(); + y_ext[VREF] = + ports_.in.template port().readSignal(); + variable_indices_ext_[VREF] = + ports_.in.template port().signalVariableIndex(); } if (auto pref_port = ports_.in.template port()) { - ws[PREF] = pref_port.readSignal(); - ws_indices_[PREF] = pref_port.signalVariableIndex(); + y_ext[PREF] = + ports_.in.template port().readSignal(); + variable_indices_ext_[PREF] = + ports_.in.template port().signalVariableIndex(); } if (auto qref_port = ports_.in.template port()) { - ws[QREF] = qref_port.readSignal(); - ws_indices_[QREF] = qref_port.signalVariableIndex(); + y_ext[QREF] = + ports_.in.template port().readSignal(); + variable_indices_ext_[QREF] = + ports_.in.template port().signalVariableIndex(); } if (auto freqref_port = ports_.in.template port()) { - ws[FREQREF] = freqref_port.readSignal(); - ws_indices_[FREQREF] = freqref_port.signalVariableIndex(); + y_ext[FREQREF] = + ports_.in.template port().readSignal(); + variable_indices_ext_[FREQREF] = + ports_.in.template port().signalVariableIndex(); } + } - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + /** + * @brief Evaluate the REPCA-owned residual rows. + */ + template + int Repca::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + /** + * @brief Evaluate internal equations and external contributions. + * + * REPCA contributes no external residual, so the base implementation + * returns zero after the internal equations are evaluated. + */ + template + int Repca::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + + /** * @brief Access the configured monitor * @@ -683,8 +727,8 @@ namespace GridKit * @param[in] y Internal variables in `RepcaInternalVariables` order and * on the bases documented by their enums. * @param[in] yp Internal derivatives in the same enum order and bases. - * @param[in] wb Regulated-bus real and imaginary voltage components. - * @param[in] ws External signals in `RepcaExternalVariables` order. + * @param[in] y_ext Regulated-bus voltage and external signals in + * `RepcaExternalVariables` order. * @param[out] f Caller-provided residual output buffer in * `RepcaInternalVariables` order. */ @@ -693,8 +737,7 @@ namespace GridKit Repca::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { const auto VMEAS = static_cast(RepcaInternalVariables::VMEAS); @@ -720,6 +763,8 @@ namespace GridKit const auto PPI = static_cast(RepcaInternalVariables::PPI); const auto PEXT = static_cast(RepcaInternalVariables::PEXT); + const auto VR = static_cast(RepcaExternalVariables::VR); + const auto VI = static_cast(RepcaExternalVariables::VI); const auto IR = static_cast(RepcaExternalVariables::IR); const auto II = static_cast(RepcaExternalVariables::II); const auto P = static_cast(RepcaExternalVariables::P); @@ -761,18 +806,18 @@ namespace GridKit const ScalarT xppi_dot = yp[XPPI]; const ScalarT pref_dot = yp[PREF_STATE]; - const ScalarT vr = wb[0]; - const ScalarT vi = wb[1]; - - const ScalarT ir = this->toComponentBase(ws[IR]); - const ScalarT ii = this->toComponentBase(ws[II]); - const ScalarT p = this->toComponentBase(ws[P]); - const ScalarT q = this->toComponentBase(ws[Q]); - const ScalarT freq = ws[FREQ]; - const ScalarT freqref = ws[FREQREF]; - const ScalarT vref = ws[VREF]; - const ScalarT qref = this->toComponentBase(ws[QREF]); - const ScalarT pref_in = this->toComponentBase(ws[PREF_INPUT]); + const ScalarT vr = y_ext[VR]; + const ScalarT vi = y_ext[VI]; + + const ScalarT ir = this->toComponentBase(y_ext[IR]); + const ScalarT ii = this->toComponentBase(y_ext[II]); + const ScalarT p = this->toComponentBase(y_ext[P]); + const ScalarT q = this->toComponentBase(y_ext[Q]); + const ScalarT freq = y_ext[FREQ]; + const ScalarT freqref = y_ext[FREQREF]; + const ScalarT vref = y_ext[VREF]; + const ScalarT qref = this->toComponentBase(y_ext[QREF]); + const ScalarT pref_in = this->toComponentBase(y_ext[PREF_INPUT]); const ScalarT vldc_r = vr - Rc_ * ir + Xc_ * ii; const ScalarT vldc_i = vi - Rc_ * ii - Xc_ * ir; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp index bb3ab989c..d9acdaf9e 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -52,6 +52,8 @@ namespace GridKit /// External variables of a `Regca` enum class RegcaExternalVariables : size_t { + VR, ///< \f$V_\mathrm{r}\f$ Terminal-bus real voltage + VI, ///< \f$V_\mathrm{i}\f$ Terminal-bus imaginary voltage IPCMD, ///< \f$I_p^\mathrm{cmd}\f$ Active-current command on system base IQCMD, ///< \f$I_q^\mathrm{cmd}\f$ Reactive-current command on system base }; @@ -70,19 +72,19 @@ namespace GridKit using Component::allocated_; using Component::abs_tol_; using Component::f_; - using Component::h_; + using Component::f_ext_; using Component::J_cols_buffer_; using Component::J_rows_buffer_; using Component::J_vals_buffer_; using Component::nnz_; using Component::residual_indices_; + using Component::residual_indices_ext_; using Component::size_; using Component::tag_; using Component::variable_indices_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -108,6 +110,8 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; + int evaluateExternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -119,15 +123,16 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT* y, const ScalarT* yp, const ScalarT* wb, const ScalarT* ws, ScalarT* f); + const ScalarT* y, const ScalarT* yp, const ScalarT* y_ext, ScalarT* f); - __attribute__((always_inline)) inline int evaluateBusResidual( - const ScalarT* y, const ScalarT* yp, const ScalarT* wb, ScalarT* h); + __attribute__((always_inline)) inline int evaluateExternalResidual( + const ScalarT* y, const ScalarT* yp, const ScalarT* y_ext, ScalarT* f_ext); private: void initializeParameters(const ModelDataT& data); void initializeMonitor(); void setDerivedParameters(); + void gatherExternalVariables(); /** * @brief Smooth approximation of the REGCA `rrpwr` rate limiter. diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp index 935c00320..668fc64b6 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp @@ -16,8 +16,6 @@ namespace GridKit { /** * @brief Assemble the sparse component Jacobian with Enzyme. - * - * @pre evaluateResidual() has run at the current state. */ template int Regca::evaluateJacobian() @@ -25,14 +23,16 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. Enzyme keeps only structural // nonzeros for each differentiated block. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + 2 * size * bus_size + size * signal_size; + auto y_ext_size = static_cast(y_ext_.getSize()); + auto f_ext_size = static_cast(f_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size + size * f_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -44,79 +44,61 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + Fn::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDwb::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + Fn::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DhDy::eval(this, - static_cast(bus_->size()), - static_cast(y_.getSize()), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); return 0; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp index 591fd614c..bc6b97103 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -312,15 +312,11 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.resize(2); - wb_.setToZero(); - h_.resize(2); - h_.setToZero(); - - auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + // Resize coupling data + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(2); + f_ext_.setToZero(); + residual_indices_ext_.assign(2, INVALID_INDEX); for (IdxT j = 0; j < size_; ++j) { @@ -549,8 +545,8 @@ namespace GridKit * * @param[in] y Internal variables. * @param[in] yp Internal variable derivatives. - * @param[in] wb Terminal-bus voltage components. - * @param[in] ws Current-command signal values on the system base. + * @param[in] y_ext External variables: terminal-bus voltage components + * followed by the current-command signals on the system base. * @param[out] f Internal residuals. */ template @@ -558,8 +554,7 @@ namespace GridKit Regca::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { const auto VM = static_cast(RegcaInternalVariables::VM); @@ -573,6 +568,8 @@ namespace GridKit const auto PBR = static_cast(RegcaInternalVariables::PBR); const auto QBR = static_cast(RegcaInternalVariables::QBR); + const auto VR = static_cast(RegcaExternalVariables::VR); + const auto VI = static_cast(RegcaExternalVariables::VI); const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); @@ -591,11 +588,11 @@ namespace GridKit const ScalarT iq_dot = yp[IQ]; const ScalarT ip_dot = yp[IP]; - const ScalarT vr = wb[0]; - const ScalarT vi = wb[1]; + const ScalarT vr = y_ext[VR]; + const ScalarT vi = y_ext[VI]; - const ScalarT ipcmd = this->toComponentBase(ws[IPCMD]); - const ScalarT iqcmd = this->toComponentBase(ws[IQCMD]); + const ScalarT ipcmd = this->toComponentBase(y_ext[IPCMD]); + const ScalarT iqcmd = this->toComponentBase(y_ext[IQCMD]); // Form the unconstrained current derivatives, then apply the REGCA // recovery rate limits in p.u./s. @@ -640,76 +637,124 @@ namespace GridKit * * @param[in] y Internal variables. * @param[in] yp Internal variable derivatives, unused. - * @param[in] wb Terminal-bus voltage components, unused. - * @param[out] h Current injected into the terminal bus. + * @param[in] y_ext External variables, unused. + * @param[out] f_ext Current injected into the terminal bus. */ template - __attribute__((always_inline)) inline int Regca::evaluateBusResidual( + __attribute__((always_inline)) inline int Regca::evaluateExternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - ScalarT* h) + [[maybe_unused]] const ScalarT* y_ext, + ScalarT* f_ext) { const auto IR = static_cast(RegcaInternalVariables::IR); const auto II = static_cast(RegcaInternalVariables::II); - h[0] = y[IR]; - h[1] = y[II]; + f_ext[0] = y[IR]; + f_ext[1] = y[II]; return 0; } /** - * @brief Evaluate model residuals and accumulate the branch current. - * - * Refreshes the bus and signal interface buffers, evaluates the internal - * and bus residuals, and accumulates the converter current into the terminal - * bus. Unattached command ports use the setpoints latched by initialize(). - * - * @pre The terminal-bus residual has been zeroed for this evaluation. + * @brief Gather external variables and index maps. */ template - int Regca::evaluateResidual() + void Regca::gatherExternalVariables() { + auto* y_ext = y_ext_.getData(); + + const auto VR = static_cast(RegcaExternalVariables::VR); + const auto VI = static_cast(RegcaExternalVariables::VI); const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); - auto* ws = ws_.getData(); - - ws[IPCMD] = ipcmd_set_; - ws[IQCMD] = iqcmd_set_; - std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + // Terminal-bus voltage + y_ext[VR] = Vr(); + y_ext[VI] = Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[VR] = bus_->getVariableIndex(0); + variable_indices_ext_[VI] = bus_->getVariableIndex(1); + residual_indices_ext_[0] = bus_->getResidualIndex(0); + residual_indices_ext_[1] = bus_->getResidualIndex(1); + } - if (auto ipcmd_port = ports_.in.template port()) + // Active-current command + y_ext[IPCMD] = ipcmd_set_; + if (ports_.in.template port().connected()) { - ws[IPCMD] = ipcmd_port.readSignal(); - ws_indices_[IPCMD] = ipcmd_port.signalVariableIndex(); + y_ext[IPCMD] = + ports_.in.template port().readSignal(); + variable_indices_ext_[IPCMD] = + ports_.in.template port().signalVariableIndex(); } - if (auto iqcmd_port = ports_.in.template port()) + // Reactive-current command + y_ext[IQCMD] = iqcmd_set_; + if (ports_.in.template port().connected()) { - ws[IQCMD] = iqcmd_port.readSignal(); - ws_indices_[IQCMD] = iqcmd_port.signalVariableIndex(); + y_ext[IQCMD] = + ports_.in.template port().readSignal(); + variable_indices_ext_[IQCMD] = + ports_.in.template port().signalVariableIndex(); } + } - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + /** + * @brief Internal residual for the converter model. + */ + template + int Regca::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - auto* h = h_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); - evaluateBusResidual(y, yp, wb, h); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); - Ir() += h[0]; - Ii() += h[1]; - bus_->getResidual().setDataUpdated(); + return 0; + } + + /** + * @brief External residual contributions to the bus. + * + * @pre The terminal-bus residual has been zeroed for this evaluation. + */ + template + int Regca::evaluateExternalResidual() + { + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + + evaluateExternalResidual(y, yp, y_ext, f_ext); + + // Regca contribution to bus algebraic equations + Ir() += f_ext[0]; + Ii() += f_ext[1]; + + if (bus_->size() > 0) + { + bus_->getResidual().setDataUpdated(); + } return 0; } + + /** + * @brief Evaluate model residuals and accumulate the branch current. + */ + template + int Regca::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } } // namespace Converter } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp index 6d70cb1a9..ed3e00142 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp @@ -43,13 +43,15 @@ namespace GridKit EFD, ///< \f$E_{\mathrm{fd}}\f$ Algebraic field-voltage output [p.u.] }; - /// External signal variables read or initialized by an `Esdc1a`. + /// External variables read by an `Esdc1a`. enum class Esdc1aExternalVariables : size_t { - OMEGA, ///< \f$\omega\f$ Known machine speed deviation [p.u.] - VREF, ///< \f$V_{\mathrm{ref}}\f$ Unknown voltage-control reference [p.u.] - VS, ///< \f$V_S\f$ Known stabilizer input signal [p.u.] - VUEL, ///< \f$V_{\mathrm{UEL}}\f$ Known under-excitation limiter input [p.u.] + VR, ///< \f$V_\mathrm{r}\f$ Terminal-bus real voltage [p.u.] + VI, ///< \f$V_\mathrm{i}\f$ Terminal-bus imaginary voltage [p.u.] + OMEGA, ///< \f$\omega\f$ Known machine speed deviation [p.u.] + VREF, ///< \f$V_{\mathrm{ref}}\f$ Unknown voltage-control reference [p.u.] + VS, ///< \f$V_S\f$ Known stabilizer input signal [p.u.] + VUEL, ///< \f$V_{\mathrm{UEL}}\f$ Known under-excitation limiter input [p.u.] }; /** @@ -74,10 +76,9 @@ namespace GridKit using Component::size_; using Component::tag_; using Component::variable_indices_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -103,6 +104,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -116,14 +118,14 @@ namespace GridKit __attribute__((always_inline)) inline int evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f); private: void initializeParameters(const ModelDataT& data); void initializeMonitor(); void setDerivedParameters(); + void gatherExternalVariables(); static __attribute__((always_inline)) inline ScalarT awmin( ScalarT x, diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp index d3ea94be6..65d325916 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp @@ -22,8 +22,6 @@ namespace GridKit * then assembles the resulting entries in COO form. * * @pre allocate() has completed. - * @pre evaluateResidual() has refreshed the interface buffers at the - * current state. * @pre Solver alpha and global variable and residual indices are set. */ template @@ -32,12 +30,13 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -48,62 +47,45 @@ namespace GridKit nnz_ = 0; - GridKit::Enzyme::Sparse::DfDy::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDyp::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDwb::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyp::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); return 0; diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 2645ed940..15ecb5722 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -108,13 +108,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.resize(2); - wb_.setToZero(); - - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); for (IdxT j = 0; j < size_; ++j) { @@ -438,65 +432,100 @@ namespace GridKit } /** - * @brief Residuals of system equations - * - * Refreshes the bus and signal interface buffers and evaluates the - * internal residual. ESDC1A injects no current, so there is no bus - * residual. An unattached input port falls back to the value latched - * by initialize(). + * @brief Gather external variables and index maps. * - * @return Zero on success. + * Unattached signal inputs retain the values latched by initialize(). */ template - int Esdc1a::evaluateResidual() + void Esdc1a::gatherExternalVariables() { - const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); - const auto VREF = static_cast(Esdc1aExternalVariables::VREF); - const auto VS = static_cast(Esdc1aExternalVariables::VS); - const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); - - auto* ws = ws_.getData(); + auto* y_ext = y_ext_.getData(); + + const auto VR_EXT = static_cast(Esdc1aExternalVariables::VR); + const auto VI_EXT = static_cast(Esdc1aExternalVariables::VI); + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VREF = static_cast(Esdc1aExternalVariables::VREF); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + + y_ext[VR_EXT] = Vr(); + y_ext[VI_EXT] = Vi(); + variable_indices_ext_[VR_EXT] = INVALID_INDEX; + variable_indices_ext_[VI_EXT] = INVALID_INDEX; + if (bus_->size() > 0) + { + variable_indices_ext_[VR_EXT] = bus_->getVariableIndex(0); + variable_indices_ext_[VI_EXT] = bus_->getVariableIndex(1); + } - ws[OMEGA] = omega_set_; - ws[VREF] = vref_set_; - ws[VS] = vs_set_; - ws[VUEL] = vuel_set_; - std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + y_ext[OMEGA] = omega_set_; + y_ext[VREF] = vref_set_; + y_ext[VS] = vs_set_; + y_ext[VUEL] = vuel_set_; + variable_indices_ext_[OMEGA] = INVALID_INDEX; + variable_indices_ext_[VREF] = INVALID_INDEX; + variable_indices_ext_[VS] = INVALID_INDEX; + variable_indices_ext_[VUEL] = INVALID_INDEX; if (auto omega_port = ports_.in.template port()) { - ws[OMEGA] = omega_port.readSignal(); - ws_indices_[OMEGA] = omega_port.signalVariableIndex(); + y_ext[OMEGA] = + ports_.in.template port().readSignal(); + variable_indices_ext_[OMEGA] = + ports_.in.template port().signalVariableIndex(); } if (auto vref_port = ports_.in.template port()) { - ws[VREF] = vref_port.readSignal(); - ws_indices_[VREF] = vref_port.signalVariableIndex(); + y_ext[VREF] = + ports_.in.template port().readSignal(); + variable_indices_ext_[VREF] = + ports_.in.template port().signalVariableIndex(); } if (auto vs_port = ports_.in.template port()) { - ws[VS] = vs_port.readSignal(); - ws_indices_[VS] = vs_port.signalVariableIndex(); + y_ext[VS] = ports_.in.template port().readSignal(); + variable_indices_ext_[VS] = + ports_.in.template port().signalVariableIndex(); } if (auto vuel_port = ports_.in.template port()) { - ws[VUEL] = vuel_port.readSignal(); - ws_indices_[VUEL] = vuel_port.signalVariableIndex(); + y_ext[VUEL] = + ports_.in.template port().readSignal(); + variable_indices_ext_[VUEL] = + ports_.in.template port().signalVariableIndex(); } + } - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + /** + * @brief Evaluate the internal ESDC1A residual equations. + */ + template + int Esdc1a::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + /** + * @brief Evaluate internal equations and external contributions. + * + * ESDC1A contributes no external residual, so the base implementation + * returns zero after the internal equations are evaluated. + */ + template + int Esdc1a::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + /** * @brief Access the monitor * @@ -520,9 +549,8 @@ namespace GridKit * * @param[in] y Internal variables in Esdc1aInternalVariables order. * @param[in] yp Internal derivatives in the same enum order. - * @param[in] wb Terminal-bus \f$(V_{\mathrm{r}},V_{\mathrm{i}})\f$ - * voltage components. - * @param[in] ws Signal values in Esdc1aExternalVariables order. + * @param[in] y_ext Terminal-bus voltage components and signal values in + * Esdc1aExternalVariables order. * @param[out] f Residuals in Esdc1aInternalVariables order. * @return Zero on success. */ @@ -531,8 +559,7 @@ namespace GridKit Esdc1a::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); @@ -547,10 +574,12 @@ namespace GridKit const auto VFE = static_cast(Esdc1aInternalVariables::VFE); const auto EFD = static_cast(Esdc1aInternalVariables::EFD); - const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); - const auto VREF = static_cast(Esdc1aExternalVariables::VREF); - const auto VS = static_cast(Esdc1aExternalVariables::VS); - const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + const auto VR_EXT = static_cast(Esdc1aExternalVariables::VR); + const auto VI_EXT = static_cast(Esdc1aExternalVariables::VI); + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VREF = static_cast(Esdc1aExternalVariables::VREF); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); const ScalarT efdp = y[EFDP]; const ScalarT vc = y[VC]; @@ -570,12 +599,13 @@ namespace GridKit const ScalarT vf_dot = yp[VF]; const ScalarT xll_dot = yp[XLL]; - const ScalarT omega = ws[OMEGA]; - const ScalarT vref = ws[VREF]; - const ScalarT vs = ws[VS]; - const ScalarT vuel = ws[VUEL]; + const ScalarT omega = y_ext[OMEGA]; + const ScalarT vref = y_ext[VREF]; + const ScalarT vs = y_ext[VS]; + const ScalarT vuel = y_ext[VUEL]; - const ScalarT ec = std::sqrt(wb[0] * wb[0] + wb[1] * wb[1]); + const ScalarT ec = std::sqrt( + y_ext[VR_EXT] * y_ext[VR_EXT] + y_ext[VI_EXT] * y_ext[VI_EXT]); const ScalarT ev_target = vref + vs + uel_on_ * vuel - vc - vf; const ScalarT vfe_target = Ke_eff_ * efdp + se; const ScalarT efdp_rate = (vr - vfe) / Te_; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index f45eef4f0..f78d000bb 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -61,9 +61,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::y_ext_; + using Component::variable_indices_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -91,6 +90,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -102,7 +102,7 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); private: static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); @@ -156,6 +156,8 @@ namespace GridKit /// Associate variable getter functions with enum values void initializeMonitor(); + + void gatherExternalVariables(); }; } // namespace Exciter diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp index adf3948a0..8676a377e 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp @@ -25,15 +25,16 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -42,65 +43,47 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus_->y().getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index ada67c9e3..19f03aada 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -104,14 +104,8 @@ namespace GridKit this->setResidualIndex(j, j); } - // Resize bus data - wb_.resize(2); - - // Resize signal variable data - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + // Resize coupling data + this->allocateExternalVectors(static_cast(Utilities::enum_size())); // Set output signals if (auto efd_port = ports_.out.template port()) @@ -366,19 +360,20 @@ namespace GridKit __attribute__((always_inline)) inline int Ieeet1::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { const auto OMEGA = static_cast(Ieeet1ExternalVariables::OMEGA); + const auto VREAL = static_cast(Ieeet1ExternalVariables::VREAL); + const auto VIMAG = static_cast(Ieeet1ExternalVariables::VIMAG); const auto VREF = static_cast(Ieeet1ExternalVariables::VREF); const auto VS = static_cast(Ieeet1ExternalVariables::VS); const auto VUEL = static_cast(Ieeet1ExternalVariables::VUEL); const auto VOEL = static_cast(Ieeet1ExternalVariables::VOEL); // Read bus voltage components - ScalarT vreal = wb[0]; - ScalarT vimag = wb[1]; + ScalarT vreal = y_ext[VREAL]; + ScalarT vimag = y_ext[VIMAG]; ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); // Read Internal Variables @@ -399,11 +394,11 @@ namespace GridKit ScalarT vfx_dot = yp[3]; // Set signal variable aliases - ScalarT omega = ws[OMEGA]; - ScalarT vref = ws[VREF]; - ScalarT vs = ws[VS]; - ScalarT vuel = ws[VUEL]; - ScalarT voel = ws[VOEL]; + ScalarT omega = y_ext[OMEGA]; + ScalarT vref = y_ext[VREF]; + ScalarT vs = y_ext[VS]; + ScalarT vuel = y_ext[VUEL]; + ScalarT voel = y_ext[VOEL]; // The 'pre-limit' derivative of Vr. ScalarT func = (-vr + Ka_ * vtr) / Ta_; @@ -425,25 +420,25 @@ namespace GridKit } /** - * @brief Residual evaluation + * @brief Gather external variables and index maps. * */ template - int Ieeet1::evaluateResidual() + void Ieeet1::gatherExternalVariables() { - auto* ws = ws_.getData(); + auto* y_ext = y_ext_.getData(); // Attached signals are read live; unattached ones keep the latched value. auto read_signal = [&](const ScalarT& latched) { const auto index = static_cast(variable); - ws[index] = latched; - ws_indices_[index] = INVALID_INDEX; + y_ext[index] = latched; + variable_indices_ext_[index] = INVALID_INDEX; if (auto port = ports_.in.template port()) { - ws[index] = port.readSignal(); - ws_indices_[index] = port.signalVariableIndex(); + y_ext[index] = port.readSignal(); + variable_indices_ext_[index] = port.signalVariableIndex(); } }; @@ -459,21 +454,44 @@ namespace GridKit Ieeet1ExternalVariables::VOEL>(voel_set_); // Bus voltages - auto* wb = wb_.getData(); - wb[0] = bus_->Vr(); - wb[1] = bus_->Vi(); + const auto VREAL = static_cast(Ieeet1ExternalVariables::VREAL); + const auto VIMAG = static_cast(Ieeet1ExternalVariables::VIMAG); + y_ext[VREAL] = bus_->Vr(); + y_ext[VIMAG] = bus_->Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[VREAL] = bus_->getVariableIndex(0); + variable_indices_ext_[VIMAG] = bus_->getVariableIndex(1); + } + } + + /** + * @brief Residual evaluation + * + */ + template + int Ieeet1::evaluateInternalResidual() + { + gatherExternalVariables(); // Residual evaluation const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + template + int Ieeet1::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + /** * @brief Initialization Exciter Parameters from data structure */ diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp index fbc3687e1..7972715bf 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp @@ -29,6 +29,8 @@ namespace GridKit /// External variables of a `SexsPti`. enum class SexsPtiExternalVariables : size_t { + VREAL, ///< Real bus voltage + VIMAG, ///< Imaginary bus voltage VREF, ///< Voltage reference VS, ///< Stabilizer output signal VUEL, ///< Under-excitation limiter signal @@ -48,9 +50,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::y_ext_; + using Component::variable_indices_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -78,6 +79,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -89,7 +91,7 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); private: BusT* bus_{nullptr}; @@ -118,6 +120,7 @@ namespace GridKit void initModelParams(const ModelDataT& data); void initializeMonitor(); + void gatherExternalVariables(); }; } // namespace Exciter diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp index a15c31c10..6820772c4 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp @@ -20,15 +20,16 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for SexsPti..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -37,65 +38,47 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus_->y().getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index 16c1e0cce..a726692c0 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -74,12 +74,7 @@ namespace GridKit this->setResidualIndex(j, j); } - wb_.resize(2); - - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); if (auto efd_port = ports_.out.template port()) { @@ -255,14 +250,15 @@ namespace GridKit __attribute__((always_inline)) inline int SexsPti::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { - const auto VREF = static_cast(SexsPtiExternalVariables::VREF); - const auto VS = static_cast(SexsPtiExternalVariables::VS); - const auto VUEL = static_cast(SexsPtiExternalVariables::VUEL); - const auto VOEL = static_cast(SexsPtiExternalVariables::VOEL); + const auto VREAL = static_cast(SexsPtiExternalVariables::VREAL); + const auto VIMAG = static_cast(SexsPtiExternalVariables::VIMAG); + const auto VREF = static_cast(SexsPtiExternalVariables::VREF); + const auto VS = static_cast(SexsPtiExternalVariables::VS); + const auto VUEL = static_cast(SexsPtiExternalVariables::VUEL); + const auto VOEL = static_cast(SexsPtiExternalVariables::VOEL); ScalarT vr = y[0]; ScalarT efd = y[1]; @@ -270,11 +266,11 @@ namespace GridKit ScalarT vr_dot = yp[0]; ScalarT efd_dot = yp[1]; - ScalarT Ec = std::sqrt(wb[0] * wb[0] + wb[1] * wb[1]); - ScalarT vref = ws[VREF]; - ScalarT vs = ws[VS]; - ScalarT vuel = ws[VUEL]; - ScalarT voel = ws[VOEL]; + ScalarT Ec = std::sqrt(y_ext[VREAL] * y_ext[VREAL] + y_ext[VIMAG] * y_ext[VIMAG]); + ScalarT vref = y_ext[VREF]; + ScalarT vs = y_ext[VS]; + ScalarT vuel = y_ext[VUEL]; + ScalarT voel = y_ext[VOEL]; ScalarT func = (-efd + (K_ / Tb_) * (-vr + Ta_ * vtr)) / Te_; @@ -285,22 +281,26 @@ namespace GridKit return 0; } + /** + * @brief Gather external variables and index maps. + * + */ template - int SexsPti::evaluateResidual() + void SexsPti::gatherExternalVariables() { - auto* ws = ws_.getData(); + auto* y_ext = y_ext_.getData(); // Attached signals are read live; unattached ones keep the latched value. auto read_signal = [&](const ScalarT& latched) { const auto index = static_cast(variable); - ws[index] = latched; - ws_indices_[index] = INVALID_INDEX; + y_ext[index] = latched; + variable_indices_ext_[index] = INVALID_INDEX; if (auto port = ports_.in.template port()) { - ws[index] = port.readSignal(); - ws_indices_[index] = port.signalVariableIndex(); + y_ext[index] = port.readSignal(); + variable_indices_ext_[index] = port.signalVariableIndex(); } }; @@ -313,20 +313,39 @@ namespace GridKit read_signal.template operator()(voel_set_); - auto* wb = wb_.getData(); - wb[0] = bus_->Vr(); - wb[1] = bus_->Vi(); + const auto VREAL = static_cast(SexsPtiExternalVariables::VREAL); + const auto VIMAG = static_cast(SexsPtiExternalVariables::VIMAG); + y_ext[VREAL] = bus_->Vr(); + y_ext[VIMAG] = bus_->Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[VREAL] = bus_->getVariableIndex(0); + variable_indices_ext_[VIMAG] = bus_->getVariableIndex(1); + } + } + + template + int SexsPti::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + template + int SexsPti::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + template void SexsPti::initModelParams(const ModelDataT& data) { diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp index 300218d25..f326ec285 100644 --- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp @@ -65,10 +65,9 @@ namespace GridKit using Component::va_component_base_; using Component::va_system_base_; using Component::variable_indices_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -92,6 +91,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -105,11 +105,11 @@ namespace GridKit [[gnu::always_inline]] inline int evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f); private: + void gatherExternalVariables(); void loadRealParameter(const ModelDataT& data, GastPtiParameters parameter, RealT& target, diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp index c41d5a7ff..fac739c81 100644 --- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp @@ -21,8 +21,6 @@ namespace GridKit * and linked signal variables, then constructs the model COO matrix. * * @pre allocate() has sized the model and Jacobian index maps. - * @pre evaluateResidual() has refreshed the current signal values and - * signal indices. * @pre The containing solver has set the current integration coefficient * and global variable/residual indices. */ @@ -32,11 +30,12 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for GastPti..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { auto size = static_cast(size_); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * signal_size; + auto buffer_size = 2 * size * size + size * static_cast(y_ext_.getSize()); J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -47,48 +46,45 @@ namespace GridKit nnz_ = 0; - GridKit::Enzyme::Sparse::DfDy::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDyp::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyp::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); return 0; diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp index d5a33b257..20b6dd9d0 100644 --- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp @@ -75,12 +75,11 @@ namespace GridKit /** * @brief Allocate the model vectors and wire the mechanical-power output * - * Sizes the state, residual, and signal-interface buffers, seeds the + * Sizes the state, residual, and external-variable buffers, seeds the * identity index maps, and points the assigned `pmech` node at the * internal state it publishes. That node aliases GASTPTI storage from * here on, which is how initialize() reads the seed the machine wrote. * Repeated allocation by this model reuses its existing link. - * GASTPTI attaches to no bus, so the bus-interface buffer stays empty. * */ template @@ -98,10 +97,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); for (IdxT j = 0; j < size_; ++j) { @@ -436,46 +432,72 @@ namespace GridKit } /** - * @brief Evaluate the seven GASTPTI-owned residual rows - * - * Refreshes the signal interface buffers and evaluates the internal - * residual. GASTPTI attaches to no bus, so there is no bus interface to - * refresh. An unattached reference port falls back to the value latched - * by initialize(); an unattached speed port reads zero deviation. + * @brief Gather external signal values and global indices. * + * An unattached reference port falls back to the value latched by + * initialize(); an unattached speed port reads zero deviation. */ template - int GastPti::evaluateResidual() + void GastPti::gatherExternalVariables() { + auto* y_ext = y_ext_.getData(); + const auto OMEGA = static_cast(GastPtiExternalVariables::OMEGA); const auto PREF = static_cast(GastPtiExternalVariables::PREF); - auto* ws = ws_.getData(); - - ws[OMEGA] = ZERO; - ws[PREF] = pref_set_; - std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + y_ext[OMEGA] = ZERO; + y_ext[PREF] = pref_set_; + std::fill(variable_indices_ext_.begin(), + variable_indices_ext_.end(), + INVALID_INDEX); if (auto omega_port = ports_.in.template port()) { - ws[OMEGA] = omega_port.readSignal(); - ws_indices_[OMEGA] = omega_port.signalVariableIndex(); + y_ext[OMEGA] = + ports_.in.template port().readSignal(); + variable_indices_ext_[OMEGA] = + ports_.in.template port().signalVariableIndex(); } if (auto pref_port = ports_.in.template port()) { - ws[PREF] = pref_port.readSignal(); - ws_indices_[PREF] = pref_port.signalVariableIndex(); + y_ext[PREF] = + ports_.in.template port().readSignal(); + variable_indices_ext_[PREF] = + ports_.in.template port().signalVariableIndex(); } + } + + /** + * @brief Evaluate the seven GASTPTI-owned residual rows. + */ + template + int GastPti::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, nullptr, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + /** + * @brief Evaluate internal equations and external contributions. + * + * GASTPTI contributes no external residual, so the base implementation + * returns zero after the internal equations are evaluated. + */ + template + int GastPti::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + + /** * @brief Access the monitor * @@ -500,10 +522,9 @@ namespace GridKit * @param[in] y Internal variables in `GastPtiInternalVariables` order; * each variable uses the base documented by its enum. * @param[in] yp Internal derivatives in the same enum order and bases. - * @param[in] wb Bus voltage components; unused, GASTPTI attaches to no bus. - * @param[in] ws External signals in `GastPtiExternalVariables` order: - * per-unit speed deviation followed by system-base - * active-power reference. + * @param[in] y_ext External signals in `GastPtiExternalVariables` order: + * per-unit speed deviation followed by system-base + * active-power reference. * @param[out] f Model-owned residuals in `GastPtiInternalVariables` order. */ template @@ -511,8 +532,7 @@ namespace GridKit GastPti::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { const auto XVALVE = static_cast(GastPtiInternalVariables::XVALVE); @@ -526,8 +546,6 @@ namespace GridKit const auto OMEGA = static_cast(GastPtiExternalVariables::OMEGA); const auto PREF = static_cast(GastPtiExternalVariables::PREF); - static_cast(wb); - const ScalarT xvalve = y[XVALVE]; const ScalarT xflow = y[XFLOW]; const ScalarT xtemp = y[XTEMP]; @@ -540,8 +558,8 @@ namespace GridKit const ScalarT xflow_dot = yp[XFLOW]; const ScalarT xtemp_dot = yp[XTEMP]; - const ScalarT omega = ws[OMEGA]; - const ScalarT pref = this->toComponentBase(ws[PREF]); + const ScalarT omega = y_ext[OMEGA]; + const ScalarT pref = this->toComponentBase(y_ext[PREF]); const ScalarT valve_target = Math::antiwindup(xvalve, vlv - xvalve, Vmin_response_, Vmax_response_); diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp index b11641318..0386ed6d0 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -76,10 +76,9 @@ namespace GridKit using Component::va_component_base_; using Component::va_system_base_; using Component::variable_indices_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -104,6 +103,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -117,11 +117,11 @@ namespace GridKit __attribute__((always_inline)) inline int evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f); private: + void gatherExternalVariables(); void initializeParameters(const ModelDataT& data); void initializeMonitor(); void setDerivedParameters(); diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp index af4cb2c90..5a1fe29bb 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp @@ -20,11 +20,12 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { auto size = static_cast(size_); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * signal_size; + auto buffer_size = 2 * size * size + size * static_cast(y_ext_.getSize()); J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -35,48 +36,45 @@ namespace GridKit nnz_ = 0; - GridKit::Enzyme::Sparse::DfDy::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDyp::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyp::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp index e85d05354..7df2e94f4 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -72,11 +72,10 @@ namespace GridKit /** * @brief Allocate the model vectors and wire the mechanical-power output * - * Sizes the state, residual, and signal-interface buffers, initializes + * Sizes the state, residual, and external-variable buffers, initializes * the identity index maps, and points the assigned `pmech` node at the * internal state it publishes. That node aliases HYGOV storage from * here on, which is how initialize() reads the machine value. - * HYGOV attaches to no bus, so the bus-interface buffer stays empty. * Repeated calls reuse the allocated vectors. * * @return int 0 on success. @@ -96,10 +95,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); for (IdxT j = 0; j < size_; ++j) { @@ -469,55 +465,77 @@ namespace GridKit } /** - * @brief Residuals of system equations + * @brief Gather external signal values and global indices. * - * Refreshes the signal interface buffers and evaluates the internal - * residual. HYGOV attaches to no bus, so there is no bus interface to - * refresh. An unattached reference or auxiliary port falls back to the - * value latched by initialize(); an unattached speed port reads zero - * deviation. - * - * @return int 0 on success. + * An unattached reference or auxiliary port falls back to the value + * latched by initialize(); an unattached speed port reads zero deviation. */ template - int Hygov::evaluateResidual() + void Hygov::gatherExternalVariables() { + auto* y_ext = y_ext_.getData(); + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); const auto PREF = static_cast(HygovExternalVariables::PREF); const auto PAUX = static_cast(HygovExternalVariables::PAUX); - auto* ws = ws_.getData(); - - ws[OMEGA] = ZERO; - ws[PREF] = pref_set_; - ws[PAUX] = paux_set_; - std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + y_ext[OMEGA] = ZERO; + y_ext[PREF] = pref_set_; + y_ext[PAUX] = paux_set_; + std::fill(variable_indices_ext_.begin(), + variable_indices_ext_.end(), + INVALID_INDEX); if (auto omega_port = ports_.in.template port()) { - ws[OMEGA] = omega_port.readSignal(); - ws_indices_[OMEGA] = omega_port.signalVariableIndex(); + y_ext[OMEGA] = ports_.in.template port().readSignal(); + variable_indices_ext_[OMEGA] = + ports_.in.template port().signalVariableIndex(); } if (auto pref_port = ports_.in.template port()) { - ws[PREF] = pref_port.readSignal(); - ws_indices_[PREF] = pref_port.signalVariableIndex(); + y_ext[PREF] = ports_.in.template port().readSignal(); + variable_indices_ext_[PREF] = + ports_.in.template port().signalVariableIndex(); } if (auto paux_port = ports_.in.template port()) { - ws[PAUX] = paux_port.readSignal(); - ws_indices_[PAUX] = paux_port.signalVariableIndex(); + y_ext[PAUX] = ports_.in.template port().readSignal(); + variable_indices_ext_[PAUX] = + ports_.in.template port().signalVariableIndex(); } + } + + /** + * @brief Evaluate the HYGOV-owned residual rows. + */ + template + int Hygov::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, nullptr, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + /** + * @brief Evaluate internal equations and external contributions. + * + * HYGOV contributes no external residual, so the base implementation + * returns zero after the internal equations are evaluated. + */ + template + int Hygov::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + /** * @brief Access the monitor * @@ -541,19 +559,17 @@ namespace GridKit * * @param[in] y Internal variables. * @param[in] yp Internal variable derivatives. - * @param[in] wb Bus voltage components; unused, HYGOV attaches to no bus. - * @param[in] ws External signal values on system base. + * @param[in] y_ext External signal values on system base. * @param[out] f Internal residuals. * @return int 0 on success. */ template __attribute__((always_inline)) inline int Hygov::evaluateInternalResidual( - const ScalarT* y, - const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - const ScalarT* ws, - ScalarT* f) + const ScalarT* y, + const ScalarT* yp, + const ScalarT* y_ext, + ScalarT* f) { const auto XN = static_cast(HygovInternalVariables::XN); const auto XF = static_cast(HygovInternalVariables::XF); @@ -591,9 +607,9 @@ namespace GridKit const ScalarT g_dot = yp[G]; const ScalarT q_dot = yp[Q]; - const ScalarT omega = ws[OMEGA]; - const ScalarT pref = ws[PREF]; - const ScalarT paux = ws[PAUX]; + const ScalarT omega = y_ext[OMEGA]; + const ScalarT pref = y_ext[PREF]; + const ScalarT paux = y_ext[PAUX]; const ScalarT yomega = xn + leadlag_gain_ * (omegadb - xn); diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index dcc5e147e..ab0ac8ae1 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -50,10 +50,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -82,6 +80,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; // Still to be implemented @@ -94,7 +93,7 @@ namespace GridKit public: __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); private: // Input parameters @@ -115,6 +114,7 @@ namespace GridKit // Parameter initialization function void initializeParameters(const ModelDataT& data); + void gatherExternalVariables(); void setDerivedParams(); static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp index 39e730e6b..f834653eb 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp @@ -25,14 +25,16 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * signal_size; + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -41,50 +43,47 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 229a58df3..7914853ad 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -196,11 +196,8 @@ namespace GridKit this->setResidualIndex(j, j); } - // Resize signal variable data - const auto signal_size = Utilities::enum_size(); - ws_.resize(static_cast(signal_size)); - ws_.setToZero(); - ws_indices_.assign(signal_size, INVALID_INDEX); + // Resize coupling data + this->allocateExternalVectors(static_cast(Utilities::enum_size())); // Set output signals if (auto pmech_port = ports_.out.template port()) @@ -349,11 +346,10 @@ namespace GridKit */ template __attribute__((always_inline)) inline int Tgov1::evaluateInternalResidual( - const ScalarT* y, - const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - const ScalarT* ws, - ScalarT* f) + const ScalarT* y, + const ScalarT* yp, + const ScalarT* y_ext, + ScalarT* f) { const auto PTX = static_cast(Tgov1InternalVariables::PTX); const auto PV = static_cast(Tgov1InternalVariables::PV); @@ -369,8 +365,8 @@ namespace GridKit const ScalarT pturb_dot = yp[PTX]; const ScalarT pv_dot = yp[PV]; - const ScalarT omega = ws[DELTAOMEGA]; - const ScalarT pref = ws[PREF]; + const ScalarT omega = y_ext[DELTAOMEGA]; + const ScalarT pref = y_ext[PREF]; f[PTX] = -pturb_dot - (pturb - pv - T2_ * pv_dot) / T3_; f[PV] = -pv_dot + Math::antiwindup(pv, -pv + (pref - omega) / R_, Pvmin_, Pvmax_) / T1_; @@ -383,40 +379,60 @@ namespace GridKit * @brief Residuals of system equations * */ + /** + * @brief Gather external variables and index maps. + * + */ template - int Tgov1::evaluateResidual() + void Tgov1::gatherExternalVariables() { + auto* y_ext = y_ext_.getData(); + const auto DELTAOMEGA = static_cast(Tgov1ExternalVariables::DELTAOMEGA); const auto PREF = static_cast(Tgov1ExternalVariables::PREF); - auto* ws = ws_.getData(); - - ws[DELTAOMEGA] = ScalarT{ZERO}; - ws[PREF] = pref_set_; - std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + y_ext[DELTAOMEGA] = ScalarT{ZERO}; + y_ext[PREF] = pref_set_; + std::fill(variable_indices_ext_.begin(), + variable_indices_ext_.end(), + INVALID_INDEX); if (auto speed_port = ports_.in.template port()) { - ws[DELTAOMEGA] = speed_port.readSignal(); - ws_indices_[DELTAOMEGA] = - speed_port.signalVariableIndex(); + y_ext[DELTAOMEGA] = ports_.in.template port().readSignal(); + variable_indices_ext_[DELTAOMEGA] = + ports_.in.template port().signalVariableIndex(); } if (auto pref_port = ports_.in.template port()) { - ws[PREF] = pref_port.readSignal(); - ws_indices_[PREF] = pref_port.signalVariableIndex(); + y_ext[PREF] = ports_.in.template port().readSignal(); + variable_indices_ext_[PREF] = + ports_.in.template port().signalVariableIndex(); } + } + + template + int Tgov1::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, nullptr, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + + template + int Tgov1::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } } // namespace Governor } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp index c27866a38..7a37202e9 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp @@ -9,6 +9,13 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `LoadZ` + enum class LoadZExternalVariables : size_t + { + VR, ///< \f$V_r\f$ + VI, ///< \f$V_i\f$ + }; + /*! * @brief Implementation of a constant load. * @@ -25,8 +32,10 @@ namespace GridKit using Component::yp_; using Component::abs_tol_; using Component::tag_; - using Component::wb_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; + using Component::residual_indices_ext_; + using Component::f_ext_; using Component::f_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; @@ -53,7 +62,9 @@ namespace GridKit virtual int initialize() override final; virtual int tagDifferentiable() override final; virtual int setAbsoluteTolerance(RealT) override final; + virtual int evaluateInternalResidual() override final; virtual int evaluateResidual() override final; + virtual int evaluateExternalResidual() override final; virtual int evaluateJacobian() override final; virtual int verify() const override final @@ -76,6 +87,7 @@ namespace GridKit private: void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -101,7 +113,7 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; public: - __attribute__((always_inline)) inline int evaluateBusResidual( + __attribute__((always_inline)) inline int evaluateExternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); __attribute__((always_inline)) inline int evaluateInternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp index 84a03966d..7f0083ced 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp @@ -23,14 +23,17 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for LoadZ..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto buffer_size = size * size + 2 * size * bus_size; + auto f_ext_size = static_cast(f_ext_.getSize()); + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = size * size + size * y_ext_size + size * f_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -46,39 +49,39 @@ namespace GridKit (this->getVariableIndices()).data(), y_.getData(), yp_.getData(), - wb_.getData(), + y_ext_.getData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz_); - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus_->y().getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, - static_cast(bus_->size()), - static_cast(y_.getSize()), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp index cea110760..17c976065 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -102,8 +103,10 @@ namespace GridKit } // Resize coupling data - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(2); + f_ext_.setToZero(); + residual_indices_ext_.assign(2, INVALID_INDEX); allocated_ = true; return 0; @@ -167,20 +170,20 @@ namespace GridKit } /** - * @brief Bus residual + * @brief External residual * */ template - __attribute__((always_inline)) int LoadZ::evaluateBusResidual( + __attribute__((always_inline)) int LoadZ::evaluateExternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - ScalarT* h) + [[maybe_unused]] const ScalarT* y_ext, + ScalarT* f_ext) { const ScalarT Ir = y[0]; const ScalarT Ii = y[1]; - h[0] = Ir; - h[1] = Ii; + f_ext[0] = Ir; + f_ext[1] = Ii; return 0; } @@ -193,11 +196,11 @@ namespace GridKit __attribute__((always_inline)) int LoadZ::evaluateInternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - const ScalarT* wb, + const ScalarT* y_ext, ScalarT* f) { - const ScalarT Vr = wb[0]; - const ScalarT Vi = wb[1]; + const ScalarT Vr = y_ext[0]; + const ScalarT Vi = y_ext[1]; const ScalarT Ir = y[0]; const ScalarT Ii = y[1]; f[0] = Ir + g_ * Vr - b_ * Vi; @@ -207,33 +210,73 @@ namespace GridKit } /** - * @brief Residual contribution of the load is pushed to the bus. + * @brief Gather external variables and index maps. * */ template - int LoadZ::evaluateResidual() + void LoadZ::gatherExternalVariables() { - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + auto* y_ext = y_ext_.getData(); + + y_ext[0] = Vr(); + y_ext[1] = Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[0] = bus_->getVariableIndex(0); + variable_indices_ext_[1] = bus_->getVariableIndex(1); + residual_indices_ext_[0] = bus_->getResidualIndex(0); + residual_indices_ext_[1] = bus_->getResidualIndex(1); + } + } + + template + int LoadZ::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - auto* h = h_.getData(); - evaluateInternalResidual(y, yp, wb, f); - evaluateBusResidual(y, yp, wb, h); - Ir() += h[0]; - Ii() += h[1]; + evaluateInternalResidual(y, yp, y_ext_.getData(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * @brief External residual contributions to the bus. + * + */ + template + int LoadZ::evaluateExternalResidual() + { + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext, f_ext); + Ir() += f_ext[0]; + Ii() += f_ext[1]; if (bus_->size() > 0) { bus_->getResidual().setDataUpdated(); } - f_.setDataUpdated(); return 0; } + /** + * @brief Residual contribution of the load is pushed to the bus. + * + */ + template + int LoadZ::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } + /** * @brief Derived parameters * diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp index 737cba8c7..8cfb774f3 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp @@ -9,6 +9,13 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `LoadZIP` + enum class LoadZIPExternalVariables : size_t + { + VR, ///< \f$V_r\f$ + VI, ///< \f$V_i\f$ + }; + /*! * @brief Implementation of a ZIP load. * @@ -25,8 +32,10 @@ namespace GridKit using Component::yp_; using Component::abs_tol_; using Component::tag_; - using Component::wb_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; + using Component::residual_indices_ext_; + using Component::f_ext_; using Component::f_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; @@ -53,7 +62,9 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; + int evaluateExternalResidual() override final; int evaluateJacobian() override final; int verify() const override final @@ -89,6 +100,7 @@ namespace GridKit private: void initializeParameters(const ModelDataT& data); void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -114,7 +126,7 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; public: - __attribute__((always_inline)) inline int evaluateBusResidual( + __attribute__((always_inline)) inline int evaluateExternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); __attribute__((always_inline)) inline int evaluateInternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp index 4e25242f7..f78a284dd 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp @@ -1,4 +1,3 @@ - #include #include "LoadZIPImpl.hpp" @@ -18,14 +17,17 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for LoadZIP..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto buffer_size = size * size + 2 * size * bus_size; + auto f_ext_size = static_cast(f_ext_.getSize()); + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = size * size + size * y_ext_size + size * f_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -41,39 +43,39 @@ namespace GridKit (this->getVariableIndices()).data(), y_.getData(), yp_.getData(), - wb_.getData(), + y_ext_.getData(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz_); - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - bus_->y().getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, - static_cast(bus_->size()), - static_cast(y_.getSize()), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp index 19bd6ed6c..b6f318d0a 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -115,8 +116,10 @@ namespace GridKit } // Resize coupling data - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(2); + f_ext_.setToZero(); + residual_indices_ext_.assign(2, INVALID_INDEX); allocated_ = true; return 0; @@ -188,52 +191,96 @@ namespace GridKit } /** - * @brief Bus residual + * @brief External residual * */ template - __attribute__((always_inline)) int LoadZIP::evaluateBusResidual( + __attribute__((always_inline)) int LoadZIP::evaluateExternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - ScalarT* h) + [[maybe_unused]] const ScalarT* y_ext, + ScalarT* f_ext) { const ScalarT Ir = y[0]; const ScalarT Ii = y[1]; - h[0] = Ir; - h[1] = Ii; + f_ext[0] = Ir; + f_ext[1] = Ii; return 0; } /** - * @brief Residual contribution of the load is pushed to the bus. + * @brief Gather external variables and index maps. * */ template - int LoadZIP::evaluateResidual() + void LoadZIP::gatherExternalVariables() + { + auto* y_ext = y_ext_.getData(); + + y_ext[0] = Vr(); + y_ext[1] = Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[0] = bus_->getVariableIndex(0); + variable_indices_ext_[1] = bus_->getVariableIndex(1); + residual_indices_ext_[0] = bus_->getResidualIndex(0); + residual_indices_ext_[1] = bus_->getResidualIndex(1); + } + } + + /** + * @brief Internal residual for the load model. + * + */ + template + int LoadZIP::evaluateInternalResidual() { - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - auto* h = h_.getData(); - evaluateInternalResidual(y, yp, wb, f); - evaluateBusResidual(y, yp, wb, h); - Ir() += h[0]; - Ii() += h[1]; + evaluateInternalResidual(y, yp, y_ext_.getData(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * @brief External residual contributions to the bus. + * + */ + template + int LoadZIP::evaluateExternalResidual() + { + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext, f_ext); + Ir() += f_ext[0]; + Ii() += f_ext[1]; if (bus_->size() > 0) { bus_->getResidual().setDataUpdated(); } - f_.setDataUpdated(); return 0; } + /** + * @brief Residual contribution of the load is pushed to the bus. + * + */ + template + int LoadZIP::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } + /** * @brief Internal residual * @@ -242,11 +289,11 @@ namespace GridKit __attribute__((always_inline)) int LoadZIP::evaluateInternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - const ScalarT* wb, + const ScalarT* y_ext, ScalarT* f) { - const ScalarT Vr = wb[0]; - const ScalarT Vi = wb[1]; + const ScalarT Vr = y_ext[0]; + const ScalarT Vi = y_ext[1]; const ScalarT Ir = y[0]; const ScalarT Ii = y[1]; const RealT Vnom2 = Vnom_ * Vnom_; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index 448e0144c..8c2321527 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -54,10 +54,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -84,6 +82,7 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -98,7 +97,6 @@ namespace GridKit const ScalarT*, const ScalarT*, const ScalarT*, - const ScalarT*, ScalarT*); private: @@ -149,6 +147,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); void initializeMonitor(); + void gatherExternalVariables(); }; } // namespace Stabilizer diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp index c814afe3e..29f9db3ca 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp @@ -27,14 +27,16 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Ieeest..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * signal_size; + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -43,50 +45,47 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - nullptr, - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index 1d7d6d9b0..d560f2c8f 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -173,10 +173,7 @@ namespace GridKit this->setResidualIndex(j, j); } - ws_.resize(1); - ws_.setToZero(); - ws_indices_.resize(1); - ws_indices_[0] = INVALID_INDEX; + this->allocateExternalVectors(static_cast(Utilities::enum_size())); if (auto output_port = ports_.out.template port()) { @@ -225,8 +222,9 @@ namespace GridKit auto input_port = ports_.in.template port(); const ScalarT u = input_port.readSignal(); - auto* y = y_.getData(); - auto* yp = yp_.getData(); + auto* y = y_.getData(); + auto* yp = yp_.getData(); + auto* y_ext = y_ext_.getData(); for (IdxT i = 0; i < size_; ++i) { @@ -234,8 +232,9 @@ namespace GridKit yp[static_cast(i)] = 0.0; } - ws_.getData()[0] = u; - ws_indices_[0] = input_port.signalVariableIndex(); + y_ext[0] = u; + variable_indices_ext_[0] = + ports_.in.template port().signalVariableIndex(); y[0] = use_notch_ * u; y[4] = u; @@ -295,11 +294,10 @@ namespace GridKit template __attribute__((always_inline)) inline int Ieeest::evaluateInternalResidual( - const ScalarT* y, - const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - const ScalarT* ws, - ScalarT* f) + const ScalarT* y, + const ScalarT* yp, + const ScalarT* y_ext, + ScalarT* f) { ScalarT x1 = y[0]; ScalarT x2 = y[1]; @@ -322,7 +320,7 @@ namespace GridKit ScalarT x6_dot = yp[5]; ScalarT x7_dot = yp[6]; - ScalarT u = ws[0]; + ScalarT u = y_ext[0]; const ScalarT x2_rhs = (use_4th_order_ + use_3rd_order_) * x3 + use_2nd_order_ * (-a0_ * x1 - a1_ * x2 + u) * safe_inv_a2_; @@ -344,26 +342,43 @@ namespace GridKit return 0; } + /** + * @brief Gather external variables and index maps. + * + */ template - int Ieeest::evaluateResidual() + void Ieeest::gatherExternalVariables() { - auto* ws = ws_.getData(); + auto* y_ext = y_ext_.getData(); if (auto input_port = ports_.in.template port()) { - ws[0] = input_port.readSignal(); - ws_indices_[0] = input_port.signalVariableIndex(); + y_ext[0] = ports_.in.template port().readSignal(); + variable_indices_ext_[0] = ports_.in.template port().signalVariableIndex(); } + } + + template + int Ieeest::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, nullptr, ws, f); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); f_.setDataUpdated(); return 0; } + template + int Ieeest::evaluateResidual() + { + evaluateInternalResidual(); + return this->evaluateExternalResidual(); + } + template const Model::VariableMonitorBase* Ieeest::getMonitor() const { diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp index e02704ca5..d902f35e4 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp @@ -63,10 +63,10 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; + using Component::residual_indices_ext_; + using Component::f_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -125,7 +125,9 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; + int evaluateExternalResidual() override final; // Still to be implemented int evaluateJacobian() override final; @@ -141,6 +143,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); /// Associate variable getter functions with enum values void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -165,8 +168,8 @@ namespace GridKit public: __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); - __attribute__((always_inline)) inline int evaluateBusResidual( + const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + __attribute__((always_inline)) inline int evaluateExternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); private: diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp index 242d92eaf..d7194d6a8 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp @@ -23,15 +23,17 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * signal_size + 2 * size * bus_size; + auto f_ext_size = static_cast(f_ext_.getSize()); + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size + size * f_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -40,79 +42,61 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, - static_cast(bus_->size()), - static_cast(y_.getSize()), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp index 2051a9214..e58105bb1 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp @@ -325,15 +325,11 @@ namespace GridKit this->setResidualIndex(j, j); } - // Resize bus data - wb_.resize(2); - h_.resize(2); - - // Resize signal variable data - ws_.resize(2); - ws_indices_.resize(2); - ws_indices_[0] = INVALID_INDEX; - ws_indices_[1] = INVALID_INDEX; + // Resize coupling data + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(2); + f_ext_.setToZero(); + residual_indices_ext_.assign(2, INVALID_INDEX); // Set output signals if (auto speed_port = ports_.out.template port()) @@ -500,8 +496,7 @@ namespace GridKit __attribute__((always_inline)) inline int Genrou::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { /* Read variables */ @@ -532,12 +527,12 @@ namespace GridKit ScalarT Edp_dot = yp[5]; // Set coupling variable aliases - ScalarT vr = wb[0]; - ScalarT vi = wb[1]; + ScalarT vr = y_ext[0]; + ScalarT vi = y_ext[1]; // Set signal variable aliases - ScalarT pmech = this->toComponentBase(ws[0]); - ScalarT efd = ws[1]; + ScalarT pmech = this->toComponentBase(y_ext[2]); + ScalarT efd = y_ext[3]; static constexpr auto pi = std::numbers::pi_v; @@ -566,77 +561,118 @@ namespace GridKit } /** - * @brief Bus residual + * @brief External residual * */ template - __attribute__((always_inline)) inline int Genrou::evaluateBusResidual( + __attribute__((always_inline)) inline int Genrou::evaluateExternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - ScalarT* h) + [[maybe_unused]] const ScalarT* y_ext, + ScalarT* f_ext) { ScalarT ir = y[15]; ScalarT ii = y[16]; // Convert current injection to system base for the network. - h[0] = this->toSystemBase(ir); - h[1] = this->toSystemBase(ii); + f_ext[0] = this->toSystemBase(ir); + f_ext[1] = this->toSystemBase(ii); return 0; } /** - * \brief Residual evaluation and contribution to the connected bus + * @brief Gather external variables and index maps. * */ template - int Genrou::evaluateResidual() + void Genrou::gatherExternalVariables() { - auto* ws = ws_.getData(); + auto* y_ext = y_ext_.getData(); + + // Bus voltages + y_ext[0] = Vr(); + y_ext[1] = Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[0] = bus_->getVariableIndex(0); + variable_indices_ext_[1] = bus_->getVariableIndex(1); + residual_indices_ext_[0] = bus_->getResidualIndex(0); + residual_indices_ext_[1] = bus_->getResidualIndex(1); + } // Mechanical Power - ws[0] = pmech_set_; - if (auto pmech_port = ports_.in.template port()) + y_ext[2] = pmech_set_; + if (ports_.in.template port().connected()) { - ws[0] = pmech_port.readSignal(); - ws_indices_[0] = pmech_port.signalVariableIndex(); + y_ext[2] = ports_.in.template port().readSignal(); + variable_indices_ext_[2] = ports_.in.template port().signalVariableIndex(); } // Exciter Efield - ws[1] = efd_set_; - if (auto efd_port = ports_.in.template port()) + y_ext[3] = efd_set_; + if (ports_.in.template port().connected()) { - ws[1] = efd_port.readSignal(); - ws_indices_[1] = efd_port.signalVariableIndex(); + y_ext[3] = ports_.in.template port().readSignal(); + variable_indices_ext_[3] = ports_.in.template port().signalVariableIndex(); } + } - // Bus voltages - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + /** + * \brief Internal residual for the generator model. + * + */ + template + int Genrou::evaluateInternalResidual() + { + gatherExternalVariables(); - // Residual evaluation const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - auto* h = h_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); - evaluateBusResidual(y, yp, wb, h); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * \brief External residual contributions to the bus. + * + */ + template + int Genrou::evaluateExternalResidual() + { + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext, f_ext); // Genrou contribution to bus algebraic equations - Ir() += h[0]; - Ii() += h[1]; + Ir() += f_ext[0]; + Ii() += f_ext[1]; if (bus_->size() > 0) { bus_->getResidual().setDataUpdated(); } - f_.setDataUpdated(); return 0; } + /** + * \brief Residual evaluation and contribution to the connected bus + * + */ + template + int Genrou::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } + template void Genrou::setDerivedParams() { diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp index cf3d379c6..b1898188e 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp @@ -59,10 +59,10 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; + using Component::residual_indices_ext_; + using Component::f_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -91,7 +91,9 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; + int evaluateExternalResidual() override final; // Still to be implemented int evaluateJacobian() override final; @@ -107,6 +109,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); /// Associate variable getter functions with enum values void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -131,8 +134,8 @@ namespace GridKit public: __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); - __attribute__((always_inline)) inline int evaluateBusResidual( + const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + __attribute__((always_inline)) inline int evaluateExternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); private: diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp index 810a725c8..bdf06cc46 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp @@ -23,15 +23,17 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Gensal..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * signal_size + 2 * size * bus_size; + auto f_ext_size = static_cast(f_ext_.getSize()); + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size + size * f_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -40,79 +42,61 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, - static_cast(bus_->size()), - static_cast(y_.getSize()), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp index e656a2f6b..8268711fe 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp @@ -201,15 +201,11 @@ namespace GridKit this->setResidualIndex(j, j); } - // Resize bus data - wb_.resize(2); - h_.resize(2); - - // Resize signal variable data - ws_.resize(2); - ws_indices_.resize(2); - ws_indices_[0] = INVALID_INDEX; - ws_indices_[1] = INVALID_INDEX; + // Resize coupling data + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(2); + f_ext_.setToZero(); + residual_indices_ext_.assign(2, INVALID_INDEX); // Set output signals if (auto speed_port = ports_.out.template port()) @@ -360,8 +356,7 @@ namespace GridKit __attribute__((always_inline)) inline int Gensal::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { /* Read variables */ @@ -388,12 +383,12 @@ namespace GridKit ScalarT psiqpp_dot = yp[4]; // Set coupling variable aliases - ScalarT vr = wb[0]; - ScalarT vi = wb[1]; + ScalarT vr = y_ext[0]; + ScalarT vi = y_ext[1]; // Set signal variable aliases - ScalarT pmech = this->toComponentBase(ws[0]); - ScalarT efd = ws[1]; + ScalarT pmech = this->toComponentBase(y_ext[2]); + ScalarT efd = y_ext[3]; static constexpr auto pi = std::numbers::pi_v; @@ -419,77 +414,118 @@ namespace GridKit } /** - * @brief Bus residual + * @brief External residual * */ template - __attribute__((always_inline)) inline int Gensal::evaluateBusResidual( + __attribute__((always_inline)) inline int Gensal::evaluateExternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - ScalarT* h) + [[maybe_unused]] const ScalarT* y_ext, + ScalarT* f_ext) { ScalarT ir = y[12]; ScalarT ii = y[13]; // Convert current injection to system base for the network. - h[0] = this->toSystemBase(ir); - h[1] = this->toSystemBase(ii); + f_ext[0] = this->toSystemBase(ir); + f_ext[1] = this->toSystemBase(ii); return 0; } /** - * \brief Residual evaluation and contribution to the connected bus + * @brief Gather external variables and index maps. * */ template - int Gensal::evaluateResidual() + void Gensal::gatherExternalVariables() { - auto* ws = ws_.getData(); + auto* y_ext = y_ext_.getData(); + + // Bus voltages + y_ext[0] = Vr(); + y_ext[1] = Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[0] = bus_->getVariableIndex(0); + variable_indices_ext_[1] = bus_->getVariableIndex(1); + residual_indices_ext_[0] = bus_->getResidualIndex(0); + residual_indices_ext_[1] = bus_->getResidualIndex(1); + } // Mechanical Power - ws[0] = pmech_set_; - if (auto pmech_port = ports_.in.template port()) + y_ext[2] = pmech_set_; + if (ports_.in.template port().connected()) { - ws[0] = pmech_port.readSignal(); - ws_indices_[0] = pmech_port.signalVariableIndex(); + y_ext[2] = ports_.in.template port().readSignal(); + variable_indices_ext_[2] = ports_.in.template port().signalVariableIndex(); } // Exciter Efield - ws[1] = efd_set_; - if (auto efd_port = ports_.in.template port()) + y_ext[3] = efd_set_; + if (ports_.in.template port().connected()) { - ws[1] = efd_port.readSignal(); - ws_indices_[1] = efd_port.signalVariableIndex(); + y_ext[3] = ports_.in.template port().readSignal(); + variable_indices_ext_[3] = ports_.in.template port().signalVariableIndex(); } + } - // Bus voltages - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + /** + * \brief Internal residual for the generator model. + * + */ + template + int Gensal::evaluateInternalResidual() + { + gatherExternalVariables(); - // Residual evaluation const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - auto* h = h_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); - evaluateBusResidual(y, yp, wb, h); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * \brief External residual contributions to the bus. + * + */ + template + int Gensal::evaluateExternalResidual() + { + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext, f_ext); // Gensal contribution to bus algebraic equations - Ir() += h[0]; - Ii() += h[1]; + Ir() += f_ext[0]; + Ii() += f_ext[1]; if (bus_->size() > 0) { bus_->getResidual().setDataUpdated(); } - f_.setDataUpdated(); return 0; } + /** + * \brief Residual evaluation and contribution to the connected bus + * + */ + template + int Gensal::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } + template void Gensal::setDerivedParams() { diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index ced9a652e..bbeb06c6e 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -50,10 +50,10 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::ws_; - using Component::ws_indices_; - using Component::h_; + using Component::y_ext_; + using Component::variable_indices_ext_; + using Component::residual_indices_ext_; + using Component::f_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -80,7 +80,9 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateInternalResidual() override final; int evaluateResidual() override final; + int evaluateExternalResidual() override final; // Still to be implemented int evaluateJacobian() override final; @@ -96,6 +98,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); /// Associate variable getter functions with enum values void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -120,8 +123,8 @@ namespace GridKit public: __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); - __attribute__((always_inline)) inline int evaluateBusResidual( + const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + __attribute__((always_inline)) inline int evaluateExternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); private: diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp index 0f8ea63ec..bd2dd64a1 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp @@ -23,15 +23,17 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + gatherExternalVariables(); + if (J_rows_buffer_ == nullptr) { // Reserve space for the dense blocks. // The size of the buffer is the sum of maximum capacities of the blocks. // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); - auto bus_size = static_cast(bus_->size()); - auto signal_size = static_cast(ws_.getSize()); - auto buffer_size = 2 * size * size + size * signal_size + 2 * size * bus_size; + auto f_ext_size = static_cast(f_ext_.getSize()); + auto y_ext_size = static_cast(y_ext_.getSize()); + auto buffer_size = 2 * size * size + size * y_ext_size + size * f_ext_size; J_rows_buffer_ = new IdxT[buffer_size]; J_cols_buffer_ = new IdxT[buffer_size]; J_vals_buffer_ = new RealT[buffer_size]; @@ -40,79 +42,61 @@ namespace GridKit nnz_ = 0; GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); - - GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(ws_.getSize()), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - ws_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + static_cast(y_ext_.getSize()), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, - static_cast(bus_->size()), - static_cast(y_.getSize()), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.getData(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + static_cast(f_ext_.getSize()), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.getData(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index cd61c29ea..abf147b68 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -142,15 +142,11 @@ namespace GridKit this->setResidualIndex(j, j); } - // Resize bus data - wb_.resize(2); - h_.resize(2); - - // Resize signal variable data - ws_.resize(2); - ws_indices_.resize(2); - ws_indices_[0] = INVALID_INDEX; - ws_indices_[1] = INVALID_INDEX; + // Resize coupling data + this->allocateExternalVectors(static_cast(Utilities::enum_size())); + f_ext_.resize(2); + f_ext_.setToZero(); + residual_indices_ext_.assign(2, INVALID_INDEX); // Set output signals if (auto speed_port = ports_.out.template port()) @@ -283,8 +279,7 @@ namespace GridKit __attribute__((always_inline)) inline int GenClassical::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, + const ScalarT* y_ext, ScalarT* f) { // Set variable aliases for better readability. @@ -299,12 +294,12 @@ namespace GridKit const ScalarT omega_dot = yp[1]; // Set coupling variable aliases - const ScalarT vr = wb[0]; - const ScalarT vi = wb[1]; + const ScalarT vr = y_ext[0]; + const ScalarT vi = y_ext[1]; // Set signal variable aliases - const ScalarT pmech = this->toComponentBase(ws[0]); - const ScalarT efd = ws[1]; + const ScalarT pmech = this->toComponentBase(y_ext[2]); + const ScalarT efd = y_ext[3]; static constexpr auto pi = std::numbers::pi_v; @@ -321,70 +316,116 @@ namespace GridKit } /** - * @brief Bus residual + * @brief External residual * */ template - __attribute__((always_inline)) inline int GenClassical::evaluateBusResidual( + __attribute__((always_inline)) inline int GenClassical::evaluateExternalResidual( const ScalarT* y, [[maybe_unused]] const ScalarT* yp, - [[maybe_unused]] const ScalarT* wb, - ScalarT* h) + [[maybe_unused]] const ScalarT* y_ext, + ScalarT* f_ext) { const ScalarT ir = y[3]; const ScalarT ii = y[4]; - h[0] = this->toSystemBase(ir); - h[1] = this->toSystemBase(ii); + f_ext[0] = this->toSystemBase(ir); + f_ext[1] = this->toSystemBase(ii); return 0; } /** - * \brief Residual for the generator model. + * @brief Gather external variables and index maps. * */ template - int GenClassical::evaluateResidual() + void GenClassical::gatherExternalVariables() { - auto* ws = ws_.getData(); + auto* y_ext = y_ext_.getData(); - ws[0] = pmech_set_; - if (auto pmech_port = ports_.in.template port()) + // Bus voltages + y_ext[0] = Vr(); + y_ext[1] = Vi(); + if (bus_->size() > 0) { - ws[0] = pmech_port.readSignal(); - ws_indices_[0] = pmech_port.signalVariableIndex(); + variable_indices_ext_[0] = bus_->getVariableIndex(0); + variable_indices_ext_[1] = bus_->getVariableIndex(1); + residual_indices_ext_[0] = bus_->getResidualIndex(0); + residual_indices_ext_[1] = bus_->getResidualIndex(1); } - ws[1] = efd_set_; - if (auto efd_port = ports_.in.template port()) + // Mechanical Power + y_ext[2] = pmech_set_; + if (ports_.in.template port().connected()) { - ws[1] = efd_port.readSignal(); - ws_indices_[1] = efd_port.signalVariableIndex(); + y_ext[2] = ports_.in.template port().readSignal(); + variable_indices_ext_[2] = ports_.in.template port().signalVariableIndex(); } - auto* wb = wb_.getData(); - wb[0] = Vr(); - wb[1] = Vi(); + // Exciter Efield + y_ext[3] = efd_set_; + if (ports_.in.template port().connected()) + { + y_ext[3] = ports_.in.template port().readSignal(); + variable_indices_ext_[3] = ports_.in.template port().signalVariableIndex(); + } + } + + /** + * \brief Internal residual for the generator model. + * + */ + template + int GenClassical::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - auto* h = h_.getData(); - evaluateInternalResidual(y, yp, wb, ws, f); - evaluateBusResidual(y, yp, wb, h); + evaluateInternalResidual(y, yp, y_ext_.getData(), f); + f_.setDataUpdated(); + + return 0; + } - Ir() += h[0]; - Ii() += h[1]; + /** + * \brief External residual contributions to the bus. + * + */ + template + int GenClassical::evaluateExternalResidual() + { + auto* y_ext = y_ext_.getData(); + auto* f_ext = f_ext_.getData(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext, f_ext); + + // GenClassical contribution to bus algebraic equations + Ir() += f_ext[0]; + Ii() += f_ext[1]; if (bus_->size() > 0) { bus_->getResidual().setDataUpdated(); } - f_.setDataUpdated(); return 0; } + /** + * \brief Residual evaluation and contribution to the connected bus + * + */ + template + int GenClassical::evaluateResidual() + { + evaluateInternalResidual(); + return evaluateExternalResidual(); + } + template void GenClassical::setDerivedParams() { diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index c39e0894b..4cf420280 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -92,6 +92,8 @@ namespace GridKit int tagDifferentiable() override; int setAbsoluteTolerance(RealT rel_tol) override; + int evaluateInternalResidual() override; + int evaluateExternalResidual() override; int evaluateResidual() override; int evaluateJacobian() override; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 2abc4f153..d44a78e23 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -673,20 +673,14 @@ namespace GridKit } /** - * @brief Compute system residual vector + * @brief Compute the residuals each bus and component owns. * * Buses and components read and write their bound system-vector slices - * directly. - * - * @warning Residuals must be computed for buses, before component - * residuals are computed. Buses own residuals for currents - * Ir and Ii, but the contributions to these residuals come - * from components. Buses assign their residual values, while components - * add to those values by in-place adition. This is why (for now) bus - * residuals need to be computed first. + * directly. Buses assign their residual values first, so that external + * residuals can accumulate into them. */ template - int SystemModel::evaluateResidual() + int SystemModel::evaluateInternalResidual() { for (const auto& bus : buses_) { @@ -695,9 +689,39 @@ namespace GridKit for (const auto& component : components_) { - component->evaluateResidual(); + component->evaluateInternalResidual(); } + return 0; + } + + /** + * @brief Accumulate component contributions to residuals owned elsewhere, + * e.g. bus current balances. + */ + template + int SystemModel::evaluateExternalResidual() + { + for (const auto& component : components_) + { + component->evaluateExternalResidual(); + } + + return 0; + } + + /** + * @brief Compute system residual vector + * + * Internal residuals assign every owned entry of the residual vector, + * then external residuals accumulate the remaining contributions. + */ + template + int SystemModel::evaluateResidual() + { + evaluateInternalResidual(); + evaluateExternalResidual(); + f_.setDataUpdated(); return 0; diff --git a/tests/UnitTests/PhasorDynamics/BranchTests.hpp b/tests/UnitTests/PhasorDynamics/BranchTests.hpp index ba09fd4e5..2d6f9e5ba 100644 --- a/tests/UnitTests/PhasorDynamics/BranchTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BranchTests.hpp @@ -269,6 +269,7 @@ namespace GridKit ref_bus2.initialize(); ref_bus2.evaluateResidual(); PhasorDynamics::Branch ref_branch(&ref_bus1, &ref_bus2, R, X, G, B, tap, phase); + ref_branch.allocate(); PhasorDynamics::Bus test_bus1(Vr1, Vi1); PhasorDynamics::Bus test_bus2(Vr2, Vi2); @@ -279,6 +280,7 @@ namespace GridKit test_bus2.initialize(); test_bus2.evaluateResidual(); PhasorDynamics::Branch test_branch(&test_bus1, &test_bus2, 1.0, 1.0, 0.0, 0.0); + test_branch.allocate(); test_branch.setR(R); test_branch.setX(X); @@ -378,6 +380,8 @@ namespace GridKit PhasorDynamics::Branch data_branch(&data_bus1, &data_bus2, data); PhasorDynamics::Branch ref_branch(&ref_bus1, &ref_bus2, R, X, G, B, 1.0, 0.0); + data_branch.allocate(); + ref_branch.allocate(); data_branch.evaluateResidual(); ref_branch.evaluateResidual(); diff --git a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp index 03528173e..0d4c4b77e 100644 --- a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -100,19 +101,31 @@ namespace GridKit // Jacobian via Enzyme auto enzyme_jacobian = EnzymeJacobian(R, X, status); - if (!status) + success *= dependency_tracking_jacobian.size() == enzyme_jacobian.size(); + + const auto remove_zeros = [](auto& jacobian) { - // HACK: Enzyme retains the fixed DfDwb/DhDy structure and masks its - // inactive values to exact zero, while DependencyTracking omits them. - for (auto& row : enzyme_jacobian) + for (auto& row : jacobian) { - std::erase_if(row, [](const auto& entry) - { return entry.second == 0.0; }); + for (auto entry = row.begin(); entry != row.end();) + { + if (entry->second == 0.0) + { + entry = row.erase(entry); + } + else + { + ++entry; + } + } } - } + }; + remove_zeros(dependency_tracking_jacobian); + remove_zeros(enzyme_jacobian); /// Compare DependencyTracking dependencies to Enzyme's - for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) + const size_t rows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); + for (size_t i = 0; i < rows; ++i) { success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i])); } @@ -157,8 +170,9 @@ namespace GridKit std::vector residual_y( residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); - auto& bus_residual_y_view = bus.getResidual(); - std::vector bus_residual_y( + auto& bus_residual_y_view = bus.getResidual(); + residual_y.insert( + residual_y.end(), bus_residual_y_view.getData(), bus_residual_y_view.getData() + bus_residual_y_view.getSize()); @@ -180,6 +194,11 @@ namespace GridKit std::vector residual_yp( residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); + auto& bus_residual_yp_view = bus.getResidual(); + residual_yp.insert( + residual_yp.end(), + bus_residual_yp_view.getData(), + bus_residual_yp_view.getData() + bus_residual_yp_view.getSize()); // Print the dependencies for (size_t i = 0; i < residual_y.size(); ++i) @@ -193,8 +212,7 @@ namespace GridKit } // Extract the dependencies and add d/dy' to d/dy - std::vector dependencies( - residual_y.size() + bus_residual_y.size()); + std::vector dependencies(residual_y.size()); for (IdxT i = 0; i < residual_y.size(); ++i) { DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); @@ -229,11 +247,6 @@ namespace GridKit } } - for (size_t i = 0; i < bus_residual_y.size(); ++i) - { - dependencies[residual_y.size() + i] = bus_residual_y[i].getDependencies(); - } - return dependencies; } diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index afc829795..04c32f433 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -219,6 +219,13 @@ namespace GridKit success *= isEqual(bus2.Ir(), Ir2); success *= isEqual(bus2.Ii(), Ii2); + system.evaluateResidual(); + + success *= isEqual(bus1.Ir(), Ir1); + success *= isEqual(bus1.Ii(), Ii1); + success *= isEqual(bus2.Ir(), Ir2); + success *= isEqual(bus2.Ii(), Ii2); + return success.report(__func__); } From 2e5f44b0fcc9d34bf470198cbb1f602263ea4bf1 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 31 Aug 2026 22:24:22 +0000 Subject: [PATCH 2/5] Apply pre-commit fixes --- .../Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp | 6 +++--- GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp | 8 ++++---- .../Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp | 6 +++--- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 19f03aada..afd2fc6bd 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -432,12 +432,12 @@ namespace GridKit auto read_signal = [&](const ScalarT& latched) { - const auto index = static_cast(variable); - y_ext[index] = latched; + const auto index = static_cast(variable); + y_ext[index] = latched; variable_indices_ext_[index] = INVALID_INDEX; if (auto port = ports_.in.template port()) { - y_ext[index] = port.readSignal(); + y_ext[index] = port.readSignal(); variable_indices_ext_[index] = port.signalVariableIndex(); } }; diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp index 7972715bf..0fa74542c 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp @@ -31,10 +31,10 @@ namespace GridKit { VREAL, ///< Real bus voltage VIMAG, ///< Imaginary bus voltage - VREF, ///< Voltage reference - VS, ///< Stabilizer output signal - VUEL, ///< Under-excitation limiter signal - VOEL, ///< Over-excitation limiter signal + VREF, ///< Voltage reference + VS, ///< Stabilizer output signal + VUEL, ///< Under-excitation limiter signal + VOEL, ///< Over-excitation limiter signal }; template diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index a726692c0..1498f5a37 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -294,12 +294,12 @@ namespace GridKit auto read_signal = [&](const ScalarT& latched) { - const auto index = static_cast(variable); - y_ext[index] = latched; + const auto index = static_cast(variable); + y_ext[index] = latched; variable_indices_ext_[index] = INVALID_INDEX; if (auto port = ports_.in.template port()) { - y_ext[index] = port.readSignal(); + y_ext[index] = port.readSignal(); variable_indices_ext_[index] = port.signalVariableIndex(); } }; From 8ca0f76bdcc01efac6c72d604ef1acf0a573df93 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 3 Sep 2026 16:03:32 +0000 Subject: [PATCH 3/5] Apply pre-commit fixes --- .../PhasorDynamics/Controller/REECB/ReecbImpl.hpp | 14 +++++++------- .../PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp | 12 ++++++------ .../PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp | 4 ++-- .../Exciter/SEXS-PTI/SexsPtiImpl.hpp | 4 ++-- .../Stabilizer/IEEEST/IeeestImpl.hpp | 2 +- .../SynchronousMachine/GENROU/GenrouImpl.hpp | 4 ++-- .../SynchronousMachine/GENSAL/GensalImpl.hpp | 4 ++-- 7 files changed, 22 insertions(+), 22 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp index 21604d7b2..9ff607981 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -799,8 +799,8 @@ namespace GridKit const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); const auto PREF = static_cast(ReecbExternalVariables::PREF); - y_ext[VR_EXT] = Vr(); - y_ext[VI_EXT] = Vi(); + y_ext[VR_EXT] = Vr(); + y_ext[VI_EXT] = Vi(); variable_indices_ext_[VR_EXT] = INVALID_INDEX; variable_indices_ext_[VI_EXT] = INVALID_INDEX; if (bus_->size() > 0) @@ -809,11 +809,11 @@ namespace GridKit variable_indices_ext_[VI_EXT] = bus_->getVariableIndex(1); } - y_ext[PE] = pe_set_; - y_ext[QGEN] = qgen_set_; - y_ext[QEXT] = qext_set_; - y_ext[PFAREF] = pfaref_set_; - y_ext[PREF] = pref_set_; + y_ext[PE] = pe_set_; + y_ext[QGEN] = qgen_set_; + y_ext[QEXT] = qext_set_; + y_ext[PFAREF] = pfaref_set_; + y_ext[PREF] = pref_set_; variable_indices_ext_[PE] = INVALID_INDEX; variable_indices_ext_[QGEN] = INVALID_INDEX; variable_indices_ext_[QEXT] = INVALID_INDEX; diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 15ecb5722..41105706d 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -448,8 +448,8 @@ namespace GridKit const auto VS = static_cast(Esdc1aExternalVariables::VS); const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); - y_ext[VR_EXT] = Vr(); - y_ext[VI_EXT] = Vi(); + y_ext[VR_EXT] = Vr(); + y_ext[VI_EXT] = Vi(); variable_indices_ext_[VR_EXT] = INVALID_INDEX; variable_indices_ext_[VI_EXT] = INVALID_INDEX; if (bus_->size() > 0) @@ -458,10 +458,10 @@ namespace GridKit variable_indices_ext_[VI_EXT] = bus_->getVariableIndex(1); } - y_ext[OMEGA] = omega_set_; - y_ext[VREF] = vref_set_; - y_ext[VS] = vs_set_; - y_ext[VUEL] = vuel_set_; + y_ext[OMEGA] = omega_set_; + y_ext[VREF] = vref_set_; + y_ext[VS] = vs_set_; + y_ext[VUEL] = vuel_set_; variable_indices_ext_[OMEGA] = INVALID_INDEX; variable_indices_ext_[VREF] = INVALID_INDEX; variable_indices_ext_[VS] = INVALID_INDEX; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index afd2fc6bd..77c4e2d8a 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -456,8 +456,8 @@ namespace GridKit // Bus voltages const auto VREAL = static_cast(Ieeet1ExternalVariables::VREAL); const auto VIMAG = static_cast(Ieeet1ExternalVariables::VIMAG); - y_ext[VREAL] = bus_->Vr(); - y_ext[VIMAG] = bus_->Vi(); + y_ext[VREAL] = bus_->Vr(); + y_ext[VIMAG] = bus_->Vi(); if (bus_->size() > 0) { variable_indices_ext_[VREAL] = bus_->getVariableIndex(0); diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index 1498f5a37..ed51739ed 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -315,8 +315,8 @@ namespace GridKit const auto VREAL = static_cast(SexsPtiExternalVariables::VREAL); const auto VIMAG = static_cast(SexsPtiExternalVariables::VIMAG); - y_ext[VREAL] = bus_->Vr(); - y_ext[VIMAG] = bus_->Vi(); + y_ext[VREAL] = bus_->Vr(); + y_ext[VIMAG] = bus_->Vi(); if (bus_->size() > 0) { variable_indices_ext_[VREAL] = bus_->getVariableIndex(0); diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index d560f2c8f..700bfbff3 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -353,7 +353,7 @@ namespace GridKit if (auto input_port = ports_.in.template port()) { - y_ext[0] = ports_.in.template port().readSignal(); + y_ext[0] = ports_.in.template port().readSignal(); variable_indices_ext_[0] = ports_.in.template port().signalVariableIndex(); } } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp index e58105bb1..ccbaad3a4 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp @@ -605,7 +605,7 @@ namespace GridKit y_ext[2] = pmech_set_; if (ports_.in.template port().connected()) { - y_ext[2] = ports_.in.template port().readSignal(); + y_ext[2] = ports_.in.template port().readSignal(); variable_indices_ext_[2] = ports_.in.template port().signalVariableIndex(); } @@ -613,7 +613,7 @@ namespace GridKit y_ext[3] = efd_set_; if (ports_.in.template port().connected()) { - y_ext[3] = ports_.in.template port().readSignal(); + y_ext[3] = ports_.in.template port().readSignal(); variable_indices_ext_[3] = ports_.in.template port().signalVariableIndex(); } } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp index 8268711fe..67f130835 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp @@ -458,7 +458,7 @@ namespace GridKit y_ext[2] = pmech_set_; if (ports_.in.template port().connected()) { - y_ext[2] = ports_.in.template port().readSignal(); + y_ext[2] = ports_.in.template port().readSignal(); variable_indices_ext_[2] = ports_.in.template port().signalVariableIndex(); } @@ -466,7 +466,7 @@ namespace GridKit y_ext[3] = efd_set_; if (ports_.in.template port().connected()) { - y_ext[3] = ports_.in.template port().readSignal(); + y_ext[3] = ports_.in.template port().readSignal(); variable_indices_ext_[3] = ports_.in.template port().signalVariableIndex(); } } From 04684ba6c18c26937d09776d1a8d265579fea5d6 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Wed, 9 Sep 2026 20:19:49 +0000 Subject: [PATCH 4/5] Apply pre-commit fixes --- .../SynchronousMachine/GenClassical/GenClassicalImpl.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index abf147b68..d136a7349 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -358,7 +358,7 @@ namespace GridKit y_ext[2] = pmech_set_; if (ports_.in.template port().connected()) { - y_ext[2] = ports_.in.template port().readSignal(); + y_ext[2] = ports_.in.template port().readSignal(); variable_indices_ext_[2] = ports_.in.template port().signalVariableIndex(); } @@ -366,7 +366,7 @@ namespace GridKit y_ext[3] = efd_set_; if (ports_.in.template port().connected()) { - y_ext[3] = ports_.in.template port().readSignal(); + y_ext[3] = ports_.in.template port().readSignal(); variable_indices_ext_[3] = ports_.in.template port().signalVariableIndex(); } } From 97cc8612fddea1bd6c4802af42cebf7c959846c9 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Fri, 11 Sep 2026 20:36:18 +0000 Subject: [PATCH 5/5] Apply pre-commit fixes --- GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp | 2 +- .../Model/PhasorDynamics/BusFault/BusFaultImpl.hpp | 2 +- .../PhasorDynamics/Controller/REECB/ReecbImpl.hpp | 1 - .../PhasorDynamics/Controller/REPCA/RepcaImpl.hpp | 1 - .../Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp | 12 ++++++------ .../PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp | 1 - .../Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp | 2 +- .../PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp | 2 +- 8 files changed, 10 insertions(+), 13 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp index 68d1000ce..09ff51094 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp @@ -11,11 +11,11 @@ #include -#include #include #include #include #include +#include namespace GridKit { diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 4c0df2de4..cba7c78b5 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -2,11 +2,11 @@ #include -#include #include #include #include #include +#include namespace GridKit { diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp index 9ff607981..12138328b 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -879,7 +879,6 @@ namespace GridKit return this->evaluateExternalResidual(); } - /** * @brief Access the optional variable monitor * diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp index 8f607a0e4..641a0a10a 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp @@ -705,7 +705,6 @@ namespace GridKit return this->evaluateExternalResidual(); } - /** * @brief Access the configured monitor * diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp index ed3e00142..251dae45a 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp @@ -46,12 +46,12 @@ namespace GridKit /// External variables read by an `Esdc1a`. enum class Esdc1aExternalVariables : size_t { - VR, ///< \f$V_\mathrm{r}\f$ Terminal-bus real voltage [p.u.] - VI, ///< \f$V_\mathrm{i}\f$ Terminal-bus imaginary voltage [p.u.] - OMEGA, ///< \f$\omega\f$ Known machine speed deviation [p.u.] - VREF, ///< \f$V_{\mathrm{ref}}\f$ Unknown voltage-control reference [p.u.] - VS, ///< \f$V_S\f$ Known stabilizer input signal [p.u.] - VUEL, ///< \f$V_{\mathrm{UEL}}\f$ Known under-excitation limiter input [p.u.] + VR, ///< \f$V_\mathrm{r}\f$ Terminal-bus real voltage [p.u.] + VI, ///< \f$V_\mathrm{i}\f$ Terminal-bus imaginary voltage [p.u.] + OMEGA, ///< \f$\omega\f$ Known machine speed deviation [p.u.] + VREF, ///< \f$V_{\mathrm{ref}}\f$ Unknown voltage-control reference [p.u.] + VS, ///< \f$V_S\f$ Known stabilizer input signal [p.u.] + VUEL, ///< \f$V_{\mathrm{UEL}}\f$ Known under-excitation limiter input [p.u.] }; /** diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp index 20b6dd9d0..ad31dd6fd 100644 --- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp @@ -497,7 +497,6 @@ namespace GridKit return this->evaluateExternalResidual(); } - /** * @brief Access the monitor * diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp index 17c976065..99378b810 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp @@ -3,11 +3,11 @@ #include #include -#include #include #include #include #include +#include namespace GridKit { diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp index b6f318d0a..cbd3bc750 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp @@ -1,10 +1,10 @@ #pragma once -#include #include #include #include #include +#include namespace GridKit {