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 7f4f4164c..c7a35a969 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -30,6 +30,16 @@ 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$ + MAXIMUM, + }; + /** * @brief Implementation of a line or off-nominal transformer branch between two buses. * @@ -49,8 +59,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_; @@ -83,7 +95,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; @@ -128,6 +142,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); @@ -142,11 +157,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(); @@ -188,13 +198,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..1827f0162 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 = f_ext_.size(); + auto y_ext_size = y_ext_.size(); + 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, + f_ext_.size(), + y_ext_.size(), + residual_indices_ext_.data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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..b04afaf4b 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchImpl.hpp @@ -127,8 +127,9 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(BranchExternalVariables::MAXIMUM)); + f_ext_.resize(4); + residual_indices_ext_.assign(4, INVALID_INDEX); allocated_ = true; return 0; @@ -197,20 +198,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 +217,83 @@ 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); - - 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}; + gatherExternalVariables(); - terminalCurrent1(ir1, ii1); - terminalCurrent2(ir2, ii2); + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); - 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 +307,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 2a4bf6365..bf1b6e720 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp @@ -20,6 +20,14 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `BusFault` + enum class BusFaultExternalVariables : size_t + { + VR, ///< \f$V_r\f$ + VI, ///< \f$V_i\f$ + MAXIMUM, + }; + template class BusFault : public Component { @@ -32,8 +40,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_; @@ -60,7 +70,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 @@ -93,6 +105,7 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; private: + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -116,7 +129,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 b96719c4f..ed07f3fb2 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 = f_ext_.size(); + auto y_ext_size = y_ext_.size(); + 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,47 +49,55 @@ namespace GridKit (this->getVariableIndices()).data(), y_.getData(), yp_.getData(), - wb_.data(), + y_ext_.data(), 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_); - if (!status_) // Value contributions from DfDwb only when status_ + const IdxT internal_external_begin = nnz_; + GridKit::Enzyme::Sparse::DfDyExt, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + static_cast(f_.getSize()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + if (!status_) // Value contributions from DfDyExt only when status_ { - for (IdxT i = nnz_tmp; i < nnz_; ++i) + for (IdxT i = internal_external_begin; i < nnz_; ++i) { J_vals_buffer_[i] = 0.0; } } - 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_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + const IdxT external_residual_begin = nnz_; + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + f_ext_.size(), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + if (!status_) // External bus-row contributions only when status_ + { + for (IdxT i = external_residual_begin; i < nnz_; ++i) + { + J_vals_buffer_[i] = 0.0; + } + } this->constructCoo(); diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 2a2953360..85428ab5a 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -125,8 +125,9 @@ namespace GridKit residual_indices_.resize(size); // Resize coupling data - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(BusFaultExternalVariables::MAXIMUM)); + f_ext_.resize(2); + 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 +206,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 +232,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,44 +246,84 @@ 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() + { + 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_.data(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * \brief External residual contributions to the bus while the fault is + * active. + * + */ + template + int BusFault::evaluateExternalResidual() { if (status_) { - wb_[0] = Vr(); - wb_[1] = Vi(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); - auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb_.data(), f); - evaluateBusResidual(y, yp, wb_.data(), h_.data()); - Ir() += h_[0]; - Ii() += h_[1]; + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); + 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_.data(), 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 a930c8ffb..61c1a8794 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_; @@ -299,6 +311,16 @@ namespace GridKit abs_tol_.resize(n); } + /** + * @brief Allocate this component's external variable vectors. + */ + void allocateExternalVectors(IdxT n) + { + y_ext_.assign(static_cast(n), ScalarT{}); + yp_ext_.assign(static_cast(n), ScalarT{}); + variable_indices_ext_.assign(static_cast(n), INVALID_INDEX); + } + int constructCoo() { if (coo_jac_ == nullptr) @@ -327,15 +349,20 @@ 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_; - - VectorT y_; - VectorT yp_; - std::vector tag_; - VectorT abs_tol_; - VectorT f_; - bool allocated_{false}; + std::vector residual_indices_ext_; + + VectorT y_; + std::vector y_ext_; + VectorT yp_; + std::vector yp_ext_; + std::vector tag_; + VectorT abs_tol_; + VectorT f_; + std::vector f_ext_; + bool allocated_{false}; std::vector g_; @@ -361,9 +388,6 @@ namespace GridKit IdxT gridkit_component_id_{0}; - std::vector wb_; - std::vector h_; - RealT time_; RealT alpha_; diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp index 73213d15d..94be0c888 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp @@ -8,7 +8,6 @@ #include #include -#include #include #include @@ -58,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.] @@ -95,8 +96,9 @@ namespace GridKit using Component::tag_; using Component::va_system_base_; using Component::variable_indices_; - using Component::wb_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -122,6 +124,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; @@ -136,8 +139,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: @@ -146,6 +148,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; @@ -213,9 +216,6 @@ namespace GridKit ComponentSignals signals_; std::unique_ptr monitor_; - - std::vector ws_; - std::vector ws_indices_; }; } // namespace Controller } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp index 0f6b1451e..3a9cd57fb 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 = ws_.size(); - const auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + 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 +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_.data(), - ws_.data(), + y_ext_.data(), 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_.data(), - ws_.data(), + y_ext_.data(), 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_.data(), - ws_.data(), - 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()), - ws_.size(), + y_ext_.size(), this->getResidualIndices().data(), - ws_indices_.data(), + variable_indices_ext_.data(), y_.getData(), yp_.getData(), - wb_.data(), - ws_.data(), + y_ext_.data(), 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 74a939ee1..d334ff1e3 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp @@ -81,7 +81,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. @@ -102,11 +102,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.assign(2, ScalarT{0}); - - const auto signal_size = static_cast(RepcaExternalVariables::MAXIMUM); - ws_.assign(signal_size, ScalarT{0}); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(RepcaExternalVariables::MAXIMUM)); for (IdxT j = 0; j < size_; ++j) { @@ -588,15 +584,17 @@ 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() { + 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); @@ -607,78 +605,108 @@ namespace GridKit const auto QREF = static_cast(RepcaExternalVariables::QREF); const auto FREQREF = static_cast(RepcaExternalVariables::FREQREF); - 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); + 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_; - ws_[IR] = + y_ext_[IR] = signals_.template readExternalVariable(); - ws_indices_[IR] = + variable_indices_ext_[IR] = signals_.template readExternalVariableIndex(); - ws_[II] = + y_ext_[II] = signals_.template readExternalVariable(); - ws_indices_[II] = + variable_indices_ext_[II] = signals_.template readExternalVariableIndex(); - ws_[P] = + y_ext_[P] = signals_.template readExternalVariable(); - ws_indices_[P] = + variable_indices_ext_[P] = signals_.template readExternalVariableIndex(); - ws_[Q] = + y_ext_[Q] = signals_.template readExternalVariable(); - ws_indices_[Q] = + variable_indices_ext_[Q] = signals_.template readExternalVariableIndex(); if (signals_.template isAttached()) { - ws_[FREQ] = + y_ext_[FREQ] = signals_.template readExternalVariable(); - ws_indices_[FREQ] = + variable_indices_ext_[FREQ] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[VREF] = + y_ext_[VREF] = signals_.template readExternalVariable(); - ws_indices_[VREF] = + variable_indices_ext_[VREF] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[PREF] = + y_ext_[PREF] = signals_.template readExternalVariable(); - ws_indices_[PREF] = + variable_indices_ext_[PREF] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[QREF] = + y_ext_[QREF] = signals_.template readExternalVariable(); - ws_indices_[QREF] = + variable_indices_ext_[QREF] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[FREQREF] = + y_ext_[FREQREF] = signals_.template readExternalVariable(); - ws_indices_[FREQREF] = + variable_indices_ext_[FREQREF] = signals_.template readExternalVariableIndex(); } + } - 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_.data(), ws_.data(), f); + evaluateInternalResidual(y, yp, y_ext_.data(), 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 REPCA signal interface * @@ -716,8 +744,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. */ @@ -726,8 +754,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); @@ -753,6 +780,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); @@ -794,18 +823,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 = toComponentBase(ws[IR]); - const ScalarT ii = toComponentBase(ws[II]); - const ScalarT p = toComponentBase(ws[P]); - const ScalarT q = toComponentBase(ws[Q]); - const ScalarT freq = ws[FREQ]; - const ScalarT freqref = ws[FREQREF]; - const ScalarT vref = ws[VREF]; - const ScalarT qref = toComponentBase(ws[QREF]); - const ScalarT pref_in = toComponentBase(ws[PREF_INPUT]); + const ScalarT vr = y_ext[VR]; + const ScalarT vi = y_ext[VI]; + + const ScalarT ir = toComponentBase(y_ext[IR]); + const ScalarT ii = toComponentBase(y_ext[II]); + const ScalarT p = toComponentBase(y_ext[P]); + const ScalarT q = 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 = toComponentBase(y_ext[QREF]); + const ScalarT pref_in = 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 2bfb937ca..442c5aab4 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 MAXIMUM, @@ -71,18 +73,20 @@ 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::va_system_base_; using Component::variable_indices_; - using Component::wb_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -107,6 +111,8 @@ namespace GridKit int tagDifferentiable() override final; int setAbsoluteTolerance(RealT rel_tol) override final; void setLvplGain(RealT KL); + int evaluateInternalResidual() override final; + int evaluateExternalResidual() override final; int evaluateResidual() override final; int evaluateJacobian() override final; @@ -122,15 +128,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(); ScalarT toComponentBase(ScalarT value) const; ScalarT toSystemBase(ScalarT value) const; @@ -239,10 +246,6 @@ namespace GridKit ComponentSignals signals_; std::unique_ptr monitor_; - - // Local copies of signal variables - std::vector ws_; - std::vector ws_indices_; }; } // namespace Converter } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp index 7627549eb..908ae0b22 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_.size()); - auto buffer_size = 2 * size * size + 2 * size * bus_size + size * signal_size; + auto y_ext_size = y_ext_.size(); + auto f_ext_size = f_ext_.size(); + 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_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), - ws_.data(), - 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_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws::eval(this, - static_cast(f_.getSize()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt::eval(this, + static_cast(f_.getSize()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy::eval(this, + f_ext_.size(), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 d5a8ae03c..b8cc94394 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -359,12 +359,10 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.assign(2, ScalarT{0}); - h_.assign(2, ScalarT{0}); - - auto signal_size = static_cast(RegcaExternalVariables::MAXIMUM); - ws_.assign(signal_size, ScalarT{0}); - ws_indices_.assign(signal_size, INVALID_INDEX); + // Resize coupling data + this->allocateExternalVectors(static_cast(RegcaExternalVariables::MAXIMUM)); + f_ext_.assign(2, ScalarT{0}); + residual_indices_ext_.assign(2, INVALID_INDEX); for (IdxT j = 0; j < size_; ++j) { @@ -606,8 +604,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 @@ -615,8 +613,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); @@ -630,6 +627,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); @@ -648,11 +647,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 = toComponentBase(ws[IPCMD]); - const ScalarT iqcmd = toComponentBase(ws[IQCMD]); + const ScalarT ipcmd = toComponentBase(y_ext[IPCMD]); + const ScalarT iqcmd = toComponentBase(y_ext[IQCMD]); // Form the unconstrained current derivatives, then apply the REGCA // recovery rate limits in p.u./s. @@ -697,74 +696,119 @@ 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() { + 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); - 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); + } + // Active-current command + y_ext_[IPCMD] = ipcmd_set_; if (signals_.template isAttached()) { - ws_[IPCMD] = signals_.template readExternalVariable(); - ws_indices_[IPCMD] = + y_ext_[IPCMD] = + signals_.template readExternalVariable(); + variable_indices_ext_[IPCMD] = signals_.template readExternalVariableIndex(); } + // Reactive-current command + y_ext_[IQCMD] = iqcmd_set_; if (signals_.template isAttached()) { - ws_[IQCMD] = signals_.template readExternalVariable(); - ws_indices_[IQCMD] = + y_ext_[IQCMD] = + signals_.template readExternalVariable(); + variable_indices_ext_[IQCMD] = signals_.template readExternalVariableIndex(); } + } - 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(); - evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); - evaluateBusResidual(y, yp, wb_.data(), h_.data()); + evaluateInternalResidual(y, yp, y_ext_.data(), 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() + { + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); + + // 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 dd42f55ac..8967c70d1 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp @@ -8,7 +8,6 @@ #include #include -#include #include #include @@ -44,14 +43,16 @@ namespace GridKit MAXIMUM, ///< Number of ESDC1A internal variables }; - /// External signal variables read or initialized by an `Esdc1a`. + /// 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.] - MAXIMUM, ///< Number of ESDC1A external signal variables + MAXIMUM, ///< Number of ESDC1A external variables }; /** @@ -76,8 +77,9 @@ namespace GridKit using Component::size_; using Component::tag_; using Component::variable_indices_; - using Component::wb_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -101,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; @@ -118,14 +121,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, @@ -177,9 +180,6 @@ namespace GridKit ComponentSignals signals_; std::unique_ptr monitor_; - - std::vector ws_; - std::vector ws_indices_; }; } // namespace Exciter } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp index 812ea4031..519ba8f0f 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 = ws_.size(); - auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + auto y_ext_size = y_ext_.size(); + 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_.data(), - ws_.data(), - 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_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws::eval(this, - static_cast(f_.getSize()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt::eval(this, + static_cast(f_.getSize()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 11906d3d9..a5acb7fad 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -106,11 +106,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.assign(2, ScalarT{0}); - - const auto signal_size = static_cast(Esdc1aExternalVariables::MAXIMUM); - ws_.assign(signal_size, ScalarT{0}); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(Esdc1aExternalVariables::MAXIMUM)); for (IdxT j = 0; j < size_; ++j) { @@ -420,66 +416,98 @@ 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); + 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 (signals_.template isAttached()) { - ws_[OMEGA] = signals_.template readExternalVariable(); - ws_indices_[OMEGA] = + y_ext_[OMEGA] = + signals_.template readExternalVariable(); + variable_indices_ext_[OMEGA] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[VREF] = signals_.template readExternalVariable(); - ws_indices_[VREF] = + y_ext_[VREF] = + signals_.template readExternalVariable(); + variable_indices_ext_[VREF] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[VS] = signals_.template readExternalVariable(); - ws_indices_[VS] = + y_ext_[VS] = signals_.template readExternalVariable(); + variable_indices_ext_[VS] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[VUEL] = signals_.template readExternalVariable(); - ws_indices_[VUEL] = + y_ext_[VUEL] = + signals_.template readExternalVariable(); + variable_indices_ext_[VUEL] = signals_.template readExternalVariableIndex(); } + } - 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_.data(), ws_.data(), f); + evaluateInternalResidual(y, yp, y_ext_.data(), 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 * @@ -503,9 +531,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. */ @@ -514,8 +541,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); @@ -530,10 +556,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]; @@ -553,12 +581,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_ + se) * efdp; 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 36d03a16f..25ca2a943 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -77,7 +77,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; + using Component::y_ext_; + using Component::variable_indices_ext_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; @@ -105,6 +106,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; @@ -121,7 +123,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); @@ -168,9 +170,7 @@ namespace GridKit /// Associate variable getter functions with enum values void initializeMonitor(); - /* Local copies of signal variables */ - std::vector ws_; - std::vector ws_indices_; + void gatherExternalVariables(); }; } // namespace Exciter diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp index 77c59abf4..7f1af558f 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_.size()); - auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + auto y_ext_size = y_ext_.size(); + 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_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), - 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_.data(), + 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()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 58d016ae6..985fc26a7 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -94,16 +94,8 @@ namespace GridKit this->setResidualIndex(j, j); } - // Resize bus data - wb_.resize(2); - - // Resize signal variable data - ws_.resize(2); - ws_indices_.resize(2); - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; - ws_[1] = 0.0; - ws_indices_[1] = INVALID_INDEX; + // Resize coupling data + this->allocateExternalVectors(static_cast(Ieeet1ExternalVariables::MAXIMUM)); // Set output signals if (signals_.template isAssigned()) @@ -303,13 +295,12 @@ 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) { // Read bus voltage components - ScalarT vreal = wb[0]; - ScalarT vimag = wb[1]; + ScalarT vreal = y_ext[1]; + ScalarT vimag = y_ext[2]; ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); // Read Internal Variables @@ -330,8 +321,8 @@ namespace GridKit ScalarT vfx_dot = yp[3]; // Set signal variable aliases - ScalarT omega = ws[0]; - ScalarT vs_signal = ws[1]; + ScalarT omega = y_ext[0]; + ScalarT vs_signal = y_ext[3]; // The 'pre-limit' derivative of Vr. ScalarT func = (-vr + Ka_ * vtr) / Ta_; @@ -353,43 +344,65 @@ namespace GridKit } /** - * @brief Residual evaluation + * @brief Gather external variables and index maps. * */ template - int Ieeet1::evaluateResidual() + void Ieeet1::gatherExternalVariables() { // Set input variables. if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + y_ext_[0] = signals_.template readExternalVariable(); + variable_indices_ext_[0] = signals_.template readExternalVariableIndex(); + } + + // Bus voltages + y_ext_[1] = bus_->Vr(); + y_ext_[2] = bus_->Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[1] = bus_->getVariableIndex(0); + variable_indices_ext_[2] = bus_->getVariableIndex(1); } // VS signal (stabilizer output, optional) - ws_[1] = 0.0; - ws_indices_[1] = INVALID_INDEX; + y_ext_[3] = 0.0; + variable_indices_ext_[3] = INVALID_INDEX; if (signals_.template isAttached()) { - ws_[1] = signals_.template readExternalVariable(); - ws_indices_[1] = signals_.template readExternalVariableIndex(); + y_ext_[3] = signals_.template readExternalVariable(); + variable_indices_ext_[3] = signals_.template readExternalVariableIndex(); } + } - // Bus voltages - wb_[0] = bus_->Vr(); - wb_[1] = bus_->Vi(); + /** + * @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_.data(), ws_.data(), f); + evaluateInternalResidual(y, yp, y_ext_.data(), 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 b1268515e..e94eb3aef 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp @@ -47,6 +47,8 @@ namespace GridKit /// External variables of a `SexsPti`. enum class SexsPtiExternalVariables : size_t { + VR, ///< Real bus voltage + VI, ///< Imaginary bus voltage VS, ///< Stabilizer output signal MAXIMUM, }; @@ -64,7 +66,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; + 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 +94,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; @@ -106,7 +110,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}; @@ -130,9 +134,7 @@ namespace GridKit void initModelParams(const ModelDataT& data); void initializeMonitor(); - - std::vector ws_; - std::vector ws_indices_; + 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 d7debf0fe..885b7f2a5 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_.size()); - auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + auto y_ext_size = y_ext_.size(); + 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_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), - 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_.data(), + 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()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 8b8e53a11..77a91621d 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); - - ws_.resize(1); - ws_indices_.resize(1); - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; + this->allocateExternalVectors(static_cast(SexsPtiExternalVariables::MAXIMUM)); if (signals_.template isAssigned()) { @@ -214,8 +209,7 @@ 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) { ScalarT vr = y[0]; @@ -224,8 +218,8 @@ 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 vs = ws[0]; + ScalarT Ec = std::sqrt(y_ext[0] * y_ext[0] + y_ext[1] * y_ext[1]); + ScalarT vs = y_ext[2]; ScalarT func = (-efd + (K_ / Tb_) * (-vr + Ta_ * vtr)) / Te_; @@ -236,30 +230,52 @@ namespace GridKit return 0; } + /** + * @brief Gather external variables and index maps. + * + */ template - int SexsPti::evaluateResidual() + void SexsPti::gatherExternalVariables() { - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; + y_ext_[0] = bus_->Vr(); + y_ext_[1] = bus_->Vi(); + if (bus_->size() > 0) + { + variable_indices_ext_[0] = bus_->getVariableIndex(0); + variable_indices_ext_[1] = bus_->getVariableIndex(1); + } + + y_ext_[2] = 0.0; + variable_indices_ext_[2] = INVALID_INDEX; if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + y_ext_[2] = signals_.template readExternalVariable(); + variable_indices_ext_[2] = signals_.template readExternalVariableIndex(); } + } - wb_[0] = bus_->Vr(); - wb_[1] = bus_->Vi(); + template + int SexsPti::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); + evaluateInternalResidual(y, yp, y_ext_.data(), 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/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp index 9844dffd3..81897ef0c 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -10,7 +10,6 @@ #include #include #include -#include #include #include @@ -77,8 +76,9 @@ namespace GridKit using Component::tag_; using Component::va_system_base_; using Component::variable_indices_; - using Component::wb_; + using Component::variable_indices_ext_; using Component::y_; + using Component::y_ext_; using Component::yp_; public: @@ -101,6 +101,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; @@ -118,11 +119,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(); @@ -190,9 +191,6 @@ namespace GridKit ComponentSignals signals_; std::unique_ptr monitor_; - - std::vector ws_; - std::vector ws_indices_; }; } // namespace Governor } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp index 174f23959..c2f701f56 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_.size()); - auto buffer_size = 2 * size * size + size * signal_size; + 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]; @@ -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(), - wb_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws::eval(this, - static_cast(f_.getSize()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyExt::eval(this, + static_cast(f_.getSize()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 9761ca9c6..706b1b58a 100644 --- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -73,11 +73,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. @@ -97,11 +96,7 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - wb_.clear(); - - const auto signal_size = static_cast(HygovExternalVariables::MAXIMUM); - ws_.assign(signal_size, ScalarT{0}); - ws_indices_.assign(signal_size, INVALID_INDEX); + this->allocateExternalVectors(static_cast(HygovExternalVariables::MAXIMUM)); for (IdxT j = 0; j < size_; ++j) { @@ -487,56 +482,75 @@ namespace GridKit } /** - * @brief Residuals of system equations - * - * 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. + * @brief Gather external signal values and global indices. * - * @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() { const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); const auto PREF = static_cast(HygovExternalVariables::PREF); const auto PAUX = static_cast(HygovExternalVariables::PAUX); - 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 (signals_.template isAttached()) { - ws_[OMEGA] = signals_.template readExternalVariable(); - ws_indices_[OMEGA] = + y_ext_[OMEGA] = signals_.template readExternalVariable(); + variable_indices_ext_[OMEGA] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[PREF] = signals_.template readExternalVariable(); - ws_indices_[PREF] = + y_ext_[PREF] = signals_.template readExternalVariable(); + variable_indices_ext_[PREF] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - ws_[PAUX] = signals_.template readExternalVariable(); - ws_indices_[PAUX] = + y_ext_[PAUX] = signals_.template readExternalVariable(); + variable_indices_ext_[PAUX] = signals_.template readExternalVariableIndex(); } + } + + /** + * @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, wb_.data(), ws_.data(), f); + evaluateInternalResidual(y, yp, y_ext_.data(), 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 * @@ -560,19 +574,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); @@ -610,9 +622,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 bb8abc988..1b6b8d17c 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -68,8 +68,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - 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_; @@ -96,6 +96,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 @@ -113,7 +114,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 @@ -137,13 +138,10 @@ namespace GridKit // Parameter initialization function void initializeParameters(const ModelDataT& data); + void gatherExternalVariables(); void setDerivedParams(); ScalarT toComponentBase(ScalarT value) const; ScalarT toSystemBase(ScalarT value) const; - - /* Local copies of signal variables */ - std::vector ws_; - std::vector ws_indices_; }; } // namespace Governor diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp index b247628a1..2af52099e 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_.size()); - auto buffer_size = 2 * size * size + size * signal_size; + auto y_ext_size = y_ext_.size(); + 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(), - wb_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), + 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()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 ee5bed430..c59d9fbfa 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -182,11 +182,8 @@ namespace GridKit this->setResidualIndex(j, j); } - // Resize signal variable data - ws_.resize(1); - ws_indices_.resize(1); - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; + // Resize coupling data + this->allocateExternalVectors(static_cast(Tgov1ExternalVariables::MAXIMUM)); // Set output signals if (signals_.template isAssigned()) @@ -296,11 +293,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) { // Read Internal Variables ScalarT ptx = y[0]; // y0 - Ptx @@ -312,7 +308,7 @@ namespace GridKit ScalarT pv_dot = yp[1]; // Set signal variable aliases - ScalarT omega = ws[0]; + ScalarT omega = y_ext[0]; // The 'pre-limit' target of Pv ScalarT func = -pv + (pref_ - omega) / R_; @@ -332,25 +328,42 @@ namespace GridKit * @brief Residuals of system equations * */ + /** + * @brief Gather external variables and index maps. + * + */ template - int Tgov1::evaluateResidual() + void Tgov1::gatherExternalVariables() { // Input Variables if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + y_ext_[0] = signals_.template readExternalVariable(); + variable_indices_ext_[0] = signals_.template readExternalVariableIndex(); } + } + + template + int Tgov1::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); + evaluateInternalResidual(y, yp, y_ext_.data(), 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 e31259e19..aa0793394 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp @@ -22,6 +22,14 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `LoadZ` + enum class LoadZExternalVariables : size_t + { + VR, ///< \f$V_r\f$ + VI, ///< \f$V_i\f$ + MAXIMUM, + }; + /*! * @brief Implementation of a constant load. * @@ -38,8 +46,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_; @@ -66,7 +76,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 @@ -89,6 +101,7 @@ namespace GridKit private: void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -114,7 +127,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 4e5304cf4..83c174794 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 = f_ext_.size(); + auto y_ext_size = y_ext_.size(); + 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_.data(), + y_ext_.data(), 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + f_ext_.size(), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 f3708e491..00e68942b 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp @@ -102,8 +102,9 @@ namespace GridKit } // Resize coupling data - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(LoadZExternalVariables::MAXIMUM)); + f_ext_.resize(2); + residual_indices_ext_.assign(2, INVALID_INDEX); allocated_ = true; return 0; @@ -167,20 +168,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 +194,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,30 +208,68 @@ 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() + { + 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() { - wb_[0] = Vr(); - wb_[1] = Vi(); + gatherExternalVariables(); + const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb_.data(), f); - evaluateBusResidual(y, yp, wb_.data(), h_.data()); - Ir() += h_[0]; - Ii() += h_[1]; + evaluateInternalResidual(y, yp, y_ext_.data(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * @brief External residual contributions to the bus. + * + */ + template + int LoadZ::evaluateExternalResidual() + { + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); + 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 5387e6f27..dc7a06341 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp @@ -22,6 +22,14 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `LoadZIP` + enum class LoadZIPExternalVariables : size_t + { + VR, ///< \f$V_r\f$ + VI, ///< \f$V_i\f$ + MAXIMUM, + }; + /*! * @brief Implementation of a ZIP load. * @@ -38,8 +46,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_; @@ -66,7 +76,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 @@ -102,6 +114,7 @@ namespace GridKit private: void initializeParameters(const ModelDataT& data); void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); ScalarT& Vr() @@ -127,7 +140,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 b067c0cbd..c410c418f 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 = f_ext_.size(); + auto y_ext_size = y_ext_.size(); + 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_.data(), + y_ext_.data(), 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + f_ext_.size(), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 58058b07c..fffa0f8e2 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp @@ -115,8 +115,9 @@ namespace GridKit } // Resize coupling data - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(LoadZIPExternalVariables::MAXIMUM)); + f_ext_.resize(2); + residual_indices_ext_.assign(2, INVALID_INDEX); allocated_ = true; return 0; @@ -188,49 +189,91 @@ 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() + { + 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() { - wb_[0] = Vr(); - wb_[1] = Vi(); + gatherExternalVariables(); + const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb_.data(), f); - evaluateBusResidual(y, yp, wb_.data(), h_.data()); - Ir() += h_[0]; - Ii() += h_[1]; + evaluateInternalResidual(y, yp, y_ext_.data(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * @brief External residual contributions to the bus. + * + */ + template + int LoadZIP::evaluateExternalResidual() + { + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); + 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 * @@ -239,11 +282,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 f77935130..4718dd5c1 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -70,8 +70,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - using Component::wb_; - 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_; @@ -97,6 +97,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,7 +117,6 @@ namespace GridKit const ScalarT*, const ScalarT*, const ScalarT*, - const ScalarT*, ScalarT*); private: @@ -167,9 +167,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); void initializeMonitor(); - - std::vector ws_; - std::vector ws_indices_; + void gatherExternalVariables(); }; } // namespace Stabilizer diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp index 4c1af73d9..ddd10a275 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_.size()); - auto buffer_size = 2 * size * size + size * signal_size; + auto y_ext_size = y_ext_.size(); + 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(), - wb_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), + 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()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 970907ab0..1ed58af6e 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -172,10 +172,7 @@ namespace GridKit this->setResidualIndex(j, j); } - ws_.resize(1); - ws_indices_.resize(1); - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; + this->allocateExternalVectors(static_cast(IeeestExternalVariables::MAXIMUM)); if (signals_.template isAssigned()) { @@ -274,11 +271,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]; @@ -301,7 +297,7 @@ namespace GridKit ScalarT x6_dot = yp[5]; ScalarT x7_dot = yp[6]; - ScalarT u = ws[0]; + ScalarT u = y_ext[0]; f[0] = -x1_dot + use_notch_ * x2; f[1] = -x2_dot + (use_4th_order_ + use_3rd_order_) * x3 @@ -321,24 +317,41 @@ namespace GridKit return 0; } + /** + * @brief Gather external variables and index maps. + * + */ template - int Ieeest::evaluateResidual() + void Ieeest::gatherExternalVariables() { if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + y_ext_[0] = signals_.template readExternalVariable(); + variable_indices_ext_[0] = signals_.template readExternalVariableIndex(); } + } + + template + int Ieeest::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); + evaluateInternalResidual(y, yp, y_ext_.data(), 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 75e765b24..e43400bb2 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp @@ -79,8 +79,10 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - 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_; @@ -137,7 +139,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; @@ -158,6 +162,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); /// Associate variable getter functions with enum values void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); /** @@ -202,8 +207,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: @@ -259,10 +264,6 @@ namespace GridKit ScalarT pmech_set_{0.0}; // TODO remove default initialization and ensure this gets set ScalarT efd_set_{0.0}; // TODO remove default initialization and ensure this gets set - /* Local copies of signal variables */ - std::vector ws_; - std::vector ws_indices_; - /// Variable monitor std::unique_ptr monitor_; }; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp index 834d027ab..501695fb8 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_.size()); - auto buffer_size = 2 * size * size + size * signal_size + 2 * size * bus_size; + auto f_ext_size = f_ext_.size(); + auto y_ext_size = y_ext_.size(); + 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_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), - ws_.data(), - 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_.data(), + 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()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + f_ext_.size(), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 2d933043d..de2e38b4e 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp @@ -337,15 +337,10 @@ 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(GenrouExternalVariables::MAXIMUM)); + f_ext_.resize(2); + residual_indices_ext_.assign(2, INVALID_INDEX); // Set output signals if (signals_.template isAssigned()) @@ -554,8 +549,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 */ @@ -588,12 +582,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 = toMachineBase(ws[0]); - ScalarT efd = ws[1]; + ScalarT pmech = toMachineBase(y_ext[2]); + ScalarT efd = y_ext[3]; static constexpr auto pi = std::numbers::pi_v; @@ -627,73 +621,113 @@ 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] = toSystemBase(ir); - h[1] = toSystemBase(ii); + f_ext[0] = toSystemBase(ir); + f_ext[1] = 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() { + // 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_; + y_ext_[2] = pmech_set_; if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + y_ext_[2] = signals_.template readExternalVariable(); + variable_indices_ext_[2] = signals_.template readExternalVariableIndex(); } // Exciter Efield - ws_[1] = efd_set_; + y_ext_[3] = efd_set_; if (signals_.template isAttached()) { - ws_[1] = signals_.template readExternalVariable(); - ws_indices_[1] = signals_.template readExternalVariableIndex(); + y_ext_[3] = signals_.template readExternalVariable(); + variable_indices_ext_[3] = signals_.template readExternalVariableIndex(); } + } - // Bus voltages - 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(); - evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); - evaluateBusResidual(y, yp, wb_.data(), h_.data()); + evaluateInternalResidual(y, yp, y_ext_.data(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * \brief External residual contributions to the bus. + * + */ + template + int Genrou::evaluateExternalResidual() + { + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); // 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 2060834ef..389de653f 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp @@ -74,8 +74,10 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - 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_; @@ -102,7 +104,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; @@ -123,6 +127,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); /// Associate variable getter functions with enum values void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); /** @@ -167,8 +172,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: @@ -215,10 +220,6 @@ namespace GridKit ScalarT pmech_set_{0.0}; // TODO remove default initialization and ensure this gets set ScalarT efd_set_{0.0}; // TODO remove default initialization and ensure this gets set - /* Local copies of signal variables */ - std::vector ws_; - std::vector ws_indices_; - /// Variable monitor std::unique_ptr monitor_; }; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp index 185eb20d7..fd99163b4 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_.size()); - auto buffer_size = 2 * size * size + size * signal_size + 2 * size * bus_size; + auto f_ext_size = f_ext_.size(); + auto y_ext_size = y_ext_.size(); + 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_.data(), - ws_.data(), - 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_.data(), + 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_.data(), - ws_.data(), - 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_.data(), - ws_.data(), - 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_.data(), + 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()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + f_ext_.size(), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 ebb3f1471..8d0c5e5ef 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp @@ -212,15 +212,10 @@ 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(GensalExternalVariables::MAXIMUM)); + f_ext_.resize(2); + residual_indices_ext_.assign(2, INVALID_INDEX); // Set output signals if (signals_.template isAssigned()) @@ -384,8 +379,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 */ @@ -414,12 +408,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 = toMachineBase(ws[0]); - ScalarT efd = ws[1]; + ScalarT pmech = toMachineBase(y_ext[2]); + ScalarT efd = y_ext[3]; static constexpr auto pi = std::numbers::pi_v; @@ -450,73 +444,113 @@ 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] = toSystemBase(ir); - h[1] = toSystemBase(ii); + f_ext[0] = toSystemBase(ir); + f_ext[1] = 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() { + // 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_; + y_ext_[2] = pmech_set_; if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + y_ext_[2] = signals_.template readExternalVariable(); + variable_indices_ext_[2] = signals_.template readExternalVariableIndex(); } // Exciter Efield - ws_[1] = efd_set_; + y_ext_[3] = efd_set_; if (signals_.template isAttached()) { - ws_[1] = signals_.template readExternalVariable(); - ws_indices_[1] = signals_.template readExternalVariableIndex(); + y_ext_[3] = signals_.template readExternalVariable(); + variable_indices_ext_[3] = signals_.template readExternalVariableIndex(); } + } - // Bus voltages - 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(); - evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); - evaluateBusResidual(y, yp, wb_.data(), h_.data()); + evaluateInternalResidual(y, yp, y_ext_.data(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * \brief External residual contributions to the bus. + * + */ + template + int Gensal::evaluateExternalResidual() + { + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); // 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 14ba8ce36..0d6a9eeb6 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -29,6 +29,13 @@ namespace GridKit { namespace PhasorDynamics { + /// External variables of a `GenClassical` + enum class GenClassicalExternalVariables : size_t + { + VR, ///< \f$V_r\f$ + VI, ///< \f$V_i\f$ + MAXIMUM, + }; template class GenClassical : public Component @@ -43,8 +50,10 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; - 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_; @@ -78,7 +87,9 @@ namespace GridKit int initialize() override final; int tagDifferentiable() override final; int setAbsoluteTolerance(RealT) override; + int evaluateInternalResidual() override final; int evaluateResidual() override final; + int evaluateExternalResidual() override final; int verify() const override final { @@ -102,6 +113,7 @@ namespace GridKit private: void initializeMonitor(); + void gatherExternalVariables(); void setDerivedParams(); /** @@ -147,7 +159,7 @@ namespace GridKit public: __attribute__((always_inline)) inline int evaluateInternalResidual( const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); - __attribute__((always_inline)) inline int evaluateBusResidual( + __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 0e9d3b8fe..cd59f5347 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp @@ -23,14 +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 buffer_size = 2 * size * size + 2 * size * bus_size; + auto f_ext_size = f_ext_.size(); + auto y_ext_size = y_ext_.size(); + 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]; @@ -46,7 +49,7 @@ namespace GridKit (this->getVariableIndices()).data(), y_.getData(), yp_.getData(), - wb_.data(), + y_ext_.data(), J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, @@ -60,40 +63,40 @@ namespace GridKit (this->getVariableIndices()).data(), y_.getData(), yp_.getData(), - wb_.data(), + y_ext_.data(), alpha_, 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()), + y_ext_.size(), + (this->getResidualIndices()).data(), + variable_indices_ext_.data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::ExternalResidual>::eval(this, + f_ext_.size(), + static_cast(y_.getSize()), + residual_indices_ext_.data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + y_ext_.data(), + 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 5bd9ab0eb..d8d957f85 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -196,8 +196,9 @@ namespace GridKit } // Resize coupling data - wb_.resize(2); - h_.resize(2); + this->allocateExternalVectors(static_cast(GenClassicalExternalVariables::MAXIMUM)); + f_ext_.resize(2); + residual_indices_ext_.assign(2, INVALID_INDEX); allocated_ = true; return 0; @@ -283,7 +284,7 @@ namespace GridKit __attribute__((always_inline)) int GenClassical::evaluateInternalResidual( const ScalarT* y, const ScalarT* yp, - const ScalarT* wb, + const ScalarT* y_ext, ScalarT* f) { // Set variable aliases for better readability. @@ -300,8 +301,8 @@ 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]; static constexpr auto pi = std::numbers::pi_v; @@ -319,52 +320,93 @@ namespace GridKit } /** - * @brief Bus residual + * @brief External residual * */ template - __attribute__((always_inline)) int GenClassical::evaluateBusResidual( + __attribute__((always_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] = toSystemBase(ir); - h[1] = toSystemBase(ii); + f_ext[0] = toSystemBase(ir); + f_ext[1] = toSystemBase(ii); return 0; } /** - * \brief Residual for the generator model. + * @brief Gather external variables and index maps. * */ template - int GenClassical::evaluateResidual() + void GenClassical::gatherExternalVariables() { - wb_[0] = Vr(); - wb_[1] = Vi(); + 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 generator model. + * + */ + template + int GenClassical::evaluateInternalResidual() + { + gatherExternalVariables(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); - evaluateInternalResidual(y, yp, wb_.data(), f); - evaluateBusResidual(y, yp, wb_.data(), h_.data()); + evaluateInternalResidual(y, yp, y_ext_.data(), f); + f_.setDataUpdated(); + + return 0; + } + + /** + * \brief External residual contributions to the bus. + * + */ + template + int GenClassical::evaluateExternalResidual() + { + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + evaluateExternalResidual(y, yp, y_ext_.data(), f_ext_.data()); - 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 for the generator model. + * + */ + 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 cfe031300..fd857ad82 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -90,6 +90,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 1edc3f62c..e17ccb658 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -961,20 +961,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_) { @@ -983,9 +977,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 6afbd6e98..5be539a2e 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); @@ -377,6 +379,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 8ac45f8d8..0d4c4b77e 100644 --- a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp @@ -1,7 +1,9 @@ #pragma once +#include #include #include +#include #include #include @@ -99,8 +101,31 @@ namespace GridKit // Jacobian via Enzyme auto enzyme_jacobian = EnzymeJacobian(R, X, status); + success *= dependency_tracking_jacobian.size() == enzyme_jacobian.size(); + + const auto remove_zeros = [](auto& jacobian) + { + for (auto& row : jacobian) + { + 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])); } @@ -145,6 +170,11 @@ 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(); + residual_y.insert( + residual_y.end(), + bus_residual_y_view.getData(), + bus_residual_y_view.getData() + bus_residual_y_view.getSize()); // Get d/dy' bus.initialize(); @@ -164,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) diff --git a/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp index 8728a23eb..371433dc6 100644 --- a/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp @@ -1284,6 +1284,17 @@ namespace GridKit return static_cast(variable); } + static constexpr std::array kSignalInputs{ + Ext::IR, + Ext::II, + Ext::P, + Ext::Q, + Ext::FREQ, + Ext::VREF, + Ext::PREF, + Ext::QREF, + Ext::FREQREF}; + struct VariableValue { Vars variable; @@ -1347,8 +1358,9 @@ namespace GridKit void attachRequiredInputs(RealT initial_value = 0.0) { const IdxT external_index_base = repca.size() + bus.size(); - for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + for (const Ext variable : kSignalInputs) { + const auto port = index(variable); input_values_[port] = static_cast(initial_value); input_indices_[port] = external_index_base + static_cast(port); input_nodes_[port].set(&input_values_[port], &input_indices_[port]); @@ -1670,9 +1682,9 @@ namespace GridKit { success = false; } - for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + for (const Ext variable : kSignalInputs) { - const auto variable = static_cast(port); + const auto port = index(variable); if (!rowMatches(implicit_defaults.input(variable), explicit_defaults.input(variable), "documented-default signal", @@ -1774,9 +1786,9 @@ namespace GridKit const auto yp_before = copyVector(fixture.repca.yp()); const auto bus_before = copyVector(fixture.bus.y()); std::array inputs_before{}; - for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + for (const Ext variable : kSignalInputs) { - inputs_before[port] = fixture.input(static_cast(port)); + inputs_before[index(variable)] = fixture.input(variable); } bool success = true; @@ -1806,9 +1818,10 @@ namespace GridKit { success = false; } - for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + for (const Ext variable : kSignalInputs) { - if (!valueUnchanged(fixture.input(static_cast(port)), + const auto port = index(variable); + if (!valueUnchanged(fixture.input(variable), inputs_before[port], "external signal", port)) @@ -2297,9 +2310,8 @@ namespace GridKit { bus_y[row].setVariableNumber(kBusVrColumn + row); } - for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + for (const Ext variable : kSignalInputs) { - const auto variable = static_cast(port); fixture.input(variable).setVariableNumber(fixture.inputIndex(variable)); } diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 968ad1462..3a4fda754 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -156,6 +156,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__); } diff --git a/tests/UnitTests/PhasorDynamics/runBusFaultTests.cpp b/tests/UnitTests/PhasorDynamics/runBusFaultTests.cpp index 5a3688fc5..28ca5b661 100644 --- a/tests/UnitTests/PhasorDynamics/runBusFaultTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runBusFaultTests.cpp @@ -12,6 +12,7 @@ int main() result += test.zeroInitialResidual(true); #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian(true); + result += test.jacobian(false); #endif return result.summary();