diff --git a/CHANGELOG.md b/CHANGELOG.md index b8dfddd8..aec17250 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,7 @@ - Added capability to print monitored variables in multiple formats, triggered from `Ida::runSimulation`. - Added IDA statistics object which can be accumulated over multiple simulations. - Minor performance improvements to residual evaluation in PowerElectronics module. +- Added full support for sparse Jacobians obtained with Enzyme in PhasorDynamics. ## v0.1 diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp index b7adffdc..b6285abc 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp @@ -14,6 +14,7 @@ #include #include +#include #include namespace GridKit @@ -266,7 +267,7 @@ namespace GridKit bool is_fixed_; ///< Constant parameter flag. mutable DependencyMap* dependencies_; - static const size_t INVALID_VAR_NUMBER = static_cast(-1); + static const size_t INVALID_VAR_NUMBER = INVALID_INDEX; }; //------------------------------------ diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp new file mode 100644 index 00000000..392a44d4 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp @@ -0,0 +1,171 @@ +/** + * @file DfDwb.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#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 + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DfDwb + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @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,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; //< @todo: Update once sparse storage format changes + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + + /** + * @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[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; //< @todo: Update once sparse storage format changes + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp new file mode 100644 index 00000000..03591b91 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp @@ -0,0 +1,112 @@ +/** + * @file DfDws.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#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 + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DfDws + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @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[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; //< @todo: Update once sparse storage format changes + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + IdxT row = res_indices[static_cast(tup.row)]; + IdxT col = var_indices[static_cast(tup.col)]; + if (col != INVALID_INDEX) + { + rtemp.push_back(row); + ctemp.push_back(col); + valtemp.push_back(tup.val); + } + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp new file mode 100644 index 00000000..3aac93d3 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp @@ -0,0 +1,259 @@ +/** + * @file DfDy.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: Internal Jacobian, df/dy + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DfDy + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @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] alpha - Time derivative jacobian coefficient + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + RealT alpha, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + // df/dy + std::vector> triplets; //< @todo: Update once sparse storage format changes + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + + // df/dy' + triplets.clear(); + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + ctemp.clear(); + rtemp.clear(); + valtemp.clear(); + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + + /** + * @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[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + RealT alpha, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + // df/dy + std::vector> triplets; //< @todo: Update once sparse storage format changes + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + + // df/dy' + triplets.clear(); + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + ctemp.clear(); + rtemp.clear(); + valtemp.clear(); + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp new file mode 100644 index 00000000..234657f4 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp @@ -0,0 +1,102 @@ +/** + * @file DhDwb.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: dh/dwb + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DhDwb + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @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,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; //< @todo: Update once sparse storage format changes + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.axpy(1.0, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp new file mode 100644 index 00000000..a0b60972 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp @@ -0,0 +1,102 @@ +/** + * @file DhDy.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#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 + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DhDy + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @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,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; //< @todo: Update once sparse storage format changes + 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, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiaation 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); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp b/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp new file mode 100644 index 00000000..4edbcd88 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp @@ -0,0 +1,42 @@ +/** + * @file EnzymeDefinitions.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include + +/** + * @brief Enzyme constants for activity analysis + * + */ +extern int enzyme_dup; +extern int enzyme_const; +extern int enzyme_dupnoneed; + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme fwddiff template + * + * @details This is a core templated intrinsic that the Enzyme pass will + * use to perform automatic differenciation. We define the template here + * so it can later be used in different places. + * + * @tparam T - return type + * @tparam ModelT - model type + */ + template + extern T __enzyme_fwddiff(void*, ModelT...) noexcept; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp new file mode 100644 index 00000000..646bbf79 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp @@ -0,0 +1,130 @@ +/** + * @file LowerSparseStorage.hpp + * + * @details This file contains functions used by Enzyme to store sparse Jacobians. + * + * @todo Replace this inner/lower storage by a more CSR-efficient approach. + * + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme todense template + * + * @details This is used by Enzyme's auto sparsity analysis. It internally maps dense storage + * to sparse ones (where structural zeros are not kept). + * + * @tparam T - return type + */ + template + extern T __enzyme_todense(void*...) noexcept; + + /** + * @brief Enzyme inner sparse storage in triplet format + * + * @tparam ScalarT - scalar data type + */ + template + struct Triple + { + size_t row; + size_t col; + ScalarT val; + Triple(Triple&&) = default; + + Triple(size_t row, size_t col, ScalarT val) + : row(row), + col(col), + val(val) + { + } + }; + + /** + * @brief Enzyme sparse accumulation for float + * + * @todo Separate sparsity analyis from value storage for a known sparsity structure. + * The current implementation always performs sparsity analysis. + */ + [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(size_t row, size_t col, float val, std::vector>& triplets) + { + triplets.emplace_back(row, col, val); + } + + /** + * @brief Enzyme sparse accumulation for double + * + * @todo Separate sparsity analyis from value storage for a known sparsity structure. + * The current implementation always performs sparsity analysis. + */ + [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(size_t row, size_t col, double val, std::vector>& triplets) + { + triplets.emplace_back(row, col, val); + } + + /** + * @brief Enzyme sparse store + * + * @details This takes in a row, column and value and stores them in a vector of triplets + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static void sparse_store(ScalarT val, size_t idx, size_t i, std::vector>& triplets) + { + if (val == 0.0) + return; + idx /= sizeof(ScalarT); + if constexpr (sizeof(ScalarT) == 4) + inner_storeflt(idx, i, val, triplets); + else + inner_storedbl(idx, i, val, triplets); + } + + /** + * @brief Enzyme sparse load + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static ScalarT sparse_load(size_t, size_t, std::vector>&) + { + return 0.0; + } + + /** + * @brief Enzyme identity store + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static void ident_store(ScalarT, size_t, size_t) + { + assert(0 && "should never store"); + } + + /** + * @brief Enzyme identity load + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static ScalarT ident_load(size_t idx, size_t i) + { + idx /= sizeof(ScalarT); + return (ScalarT) (idx == i); + } + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp new file mode 100644 index 00000000..09d67f08 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp @@ -0,0 +1,188 @@ +/** + * @file ModelWrappers.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Model member function parameter keys + * + */ + 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 + }; + + /** + * @brief Template definition for wrapper around residual methods inside model classes + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * + */ + template + struct ModelWrapper + { + }; + + /** + * @brief Residual wrapper partial template specialization for InternalResidual + * + */ + template + struct ModelWrapper + { + /** + * @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] f - Internal residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* f) + { + model->evaluateInternalResidual(y, yp, wb, f); + } + }; + + /** + * @brief Residual wrapper partial template specialization for InternalResidualWithSignal + * + */ + template + struct ModelWrapper + { + /** + * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* ws, ScalarT* f) + { + model->evaluateInternalResidual(y, yp, wb, ws, f); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual + * + */ + template + struct ModelWrapper + { + /** + * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual11 (branch member function) + * + */ + template + struct ModelWrapper + { + /** + * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual11(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual12 (branch member function) + * + */ + template + struct ModelWrapper + { + /** + * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual12(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual21 (branch member function) + * + */ + template + struct ModelWrapper + { + /** + * @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 + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual21(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual22 (branch member function) + * + */ + template + struct ModelWrapper + { + /** + * @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, ScalarT* y, ScalarT* yp, 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 new file mode 100644 index 00000000..9a9df277 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp @@ -0,0 +1,13 @@ +/** + * @file SparseWrapper.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include diff --git a/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp b/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp deleted file mode 100644 index 3bba002f..00000000 --- a/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp +++ /dev/null @@ -1,721 +0,0 @@ -/** - * @file SparseWrapper.hpp - * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) - * - */ - -#pragma once - -#include -#include -#include -#include - -#include -#include - -/** - * @brief Enzyme constants for activity analysis - * - */ -extern int enzyme_dup; -extern int enzyme_const; -extern int enzyme_dupnoneed; - -namespace GridKit -{ - namespace Enzyme - { - namespace Sparse - { - /** - * @brief Model member function parameter keys - * - */ - 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 - }; - - /** - * @brief Template definition for wrapper around residual methods inside model classes - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * - */ - template - struct Wrapper - { - }; - - /** - * @brief Residual wrapper partial template specialization for InternalResidual - * - */ - template - struct Wrapper - { - /** - * @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] f - Internal residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* f) - { - model->evaluateInternalResidual(y, yp, wb, f); - } - }; - - /** - * @brief Residual wrapper partial template specialization for InternalResidualWithSignal - * - */ - template - struct Wrapper - { - /** - * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* ws, ScalarT* f) - { - model->evaluateInternalResidual(y, yp, wb, ws, f); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual - * - */ - template - struct Wrapper - { - /** - * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual11 - * - */ - template - struct Wrapper - { - /** - * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual11(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual12 - * - */ - template - struct Wrapper - { - /** - * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual12(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual21 - * - */ - template - struct Wrapper - { - /** - * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual21(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual22 - * - */ - template - struct Wrapper - { - /** - * @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, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual22(y, yp, wb, h); - } - }; - - /** - * @brief Enzyme fwddiff template - * - * @tparam T - return type - * @tparam ModelT - model type - */ - template - extern T __enzyme_fwddiff(void*, ModelT...) noexcept; - - /** - * @brief Enzyme todense template - * - * @tparam T - return type - */ - template - extern T __enzyme_todense(void*...) noexcept; - - /** - * @brief Enzyme sparse storage in triplet format - * - * @tparam ScalarT - scalar data type - */ - template - struct Triple - { - size_t row; - size_t col; - ScalarT val; - Triple(Triple&&) = default; - - Triple(size_t row, size_t col, ScalarT val) - : row(row), - col(col), - val(val) - { - } - }; - - /** - * @brief Enzyme sparse accumulation for float - * - */ - [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(size_t row, size_t col, float val, std::vector>& triplets) - { - triplets.emplace_back(row, col, val); - } - - /** - * @brief Enzyme sparse accumulation for double - * - */ - [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(size_t row, size_t col, double val, std::vector>& triplets) - { - triplets.emplace_back(row, col, val); - } - - /** - * @brief Enzyme sparse store - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static void sparse_store(ScalarT val, size_t idx, size_t i, std::vector>& triplets) - { - if (val == 0.0) - return; - idx /= sizeof(ScalarT); - if constexpr (sizeof(ScalarT) == 4) - inner_storeflt(idx, i, val, triplets); - else - inner_storedbl(idx, i, val, triplets); - } - - /** - * @brief Enzyme sparse load - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static ScalarT sparse_load(size_t, size_t, std::vector>&) - { - return 0.0; - } - - /** - * @brief Enzyme identity store - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static void ident_store(ScalarT, size_t, size_t) - { - assert(0 && "should never store"); - } - - /** - * @brief Enzyme identity load - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static ScalarT ident_load(size_t idx, size_t i) - { - idx /= sizeof(ScalarT); - return (ScalarT) (idx == i); - } - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Internal Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct InternalJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @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,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_res); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - 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, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_dup, - y, - output, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Internal Jacobian with signal variables - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct InternalJacobianWithSignal - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @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,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_res); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - 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, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_dup, - y, - output, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_const, - ws, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: External Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct ExternalJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @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,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_res); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - 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, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_dup, - ws, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - IdxT row = res_indices.at(static_cast(tup.row)); - IdxT col = var_indices.at(static_cast(tup.col)); - if (col != static_cast(-1)) - { - rtemp.push_back(row); - ctemp.push_back(col); - valtemp.push_back(tup.val); - } - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Bus Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct BusJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @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,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_res); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - 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, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.axpy(1.0, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Branch Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct BranchJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @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,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_res); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - 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, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - } // namespace Sparse - } // namespace Enzyme -} // namespace GridKit diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index 7866d8f1..f896068e 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -20,7 +20,7 @@ namespace GridKit * * @tparam ScalarT - scalar data type * - * @param[in] x + * @param[in] x - expected to be of order 1 * @return value of the sigmoid function */ template diff --git a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp index 80f80ecb..1cd2f6e9 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp @@ -292,11 +292,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, COO_Matrix& a) { - if (alpha == 0) - { - return; - } - if (!this->sorted_) { this->sortSparse(); @@ -372,9 +367,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, std::vector r, std::vector c, std::vector v) { - if (alpha == 0) - return; - if (!this->sorted_) { this->sortSparse(); diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index 9bc9f010..93c812ff 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -2,7 +2,6 @@ #include -#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.cpp b/GridKit/Model/PhasorDynamics/Branch/Branch.cpp index cc1dd85f..b2aff31f 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.cpp @@ -20,8 +20,8 @@ namespace GridKit template int Branch::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Branch..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index c917e431..123eeb82 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -52,6 +52,8 @@ namespace GridKit using Component::wb_; using Component::h_; using Component::J_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; @@ -76,10 +78,6 @@ namespace GridKit return 0; } - virtual void updateTime(RealT /* t */, RealT /* a */) override - { - } - void setR(RealT R) { R_ = R; diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp index 630df803..5f384af2 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp @@ -20,8 +20,8 @@ namespace GridKit template int Branch::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Branch..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp index d66948c4..48582956 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "BranchImpl.hpp" @@ -22,68 +22,66 @@ namespace GridKit template int Branch::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Branch..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); // Bus 1 diagonal Jacobian block owned by the bus - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11, - ScalarT, - IdxT>::eval(this, - static_cast(bus1_->size()), - static_cast((bus1_->y()).size()), - bus1_->getResidualIndices(), - bus1_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus1_->y()).data(), - bus1_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11, + ScalarT, + IdxT>::eval(this, + static_cast(bus1_->size()), + static_cast((bus1_->y()).size()), + bus1_->getResidualIndices(), + bus1_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus1_->y()).data(), + bus1_->getJacobian()); // Bus 2 diagonal Jacobian block owned by the bus - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22, - ScalarT, - IdxT>::eval(this, - static_cast(bus2_->size()), - static_cast((bus2_->y()).size()), - bus2_->getResidualIndices(), - bus2_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus2_->y()).data(), - bus2_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22, + ScalarT, + IdxT>::eval(this, + static_cast(bus2_->size()), + static_cast((bus2_->y()).size()), + bus2_->getResidualIndices(), + bus2_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus2_->y()).data(), + bus2_->getJacobian()); // Off-diagonal Jacobian block (Bus2 variables) owned by the branch - GridKit::Enzyme::Sparse::BranchJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12, - ScalarT, - IdxT>::eval(this, - static_cast(bus1_->size()), - static_cast((bus2_->y()).size()), - bus1_->getResidualIndices(), - bus2_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus2_->y()).data(), - J_); - - J_.printMatrix("Branch after BusJacobian12 call"); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12, + ScalarT, + IdxT>::eval(this, + static_cast(bus1_->size()), + static_cast((bus2_->y()).size()), + bus1_->getResidualIndices(), + bus2_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus2_->y()).data(), + J_); // Off-diagonal Jacobian block (Bus1 variables) owned by the branch - GridKit::Enzyme::Sparse::BranchJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21, - ScalarT, - IdxT>::eval(this, - static_cast(bus2_->size()), - static_cast((bus1_->y()).size()), - bus2_->getResidualIndices(), - bus1_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus1_->y()).data(), - J_); - - J_.printMatrix("Branch after BusJacobian21 call"); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21, + ScalarT, + IdxT>::eval(this, + static_cast(bus2_->size()), + static_cast((bus1_->y()).size()), + bus2_->getResidualIndices(), + bus1_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus1_->y()).data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt index 10761d6f..3325ea60 100644 --- a/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt @@ -19,6 +19,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -27,13 +28,19 @@ else() Branch.cpp HEADERS ${_install_headers} + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) endif() gridkit_add_library(phasor_dynamics_branch_dependency_tracking - SOURCES BranchDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + BranchDependencyTracking.cpp + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) # Link to interface target for all components target_link_libraries(phasor_dynamics_components diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.cpp b/GridKit/Model/PhasorDynamics/Bus/Bus.cpp index 7e633dae..0971e250 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Bus::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Bus..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Bus..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp index 3d51ebd1..fe116cf7 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Bus::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Bus..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Bus..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp index 2e3f5a37..d9aa774b 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp @@ -24,14 +24,15 @@ namespace GridKit template int Bus::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Bus..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Bus..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); std::vector rtemp = {residual_indices_.at(0), residual_indices_.at(0), residual_indices_.at(1), residual_indices_.at(1)}; std::vector ctemp = {variable_indices_.at(0), variable_indices_.at(1), variable_indices_.at(0), variable_indices_.at(1)}; std::vector valtemp = {0.0, 0.0, 0.0, 0.0}; J_.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - J_.printMatrix("Bus Jacobian initialization... should be 2x2 with zero values."); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 37b70c01..77fe6ecf 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -86,6 +86,8 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp index ad97593f..a186f02c 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp @@ -22,6 +22,8 @@ namespace GridKit using BusBase::yp_; using BusBase::f_; using BusBase::J_; + using BusBase::variable_indices_; + using BusBase::residual_indices_; public: using RealT = typename BusBase::RealT; diff --git a/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt index b5086ed7..ca318a3b 100644 --- a/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt @@ -38,8 +38,8 @@ gridkit_add_library(phasor_dynamics_bus_dependency_tracking SOURCES BusDependencyTracking.cpp BusInfiniteDependencyTracking.cpp - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PUBLIC GridKit::phasor_dynamics_core) diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 603f9699..12dc04ed 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -5,14 +5,18 @@ #include #include +#include #include #include #include +#include namespace GridKit { namespace PhasorDynamics { + using Log = ::GridKit::Utilities::Logger; + /*! * @brief BusBase model implementation base class. * @@ -127,32 +131,32 @@ namespace GridKit int setVariableIndex(IdxT local_index, IdxT global_index) { - variable_indices_[local_index] = global_index; + variable_indices_[static_cast(local_index)] = global_index; return 0; } IdxT getVariableIndex(IdxT local_index) const { - return variable_indices_.at(local_index); + return variable_indices_[static_cast(local_index)]; } - const std::map& getVariableIndices() const + const std::vector& getVariableIndices() const { return variable_indices_; } int setResidualIndex(IdxT local_index, IdxT global_index) { - residual_indices_[local_index] = global_index; + residual_indices_[static_cast(local_index)] = global_index; return 0; } IdxT getResidualIndex(IdxT local_index) const { - return residual_indices_.at(local_index); + return residual_indices_[static_cast(local_index)]; } - const std::map& getResidualIndices() const + const std::vector& getResidualIndices() const { return residual_indices_; } @@ -160,14 +164,12 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; protected: - IdxT bus_id_{static_cast(-1)}; + IdxT bus_id_{INVALID_INDEX}; - IdxT size_{0}; - IdxT nnz_{0}; - std::map variable_indices_; ///< Map between local and global (system-level) - /// variable indices - std::map residual_indices_; ///< Map between local and global (system-level) - /// residual indices + IdxT size_{0}; + IdxT nnz_{0}; + std::vector variable_indices_; ///< Global (system-level) variable indices + std::vector residual_indices_; ///< Global (system-level) residual indices /// Variable monitor std::unique_ptr monitor_; @@ -179,9 +181,6 @@ namespace GridKit MatrixT J_; - RealT time_; - RealT alpha_; - RealT rtol_; RealT atol_; diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp index 33b46987..9f8124c5 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp @@ -22,8 +22,8 @@ namespace GridKit template int BusFault::evaluateJacobian() { - std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp index 342e932d..6dc1173c 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp @@ -14,8 +14,8 @@ namespace GridKit template int BusFault::evaluateJacobian() { - std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp index 369bf74f..8defda6d 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "BusFaultImpl.hpp" @@ -22,23 +22,23 @@ namespace GridKit template int BusFault::evaluateJacobian() { - std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; if (status_) { - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - bus_->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - bus_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + bus_->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + bus_->getJacobian()); } return 0; diff --git a/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt b/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt index 40f438fb..53b0de98 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt @@ -13,6 +13,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -22,12 +23,18 @@ else() HEADERS ${_install_headers} INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core) endif() gridkit_add_library(phasor_dynamics_bus_fault_dependency_tracking - SOURCES BusFaultDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + BusFaultDependencyTracking.cpp + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core) # Link to interface target for all components target_link_libraries(phasor_dynamics_components diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index cab9f03a..9545732c 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -18,6 +18,7 @@ gridkit_add_library(phasor_dynamics_core SystemModel.hpp SystemModelData.hpp LINK_LIBRARIES + PUBLIC GridKit::definitions PUBLIC GridKit::utilities_logger INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/nlohmann-json/include diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 7fdf3594..882ad575 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -3,12 +3,16 @@ #include #include +#include #include +#include namespace GridKit { namespace PhasorDynamics { + using Log = ::GridKit::Utilities::Logger; + /** * @brief Component model implementation base class. */ @@ -39,15 +43,15 @@ namespace GridKit /// @todo Remove this method. It should be part of DynamicSolver class. virtual bool hasJacobian() override { - return false; + return true; } - // virtual void updateTime(RealT t, RealT a) - // { - // time_ = t; - // alpha_ = a; - // std::cout << "updateTime: t = " << time_ << "\n"; - // } + virtual void updateTime(RealT t, RealT a) override + { + time_ = t; + alpha_ = a; + // std::cout << "updateTime: t = " << time_ << ", alpha = " << alpha_ << "\n"; + } virtual void setTolerances(RealT& rtol, RealT& atol) const override { @@ -119,44 +123,42 @@ namespace GridKit int setVariableIndex(IdxT local_index, IdxT global_index) { - variable_indices_[local_index] = global_index; + variable_indices_[static_cast(local_index)] = global_index; return 0; } IdxT& getVariableIndex(IdxT local_index) { - return variable_indices_.at(local_index); + return variable_indices_[static_cast(local_index)]; } - const std::map& getVariableIndices() const + const std::vector& getVariableIndices() const { return variable_indices_; } int setResidualIndex(IdxT local_index, IdxT global_index) { - residual_indices_[local_index] = global_index; + residual_indices_[static_cast(local_index)] = global_index; return 0; } IdxT& getResidualIndex(IdxT local_index) { - return residual_indices_.at(local_index); + return residual_indices_[static_cast(local_index)]; } - const std::map& getResidualIndices() const + const std::vector& getResidualIndices() const { return residual_indices_; } protected: - IdxT gridkit_component_id_{0}; - IdxT size_{0}; - IdxT nnz_{0}; - std::map variable_indices_; ///< Map between local and global (system-level) - /// variable indices - std::map residual_indices_; ///< Map between local and global (system-level) - /// residual indices + IdxT gridkit_component_id_{0}; + IdxT size_{0}; + IdxT nnz_{0}; + std::vector variable_indices_; ///< Global (system-level) variable indices + std::vector residual_indices_; ///< Global (system-level) residual indices std::vector y_; std::vector yp_; @@ -180,7 +182,7 @@ namespace GridKit ------ WARNING: Temporary ------ - The protexted variable mva_system_base_ is temporarily + The protected variable mva_system_base_ is temporarily hard coded. This eventually needs to be configured from the input JSON format, which specifies the system MVA base. diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt index b043b66a..52d50ded 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt @@ -8,16 +8,32 @@ set(_install_headers Ieeet1.hpp Ieeet1Data.hpp) -gridkit_add_library(phasor_dynamics_exciter_ieeet1 - SOURCES - Ieeet1.cpp - HEADERS - ${_install_headers} - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include - LINK_LIBRARIES - GridKit::phasor_dynamics_core - GridKit::phasor_dynamics_signal) +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library(phasor_dynamics_exciter_ieeet1 + SOURCES + Ieeet1Enzyme.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + PUBLIC GridKit::phasor_dynamics_signal + PRIVATE ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) +else() + gridkit_add_library(phasor_dynamics_exciter_ieeet1 + SOURCES + Ieeet1.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + GridKit::phasor_dynamics_core + GridKit::phasor_dynamics_signal) +endif() gridkit_add_library(phasor_dynamics_exciter_ieeet1_dependency_tracking SOURCES diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp index f626ddea..da21cd77 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp @@ -15,6 +15,21 @@ namespace GridKit { namespace Exciter { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeet1::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + // Available template instantiations template class Ieeet1; template class Ieeet1; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index fffb5674..95066f4d 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -75,6 +75,10 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; + using Component::wb_; + using Component::J_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; @@ -100,10 +104,6 @@ namespace GridKit int evaluateResidual() override; int evaluateJacobian() override; - void updateTime(RealT /* t */, RealT /* a */) override - { - } - /// Get the `ComponentSignals` from this `Ieeet1` auto getSignals() -> ComponentSignals signals_; @@ -162,6 +164,10 @@ namespace GridKit /// Associate variable getter functions with enum values void initializeMonitor(); + + /* Local copies of signal variables */ + std::vector ws_; + std::vector ws_indices_; }; } // namespace Exciter diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp index 7b5c1109..9346fabb 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp @@ -15,6 +15,21 @@ namespace GridKit { namespace Exciter { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeet1::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + // Available template instantiations template class Ieeet1; template class Ieeet1; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp new file mode 100644 index 00000000..8ffead08 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp @@ -0,0 +1,84 @@ +/** + * @file Ieeet1Enzyme.cpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#include + +#include "Ieeet1Impl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeet1::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_); + + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + ws_.data(), + J_); + + GridKit::Enzyme::Sparse::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + this->getResidualIndices(), + ws_indices_, + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); + + return 0; + } + + // Available template instantiations + template class Ieeet1; + template class Ieeet1; + + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 6e935124..f2cffa30 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -122,6 +122,17 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); + + // Resize bus data + wb_.resize(2); + + // Resize signal variable data + ws_.resize(1); + ws_indices_.resize(1); + ws_[0] = 0.0; + ws_indices_[0] = INVALID_INDEX; // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) @@ -189,7 +200,7 @@ namespace GridKit // Terminal Voltage ScalarT vreal = bus_->Vr(); ScalarT vimag = bus_->Vi(); - Ec_ = std::sqrt(vreal * vreal + vimag * vimag); + ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); // Derived from External initial values ScalarT vr = Ke_ * efd0; @@ -197,10 +208,10 @@ namespace GridKit ScalarT vtr = Ke_ / Ka_ * efd0; // Vref (setpoint = terminal + error) - vref_ = Ec_ + vtr; + vref_ = Ec + vtr; // IVP for Internal Variables - y_[0] = Ec_; // y0 - vts - Sensed term volt + y_[0] = Ec; // y0 - vts - Sensed term volt y_[1] = vr; // y1 - vr - Voltage reg y_[2] = efd0; // y2 - efdp - Efd pre mult y_[3] = vfx; // y3 - vfx - Exciter feedback @@ -249,84 +260,94 @@ namespace GridKit } /** - * @brief Residuals of system equations + * @brief Internal Residual * */ template - int Ieeet1::evaluateResidual() + __attribute__((always_inline)) inline int Ieeet1::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) { - // Input Variables - ScalarT omega{0}; - - // Set Input Variables - // Meta PR Note: This seems to be very slow, - // but I see how read/write ownership may require this - // - // I believe implementing the equivalent to signal->read() - // at the system level would address this, by routing - // external signals into a generic inputs_ vector - // at the same time as the internal state values y_ - // are recieved from IDA. - if (signals_.template isAttached()) - { - omega = signals_.template readExternalVariable(); - } - // Read E comp (terminal voltage, unless compensation impedance) - ScalarT vreal = bus_->Vr(); - ScalarT vimag = bus_->Vi(); - Ec_ = std::sqrt(vreal * vreal + vimag * vimag); + ScalarT vreal = wb[0]; + ScalarT vimag = wb[1]; + ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); // Read Internal Variables - ScalarT vts = y_[0]; // y0 - Sensed term volt - ScalarT vr = y_[1]; // y1 - Voltage reg - ScalarT efdp = y_[2]; // y2 - Efd pre mult - ScalarT vfx = y_[3]; // y3 - Exciter feedback - ScalarT vtr = y_[4]; // y4 - Term Volt Err - ScalarT vf = y_[5]; // y5 - Feedback volt - ScalarT ve = y_[6]; // y6 - Excit. Cntrl Volt - ScalarT efd = y_[7]; // y7 - Efd - ScalarT ksat = y_[8]; // y8 - Saturation + ScalarT vts = y[0]; // y0 - Sensed term volt + ScalarT vr = y[1]; // y1 - Voltage reg + ScalarT efdp = y[2]; // y2 - Efd pre mult + ScalarT vfx = y[3]; // y3 - Exciter feedback + ScalarT vtr = y[4]; // y4 - Term Volt Err + ScalarT vf = y[5]; // y5 - Feedback volt + ScalarT ve = y[6]; // y6 - Excit. Cntrl Volt + ScalarT efd = y[7]; // y7 - Efd + ScalarT ksat = y[8]; // y8 - Saturation // Read Internal Derivatives - ScalarT vts_dot = yp_[0]; - ScalarT vr_dot = yp_[1]; - ScalarT efdp_dot = yp_[2]; - ScalarT vfx_dot = yp_[3]; + ScalarT vts_dot = yp[0]; + ScalarT vr_dot = yp[1]; + ScalarT efdp_dot = yp[2]; + ScalarT vfx_dot = yp[3]; + + // Set signal variable aliases + ScalarT omega = ws[0]; // The 'pre-limit' derivative of Pv - ScalarT f = -vr + Ka_ * vtr; - ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, f); + ScalarT func = -vr + Ka_ * vtr; + ScalarT func_normalized = func * 12.0 / 240.0; // Arbitrary normalization to use standardized sigmoid + ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, func_normalized); // Internal Differential Equations - f_[0] = -vts_dot + (Ec_ - vts) / Tr_; - f_[1] = -vr_dot + vr_ind * f / Ta_; - f_[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_; - f_[3] = -vfx_dot + vf / Tf_; + f[0] = -vts_dot + (Ec - vts) / Tr_; + f[1] = -vr_dot + vr_ind * func / Ta_; + f[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_; + f[3] = -vfx_dot + vf / Tf_; // Internal Algebraic Equations - f_[4] = -vtr + vref_ + vUEL_ + vOEL_ + vS_ - vf - vts; - f_[5] = -vf + (efdp * Kf_) / Tf_ - vfx; - f_[6] = -ve + ksat * efdp; - f_[7] = -efd + efdp + omega * efdp * Ispdlim_; + f[4] = -vts + vref_ + vUEL_ + vOEL_ + vS_ - vtr - vf; + f[5] = -vf + (efdp * Kf_) / Tf_ - vfx; + f[6] = -ve + ksat * efdp; + f[7] = -efd + efdp + omega * efdp * Ispdlim_; ScalarT efd_sat = (efdp - SA_) * (Math::sigmoid(efdp - SA_)); - f_[8] = -ksat + SB_ * efd_sat * efd_sat; + f[8] = -ksat + SB_ * efd_sat * efd_sat; return 0; } /** - * @brief Jacobian evaluation not implemented yet + * @brief Residual evaluation * - * @tparam ScalarT - Scalar data type - * @tparam IdxT - Index data type - * @return int - error code, 0 = success */ template - int Ieeet1::evaluateJacobian() + int Ieeet1::evaluateResidual() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + // Set Input Variables + // Meta PR Note: This seems to be very slow, + // but I see how read/write ownership may require this + // + // I believe implementing the equivalent to signal->read() + // at the system level would address this, by routing + // external signals into a generic inputs_ vector + // at the same time as the internal state values y_ + // are recieved from IDA. + if (signals_.template isAttached()) + { + ws_[0] = signals_.template readExternalVariable(); + ws_indices_[0] = signals_.template readExternalVariableIndex(); + } + + // Bus voltages + wb_[0] = bus_->Vr(); + wb_[1] = bus_->Vi(); + + // Residual evaluation + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp index 22743387..048f058a 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp @@ -21,7 +21,8 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index 95bb6fda..320e8c29 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -70,6 +70,8 @@ namespace GridKit using Component::wb_; using Component::h_; using Component::J_; + using Component::variable_indices_; + using Component::residual_indices_; using RealT = typename Component::RealT; using model_data_type = Tgov1Data; @@ -91,10 +93,6 @@ namespace GridKit // Still to be implemented int evaluateJacobian() override; - void updateTime(RealT /* t */, RealT /* a */) override - { - } - /// Get the `ComponentSignals` from this `Tgov1` auto getSignals() -> ComponentSignals ws_; - std::map ws_indices_; + std::vector ws_indices_; }; } // namespace Governor diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp index 62dc3af8..63476a79 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp @@ -21,7 +21,8 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp index 57ecbb12..aee05e06 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "Tgov1Impl.hpp" @@ -24,40 +24,39 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Tgov1..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - GridKit::Enzyme::Sparse::InternalJacobianWithSignal, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - this->getResidualIndices(), - this->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); + J_.zeroMatrix(); - J_.printMatrix("Tgov1 internal Jacobian"); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_); - GridKit::Enzyme::Sparse::ExternalJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - this->getResidualIndices(), - ws_indices_, - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); - - J_.printMatrix("Tgov1 Jacobian after signal evaluation"); + GridKit::Enzyme::Sparse::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + this->getResidualIndices(), + ws_indices_, + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index f6add0df..2d7c7d87 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -140,11 +140,14 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); - // Resize external variable data + // Resize signal variable data ws_.resize(1); + ws_indices_.resize(1); ws_[0] = 0.0; - ws_indices_[0] = static_cast(-1); + ws_indices_[0] = INVALID_INDEX; // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) @@ -238,7 +241,7 @@ namespace GridKit ScalarT* y, ScalarT* yp, [[maybe_unused]] ScalarT* wb, - [[maybe_unused]] ScalarT* ws, + ScalarT* ws, ScalarT* f) { // Read Internal Variables @@ -250,7 +253,7 @@ namespace GridKit ScalarT ptx_dot = yp[0]; ScalarT pv_dot = yp[1]; - // Set external variable aliases + // Set signal variable aliases ScalarT omega = ws[0]; // The 'pre-limit' derivative of Pv diff --git a/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt index ef86891b..91fb5164 100644 --- a/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt @@ -19,6 +19,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -27,13 +28,19 @@ else() Load.cpp INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core HEADERS ${_install_headers}) endif() gridkit_add_library(phasor_dynamics_load_dependency_tracking - SOURCES LoadDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + LoadDependencyTracking.cpp + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core) # Link to interface target for all components target_link_libraries(phasor_dynamics_components diff --git a/GridKit/Model/PhasorDynamics/Load/Load.cpp b/GridKit/Model/PhasorDynamics/Load/Load.cpp index 8f130cba..324b0af8 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.cpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Load::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Load..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Load..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/Load.hpp b/GridKit/Model/PhasorDynamics/Load/Load.hpp index 091853b5..8ffe2948 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.hpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.hpp @@ -39,6 +39,8 @@ namespace GridKit using Component::tag_; using Component::wb_; using Component::h_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; @@ -63,10 +65,6 @@ namespace GridKit return 0; } - virtual void updateTime(RealT /* t */, RealT /* a */) override - { - } - public: void setR(RealT R) { diff --git a/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp index 8219abca..37e9f023 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Load::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Load..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Load..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp index b5c1904a..b5309c6f 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "LoadImpl.hpp" @@ -22,21 +22,21 @@ namespace GridKit template int Load::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Load..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Load..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - bus_->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - bus_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + bus_->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + bus_->getJacobian()); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp index 7bd6e91b..861b71dc 100644 --- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp +++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp @@ -67,7 +67,7 @@ namespace GridKit virtual void updateTime(RealT /* t */, RealT /* a */) override { - // No time to update in bus models + // No time to update in signal nodes } virtual void setTolerances(RealT& rtol, RealT& atol) const override @@ -131,7 +131,7 @@ namespace GridKit IdxT signal_id_{0}; protected: - const IdxT bus_id_{static_cast(-1)}; + const IdxT bus_id_{INVALID_INDEX}; IdxT size_{0}; IdxT nnz_{0}; @@ -144,9 +144,6 @@ namespace GridKit MatrixT J_; - RealT time_; - RealT alpha_; - RealT rtol_; RealT atol_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp index f9f38948..efcea4fc 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp @@ -21,8 +21,8 @@ namespace GridKit template int Genrou::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp index 8c9aaa61..1ea500dd 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp @@ -82,6 +82,8 @@ namespace GridKit using Component::h_; using Component::J_; using Component::mva_system_base_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; @@ -133,10 +135,6 @@ namespace GridKit // Still to be implemented int evaluateJacobian() override; - void updateTime(RealT /* t */, RealT /* a */) override - { - } - // Temporary access functions for governor // Should be abstracted ScalarT getSpeed(); @@ -237,9 +235,9 @@ 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 external variables */ + /* Local copies of signal variables */ std::vector ws_; - std::map ws_indices_; + std::vector ws_indices_; /// Variable monitor std::unique_ptr monitor_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp index 0419269d..8d4c004c 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp @@ -14,8 +14,8 @@ namespace GridKit template int Genrou::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp index b6984666..8360f1c8 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "GenrouImpl.hpp" @@ -22,53 +22,79 @@ namespace GridKit template int Genrou::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - GridKit::Enzyme::Sparse::InternalJacobianWithSignal, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - this->getResidualIndices(), - this->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); + J_.zeroMatrix(); - J_.printMatrix("Genrou internal Jacobian"); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_); - GridKit::Enzyme::Sparse::ExternalJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - this->getResidualIndices(), - ws_indices_, - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); - J_.printMatrix("Genrou Jacobian after signal evaluation"); + GridKit::Enzyme::Sparse::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + this->getResidualIndices(), + ws_indices_, + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - bus_->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - bus_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + bus_->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + bus_->getJacobian()); + + GridKit::Enzyme::Sparse::DhDy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + y_.size(), + bus_->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index 7f3ae05f..56fae73b 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -310,19 +310,23 @@ namespace GridKit int Genrou::allocate() { // Resize component model data - f_.resize(static_cast(size_)); - y_.resize(static_cast(size_)); - yp_.resize(static_cast(size_)); - tag_.resize(static_cast(size_)); - - // Resize coupling data + auto size = static_cast(size_); + f_.resize(size); + y_.resize(size); + yp_.resize(size); + tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); + + // Resize bus data wb_.resize(2); h_.resize(2); - // Resize external variable data + // Resize signal variable data ws_.resize(2); - ws_indices_[0] = static_cast(-1); - ws_indices_[1] = static_cast(-1); + ws_indices_.resize(2); + ws_indices_[0] = INVALID_INDEX; + ws_indices_[1] = INVALID_INDEX; // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) @@ -513,7 +517,7 @@ namespace GridKit ScalarT vr = wb[0]; ScalarT vi = wb[1]; - // Set external variable aliases + // Set signal variable aliases ScalarT pmech = ws[0]; ScalarT efd = ws[1]; @@ -551,7 +555,7 @@ namespace GridKit */ template __attribute__((always_inline)) inline int Genrou::evaluateBusResidual( - [[maybe_unused]] ScalarT* y, + ScalarT* y, [[maybe_unused]] ScalarT* yp, ScalarT* wb, ScalarT* h) diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt index ef1c782b..ef03af11 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt @@ -18,6 +18,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -26,13 +27,19 @@ else() GenClassical.cpp HEADERS ${_install_headers} + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) endif() gridkit_add_library(phasor_dynamics_gen_classical_dependency_tracking - SOURCES GenClassicalDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + GenClassicalDependencyTracking.cpp + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) # Link to interface target for all components target_link_libraries(phasor_dynamics_components diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index 52e06f95..b13b6624 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -15,8 +15,8 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index 7a6d69df..d1ca4c29 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -46,6 +46,8 @@ namespace GridKit using Component::h_; using Component::J_; using Component::mva_system_base_; + using Component::variable_indices_; + using Component::residual_indices_; public: using bus_type = BusBase; @@ -79,10 +81,6 @@ namespace GridKit // Still to be implemented int evaluateJacobian() override; - void updateTime(RealT /* t */, RealT /* a */) override - { - } - void setPmech(RealT pmech) { pmech_set_ = pmech; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp index d6f4885f..1fa1a703 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp index c326316e..81f8d0a0 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "GenClassicalImpl.hpp" @@ -22,23 +22,50 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; - - GridKit::Enzyme::Sparse::InternalJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - this->getResidualIndices(), - this->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - J_); - - J_.printMatrix("GenClassical internal Jacobian"); + Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + alpha_, + J_); + + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_); + + GridKit::Enzyme::Sparse::DhDy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + y_.size(), + bus_->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index 58ba7039..bec605fc 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -170,6 +170,8 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); // Resize coupling data wb_.resize(2); diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 2d787a02..a2ce708e 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -4,13 +4,13 @@ #include #include +#include #include #include #include #include #include #include -#include // Include all components #include @@ -19,8 +19,6 @@ namespace GridKit { namespace PhasorDynamics { - using Log = ::GridKit::Utilities::Logger; - /** * @brief Prototype for a system model class * @@ -51,6 +49,8 @@ namespace GridKit using PhasorDynamics::Component::J_; using PhasorDynamics::Component::rel_tol_; using PhasorDynamics::Component::abs_tol_; + using PhasorDynamics::Component::variable_indices_; + using PhasorDynamics::Component::residual_indices_; public: /** @@ -327,6 +327,8 @@ namespace GridKit yp_.resize(size_); f_.resize(size_); tag_.resize(size_); + variable_indices_.resize(size_); + residual_indices_.resize(size_); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) @@ -335,6 +337,7 @@ namespace GridKit this->setResidualIndex(j, j); } + // Verify component configuration int errorCount = this->verify(); if (errorCount > 0) { @@ -342,6 +345,17 @@ namespace GridKit throw std::runtime_error("SystemModel allocation failed"); } + // Perform an initial Jacobian evaluation for sparse Jacobians, such that + // the dynamic solver can querry the NNZ value when it is configured. + // @todo Replace with a sparsity analysis that sets the NNZ and allocates the Jacobian + // without needing the Jacobian values. + if (hasJacobian()) + { + initialize(); + evaluateResidual(); + evaluateJacobian(); + } + return 0; } @@ -365,14 +379,31 @@ namespace GridKit } /** - * @brief Assume that jacobian is not available + * @brief Check components for Jacobian availability * * @return true * @return false */ bool hasJacobian() override { - return false; + bool has_jacobian = false; +#ifdef GRIDKIT_ENABLE_ENZYME + has_jacobian = true; + for (const auto& component : components_) + { + has_jacobian = has_jacobian && component->hasJacobian(); + } + + if (!has_jacobian) + { + Log::warning() << "GritKit was built with Enzyme, but some models don't have Jacobians available. " + << "Falling back to dense Jacobians in PhasorDynamics.\n"; + } +#else + Log::warning() << "GritKit was not built with Enzyme. " + << "Falling back to dense Jacobians in PhasorDynamics.\n"; +#endif + return has_jacobian; } /** @@ -573,6 +604,7 @@ namespace GridKit */ int evaluateJacobian() override { + J_.zeroMatrix(); std::vector ctemp{}; std::vector rtemp{}; std::vector valtemp{}; @@ -617,6 +649,8 @@ namespace GridKit J_.setValues(rtemp, ctemp, valtemp); + // J_.printMatrix("System Jacobian"); + return 0; } @@ -659,7 +693,8 @@ namespace GridKit */ void updateTime(RealT t, RealT a) override { - this->time_ = t; + time_ = t; + alpha_ = a; for (const auto& component : components_) { component->updateTime(t, a); diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index 55bcf3f1..cbf8a613 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -40,7 +40,7 @@ namespace GridKit this->alpha_ = a; } - bool hasJacobian() + bool hasJacobian() override { return true; } diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 1bffd0d8..f5796fc8 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -409,7 +410,7 @@ namespace GridKit } private: - static constexpr IdxT neg1_ = static_cast(-1); + static constexpr IdxT neg1_ = INVALID_INDEX; std::vector components_; diff --git a/GridKit/Model/VariableMonitor.hpp b/GridKit/Model/VariableMonitor.hpp index ce7df4f4..c50654ae 100644 --- a/GridKit/Model/VariableMonitor.hpp +++ b/GridKit/Model/VariableMonitor.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/VariableMonitorController.hpp b/GridKit/Model/VariableMonitorController.hpp index 08bc9a77..79db154c 100644 --- a/GridKit/Model/VariableMonitorController.hpp +++ b/GridKit/Model/VariableMonitorController.hpp @@ -164,7 +164,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { printFullHeader(*sink.os, sink.format); }, + { this->printFullHeader(*sink.os, sink.format); }, sink); } } @@ -256,7 +256,7 @@ namespace GridKit { std::visit([this](auto&& sink) { - printFull(*sink.os, sink.format); + this->printFull(*sink.os, sink.format); using T = std::remove_cvref_t; if constexpr (std::is_same_v>) { @@ -295,7 +295,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { printFullFooter(*sink.os, sink.format); }, + { this->printFullFooter(*sink.os, sink.format); }, sink); } } @@ -403,7 +403,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { sink.start(); }, + { (void) this; sink.start(); }, sink); } } @@ -416,7 +416,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { sink.stop(); }, + { (void) this; sink.stop(); }, sink); } } diff --git a/GridKit/Model/VariableMonitorImpl.hpp b/GridKit/Model/VariableMonitorImpl.hpp index fdcca543..eda73332 100644 --- a/GridKit/Model/VariableMonitorImpl.hpp +++ b/GridKit/Model/VariableMonitorImpl.hpp @@ -97,7 +97,7 @@ namespace GridKit template void set(VariableEnum v, FuncT f) { - f_[enum_integer(v)] = ValuePrinter{f}; + f_[static_cast(enum_integer(v))] = ValuePrinter{f}; } private: @@ -166,7 +166,7 @@ namespace GridKit { for (auto v : variables_) { - os << csv.delim << f_[enum_integer(v)]; + os << csv.delim << f_[static_cast(enum_integer(v))]; } } @@ -176,7 +176,7 @@ namespace GridKit void print(std::ostream& os, VariableEnum v, Json) const { os << indent_ << std::quoted(enum_name(v)) << ": " - << f_[enum_integer(v)] << ",\n"; + << f_[static_cast(enum_integer(v))] << ",\n"; } void print(std::ostream& os, Json) const override @@ -209,7 +209,7 @@ namespace GridKit */ void print(std::ostream& os, VariableEnum v, Yaml) const { - os << indent_ << enum_name(v) << ": " << f_[enum_integer(v)] + os << indent_ << enum_name(v) << ": " << f_[static_cast(enum_integer(v))] << '\n'; } diff --git a/GridKit/Solver/Dynamic/CMakeLists.txt b/GridKit/Solver/Dynamic/CMakeLists.txt index c8e77529..0cceed1f 100644 --- a/GridKit/Solver/Dynamic/CMakeLists.txt +++ b/GridKit/Solver/Dynamic/CMakeLists.txt @@ -19,7 +19,8 @@ if (GRIDKIT_ENABLE_SUNDIALS_SPARSE) PUBLIC SUNDIALS::nvecserial PUBLIC SUNDIALS::idas PUBLIC SUNDIALS::sunlinsolklu - PUBLIC GridKit::definitions) + PUBLIC GridKit::definitions + PUBLIC GridKit::utilities_logger) else() gridkit_add_library(solvers_dyn SOURCES @@ -29,5 +30,6 @@ else() LINK_LIBRARIES PUBLIC SUNDIALS::nvecserial PUBLIC SUNDIALS::idas - PUBLIC GridKit::definitions) + PUBLIC GridKit::definitions + PUBLIC GridKit::utilities_logger) endif() diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index a05ba1ca..aec2d9b0 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -142,15 +142,13 @@ namespace AnalysisManager this->configureLinearSolverDense(); } #else + /// Todo - Improve error handling capabilities and hasJacobian_ ownership if (model_->hasJacobian()) { - /// Todo - Improve error handling capabilities and hasJacobian_ ownership - throw std::runtime_error("SUNDIALS is not configured with KLU, but the model has a (sparse) Jacobian."); - } - else - { - this->configureLinearSolverDense(); + Log::warning() << "SUNDIALS is not configured with KLU, but the model has a (sparse) Jacobian. " + << "Falling back to dense Jacobian.\n"; } + this->configureLinearSolverDense(); #endif return retval; @@ -275,9 +273,12 @@ namespace AnalysisManager if (tag_) initType = IDA_YA_YDP_INIT; - retval = IDACalcIC(solver_, initType, 0.1); + retval = IDACalcIC(solver_, initType, t0 + 0.1); checkOutput(retval, "IDACalcIC"); + retval = IDAGetConsistentIC(solver_, yy_, yp_); + checkOutput(retval, "IDAGetConsistentIC"); + copyVec(yy_, model_->y()); copyVec(yp_, model_->yp()); } diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 650e2885..20ff56aa 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -19,11 +19,14 @@ #include #include +#include namespace AnalysisManager { namespace Sundials { + using Log = ::GridKit::Utilities::Logger; + struct IdaStats { long int num_steps_ = 0; diff --git a/buildsystem/spack_repo/gridkit/packages/enzyme/package.py b/buildsystem/spack_repo/gridkit/packages/enzyme/package.py index a57c5cd9..930848d4 100644 --- a/buildsystem/spack_repo/gridkit/packages/enzyme/package.py +++ b/buildsystem/spack_repo/gridkit/packages/enzyme/package.py @@ -3,6 +3,7 @@ from spack.package import * class Enzyme(BuiltinEnzyme): + version("0.0.216", sha256="c41f1a1286934f5155614fde0f3c63e0f07bc3b67b8f2ff728540bc862cb8f7a") version("0.0.206", sha256="600fd2db370fb40abb6411e0e80df524aea03f2c1ad50a2765ecaab9e1115c77") version("0.0.196", sha256="2b9cfcb7c34e56fc8191423042df06241cf32928eefbb113ac3c5199e3361cb2") version("0.0.186", sha256="125e612df0b6b82b07e1e13218c515bc54e04aa1407e57f4f31d3abe995f4714") diff --git a/cmake/FindEnzyme.cmake b/cmake/FindEnzyme.cmake index fbfaa440..ad19c20a 100644 --- a/cmake/FindEnzyme.cmake +++ b/cmake/FindEnzyme.cmake @@ -12,6 +12,15 @@ Author(s): ]] + +# Error if CMAKE_BUILD_TYPE is not Release or RelWithDebInfo due to an Enzyme bug +if (NOT ((CMAKE_BUILD_TYPE STREQUAL "Release") OR (CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo"))) + message(STATUS "CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}") + message(FATAL_ERROR "Enzyme builds currently only support Release and RelWithDebInfo as \ + CMAKE_BUILD_TYPE due to a bug within the Enzyme library.") +endif() + +# Find Enzyme and necessary programs find_package(Enzyme REQUIRED CONFIG PATHS ${ENZYME_DIR} diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt b/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt index 0322107b..72a04361 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt +++ b/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt @@ -5,4 +5,4 @@ target_link_libraries(TenGenClassical install(TARGETS TenGenClassical RUNTIME DESTINATION bin) # Not used for tesing for now. -# add_test(NAME GenrouTest3_genclassical COMMAND $) +add_test(NAME TenGenClassical COMMAND $) diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp index 05c29f6b..50c9dd2d 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp @@ -123,7 +123,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& yval = sys.y(); + std::vector& yval = sys.y(); // Output time out << t << ","; @@ -161,7 +161,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -179,7 +179,7 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); fileout.close(); diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt b/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt index ef1fc7a9..c5ee2d80 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt @@ -5,4 +5,4 @@ target_link_libraries(TenGenGenrou install(TARGETS TenGenGenrou RUNTIME DESTINATION bin) # Not used for tesing for now. -# add_test(NAME GenrouTest3 COMMAND $) +add_test(NAME TenGenGenrou COMMAND $) diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp index 4f808c69..b3be48eb 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp @@ -118,7 +118,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& yval = sys.y(); + std::vector& yval = sys.y(); // Output time out << t << ","; @@ -147,7 +147,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -165,7 +165,7 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); fileout.close(); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt index ee4ffb76..aecb3340 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt @@ -13,9 +13,9 @@ target_link_libraries(ThreeBusBasicJson install(TARGETS ThreeBusBasicJson RUNTIME DESTINATION ${_install_path}) gridkit_example_add_file(ThreeBusBasic.json) -add_test(NAME GenrouTest2 COMMAND ThreeBusBasic) -add_test(NAME GenrouTest2Json +add_test(NAME ThreeBusBasic COMMAND ThreeBusBasic) +add_test(NAME ThreeBusBasicJson COMMAND ThreeBusBasicJson ${CMAKE_CURRENT_BINARY_DIR}/ThreeBusBasic.json) -add_test(NAME GenrouTest2Json_no_arg +add_test(NAME ThreeBusBasicJson_no_arg COMMAND ThreeBusBasicJson WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp index 55a07472..baad1107 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp @@ -27,11 +27,11 @@ using index_type = size_t; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -43,7 +43,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -73,7 +73,7 @@ int main() { using namespace GridKit::PhasorDynamics; using namespace AnalysisManager::Sundials; - using BusType = BusData::BusType; + using BusType = BusData::BusType; auto error_allowed = static_cast(1e-4); @@ -83,7 +83,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(3); @@ -210,7 +210,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); // Bus 1 -> +2 // Bus 2 -> +2 @@ -219,10 +219,10 @@ int main() // output.push_back(OutputData{t, - 1.0 + y_val[5], // Gen 1 Speed -> 4 + 1 - 1.0 + y_val[24], // Gen 2 Speed -> 23 + 1 - std::sqrt(y_val[0] * y_val[0] + y_val[1] * y_val[1]), // Bus 1 Vmag - std::sqrt(y_val[2] * y_val[2] + y_val[3] * y_val[3])}); // Bus 2 Vmag + 1.0 + static_cast(y_val[5]), // Gen 1 Speed -> 4 + 1 + 1.0 + static_cast(y_val[24]), // Gen 2 Speed -> 23 + 1 + std::sqrt(static_cast(y_val[0]) * static_cast(y_val[0]) + static_cast(y_val[1]) * static_cast(y_val[1])), // Bus 1 Vmag + std::sqrt(static_cast(y_val[2]) * static_cast(y_val[2]) + static_cast(y_val[3]) * static_cast(y_val[3]))}); // Bus 2 Vmag }; // Set up simulation @@ -230,7 +230,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -248,19 +248,19 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); /* Check worst-case error */ real_type worst_error = 0; real_type worst_error_time = 0; - std::ostream nullout(nullptr); - std::ostream& out = nullout; + // std::ostream nullout(nullptr); + // std::ostream& out = nullout; - // // Uncomment code below to print output to a file: - // std::ofstream fileout; - // fileout.open("Example_ThreeBus_Basic_results.csv"); - // std::ostream& out = fileout; + // Uncomment code below to print output to a file: + std::ofstream fileout; + fileout.open("Example_ThreeBus_Basic_results.csv"); + std::ostream& out = fileout; out << "Time,gen2speed,gen3speed,v2mag,v3mag\n"; out << 0. << "," << 1. << "," << 1. << "," << 1. << "," << 1. << "\n"; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp index e01906ff..225fb919 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp @@ -27,11 +27,11 @@ using index_type = size_t; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -43,7 +43,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -134,7 +134,7 @@ int main(int argc, const char* argv[]) auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); // Bus 1 -> +2 // Bus 2 -> +2 @@ -143,10 +143,10 @@ int main(int argc, const char* argv[]) // output.push_back(OutputData{t, - 1.0 + y_val[5], // Gen 1 Speed -> 4 + 1 - 1.0 + y_val[24], // Gen 2 Speed -> 23 + 1 - std::sqrt(y_val[0] * y_val[0] + y_val[1] * y_val[1]), // Bus 1 Vmag - std::sqrt(y_val[2] * y_val[2] + y_val[3] * y_val[3])}); // Bus 2 Vmag + 1.0 + static_cast(y_val[5]), // Gen 1 Speed -> 4 + 1 + 1.0 + static_cast(y_val[24]), // Gen 2 Speed -> 23 + 1 + std::sqrt(static_cast(y_val[0]) * static_cast(y_val[0]) + static_cast(y_val[1]) * static_cast(y_val[1])), // Bus 1 Vmag + std::sqrt(static_cast(y_val[2]) * static_cast(y_val[2]) + static_cast(y_val[3]) * static_cast(y_val[3]))}); // Bus 2 Vmag }; // Set up simulation @@ -154,7 +154,7 @@ int main(int argc, const char* argv[]) ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -172,7 +172,7 @@ int main(int argc, const char* argv[]) ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); sys.stopMonitor(); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp index 5b6d073b..c7660de6 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp @@ -25,15 +25,15 @@ using scalar_type = double; using real_type = double; using index_type = size_t; -using BusType = GridKit::PhasorDynamics::BusData::BusType; +using BusType = GridKit::PhasorDynamics::BusData::BusType; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -45,7 +45,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -84,7 +84,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(3); @@ -103,7 +103,7 @@ int main() // Bus 2 data.bus[2].bus_id = 2; - data.bus[2].bus_type = BusData::BusType::DEFAULT; + data.bus[2].bus_type = BusData::BusType::DEFAULT; data.bus[2].Vr0 = 0.9610827543495831; data.bus[2].Vi0 = -0.13122476630506485; @@ -187,13 +187,13 @@ int main() auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); output.push_back(OutputData{t, - 1 + y_val[5], - 1 + y_val[10], - std::hypot(y_val[0], y_val[1]), - std::hypot(y_val[2], y_val[3])}); + 1 + static_cast(y_val[5]), + 1 + static_cast(y_val[10]), + std::hypot(static_cast(y_val[0]), static_cast(y_val[1])), + std::hypot(static_cast(y_val[2]), static_cast(y_val[3]))}); }; // Set up simulation @@ -201,7 +201,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -219,19 +219,19 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); /* Check worst-case error */ real_type worst_error = 0; real_type worst_error_time = 0; - std::ostream nullout(nullptr); - std::ostream& out = nullout; + // std::ostream nullout(nullptr); + // std::ostream& out = nullout; - // // Uncomment code below to print output to a file: - // std::ofstream fileout; - // fileout.open("Example_ThreeBus_Classical_results.csv"); - // std::ostream& out = fileout; + // Uncomment code below to print output to a file: + std::ofstream fileout; + fileout.open("Example_ThreeBus_Classical_results.csv"); + std::ostream& out = fileout; out << "Time,gen2speed,gen3speed,v2mag,v3mag\n"; out << 0. << "," << 1. << "," << 1. << "," << 1. << "," << 1. << "\n"; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp index 1e69f8f9..98ddc735 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp @@ -27,11 +27,11 @@ using index_type = size_t; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -43,7 +43,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -134,13 +134,13 @@ int main(int argc, const char* argv[]) auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); output.push_back(OutputData{t, - 1 + y_val[5], - 1 + y_val[10], - std::hypot(y_val[0], y_val[1]), - std::hypot(y_val[2], y_val[3])}); + 1 + static_cast(y_val[5]), + 1 + static_cast(y_val[10]), + std::hypot(static_cast(y_val[0]), static_cast(y_val[1])), + std::hypot(static_cast(y_val[2]), static_cast(y_val[3]))}); }; // Set up simulation @@ -148,7 +148,7 @@ int main(int argc, const char* argv[]) ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -166,7 +166,7 @@ int main(int argc, const char* argv[]) ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); /* Check worst-case error */ real_type worst_error = 0; diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp index d7f77998..117e6e06 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp @@ -28,7 +28,7 @@ int main() using real_type = double; using index_type = size_t; - using BusType = BusData::BusType; + using BusType = BusData::BusType; std::cout << "Example: TwoBusBasic\n"; @@ -36,7 +36,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(2); @@ -135,7 +135,10 @@ int main() { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp index 1bdd42ff..a2961c70 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp @@ -115,7 +115,10 @@ int main(int argc, const char* argv[]) { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt index a1b3895d..8e514ab8 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt @@ -15,9 +15,9 @@ target_link_libraries(TwoBusIeeet1Json install(TARGETS TwoBusIeeet1Json RUNTIME DESTINATION ${_install_path}) gridkit_example_add_file(TwoBusIeeet1.json) -add_test(NAME GenrouTest1_Ieeet1 COMMAND TwoBusIeeet1) -add_test(NAME GenrouTest1_Ieeet1_Json +add_test(NAME TwoBusIeeet1 COMMAND TwoBusIeeet1) +add_test(NAME TwoBusIeeet1Json COMMAND TwoBusIeeet1Json ${CMAKE_CURRENT_BINARY_DIR}/TwoBusIeeet1.json) -add_test(NAME GenrouTest1_Ieeet1_Json_no_arg +add_test(NAME TwoBusIeeet1Json_no_arg COMMAND TwoBusIeeet1Json WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp index aa8f478e..886aaae1 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp @@ -27,7 +27,7 @@ int main() using real_type = double; using index_type = size_t; - using BusType = BusData::BusType; + using BusType = BusData::BusType; std::cout << "Example: TwoBusTgov1 + IEEET1 Exciter\n"; @@ -35,7 +35,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(2); @@ -189,11 +189,11 @@ int main() // Exc -> 9 States -> Start Idx 24 output.push_back(OutputData{ t, - y_val[0], // Bus Vr - y_val[1], // Bus Vi - y_val[3], // Gen Speed - y_val[23], // Gov Pmech - y_val[26], // Exc Efd + static_cast(y_val[0]), // Bus Vr + static_cast(y_val[1]), // Bus Vi + static_cast(y_val[3]), // Gen Speed + static_cast(y_val[23]), // Gov Pmech + static_cast(y_val[26]), // Exc Efd }); }; diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp index 2e546f0c..4c17ae62 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp @@ -122,11 +122,11 @@ int main(int argc, const char* argv[]) // Exc -> 9 States -> Start Idx 24 output.push_back(OutputData{ t, - y_val[0], // Bus Vr - y_val[1], // Bus Vi - y_val[3], // Gen Speed - y_val[23], // Gov Pmech - y_val[26], // Exc Efd + static_cast(y_val[0]), // Bus Vr + static_cast(y_val[1]), // Bus Vi + static_cast(y_val[3]), // Gen Speed + static_cast(y_val[23]), // Gov Pmech + static_cast(y_val[26]), // Exc Efd }); }; diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt index c0711d20..c6c8008e 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt @@ -15,9 +15,9 @@ target_link_libraries(TwoBusTgov1Json install(TARGETS TwoBusTgov1Json RUNTIME DESTINATION ${_install_path}) gridkit_example_add_file(TwoBusTgov1.json) -add_test(NAME GenrouTest1_tgov1 COMMAND TwoBusTgov1) -add_test(NAME GenrouTest1_tgov1_json +add_test(NAME TwoBusTgov1 COMMAND TwoBusTgov1) +add_test(NAME TwoBusTgov1Json COMMAND TwoBusTgov1Json ${CMAKE_CURRENT_BINARY_DIR}/TwoBusTgov1.json) -add_test(NAME GenrouTest1_tgov1_json_no_arg +add_test(NAME TwoBusTgov1Json_no_arg COMMAND TwoBusTgov1Json WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp index 4f33d3d3..b4893518 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp @@ -156,7 +156,10 @@ int main() { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp index ad3f7ff7..2d825a71 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp @@ -119,7 +119,10 @@ int main(int argc, const char* argv[]) { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index c6bd98a7..a086c754 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -236,24 +236,22 @@ namespace GridKit RealT Xdp{0.5}; // Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(H, D, Ra, Xdp); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(H, D, Ra, Xdp); // Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(H, D, Ra, Xdp); /// Compare DependencyTracking dependencies to Enzyme's - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i])); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i])); } return success.report(__func__); } private: - std::vector DependencyTrackingJacobian( + std::vector DependencyTrackingJacobian( const RealT H, const RealT D, const RealT Ra, const RealT Xdp) { DependencyTracking::Variable Vr1{1.0}; ///< Bus-1 real voltage @@ -261,29 +259,90 @@ namespace GridKit PhasorDynamics::Bus bus(Vr1, Vi1); PhasorDynamics::GenClassical gen(&bus, 1, 1.0, 1.0, H, D, Ra, Xdp); + bus.allocate(); - bus.initialize(); gen.allocate(); + + // Get d/dy + bus.initialize(); + gen.initialize(); + + for (size_t i = 0; i < gen.size(); ++i) + { + gen.y()[i].setVariableNumber(i); ///< Generator independent variables + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + gen.size()); // Bus independent variables + } + + bus.evaluateResidual(); + gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking + ///< the dependencies + std::vector residual_y = gen.getResidual(); + + // Get d/dy' + bus.initialize(); gen.initialize(); for (size_t i = 0; i < gen.size(); ++i) { - gen.y()[i].setVariableNumber(i); ///< Independent variables + gen.yp()[i].setVariableNumber(i); ///< Generator independent variables } + bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual = gen.getResidual(); + std::vector residual_yp = gen.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residual.size(); ++i) + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) { - std::cout << i << "th residual: "; - (residual[i]).print(std::cout); + std::cout << i << "th residual, y: "; + (residual_y[i]).print(std::cout); + std::cout << "\n"; + std::cout << i << "th residual, yp: "; + (residual_yp[i]).print(std::cout); std::cout << "\n"; } - return residual; + // Extract the dependencies and add d/dy' to d/dy + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); + DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + auto value_yp = it_yp->second; + dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); + } + else + { + dependencies[i].insert(std::make_pair(index_y, value_y)); + } + } + + // Insert yp dependencies that did not exist in the y dependencies + for (const auto& pair_yp : dependency_yp) + { + auto index_yp = pair_yp.first; + auto value_yp = pair_yp.second; + auto it_y = dependency_y.find(index_yp); + if (it_y == dependency_y.end()) + { + dependencies[i].insert(std::make_pair(index_yp, value_yp)); + } + } + } + + return dependencies; } std::vector EnzymeJacobian( @@ -294,12 +353,24 @@ namespace GridKit PhasorDynamics::Bus bus(Vr1, Vi1); PhasorDynamics::GenClassical gen(&bus, 1, 1.0, 1.0, H, D, Ra, Xdp); + bus.allocate(); - bus.initialize(); gen.allocate(); + + bus.initialize(); gen.initialize(); + gen.updateTime(0.0, 1.0); + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + gen.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + gen.size()); // Reset bus residual indices + } + + bus.evaluateResidual(); gen.evaluateResidual(); + gen.evaluateJacobian(); GridKit::LinearAlgebra::COO_Matrix model_jacobian = gen.getJacobian(); model_jacobian.printMatrix("Model Jacobian"); diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 90a160fa..081bd5f0 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -231,24 +231,22 @@ namespace GridKit auto tol = 10 * std::numeric_limits::epsilon(); // Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(); // Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(); /// Compare DependencyTracking dependencies to Enzyme's - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i], tol)); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol)); } return success.report(__func__); } private: - std::vector DependencyTrackingJacobian() + std::vector DependencyTrackingJacobian() { DependencyTracking::Variable Vr1{1.0}; ///< Bus real voltage DependencyTracking::Variable Vi1{0.0}; ///< Bus imaginary voltage @@ -275,30 +273,88 @@ namespace GridKit 0); bus.allocate(); + gen.allocate(); + + // Get d/dy bus.initialize(); + gen.initialize(); - gen.allocate(); + for (size_t i = 0; i < gen.size(); ++i) + { + gen.y()[i].setVariableNumber(i); ///< Generator independent variables + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + gen.size()); // Bus independent variables + } + + bus.evaluateResidual(); + gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking + ///< the dependencies + std::vector residual_y = gen.getResidual(); + + // Get d/dy' + bus.initialize(); gen.initialize(); for (size_t i = 0; i < gen.size(); ++i) { - gen.y()[i].setVariableNumber(i); ///< Independent variables + gen.yp()[i].setVariableNumber(i); ///< Generator independent variables } bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual = gen.getResidual(); + std::vector residual_yp = gen.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residual.size(); ++i) + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) { - std::cout << i << "th residual: "; - (residual[i]).print(std::cout); + std::cout << i << "th residual, y: "; + (residual_y[i]).print(std::cout); + std::cout << "\n"; + std::cout << i << "th residual, yp: "; + (residual_yp[i]).print(std::cout); std::cout << "\n"; } - return residual; + // Extract the dependencies and add d/dy' to d/dy + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); + DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + auto value_yp = it_yp->second; + dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); + } + else + { + dependencies[i].insert(std::make_pair(index_y, value_y)); + } + } + + // Insert yp dependencies that did not exist in the y dependencies + for (const auto& pair_yp : dependency_yp) + { + auto index_yp = pair_yp.first; + auto value_yp = pair_yp.second; + auto it_y = dependency_y.find(index_yp); + if (it_y == dependency_y.end()) + { + dependencies[i].insert(std::make_pair(index_yp, value_yp)); + } + } + } + + return dependencies; } std::vector EnzymeJacobian() @@ -328,11 +384,20 @@ namespace GridKit 0); bus.allocate(); - bus.initialize(); - bus.evaluateResidual(); - gen.allocate(); + + bus.initialize(); gen.initialize(); + + gen.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + gen.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + gen.size()); // Reset bus residual indices + } + + bus.evaluateResidual(); gen.evaluateResidual(); gen.evaluateJacobian(); diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 99d0e84c..9054bbda 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -272,18 +272,16 @@ namespace GridKit gendata.parameters[GenrouParameters::S12] = 0.; /// Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(busdata, gendata); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(busdata, gendata); /// Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(busdata, gendata); /// Compare DependencyTracking dependencies to Enzyme's auto tol = 10 * std::numeric_limits::epsilon(); - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i], tol)); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol)); } success.expectFailure(); // TODO: Resolve mismatch from sigmoid approximation @@ -292,7 +290,7 @@ namespace GridKit } private: - std::vector DependencyTrackingJacobian( + std::vector DependencyTrackingJacobian( PhasorDynamics::BusData busdata, PhasorDynamics::GenrouData gendata) { @@ -307,6 +305,7 @@ namespace GridKit gov.allocate(); gen.allocate(); + // Get d/dy bus.initialize(); gen.initialize(); gov.initialize(); @@ -321,17 +320,72 @@ namespace GridKit gen.evaluateResidual(); gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking // the dependencies - std::vector residual = gov.getResidual(); + std::vector residual_y = gov.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residual.size(); ++i) + // Get d/dy' + bus.initialize(); + gen.initialize(); + gov.initialize(); + + for (size_t i = 0; i < gov.size(); ++i) + { + gov.yp()[i].setVariableNumber(i); ///< Governor independent variables + } + + bus.evaluateResidual(); + gen.evaluateResidual(); + gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking + // the dependencies + std::vector residual_yp = gov.getResidual(); + + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) { - std::cout << i << "th residual: "; - (residual[i]).print(std::cout); + std::cout << i << "th residual, y: "; + (residual_y[i]).print(std::cout); + std::cout << "\n"; + std::cout << i << "th residual, yp: "; + (residual_yp[i]).print(std::cout); std::cout << "\n"; } - return residual; + // Extract the dependencies and add d/dy' to d/dy + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); + DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + auto value_yp = it_yp->second; + dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); + } + else + { + dependencies[i].insert(std::make_pair(index_y, value_y)); + } + } + + // Insert yp dependencies that did not exist in the y dependencies + for (const auto& pair_yp : dependency_yp) + { + auto index_yp = pair_yp.first; + auto value_yp = pair_yp.second; + auto it_y = dependency_y.find(index_yp); + if (it_y == dependency_y.end()) + { + dependencies[i].insert(std::make_pair(index_yp, value_yp)); + } + } + } + + return dependencies; } std::vector EnzymeJacobian( @@ -355,6 +409,8 @@ namespace GridKit gen.initialize(); gov.initialize(); + gov.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term + bus.evaluateResidual(); gen.evaluateResidual(); gov.evaluateResidual(); diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 727d2656..1993c94a 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -204,24 +204,22 @@ namespace GridKit data.branch[0].parameters[BranchParameters::B] = 1.2; // Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(data); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(data); // Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(data); /// Compare DependencyTracking dependencies to Enzyme's - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i])); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i])); } return success.report(__func__); } private: - std::vector DependencyTrackingJacobian( + std::vector DependencyTrackingJacobian( PhasorDynamics::SystemModelData data) { // Create an empty system model @@ -239,17 +237,24 @@ namespace GridKit // Evaluate and get the system residuals system.evaluateResidual(); - std::vector residuals = system.getResidual(); + std::vector residual = system.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residuals.size(); ++i) + // Print the dependencies + for (size_t i = 0; i < residual.size(); ++i) { std::cout << i << "th residual: "; - (residuals[i]).print(std::cout); + (residual[i]).print(std::cout); std::cout << "\n"; } - return residuals; + // Extract the dependencies + std::vector dependencies(residual.size()); + for (IdxT i = 0; i < residual.size(); ++i) + { + dependencies[i] = (residual[i]).getDependencies(); + } + + return dependencies; } std::vector EnzymeJacobian(