diff --git a/.github/workflows/pre_commit.yaml b/.github/workflows/pre_commit.yaml index a6ecfed8a..a56c74938 100644 --- a/.github/workflows/pre_commit.yaml +++ b/.github/workflows/pre_commit.yaml @@ -19,4 +19,4 @@ jobs: fetch: false committer_name: GitHub Actions committer_email: actions@github.com - message: Apply pre-commmit fixes + message: Apply pre-commit fixes diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp index 88b62a0d1..a5293ddf9 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp @@ -9,8 +9,6 @@ #include #include #include -#include -#include namespace GridKit { @@ -23,14 +21,13 @@ namespace GridKit * * @tparam ModelT - model type * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type */ - template + template struct DfDwb { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + using ScalarT = typename ModelT::ScalarT; + using IdxT = typename ModelT::IdxT; + using RealT = typename ModelT::RealT; /** * @param[in] model - Pointer to the model to be differentiated @@ -41,7 +38,7 @@ namespace GridKit * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian + * @param[out] nnz - Number of nonzeros */ static void eval(ModelT* model, size_t n_res, @@ -54,12 +51,11 @@ namespace GridKit IdxT* rows, IdxT* cols, RealT* vals, - MatrixT& jac) + IdxT& nnz) { if (n_res > 0 && n_var > 0) { std::vector elementary_v(n_var); - IdxT nnz = 0; for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage. @see LowerSparseStorage.hpp @@ -69,6 +65,7 @@ namespace GridKit ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, + 1.0, // value scaling res_indices, var_indices, rows, @@ -80,8 +77,8 @@ namespace GridKit 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, + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, enzyme_const, @@ -95,9 +92,6 @@ namespace GridKit elementary_v.data(), d_output); } - - // Store result - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes } } @@ -111,7 +105,7 @@ namespace GridKit * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in] ws - Signal variables - * @param[in,out] jac - Jacobian + * @param[out] nnz - Number of nonzeros */ static void eval(ModelT* model, size_t n_res, @@ -125,12 +119,11 @@ namespace GridKit IdxT* rows, IdxT* cols, RealT* vals, - MatrixT& jac) + IdxT& nnz) { if (n_res > 0 && n_var > 0) { std::vector elementary_v(n_var); - IdxT nnz = 0; for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage. @see LowerSparseStorage.hpp @@ -140,6 +133,7 @@ namespace GridKit ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, + 1.0, // value scaling res_indices, var_indices, rows, @@ -151,8 +145,8 @@ namespace GridKit 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, + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, enzyme_const, @@ -168,9 +162,6 @@ namespace GridKit elementary_v.data(), d_output); } - - // Store result - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes } } }; diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp index 7124ad0b0..ad3e3e568 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp @@ -9,8 +9,6 @@ #include #include #include -#include -#include namespace GridKit { @@ -23,14 +21,13 @@ namespace GridKit * * @tparam ModelT - model type * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type */ - template + template struct DfDws { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + using ScalarT = typename ModelT::ScalarT; + using IdxT = typename ModelT::IdxT; + using RealT = typename ModelT::RealT; /** * @param[in] model - Pointer to the model to be differentiated @@ -42,7 +39,7 @@ namespace GridKit * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in] ws - Signal variables - * @param[in,out] jac - Jacobian + * @param[out] nnz - Number of nonzeros */ static void eval(ModelT* model, size_t n_res, @@ -56,12 +53,11 @@ namespace GridKit IdxT* rows, IdxT* cols, RealT* vals, - MatrixT& jac) + IdxT& nnz) { if (n_res > 0 && n_var > 0) { std::vector elementary_v(n_var); - IdxT nnz = 0; for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage. @see LowerSparseStorage.hpp @@ -71,6 +67,7 @@ namespace GridKit ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, + 1.0, // value scaling res_indices, var_indices, rows, @@ -82,8 +79,8 @@ namespace GridKit 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, + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, enzyme_const, @@ -99,9 +96,6 @@ namespace GridKit elementary_v.data(), d_output); } - - // Store result - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes } } }; diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp index 51142932c..4b1b8f87b 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp @@ -9,8 +9,6 @@ #include #include #include -#include -#include namespace GridKit { @@ -23,14 +21,13 @@ namespace GridKit * * @tparam ModelT - model type * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type */ - template + template struct DfDy { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + using ScalarT = typename ModelT::ScalarT; + using IdxT = typename ModelT::IdxT; + using RealT = typename ModelT::RealT; /** * @param[in] model - Pointer to the model to be differentiated @@ -41,7 +38,7 @@ namespace GridKit * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian + * @param[out] nnz - Number of nonzeros */ static void eval(ModelT* model, size_t n_res, @@ -54,13 +51,11 @@ namespace GridKit IdxT* rows, IdxT* cols, RealT* vals, - MatrixT& jac) + IdxT& nnz) { if (n_res > 0 && n_var > 0) { - // df/dy std::vector elementary_v(n_var); - IdxT nnz = 0; for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage. @see LowerSparseStorage.hpp @@ -70,6 +65,7 @@ namespace GridKit ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, + 1.0, // value scaling res_indices, var_indices, rows, @@ -81,8 +77,8 @@ namespace GridKit 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, + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, enzyme_dup, @@ -96,125 +92,6 @@ namespace GridKit elementary_v.data(), d_output); } - - // Store result - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes - - // There is no df/dy' when alpha is not passed as an argument - // @todo: Implement a generic way to identify these cases at compile time - } - } - - /** - * @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 IdxT* res_indices, - const IdxT* var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - RealT alpha, - IdxT* rows, - IdxT* cols, - RealT* vals, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - // df/dy - std::vector elementary_v(n_var); - IdxT nnz = 0; - 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, - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic 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 - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes - - // df/dy' - nnz = 0; - 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, - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic 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 - jac.setValues(alpha, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes } } @@ -228,8 +105,7 @@ namespace GridKit * @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 + * @param[out] nnz - Number of nonzeros */ static void eval(ModelT* model, size_t n_res, @@ -240,17 +116,14 @@ namespace GridKit ScalarT* yp, ScalarT* wb, ScalarT* ws, - RealT alpha, IdxT* rows, IdxT* cols, RealT* vals, - MatrixT& jac) + IdxT& nnz) { if (n_res > 0 && n_var > 0) { - // df/dy std::vector elementary_v(n_var); - IdxT nnz = 0; for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage. @see LowerSparseStorage.hpp @@ -260,6 +133,7 @@ namespace GridKit ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, + 1.0, // value scaling res_indices, var_indices, rows, @@ -271,8 +145,8 @@ namespace GridKit 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, + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, enzyme_dup, @@ -288,52 +162,6 @@ namespace GridKit elementary_v.data(), d_output); } - - // Store result - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes - - // df/dy' - nnz = 0; - 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, - res_indices, - var_indices, - rows, - cols, - vals, - &nnz); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Core automatic 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 - jac.setValues(alpha, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes } } }; diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDyp.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDyp.hpp new file mode 100644 index 000000000..83baff9a6 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDyp.hpp @@ -0,0 +1,174 @@ +/** + * @file DfDyp.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: Internal Jacobian, alpha*df/dyp + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + */ + template + struct DfDyp + { + using ScalarT = typename ModelT::ScalarT; + using IdxT = typename ModelT::IdxT; + using RealT = typename ModelT::RealT; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Global residual indices + * @param[in] var_indices - Global variable indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in] alpha - Time derivative jacobian coefficient + * @param[out] nnz - Number of nonzeros + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const IdxT* res_indices, + const IdxT* var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + RealT alpha, + IdxT* rows, + IdxT* cols, + RealT* vals, + IdxT& nnz) + { + if (n_res > 0 && n_var > 0) + { + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage. @see LowerSparseStorage.hpp + ScalarT* output = __enzyme_todense((void*) ident_load, + (void*) ident_store, + var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, + (void*) sparse_store, + var_i, + alpha, // value scaling + res_indices, + var_indices, + rows, + cols, + vals, + &nnz); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_dup, + yp, + output, + enzyme_const, + wb, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + } + } + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in] ws - Signal variables + * @param[in] alpha - Time derivative jacobian coefficient + * @param[out] nnz - Number of nonzeros + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const IdxT* res_indices, + const IdxT* var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + RealT alpha, + IdxT* rows, + IdxT* cols, + RealT* vals, + IdxT& nnz) + { + if (n_res > 0 && n_var > 0) + { + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage. @see LowerSparseStorage.hpp + ScalarT* output = __enzyme_todense((void*) ident_load, + (void*) ident_store, + var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, + (void*) sparse_store, + var_i, + alpha, // value scaling + res_indices, + var_indices, + rows, + cols, + vals, + &nnz); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_dup, + yp, + output, + enzyme_const, + wb, + enzyme_const, + ws, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp index 11ea1dcdc..7bad3d0c8 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp @@ -9,8 +9,6 @@ #include #include #include -#include -#include namespace GridKit { @@ -23,14 +21,13 @@ namespace GridKit * * @tparam ModelT - model type * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type */ - template + template struct DhDwb { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + using ScalarT = typename ModelT::ScalarT; + using IdxT = typename ModelT::IdxT; + using RealT = typename ModelT::RealT; /** * @param[in] model - Pointer to the model to be differentiated @@ -41,8 +38,7 @@ namespace GridKit * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian - * @param[in] flag - To append values or deduplicate right away (bus-owned values) + * @param[out] nnz - Number of nonzeros */ static void eval(ModelT* model, size_t n_res, @@ -55,13 +51,11 @@ namespace GridKit IdxT* rows, IdxT* cols, RealT* vals, - MatrixT& jac, - const bool flag = false) + IdxT& nnz) { if (n_res > 0 && n_var > 0) { std::vector elementary_v(n_var); - IdxT nnz = 0; for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage. @see LowerSparseStorage.hpp @@ -71,6 +65,7 @@ namespace GridKit ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, + 1.0, // value scaling res_indices, var_indices, rows, @@ -82,8 +77,8 @@ namespace GridKit 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, + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, enzyme_const, @@ -97,16 +92,6 @@ namespace GridKit elementary_v.data(), d_output); } - - // Store result - if (flag) - { - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes - } - else - { - jac.axpy(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes - } } } }; diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp index cddcf93a3..0bc62301f 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp @@ -9,8 +9,6 @@ #include #include #include -#include -#include namespace GridKit { @@ -23,14 +21,13 @@ namespace GridKit * * @tparam ModelT - model type * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type */ - template + template struct DhDy { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + using ScalarT = typename ModelT::ScalarT; + using IdxT = typename ModelT::IdxT; + using RealT = typename ModelT::RealT; /** * @param[in] model - Pointer to the model to be differentiated @@ -41,7 +38,7 @@ namespace GridKit * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian + * @param[out] nnz - Number of nonzeros */ static void eval(ModelT* model, size_t n_res, @@ -54,12 +51,11 @@ namespace GridKit IdxT* rows, IdxT* cols, RealT* vals, - MatrixT& jac) + IdxT& nnz) { if (n_res > 0 && n_var > 0) { std::vector elementary_v(n_var); - IdxT nnz = 0; for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage. @see LowerSparseStorage.hpp @@ -69,6 +65,7 @@ namespace GridKit ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, + 1.0, // value scaling res_indices, var_indices, rows, @@ -80,8 +77,8 @@ namespace GridKit 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, + // Core automatic differentiation intrinsic that will be replaced by a derivative + __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, enzyme_dup, @@ -95,9 +92,6 @@ namespace GridKit elementary_v.data(), d_output); } - - // Store result - jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes } } }; diff --git a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp index 4f9cfabe7..6714a28a7 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp @@ -38,6 +38,7 @@ namespace GridKit * @param[in] row - row to be stored * @param[in] col - column to be stored * @param[in] val - value to be stored + * @param[in] scaling - scaling factor for values * @param[in] res_indices - Global residual indices * @param[in] var_indices - Global variable indices * @param[in,out] rows - buffer where row will be stored @@ -49,6 +50,7 @@ namespace GridKit size_t row, size_t col, float val, + float scaling, const size_t* row_indices, const size_t* col_indices, size_t* rows, @@ -62,7 +64,7 @@ namespace GridKit { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; - vals[static_cast(nnz)] = val; + vals[static_cast(nnz)] = scaling * val; nnz++; } } @@ -75,6 +77,7 @@ namespace GridKit * @param[in] row - row to be stored * @param[in] col - column to be stored * @param[in] val - value to be stored + * @param[in] scaling - scaling factor for values * @param[in] res_indices - Global residual indices * @param[in] var_indices - Global variable indices * @param[in,out] rows - buffer where row will be stored @@ -86,6 +89,7 @@ namespace GridKit long int row, long int col, float val, + float scaling, const long int* row_indices, const long int* col_indices, long int* rows, @@ -99,7 +103,7 @@ namespace GridKit { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; - vals[static_cast(nnz)] = val; + vals[static_cast(nnz)] = scaling * val; nnz++; } } @@ -112,6 +116,7 @@ namespace GridKit * @param[in] row - row to be stored * @param[in] col - column to be stored * @param[in] val - value to be stored + * @param[in] scaling - scaling factor for values * @param[in] res_indices - Global residual indices * @param[in] var_indices - Global variable indices * @param[in,out] rows - buffer where row will be stored @@ -123,6 +128,7 @@ namespace GridKit size_t row, size_t col, double val, + double scaling, const size_t* row_indices, const size_t* col_indices, size_t* rows, @@ -136,7 +142,7 @@ namespace GridKit { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; - vals[static_cast(nnz)] = val; + vals[static_cast(nnz)] = scaling * val; nnz++; } } @@ -149,6 +155,7 @@ namespace GridKit * @param[in] row - row to be stored * @param[in] col - column to be stored * @param[in] val - value to be stored + * @param[in] scaling - scaling factor for values * @param[in] res_indices - Global residual indices * @param[in] var_indices - Global variable indices * @param[in,out] rows - buffer where row will be stored @@ -160,6 +167,7 @@ namespace GridKit long int row, long int col, double val, + double scaling, const long int* row_indices, const long int* col_indices, long int* rows, @@ -173,7 +181,7 @@ namespace GridKit { rows[static_cast(nnz)] = row_mapped; cols[static_cast(nnz)] = col_mapped; - vals[static_cast(nnz)] = val; + vals[static_cast(nnz)] = scaling * val; nnz++; } } @@ -189,6 +197,7 @@ namespace GridKit * @param[in] val - value to be stored * @param[in] row - row to be stored * @param[in] col - column to be stored + * @param[in] scaling - scaling factor for values * @param[in] res_indices - Global residual indices * @param[in] var_indices - Global variable indices * @param[in,out] rows - buffer where row will be stored @@ -201,6 +210,7 @@ namespace GridKit ScalarT val, IdxT row, IdxT col, + ScalarT scaling, const IdxT* row_indices, const IdxT* col_indices, IdxT* rows, @@ -217,16 +227,16 @@ namespace GridKit if constexpr (std::is_same::value) { if constexpr (std::is_same::value) - inner_store_float_size_t(row, col, val, row_indices, col_indices, rows, cols, vals, nnz); + inner_store_float_size_t(row, col, val, scaling, row_indices, col_indices, rows, cols, vals, nnz); else - inner_store_double_size_t(row, col, val, row_indices, col_indices, rows, cols, vals, nnz); + inner_store_double_size_t(row, col, val, scaling, row_indices, col_indices, rows, cols, vals, nnz); } else if constexpr (std::is_same::value) { if constexpr (std::is_same::value) - inner_store_float_long_int(row, col, val, row_indices, col_indices, rows, cols, vals, nnz); + inner_store_float_long_int(row, col, val, scaling, row_indices, col_indices, rows, cols, vals, nnz); else - inner_store_double_long_int(row, col, val, row_indices, col_indices, rows, cols, vals, nnz); + inner_store_double_long_int(row, col, val, scaling, row_indices, col_indices, rows, cols, vals, nnz); } else { diff --git a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp index 09d67f08d..efbb0eaa6 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp @@ -34,10 +34,9 @@ namespace GridKit * * @tparam ModelT - model type * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type * */ - template + template struct ModelWrapper { }; @@ -46,9 +45,11 @@ namespace GridKit * @brief Residual wrapper partial template specialization for InternalResidual * */ - template - struct ModelWrapper + template + struct ModelWrapper { + using ScalarT = typename ModelT::ScalarT; + /** * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables @@ -66,9 +67,11 @@ namespace GridKit * @brief Residual wrapper partial template specialization for InternalResidualWithSignal * */ - template - struct ModelWrapper + template + struct ModelWrapper { + using ScalarT = typename ModelT::ScalarT; + /** * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables @@ -87,9 +90,11 @@ namespace GridKit * @brief Residual wrapper partial template specialization for BusResidual * */ - template - struct ModelWrapper + template + struct ModelWrapper { + using ScalarT = typename ModelT::ScalarT; + /** * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables @@ -107,9 +112,11 @@ namespace GridKit * @brief Residual wrapper partial template specialization for BusResidual11 (branch member function) * */ - template - struct ModelWrapper + template + struct ModelWrapper { + using ScalarT = typename ModelT::ScalarT; + /** * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables @@ -127,9 +134,11 @@ namespace GridKit * @brief Residual wrapper partial template specialization for BusResidual12 (branch member function) * */ - template - struct ModelWrapper + template + struct ModelWrapper { + using ScalarT = typename ModelT::ScalarT; + /** * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables @@ -147,12 +156,14 @@ namespace GridKit * @brief Residual wrapper partial template specialization for BusResidual21 (branch member function) * */ - template - struct ModelWrapper + template + struct ModelWrapper { + using ScalarT = typename ModelT::ScalarT; + /** * @param[in] model - Pointer to the model to be differentiated - __enzyme_fwddiff((void*) ModelWrapper::eval, + __enzyme_fwddiff((void*) ModelWrapper::eval, * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables @@ -168,9 +179,11 @@ namespace GridKit * @brief Residual wrapper partial template specialization for BusResidual22 (branch member function) * */ - template - struct ModelWrapper + template + struct ModelWrapper { + using ScalarT = typename ModelT::ScalarT; + /** * @param[in] model - Pointer to the model to be differentiated * @param[in] y - Internal variables diff --git a/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp b/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp index 9a9df2777..6be78e42f 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp @@ -9,5 +9,6 @@ #include #include #include +#include #include #include diff --git a/GridKit/LinearAlgebra/DenseMatrix/DenseMatrix.hpp b/GridKit/LinearAlgebra/DenseMatrix/DenseMatrix.hpp index 364fab55e..7a255e7de 100644 --- a/GridKit/LinearAlgebra/DenseMatrix/DenseMatrix.hpp +++ b/GridKit/LinearAlgebra/DenseMatrix/DenseMatrix.hpp @@ -6,8 +6,6 @@ #include #include -#include - /** * @brief Class to provide dense matrices. * @@ -22,12 +20,10 @@ namespace GridKit class DenseMatrix { private: - IdxT rows_size_; - IdxT columns_size_; - std::vector values_; - COO_Matrix values_COO_; - bool values_changed_ = false; - bool sparsified_ = false; + IdxT rows_size_; + IdxT columns_size_; + std::vector values_; + bool values_changed_ = false; public: // Constructors and destructors @@ -35,16 +31,13 @@ namespace GridKit ~DenseMatrix(); // Getters and setters - RealT getValue(const IdxT i, const IdxT j) const; - void setValue(const IdxT i, const IdxT j, const RealT value); - void setValues(COO_Matrix values_COO); - void setValues(size_t nnz, const IdxT* rows_coo, const IdxT* cols_coo, const RealT* vals_coo); - std::vector* getValues(); - COO_Matrix* getValuesCOO(); + RealT getValue(const IdxT i, const IdxT j) const; + void setValue(const IdxT i, const IdxT j, const RealT value); + void setValues(size_t nnz, const IdxT* rows_coo, const IdxT* cols_coo, const RealT* vals_coo); + std::vector* getValues(); // Utilities - void toCOO(); - void printMatrix(std::string name = ""); + void print(std::string name = ""); // Purposefully not defining BLAS operations. This class should not be used // for compute. @@ -53,7 +46,7 @@ namespace GridKit /** * @brief DenseMatrix constructor * - * @tparam RealT - Real type for Jacobian entries + * @tparam RealT - Real type for matrix values * @tparam IdxT - Integer data type for matrix indices * * @param[in] IdxT - rows_size @@ -63,15 +56,14 @@ namespace GridKit DenseMatrix::DenseMatrix(const IdxT rows_size, const IdxT columns_size) : rows_size_(rows_size), columns_size_(columns_size), - values_(rows_size * columns_size, 0), - values_COO_(rows_size, columns_size) + values_(rows_size * columns_size, 0) { } /** * @brief DenseMatrix single value getter * - * @tparam RealT - Real type for Jacobian entries + * @tparam RealT - Real type for matrix values * @tparam IdxT - Integer data type for matrix indices * * @param[in] IdxT - i row index @@ -81,15 +73,15 @@ namespace GridKit template inline RealT DenseMatrix::getValue(const IdxT i, const IdxT j) const { - assert(i < this->columns_size_); - assert(j < this->rows_size_); - return this->values_[j * rows_size_ + i]; + assert(i < columns_size_); + assert(j < rows_size_); + return values_[j * rows_size_ + i]; } /** * @brief DenseMatrix single value setter * - * @tparam RealT - Real type for Jacobian entries + * @tparam RealT - Real type for matrix values * @tparam IdxT - Integer data type for matrix indices * * @param[in] IdxT - i row index @@ -99,29 +91,10 @@ namespace GridKit template inline void DenseMatrix::setValue(const IdxT i, const IdxT j, const RealT value) { - assert(i < this->columns_size_); - assert(j < this->rows_size_); - this->values_[j * rows_size_ + i] = value; - values_changed_ = true; - } - - /** - * @brief DenseMatrix value setter from COO - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] COO_Matrix - values_COO - */ - template - inline void DenseMatrix::setValues(COO_Matrix values_COO) - { - std::tuple&, std::vector&, std::vector&> entries = values_COO.getEntries(); - const auto [rcord, ccord, vals] = entries; - for (IdxT idx = 0; idx < values_COO.nnz(); ++idx) - { - this->setValue(rcord[idx], ccord[idx], vals[idx]); - } + assert(i < columns_size_); + assert(j < rows_size_); + values_[j * rows_size_ + i] = value; + values_changed_ = true; } /** @@ -137,14 +110,14 @@ namespace GridKit { for (size_t idx = 0; idx < nnz; ++idx) { - this->setValue(rows_coo[idx], cols_coo[idx], vals_coo[idx]); + setValue(rows_coo[idx], cols_coo[idx], vals_coo[idx]); } } /** * @brief DenseMatrix getter for all values stored as a vector * - * @tparam RealT - Real type for Jacobian entries + * @tparam RealT - Real type for matrix values * @tparam IdxT - Integer data type for matrix indices * * @return Address of the vector containing matrix values @@ -152,79 +125,26 @@ namespace GridKit template inline std::vector* DenseMatrix::getValues() { - return &(this->values_); - } - - /** - * @brief DenseMatrix getter for all values stored as a COO sparse matrix - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @return Address of the COO matrix containing the sparsified matrix values - */ - template - inline COO_Matrix* DenseMatrix::getValuesCOO() - { - if (!sparsified_ || values_changed_) - { - this->toCOO(); - } - return &(this->values_COO_); - } - - /** - * @brief Dense matrix conversion to COO form - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - */ - template - inline void DenseMatrix::toCOO() - { - if (!sparsified_ || values_changed_) - { - IdxT nnz = 0; - std::vector rcord; - std::vector ccord; - std::vector vals; - for (IdxT j = 0; j < this->columns_size_; ++j) - { - for (IdxT i = 0; i < this->rows_size_; ++i) - { - RealT value = this->values_[j * rows_size_ + i]; - if (std::abs(value) > std::numeric_limits::epsilon()) - { - nnz++; - rcord.push_back(i); - ccord.push_back(j); - vals.push_back(value); - } - } - } - values_COO_.setValues(rcord, ccord, vals); - sparsified_ = true; - values_changed_ = false; - } + return &(values_); } /** * @brief Print matrix * - * @tparam RealT - Real type for Jacobian entries + * @tparam RealT - Real type for matrix values * @tparam IdxT - Integer data type for matrix indices * * @param[in] name to identify the specific matrix printed */ template - inline void DenseMatrix::printMatrix(std::string name) + inline void DenseMatrix::print(std::string name) { std::cout << "Dense matrix: " << name << "\n"; - for (IdxT i = 0; i < this->rows_size_; ++i) + for (IdxT i = 0; i < rows_size_; ++i) { - for (IdxT j = 0; j < this->columns_size_; ++j) + for (IdxT j = 0; j < columns_size_; ++j) { - std::cout << this->values_[j * rows_size_ + i] << " "; + std::cout << values_[j * rows_size_ + i] << " "; } std::cout << "\n"; } @@ -233,7 +153,7 @@ namespace GridKit /** * @brief DenseMatrix destructor * - * @tparam RealT - Real type for Jacobian entries + * @tparam RealT - Real type for matrix values * @tparam IdxT - Integer data type for matrix indices */ template diff --git a/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt b/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt index 443baa6fb..c7f5fb4f5 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt +++ b/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt @@ -1,4 +1,4 @@ gridkit_add_library( sparse_matrix SOURCES CooMatrix.cpp CsrMatrix.cpp - HEADERS COO_Matrix.hpp CooMatrix.hpp CsrMatrix.hpp) + HEADERS CooMatrix.hpp CsrMatrix.hpp) diff --git a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp deleted file mode 100644 index b699dedcb..000000000 --- a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp +++ /dev/null @@ -1,1105 +0,0 @@ -#pragma once - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace GridKit -{ - namespace LinearAlgebra - { - /** - * @brief Quick class to provide sparse matrices of COO type. - * - * Simplifies data movement - * - * @todo add functionality to keep track of multiple sorted lists. Faster - * adding of new entries and will have a threshold to sort completely. - * - * m x n sparse matrix - */ - template - class COO_Matrix - { - private: - std::vector values_; - std::vector row_indices_; - std::vector column_indices_; - IdxT rows_size_; - IdxT columns_size_; - bool sorted_; - - public: - // Constructors - COO_Matrix(std::vector r, std::vector c, std::vector v, IdxT m, IdxT n); - COO_Matrix(IdxT m, IdxT n); - COO_Matrix(); - ~COO_Matrix(); - - // Operations - - // --- Functions which call sort --- - std::tuple, std::vector> getRowCopy(IdxT r); - std::tuple&, std::vector&, std::vector&> getEntries(const bool sort = true); - std::tuple, std::vector, std::vector> getEntryCopies(); - std::tuple, std::vector, std::vector> getEntryCopiesSubMatrix(std::vector submap); - - std::tuple, std::vector, std::vector> setDataToCSR(); - std::vector getCSRRowData(); - - // Set values from vector storage. Will sort before storing - void setValues(std::vector r, std::vector c, std::vector v); - void setValues(RealT alpha, IdxT* r, IdxT* c, RealT* v, IdxT nnz); - - // BLAS. Will sort before running - void axpy(RealT alpha, COO_Matrix& a, const bool sort = true); - void axpy(RealT alpha, std::vector r, std::vector c, std::vector v, const bool sort = true); - void axpy(RealT alpha, IdxT* r, IdxT* c, RealT* v, IdxT nnz); - void scal(RealT alpha); - RealT frobNorm(); - - // --- Permutation Operations --- - // Sorting is only done if not already sorted. - void permutation(std::vector row_perm, std::vector col_perm); - void permutationSizeMap(std::vector row_perm, std::vector col_perm, IdxT m, IdxT n); - - // Special matrices (zero and identity) - void zeroMatrix(); // Actually null matrix - void zeroValuedMatrix(); - void identityMatrix(IdxT n); - - // Resort values_ - void sortSparse(); - bool isSorted(); - IdxT nnz(); - void deduplicate(); - - std::tuple getDimensions(); - - void printMatrix(std::string name = ""); - - void printMatrixMarket(const std::string& filename, const std::string& comment); - - static void sortSparseCOO(std::vector& rows, std::vector& columns, std::vector& values, size_t nnz); - - private: - IdxT indexStartRow(const std::vector& rows, IdxT r); - IdxT sparseCordBinarySearch(const std::vector& rows, const std::vector& columns, IdxT ri, IdxT ci); - bool checkIncreaseSize(IdxT r, IdxT c); - }; - - /** - * @brief Get copy of row index - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] r row index - * @return std::tuple, std::vector> - */ - template - inline std::tuple, std::vector> COO_Matrix::getRowCopy(IdxT r) - { - if (!this->sorted_) - { - this->sortSparse(); - } - IdxT row_index = this->indexStartRow(r); - - if (row_index == -1) - { - return {std::vector(), std::vector()}; - } - - IdxT rsize = row_index; - do - { - rsize++; - } while (rsize < this->values_.size() && this->row_indices_[rsize] == r); - - return {{this->column_indices_.begin() + row_index, this->column_indices_.begin() + rsize}, - {this->values_.begin() + row_index, this->values_.begin() + rsize}}; - } - - /** - * @brief Get all entry pointers. Will sort before returning - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @return std::tuple, std::vector, std::vector> - */ - template - inline std::tuple&, std::vector&, std::vector&> COO_Matrix::getEntries(const bool sort) - { - if (!this->sorted_ && sort) - { - this->sortSparse(); - } - return {this->row_indices_, this->column_indices_, this->values_}; - } - - /** - * @brief Sorts the data if it's not already sorted - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @return std::tuple, std::vector, std::vector> - */ - template - inline std::tuple, std::vector, std::vector> COO_Matrix::getEntryCopies() - { - if (!this->sorted_) - { - this->sortSparse(); - } - return {this->row_indices_, this->column_indices_, this->values_}; - } - - /** - * @brief Returns the data in CSR Format - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @return std::tuple, std::vector, std::vector> - */ - template - inline std::tuple, std::vector, std::vector> COO_Matrix::setDataToCSR() - { - if (!this->isSorted()) - this->sortSparse(); - std::vector row_size_vec(this->rows_size_ + 1, 0); - IdxT counter = 0; - for (IdxT i = 0; i < static_cast(row_size_vec.size() - 1); i++) - { - row_size_vec[i + 1] = row_size_vec[i]; - while (counter < static_cast(this->row_indices_.size()) && i == this->row_indices_[counter]) - { - row_size_vec[i + 1]++; - counter++; - } - } - return {row_size_vec, this->column_indices_, this->values_}; - } - - /** - * @brief Only creates the row data - * - * @todo swap this with having the matrix store the data and updates. This can then be passed by reference - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @return std::vector - */ - template - inline std::vector COO_Matrix::getCSRRowData() - { - if (!this->isSorted()) - this->sortSparse(); - std::vector row_size_vec(static_cast(this->rows_size_ + 1), 0); - size_t counter = 0; - for (size_t i = 0; i < row_size_vec.size() - 1; i++) - { - row_size_vec[i + 1] = row_size_vec[i]; - while (counter < this->row_indices_.size() && i == static_cast(this->row_indices_[counter])) - { - row_size_vec[i + 1]++; - counter++; - } - } - return row_size_vec; - } - - /** - * @brief Set coordinates and values of the matrix. - * - * Matrix entries will be sorted in row-major order before the method returns. - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] r row indices of the matrix - * @param[in] c column indices of the matrix - * @param[in] v values of the matrix - * - * @pre r.size() == c.size() == v.size() - * @pre r,c,v represent an array in COO format - * - * @post Coordinates and values are set in the matrix. - */ - template - inline void COO_Matrix::setValues(std::vector r, std::vector c, std::vector v) - { - // sort input - this->sortSparseCOO(r, c, v, r.size()); - - // Duplicated with axpy. Could replace with function depdent on lambda expression - size_t a_iter = 0; - // iterate for all current values_ in matrix - for (size_t i = 0; i < this->row_indices_.size(); i++) - { - // pushback values_ when they are not in current matrix - while (a_iter < r.size() && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) - { - this->row_indices_.push_back(r[a_iter]); - this->column_indices_.push_back(c[a_iter]); - this->values_.push_back(v[a_iter]); - this->checkIncreaseSize(r[a_iter], c[a_iter]); - a_iter++; - } - if (a_iter >= r.size()) - { - break; - } - - if (r[a_iter] == this->row_indices_[i] && c[a_iter] == this->column_indices_[i]) - { - this->values_[i] = v[a_iter]; - a_iter++; - } - } - // push back rest that was not found sorted - for (size_t i = a_iter; i < r.size(); i++) - { - this->row_indices_.push_back(r[i]); - this->column_indices_.push_back(c[i]); - this->values_.push_back(v[i]); - - this->checkIncreaseSize(r[i], c[i]); - } - - this->sorted_ = false; - } - - /** - * @brief Append coordinates and values of the matrix without sorting or deduplicating. - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] r row indices to be stored - * @param[in] c column indices to be stored - * @param[in] v values to be stored - * @param[in] nnz to be stored - * - */ - template - inline void COO_Matrix::setValues(RealT alpha, IdxT* r, IdxT* c, RealT* v, IdxT nnz) - { - for (size_t i = 0; i < static_cast(nnz); i++) - { - this->row_indices_.emplace_back(r[i]); - this->column_indices_.emplace_back(c[i]); - this->values_.emplace_back(alpha * v[i]); - - this->checkIncreaseSize(r[i], c[i]); - } - this->sorted_ = false; - } - - /** - * @brief Implements axpy this += alpha * a. Will sort before running - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] alpha matrix to be added - * @param[in] a scalar to multiply by - * - * @post this = this + alpha * a - */ - template - inline void COO_Matrix::axpy(RealT alpha, COO_Matrix& a, const bool sort) - { - if (!this->sorted_ && sort) - { - this->sortSparse(); - } - if (!a.isSorted()) - { - a.sortSparse(); - } - IdxT m = 0; - IdxT n = 0; - std::tuple&, std::vector&, std::vector&> tpm = a.getEntries(); - const auto& [r, c, val] = tpm; - std::tie(m, n) = a.getDimensions(); - - // Increase size as necessary - this->rows_size_ = this->rows_size_ > m ? this->rows_size_ : m; - this->columns_size_ = this->columns_size_ > n ? this->columns_size_ : n; - - size_t a_iter = 0; - // iterate for all current values in matrix - for (size_t i = 0; i < this->row_indices_.size(); i++) - { - if (sort) // highjacking sort variable to signify that the sparsity pattern can change - { - // pushback values when they are not in current matrix - while (a_iter < r.size() && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) - { - this->row_indices_.push_back(r[a_iter]); - this->column_indices_.push_back(c[a_iter]); - this->values_.push_back(alpha * val[a_iter]); - - this->checkIncreaseSize(r[a_iter], c[a_iter]); - a_iter++; - } - } - - if (a_iter >= r.size()) - { - break; - } - - if (r[a_iter] == this->row_indices_[i] && c[a_iter] == this->column_indices_[i]) - { - this->values_[i] += alpha * val[a_iter]; - a_iter++; - } - } - // push back rest that was not found sorted_ - for (size_t i = a_iter; i < r.size(); i++) - { - this->row_indices_.push_back(r[i]); - this->column_indices_.push_back(c[i]); - this->values_.push_back(alpha * val[i]); - - this->checkIncreaseSize(r[i], c[i]); - } - - this->sorted_ = false; - } - - /** - * @brief axpy on a COO representation of a matrix. Will sort before running - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param alpha scalar to multiply by - * @param r row indices - * @param c column indices - * @param v values - * - * @pre r.size() == c.size() == v.size() - * @pre r,c,v represent an array a in COO format - * - * @post this = this + alpha * a - */ - template - inline void COO_Matrix::axpy(RealT alpha, std::vector r, std::vector c, std::vector v, const bool sort) - { - if (!this->sorted_ && sort) - { - this->sortSparse(); - } - - // sort input - this->sortSparseCOO(r, c, v, r.size()); - - size_t a_iter = 0; - // iterate for all current values_ in matrix - for (size_t i = 0; i < this->row_indices_.size(); i++) - { - if (sort) // highjacking sort variable to signify that the sparsity pattern can change - { - // pushback values_ when they are not in current matrix - while (a_iter < r.size() && (r[a_iter] < this->row_indices_[i] || (r[a_iter] == this->row_indices_[i] && c[a_iter] < this->column_indices_[i]))) - { - this->row_indices_.push_back(r[a_iter]); - this->column_indices_.push_back(c[a_iter]); - this->values_.push_back(alpha * v[a_iter]); - - this->checkIncreaseSize(r[a_iter], c[a_iter]); - a_iter++; - } - } - - if (a_iter >= r.size()) - { - break; - } - - if (r[a_iter] == this->row_indices_[i] && c[a_iter] == this->column_indices_[i]) - { - this->values_[i] += alpha * v[a_iter]; - a_iter++; - } - } - // push back rest that was not found sorted_ - for (size_t i = a_iter; i < r.size(); i++) - { - this->row_indices_.push_back(r[i]); - this->column_indices_.push_back(c[i]); - this->values_.push_back(alpha * v[i]); - - this->checkIncreaseSize(r[i], c[i]); - } - - this->sorted_ = false; - } - - /** - * @brief Append coordinates and values of the matrix without sorting or deduplicating. - * - * Need to deduplicate bus entries right away to get around the new mappings introduced by bus faults. - * - * Only use for buses, which are 2x2, so the cost of loops is not uncontrolled. Don't allow resizing. - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] r row indices to be stored - * @param[in] c column indices to be stored - * @param[in] v values to be stored - * @param[in] nnz to be stored - * - */ - template - inline void COO_Matrix::axpy(RealT alpha, IdxT* r, IdxT* c, RealT* v, IdxT nnz) - { - if (this->row_indices_.size() == 0) // Do nothing for infinite bus - { - } - else if (this->row_indices_.size() == 4) - { - for (size_t i = 0; i < this->row_indices_.size(); i++) - { - for (size_t j = 0; j < static_cast(nnz); j++) - { - if (this->row_indices_[i] == r[j] && this->column_indices_[i] == c[j]) - { - this->values_[i] += alpha * v[j]; - } - } - } - } - else - { - std::cout << "Warning: Unexpected size in axpy\n"; - } - - this->sorted_ = false; - } - - /** - * @brief Deduplicate matrix entries - * - * @node This is currently only used in component-level Jacobian tests, - * which don't use CooMatrix or CsrMatrix yet, to compare with Dependency::Tracking maps. - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - */ - template - inline void COO_Matrix::deduplicate() - { - std::vector row_temp = this->row_indices_; - std::vector col_temp = this->column_indices_; - std::vector val_temp = this->values_; - - this->row_indices_.clear(); - this->column_indices_.clear(); - this->values_.clear(); - - IdxT nnz = 0; - for (size_t i = 0; i < row_temp.size(); i++) - { - IdxT row = row_temp[i]; - IdxT col = col_temp[i]; - RealT val = val_temp[i]; - bool exists = false; - for (size_t j = 0; j < nnz; j++) - { - if (this->row_indices_[j] == row && this->column_indices_[j] == col) - { - this->values_[j] += val; - exists = true; - } - } - if (!exists) - { - this->row_indices_.emplace_back(row); - this->column_indices_.emplace_back(col); - this->values_.emplace_back(val); - nnz++; - } - } - } - - /** - * @brief Scale all values by alpha - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] alpha scalar to scale by - */ - template - inline void COO_Matrix::scal(RealT alpha) - { - for (auto i = this->values_.begin(); i < this->values_.end(); i++) - { - *i *= alpha; - } - } - - /** - * @brief Calculates the Frobenius Norm of the matrix - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @return RealT - Frobenius Norm of the matrix - */ - template - inline RealT COO_Matrix::frobNorm() - { - RealT totsum = 0.0; - for (auto i = this->values_.begin(); i < this->values_.end(); i++) - { - totsum += abs(*i) ^ 2; - } - return totsum; - } - - /** - * @brief Permutate the matrix to a different one. Only changes the coordinates - * - * @tparam RealT - Real type for Jacobian entries - * - * @param[in] row_perm - * @param[out] col_perm - * - * @pre row_perm.size() == this->rows_size_ = col_perm.size() == this->columns_size_ - * - * @post this = this(row_perm, col_perm) - */ - template - inline void COO_Matrix::permutation(std::vector row_perm, std::vector col_perm) - { - assert(row_perm.size() = this->rows_size_); - assert(col_perm.size() = this->columns_size_); - - for (int i = 0; i < this->values_.size(); i++) - { - this->row_indices_[i] = row_perm[this->row_indices_[i]]; - this->column_indices_[i] = col_perm[this->column_indices_[i]]; - } - this->sorted_ = false; - // cycle sorting maybe useful since permutations are already known - } - - /** - * @brief Permutes the matrix and can change its size efficently - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] row_perm row permutation - * @param[in] col_perm column permutation - * @param[in] m number of rows - * @param[in] n number of columns - * - * @pre row_perm.size() == this->rows_size_ - * @pre col_perm.size() == this->columns_size_ - * @pre indices are set to -1 if they are to be removed - * - * @post this = this(row_perm, col_perm) and removed indices have corresponding values set to 0 - */ - template - inline void COO_Matrix::permutationSizeMap(std::vector row_perm, std::vector col_perm, IdxT m, IdxT n) - { - assert(row_perm.size() == this->rows_size_); - assert(col_perm.size() == this->columns_size_); - - this->rows_size_ = m; - this->columns_size_ = n; - - for (size_t i = 0; i < this->values_.size(); i++) - { - if (row_perm[this->row_indices_[i]] == -1 || col_perm[this->column_indices_[i]] == -1) - { - this->values_[i] = 0; - } - else - { - this->row_indices_[i] = row_perm[this->row_indices_[i]]; - this->column_indices_[i] = col_perm[this->column_indices_[i]]; - } - } - this->sorted_ = false; - } - - /** - * @brief Turn matrix into the null matrix. Does not actually delete memory - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - */ - template - inline void COO_Matrix::zeroMatrix() - { - // resize doesn't affect capacity if smaller - this->column_indices_.resize(0); - this->row_indices_.resize(0); - this->values_.resize(0); - this->sorted_ = true; - } - - /** - * @brief Initializes matrix values to 0 without changing the sparsity pattern - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - */ - template - inline void COO_Matrix::zeroValuedMatrix() - { - for (size_t i = 0; i < this->row_indices_.size(); i++) - { - this->values_[i] = 0.0; - } - } - - /** - * @brief Turn matrix into the identity matrix - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] n size of the identity matrix - * - * @post this = I_n - * - * @todo - it might be better to explicitly zero out the matrix and require to be so in preconditions - */ - - template - inline void COO_Matrix::identityMatrix(IdxT n) - { - // Reset Matrix - this->zeroMatrix(); - for (IdxT i = 0; i < n; i++) - { - this->column_indices_[i] = i; - this->row_indices_[i] = i; - this->values_[i] = 1.0; - } - this->sorted_ = true; - } - - /** - * @brief Restructure the sparse matrix for faster accesses and modifications - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - */ - template - inline void COO_Matrix::sortSparse() - { - this->sortSparseCOO(this->row_indices_, this->column_indices_, this->values_, (this->row_indices_).size()); - this->sorted_ = true; - } - - /** - * @brief Check if the matrix is sorted - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[out] bool - true if sorted, false otherwise - */ - - template - inline bool COO_Matrix::isSorted() - { - return this->sorted_; - } - - /** - * @brief Get the number of non-zero elements in the matrix - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[out] IdxT - number of non-zero elements in the matrix - */ - template - inline IdxT COO_Matrix::nnz() - { - return static_cast(this->values_.size()); - } - - template - inline std::tuple COO_Matrix::getDimensions() - { - return std::tuple(this->rows_size_, this->columns_size_); - } - - /** - * @brief Print matrix in sorted order - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] name to identify the specific matrix printed - */ - template - inline void COO_Matrix::printMatrix(std::string name) - { - if (this->sorted_ == false) - { - this->sortSparse(); - } - - std::cout << "Sparse COO Matrix: " << name << "\n"; - std::cout << "(x , y, value)\n"; - for (size_t i = 0; i < this->values_.size(); i++) - { - std::cout << "(" << this->row_indices_[i] - << ", " << this->column_indices_[i] - << ", " << this->values_[i] << ")\n"; - } - std::cout << std::flush; - } - - /** - * @brief Print matrix to file in Matrix Market format - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] filename Output file path - * @param[in] comment Optional comment for the Matrix Market file - */ - template - void COO_Matrix::printMatrixMarket(const std::string& filename, const std::string& comment) - { - if (this->sorted_ == false) - { - this->sortSparse(); - } - - std::ofstream outfile(filename); - if (!outfile.is_open()) - { - std::cerr << "Error: Unable to open file " << filename << " for writing." << std::endl; - return; - } - - // Write Matrix Market header - outfile << "%%MatrixMarket matrix coordinate "; - - if (!std::is_floating_point::value) - { - // error message - return; - } - - outfile << "general" << std::endl; - - // Write comment if provided - if (!comment.empty()) - outfile << "% " << comment << std::endl; - - // Write dimensions and number of entries - // Get max row and column indices to determine matrix dimensions - IdxT max_row = 0; - IdxT max_col = 0; - for (size_t i = 0; i < this->values_.size(); i++) - { - max_row = std::max(max_row, this->row_indices_[i]); - max_col = std::max(max_col, this->column_indices_[i]); - } - - // Matrix Market uses 1-based indexing, so add 1 to dimensions - outfile << (max_row + 1) << " " << (max_col + 1) << " " << this->values_.size() << std::endl; - - // Write the matrix entries - for (size_t i = 0; i < this->values_.size(); i++) - { - // Matrix Market uses 1-based indexing, so add 1 to indices - outfile << (this->row_indices_[i] + 1) << " " - << (this->column_indices_[i] + 1) << " " - << this->values_[i] << std::endl; - } - - outfile.close(); - - std::cout << "Matrix Market file written to: " << filename << std::endl; - } - - /** - * @brief Find the lowest row cordinate from set of provided coordinates - * - * Assumes rows and columns are sorted - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] rows - row indices - * @param[in] r - row index - * - * @return IdxT - index of lowest row - */ - template - inline IdxT COO_Matrix::indexStartRow(const std::vector& rows, IdxT r) - { - // Specialized Binary Search for Lowest Row - IdxT i1 = 0; - IdxT i2 = rows->size() - 1; - IdxT m_smallest = -1; - IdxT m = -1; - while (i1 <= i2) - { - m = (i2 + i1) / 2; - // rows - if (rows[m] < r) - { - i1 = m + 1; - } - else if (r < rows[m]) - { - i2 = m - 1; - } - else - { - if (i1 == i2) - { - return m_smallest; - } - - // Keep track of smallest cordinate - m_smallest = m; - i2 = m - 1; - } - } - return m_smallest; - } - - /** - * @brief Basic binary search - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] rows - row indices - * @param[in] columns - column indices - * @param[in] ri - row index - * @param[in] ci - column index - * @return IdxT - returns the index of the coordinate - */ - template - inline IdxT COO_Matrix::sparseCordBinarySearch(const std::vector& rows, const std::vector& columns, IdxT ri, IdxT ci) - { - assert(rows.size() == columns.size()); - // basic binary search - IdxT i1 = 0; - IdxT i2 = rows.size() - 1; - IdxT m = 0; - while (i1 <= i2) - { - m = (i2 + i1) / 2; - // rows - if (rows[m] < ri) - { - i1 = m + 1; - } - else if (ri < rows[m]) - { - i2 = m - 1; - } - else - { - if (columns[m] < ci) - { - i1 = m + 1; - } - else if (ci < columns[m]) - { - i2 = m - 1; - } - break; - } - } - - return m; - } - - /** - * @brief Check if the size of the matrix needs to be increased - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] r row index - * @param[in] c column index - * @return true if size was increased - */ - - template - inline bool COO_Matrix::checkIncreaseSize(IdxT r, IdxT c) - { - bool changed = false; - if (r + 1 > this->rows_size_) - { - this->rows_size_ = r + 1; - changed = true; - } - if (c + 1 > this->columns_size_) - { - this->columns_size_ = c + 1; - changed = true; - } - - return changed; - } - - /** - * @brief Sorts unordered COO matrix - * - * Matrix entries can appear in arbitrary order and will be sorted in - * row-major order before the method returns. - * Duplicate entries are not allowed and should be pre-summed. - * - * @pre rows, columns, and values are of the same size and represent a COO matrix with no duplicates - * @post Matrix entries are sorted in row-major order - * - * @todo simple setup. Should add stable sorting since lists are pre-sorted_ - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param rows - * @param columns - * @param values - */ - template - inline void COO_Matrix::sortSparseCOO(std::vector& rows, std::vector& columns, std::vector& values, size_t nnz) - { - - // index based sort code - // https://stackoverflow.com/questions/25921706/creating-a-vector-of-indices-of-a-sorted_-vector - // cannot call sort since two arrays are used instead - std::vector ordervec(nnz); - std::size_t n(0); - std::generate(std::begin(ordervec), std::end(ordervec), [&] - { return n++; }); - - // Sort by row first then column. - std::sort(std::begin(ordervec), - std::end(ordervec), - [&](auto i1, auto i2) - { return (rows[i1] < rows[i2]) || (rows[i1] == rows[i2] && columns[i1] < columns[i2]); }); - - // reorder based of index-sorting. Only swap cost no extra memory. - // @todo see if extra memory creation is fine - // https://stackoverflow.com/a/22183350 - for (size_t i = 0; i < ordervec.size(); i++) - { - // permutation swap - while (ordervec[i] != ordervec[ordervec[i]]) - { - std::swap(rows[ordervec[i]], rows[ordervec[ordervec[i]]]); - std::swap(columns[ordervec[i]], columns[ordervec[ordervec[i]]]); - std::swap(values[ordervec[i]], values[ordervec[ordervec[i]]]); - - // swap orderings - std::swap(ordervec[i], ordervec[ordervec[i]]); - } - } - } - - /** - * @brief Constructor for COO Matrix with given cooridnates and values - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] r row indices - * @param[in] c column indices - * @param[in] v values - * @param[in] m number of rows - * @param[in] n number of columns - * - * @pre r.size() == c.size() == v.size() - * @pre r,c,v represent an array in COO format - * - * @post COO_Matrix is created with given coordinates and values - */ - - template - inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, IdxT m, IdxT n) - { - this->values_ = v; - this->row_indices_ = r; - this->column_indices_ = c; - this->rows_size_ = m; - this->columns_size_ = n; - this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. - } - - /** - * @brief Constructor for empty COO Matrix of a given size - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @param[in] m number of rows - * @param[in] n number of columns - * - * @post empty COO Matrix is created with given size - */ - - template - inline COO_Matrix::COO_Matrix(IdxT m, IdxT n) - { - this->rows_size_ = m; - this->columns_size_ = n; - this->values_ = std::vector(); - this->row_indices_ = std::vector(); - this->column_indices_ = std::vector(); - this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. - } - - /** - * @brief Constructor for empty COO Matrix of size 0 - * - * @tparam RealT - Real type for Jacobian entries - * @tparam IdxT - Integer data type for matrix indices - * - * @post empty COO Matrix of size 0 is created - */ - - template - inline COO_Matrix::COO_Matrix() - { - this->rows_size_ = 0; - this->columns_size_ = 0; - this->values_ = std::vector(); - this->row_indices_ = std::vector(); - this->column_indices_ = std::vector(); - this->sorted_ = false; // Set to false until explicitly sorted, though logically it is sorted. - } - - template - COO_Matrix::~COO_Matrix() - { - } - } // namespace LinearAlgebra -} // namespace GridKit diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index 3635c7027..b7953a277 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -3,7 +3,6 @@ #include #include -#include #include #include #include @@ -24,8 +23,8 @@ namespace GridKit using ScalarT = scalar_type; using IdxT = index_type; using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; //\todo Use CsrMatrix using CsrMatrixT = GridKit::LinearAlgebra::CsrMatrix; + using CooMatrixT = GridKit::LinearAlgebra::CooMatrix; Evaluator() { @@ -138,9 +137,6 @@ namespace GridKit virtual std::vector& getResidual() = 0; virtual const std::vector& getResidual() const = 0; - virtual MatrixT& getJacobian() = 0; - virtual const MatrixT& getJacobian() const = 0; - virtual std::vector& getIntegrand() = 0; virtual const std::vector& getIntegrand() const = 0; diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index 84149f8a4..492c92930 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -51,7 +51,6 @@ namespace GridKit using Component::f_; using Component::wb_; using Component::h_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp index a6cd306b2..72aa13751 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp @@ -23,83 +23,79 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - J_rows_buffer_ = new IdxT[4]; - J_cols_buffer_ = new IdxT[4]; - J_vals_buffer_ = new RealT[4]; + auto bus1_size = static_cast(bus1_->size()); + auto bus2_size = static_cast(bus2_->size()); + auto buffer_size = (bus1_size + bus2_size) * (bus1_size + bus2_size); + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + // Bus 1 diagonal Jacobian block owned by the bus GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11, - ScalarT, - IdxT>::eval(this, - static_cast(bus1_->size()), - static_cast((bus1_->y()).size()), - (bus1_->getResidualIndices()).data(), - (bus1_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus1_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - bus1_->getJacobian()); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11>::eval(this, + static_cast(bus1_->size()), + static_cast((bus1_->y()).size()), + (bus1_->getResidualIndices()).data(), + (bus1_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus1_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); // Bus 2 diagonal Jacobian block owned by the bus GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22, - ScalarT, - IdxT>::eval(this, - static_cast(bus2_->size()), - static_cast((bus2_->y()).size()), - (bus2_->getResidualIndices()).data(), - (bus2_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus2_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - bus2_->getJacobian()); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22>::eval(this, + static_cast(bus2_->size()), + static_cast((bus2_->y()).size()), + (bus2_->getResidualIndices()).data(), + (bus2_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus2_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); // Off-diagonal Jacobian block (Bus2 variables) owned by the branch GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12, - ScalarT, - IdxT>::eval(this, - static_cast(bus1_->size()), - static_cast((bus2_->y()).size()), - (bus1_->getResidualIndices()).data(), - (bus2_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus2_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_, - true); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12>::eval(this, + static_cast(bus1_->size()), + static_cast((bus2_->y()).size()), + (bus1_->getResidualIndices()).data(), + (bus2_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus2_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); // Off-diagonal Jacobian block (Bus1 variables) owned by the branch GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21, - ScalarT, - IdxT>::eval(this, - static_cast(bus2_->size()), - static_cast((bus1_->y()).size()), - (bus2_->getResidualIndices()).data(), - (bus1_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus1_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_, - true); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21>::eval(this, + static_cast(bus2_->size()), + static_cast((bus1_->y()).size()), + (bus2_->getResidualIndices()).data(), + (bus1_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus1_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp index 0da7765f7..13c24c0be 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp @@ -20,22 +20,21 @@ namespace GridKit { using BusBase::bus_id_; using BusBase::size_; + using BusBase::nnz_; using BusBase::y_; using BusBase::yp_; using BusBase::f_; - using BusBase::J_; - using BusBase::J_rows_buffer_; - using BusBase::J_cols_buffer_; - using BusBase::J_vals_buffer_; using BusBase::tag_; using BusBase::variable_indices_; using BusBase::residual_indices_; + using BusBase::coo_jac_; using BusBase::monitor_; public: using ScalarT = scalar_type; using IdxT = index_type; using RealT = typename BusBase::RealT; + using CooMatrixT = typename BusBase::CooMatrixT; using MonitorT = typename BusBase::MonitorT; using ModelDataT = BusData; using BusTypeT = typename BusData::BusType; @@ -97,6 +96,35 @@ namespace GridKit return f_[1]; } + protected: + int constructCoo() + { + if (coo_jac_ == nullptr) + { + IdxT num_rows = 0; + IdxT num_cols = 0; + for (IdxT i = 0; i < nnz_; ++i) + { + if (J_rows_buffer_[i] > num_rows) + { + num_rows = J_rows_buffer_[i]; + } + if (J_cols_buffer_[i] > num_cols) + { + num_cols = J_cols_buffer_[i]; + } + } + coo_jac_ = new CooMatrixT(num_rows, num_cols, nnz_); + coo_jac_->setDataPointers(J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, LinearAlgebra::memory::HOST); + } + + return 0; + } + + IdxT* J_rows_buffer_{nullptr}; + IdxT* J_cols_buffer_{nullptr}; + RealT* J_vals_buffer_{nullptr}; + private: ScalarT Vr0_{0.0}; ScalarT Vi0_{0.0}; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp index bcb02889c..4af2d5a35 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp @@ -11,11 +11,11 @@ namespace GridKit namespace PhasorDynamics { /** - * @brief Jacobian evaluation experimental. This sets values to 0, for other - * components to add their contributions. + * @brief Jacobian evaluation experimental. * - * @warning This implementation assumes bus Jacobians are always evaluated - * _before_ component model Jacobians. + * This sets values to 0, and these remain unchanged. It is needed to get + * the indices into the list of entries that will later be deduplicated. + * Contributions to bus Jacobians from other components are stored in those components. * * @return int - error code */ @@ -27,10 +27,6 @@ namespace GridKit if (J_rows_buffer_ == nullptr) { - J_.zeroMatrix(); - - // Reserve space for a dense matrix of size_*size_. - // Enyme will compute the appropriate nnz from sparsification. J_rows_buffer_ = new IdxT[4]; J_cols_buffer_ = new IdxT[4]; J_vals_buffer_ = new RealT[4]; @@ -47,11 +43,9 @@ namespace GridKit J_vals_buffer_[1] = 0.0; J_vals_buffer_[2] = 0.0; J_vals_buffer_[3] = 0.0; - J_.setValues(1.0, J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, 4); //< @todo: Update once sparse storage format changes - } - else - { - J_.zeroValuedMatrix(); + + nnz_ = 4; + this->constructCoo(); } return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 534af06ed..05a3f7808 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -83,6 +83,21 @@ namespace GridKit Bus::~Bus() { // std::cout << "Destroy PQ bus ..." << std::endl; + if (J_rows_buffer_ != nullptr) + { + delete[] J_rows_buffer_; + delete[] J_cols_buffer_; + delete[] J_vals_buffer_; + J_rows_buffer_ = nullptr; + J_cols_buffer_ = nullptr; + J_vals_buffer_ = nullptr; + } + + if (coo_jac_ != nullptr) + { + delete coo_jac_; + coo_jac_ = nullptr; + } } /*! diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp index 5e1ff1d3a..3abec94ff 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp @@ -18,18 +18,20 @@ namespace GridKit { using BusBase::bus_id_; using BusBase::size_; + using BusBase::nnz_; using BusBase::y_; using BusBase::yp_; using BusBase::f_; - using BusBase::J_; using BusBase::variable_indices_; using BusBase::residual_indices_; + using BusBase::coo_jac_; using BusBase::monitor_; public: using ScalarT = scalar_type; using IdxT = index_type; using RealT = typename BusBase::RealT; + using CooMatrixT = typename BusBase::CooMatrixT; using MonitorT = typename BusBase::MonitorT; using ModelDataT = BusData; using BusTypeT = typename BusData::BusType; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp index f84e99def..1f3a151b8 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfiniteImpl.hpp @@ -81,6 +81,11 @@ namespace GridKit template BusInfinite::~BusInfinite() { + if (coo_jac_ != nullptr) + { + delete coo_jac_; + coo_jac_ = nullptr; + } } /** @@ -148,6 +153,11 @@ namespace GridKit template int BusInfinite::evaluateJacobian() { + if (coo_jac_ == nullptr) + { + nnz_ = 0; + coo_jac_ = new CooMatrixT(0, 0, 0); + } return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index f6591322c..05d8fcd3b 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -27,12 +27,13 @@ namespace GridKit class BusBase : public Model::Evaluator { public: - using ScalarT = scalar_type; - using IdxT = index_type; - using RealT = typename Model::Evaluator::RealT; - using MatrixT = typename Model::Evaluator::MatrixT; - using BusTypeT = typename BusData::BusType; - using MonitorT = Model::VariableMonitor; + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Model::Evaluator::RealT; + using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; + using CooMatrixT = typename Model::Evaluator::CooMatrixT; + using BusTypeT = typename BusData::BusType; + using MonitorT = Model::VariableMonitor; BusBase() = default; @@ -104,16 +105,6 @@ namespace GridKit return f_; } - MatrixT& getJacobian() override - { - return J_; - } - - const MatrixT& getJacobian() const override - { - return J_; - } - int setVariableIndex(IdxT local_index, IdxT global_index) { variable_indices_[static_cast(local_index)] = global_index; @@ -152,9 +143,19 @@ namespace GridKit return BusTypeT::DEFAULT; } + CsrMatrixT* getCsrJacobian() const override + { + return csr_jac_; + } + + CooMatrixT* getCooJacobian() const + { + return coo_jac_; + } + bool hasJacobian() override { - return false; + return true; } void updateTime(RealT /* t */, RealT /* a */) override @@ -194,10 +195,8 @@ namespace GridKit std::vector f_; std::vector g_; - MatrixT J_; - IdxT* J_rows_buffer_{nullptr}; - IdxT* J_cols_buffer_{nullptr}; - RealT* J_vals_buffer_{nullptr}; + CsrMatrixT* csr_jac_{nullptr}; + CooMatrixT* coo_jac_{nullptr}; RealT rel_tol_; RealT abs_tol_; diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp index 3b559fd1b..62c169581 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp @@ -24,18 +24,21 @@ namespace GridKit class BusFault : public Component { using Component::gridkit_component_id_; - using Component::alpha_; - using Component::nnz_; using Component::size_; - using Component::tag_; + using Component::nnz_; using Component::time_; + using Component::alpha_; using Component::y_; using Component::yp_; + using Component::tag_; using Component::wb_; using Component::h_; + using Component::f_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; + using Component::variable_indices_; + using Component::residual_indices_; public: using ScalarT = scalar_type; @@ -111,6 +114,7 @@ namespace GridKit public: __attribute__((always_inline)) inline int evaluateBusResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*); + __attribute__((always_inline)) inline int evaluateInternalResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*); private: BusT* bus_; diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp index 9d3c45854..0a4b72a38 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp @@ -23,31 +23,70 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - if (status_) + if (J_rows_buffer_ == nullptr) { - if (J_rows_buffer_ == nullptr) + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto buffer_size = size * size + 2 * size * bus_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; + } + + nnz_ = 0; + + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + IdxT nnz_tmp = nnz_; + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + if (!status_) // Value contributions from DfDwb only when status_ + { + for (IdxT i = nnz_tmp; i < nnz_; ++i) { - J_rows_buffer_ = new IdxT[4]; - J_cols_buffer_ = new IdxT[4]; - J_vals_buffer_ = new RealT[4]; + J_vals_buffer_[i] = 0.0; } - GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - (bus_->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - bus_->getJacobian()); } + GridKit::Enzyme::Sparse::DhDy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp index 8c1212d4a..2521a9d34 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp @@ -24,7 +24,7 @@ namespace GridKit : bus_(bus), R_(0), X_(0.01), status_(0), bus_id_(0) { (void) bus_id_; - size_ = 0; + size_ = 2; setDerivedParams(); } @@ -42,7 +42,7 @@ namespace GridKit BusFault::BusFault(BusT* bus, RealT R, RealT X, int status) : bus_(bus), R_(R), X_(X), status_(status), bus_id_(0) { - size_ = 0; + size_ = 2; setDerivedParams(); } @@ -81,11 +81,11 @@ namespace GridKit monitor_->set(Variable::state, [this] { return status_; }); monitor_->set(Variable::ir, [this] - { return Ir(); }); + { return y_[0]; }); monitor_->set(Variable::ii, [this] - { return Ii(); }); + { return y_[1]; }); - size_ = 0; + size_ = 2; setDerivedParams(); } @@ -112,9 +112,25 @@ namespace GridKit { // std::cout << "Allocate BusFault..." << std::endl; + auto size = static_cast(size_); // avoid compiler warnings + f_.resize(size); + y_.resize(size); + yp_.resize(size); + tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); + + // Resize coupling data wb_.resize(2); h_.resize(2); + // Default variable and residual index mapping to local index + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + return 0; } @@ -125,6 +141,24 @@ namespace GridKit template int BusFault::initialize() { + if (status_) + { + ScalarT vr = Vr(); + ScalarT vi = Vi(); + ScalarT ir = -(vr * G_ - vi * B_); + ScalarT ii = -(vr * B_ + vi * G_); + y_[0] = ir; + y_[1] = ii; + } + else + { + y_[0] = 0.0; + y_[1] = 0.0; + } + + yp_[0] = 0.0; + yp_[1] = 0.0; + return 0; } @@ -134,6 +168,9 @@ namespace GridKit template int BusFault::tagDifferentiable() { + tag_[0] = false; + tag_[1] = false; + return 0; } @@ -143,18 +180,40 @@ namespace GridKit */ template __attribute__((always_inline)) int BusFault::evaluateBusResidual( - [[maybe_unused]] ScalarT* y, [[maybe_unused]] ScalarT* yp, ScalarT* wb, ScalarT* h) + ScalarT* y, + [[maybe_unused]] ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* h) { - ScalarT Vr = wb[0]; - ScalarT Vi = wb[1]; - ScalarT Ir = -Vr * G_ + Vi * B_; - ScalarT Ii = -Vr * B_ - Vi * G_; + ScalarT Ir = y[0]; + ScalarT Ii = y[1]; h[0] = Ir; h[1] = Ii; return 0; } + /** + * @brief Internal residual + * + */ + template + __attribute__((always_inline)) int BusFault::evaluateInternalResidual( + ScalarT* y, + [[maybe_unused]] ScalarT* yp, + ScalarT* wb, + ScalarT* f) + { + ScalarT Vr = wb[0]; + ScalarT Vi = wb[1]; + ScalarT Ir = y[0]; + ScalarT Ii = y[1]; + f[0] = Ir + Vr * G_ - Vi * B_; + f[1] = Ii + Vr * B_ + Vi * G_; + + return 0; + } + /** * \brief Residual contribution of the branch is pushed to the * two terminal buses. @@ -167,10 +226,18 @@ namespace GridKit { wb_[0] = Vr(); wb_[1] = Vi(); + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), f_.data()); evaluateBusResidual(y_.data(), yp_.data(), wb_.data(), h_.data()); Ir() += h_[0]; Ii() += h_[1]; } + else + { + wb_[0] = 0.0; + wb_[1] = 0.0; + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), f_.data()); + } + return 0; } diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index 90ddb2abc..3a44ebb24 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -37,11 +37,10 @@ add_subdirectory(Exciter) add_subdirectory(Governor) add_subdirectory(Load) add_subdirectory(LoadZIP) +add_subdirectory(SignalNode) add_subdirectory(Stabilizer) add_subdirectory(SynchronousMachine) -add_subdirectory(SignalNode) - add_library(GridKit::phasor_dynamics_components ALIAS phasor_dynamics_components) add_library(GridKit::phasor_dynamics_components_dependency_tracking ALIAS phasor_dynamics_components_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 985c4eae1..8557dc4c8 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -17,12 +17,15 @@ namespace GridKit /** * @brief Component model implementation base class. */ - template - class Component : public Model::Evaluator + template + class Component : public Model::Evaluator { public: - using RealT = typename Model::Evaluator::RealT; - using MatrixT = typename Model::Evaluator::MatrixT; + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Model::Evaluator::RealT; + using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; + using CooMatrixT = typename Model::Evaluator::CooMatrixT; Component() = default; @@ -37,6 +40,24 @@ namespace GridKit J_cols_buffer_ = nullptr; J_vals_buffer_ = nullptr; } + + if (coo_jac_ != nullptr) + { + delete coo_jac_; + coo_jac_ = nullptr; + } + + if (csr_jac_ != nullptr) + { + delete csr_jac_; + csr_jac_ = nullptr; + } + + if (map_to_csr_ != nullptr) + { + delete[] map_to_csr_; + map_to_csr_ = nullptr; + } } virtual int verify() const = 0; @@ -102,16 +123,6 @@ namespace GridKit return f_; } - MatrixT& getJacobian() override - { - return J_; - } - - const MatrixT& getJacobian() const override - { - return J_; - } - int setVariableIndex(IdxT local_index, IdxT global_index) { variable_indices_[static_cast(local_index)] = global_index; @@ -144,6 +155,16 @@ namespace GridKit return residual_indices_; } + CsrMatrixT* getCsrJacobian() const override + { + return csr_jac_; + } + + CooMatrixT* getCooJacobian() const + { + return coo_jac_; + } + /// @todo Remove this method. It should be part of DynamicSolver class. bool hasJacobian() override { @@ -175,7 +196,56 @@ namespace GridKit return gridkit_component_id_; } + int constructCsr() + { + if (coo_jac_ == nullptr) + { + constructCoo(); + } + + if (csr_jac_ == nullptr) + { + IdxT* row_ptrs = coo_jac_->getCsrRowData(); + + nnz_ = coo_jac_->getNnz(); + + IdxT* cols = new IdxT[static_cast(nnz_)]; + RealT* vals = new RealT[static_cast(nnz_)]; + + std::copy(coo_jac_->getColData(), coo_jac_->getColData() + nnz_, cols); + std::copy(coo_jac_->getValues(), coo_jac_->getValues() + nnz_, vals); + + csr_jac_ = new CsrMatrixT(coo_jac_->getNumRows(), coo_jac_->getNumColumns(), nnz_, &row_ptrs, &cols, &vals); + } + + return 0; + } + protected: + int constructCoo() + { + if (coo_jac_ == nullptr) + { + IdxT num_rows = 0; + IdxT num_cols = 0; + for (IdxT i = 0; i < nnz_; ++i) + { + if (J_rows_buffer_[i] + 1 > num_rows) + { + num_rows = J_rows_buffer_[i] + 1; + } + if (J_cols_buffer_[i] + 1 > num_cols) + { + num_cols = J_cols_buffer_[i] + 1; + } + } + coo_jac_ = new CooMatrixT(num_rows, num_cols, nnz_); + coo_jac_->setDataPointers(J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, LinearAlgebra::memory::HOST); + } + + return 0; + } + IdxT size_{0}; IdxT nnz_{0}; /// Global (system-level) variable indices @@ -189,10 +259,12 @@ namespace GridKit std::vector f_; std::vector g_; - MatrixT J_; - IdxT* J_rows_buffer_{nullptr}; - IdxT* J_cols_buffer_{nullptr}; - RealT* J_vals_buffer_{nullptr}; + IdxT* J_rows_buffer_{nullptr}; + IdxT* J_cols_buffer_{nullptr}; + RealT* J_vals_buffer_{nullptr}; + IdxT* map_to_csr_{nullptr}; + CsrMatrixT* csr_jac_{nullptr}; + CooMatrixT* coo_jac_{nullptr}; RealT rel_tol_; RealT abs_tol_; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index add41ad66..f0b93e343 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -77,7 +77,6 @@ namespace GridKit using Component::y_; using Component::yp_; using Component::wb_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp index 31ef04039..1542ba8e3 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp @@ -25,67 +25,83 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - // Reserve space for a dense matrix of size_*size_. + // Reserve space for the dense blocks. // Enyme will compute the appropriate nnz from sparsification. - J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp index 56f5667d0..5406b89a1 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp @@ -64,7 +64,6 @@ namespace GridKit using Component::y_; using Component::yp_; using Component::wb_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp index 2cc959481..f73f6762a 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp @@ -20,65 +20,83 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for SexsPti..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + // Reserve space for the dense blocks. + // Enyme will compute the appropriate nnz from sparsification. + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index f8ceb7207..a919be9cd 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -69,7 +69,6 @@ namespace GridKit using Component::yp_; using Component::wb_; using Component::h_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp index 115670d4f..2202730a7 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp @@ -25,50 +25,67 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - // Reserve space for a dense matrix of size_*size_. + // Reserve space for the dense blocks. // Enyme will compute the appropriate nnz from sparsification. - J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + auto size = static_cast(size_); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/Load.hpp b/GridKit/Model/PhasorDynamics/Load/Load.hpp index 33ce3cc1e..e21bc8500 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.hpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.hpp @@ -40,12 +40,13 @@ namespace GridKit using Component::wb_; using Component::h_; using Component::f_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; using Component::variable_indices_; using Component::residual_indices_; + using Component::csr_jac_; + using Component::map_to_csr_; public: using ScalarT = scalar_type; diff --git a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp index 58ea9c1f5..df740f88d 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp @@ -23,63 +23,61 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Load..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - J_rows_buffer_ = new IdxT[4]; - J_cols_buffer_ = new IdxT[4]; - J_vals_buffer_ = new RealT[4]; + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto buffer_size = size * size + 2 * size * bus_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } - // DfDy call without alpha_ to indicate that df/dy' is null - // @todo: deduce from a compile-time differential tag + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - y_.size(), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp index ebe3a9be2..8b0f0adbf 100644 --- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp +++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp @@ -40,7 +40,6 @@ namespace GridKit using Component::wb_; using Component::h_; using Component::f_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp index a5f7e6dfc..93a345fe9 100644 --- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp @@ -18,63 +18,61 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for LoadZIP..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - J_rows_buffer_ = new IdxT[4]; - J_cols_buffer_ = new IdxT[4]; - J_vals_buffer_ = new RealT[4]; + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto buffer_size = size * size + 2 * size * bus_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } - // DfDy call without alpha_ to indicate that df/dy' is null - // @todo: deduce from a compile-time differential tag + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - y_.size(), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index 486f1820b..09051d7af 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -71,7 +71,6 @@ namespace GridKit using Component::yp_; using Component::wb_; using Component::h_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp index 1ea9e6434..8992b7995 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp @@ -27,50 +27,67 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Ieeest..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - // Reserve space for a dense matrix of size_*size_. - // Enzyme will compute the appropriate nnz from sparsification. - J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + // Reserve space for the dense blocks. + // Enyme will compute the appropriate nnz from sparsification. + auto size = static_cast(size_); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp index 73b6f764d..0703c3b22 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp @@ -80,7 +80,6 @@ namespace GridKit using Component::yp_; using Component::wb_; using Component::h_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp index 2fb91a858..95f17da98 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp @@ -23,99 +23,111 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - // Reserve space for a dense matrix of size_*size_. + // Reserve space for the dense blocks. // Enyme will compute the appropriate nnz from sparsification. - J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * signal_size + bus_size * bus_size + 2 * size * bus_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - (bus_->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - bus_->getJacobian()); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + (bus_->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - y_.size(), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp index 8e00fc2e8..d9f90d2be 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp @@ -75,7 +75,6 @@ namespace GridKit using Component::yp_; using Component::wb_; using Component::h_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp index deb4b1f5a..0c4a882c7 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp @@ -23,99 +23,111 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for Gensal..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - // Reserve space for a dense matrix of size_*size_. - // Enzyme will compute the appropriate nnz from sparsification. - J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + // Reserve space for the dense blocks. + // Enyme will compute the appropriate nnz from sparsification. + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * signal_size + bus_size * bus_size + 2 * size * bus_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DhDwb, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - (bus_->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - bus_->getJacobian()); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + (bus_->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - y_.size(), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index a2696fc0d..301a6621b 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -44,7 +44,6 @@ namespace GridKit using Component::yp_; using Component::wb_; using Component::h_; - using Component::J_; using Component::J_rows_buffer_; using Component::J_cols_buffer_; using Component::J_vals_buffer_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp index f5f664035..e6bd63484 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp @@ -23,64 +23,78 @@ namespace GridKit Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); if (J_rows_buffer_ == nullptr) { - // Reserve space for a dense matrix of size_*size_. + // Reserve space for the dense blocks. // Enyme will compute the appropriate nnz from sparsification. - J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; - J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto buffer_size = 2 * size * size + 2 * size * bus_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; } + nnz_ = 0; + GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DfDwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - (this->getResidualIndices()).data(), - (bus_->getVariableIndices()).data(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); GridKit::Enzyme::Sparse::DhDy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - y_.size(), - (bus_->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.data(), - yp_.data(), - wb_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - J_); + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual>::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 221b2c019..a9f0528d9 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -39,31 +39,33 @@ namespace GridKit * */ template - class SystemModel : public PhasorDynamics::Component + class SystemModel : public Component { - using PhasorDynamics::Component::gridkit_component_id_; - using PhasorDynamics::Component::size_; - using PhasorDynamics::Component::nnz_; - using PhasorDynamics::Component::time_; - using PhasorDynamics::Component::alpha_; - using PhasorDynamics::Component::y_; - using PhasorDynamics::Component::yp_; - using PhasorDynamics::Component::tag_; - using PhasorDynamics::Component::f_; - using PhasorDynamics::Component::J_; - using PhasorDynamics::Component::rel_tol_; - using PhasorDynamics::Component::abs_tol_; - using PhasorDynamics::Component::variable_indices_; - using PhasorDynamics::Component::residual_indices_; + using Component::gridkit_component_id_; + using Component::size_; + using Component::nnz_; + using Component::time_; + using Component::alpha_; + using Component::y_; + using Component::yp_; + using Component::tag_; + using Component::f_; + using Component::rel_tol_; + using Component::abs_tol_; + using Component::variable_indices_; + using Component::residual_indices_; + using Component::csr_jac_; + using Component::map_to_csr_; public: using ScalarT = scalar_type; using IdxT = index_type; using RealT = typename Model::Evaluator::RealT; using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; - using BusT = PhasorDynamics::BusBase; - using SignalT = PhasorDynamics::SignalNode; - using ComponentT = PhasorDynamics::Component; + using CooMatrixT = typename Model::Evaluator::CooMatrixT; + using BusT = BusBase; + using SignalT = SignalNode; + using ComponentT = Component; using MonitorT = Model::VariableMonitorController; SystemModel(); @@ -90,11 +92,6 @@ namespace GridKit int evaluateResidual() override; int evaluateJacobian() override; - CsrMatrixT* getCsrJacobian() const override - { - return csr_jac_; - } - void updateVariables(); void updateTime(RealT t, RealT a) override; @@ -121,9 +118,6 @@ namespace GridKit bool owns_components_{false}; - IdxT* map_to_csr_{nullptr}; - CsrMatrixT* csr_jac_{nullptr}; - /// Variable monitor std::unique_ptr monitor_; }; // class SystemModel diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 891c2c57d..6761422e2 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -395,16 +395,6 @@ namespace GridKit delete signal; } } - if (csr_jac_ != nullptr) - { - delete csr_jac_; - csr_jac_ = nullptr; - } - if (map_to_csr_ != nullptr) - { - delete[] map_to_csr_; - map_to_csr_ = nullptr; - } } /** @@ -537,6 +527,11 @@ namespace GridKit has_jacobian = has_jacobian && component->hasJacobian(); } + for (const auto& bus : buses_) + { + has_jacobian = has_jacobian && bus->hasJacobian(); + } + if (!has_jacobian) { Log::warning() << "GridKit was built with Enzyme, but some models " @@ -746,9 +741,6 @@ namespace GridKit * Finally, store bus Jacobians into the system Jacobian after all component have added their * contributions. * - * @todo split the initial assembly from updating values. This will the - * slow otherwise. - * */ template int SystemModel::evaluateJacobian() @@ -759,7 +751,7 @@ namespace GridKit bus->evaluateJacobian(); } - // Evaluate component Jacobians and update bus Jacobians + // Evaluate component Jacobians, including contribution to the bus Jacobians for (const auto& component : components_) { component->evaluateJacobian(); @@ -772,64 +764,86 @@ namespace GridKit IdxT nnz_dup = 0; for (const auto& component : components_) { - auto component_jacobian = component->getJacobian(); + auto component_jacobian = component->getCooJacobian(); - std::tuple&, std::vector&, std::vector&> component_jacobian_entries = component_jacobian.getEntries(false); - const auto [rows, columns, values] = component_jacobian_entries; - for (size_t i = 0; i < rows.size(); ++i) + if (component_jacobian != nullptr) { - ++nnz_dup; + nnz_dup += component_jacobian->getNnz(); + } + else + { + Log::warning() << "A component has returned a nullptr Jacobian.\n"; } } + for (const auto& bus : buses_) { - auto bus_jacobian = bus->getJacobian(); + auto bus_jacobian = bus->getCooJacobian(); - std::tuple&, std::vector&, std::vector&> bus_jacobian_entries = bus_jacobian.getEntries(false); - const auto [rows, columns, values] = bus_jacobian_entries; - for (size_t i = 0; i < rows.size(); ++i) + if (bus_jacobian != nullptr) { - ++nnz_dup; + nnz_dup += bus_jacobian->getNnz(); + } + else + { + Log::warning() << "A bus has returned a nullptr Jacobian.\n"; } } // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix) - IdxT* rows_dup = new IdxT[nnz_dup]; - IdxT* cols_dup = new IdxT[nnz_dup]; - RealT* vals_dup = new RealT[nnz_dup]; + IdxT* rows_dup = new IdxT[static_cast(nnz_dup)]; + IdxT* cols_dup = new IdxT[static_cast(nnz_dup)]; + RealT* vals_dup = new RealT[static_cast(nnz_dup)]; IdxT counter = 0; for (const auto& component : components_) { - auto component_jacobian = component->getJacobian(); + auto component_jacobian = component->getCooJacobian(); - std::tuple&, std::vector&, std::vector&> component_jacobian_entries = component_jacobian.getEntries(false); - const auto [rows, columns, values] = component_jacobian_entries; - for (size_t i = 0; i < rows.size(); ++i) + if (component_jacobian != nullptr) + { + const IdxT* rows = component_jacobian->getRowData(); + const IdxT* columns = component_jacobian->getColData(); + const RealT* values = component_jacobian->getValues(); + for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + { + rows_dup[counter] = rows[i]; + cols_dup[counter] = columns[i]; + vals_dup[counter] = values[i]; + counter++; + } + } + else { - rows_dup[counter] = rows[i]; - cols_dup[counter] = columns[i]; - vals_dup[counter] = values[i]; - counter++; + Log::warning() << "A component has returned a nullptr Jacobian.\n"; } } + for (const auto& bus : buses_) { - auto bus_jacobian = bus->getJacobian(); + auto bus_jacobian = bus->getCooJacobian(); - std::tuple&, std::vector&, std::vector&> bus_jacobian_entries = bus_jacobian.getEntries(false); - const auto [rows, columns, values] = bus_jacobian_entries; - for (size_t i = 0; i < rows.size(); ++i) + if (bus_jacobian != nullptr) + { + const IdxT* rows = bus_jacobian->getRowData(); + const IdxT* columns = bus_jacobian->getColData(); + const RealT* values = bus_jacobian->getValues(); + for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + { + rows_dup[counter] = rows[i]; + cols_dup[counter] = columns[i]; + vals_dup[counter] = values[i]; + counter++; + } + } + else { - rows_dup[counter] = rows[i]; - cols_dup[counter] = columns[i]; - vals_dup[counter] = values[i]; - counter++; + Log::warning() << "A bus has returned a nullptr Jacobian.\n"; } } // Build the system COO Jacobian - LinearAlgebra::CooMatrix jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); + CooMatrixT jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); // Populate CSR data with sort and deduplicate IdxT* row_ptrs = jac.getCsrRowData(); @@ -838,8 +852,8 @@ namespace GridKit nnz_ = jac.getNnz(); // Allocate cols/vals with deduplicated nnz - IdxT* cols = new IdxT[nnz_]; - RealT* vals = new RealT[nnz_]; + IdxT* cols = new IdxT[static_cast(nnz_)]; + RealT* vals = new RealT[static_cast(nnz_)]; std::copy(jac.getColData(), jac.getColData() + nnz_, cols); std::copy(jac.getValues(), jac.getValues() + nnz_, vals); @@ -851,7 +865,7 @@ namespace GridKit const IdxT* map_to_dedup = jac.getMapToDeduplicated(); // Build a mappping from original COO index to CSR index - map_to_csr_ = new IdxT[nnz_dup]; + map_to_csr_ = new IdxT[static_cast(nnz_dup)]; for (IdxT i = 0; i < nnz_dup; ++i) { map_to_csr_[map_to_sorted[i]] = map_to_dedup[i]; @@ -870,26 +884,39 @@ namespace GridKit IdxT counter = 0; for (const auto& component : components_) { - auto component_jacobian = component->getJacobian(); + auto component_jacobian = component->getCooJacobian(); - std::tuple&, std::vector&, std::vector&> component_jacobian_entries = component_jacobian.getEntries(false); - const auto [rows, columns, values] = component_jacobian_entries; - for (size_t i = 0; i < rows.size(); ++i) + if (component_jacobian != nullptr) + { + const RealT* values = component_jacobian->getValues(); + for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + { + vals[map_to_csr_[counter]] += values[i]; + counter++; + } + } + else { - vals[map_to_csr_[counter]] += values[i]; - ++counter; + Log::warning() << "A component has returned a nullptr Jacobian.\n"; } } + for (const auto& bus : buses_) { - auto bus_jacobian = bus->getJacobian(); + auto bus_jacobian = bus->getCooJacobian(); - std::tuple&, std::vector&, std::vector&> bus_jacobian_entries = bus_jacobian.getEntries(false); - const auto [rows, columns, values] = bus_jacobian_entries; - for (size_t i = 0; i < rows.size(); ++i) + if (bus_jacobian != nullptr) + { + const RealT* values = bus_jacobian->getValues(); + for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + { + vals[map_to_csr_[counter]] += values[i]; + counter++; + } + } + else { - vals[map_to_csr_[counter]] += values[i]; - ++counter; + Log::warning() << "A bus has returned a nullptr Jacobian.\n"; } } } diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index a99257454..017713b3b 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -19,7 +19,6 @@ namespace GridKit { public: using RealT = typename Model::Evaluator::RealT; - using MatrixT = typename Model::Evaluator::MatrixT; using CsrMatrixT = typename Model::Evaluator::CsrMatrixT; CircuitComponent() = default; @@ -360,16 +359,6 @@ namespace GridKit return f_; } - MatrixT& getJacobian() final - { - throw "Not Implemented"; - } - - const MatrixT& getJacobian() const final - { - throw "Not Implemented"; - } - std::vector& getIntegrand() final { return g_; diff --git a/GridKit/Model/PowerElectronics/CircuitNode.hpp b/GridKit/Model/PowerElectronics/CircuitNode.hpp index 5fd096f62..197ec9cf2 100644 --- a/GridKit/Model/PowerElectronics/CircuitNode.hpp +++ b/GridKit/Model/PowerElectronics/CircuitNode.hpp @@ -15,8 +15,7 @@ namespace GridKit template class CircuitNode : public Model::Evaluator { - using RealT = typename Model::Evaluator::RealT; - using MatrixT = typename Model::Evaluator::MatrixT; + using RealT = typename Model::Evaluator::RealT; public: CircuitNode() @@ -165,8 +164,6 @@ namespace GridKit std::vector tag_; std::vector f_; - MatrixT J_; - std::vector g_{}; std::vector param_{}; std::vector param_up_{}; @@ -312,16 +309,6 @@ namespace GridKit return f_; } - MatrixT& getJacobian() final - { - return J_; - } - - const MatrixT& getJacobian() const final - { - return J_; - } - std::vector& getIntegrand() final { return g_; diff --git a/GridKit/Model/PowerElectronics/NodeBase.hpp b/GridKit/Model/PowerElectronics/NodeBase.hpp index 9e9748529..612916ea8 100644 --- a/GridKit/Model/PowerElectronics/NodeBase.hpp +++ b/GridKit/Model/PowerElectronics/NodeBase.hpp @@ -12,8 +12,7 @@ namespace GridKit class NodeBase : public Model::Evaluator { public: - using RealT = typename Model::Evaluator::RealT; - using MatrixT = typename Model::Evaluator::MatrixT; + using RealT = typename Model::Evaluator::RealT; NodeBase(size_t n_intern, size_t n_extern) : n_intern_(n_intern), n_extern_(n_extern) @@ -90,16 +89,6 @@ namespace GridKit return tag_; } - MatrixT& getJacobian() final - { - return J_; - } - - const MatrixT& getJacobian() const final - { - return J_; - } - int setBusID(IdxT bus_id) { bus_id_ = bus_id; @@ -191,10 +180,9 @@ namespace GridKit std::vector tag_; std::vector f_; - MatrixT J_; - IdxT* J_rows_buffer_{nullptr}; - IdxT* J_cols_buffer_{nullptr}; - RealT* J_vals_buffer_{nullptr}; + IdxT* J_rows_buffer_{nullptr}; + IdxT* J_cols_buffer_{nullptr}; + RealT* J_vals_buffer_{nullptr}; RealT rtol_; RealT atol_; diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 3b44f4ec6..84b543168 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -16,48 +16,6 @@ namespace GridKit { - /** - * Writes a vector to a file in Matrix Market format - * - * @param vec The vector to write - * @param filename The name of the output file - * @param header Additional header information/comments - * @return true if the write was successful, false otherwise - */ - template - void writeVectorToMatrixMarket(const std::vector& vec, const std::string& filename, const std::string& header) - { - std::ofstream outFile(filename); - - if (!outFile.is_open()) - { - std::cerr << "Error: Could not open file " << filename << " for writing." << std::endl; - return; - } - - // Uncomment to write Matrix Market header - // outFile << "%%MatrixMarket vector array real general" << std::endl; - - // Write additional header information as comments - if (!header.empty()) - { - outFile << "% " << header << std::endl; - } - - // Write the vector size - outFile << vec.size() << std::endl; - - // Write the vector elements - outFile << std::scientific << std::setprecision(16); - for (const auto& val : vec) - { - outFile << val << std::endl; - } - - outFile.close(); - return; - } - template class PowerElectronicsModel : public CircuitComponent { @@ -522,17 +480,6 @@ namespace GridKit alpha_ = a; } - /** - * @brief print the system residual in COO format - * - * @param[in] filename - * @param[in] title - */ - void printResidualMatrixMarket(std::string filename, std::string title) - { - writeVectorToMatrixMarket(f_, filename, title); - } - CsrMatrixT* getCsrJacobian() const override { return csr_jac_; diff --git a/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp b/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp index 2457b16b9..66f31d7a6 100644 --- a/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp +++ b/GridKit/Model/PowerFlow/ModelEvaluatorImpl.hpp @@ -16,8 +16,7 @@ namespace GridKit class ModelEvaluatorImpl : public Model::Evaluator { public: - using RealT = typename Model::Evaluator::RealT; - using MatrixT = typename Model::Evaluator::MatrixT; + using RealT = typename Model::Evaluator::RealT; ModelEvaluatorImpl() : size_(0), @@ -38,7 +37,6 @@ namespace GridKit ypB_(static_cast(size_)), fB_(static_cast(size_)), gB_(static_cast(size_opt_)), - jac_(MatrixT()), param_(static_cast(size_opt_)), param_up_(static_cast(size_opt_)), param_lo_(static_cast(size_opt_)) @@ -178,16 +176,6 @@ namespace GridKit return f_; } - MatrixT& getJacobian() - { - return jac_; - } - - const MatrixT& getJacobian() const - { - return jac_; - } - std::vector& getIntegrand() { return g_; @@ -241,8 +229,6 @@ namespace GridKit std::vector fB_; std::vector gB_; - MatrixT jac_; - std::vector param_; std::vector param_up_; std::vector param_lo_; diff --git a/GridKit/Utilities/CMakeLists.txt b/GridKit/Utilities/CMakeLists.txt index 478eb11e8..7a94598a6 100644 --- a/GridKit/Utilities/CMakeLists.txt +++ b/GridKit/Utilities/CMakeLists.txt @@ -13,5 +13,5 @@ install( FILES Colors.hpp Errors.hpp FileIO.hpp - MapFromCOO.hpp + MapFromCsr.hpp DESTINATION include/GridKit/Utilities) diff --git a/GridKit/Utilities/MapFromCOO.hpp b/GridKit/Utilities/MapFromCOO.hpp deleted file mode 100644 index d4fa91c29..000000000 --- a/GridKit/Utilities/MapFromCOO.hpp +++ /dev/null @@ -1,37 +0,0 @@ - -/** - * @file MapFromCOO.hpp - * - * @author Nicholson Koukpaizan , ORNL - * @todo This should go away once we settle on a sparse format - * - */ -#pragma once - -#include -#include - -namespace GridKit -{ - namespace Testing - { - template - std::vector MapFromCOO(LinearAlgebra::COO_Matrix matrix) - { - std::tuple&, std::vector&, std::vector&> matrix_entries = matrix.getEntries(); - const auto [rows, columns, values] = matrix_entries; - - std::tuple matrix_dimensions = matrix.getDimensions(); - const auto [n_rows, n_columns] = matrix_dimensions; - - std::vector dependencies(n_rows); - - for (IdxT i = 0; i < rows.size(); ++i) - { - dependencies[rows[i]].insert(std::make_pair(columns[i], values[i])); - } - - return dependencies; - } - } // namespace Testing -} // namespace GridKit diff --git a/buildsystem/spack_repo/gridkit/packages/gridkit/package.py b/buildsystem/spack_repo/gridkit/packages/gridkit/package.py index c0f5a6491..f97c416c8 100644 --- a/buildsystem/spack_repo/gridkit/packages/gridkit/package.py +++ b/buildsystem/spack_repo/gridkit/packages/gridkit/package.py @@ -39,11 +39,11 @@ def cmake_args(self): args.extend( [ - self.define_from_variant("GRIDKIT_ENABLE_IPOPT", "ipopt"), - self.define_from_variant("GRIDKIT_ENABLE_SUNDIALS", "sundials"), - self.define_from_variant("GRIDKIT_ENABLE_ENZYME", "enzyme"), - self.define_from_variant("GRIDKIT_ENABLE_ASAN", "asan"), - self.define_from_variant("GRIDKIT_ENABLE_UBSAN", "ubsan"), + self.define_from_variant("GridKit_ENABLE_IPOPT", "ipopt"), + self.define_from_variant("GridKit_ENABLE_SUNDIALS", "sundials"), + self.define_from_variant("GridKit_ENABLE_ENZYME", "enzyme"), + self.define_from_variant("GridKit_ENABLE_ASAN", "asan"), + self.define_from_variant("GridKit_ENABLE_UBSAN", "ubsan"), ] ) diff --git a/examples/Enzyme/Library/Vector/EnzymeVector.cpp b/examples/Enzyme/Library/Vector/EnzymeVector.cpp index dd8a82b31..c16b1e6cf 100644 --- a/examples/Enzyme/Library/Vector/EnzymeVector.cpp +++ b/examples/Enzyme/Library/Vector/EnzymeVector.cpp @@ -77,8 +77,8 @@ int main() } if (verbose) { - jac.printMatrix("Autodiff Jacobian"); - jac_ref.printMatrix("Reference Jacobian"); + jac.print("Autodiff Jacobian"); + jac_ref.print("Reference Jacobian"); } std::cout << "Status: " << fail << "\n"; diff --git a/examples/Enzyme/PowerElectronics/main.cpp b/examples/Enzyme/PowerElectronics/main.cpp index 246cb92b6..351378d51 100644 --- a/examples/Enzyme/PowerElectronics/main.cpp +++ b/examples/Enzyme/PowerElectronics/main.cpp @@ -17,7 +17,6 @@ */ using DenseMatrix = GridKit::LinearAlgebra::DenseMatrix; -using SparseMatrix = GridKit::LinearAlgebra::COO_Matrix; using DG = GridKit::DistributedGenerator; using DGParameters = GridKit::DistributedGeneratorParameters; @@ -222,8 +221,8 @@ int main() } if (verbose) { - jac_autodiff.printMatrix("Autodiff Jacobian"); - jac_ref_dense.printMatrix("Reference Jacobian"); + jac_autodiff.print("Autodiff Jacobian"); + jac_ref_dense.print("Reference Jacobian"); } std::cout << "Status: " << fail << "\n"; diff --git a/examples/Enzyme/Standalone/EnzymeSparse.cpp b/examples/Enzyme/Standalone/EnzymeSparse.cpp index 1b42bb635..4f1818c4e 100644 --- a/examples/Enzyme/Standalone/EnzymeSparse.cpp +++ b/examples/Enzyme/Standalone/EnzymeSparse.cpp @@ -84,6 +84,7 @@ __attribute__((noinline)) SparseMatrix* jac_f(size_t N, ScalarT* input) ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, i, + 1.0, index_maps, index_maps, rows_buffer, diff --git a/examples/Enzyme/Standalone/EnzymeVector.cpp b/examples/Enzyme/Standalone/EnzymeVector.cpp index 7b7b13815..f25239028 100644 --- a/examples/Enzyme/Standalone/EnzymeVector.cpp +++ b/examples/Enzyme/Standalone/EnzymeVector.cpp @@ -123,8 +123,8 @@ int main() } if (verbose) { - dsq.printMatrix("Autodiff Jacobian"); - dsq_ref.printMatrix("Reference Jacobian"); + dsq.print("Autodiff Jacobian"); + dsq_ref.print("Reference Jacobian"); } std::cout << "Status: " << fail << "\n"; return fail; diff --git a/examples/LinearAlgebra/CMakeLists.txt b/examples/LinearAlgebra/CMakeLists.txt index aaad9baf7..6c3edab8c 100644 --- a/examples/LinearAlgebra/CMakeLists.txt +++ b/examples/LinearAlgebra/CMakeLists.txt @@ -1,2 +1 @@ -add_subdirectory(SparseTest) add_subdirectory(DenseTest) diff --git a/examples/LinearAlgebra/DenseTest/DenseTest.cpp b/examples/LinearAlgebra/DenseTest/DenseTest.cpp index 0efd5e94a..22a537d42 100644 --- a/examples/LinearAlgebra/DenseTest/DenseTest.cpp +++ b/examples/LinearAlgebra/DenseTest/DenseTest.cpp @@ -5,10 +5,9 @@ int main() { int fail = 0; - size_t m = 4; - size_t n = 4; - GridKit::LinearAlgebra::DenseMatrix A = - GridKit::LinearAlgebra::DenseMatrix(m, n); + size_t m = 4; + size_t n = 4; + auto A = GridKit::LinearAlgebra::DenseMatrix(m, n); double val = 0.0; for (size_t j = 0; j < n; ++j) @@ -19,8 +18,7 @@ int main() val += 1.0; } } - A.printMatrix("Dense matrix test output"); - (A.getValuesCOO())->printMatrix(); + A.print("Dense matrix test output"); return fail; } diff --git a/examples/LinearAlgebra/SparseTest/CMakeLists.txt b/examples/LinearAlgebra/SparseTest/CMakeLists.txt deleted file mode 100644 index abda11b7f..000000000 --- a/examples/LinearAlgebra/SparseTest/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -add_executable(spmattest SparseTest.cpp) -target_link_libraries( - spmattest GridKit::sparse_matrix) - -add_test(NAME SparseMatrixTest COMMAND spmattest) -install(TARGETS spmattest RUNTIME DESTINATION bin) diff --git a/examples/LinearAlgebra/SparseTest/SparseTest.cpp b/examples/LinearAlgebra/SparseTest/SparseTest.cpp deleted file mode 100644 index 2e5224c7b..000000000 --- a/examples/LinearAlgebra/SparseTest/SparseTest.cpp +++ /dev/null @@ -1,102 +0,0 @@ - - -#include -#include -#include -#include -#include -#include -#include - -#include - -int main() -{ - std::vector val{0.1, 0.2, 0.3, 0.4}; - std::vector x{2, 1, 3, 1}; - std::vector y{1, 3, 2, 2}; - size_t n = 4; - size_t m = 4; - - GridKit::LinearAlgebra::COO_Matrix A = GridKit::LinearAlgebra::COO_Matrix(x, y, val, m, n); - - std::vector valn(4); - std::vector xn(4); - std::vector yn(4); - - std::tie(xn, yn, valn) = A.getEntries(); - - for (size_t i = 0; i < valn.size(); i++) - { - std::cout << valn[i] << "\n"; - } - - std::cout << "A:\n"; - A.printMatrix(); - - std::vector val2{0.5, 0.6, 0.7, 0.8, 1.0}; - std::vector x2{0, 2, 0, 2, 1}; - std::vector y2{3, 3, 2, 2, 3}; - GridKit::LinearAlgebra::COO_Matrix B = GridKit::LinearAlgebra::COO_Matrix(x2, y2, val2, m, n); - - std::cout << "B:\n"; - B.printMatrix(); - - A.axpy(2.0, B); - - std::cout << "A + 2B:\n"; - A.printMatrix(); - - std::vector r; - std::vector c; - std::vector v; - std::tie(r, c, v) = A.setDataToCSR(); - - for (size_t i = 0; i < r.size() - 1; i++) - { - std::cout << r[i] << std::endl; - size_t rdiff = r[i + 1] - r[i]; - for (size_t j = 0; j < rdiff; j++) - { - std::cout << c[j + r[i]] << ", " << v[j + r[i]] << std::endl; - } - } - std::cout << r[r.size() - 1] << std::endl; - - // Basic Verification test - std::vector rtest = {0, 2, 4, 7, 8}; - std::vector ctest = {2, 3, 2, 3, 1, 2, 3, 2}; - std::vector valtest = {1.4, 1.0, 0.4, 2.2, 0.1, 1.6, 1.2, 0.3}; - - assert(rtest.size() == r.size()); - assert(ctest.size() == c.size()); - assert(valtest.size() == v.size()); - - int failval = 0; - for (size_t i = 0; i < rtest.size(); i++) - { - if (r[i] != rtest[i]) - { - failval--; - } - } - for (size_t i = 0; i < ctest.size(); i++) - { - double vdiff = v[i] - valtest[i]; - if (c[i] != ctest[i] || -1e-14 > vdiff || vdiff > 1e-14) - { - failval--; - } - } - - if (failval == 0) - { - std::cout << "Success!" << std::endl; - } - else - { - std::cout << "Failed!" << std::endl; - } - - return failval; -} diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp index 07246e83a..2cb9220ac 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp +++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp @@ -261,8 +261,6 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) sys_model->initialize(); sys_model->evaluateResidual(); - // print the residual in matrix market format - sys_model->printResidualMatrixMarket("ScaleMicrogrid_Residual_N" + std::to_string(Nsize) + ".mtx", "ScaleMicrogrid Residual N" + std::to_string(Nsize)); std::vector& fres = sys_model->getResidual(); if (debug_output) { @@ -276,8 +274,6 @@ int test(index_type Nsize, real_type error_tol, bool debug_output) sys_model->updateTime(0.0, 1.0e-8); sys_model->evaluateJacobian(); - // sys_model->printJacobianMatrixMarket("ScaleMicrogrid_Jacobian_N" + std::to_string(Nsize) + ".mtx", "ScaleMicrogrid Jacobian N" + std::to_string(Nsize)); - // print the jacobian in matrix market format if (debug_output) { diff --git a/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt b/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt index a7075f92f..16e3e85ab 100644 --- a/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt +++ b/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt @@ -2,7 +2,7 @@ add_executable(test_sparse_csr runSparseCsrTests.cpp) target_link_libraries( test_sparse_csr GridKit::sparse_matrix GridKit::testing) -add_test(NAME SparseTest COMMAND $) +add_test(NAME SparseCsrTest COMMAND $) add_executable(test_sparse_coo runSparseCooTests.cpp) target_link_libraries( diff --git a/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp new file mode 100644 index 000000000..3e5ea1083 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/BusFaultTests.hpp @@ -0,0 +1,244 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + + template + class BusFaultTests + { + private: + using RealT = typename PhasorDynamics::Component::RealT; + + public: + BusFaultTests() = default; + ~BusFaultTests() = default; + + TestOutcome constructor() + { + TestStatus success = true; + + auto* bus = new PhasorDynamics::Bus(1.0, 0.0); + + PhasorDynamics::Component* fault = + new PhasorDynamics::BusFault(bus); + + success *= (fault != nullptr); + + if (fault) + { + delete fault; + } + delete bus; + + return success.report(__func__); + } + + /** + * Verifies the residual evaluates to zero for the initial conditions + */ + TestOutcome zeroInitialResidual(bool status = false) + { + TestStatus success = true; + + ScalarT Vr1{1.0}; ///< Bus real voltage + ScalarT Vi1{1.0}; ///< Bus imaginary voltage + + PhasorDynamics::Bus bus(Vr1, Vi1); + PhasorDynamics::BusFault fault(&bus, 0.0, 1e-3, status); + bus.allocate(); + bus.initialize(); + fault.allocate(); + fault.initialize(); + fault.evaluateResidual(); + std::vector res = fault.getResidual(); + + for (size_t i = 0; i < res.size(); ++i) + { + if (!isEqual(res[i], 0.0)) + { + std::cout << "Incorrect result: " + << fault.yp()[i] << " != 0\n"; + success = false; + break; + } + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /** + * A test case to verify Jacobian values + */ + TestOutcome jacobian(bool status = false) + { + TestStatus success = true; + + RealT R = 0.0; + RealT X = 1e-3; + + // Jacobian via DependencyTracking + auto dependency_tracking_jacobian = DependencyTrackingJacobian(R, X, status); + + // Jacobian via Enzyme + auto enzyme_jacobian = EnzymeJacobian(R, X, status); + + /// Compare DependencyTracking dependencies to Enzyme's + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) + { + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i])); + } + + return success.report(__func__); + } + + private: + std::vector DependencyTrackingJacobian( + const RealT R, const RealT X, const bool status) + { + DependencyTracking::Variable Vr1{1.0}; ///< Bus-1 real voltage + DependencyTracking::Variable Vi1{1.0}; ///< Bus-1 imaginary voltage + + PhasorDynamics::Bus bus(Vr1, Vi1); + PhasorDynamics::BusFault fault(&bus, R, X, status); + + bus.allocate(); + fault.allocate(); + + // Get d/dy + bus.initialize(); + fault.initialize(); + + for (size_t i = 0; i < fault.size(); ++i) + { + fault.y()[i].setVariableNumber(i); ///< fault independent variables + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + fault.size()); // Bus independent variables + } + + bus.evaluateResidual(); + fault.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking + ///< the dependencies + std::vector residual_y = fault.getResidual(); + + // Get d/dy' + bus.initialize(); + fault.initialize(); + + for (size_t i = 0; i < fault.size(); ++i) + { + fault.yp()[i].setVariableNumber(i); ///< fault independent variables + } + + bus.evaluateResidual(); + fault.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking + ///< the dependencies + std::vector residual_yp = fault.getResidual(); + + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) + { + 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"; + } + + // 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( + const RealT R, const RealT X, const bool status) + { + ScalarT Vr1{1.0}; ///< Bus-1 real voltage + ScalarT Vi1{1.0}; ///< Bus-1 imaginary voltage + + PhasorDynamics::Bus bus(Vr1, Vi1); + PhasorDynamics::BusFault fault(&bus, R, X, status); + + bus.allocate(); + fault.allocate(); + + bus.initialize(); + fault.initialize(); + + fault.updateTime(0.0, 1.0); + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + fault.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + fault.size()); // Reset bus residual indices + } + + bus.evaluateResidual(); + fault.evaluateResidual(); + + bus.evaluateJacobian(); + fault.evaluateJacobian(); + fault.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = fault.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: BusFault Jacobian\n"; + model_jacobian->print(); + + return GridKit::Testing::MapFromCsr(model_jacobian); + } +#endif + + }; // class BusFaultTests + + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index bac53ec7b..15fde1d52 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -7,6 +7,16 @@ add_executable(test_phasor_bus runBusTests.cpp) target_link_libraries( test_phasor_bus GridKit::phasor_dynamics_bus GridKit::testing) +add_executable(test_phasor_bus_fault runBusFaultTests.cpp) +target_link_libraries( + test_phasor_bus_fault + GridKit::definitions + GridKit::phasor_dynamics_bus_fault + GridKit::phasor_dynamics_bus_fault_dependency_tracking + GridKit::phasor_dynamics_bus + GridKit::phasor_dynamics_bus_dependency_tracking + GridKit::testing) + add_executable(test_phasor_bustosignaladapter runBusToSignalAdapterTests.cpp) target_link_libraries( test_phasor_bustosignaladapter @@ -131,6 +141,7 @@ target_link_libraries( GridKit::testing) add_test(NAME PhasorDynamicsBusTest COMMAND test_phasor_bus) +add_test(NAME PhasorDynamicsBusFaultTest COMMAND test_phasor_bus_fault) add_test(NAME PhasorDynamicsBusToSignalAdapterTest COMMAND test_phasor_bustosignaladapter) add_test(NAME PhasorDynamicsBranchTest COMMAND test_phasor_branch) add_test(NAME PhasorDynamicsGenrouTest COMMAND test_phasor_genrou) @@ -150,6 +161,7 @@ add_test(NAME PhasorDynamicsSystemSingleComponentTest COMMAND test_phasor_system install( TARGETS test_phasor_bus + test_phasor_bus_fault test_phasor_bustosignaladapter test_phasor_branch test_phasor_load diff --git a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp index e514d7f3a..fa440b3b4 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp @@ -12,7 +12,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -220,11 +220,12 @@ namespace GridKit bus.evaluateJacobian(); exciter.evaluateJacobian(); - GridKit::LinearAlgebra::COO_Matrix& model_jacobian = exciter.getJacobian(); - model_jacobian.deduplicate(); - model_jacobian.printMatrix("Model Jacobian"); + exciter.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = exciter.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Ieeet1 Jacobian\n"; + model_jacobian->print(); - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif diff --git a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp index 085ba7416..d6e2d7369 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -13,7 +14,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -436,11 +437,12 @@ namespace GridKit bus.evaluateJacobian(); exciter.evaluateJacobian(); - GridKit::LinearAlgebra::COO_Matrix& model_jacobian = exciter.getJacobian(); - model_jacobian.deduplicate(); - model_jacobian.printMatrix("Model Jacobian"); + exciter.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = exciter.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: SexsPti Jacobian\n"; + model_jacobian->print(); - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index d956ba9f3..73b9cd3dc 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -20,7 +20,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -445,11 +445,12 @@ namespace GridKit bus.evaluateJacobian(); gen.evaluateJacobian(); - GridKit::LinearAlgebra::COO_Matrix& model_jacobian = gen.getJacobian(); - model_jacobian.deduplicate(); - model_jacobian.printMatrix("Model Jacobian"); + gen.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: GenClassical Jacobian\n"; + model_jacobian->print(); - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 62bf255b2..84a888b7c 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -11,7 +11,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -486,11 +486,12 @@ namespace GridKit bus.evaluateJacobian(); gen.evaluateJacobian(); - GridKit::LinearAlgebra::COO_Matrix& model_jacobian = gen.getJacobian(); - model_jacobian.deduplicate(); - model_jacobian.printMatrix("Model Jacobian"); + gen.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Genrou Jacobian\n"; + model_jacobian->print(); - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif }; // class GenrouTest diff --git a/tests/UnitTests/PhasorDynamics/GensalTests.hpp b/tests/UnitTests/PhasorDynamics/GensalTests.hpp index a382e3172..a4b3982d8 100644 --- a/tests/UnitTests/PhasorDynamics/GensalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GensalTests.hpp @@ -12,7 +12,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -384,6 +384,17 @@ namespace GridKit gen.evaluateResidual(); std::vector residual_yp = gen.getResidual(); + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) + { + 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"; + } + std::vector dependencies(residual_y.size()); for (IdxT i = 0; i < residual_y.size(); ++i) { @@ -448,10 +459,12 @@ namespace GridKit bus.evaluateJacobian(); gen.evaluateJacobian(); - GridKit::LinearAlgebra::COO_Matrix& model_jacobian = gen.getJacobian(); - model_jacobian.deduplicate(); + gen.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gen.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Gensal Jacobian\n"; + model_jacobian->print(); - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif }; // class GensalTests diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 3d9d40eb0..2e4004975 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -14,7 +14,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -416,11 +416,12 @@ namespace GridKit gov.evaluateResidual(); gov.evaluateJacobian(); - GridKit::LinearAlgebra::COO_Matrix model_jacobian = gov.getJacobian(); - model_jacobian.deduplicate(); - model_jacobian.printMatrix("Model Jacobian"); + gov.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = gov.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Tgov1 Jacobian\n"; + model_jacobian->print(); - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif }; // class GovernorTgov1Tests diff --git a/tests/UnitTests/PhasorDynamics/LoadTests.hpp b/tests/UnitTests/PhasorDynamics/LoadTests.hpp index c50f60a16..7229d48fa 100644 --- a/tests/UnitTests/PhasorDynamics/LoadTests.hpp +++ b/tests/UnitTests/PhasorDynamics/LoadTests.hpp @@ -10,7 +10,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -123,7 +123,7 @@ namespace GridKit } #ifdef GRIDKIT_ENABLE_ENZYME - TestOutcome enzyme_jacobian() + TestOutcome enzymeJacobian() { TestStatus success = true; @@ -150,13 +150,14 @@ namespace GridKit bus.evaluateJacobian(); load.evaluateJacobian(); - - GridKit::LinearAlgebra::COO_Matrix model_jacobian = load.getJacobian(); - model_jacobian.printMatrix("Model Jacobian"); + load.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = load.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Load Jacobian\n"; + model_jacobian->print(); /// Compare model Jacobian wih dependencies computed analytically std::vector ref = analyticalJacobian(R, X); - std::vector model_dependencies = GridKit::Testing::MapFromCOO(model_jacobian); + std::vector model_dependencies = GridKit::Testing::MapFromCsr(model_jacobian); for (size_t i = 0; i < ref.size(); ++i) { success *= (GridKit::Testing::isEqual(model_dependencies[i], ref[i])); diff --git a/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp b/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp index 15b68d535..1c374c04c 100644 --- a/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp +++ b/tests/UnitTests/PhasorDynamics/LoadZIPTests.hpp @@ -10,7 +10,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -237,11 +237,12 @@ namespace GridKit bus.evaluateJacobian(); load.evaluateJacobian(); - GridKit::LinearAlgebra::COO_Matrix& model_jacobian = load.getJacobian(); - model_jacobian.deduplicate(); - model_jacobian.printMatrix("Model Jacobian"); + load.constructCsr(); + GridKit::LinearAlgebra::CsrMatrix* model_jacobian = load.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: LoadZIP Jacobian\n"; + model_jacobian->print(); - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif }; // class LoadZIPTests diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 0c1c77d66..15b64876f 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -11,7 +11,7 @@ #include #include #include -#include +#include namespace GridKit { @@ -307,12 +307,12 @@ namespace GridKit stab.evaluateResidual(); stab.evaluateJacobian(); + stab.constructCsr(); + auto model_jacobian = stab.getCsrJacobian(); + std::cout << "Sparse Csr Matrix: Ieeest Jacobian\n"; + model_jacobian->print(); - auto model_jacobian = stab.getJacobian(); - model_jacobian.deduplicate(); - model_jacobian.printMatrix("IEEEST Jacobian"); - - return GridKit::Testing::MapFromCOO(model_jacobian); + return GridKit::Testing::MapFromCsr(model_jacobian); } #endif diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index c32df9f58..985bae0f6 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -85,7 +85,6 @@ namespace GridKit success *= system->initialize() == 0; success *= system->evaluateResidual() == 0; success *= system->evaluateJacobian() == 0; - success *= system->size() == 0; success *= system->size() == fault.size(); delete system; diff --git a/tests/UnitTests/PhasorDynamics/runBusFaultTests.cpp b/tests/UnitTests/PhasorDynamics/runBusFaultTests.cpp new file mode 100644 index 000000000..5a3688fc5 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runBusFaultTests.cpp @@ -0,0 +1,18 @@ +#include "BusFaultTests.hpp" + +int main() +{ + using namespace GridKit; + using namespace GridKit::Testing; + + GridKit::Testing::TestingResults result; + GridKit::Testing::BusFaultTests test; + + result += test.constructor(); + result += test.zeroInitialResidual(true); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(true); +#endif + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runLoadTests.cpp b/tests/UnitTests/PhasorDynamics/runLoadTests.cpp index 1dac30e27..e020a530e 100644 --- a/tests/UnitTests/PhasorDynamics/runLoadTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runLoadTests.cpp @@ -12,7 +12,7 @@ int main() result += test.residual(); result += test.jacobian(); #ifdef GRIDKIT_ENABLE_ENZYME - result += test.enzyme_jacobian(); + result += test.enzymeJacobian(); #endif return result.summary(); diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index e3592bcf2..63196214e 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -200,14 +200,9 @@ namespace GridKit return f_; } - GridKit::LinearAlgebra::COO_Matrix& getJacobian() override + GridKit::LinearAlgebra::CsrMatrix* getCsrJacobian() const override { - return jac_; - } - - const GridKit::LinearAlgebra::COO_Matrix& getJacobian() const override - { - return jac_; + return csr_jac_; } std::vector& getIntegrand() override @@ -257,7 +252,7 @@ namespace GridKit std::vector fB_; std::vector gB_; - GridKit::LinearAlgebra::COO_Matrix jac_; + GridKit::LinearAlgebra::CsrMatrix* csr_jac_; std::vector param_; std::vector param_up_;