From cf71e3b861e0ad0eed25b28087457e37ced924ea Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 24 Nov 2025 14:42:00 -0500 Subject: [PATCH 01/14] Functional PhasorDynamics simulations with sparse Jacobians. --- .../Enzyme/SparseWrapper.hpp | 353 +++++++++++++++++- .../LinearAlgebra/SparseMatrix/COO_Matrix.hpp | 8 - .../Model/PhasorDynamics/Branch/Branch.cpp | 4 +- .../Model/PhasorDynamics/Branch/Branch.hpp | 4 - .../Branch/BranchDependencyTracking.cpp | 4 +- .../PhasorDynamics/Branch/BranchEnzyme.cpp | 10 +- GridKit/Model/PhasorDynamics/Bus/Bus.cpp | 4 +- .../Bus/BusDependencyTracking.cpp | 4 +- .../Model/PhasorDynamics/Bus/BusEnzyme.cpp | 7 +- GridKit/Model/PhasorDynamics/BusBase.hpp | 3 - .../PhasorDynamics/BusFault/BusFault.cpp | 4 +- .../BusFault/BusFaultDependencyTracking.cpp | 4 +- .../BusFault/BusFaultEnzyme.cpp | 4 +- GridKit/Model/PhasorDynamics/Component.hpp | 16 +- .../Exciter/IEEET1/CMakeLists.txt | 36 +- .../PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp | 14 + .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 14 +- .../IEEET1/Ieeet1DependencyTracking.cpp | 14 + .../Exciter/IEEET1/Ieeet1Enzyme.cpp | 88 +++++ .../Exciter/IEEET1/Ieeet1Impl.hpp | 117 +++--- .../PhasorDynamics/Governor/Tgov1/Tgov1.cpp | 2 +- .../PhasorDynamics/Governor/Tgov1/Tgov1.hpp | 6 +- .../Tgov1/Tgov1DependencyTracking.cpp | 2 +- .../Governor/Tgov1/Tgov1Enzyme.cpp | 11 +- .../Governor/Tgov1/Tgov1Impl.hpp | 6 +- GridKit/Model/PhasorDynamics/Load/Load.cpp | 4 +- GridKit/Model/PhasorDynamics/Load/Load.hpp | 4 - .../Load/LoadDependencyTracking.cpp | 4 +- .../Model/PhasorDynamics/Load/LoadEnzyme.cpp | 4 +- .../PhasorDynamics/SignalNode/SignalNode.hpp | 5 +- .../SynchronousMachine/GENROUwS/Genrou.cpp | 4 +- .../SynchronousMachine/GENROUwS/Genrou.hpp | 6 +- .../GENROUwS/GenrouDependencyTracking.cpp | 4 +- .../GENROUwS/GenrouEnzyme.cpp | 36 +- .../GENROUwS/GenrouImpl.hpp | 8 +- .../GenClassical/GenClassical.cpp | 4 +- .../GenClassical/GenClassical.hpp | 4 - .../GenClassicalDependencyTracking.cpp | 4 +- .../GenClassical/GenClassicalEnzyme.cpp | 33 +- GridKit/Model/PhasorDynamics/SystemModel.hpp | 16 +- GridKit/Solver/Dynamic/Ida.cpp | 5 +- .../gridkit/packages/enzyme/package.py | 1 + .../Small/TenGen/Classical/CMakeLists.txt | 2 +- .../TenGen/Classical/TenGenClassical.cpp | 9 +- .../Small/TenGen/Genrou/CMakeLists.txt | 2 +- .../Small/TenGen/Genrou/TenGenGenrou.cpp | 9 +- .../Tiny/ThreeBus/Basic/ThreeBusBasic.cpp | 45 +-- .../Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp | 29 +- .../ThreeBus/Classical/ThreeBusClassical.cpp | 47 +-- .../Classical/ThreeBusClassicalJson.cpp | 29 +- .../Tiny/TwoBus/Basic/TwoBusBasic.cpp | 12 +- .../Tiny/TwoBus/Basic/TwoBusBasicJson.cpp | 8 +- .../Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp | 17 +- .../Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp | 13 +- .../Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp | 8 +- .../Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp | 8 +- .../PhasorDynamics/GenClassicalTests.hpp | 101 ++++- .../UnitTests/PhasorDynamics/GenrouTests.hpp | 100 ++++- .../PhasorDynamics/GovernorTgov1Tests.hpp | 81 +++- .../UnitTests/PhasorDynamics/SystemTests.hpp | 27 +- 60 files changed, 1106 insertions(+), 326 deletions(-) create mode 100644 GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp diff --git a/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp b/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp index 3bba002f0..6beb0a565 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp @@ -43,6 +43,17 @@ namespace GridKit BusResidual22 //< Special case for branches that are connected to two buses }; + /** + * @brief Model independent variable parameter keys + * + */ + enum class IndependentVariables + { + Internal, + Bus, + Signal + }; + /** * @brief Template definition for wrapper around residual methods inside model classes * @@ -327,6 +338,7 @@ namespace GridKit * @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, @@ -337,12 +349,14 @@ namespace GridKit ScalarT* y, ScalarT* yp, ScalarT* wb, + RealT alpha, MatrixT& jac) { if (n_res > 0 && n_var > 0) { + // df/dy std::vector> triplets; - std::vector elementary_v(n_res); + std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage @@ -380,6 +394,126 @@ namespace GridKit valtemp.push_back(tup.val); } jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + + // df/dy' + triplets.clear(); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) Wrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_dup, + yp, + output, + enzyme_const, + wb, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + ctemp.clear(); + rtemp.clear(); + valtemp.clear(); + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + + /** + * @brief Enzyme automatic differentiation Jacobian evaluator : df/dwb + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct df_dwb + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) Wrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_const, + yp, + enzyme_dup, + wb, + output, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes } } }; @@ -408,6 +542,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 */ static void eval(ModelT* model, @@ -419,12 +554,14 @@ namespace GridKit ScalarT* yp, ScalarT* wb, ScalarT* ws, + RealT alpha, MatrixT& jac) { if (n_res > 0 && n_var > 0) { + // df/dy std::vector> triplets; - std::vector elementary_v(n_res); + std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage @@ -464,6 +601,132 @@ namespace GridKit valtemp.push_back(tup.val); } jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + + // df/dy' + triplets.clear(); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) Wrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_dup, + yp, + output, + enzyme_const, + wb, + enzyme_const, + ws, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + ctemp.clear(); + rtemp.clear(); + valtemp.clear(); + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: df/dwb with signal variables + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct df_dwbWithSignal + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in] ws - Signal variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) Wrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_const, + yp, + enzyme_dup, + wb, + output, + enzyme_const, + ws, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes } } }; @@ -508,7 +771,7 @@ namespace GridKit if (n_res > 0 && n_var > 0) { std::vector> triplets; - std::vector elementary_v(n_res); + std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage @@ -595,7 +858,7 @@ namespace GridKit if (n_res > 0 && n_var > 0) { std::vector> triplets; - std::vector elementary_v(n_res); + std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage @@ -637,6 +900,86 @@ namespace GridKit } }; + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: dh/dy + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct dh_dy + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) Wrapper::eval, + enzyme_const, + model, + enzyme_dup, + y, + output, + enzyme_const, + yp, + enzyme_const, + wb, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + /** * @brief Enzyme automatic differentiation Jacobian evaluator: Branch Jacobian * @@ -675,7 +1018,7 @@ namespace GridKit if (n_res > 0 && n_var > 0) { std::vector> triplets; - std::vector elementary_v(n_res); + std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { // Sparse storage diff --git a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp index 80f80ecb1..1cd2f6e9f 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp @@ -292,11 +292,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, COO_Matrix& a) { - if (alpha == 0) - { - return; - } - if (!this->sorted_) { this->sortSparse(); @@ -372,9 +367,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, std::vector r, std::vector c, std::vector v) { - if (alpha == 0) - return; - if (!this->sorted_) { this->sortSparse(); diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.cpp b/GridKit/Model/PhasorDynamics/Branch/Branch.cpp index cc1dd85fa..527216e92 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.cpp @@ -20,8 +20,8 @@ namespace GridKit template int Branch::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Branch..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Branch..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index c917e4311..a45391ebb 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -76,10 +76,6 @@ namespace GridKit return 0; } - virtual void updateTime(RealT /* t */, RealT /* a */) override - { - } - void setR(RealT R) { R_ = R; diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp index 630df803d..4c5780ffb 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp @@ -20,8 +20,8 @@ namespace GridKit template int Branch::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Branch..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Branch..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp index d66948c42..8ba9da1c4 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp @@ -22,8 +22,10 @@ namespace GridKit template int Branch::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Branch..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + //std::cout << "Evaluate Jacobian for Branch..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); // Bus 1 diagonal Jacobian block owned by the bus GridKit::Enzyme::Sparse::BusJacobian, @@ -67,8 +69,6 @@ namespace GridKit (bus2_->y()).data(), J_); - J_.printMatrix("Branch after BusJacobian12 call"); - // Off-diagonal Jacobian block (Bus1 variables) owned by the branch GridKit::Enzyme::Sparse::BranchJacobian, GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21, @@ -83,8 +83,6 @@ namespace GridKit (bus1_->y()).data(), J_); - J_.printMatrix("Branch after BusJacobian21 call"); - return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.cpp b/GridKit/Model/PhasorDynamics/Bus/Bus.cpp index 7e633dae7..261d11038 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Bus::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Bus..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Bus..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp index 3d51ebd16..c010cb816 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Bus::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Bus..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Bus..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp index 2e3f5a372..5618ffb39 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp @@ -24,14 +24,15 @@ namespace GridKit template int Bus::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Bus..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + //std::cout << "Evaluate Jacobian for Bus..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); std::vector rtemp = {residual_indices_.at(0), residual_indices_.at(0), residual_indices_.at(1), residual_indices_.at(1)}; std::vector ctemp = {variable_indices_.at(0), variable_indices_.at(1), variable_indices_.at(0), variable_indices_.at(1)}; std::vector valtemp = {0.0, 0.0, 0.0, 0.0}; J_.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - J_.printMatrix("Bus Jacobian initialization... should be 2x2 with zero values."); return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 603f96992..6bb55707e 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -179,9 +179,6 @@ namespace GridKit MatrixT J_; - RealT time_; - RealT alpha_; - RealT rtol_; RealT atol_; diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp index 33b469879..d14a5b8f2 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp @@ -22,8 +22,8 @@ namespace GridKit template int BusFault::evaluateJacobian() { - std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for BusFault..." << std::endl; + //std::cout << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp index 342e932d1..5ea31c626 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp @@ -14,8 +14,8 @@ namespace GridKit template int BusFault::evaluateJacobian() { - std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for BusFault..." << std::endl; + //std::cout << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp index 369bf74ff..09ea82f1e 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp @@ -22,8 +22,8 @@ namespace GridKit template int BusFault::evaluateJacobian() { - std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + //std::cout << "Evaluate Jacobian for BusFault..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; if (status_) { diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 7fdf3594b..751deb4cb 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -39,15 +39,15 @@ namespace GridKit /// @todo Remove this method. It should be part of DynamicSolver class. virtual bool hasJacobian() override { - return false; + return true; } - // virtual void updateTime(RealT t, RealT a) - // { - // time_ = t; - // alpha_ = a; - // std::cout << "updateTime: t = " << time_ << "\n"; - // } + virtual void updateTime(RealT t, RealT a) override + { + time_ = t; + alpha_ = a; + //std::cout << "updateTime: t = " << time_ << ", alpha = " << alpha_ << "\n"; + } virtual void setTolerances(RealT& rtol, RealT& atol) const override { @@ -180,7 +180,7 @@ namespace GridKit ------ WARNING: Temporary ------ - The protexted variable mva_system_base_ is temporarily + The protected variable mva_system_base_ is temporarily hard coded. This eventually needs to be configured from the input JSON format, which specifies the system MVA base. diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt index b043b66a4..52d50ded9 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/CMakeLists.txt @@ -8,16 +8,32 @@ set(_install_headers Ieeet1.hpp Ieeet1Data.hpp) -gridkit_add_library(phasor_dynamics_exciter_ieeet1 - SOURCES - Ieeet1.cpp - HEADERS - ${_install_headers} - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include - LINK_LIBRARIES - GridKit::phasor_dynamics_core - GridKit::phasor_dynamics_signal) +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library(phasor_dynamics_exciter_ieeet1 + SOURCES + Ieeet1Enzyme.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + PUBLIC GridKit::phasor_dynamics_signal + PRIVATE ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) +else() + gridkit_add_library(phasor_dynamics_exciter_ieeet1 + SOURCES + Ieeet1.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + GridKit::phasor_dynamics_core + GridKit::phasor_dynamics_signal) +endif() gridkit_add_library(phasor_dynamics_exciter_ieeet1_dependency_tracking SOURCES diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp index f626ddea6..4dffe256b 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp @@ -15,6 +15,20 @@ namespace GridKit { namespace Exciter { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeet1::evaluateJacobian() + { + std::cout << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + // Available template instantiations template class Ieeet1; template class Ieeet1; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index fffb5674a..9e8cfa0d8 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -75,6 +75,8 @@ namespace GridKit using Component::time_; using Component::y_; using Component::yp_; + using Component::wb_; + using Component::J_; public: using RealT = typename Component::RealT; @@ -99,9 +101,9 @@ namespace GridKit int tagDifferentiable() override; int evaluateResidual() override; int evaluateJacobian() override; - - void updateTime(RealT /* t */, RealT /* a */) override + bool hasJacobian() override { + return false; } /// Get the `ComponentSignals` from this `Ieeet1` @@ -114,7 +116,11 @@ namespace GridKit return signals_; } +<<<<<<< HEAD const Model::VariableMonitorBase* getMonitor() const override; +======= + __attribute__((always_inline)) inline int evaluateInternalResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*, ScalarT*); +>>>>>>> 0273332f (Functional PhasorDynamics simulations with sparse Jacobians.) private: // Signal pointers @@ -162,6 +168,10 @@ namespace GridKit /// Associate variable getter functions with enum values void initializeMonitor(); + + /* Local copies of signal variables */ + std::vector ws_; + std::map ws_indices_; }; } // namespace Exciter diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp index 7b5c1109a..9068f1cca 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp @@ -15,6 +15,20 @@ namespace GridKit { namespace Exciter { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeet1::evaluateJacobian() + { + std::cout << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + // Available template instantiations template class Ieeet1; template class Ieeet1; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp new file mode 100644 index 000000000..b142927f9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp @@ -0,0 +1,88 @@ +/** + * @file Ieeet1Enzyme.cpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#include + +#include "Ieeet1Impl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeet1::evaluateJacobian() + { + //std::cout << "Evaluate Jacobian for Ieeet1..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; + +// The following will fail to build without a smooth approximation for +// the saturation +#if 0 + J_.zeroMatrix(); + + GridKit::Enzyme::Sparse::InternalJacobianWithSignal, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_); + + GridKit::Enzyme::Sparse::df_dwbWithSignal, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); + + GridKit::Enzyme::Sparse::ExternalJacobian, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + this->getResidualIndices(), + ws_indices_, + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); +#endif + + return 0; + } + + // Available template instantiations + template class Ieeet1; + template class Ieeet1; + + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 6e9351242..6858ef8e2 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -123,6 +123,14 @@ namespace GridKit yp_.resize(size); tag_.resize(size); + // Resize bus data + wb_.resize(2); + + // Resize signal variable data + ws_.resize(1); + ws_[0] = 0.0; + ws_indices_[0] = static_cast(-1); + // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) { @@ -249,66 +257,57 @@ namespace GridKit } /** - * @brief Residuals of system equations + * @brief Internal Residual * */ template - int Ieeet1::evaluateResidual() + __attribute__((always_inline)) inline int Ieeet1::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) { - // Input Variables - ScalarT omega{0}; - - // Set Input Variables - // Meta PR Note: This seems to be very slow, - // but I see how read/write ownership may require this - // - // I believe implementing the equivalent to signal->read() - // at the system level would address this, by routing - // external signals into a generic inputs_ vector - // at the same time as the internal state values y_ - // are recieved from IDA. - if (signals_.template isAttached()) - { - omega = signals_.template readExternalVariable(); - } - // Read E comp (terminal voltage, unless compensation impedance) - ScalarT vreal = bus_->Vr(); - ScalarT vimag = bus_->Vi(); + ScalarT vreal = wb[0]; + ScalarT vimag = wb[1]; Ec_ = std::sqrt(vreal * vreal + vimag * vimag); // Read Internal Variables - ScalarT vts = y_[0]; // y0 - Sensed term volt - ScalarT vr = y_[1]; // y1 - Voltage reg - ScalarT efdp = y_[2]; // y2 - Efd pre mult - ScalarT vfx = y_[3]; // y3 - Exciter feedback - ScalarT vtr = y_[4]; // y4 - Term Volt Err - ScalarT vf = y_[5]; // y5 - Feedback volt - ScalarT ve = y_[6]; // y6 - Excit. Cntrl Volt - ScalarT efd = y_[7]; // y7 - Efd - ScalarT ksat = y_[8]; // y8 - Saturation + ScalarT vts = y[0]; // y0 - Sensed term volt + ScalarT vr = y[1]; // y1 - Voltage reg + ScalarT efdp = y[2]; // y2 - Efd pre mult + ScalarT vfx = y[3]; // y3 - Exciter feedback + ScalarT vtr = y[4]; // y4 - Term Volt Err + ScalarT vf = y[5]; // y5 - Feedback volt + ScalarT ve = y[6]; // y6 - Excit. Cntrl Volt + ScalarT efd = y[7]; // y7 - Efd + ScalarT ksat = y[8]; // y8 - Saturation // Read Internal Derivatives - ScalarT vts_dot = yp_[0]; - ScalarT vr_dot = yp_[1]; - ScalarT efdp_dot = yp_[2]; - ScalarT vfx_dot = yp_[3]; + ScalarT vts_dot = yp[0]; + ScalarT vr_dot = yp[1]; + ScalarT efdp_dot = yp[2]; + ScalarT vfx_dot = yp[3]; + + // Set signal variable aliases + ScalarT omega = ws[0]; // The 'pre-limit' derivative of Pv - ScalarT f = -vr + Ka_ * vtr; - ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, f); + ScalarT func = -vr + Ka_ * vtr; + ScalarT vr_ind = this->indicator(vr, func); // Internal Differential Equations - f_[0] = -vts_dot + (Ec_ - vts) / Tr_; - f_[1] = -vr_dot + vr_ind * f / Ta_; - f_[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_; - f_[3] = -vfx_dot + vf / Tf_; + f[0] = -vts_dot + (Ec_ - vts) / Tr_; + f[1] = -vr_dot + vr_ind * func / Ta_; + f[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_; + f[3] = -vfx_dot + vf / Tf_; // Internal Algebraic Equations - f_[4] = -vtr + vref_ + vUEL_ + vOEL_ + vS_ - vf - vts; - f_[5] = -vf + (efdp * Kf_) / Tf_ - vfx; - f_[6] = -ve + ksat * efdp; - f_[7] = -efd + efdp + omega * efdp * Ispdlim_; + f[4] = -vtr + vref_ + vUEL_ + vOEL_ + vS_ - vf - vts; + f[5] = -vf + (efdp * Kf_) / Tf_ - vfx; + f[6] = -ve + ksat * efdp; + f[7] = -efd + efdp + omega * efdp * Ispdlim_; ScalarT efd_sat = (efdp - SA_) * (Math::sigmoid(efdp - SA_)); f_[8] = -ksat + SB_ * efd_sat * efd_sat; @@ -317,16 +316,34 @@ namespace GridKit } /** - * @brief Jacobian evaluation not implemented yet + * @brief Residual evaluation * - * @tparam ScalarT - Scalar data type - * @tparam IdxT - Index data type - * @return int - error code, 0 = success */ template - int Ieeet1::evaluateJacobian() + int Ieeet1::evaluateResidual() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + // Set Input Variables + // Meta PR Note: This seems to be very slow, + // but I see how read/write ownership may require this + // + // I believe implementing the equivalent to signal->read() + // at the system level would address this, by routing + // external signals into a generic inputs_ vector + // at the same time as the internal state values y_ + // are recieved from IDA. + if (signals_.template isAttached()) + { + ws_[0] = signals_.template readExternalVariable(); + ws_indices_[0] = signals_.template readExternalVariableIndex(); + } + + // Bus voltages + wb_[0] = bus_->Vr(); + wb_[1] = bus_->Vi(); + + // Residual evaluation + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp index 22743387c..6c90f46cd 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp @@ -21,7 +21,7 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + //std::cout << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index 95bb6fda5..ab146abb1 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -91,10 +91,6 @@ namespace GridKit // Still to be implemented int evaluateJacobian() override; - void updateTime(RealT /* t */, RealT /* a */) override - { - } - /// Get the `ComponentSignals` from this `Tgov1` auto getSignals() -> ComponentSignals ws_; std::map ws_indices_; }; diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp index 62dc3af8f..1e46da88a 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp @@ -21,7 +21,7 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + //std::cout << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp index 57ecbb121..57307a81a 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp @@ -24,8 +24,10 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Tgov1..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + //std::cout << "Evaluate Jacobian for Tgov1..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); GridKit::Enzyme::Sparse::InternalJacobianWithSignal, GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, @@ -39,10 +41,9 @@ namespace GridKit yp_.data(), wb_.data(), ws_.data(), + alpha_, J_); - J_.printMatrix("Tgov1 internal Jacobian"); - GridKit::Enzyme::Sparse::ExternalJacobian, GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, ScalarT, @@ -57,8 +58,6 @@ namespace GridKit ws_.data(), J_); - J_.printMatrix("Tgov1 Jacobian after signal evaluation"); - return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index f6add0df1..4e304071f 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -141,7 +141,7 @@ namespace GridKit yp_.resize(size); tag_.resize(size); - // Resize external variable data + // Resize signal variable data ws_.resize(1); ws_[0] = 0.0; ws_indices_[0] = static_cast(-1); @@ -238,7 +238,7 @@ namespace GridKit ScalarT* y, ScalarT* yp, [[maybe_unused]] ScalarT* wb, - [[maybe_unused]] ScalarT* ws, + ScalarT* ws, ScalarT* f) { // Read Internal Variables @@ -250,7 +250,7 @@ namespace GridKit ScalarT ptx_dot = yp[0]; ScalarT pv_dot = yp[1]; - // Set external variable aliases + // Set signal variable aliases ScalarT omega = ws[0]; // The 'pre-limit' derivative of Pv diff --git a/GridKit/Model/PhasorDynamics/Load/Load.cpp b/GridKit/Model/PhasorDynamics/Load/Load.cpp index 8f130cbae..eb2e91ca4 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.cpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Load::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Load..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Load..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/Load.hpp b/GridKit/Model/PhasorDynamics/Load/Load.hpp index 091853b50..0510aabd8 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.hpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.hpp @@ -63,10 +63,6 @@ namespace GridKit return 0; } - virtual void updateTime(RealT /* t */, RealT /* a */) override - { - } - public: void setR(RealT R) { diff --git a/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp index 8219abcad..672998408 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Load::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Load..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Load..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp index b5c1904ac..4f5724b15 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp @@ -22,8 +22,8 @@ namespace GridKit template int Load::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Load..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + //std::cout << "Evaluate Jacobian for Load..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; GridKit::Enzyme::Sparse::BusJacobian, GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp index 7bd6e91bb..1336406a4 100644 --- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp +++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp @@ -67,7 +67,7 @@ namespace GridKit virtual void updateTime(RealT /* t */, RealT /* a */) override { - // No time to update in bus models + // No time to update in signal nodes } virtual void setTolerances(RealT& rtol, RealT& atol) const override @@ -144,9 +144,6 @@ namespace GridKit MatrixT J_; - RealT time_; - RealT alpha_; - RealT rtol_; RealT atol_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp index f9f389486..07825e666 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp @@ -21,8 +21,8 @@ namespace GridKit template int Genrou::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Genrou..." << std::endl; + //std::cout << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp index 8c9aaa61c..8c94a259d 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp @@ -133,10 +133,6 @@ namespace GridKit // Still to be implemented int evaluateJacobian() override; - void updateTime(RealT /* t */, RealT /* a */) override - { - } - // Temporary access functions for governor // Should be abstracted ScalarT getSpeed(); @@ -237,7 +233,7 @@ namespace GridKit ScalarT pmech_set_{0.0}; // TODO remove default initialization and ensure this gets set ScalarT efd_set_{0.0}; // TODO remove default initialization and ensure this gets set - /* Local copies of external variables */ + /* Local copies of signal variables */ std::vector ws_; std::map ws_indices_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp index 0419269d4..4cc9ea806 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp @@ -14,8 +14,8 @@ namespace GridKit template int Genrou::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - std::cout << "Jacobian evaluation not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for Genrou..." << std::endl; + //std::cout << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp index b69846666..d9251be4d 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp @@ -22,8 +22,10 @@ namespace GridKit template int Genrou::evaluateJacobian() { - std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + //std::cout << "Evaluate Jacobian for Genrou..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); GridKit::Enzyme::Sparse::InternalJacobianWithSignal, GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, @@ -37,9 +39,22 @@ namespace GridKit yp_.data(), wb_.data(), ws_.data(), + alpha_, J_); - J_.printMatrix("Genrou internal Jacobian"); + GridKit::Enzyme::Sparse::df_dwbWithSignal, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); GridKit::Enzyme::Sparse::ExternalJacobian, GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, @@ -55,8 +70,6 @@ namespace GridKit ws_.data(), J_); - J_.printMatrix("Genrou Jacobian after signal evaluation"); - GridKit::Enzyme::Sparse::BusJacobian, GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, ScalarT, @@ -70,6 +83,19 @@ namespace GridKit (bus_->y()).data(), bus_->getJacobian()); + GridKit::Enzyme::Sparse::dh_dy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + y_.size(), + bus_->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + J_); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index 7f3ae05fd..4bc02cbed 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -315,11 +315,11 @@ namespace GridKit yp_.resize(static_cast(size_)); tag_.resize(static_cast(size_)); - // Resize coupling data + // Resize bus data wb_.resize(2); h_.resize(2); - // Resize external variable data + // Resize signal variable data ws_.resize(2); ws_indices_[0] = static_cast(-1); ws_indices_[1] = static_cast(-1); @@ -513,7 +513,7 @@ namespace GridKit ScalarT vr = wb[0]; ScalarT vi = wb[1]; - // Set external variable aliases + // Set signal variable aliases ScalarT pmech = ws[0]; ScalarT efd = ws[1]; @@ -551,7 +551,7 @@ namespace GridKit */ template __attribute__((always_inline)) inline int Genrou::evaluateBusResidual( - [[maybe_unused]] ScalarT* y, + ScalarT* y, [[maybe_unused]] ScalarT* yp, ScalarT* wb, ScalarT* h) diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index 52e06f956..bd4c6e8e7 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -15,8 +15,8 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index 7a6d69df3..63139e9ec 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -79,10 +79,6 @@ namespace GridKit // Still to be implemented int evaluateJacobian() override; - void updateTime(RealT /* t */, RealT /* a */) override - { - } - void setPmech(RealT pmech) { pmech_set_ = pmech; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp index d6f4885fa..8e91da2e2 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - std::cout << "Jacobian evaluation is not implemented!" << std::endl; + //std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; + //std::cout << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp index c326316e6..83f80f724 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp @@ -22,8 +22,10 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - std::cout << "Jacobian evaluation is experimental!" << std::endl; + //std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; + //std::cout << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); GridKit::Enzyme::Sparse::InternalJacobian, GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, @@ -36,9 +38,34 @@ namespace GridKit y_.data(), yp_.data(), wb_.data(), + alpha_, J_); - J_.printMatrix("GenClassical internal Jacobian"); + GridKit::Enzyme::Sparse::df_dwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_); + + GridKit::Enzyme::Sparse::dh_dy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + y_.size(), + bus_->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 2d787a02a..7cb331536 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -335,6 +335,7 @@ namespace GridKit this->setResidualIndex(j, j); } + // Verify component configuration int errorCount = this->verify(); if (errorCount > 0) { @@ -365,14 +366,19 @@ namespace GridKit } /** - * @brief Assume that jacobian is not available + * @brief Check components for Jacobian availability * * @return true * @return false */ bool hasJacobian() override { - return false; + bool has_jacobian = true; + for (const auto& component : components_) + { + has_jacobian *= component->hasJacobian(); + } + return has_jacobian; } /** @@ -573,6 +579,7 @@ namespace GridKit */ int evaluateJacobian() override { + J_.zeroMatrix(); std::vector ctemp{}; std::vector rtemp{}; std::vector valtemp{}; @@ -617,6 +624,8 @@ namespace GridKit J_.setValues(rtemp, ctemp, valtemp); + //J_.printMatrix("System Jacobian"); + return 0; } @@ -659,7 +668,8 @@ namespace GridKit */ void updateTime(RealT t, RealT a) override { - this->time_ = t; + time_ = t; + alpha_ = a; for (const auto& component : components_) { component->updateTime(t, a); diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index a05ba1ca8..83ee30d47 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -275,9 +275,12 @@ namespace AnalysisManager if (tag_) initType = IDA_YA_YDP_INIT; - retval = IDACalcIC(solver_, initType, 0.1); + retval = IDACalcIC(solver_, initType, t0 + 0.1); checkOutput(retval, "IDACalcIC"); + retval = IDAGetConsistentIC(solver_, yy_, yp_); + checkOutput(retval, "IDAGetConsistentIC"); + copyVec(yy_, model_->y()); copyVec(yp_, model_->yp()); } diff --git a/buildsystem/spack_repo/gridkit/packages/enzyme/package.py b/buildsystem/spack_repo/gridkit/packages/enzyme/package.py index a57c5cd98..930848d4e 100644 --- a/buildsystem/spack_repo/gridkit/packages/enzyme/package.py +++ b/buildsystem/spack_repo/gridkit/packages/enzyme/package.py @@ -3,6 +3,7 @@ from spack.package import * class Enzyme(BuiltinEnzyme): + version("0.0.216", sha256="c41f1a1286934f5155614fde0f3c63e0f07bc3b67b8f2ff728540bc862cb8f7a") version("0.0.206", sha256="600fd2db370fb40abb6411e0e80df524aea03f2c1ad50a2765ecaab9e1115c77") version("0.0.196", sha256="2b9cfcb7c34e56fc8191423042df06241cf32928eefbb113ac3c5199e3361cb2") version("0.0.186", sha256="125e612df0b6b82b07e1e13218c515bc54e04aa1407e57f4f31d3abe995f4714") diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt b/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt index 0322107bb..c000e4446 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt +++ b/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt @@ -5,4 +5,4 @@ target_link_libraries(TenGenClassical install(TARGETS TenGenClassical RUNTIME DESTINATION bin) # Not used for tesing for now. -# add_test(NAME GenrouTest3_genclassical COMMAND $) +add_test(NAME GenrouTest3_genclassical COMMAND $) diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp index 05c29f6bb..374ca81b2 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp @@ -97,6 +97,9 @@ int main() sys.addComponent(&gen10); sys.addComponent(&fault); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); real_type dt = 1.0 / 4.0 / 60.0; @@ -123,7 +126,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& yval = sys.y(); + std::vector& yval = sys.y(); // Output time out << t << ","; @@ -161,7 +164,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -179,7 +182,7 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); fileout.close(); diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt b/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt index ef1fc7a96..b9baf8795 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt @@ -5,4 +5,4 @@ target_link_libraries(TenGenGenrou install(TARGETS TenGenGenrou RUNTIME DESTINATION bin) # Not used for tesing for now. -# add_test(NAME GenrouTest3 COMMAND $) +add_test(NAME GenrouTest3 COMMAND $) diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp index 4f808c69f..a3ae0d41f 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp @@ -96,6 +96,9 @@ int main() sys.addComponent(&gen10); sys.addComponent(&fault); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); real_type dt = 1.0 / 4.0 / 60.0; @@ -118,7 +121,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& yval = sys.y(); + std::vector& yval = sys.y(); // Output time out << t << ","; @@ -147,7 +150,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -165,7 +168,7 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); fileout.close(); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp index 55a07472b..d9ee49b11 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp @@ -27,11 +27,11 @@ using index_type = size_t; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -43,7 +43,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -73,7 +73,7 @@ int main() { using namespace GridKit::PhasorDynamics; using namespace AnalysisManager::Sundials; - using BusType = BusData::BusType; + using BusType = BusData::BusType; auto error_allowed = static_cast(1e-4); @@ -83,7 +83,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(3); @@ -200,6 +200,9 @@ int main() SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); @@ -210,7 +213,7 @@ int main() auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); // Bus 1 -> +2 // Bus 2 -> +2 @@ -219,10 +222,10 @@ int main() // output.push_back(OutputData{t, - 1.0 + y_val[5], // Gen 1 Speed -> 4 + 1 - 1.0 + y_val[24], // Gen 2 Speed -> 23 + 1 - std::sqrt(y_val[0] * y_val[0] + y_val[1] * y_val[1]), // Bus 1 Vmag - std::sqrt(y_val[2] * y_val[2] + y_val[3] * y_val[3])}); // Bus 2 Vmag + 1.0 + static_cast(y_val[5]), // Gen 1 Speed -> 4 + 1 + 1.0 + static_cast(y_val[24]), // Gen 2 Speed -> 23 + 1 + std::sqrt(static_cast(y_val[0]) * static_cast(y_val[0]) + static_cast(y_val[1]) * static_cast(y_val[1])), // Bus 1 Vmag + std::sqrt(static_cast(y_val[2]) * static_cast(y_val[2]) + static_cast(y_val[3]) * static_cast(y_val[3]))}); // Bus 2 Vmag }; // Set up simulation @@ -230,7 +233,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -248,19 +251,19 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); /* Check worst-case error */ real_type worst_error = 0; real_type worst_error_time = 0; - std::ostream nullout(nullptr); - std::ostream& out = nullout; + //std::ostream nullout(nullptr); + //std::ostream& out = nullout; - // // Uncomment code below to print output to a file: - // std::ofstream fileout; - // fileout.open("Example_ThreeBus_Basic_results.csv"); - // std::ostream& out = fileout; + // Uncomment code below to print output to a file: + std::ofstream fileout; + fileout.open("Example_ThreeBus_Basic_results.csv"); + std::ostream& out = fileout; out << "Time,gen2speed,gen3speed,v2mag,v3mag\n"; out << 0. << "," << 1. << "," << 1. << "," << 1. << "," << 1. << "\n"; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp index e01906ffe..952c7f326 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp @@ -27,11 +27,11 @@ using index_type = size_t; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -43,7 +43,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -124,6 +124,9 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); @@ -134,7 +137,7 @@ int main(int argc, const char* argv[]) auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); // Bus 1 -> +2 // Bus 2 -> +2 @@ -143,10 +146,10 @@ int main(int argc, const char* argv[]) // output.push_back(OutputData{t, - 1.0 + y_val[5], // Gen 1 Speed -> 4 + 1 - 1.0 + y_val[24], // Gen 2 Speed -> 23 + 1 - std::sqrt(y_val[0] * y_val[0] + y_val[1] * y_val[1]), // Bus 1 Vmag - std::sqrt(y_val[2] * y_val[2] + y_val[3] * y_val[3])}); // Bus 2 Vmag + 1.0 + static_cast(y_val[5]), // Gen 1 Speed -> 4 + 1 + 1.0 + static_cast(y_val[24]), // Gen 2 Speed -> 23 + 1 + std::sqrt(static_cast(y_val[0]) * static_cast(y_val[0]) + static_cast(y_val[1]) * static_cast(y_val[1])), // Bus 1 Vmag + std::sqrt(static_cast(y_val[2]) * static_cast(y_val[2]) + static_cast(y_val[3]) * static_cast(y_val[3]))}); // Bus 2 Vmag }; // Set up simulation @@ -154,7 +157,7 @@ int main(int argc, const char* argv[]) ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -172,7 +175,7 @@ int main(int argc, const char* argv[]) ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); sys.stopMonitor(); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp index 5b6d073bf..6825eaeb3 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp @@ -25,15 +25,15 @@ using scalar_type = double; using real_type = double; using index_type = size_t; -using BusType = GridKit::PhasorDynamics::BusData::BusType; +using BusType = GridKit::PhasorDynamics::BusData::BusType; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -45,7 +45,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -84,7 +84,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(3); @@ -103,7 +103,7 @@ int main() // Bus 2 data.bus[2].bus_id = 2; - data.bus[2].bus_type = BusData::BusType::DEFAULT; + data.bus[2].bus_type = BusData::BusType::DEFAULT; data.bus[2].Vr0 = 0.9610827543495831; data.bus[2].Vi0 = -0.13122476630506485; @@ -177,6 +177,9 @@ int main() SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); @@ -187,13 +190,13 @@ int main() auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); output.push_back(OutputData{t, - 1 + y_val[5], - 1 + y_val[10], - std::hypot(y_val[0], y_val[1]), - std::hypot(y_val[2], y_val[3])}); + 1 + static_cast(y_val[5]), + 1 + static_cast(y_val[10]), + std::hypot(static_cast(y_val[0]), static_cast(y_val[1])), + std::hypot(static_cast(y_val[2]), static_cast(y_val[3]))}); }; // Set up simulation @@ -201,7 +204,7 @@ int main() ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -219,19 +222,19 @@ int main() ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); /* Check worst-case error */ real_type worst_error = 0; real_type worst_error_time = 0; - std::ostream nullout(nullptr); - std::ostream& out = nullout; + //std::ostream nullout(nullptr); + //std::ostream& out = nullout; - // // Uncomment code below to print output to a file: - // std::ofstream fileout; - // fileout.open("Example_ThreeBus_Classical_results.csv"); - // std::ostream& out = fileout; + // Uncomment code below to print output to a file: + std::ofstream fileout; + fileout.open("Example_ThreeBus_Classical_results.csv"); + std::ostream& out = fileout; out << "Time,gen2speed,gen3speed,v2mag,v3mag\n"; out << 0. << "," << 1. << "," << 1. << "," << 1. << "," << 1. << "\n"; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp index 1e69f8f95..774586b22 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp @@ -27,11 +27,11 @@ using index_type = size_t; struct OutputData { - real_type t; - scalar_type gen2speed; - scalar_type gen3speed; - scalar_type v2mag; - scalar_type v3mag; + real_type t; + real_type gen2speed; + real_type gen3speed; + real_type v2mag; + real_type v3mag; OutputData& operator-=(const OutputData& other) { @@ -43,7 +43,7 @@ struct OutputData return *this; } - double norm() const + real_type norm() const { return std::max({ std::abs(gen2speed), @@ -124,6 +124,9 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); @@ -134,13 +137,13 @@ int main(int argc, const char* argv[]) auto output_cb = [&](real_type t) { - std::vector& y_val = sys.y(); + std::vector& y_val = sys.y(); output.push_back(OutputData{t, - 1 + y_val[5], - 1 + y_val[10], - std::hypot(y_val[0], y_val[1]), - std::hypot(y_val[2], y_val[3])}); + 1 + static_cast(y_val[5]), + 1 + static_cast(y_val[10]), + std::hypot(static_cast(y_val[0]), static_cast(y_val[1])), + std::hypot(static_cast(y_val[2]), static_cast(y_val[3]))}); }; // Set up simulation @@ -148,7 +151,7 @@ int main(int argc, const char* argv[]) ida.configureSimulation(); // Run simulation, output each `dt` interval - scalar_type start = static_cast(clock()); + real_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); // Run for 1s @@ -166,7 +169,7 @@ int main(int argc, const char* argv[]) ida.initializeSimulation(1.1, false); nout = static_cast(std::round((10.0 - 1.1) / dt)); ida.runSimulation(10.0, nout, output_cb); - double stop = static_cast(clock()); + real_type stop = static_cast(clock()); /* Check worst-case error */ real_type worst_error = 0; diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp index d7f779987..86f9ccaf7 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp @@ -28,7 +28,7 @@ int main() using real_type = double; using index_type = size_t; - using BusType = BusData::BusType; + using BusType = BusData::BusType; std::cout << "Example: TwoBusBasic\n"; @@ -36,7 +36,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(2); @@ -96,6 +96,9 @@ int main() SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); @@ -135,7 +138,10 @@ int main() { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp index 1bdd42fff..021dd0650 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp @@ -76,6 +76,9 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); @@ -115,7 +118,10 @@ int main(int argc, const char* argv[]) { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp index aa8f478ef..5f4d387c1 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp @@ -27,7 +27,7 @@ int main() using real_type = double; using index_type = size_t; - using BusType = BusData::BusType; + using BusType = BusData::BusType; std::cout << "Example: TwoBusTgov1 + IEEET1 Exciter\n"; @@ -35,7 +35,7 @@ int main() // Create model data // - SystemModelData data; + SystemModelData data; // Set bus data data.bus.resize(2); @@ -143,6 +143,9 @@ int main() SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); @@ -189,11 +192,11 @@ int main() // Exc -> 9 States -> Start Idx 24 output.push_back(OutputData{ t, - y_val[0], // Bus Vr - y_val[1], // Bus Vi - y_val[3], // Gen Speed - y_val[23], // Gov Pmech - y_val[26], // Exc Efd + static_cast(y_val[0]), // Bus Vr + static_cast(y_val[1]), // Bus Vi + static_cast(y_val[3]), // Gen Speed + static_cast(y_val[23]), // Gov Pmech + static_cast(y_val[26]), // Exc Efd }); }; diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp index 2e546f0cd..61289c97b 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp @@ -76,6 +76,9 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); @@ -122,11 +125,11 @@ int main(int argc, const char* argv[]) // Exc -> 9 States -> Start Idx 24 output.push_back(OutputData{ t, - y_val[0], // Bus Vr - y_val[1], // Bus Vi - y_val[3], // Gen Speed - y_val[23], // Gov Pmech - y_val[26], // Exc Efd + static_cast(y_val[0]), // Bus Vr + static_cast(y_val[1]), // Bus Vi + static_cast(y_val[3]), // Gen Speed + static_cast(y_val[23]), // Gov Pmech + static_cast(y_val[26]), // Exc Efd }); }; diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp index 4f33d3d3f..538d93076 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp @@ -117,6 +117,9 @@ int main() SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); @@ -156,7 +159,10 @@ int main() { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp index ad3f7ff76..046f5a606 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp @@ -80,6 +80,9 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); + sys.initialize(); + sys.evaluateResidual(); + sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); @@ -119,7 +122,10 @@ int main(int argc, const char* argv[]) { std::vector& y_val = sys.y(); - output.push_back(OutputData{t, y_val[0], y_val[1], y_val[3]}); + output.push_back(OutputData{t, + static_cast(y_val[0]), + static_cast(y_val[1]), + static_cast(y_val[3])}); }; // Set up simulation diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index c6bd98a7c..3151ba9e3 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -236,24 +236,22 @@ namespace GridKit RealT Xdp{0.5}; // Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(H, D, Ra, Xdp); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(H, D, Ra, Xdp); // Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(H, D, Ra, Xdp); /// Compare DependencyTracking dependencies to Enzyme's - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i])); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i])); } return success.report(__func__); } private: - std::vector DependencyTrackingJacobian( + std::vector DependencyTrackingJacobian( const RealT H, const RealT D, const RealT Ra, const RealT Xdp) { DependencyTracking::Variable Vr1{1.0}; ///< Bus-1 real voltage @@ -261,29 +259,90 @@ namespace GridKit PhasorDynamics::Bus bus(Vr1, Vi1); PhasorDynamics::GenClassical gen(&bus, 1, 1.0, 1.0, H, D, Ra, Xdp); + bus.allocate(); - bus.initialize(); gen.allocate(); + + // Get d/dy + bus.initialize(); + gen.initialize(); + + for (size_t i = 0; i < gen.size(); ++i) + { + gen.y()[i].setVariableNumber(i); ///< Generator independent variables + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + gen.size()); // Bus independent variables + } + + bus.evaluateResidual(); + gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking + ///< the dependencies + std::vector residual_y = gen.getResidual(); + + // Get d/dy' + bus.initialize(); gen.initialize(); for (size_t i = 0; i < gen.size(); ++i) { - gen.y()[i].setVariableNumber(i); ///< Independent variables + gen.yp()[i].setVariableNumber(i); ///< Generator independent variables } + bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual = gen.getResidual(); + std::vector residual_yp = gen.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residual.size(); ++i) + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) { - std::cout << i << "th residual: "; - (residual[i]).print(std::cout); + std::cout << i << "th residual, y: "; + (residual_y[i]).print(std::cout); + std::cout << "\n"; + std::cout << i << "th residual, yp: "; + (residual_yp[i]).print(std::cout); std::cout << "\n"; } - return residual; + // Extract the dependencies and add d/dy' to d/dy + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); + DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + auto value_yp = it_yp->second; + dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); + } + else + { + dependencies[i].insert(std::make_pair(index_y, value_y)); + } + } + + // Insert yp dependencies that did not exist in the y dependencies + for (const auto& pair_yp : dependency_yp) + { + auto index_yp = pair_yp.first; + auto value_yp = pair_yp.second; + auto it_y = dependency_y.find(index_yp); + if (it_y == dependency_y.end()) + { + dependencies[i].insert(std::make_pair(index_yp, value_yp)); + } + } + } + + return dependencies; } std::vector EnzymeJacobian( @@ -294,12 +353,24 @@ namespace GridKit PhasorDynamics::Bus bus(Vr1, Vi1); PhasorDynamics::GenClassical gen(&bus, 1, 1.0, 1.0, H, D, Ra, Xdp); + bus.allocate(); - bus.initialize(); gen.allocate(); + + bus.initialize(); gen.initialize(); + gen.updateTime(0.0, 1.0); + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + gen.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + gen.size()); // Reset bus residual indices + } + + bus.evaluateResidual(); gen.evaluateResidual(); + gen.evaluateJacobian(); GridKit::LinearAlgebra::COO_Matrix model_jacobian = gen.getJacobian(); model_jacobian.printMatrix("Model Jacobian"); diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 90a160faa..4a2d53995 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -231,24 +231,22 @@ namespace GridKit auto tol = 10 * std::numeric_limits::epsilon(); // Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(); // Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(); /// Compare DependencyTracking dependencies to Enzyme's - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i], tol)); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol)); } return success.report(__func__); } private: - std::vector DependencyTrackingJacobian() + std::vector DependencyTrackingJacobian() { DependencyTracking::Variable Vr1{1.0}; ///< Bus real voltage DependencyTracking::Variable Vi1{0.0}; ///< Bus imaginary voltage @@ -275,30 +273,89 @@ namespace GridKit 0); bus.allocate(); + gen.allocate(); + + // Get d/dy bus.initialize(); + gen.initialize(); - gen.allocate(); + for (size_t i = 0; i < gen.size(); ++i) + { + gen.y()[i].setVariableNumber(i); ///< Generator independent variables + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + gen.size()); // Bus independent variables + } + + bus.evaluateResidual(); + gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking + ///< the dependencies + std::vector residual_y = gen.getResidual(); + + // Get d/dy' + bus.initialize(); gen.initialize(); for (size_t i = 0; i < gen.size(); ++i) { - gen.y()[i].setVariableNumber(i); ///< Independent variables + gen.yp()[i].setVariableNumber(i); ///< Generator independent variables } bus.evaluateResidual(); gen.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking ///< the dependencies - std::vector residual = gen.getResidual(); + std::vector residual_yp = gen.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residual.size(); ++i) + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) { - std::cout << i << "th residual: "; - (residual[i]).print(std::cout); + std::cout << i << "th residual, y: "; + (residual_y[i]).print(std::cout); + std::cout << "\n"; + std::cout << i << "th residual, yp: "; + (residual_yp[i]).print(std::cout); std::cout << "\n"; } - return residual; + // Extract the dependencies and add d/dy' to d/dy + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); + DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + auto value_yp = it_yp->second; + dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); + } + else + { + dependencies[i].insert(std::make_pair(index_y, value_y)); + } + } + + // Insert yp dependencies that did not exist in the y dependencies + for (const auto& pair_yp : dependency_yp) + { + auto index_yp = pair_yp.first; + auto value_yp = pair_yp.second; + auto it_y = dependency_y.find(index_yp); + if (it_y == dependency_y.end()) + { + dependencies[i].insert(std::make_pair(index_yp, value_yp)); + } + } + } + + return dependencies; + } std::vector EnzymeJacobian() @@ -328,11 +385,20 @@ namespace GridKit 0); bus.allocate(); - bus.initialize(); - bus.evaluateResidual(); - gen.allocate(); + + bus.initialize(); gen.initialize(); + + gen.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + gen.size()); // Reset bus variable indices + bus.setResidualIndex(i, i + gen.size()); // Reset bus residual indices + } + + bus.evaluateResidual(); gen.evaluateResidual(); gen.evaluateJacobian(); diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 99d0e84c8..026128a46 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -272,18 +272,16 @@ namespace GridKit gendata.parameters[GenrouParameters::S12] = 0.; /// Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(busdata, gendata); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(busdata, gendata); /// Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(busdata, gendata); /// Compare DependencyTracking dependencies to Enzyme's auto tol = 10 * std::numeric_limits::epsilon(); - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i], tol)); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol)); } success.expectFailure(); // TODO: Resolve mismatch from sigmoid approximation @@ -292,7 +290,7 @@ namespace GridKit } private: - std::vector DependencyTrackingJacobian( + std::vector DependencyTrackingJacobian( PhasorDynamics::BusData busdata, PhasorDynamics::GenrouData gendata) { @@ -307,6 +305,7 @@ namespace GridKit gov.allocate(); gen.allocate(); + // Get d/dy bus.initialize(); gen.initialize(); gov.initialize(); @@ -321,17 +320,73 @@ namespace GridKit gen.evaluateResidual(); gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking // the dependencies - std::vector residual = gov.getResidual(); + std::vector residual_y = gov.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residual.size(); ++i) + // Get d/dy' + bus.initialize(); + gen.initialize(); + gov.initialize(); + + for (size_t i = 0; i < gov.size(); ++i) + { + gov.yp()[i].setVariableNumber(i); ///< Governor independent variables + } + + bus.evaluateResidual(); + gen.evaluateResidual(); + gov.evaluateResidual(); // Computes the residual and the Jacobian values by tracking + // the dependencies + std::vector residual_yp = gov.getResidual(); + + // Print the dependencies + for (size_t i = 0; i < residual_y.size(); ++i) { - std::cout << i << "th residual: "; - (residual[i]).print(std::cout); + std::cout << i << "th residual, y: "; + (residual_y[i]).print(std::cout); std::cout << "\n"; + std::cout << i << "th residual, yp: "; + (residual_yp[i]).print(std::cout); + std::cout << "\n"; + } + + // 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 residual; + return dependencies; + } std::vector EnzymeJacobian( @@ -355,6 +410,8 @@ namespace GridKit gen.initialize(); gov.initialize(); + gov.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term + bus.evaluateResidual(); gen.evaluateResidual(); gov.evaluateResidual(); diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 727d2656a..1993c94a6 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -204,24 +204,22 @@ namespace GridKit data.branch[0].parameters[BranchParameters::B] = 1.2; // Jacobian via DependencyTracking - std::vector dependency_tracking_residuals = DependencyTrackingJacobian(data); + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(data); // Jacobian via Enzyme std::vector enzyme_jacobian = EnzymeJacobian(data); /// Compare DependencyTracking dependencies to Enzyme's - for (size_t i = 0; i < dependency_tracking_residuals.size(); ++i) + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) { - DependencyTracking::Variable res = dependency_tracking_residuals[i]; - const DependencyTracking::Variable::DependencyMap& dependencies = res.getDependencies(); - success *= (GridKit::Testing::isEqual(dependencies, enzyme_jacobian[i])); + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i])); } return success.report(__func__); } private: - std::vector DependencyTrackingJacobian( + std::vector DependencyTrackingJacobian( PhasorDynamics::SystemModelData data) { // Create an empty system model @@ -239,17 +237,24 @@ namespace GridKit // Evaluate and get the system residuals system.evaluateResidual(); - std::vector residuals = system.getResidual(); + std::vector residual = system.getResidual(); - /// Print the dependencies - for (size_t i = 0; i < residuals.size(); ++i) + // Print the dependencies + for (size_t i = 0; i < residual.size(); ++i) { std::cout << i << "th residual: "; - (residuals[i]).print(std::cout); + (residual[i]).print(std::cout); std::cout << "\n"; } - return residuals; + // Extract the dependencies + std::vector dependencies(residual.size()); + for (IdxT i = 0; i < residual.size(); ++i) + { + dependencies[i] = (residual[i]).getDependencies(); + } + + return dependencies; } std::vector EnzymeJacobian( From e49d7b38ac421646c992d411be95e5d9f18494ce Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sun, 21 Dec 2025 16:29:10 -0500 Subject: [PATCH 02/14] Cleaner Enzyme sparse wrapper implementation. --- CHANGELOG.md | 1 + .../AutomaticDifferentiation/Enzyme/DfDwb.hpp | 171 +++ .../AutomaticDifferentiation/Enzyme/DfDws.hpp | 111 ++ .../AutomaticDifferentiation/Enzyme/DfDy.hpp | 259 ++++ .../AutomaticDifferentiation/Enzyme/DhDwb.hpp | 102 ++ .../AutomaticDifferentiation/Enzyme/DhDy.hpp | 102 ++ .../Enzyme/EnzymeDefinitions.hpp | 38 + .../Enzyme/LowerSparseStorage.hpp | 116 ++ .../Enzyme/ModelWrappers.hpp | 188 +++ .../Enzyme/SparseJacobians.hpp | 13 + .../Enzyme/SparseWrapper.hpp | 1064 ----------------- GridKit/CommonMath.hpp | 2 +- GridKit/Model/Evaluator.hpp | 1 - .../Model/PhasorDynamics/Branch/Branch.cpp | 4 +- .../Branch/BranchDependencyTracking.cpp | 4 +- .../PhasorDynamics/Branch/BranchEnzyme.cpp | 102 +- GridKit/Model/PhasorDynamics/Bus/Bus.cpp | 4 +- .../Bus/BusDependencyTracking.cpp | 4 +- .../Model/PhasorDynamics/Bus/BusEnzyme.cpp | 4 +- GridKit/Model/PhasorDynamics/BusBase.hpp | 3 + .../PhasorDynamics/BusFault/BusFault.cpp | 4 +- .../BusFault/BusFaultDependencyTracking.cpp | 4 +- .../BusFault/BusFaultEnzyme.cpp | 30 +- GridKit/Model/PhasorDynamics/CMakeLists.txt | 1 + GridKit/Model/PhasorDynamics/Component.hpp | 10 +- .../PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp | 3 +- .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 10 +- .../IEEET1/Ieeet1DependencyTracking.cpp | 5 +- .../Exciter/IEEET1/Ieeet1Enzyme.cpp | 98 +- .../Exciter/IEEET1/Ieeet1Impl.hpp | 23 +- .../PhasorDynamics/Governor/Tgov1/Tgov1.cpp | 3 +- .../Tgov1/Tgov1DependencyTracking.cpp | 3 +- .../Governor/Tgov1/Tgov1Enzyme.cpp | 64 +- GridKit/Model/PhasorDynamics/Load/Load.cpp | 4 +- .../Load/LoadDependencyTracking.cpp | 4 +- .../Model/PhasorDynamics/Load/LoadEnzyme.cpp | 30 +- .../SynchronousMachine/GENROUwS/Genrou.cpp | 4 +- .../GENROUwS/GenrouDependencyTracking.cpp | 4 +- .../GENROUwS/GenrouEnzyme.cpp | 124 +- .../GenClassical/GenClassical.cpp | 4 +- .../GenClassicalDependencyTracking.cpp | 4 +- .../GenClassical/GenClassicalEnzyme.cpp | 72 +- GridKit/Model/PhasorDynamics/SystemModel.hpp | 22 +- GridKit/Solver/Dynamic/CMakeLists.txt | 6 +- GridKit/Solver/Dynamic/Ida.cpp | 10 +- GridKit/Solver/Dynamic/Ida.hpp | 3 + .../Tiny/ThreeBus/Basic/ThreeBusBasic.cpp | 4 +- .../ThreeBus/Classical/ThreeBusClassical.cpp | 4 +- .../PhasorDynamics/GenClassicalTests.hpp | 12 +- .../UnitTests/PhasorDynamics/GenrouTests.hpp | 15 +- .../PhasorDynamics/GovernorTgov1Tests.hpp | 13 +- 51 files changed, 1470 insertions(+), 1420 deletions(-) create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp create mode 100644 GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp delete mode 100644 GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp diff --git a/CHANGELOG.md b/CHANGELOG.md index b8dfddd8f..aec17250f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,6 +42,7 @@ - Added capability to print monitored variables in multiple formats, triggered from `Ida::runSimulation`. - Added IDA statistics object which can be accumulated over multiple simulations. - Minor performance improvements to residual evaluation in PowerElectronics module. +- Added full support for sparse Jacobians obtained with Enzyme in PhasorDynamics. ## v0.1 diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp new file mode 100644 index 000000000..c213216ed --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp @@ -0,0 +1,171 @@ +/** + * @file DfDwb.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator : df/dwb + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DfDwb + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_const, + yp, + enzyme_dup, + wb, + output, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in] ws - Signal variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_const, + yp, + enzyme_dup, + wb, + output, + enzyme_const, + ws, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp new file mode 100644 index 000000000..f0811e769 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp @@ -0,0 +1,111 @@ +/** + * @file DfDws.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: df/dws + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DfDws + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in] ws - Signal variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_const, + yp, + enzyme_const, + wb, + enzyme_dup, + ws, + output, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + IdxT row = res_indices.at(static_cast(tup.row)); + IdxT col = var_indices.at(static_cast(tup.col)); + if (col != static_cast(-1)) + { + rtemp.push_back(row); + ctemp.push_back(col); + valtemp.push_back(tup.val); + } + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp new file mode 100644 index 000000000..9b86cc482 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp @@ -0,0 +1,259 @@ +/** + * @file DfDy.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: Internal Jacobian, df/dy + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DfDy + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - 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] alpha - Time derivative jacobian coefficient + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + RealT alpha, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + // df/dy + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_dup, + y, + output, + enzyme_const, + yp, + enzyme_const, + wb, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + + // df/dy' + triplets.clear(); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_dup, + yp, + output, + enzyme_const, + wb, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + ctemp.clear(); + rtemp.clear(); + valtemp.clear(); + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in] ws - Signal variables + * @param[in] alpha - Time derivative jacobian coefficient + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + RealT alpha, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + // df/dy + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_dup, + y, + output, + enzyme_const, + yp, + enzyme_const, + wb, + enzyme_const, + ws, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + + // df/dy' + triplets.clear(); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_dup, + yp, + output, + enzyme_const, + wb, + enzyme_const, + ws, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + ctemp.clear(); + rtemp.clear(); + valtemp.clear(); + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp new file mode 100644 index 000000000..453ed256e --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp @@ -0,0 +1,102 @@ +/** + * @file DhDwb.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: dh/dwb + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DhDwb + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_const, + y, + enzyme_const, + yp, + enzyme_dup, + wb, + output, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.axpy(1.0, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp new file mode 100644 index 000000000..4a78d7073 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp @@ -0,0 +1,102 @@ +/** + * @file DhDy.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme automatic differentiation Jacobian evaluator: dh/dy + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + */ + template + struct DhDy + { + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; + + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] n_res - Number of residual functions + * @param[in] n_var - Number of independent variables + * @param[in] res_indices - Map from local residual indices to global indices + * @param[in] var_indices - Map from local variable indices to global indices + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in,out] jac - Jacobian + */ + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::map& res_indices, + const std::map& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) + { + if (n_res > 0 && n_var > 0) + { + std::vector> triplets; + std::vector elementary_v(n_var); + for (size_t var_i = 0; var_i < n_var; ++var_i) + { + // Sparse storage + ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); + ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); + + // Elementary vector for Jacobian-vector product + std::ranges::fill(elementary_v, 0.0); + elementary_v[var_i] = 1.0; + + // Autodiff + __enzyme_fwddiff((void*) ModelWrapper::eval, + enzyme_const, + model, + enzyme_dup, + y, + output, + enzyme_const, + yp, + enzyme_const, + wb, + enzyme_dupnoneed, + elementary_v.data(), + d_output); + } + + // Store result + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + for (auto& tup : triplets) + { + rtemp.push_back(res_indices.at(static_cast(tup.row))); + ctemp.push_back(var_indices.at(static_cast(tup.col))); + valtemp.push_back(tup.val); + } + jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes + } + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp b/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp new file mode 100644 index 000000000..858385ce7 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp @@ -0,0 +1,38 @@ +/** + * @file EnzymeDefinitions.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include + +/** + * @brief Enzyme constants for activity analysis + * + */ +extern int enzyme_dup; +extern int enzyme_const; +extern int enzyme_dupnoneed; + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme fwddiff template + * + * @tparam T - return type + * @tparam ModelT - model type + */ + template + extern T __enzyme_fwddiff(void*, ModelT...) noexcept; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp new file mode 100644 index 000000000..564584077 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp @@ -0,0 +1,116 @@ +/** + * @file LowerSparseStorage.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Enzyme todense template + * + * @tparam T - return type + */ + template + extern T __enzyme_todense(void*...) noexcept; + + /** + * @brief Enzyme sparse storage in triplet format + * + * @tparam ScalarT - scalar data type + */ + template + struct Triple + { + size_t row; + size_t col; + ScalarT val; + Triple(Triple&&) = default; + + Triple(size_t row, size_t col, ScalarT val) + : row(row), + col(col), + val(val) + { + } + }; + + /** + * @brief Enzyme sparse accumulation for float + * + */ + [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(size_t row, size_t col, float val, std::vector>& triplets) + { + triplets.emplace_back(row, col, val); + } + + /** + * @brief Enzyme sparse accumulation for double + * + */ + [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(size_t row, size_t col, double val, std::vector>& triplets) + { + triplets.emplace_back(row, col, val); + } + + /** + * @brief Enzyme sparse store + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static void sparse_store(ScalarT val, size_t idx, size_t i, std::vector>& triplets) + { + if (val == 0.0) + return; + idx /= sizeof(ScalarT); + if constexpr (sizeof(ScalarT) == 4) + inner_storeflt(idx, i, val, triplets); + else + inner_storedbl(idx, i, val, triplets); + } + + /** + * @brief Enzyme sparse load + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static ScalarT sparse_load(size_t, size_t, std::vector>&) + { + return 0.0; + } + + /** + * @brief Enzyme identity store + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static void ident_store(ScalarT, size_t, size_t) + { + assert(0 && "should never store"); + } + + /** + * @brief Enzyme identity load + * + * @tparam ScalarT - scalar data type + */ + template + __attribute__((always_inline)) static ScalarT ident_load(size_t idx, size_t i) + { + idx /= sizeof(ScalarT); + return (ScalarT) (idx == i); + } + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp new file mode 100644 index 000000000..aff23d92a --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp @@ -0,0 +1,188 @@ +/** + * @file ModelWrappers.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace Enzyme + { + namespace Sparse + { + /** + * @brief Model member function parameter keys + * + */ + enum class MemberFunctions + { + InternalResidual, + InternalResidualWithSignal, + BusResidual, + BusResidual11, //< Special case for branches that are connected to two buses + BusResidual12, //< Special case for branches that are connected to two buses + BusResidual21, //< Special case for branches that are connected to two buses + BusResidual22 //< Special case for branches that are connected to two buses + }; + + /** + * @brief Template definition for wrapper around residual methods inside model classes + * + * @tparam ModelT - model type + * @tparam MemberFunctions - member function parameter key + * @tparam ScalarT - scalar data type + * + */ + template + struct ModelWrapper + { + }; + + /** + * @brief Residual wrapper partial template specialization for InternalResidual + * + */ + template + struct ModelWrapper + { + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[out] f - Internal residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* f) + { + model->evaluateInternalResidual(y, yp, wb, f); + } + }; + + /** + * @brief Residual wrapper partial template specialization for InternalResidualWithSignal + * + */ + template + struct ModelWrapper + { + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[in] ws - Signal variables + * @param[out] f - Internal residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* ws, ScalarT* f) + { + model->evaluateInternalResidual(y, yp, wb, ws, f); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual + * + */ + template + struct ModelWrapper + { + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[out] h - Bus residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual11 + * + */ + template + struct ModelWrapper + { + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[out] h - Bus residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual11(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual12 + * + */ + template + struct ModelWrapper + { + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[out] h - Bus residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual12(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual21 + * + */ + template + struct ModelWrapper + { + /** + * @param[in] model - Pointer to the model to be differentiated + __enzyme_fwddiff((void*) ModelWrapper::eval, + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[out] h - Bus residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual21(y, yp, wb, h); + } + }; + + /** + * @brief Residual wrapper partial template specialization for BusResidual22 + * + */ + template + struct ModelWrapper + { + /** + * @param[in] model - Pointer to the model to be differentiated + * @param[in] y - Internal variables + * @param[in] yp - Internal variable derivatives + * @param[in] wb - Bus variables + * @param[out] h - Bus residual + */ + static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) + { + model->evaluateBusResidual22(y, yp, wb, h); + } + }; + } // namespace Sparse + } // namespace Enzyme +} // namespace GridKit diff --git a/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp b/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp new file mode 100644 index 000000000..9a9df2777 --- /dev/null +++ b/GridKit/AutomaticDifferentiation/Enzyme/SparseJacobians.hpp @@ -0,0 +1,13 @@ +/** + * @file SparseWrapper.hpp + * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) + * + */ + +#pragma once + +#include +#include +#include +#include +#include diff --git a/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp b/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp deleted file mode 100644 index 6beb0a565..000000000 --- a/GridKit/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp +++ /dev/null @@ -1,1064 +0,0 @@ -/** - * @file SparseWrapper.hpp - * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) - * - */ - -#pragma once - -#include -#include -#include -#include - -#include -#include - -/** - * @brief Enzyme constants for activity analysis - * - */ -extern int enzyme_dup; -extern int enzyme_const; -extern int enzyme_dupnoneed; - -namespace GridKit -{ - namespace Enzyme - { - namespace Sparse - { - /** - * @brief Model member function parameter keys - * - */ - enum class MemberFunctions - { - InternalResidual, - InternalResidualWithSignal, - BusResidual, - BusResidual11, //< Special case for branches that are connected to two buses - BusResidual12, //< Special case for branches that are connected to two buses - BusResidual21, //< Special case for branches that are connected to two buses - BusResidual22 //< Special case for branches that are connected to two buses - }; - - /** - * @brief Model independent variable parameter keys - * - */ - enum class IndependentVariables - { - Internal, - Bus, - Signal - }; - - /** - * @brief Template definition for wrapper around residual methods inside model classes - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * - */ - template - struct Wrapper - { - }; - - /** - * @brief Residual wrapper partial template specialization for InternalResidual - * - */ - template - struct Wrapper - { - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] f - Internal residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* f) - { - model->evaluateInternalResidual(y, yp, wb, f); - } - }; - - /** - * @brief Residual wrapper partial template specialization for InternalResidualWithSignal - * - */ - template - struct Wrapper - { - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[out] f - Internal residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* ws, ScalarT* f) - { - model->evaluateInternalResidual(y, yp, wb, ws, f); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual - * - */ - template - struct Wrapper - { - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual11 - * - */ - template - struct Wrapper - { - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual11(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual12 - * - */ - template - struct Wrapper - { - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual12(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual21 - * - */ - template - struct Wrapper - { - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual21(y, yp, wb, h); - } - }; - - /** - * @brief Residual wrapper partial template specialization for BusResidual22 - * - */ - template - struct Wrapper - { - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[out] h - Bus residual - */ - static void eval(ModelT* model, ScalarT* y, ScalarT* yp, ScalarT* wb, ScalarT* h) - { - model->evaluateBusResidual22(y, yp, wb, h); - } - }; - - /** - * @brief Enzyme fwddiff template - * - * @tparam T - return type - * @tparam ModelT - model type - */ - template - extern T __enzyme_fwddiff(void*, ModelT...) noexcept; - - /** - * @brief Enzyme todense template - * - * @tparam T - return type - */ - template - extern T __enzyme_todense(void*...) noexcept; - - /** - * @brief Enzyme sparse storage in triplet format - * - * @tparam ScalarT - scalar data type - */ - template - struct Triple - { - size_t row; - size_t col; - ScalarT val; - Triple(Triple&&) = default; - - Triple(size_t row, size_t col, ScalarT val) - : row(row), - col(col), - val(val) - { - } - }; - - /** - * @brief Enzyme sparse accumulation for float - * - */ - [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(size_t row, size_t col, float val, std::vector>& triplets) - { - triplets.emplace_back(row, col, val); - } - - /** - * @brief Enzyme sparse accumulation for double - * - */ - [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(size_t row, size_t col, double val, std::vector>& triplets) - { - triplets.emplace_back(row, col, val); - } - - /** - * @brief Enzyme sparse store - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static void sparse_store(ScalarT val, size_t idx, size_t i, std::vector>& triplets) - { - if (val == 0.0) - return; - idx /= sizeof(ScalarT); - if constexpr (sizeof(ScalarT) == 4) - inner_storeflt(idx, i, val, triplets); - else - inner_storedbl(idx, i, val, triplets); - } - - /** - * @brief Enzyme sparse load - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static ScalarT sparse_load(size_t, size_t, std::vector>&) - { - return 0.0; - } - - /** - * @brief Enzyme identity store - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static void ident_store(ScalarT, size_t, size_t) - { - assert(0 && "should never store"); - } - - /** - * @brief Enzyme identity load - * - * @tparam ScalarT - scalar data type - */ - template - __attribute__((always_inline)) static ScalarT ident_load(size_t idx, size_t i) - { - idx /= sizeof(ScalarT); - return (ScalarT) (idx == i); - } - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Internal Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct InternalJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] alpha - Time derivative jacobian coefficient - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - RealT alpha, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - // df/dy - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_dup, - y, - output, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - - // df/dy' - triplets.clear(); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_dup, - yp, - output, - enzyme_const, - wb, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - ctemp.clear(); - rtemp.clear(); - valtemp.clear(); - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator : df/dwb - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct df_dwb - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Internal Jacobian with signal variables - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct InternalJacobianWithSignal - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[in] alpha - Time derivative jacobian coefficient - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - RealT alpha, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - // df/dy - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_dup, - y, - output, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_const, - ws, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - - // df/dy' - triplets.clear(); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_dup, - yp, - output, - enzyme_const, - wb, - enzyme_const, - ws, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - ctemp.clear(); - rtemp.clear(); - valtemp.clear(); - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: df/dwb with signal variables - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct df_dwbWithSignal - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_const, - ws, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: External Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct ExternalJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in] ws - Signal variables - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_dup, - ws, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - IdxT row = res_indices.at(static_cast(tup.row)); - IdxT col = var_indices.at(static_cast(tup.col)); - if (col != static_cast(-1)) - { - rtemp.push_back(row); - ctemp.push_back(col); - valtemp.push_back(tup.val); - } - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Bus Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct BusJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.axpy(1.0, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: dh/dy - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct dh_dy - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_dup, - y, - output, - enzyme_const, - yp, - enzyme_const, - wb, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - - /** - * @brief Enzyme automatic differentiation Jacobian evaluator: Branch Jacobian - * - * @tparam ModelT - model type - * @tparam MemberFunctions - member function parameter key - * @tparam ScalarT - scalar data type - * @tparam IdxT - matrix index data type - */ - template - struct BranchJacobian - { - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; - - /** - * @param[in] model - Pointer to the model to be differentiated - * @param[in] n_res - Number of residual functions - * @param[in] n_var - Number of independent variables - * @param[in] res_indices - Map from local residual indices to global indices - * @param[in] var_indices - Map from local variable indices to global indices - * @param[in] y - Internal variables - * @param[in] yp - Internal variable derivatives - * @param[in] wb - Bus variables - * @param[in,out] jac - Jacobian - */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) - { - if (n_res > 0 && n_var > 0) - { - std::vector> triplets; - std::vector elementary_v(n_var); - for (size_t var_i = 0; var_i < n_var; ++var_i) - { - // Sparse storage - ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); - ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); - - // Elementary vector for Jacobian-vector product - std::ranges::fill(elementary_v, 0.0); - elementary_v[var_i] = 1.0; - - // Autodiff - __enzyme_fwddiff((void*) Wrapper::eval, - enzyme_const, - model, - enzyme_const, - y, - enzyme_const, - yp, - enzyme_dup, - wb, - output, - enzyme_dupnoneed, - elementary_v.data(), - d_output); - } - - // Store result - std::vector ctemp{}; - std::vector rtemp{}; - std::vector valtemp{}; - for (auto& tup : triplets) - { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); - valtemp.push_back(tup.val); - } - jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes - } - } - }; - } // namespace Sparse - } // namespace Enzyme -} // namespace GridKit diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index 7866d8f17..f896068e5 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -20,7 +20,7 @@ namespace GridKit * * @tparam ScalarT - scalar data type * - * @param[in] x + * @param[in] x - expected to be of order 1 * @return value of the sigmoid function */ template diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index 9bc9f010e..93c812ffc 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -2,7 +2,6 @@ #include -#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.cpp b/GridKit/Model/PhasorDynamics/Branch/Branch.cpp index 527216e92..b2aff31fe 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.cpp @@ -20,8 +20,8 @@ namespace GridKit template int Branch::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Branch..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp index 4c5780ffb..5f384af29 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp @@ -20,8 +20,8 @@ namespace GridKit template int Branch::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Branch..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp index 8ba9da1c4..485829569 100644 --- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "BranchImpl.hpp" @@ -22,66 +22,66 @@ namespace GridKit template int Branch::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Branch..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Branch..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; J_.zeroMatrix(); // Bus 1 diagonal Jacobian block owned by the bus - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11, - ScalarT, - IdxT>::eval(this, - static_cast(bus1_->size()), - static_cast((bus1_->y()).size()), - bus1_->getResidualIndices(), - bus1_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus1_->y()).data(), - bus1_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11, + ScalarT, + IdxT>::eval(this, + static_cast(bus1_->size()), + static_cast((bus1_->y()).size()), + bus1_->getResidualIndices(), + bus1_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus1_->y()).data(), + bus1_->getJacobian()); // Bus 2 diagonal Jacobian block owned by the bus - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22, - ScalarT, - IdxT>::eval(this, - static_cast(bus2_->size()), - static_cast((bus2_->y()).size()), - bus2_->getResidualIndices(), - bus2_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus2_->y()).data(), - bus2_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22, + ScalarT, + IdxT>::eval(this, + static_cast(bus2_->size()), + static_cast((bus2_->y()).size()), + bus2_->getResidualIndices(), + bus2_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus2_->y()).data(), + bus2_->getJacobian()); // Off-diagonal Jacobian block (Bus2 variables) owned by the branch - GridKit::Enzyme::Sparse::BranchJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12, - ScalarT, - IdxT>::eval(this, - static_cast(bus1_->size()), - static_cast((bus2_->y()).size()), - bus1_->getResidualIndices(), - bus2_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus2_->y()).data(), - J_); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12, + ScalarT, + IdxT>::eval(this, + static_cast(bus1_->size()), + static_cast((bus2_->y()).size()), + bus1_->getResidualIndices(), + bus2_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus2_->y()).data(), + J_); // Off-diagonal Jacobian block (Bus1 variables) owned by the branch - GridKit::Enzyme::Sparse::BranchJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21, - ScalarT, - IdxT>::eval(this, - static_cast(bus2_->size()), - static_cast((bus1_->y()).size()), - bus2_->getResidualIndices(), - bus1_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus1_->y()).data(), - J_); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21, + ScalarT, + IdxT>::eval(this, + static_cast(bus2_->size()), + static_cast((bus1_->y()).size()), + bus2_->getResidualIndices(), + bus1_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus1_->y()).data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.cpp b/GridKit/Model/PhasorDynamics/Bus/Bus.cpp index 261d11038..0971e2504 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Bus::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Bus..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Bus..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp index c010cb816..fe116cf7c 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Bus::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Bus..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Bus..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp index 5618ffb39..d9aa774b0 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp @@ -24,8 +24,8 @@ namespace GridKit template int Bus::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Bus..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Bus..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; J_.zeroMatrix(); diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 6bb55707e..cea918fa3 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -8,11 +8,14 @@ #include #include #include +#include namespace GridKit { namespace PhasorDynamics { + using Log = ::GridKit::Utilities::Logger; + /*! * @brief BusBase model implementation base class. * diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp index d14a5b8f2..9f8124c53 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.cpp @@ -22,8 +22,8 @@ namespace GridKit template int BusFault::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - //std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp index 5ea31c626..6dc1173cb 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp @@ -14,8 +14,8 @@ namespace GridKit template int BusFault::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - //std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp index 09ea82f1e..8defda6df 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "BusFaultImpl.hpp" @@ -22,23 +22,23 @@ namespace GridKit template int BusFault::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for BusFault..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; if (status_) { - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - bus_->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - bus_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + bus_->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + bus_->getJacobian()); } return 0; diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index cab9f03a6..9545732cb 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -18,6 +18,7 @@ gridkit_add_library(phasor_dynamics_core SystemModel.hpp SystemModelData.hpp LINK_LIBRARIES + PUBLIC GridKit::definitions PUBLIC GridKit::utilities_logger INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/nlohmann-json/include diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 751deb4cb..ea78a4a76 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -3,12 +3,16 @@ #include #include +#include #include +#include namespace GridKit { namespace PhasorDynamics { + using Log = ::GridKit::Utilities::Logger; + /** * @brief Component model implementation base class. */ @@ -44,9 +48,9 @@ namespace GridKit virtual void updateTime(RealT t, RealT a) override { - time_ = t; - alpha_ = a; - //std::cout << "updateTime: t = " << time_ << ", alpha = " << alpha_ << "\n"; + time_ = t; + alpha_ = a; + // std::cout << "updateTime: t = " << time_ << ", alpha = " << alpha_ << "\n"; } virtual void setTolerances(RealT& rtol, RealT& atol) const override diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp index 4dffe256b..da21cd779 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.cpp @@ -25,7 +25,8 @@ namespace GridKit template int Ieeet1::evaluateJacobian() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index 9e8cfa0d8..ca2974071 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -101,10 +101,6 @@ namespace GridKit int tagDifferentiable() override; int evaluateResidual() override; int evaluateJacobian() override; - bool hasJacobian() override - { - return false; - } /// Get the `ComponentSignals` from this `Ieeet1` auto getSignals() @@ -116,11 +112,9 @@ namespace GridKit return signals_; } -<<<<<<< HEAD const Model::VariableMonitorBase* getMonitor() const override; -======= + __attribute__((always_inline)) inline int evaluateInternalResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*, ScalarT*); ->>>>>>> 0273332f (Functional PhasorDynamics simulations with sparse Jacobians.) private: // Signal pointers @@ -155,7 +149,7 @@ namespace GridKit ScalarT vUEL_{0}; ScalarT vOEL_{0}; ScalarT vS_{0}; - ScalarT Ec_{0}; // "Compensated" terminal measurment + ScalarT Ec_{0}; // "Compensated" terminal measurment, currently unused /// Component signal extension ComponentSignals signals_; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp index 9068f1cca..9346fabb3 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp @@ -25,10 +25,11 @@ namespace GridKit template int Ieeet1::evaluateJacobian() { - std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } - + // Available template instantiations template class Ieeet1; template class Ieeet1; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp index b142927f9..8ffead080 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "Ieeet1Impl.hpp" @@ -24,58 +24,54 @@ namespace GridKit template int Ieeet1::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Ieeet1..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; - -// The following will fail to build without a smooth approximation for -// the saturation -#if 0 - J_.zeroMatrix(); + Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - GridKit::Enzyme::Sparse::InternalJacobianWithSignal, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - this->getResidualIndices(), - this->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_); + J_.zeroMatrix(); + + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_); + + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + ws_.data(), + J_); + + GridKit::Enzyme::Sparse::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + this->getResidualIndices(), + ws_indices_, + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); - GridKit::Enzyme::Sparse::df_dwbWithSignal, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - this->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); - - GridKit::Enzyme::Sparse::ExternalJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - this->getResidualIndices(), - ws_indices_, - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); -#endif - return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 6858ef8e2..0e3048127 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -128,7 +128,7 @@ namespace GridKit // Resize signal variable data ws_.resize(1); - ws_[0] = 0.0; + ws_[0] = 0.0; ws_indices_[0] = static_cast(-1); // Default variable and residual index mapping to local index @@ -197,7 +197,7 @@ namespace GridKit // Terminal Voltage ScalarT vreal = bus_->Vr(); ScalarT vimag = bus_->Vi(); - Ec_ = std::sqrt(vreal * vreal + vimag * vimag); + ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); // Derived from External initial values ScalarT vr = Ke_ * efd0; @@ -205,10 +205,10 @@ namespace GridKit ScalarT vtr = Ke_ / Ka_ * efd0; // Vref (setpoint = terminal + error) - vref_ = Ec_ + vtr; + vref_ = Ec + vtr; // IVP for Internal Variables - y_[0] = Ec_; // y0 - vts - Sensed term volt + y_[0] = Ec; // y0 - vts - Sensed term volt y_[1] = vr; // y1 - vr - Voltage reg y_[2] = efd0; // y2 - efdp - Efd pre mult y_[3] = vfx; // y3 - vfx - Exciter feedback @@ -271,7 +271,7 @@ namespace GridKit // Read E comp (terminal voltage, unless compensation impedance) ScalarT vreal = wb[0]; ScalarT vimag = wb[1]; - Ec_ = std::sqrt(vreal * vreal + vimag * vimag); + ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); // Read Internal Variables ScalarT vts = y[0]; // y0 - Sensed term volt @@ -289,28 +289,29 @@ namespace GridKit ScalarT vr_dot = yp[1]; ScalarT efdp_dot = yp[2]; ScalarT vfx_dot = yp[3]; - + // Set signal variable aliases ScalarT omega = ws[0]; // The 'pre-limit' derivative of Pv - ScalarT func = -vr + Ka_ * vtr; - ScalarT vr_ind = this->indicator(vr, func); + ScalarT func = -vr + Ka_ * vtr; + ScalarT func_normalized = func * 12.0 / 240.0; // Arbitrary normalization to use standardized sigmoid + ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, func_normalized); // Internal Differential Equations - f[0] = -vts_dot + (Ec_ - vts) / Tr_; + f[0] = -vts_dot + (Ec - vts) / Tr_; f[1] = -vr_dot + vr_ind * func / Ta_; f[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_; f[3] = -vfx_dot + vf / Tf_; // Internal Algebraic Equations - f[4] = -vtr + vref_ + vUEL_ + vOEL_ + vS_ - vf - vts; + f[4] = -vts + vref_ + vUEL_ + vOEL_ + vS_ - vtr - vf; f[5] = -vf + (efdp * Kf_) / Tf_ - vfx; f[6] = -ve + ksat * efdp; f[7] = -efd + efdp + omega * efdp * Ispdlim_; ScalarT efd_sat = (efdp - SA_) * (Math::sigmoid(efdp - SA_)); - f_[8] = -ksat + SB_ * efd_sat * efd_sat; + f[8] = -ksat + SB_ * efd_sat * efd_sat; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp index 6c90f46cd..048f058a8 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.cpp @@ -21,7 +21,8 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - //std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp index 1e46da88a..63476a79b 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1DependencyTracking.cpp @@ -21,7 +21,8 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - //std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp index 57307a81a..aee05e06d 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "Tgov1Impl.hpp" @@ -24,39 +24,39 @@ namespace GridKit template int Tgov1::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Tgov1..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; - - J_.zeroMatrix(); + Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - GridKit::Enzyme::Sparse::InternalJacobianWithSignal, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - this->getResidualIndices(), - this->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_); + J_.zeroMatrix(); - GridKit::Enzyme::Sparse::ExternalJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - this->getResidualIndices(), - ws_indices_, - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_); + + GridKit::Enzyme::Sparse::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + this->getResidualIndices(), + ws_indices_, + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/Load.cpp b/GridKit/Model/PhasorDynamics/Load/Load.cpp index eb2e91ca4..324b0af88 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.cpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Load::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Load..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Load..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp index 672998408..37e9f0238 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int Load::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Load..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Load..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp index 4f5724b15..b5309c6f6 100644 --- a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "LoadImpl.hpp" @@ -22,21 +22,21 @@ namespace GridKit template int Load::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Load..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Load..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - bus_->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - bus_->getJacobian()); + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + bus_->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + bus_->getJacobian()); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp index 07825e666..efcea4fc6 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp @@ -21,8 +21,8 @@ namespace GridKit template int Genrou::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - //std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp index 4cc9ea806..8d4c004ca 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouDependencyTracking.cpp @@ -14,8 +14,8 @@ namespace GridKit template int Genrou::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - //std::cout << "Jacobian evaluation not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp index d9251be4d..8360f1c84 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "GenrouImpl.hpp" @@ -22,79 +22,79 @@ namespace GridKit template int Genrou::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for Genrou..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); + J_.zeroMatrix(); - GridKit::Enzyme::Sparse::InternalJacobianWithSignal, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - this->getResidualIndices(), - this->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - alpha_, - J_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_); - GridKit::Enzyme::Sparse::df_dwbWithSignal, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - this->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); - - GridKit::Enzyme::Sparse::ExternalJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, - ScalarT, - IdxT>::eval(this, - f_.size(), - ws_.size(), - this->getResidualIndices(), - ws_indices_, - y_.data(), - yp_.data(), - wb_.data(), - ws_.data(), - J_); + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + this->getResidualIndices(), + bus_->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); - GridKit::Enzyme::Sparse::BusJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, - ScalarT, - IdxT>::eval(this, - static_cast(bus_->size()), - static_cast(bus_->size()), - bus_->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - bus_->getJacobian()); + GridKit::Enzyme::Sparse::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + this->getResidualIndices(), + ws_indices_, + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_); - GridKit::Enzyme::Sparse::dh_dy, + GridKit::Enzyme::Sparse::DhDwb, GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, ScalarT, IdxT>::eval(this, static_cast(bus_->size()), - y_.size(), + static_cast(bus_->size()), bus_->getResidualIndices(), - this->getVariableIndices(), + bus_->getVariableIndices(), y_.data(), yp_.data(), - wb_.data(), - J_); + (bus_->y()).data(), + bus_->getJacobian()); + + GridKit::Enzyme::Sparse::DhDy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + y_.size(), + bus_->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp index bd4c6e8e7..b13b66247 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.cpp @@ -15,8 +15,8 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp index 8e91da2e2..1fa1a7036 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp @@ -15,8 +15,8 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - //std::cout << "Jacobian evaluation is not implemented!" << std::endl; + Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; return 0; } diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp index 83f80f724..81f8d0a09 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp @@ -4,7 +4,7 @@ * */ -#include +#include #include "GenClassicalImpl.hpp" @@ -22,51 +22,51 @@ namespace GridKit template int GenClassical::evaluateJacobian() { - //std::cout << "Evaluate Jacobian for GenClassical..." << std::endl; - //std::cout << "Jacobian evaluation is experimental!" << std::endl; + Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; - J_.zeroMatrix(); + J_.zeroMatrix(); - GridKit::Enzyme::Sparse::InternalJacobian, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - y_.size(), - this->getResidualIndices(), - this->getVariableIndices(), - y_.data(), - yp_.data(), - wb_.data(), - alpha_, - J_); + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + this->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + alpha_, + J_); - GridKit::Enzyme::Sparse::df_dwb, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, - ScalarT, - IdxT>::eval(this, - f_.size(), - static_cast(bus_->size()), - this->getResidualIndices(), - bus_->getVariableIndices(), - y_.data(), - yp_.data(), - (bus_->y()).data(), - J_); - - GridKit::Enzyme::Sparse::dh_dy, - GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidual, ScalarT, IdxT>::eval(this, + f_.size(), static_cast(bus_->size()), - y_.size(), - bus_->getResidualIndices(), - this->getVariableIndices(), + this->getResidualIndices(), + bus_->getVariableIndices(), y_.data(), yp_.data(), - wb_.data(), + (bus_->y()).data(), J_); + GridKit::Enzyme::Sparse::DhDy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + y_.size(), + bus_->getResidualIndices(), + this->getVariableIndices(), + y_.data(), + yp_.data(), + wb_.data(), + J_); + return 0; } diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 7cb331536..17ad5e2c0 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -4,13 +4,13 @@ #include #include +#include #include #include #include #include #include #include -#include // Include all components #include @@ -19,8 +19,6 @@ namespace GridKit { namespace PhasorDynamics { - using Log = ::GridKit::Utilities::Logger; - /** * @brief Prototype for a system model class * @@ -373,11 +371,23 @@ namespace GridKit */ bool hasJacobian() override { - bool has_jacobian = true; + bool has_jacobian = false; +#ifdef GRIDKIT_ENABLE_ENZYME + has_jacobian = true; for (const auto& component : components_) { - has_jacobian *= component->hasJacobian(); + has_jacobian = has_jacobian && component->hasJacobian(); + } + + if (!has_jacobian) + { + Log::warning() << "GritKit was built with Enzyme, but some models don't have Jacobians available. " + << "Falling back to dense Jacobians in PhasorDynamics.\n"; } +#else + Log::warning() << "GritKit was not built with Enzyme. " + << "Falling back to dense Jacobians in PhasorDynamics.\n"; +#endif return has_jacobian; } @@ -624,7 +634,7 @@ namespace GridKit J_.setValues(rtemp, ctemp, valtemp); - //J_.printMatrix("System Jacobian"); + // J_.printMatrix("System Jacobian"); return 0; } diff --git a/GridKit/Solver/Dynamic/CMakeLists.txt b/GridKit/Solver/Dynamic/CMakeLists.txt index c8e77529a..0cceed1fd 100644 --- a/GridKit/Solver/Dynamic/CMakeLists.txt +++ b/GridKit/Solver/Dynamic/CMakeLists.txt @@ -19,7 +19,8 @@ if (GRIDKIT_ENABLE_SUNDIALS_SPARSE) PUBLIC SUNDIALS::nvecserial PUBLIC SUNDIALS::idas PUBLIC SUNDIALS::sunlinsolklu - PUBLIC GridKit::definitions) + PUBLIC GridKit::definitions + PUBLIC GridKit::utilities_logger) else() gridkit_add_library(solvers_dyn SOURCES @@ -29,5 +30,6 @@ else() LINK_LIBRARIES PUBLIC SUNDIALS::nvecserial PUBLIC SUNDIALS::idas - PUBLIC GridKit::definitions) + PUBLIC GridKit::definitions + PUBLIC GridKit::utilities_logger) endif() diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index 83ee30d47..aec2d9b0c 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -142,15 +142,13 @@ namespace AnalysisManager this->configureLinearSolverDense(); } #else + /// Todo - Improve error handling capabilities and hasJacobian_ ownership if (model_->hasJacobian()) { - /// Todo - Improve error handling capabilities and hasJacobian_ ownership - throw std::runtime_error("SUNDIALS is not configured with KLU, but the model has a (sparse) Jacobian."); - } - else - { - this->configureLinearSolverDense(); + Log::warning() << "SUNDIALS is not configured with KLU, but the model has a (sparse) Jacobian. " + << "Falling back to dense Jacobian.\n"; } + this->configureLinearSolverDense(); #endif return retval; diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 650e28852..20ff56aa8 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -19,11 +19,14 @@ #include #include +#include namespace AnalysisManager { namespace Sundials { + using Log = ::GridKit::Utilities::Logger; + struct IdaStats { long int num_steps_ = 0; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp index d9ee49b11..467dc7543 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp @@ -257,8 +257,8 @@ int main() real_type worst_error = 0; real_type worst_error_time = 0; - //std::ostream nullout(nullptr); - //std::ostream& out = nullout; + // std::ostream nullout(nullptr); + // std::ostream& out = nullout; // Uncomment code below to print output to a file: std::ofstream fileout; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp index 6825eaeb3..6cfb90456 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp @@ -228,8 +228,8 @@ int main() real_type worst_error = 0; real_type worst_error_time = 0; - //std::ostream nullout(nullptr); - //std::ostream& out = nullout; + // std::ostream nullout(nullptr); + // std::ostream& out = nullout; // Uncomment code below to print output to a file: std::ofstream fileout; diff --git a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp index 3151ba9e3..a086c7545 100644 --- a/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenClassicalTests.hpp @@ -310,31 +310,31 @@ namespace GridKit 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_y = (residual_y[i]).getDependencies(); DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - for (const auto& pair_y : dependency_y) + 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); + 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 + 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) + 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); + auto it_y = dependency_y.find(index_yp); if (it_y == dependency_y.end()) { dependencies[i].insert(std::make_pair(index_yp, value_yp)); diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 4a2d53995..081bd5f0e 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -322,31 +322,31 @@ namespace GridKit 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_y = (residual_y[i]).getDependencies(); DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - for (const auto& pair_y : dependency_y) + 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); + 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 + 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) + 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); + auto it_y = dependency_y.find(index_yp); if (it_y == dependency_y.end()) { dependencies[i].insert(std::make_pair(index_yp, value_yp)); @@ -355,7 +355,6 @@ namespace GridKit } return dependencies; - } std::vector EnzymeJacobian() @@ -391,7 +390,7 @@ namespace GridKit gen.initialize(); gen.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term - + for (size_t i = 0; i < bus.size(); ++i) { bus.setVariableIndex(i, i + gen.size()); // Reset bus variable indices diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 026128a46..9054bbda9 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -353,31 +353,31 @@ namespace GridKit 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_y = (residual_y[i]).getDependencies(); DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); - for (const auto& pair_y : dependency_y) + 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); + 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 + 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) + 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); + auto it_y = dependency_y.find(index_yp); if (it_y == dependency_y.end()) { dependencies[i].insert(std::make_pair(index_yp, value_yp)); @@ -386,7 +386,6 @@ namespace GridKit } return dependencies; - } std::vector EnzymeJacobian( From d70afcea66925a905491d1f7243f391f63e8c7d8 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 22 Dec 2025 12:37:59 -0500 Subject: [PATCH 03/14] Replace int-int maps with vector of integers for residual/variable indices. --- .../AutomaticDifferentiation/Enzyme/DfDwb.hpp | 54 ++++++++-------- .../AutomaticDifferentiation/Enzyme/DfDws.hpp | 28 ++++----- .../AutomaticDifferentiation/Enzyme/DfDy.hpp | 62 +++++++++---------- .../AutomaticDifferentiation/Enzyme/DhDwb.hpp | 26 ++++---- .../AutomaticDifferentiation/Enzyme/DhDy.hpp | 26 ++++---- .../Model/PhasorDynamics/Branch/Branch.hpp | 2 + GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp | 2 + .../Model/PhasorDynamics/Bus/BusInfinite.hpp | 2 + GridKit/Model/PhasorDynamics/BusBase.hpp | 18 +++--- GridKit/Model/PhasorDynamics/Component.hpp | 18 +++--- .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 4 +- .../Exciter/IEEET1/Ieeet1Impl.hpp | 3 + .../PhasorDynamics/Governor/Tgov1/Tgov1.hpp | 4 +- .../Governor/Tgov1/Tgov1Impl.hpp | 3 + GridKit/Model/PhasorDynamics/Load/Load.hpp | 2 + .../SynchronousMachine/GENROUwS/Genrou.hpp | 4 +- .../GENROUwS/GenrouImpl.hpp | 12 ++-- .../GenClassical/GenClassical.hpp | 2 + .../GenClassical/GenClassicalImpl.hpp | 2 + GridKit/Model/PhasorDynamics/SystemModel.hpp | 4 ++ 20 files changed, 153 insertions(+), 125 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp index c213216ed..f2a6abacf 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp @@ -36,22 +36,22 @@ namespace GridKit * @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] res_indices - Global residual indices + * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in,out] jac - Jacobian */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) { if (n_res > 0 && n_var > 0) { @@ -89,8 +89,8 @@ namespace GridKit std::vector valtemp{}; for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes @@ -101,24 +101,24 @@ namespace GridKit * @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] res_indices - Global residual indices + * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in] ws - Signal variables * @param[in,out] jac - Jacobian */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - MatrixT& jac) + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + MatrixT& jac) { if (n_res > 0 && n_var > 0) { @@ -158,8 +158,8 @@ namespace GridKit std::vector valtemp{}; for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp index f0811e769..d68a12780 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp @@ -36,24 +36,24 @@ namespace GridKit * @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] res_indices - Global residual indices + * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in] ws - Signal variables * @param[in,out] jac - Jacobian */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - MatrixT& jac) + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + MatrixT& jac) { if (n_res > 0 && n_var > 0) { @@ -93,8 +93,8 @@ namespace GridKit std::vector valtemp{}; for (auto& tup : triplets) { - IdxT row = res_indices.at(static_cast(tup.row)); - IdxT col = var_indices.at(static_cast(tup.col)); + IdxT row = res_indices[static_cast(tup.row)]; + IdxT col = var_indices[static_cast(tup.col)]; if (col != static_cast(-1)) { rtemp.push_back(row); diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp index 9b86cc482..ee2575ef9 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp @@ -36,24 +36,24 @@ namespace GridKit * @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] res_indices - Global residual indices + * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in] alpha - Time derivative jacobian coefficient * @param[in,out] jac - Jacobian */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - RealT alpha, - MatrixT& jac) + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + RealT alpha, + MatrixT& jac) { if (n_res > 0 && n_var > 0) { @@ -92,8 +92,8 @@ namespace GridKit std::vector valtemp{}; for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes @@ -132,8 +132,8 @@ namespace GridKit valtemp.clear(); for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes @@ -153,17 +153,17 @@ namespace GridKit * @param[in] alpha - Time derivative jacobian coefficient * @param[in,out] jac - Jacobian */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - ScalarT* ws, - RealT alpha, - MatrixT& jac) + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + RealT alpha, + MatrixT& jac) { if (n_res > 0 && n_var > 0) { @@ -204,8 +204,8 @@ namespace GridKit std::vector valtemp{}; for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes @@ -246,8 +246,8 @@ namespace GridKit valtemp.clear(); for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.axpy(alpha, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp index 453ed256e..87408de16 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp @@ -36,22 +36,22 @@ namespace GridKit * @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] res_indices - Global residual indices + * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in,out] jac - Jacobian */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) { if (n_res > 0 && n_var > 0) { @@ -89,8 +89,8 @@ namespace GridKit std::vector valtemp{}; for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.axpy(1.0, rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp index 4a78d7073..9e05f30d2 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp @@ -36,22 +36,22 @@ namespace GridKit * @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] res_indices - Global residual indices + * @param[in] var_indices - Global variable indices * @param[in] y - Internal variables * @param[in] yp - Internal variable derivatives * @param[in] wb - Bus variables * @param[in,out] jac - Jacobian */ - static void eval(ModelT* model, - size_t n_res, - size_t n_var, - const std::map& res_indices, - const std::map& var_indices, - ScalarT* y, - ScalarT* yp, - ScalarT* wb, - MatrixT& jac) + static void eval(ModelT* model, + size_t n_res, + size_t n_var, + const std::vector& res_indices, + const std::vector& var_indices, + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + MatrixT& jac) { if (n_res > 0 && n_var > 0) { @@ -89,8 +89,8 @@ namespace GridKit std::vector valtemp{}; for (auto& tup : triplets) { - rtemp.push_back(res_indices.at(static_cast(tup.row))); - ctemp.push_back(var_indices.at(static_cast(tup.col))); + rtemp.push_back(res_indices[static_cast(tup.row)]); + ctemp.push_back(var_indices[static_cast(tup.col)]); valtemp.push_back(tup.val); } jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes diff --git a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp index a45391ebb..123eeb821 100644 --- a/GridKit/Model/PhasorDynamics/Branch/Branch.hpp +++ b/GridKit/Model/PhasorDynamics/Branch/Branch.hpp @@ -52,6 +52,8 @@ namespace GridKit using Component::wb_; using Component::h_; using Component::J_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 37b70c01d..77fe6ecf0 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -86,6 +86,8 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) diff --git a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp index ad97593f1..a186f02c5 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp @@ -22,6 +22,8 @@ namespace GridKit using BusBase::yp_; using BusBase::f_; using BusBase::J_; + using BusBase::variable_indices_; + using BusBase::residual_indices_; public: using RealT = typename BusBase::RealT; diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index cea918fa3..04ea16eeb 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -136,10 +136,10 @@ namespace GridKit IdxT getVariableIndex(IdxT local_index) const { - return variable_indices_.at(local_index); + return variable_indices_[local_index]; } - const std::map& getVariableIndices() const + const std::vector& getVariableIndices() const { return variable_indices_; } @@ -152,10 +152,10 @@ namespace GridKit IdxT getResidualIndex(IdxT local_index) const { - return residual_indices_.at(local_index); + return residual_indices_[local_index]; } - const std::map& getResidualIndices() const + const std::vector& getResidualIndices() const { return residual_indices_; } @@ -165,12 +165,10 @@ namespace GridKit protected: IdxT bus_id_{static_cast(-1)}; - IdxT size_{0}; - IdxT nnz_{0}; - std::map variable_indices_; ///< Map between local and global (system-level) - /// variable indices - std::map residual_indices_; ///< Map between local and global (system-level) - /// residual indices + IdxT size_{0}; + IdxT nnz_{0}; + std::vector variable_indices_; ///< Global (system-level) variable indices + std::vector residual_indices_; ///< Global (system-level) residual indices /// Variable monitor std::unique_ptr monitor_; diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index ea78a4a76..3c355e2e5 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -129,10 +129,10 @@ namespace GridKit IdxT& getVariableIndex(IdxT local_index) { - return variable_indices_.at(local_index); + return variable_indices_[local_index]; } - const std::map& getVariableIndices() const + const std::vector& getVariableIndices() const { return variable_indices_; } @@ -148,19 +148,17 @@ namespace GridKit return residual_indices_.at(local_index); } - const std::map& getResidualIndices() const + const std::vector& getResidualIndices() const { return residual_indices_; } protected: - IdxT gridkit_component_id_{0}; - IdxT size_{0}; - IdxT nnz_{0}; - std::map variable_indices_; ///< Map between local and global (system-level) - /// variable indices - std::map residual_indices_; ///< Map between local and global (system-level) - /// residual indices + IdxT gridkit_component_id_{0}; + IdxT size_{0}; + IdxT nnz_{0}; + std::vector variable_indices_; ///< Global (system-level) variable indices + std::vector residual_indices_; ///< Global (system-level) residual indices std::vector y_; std::vector yp_; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index ca2974071..95066f4d6 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -77,6 +77,8 @@ namespace GridKit using Component::yp_; using Component::wb_; using Component::J_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; @@ -165,7 +167,7 @@ namespace GridKit /* Local copies of signal variables */ std::vector ws_; - std::map ws_indices_; + std::vector ws_indices_; }; } // namespace Exciter diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 0e3048127..1b8634b8c 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -122,12 +122,15 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); // Resize bus data wb_.resize(2); // Resize signal variable data ws_.resize(1); + ws_indices_.resize(1); ws_[0] = 0.0; ws_indices_[0] = static_cast(-1); diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index ab146abb1..320e8c290 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -70,6 +70,8 @@ namespace GridKit using Component::wb_; using Component::h_; using Component::J_; + using Component::variable_indices_; + using Component::residual_indices_; using RealT = typename Component::RealT; using model_data_type = Tgov1Data; @@ -125,7 +127,7 @@ namespace GridKit /* Local copies of signal variables */ std::vector ws_; - std::map ws_indices_; + std::vector ws_indices_; }; } // namespace Governor diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 4e304071f..0cc282e21 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -140,9 +140,12 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); // Resize signal variable data ws_.resize(1); + ws_indices_.resize(1); ws_[0] = 0.0; ws_indices_[0] = static_cast(-1); diff --git a/GridKit/Model/PhasorDynamics/Load/Load.hpp b/GridKit/Model/PhasorDynamics/Load/Load.hpp index 0510aabd8..8ffe29482 100644 --- a/GridKit/Model/PhasorDynamics/Load/Load.hpp +++ b/GridKit/Model/PhasorDynamics/Load/Load.hpp @@ -39,6 +39,8 @@ namespace GridKit using Component::tag_; using Component::wb_; using Component::h_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp index 8c94a259d..1ea500dda 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp @@ -82,6 +82,8 @@ namespace GridKit using Component::h_; using Component::J_; using Component::mva_system_base_; + using Component::variable_indices_; + using Component::residual_indices_; public: using RealT = typename Component::RealT; @@ -235,7 +237,7 @@ namespace GridKit /* Local copies of signal variables */ std::vector ws_; - std::map ws_indices_; + std::vector ws_indices_; /// Variable monitor std::unique_ptr monitor_; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index 4bc02cbed..7d20070f2 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -310,10 +310,13 @@ namespace GridKit int Genrou::allocate() { // Resize component model data - f_.resize(static_cast(size_)); - y_.resize(static_cast(size_)); - yp_.resize(static_cast(size_)); - tag_.resize(static_cast(size_)); + auto size = static_cast(size_); + f_.resize(size); + y_.resize(size); + yp_.resize(size); + tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); // Resize bus data wb_.resize(2); @@ -321,6 +324,7 @@ namespace GridKit // Resize signal variable data ws_.resize(2); + ws_indices_.resize(2); ws_indices_[0] = static_cast(-1); ws_indices_[1] = static_cast(-1); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp index 63139e9ec..d1ca4c29e 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp @@ -46,6 +46,8 @@ namespace GridKit using Component::h_; using Component::J_; using Component::mva_system_base_; + using Component::variable_indices_; + using Component::residual_indices_; public: using bus_type = BusBase; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp index 58ba7039a..bec605fcc 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp @@ -170,6 +170,8 @@ namespace GridKit y_.resize(size); yp_.resize(size); tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); // Resize coupling data wb_.resize(2); diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 17ad5e2c0..92cbac041 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -49,6 +49,8 @@ namespace GridKit using PhasorDynamics::Component::J_; using PhasorDynamics::Component::rel_tol_; using PhasorDynamics::Component::abs_tol_; + using PhasorDynamics::Component::variable_indices_; + using PhasorDynamics::Component::residual_indices_; public: /** @@ -325,6 +327,8 @@ namespace GridKit yp_.resize(size_); f_.resize(size_); tag_.resize(size_); + variable_indices_.resize(size_); + residual_indices_.resize(size_); // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) From d0763c90639c6a1e7ced8091f19871d8a7d7b868 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 22 Dec 2025 12:58:30 -0500 Subject: [PATCH 04/14] Make use of INVALID_INDEX constant. --- .../AutomaticDifferentiation/DependencyTracking/Variable.hpp | 3 ++- GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp | 3 ++- GridKit/Model/PhasorDynamics/BusBase.hpp | 3 ++- GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp | 2 +- GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp | 2 +- GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp | 2 +- .../PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp | 4 ++-- .../Model/PowerElectronics/SystemModelPowerElectronics.hpp | 3 ++- 8 files changed, 13 insertions(+), 9 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp index b7adffdc9..b6285abc0 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp @@ -14,6 +14,7 @@ #include #include +#include #include namespace GridKit @@ -266,7 +267,7 @@ namespace GridKit bool is_fixed_; ///< Constant parameter flag. mutable DependencyMap* dependencies_; - static const size_t INVALID_VAR_NUMBER = static_cast(-1); + static const size_t INVALID_VAR_NUMBER = INVALID_INDEX; }; //------------------------------------ diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp index d68a12780..7e8aa78c9 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp @@ -6,6 +6,7 @@ #pragma once +#include #include #include #include @@ -95,7 +96,7 @@ namespace GridKit { IdxT row = res_indices[static_cast(tup.row)]; IdxT col = var_indices[static_cast(tup.col)]; - if (col != static_cast(-1)) + if (col != INVALID_INDEX) { rtemp.push_back(row); ctemp.push_back(col); diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 04ea16eeb..e6120114e 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -4,6 +4,7 @@ #include #include +#include #include #include #include @@ -163,7 +164,7 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; protected: - IdxT bus_id_{static_cast(-1)}; + IdxT bus_id_{INVALID_INDEX}; IdxT size_{0}; IdxT nnz_{0}; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 1b8634b8c..f2cffa300 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -132,7 +132,7 @@ namespace GridKit ws_.resize(1); ws_indices_.resize(1); ws_[0] = 0.0; - ws_indices_[0] = static_cast(-1); + ws_indices_[0] = INVALID_INDEX; // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 0cc282e21..2d7c7d87a 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -147,7 +147,7 @@ namespace GridKit ws_.resize(1); ws_indices_.resize(1); ws_[0] = 0.0; - ws_indices_[0] = static_cast(-1); + ws_indices_[0] = INVALID_INDEX; // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp index 1336406a4..861b71dc6 100644 --- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp +++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp @@ -131,7 +131,7 @@ namespace GridKit IdxT signal_id_{0}; protected: - const IdxT bus_id_{static_cast(-1)}; + const IdxT bus_id_{INVALID_INDEX}; IdxT size_{0}; IdxT nnz_{0}; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index 7d20070f2..56fae73b6 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -325,8 +325,8 @@ namespace GridKit // Resize signal variable data ws_.resize(2); ws_indices_.resize(2); - ws_indices_[0] = static_cast(-1); - ws_indices_[1] = static_cast(-1); + ws_indices_[0] = INVALID_INDEX; + ws_indices_[1] = INVALID_INDEX; // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 1bffd0d8d..f5796fc88 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -409,7 +410,7 @@ namespace GridKit } private: - static constexpr IdxT neg1_ = static_cast(-1); + static constexpr IdxT neg1_ = INVALID_INDEX; std::vector components_; From 0429a1048673dfab4a2f022502f2fef21405a7fe Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 22 Dec 2025 13:06:55 -0500 Subject: [PATCH 05/14] Suppress warnings. --- GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp | 2 +- GridKit/Model/PhasorDynamics/BusBase.hpp | 8 ++++---- GridKit/Model/PhasorDynamics/Component.hpp | 8 ++++---- GridKit/Model/VariableMonitorController.hpp | 10 +++++----- GridKit/Model/VariableMonitorImpl.hpp | 8 ++++---- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp index 7e8aa78c9..0155193df 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp @@ -6,10 +6,10 @@ #pragma once -#include #include #include #include +#include #include #include diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index e6120114e..9bb31d747 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -131,13 +131,13 @@ namespace GridKit int setVariableIndex(IdxT local_index, IdxT global_index) { - variable_indices_[local_index] = global_index; + variable_indices_[static_cast(local_index)] = global_index; return 0; } IdxT getVariableIndex(IdxT local_index) const { - return variable_indices_[local_index]; + return variable_indices_[static_cast(local_index)]; } const std::vector& getVariableIndices() const @@ -147,13 +147,13 @@ namespace GridKit int setResidualIndex(IdxT local_index, IdxT global_index) { - residual_indices_[local_index] = global_index; + residual_indices_[static_cast(local_index)] = global_index; return 0; } IdxT getResidualIndex(IdxT local_index) const { - return residual_indices_[local_index]; + return residual_indices_[static_cast(local_index)]; } const std::vector& getResidualIndices() const diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 3c355e2e5..882ad5759 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -123,13 +123,13 @@ namespace GridKit int setVariableIndex(IdxT local_index, IdxT global_index) { - variable_indices_[local_index] = global_index; + variable_indices_[static_cast(local_index)] = global_index; return 0; } IdxT& getVariableIndex(IdxT local_index) { - return variable_indices_[local_index]; + return variable_indices_[static_cast(local_index)]; } const std::vector& getVariableIndices() const @@ -139,13 +139,13 @@ namespace GridKit int setResidualIndex(IdxT local_index, IdxT global_index) { - residual_indices_[local_index] = global_index; + residual_indices_[static_cast(local_index)] = global_index; return 0; } IdxT& getResidualIndex(IdxT local_index) { - return residual_indices_.at(local_index); + return residual_indices_[static_cast(local_index)]; } const std::vector& getResidualIndices() const diff --git a/GridKit/Model/VariableMonitorController.hpp b/GridKit/Model/VariableMonitorController.hpp index 08bc9a779..79db154ca 100644 --- a/GridKit/Model/VariableMonitorController.hpp +++ b/GridKit/Model/VariableMonitorController.hpp @@ -164,7 +164,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { printFullHeader(*sink.os, sink.format); }, + { this->printFullHeader(*sink.os, sink.format); }, sink); } } @@ -256,7 +256,7 @@ namespace GridKit { std::visit([this](auto&& sink) { - printFull(*sink.os, sink.format); + this->printFull(*sink.os, sink.format); using T = std::remove_cvref_t; if constexpr (std::is_same_v>) { @@ -295,7 +295,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { printFullFooter(*sink.os, sink.format); }, + { this->printFullFooter(*sink.os, sink.format); }, sink); } } @@ -403,7 +403,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { sink.start(); }, + { (void) this; sink.start(); }, sink); } } @@ -416,7 +416,7 @@ namespace GridKit for (auto&& sink : sinks_) { std::visit([this](auto&& sink) - { sink.stop(); }, + { (void) this; sink.stop(); }, sink); } } diff --git a/GridKit/Model/VariableMonitorImpl.hpp b/GridKit/Model/VariableMonitorImpl.hpp index fdcca543b..eda73332e 100644 --- a/GridKit/Model/VariableMonitorImpl.hpp +++ b/GridKit/Model/VariableMonitorImpl.hpp @@ -97,7 +97,7 @@ namespace GridKit template void set(VariableEnum v, FuncT f) { - f_[enum_integer(v)] = ValuePrinter{f}; + f_[static_cast(enum_integer(v))] = ValuePrinter{f}; } private: @@ -166,7 +166,7 @@ namespace GridKit { for (auto v : variables_) { - os << csv.delim << f_[enum_integer(v)]; + os << csv.delim << f_[static_cast(enum_integer(v))]; } } @@ -176,7 +176,7 @@ namespace GridKit void print(std::ostream& os, VariableEnum v, Json) const { os << indent_ << std::quoted(enum_name(v)) << ": " - << f_[enum_integer(v)] << ",\n"; + << f_[static_cast(enum_integer(v))] << ",\n"; } void print(std::ostream& os, Json) const override @@ -209,7 +209,7 @@ namespace GridKit */ void print(std::ostream& os, VariableEnum v, Yaml) const { - os << indent_ << enum_name(v) << ": " << f_[enum_integer(v)] + os << indent_ << enum_name(v) << ": " << f_[static_cast(enum_integer(v))] << '\n'; } From ffe01f5687a283f6d72d885de3ab00cb222ff83a Mon Sep 17 00:00:00 2001 From: nkoukpaizan Date: Mon, 22 Dec 2025 18:35:06 +0000 Subject: [PATCH 06/14] Apply pre-commmit fixes --- GridKit/Model/PhasorDynamics/BusBase.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 9bb31d747..12dc04ed0 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -4,8 +4,8 @@ #include #include -#include #include +#include #include #include #include From 02410e00f82334978d65b5f7a0b66fbbf4d6e58f Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 22 Dec 2025 17:58:58 -0500 Subject: [PATCH 07/14] More intuitive PhasorDynamics example test names. --- .../PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt | 2 +- examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt | 2 +- examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt | 6 +++--- examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt | 6 +++--- examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt | 6 +++--- 5 files changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt b/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt index c000e4446..72a043616 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt +++ b/examples/PhasorDynamics/Small/TenGen/Classical/CMakeLists.txt @@ -5,4 +5,4 @@ target_link_libraries(TenGenClassical install(TARGETS TenGenClassical RUNTIME DESTINATION bin) # Not used for tesing for now. -add_test(NAME GenrouTest3_genclassical COMMAND $) +add_test(NAME TenGenClassical COMMAND $) diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt b/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt index b9baf8795..c5ee2d800 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/CMakeLists.txt @@ -5,4 +5,4 @@ target_link_libraries(TenGenGenrou install(TARGETS TenGenGenrou RUNTIME DESTINATION bin) # Not used for tesing for now. -add_test(NAME GenrouTest3 COMMAND $) +add_test(NAME TenGenGenrou COMMAND $) diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt index ee4ffb76e..aecb33407 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/CMakeLists.txt @@ -13,9 +13,9 @@ target_link_libraries(ThreeBusBasicJson install(TARGETS ThreeBusBasicJson RUNTIME DESTINATION ${_install_path}) gridkit_example_add_file(ThreeBusBasic.json) -add_test(NAME GenrouTest2 COMMAND ThreeBusBasic) -add_test(NAME GenrouTest2Json +add_test(NAME ThreeBusBasic COMMAND ThreeBusBasic) +add_test(NAME ThreeBusBasicJson COMMAND ThreeBusBasicJson ${CMAKE_CURRENT_BINARY_DIR}/ThreeBusBasic.json) -add_test(NAME GenrouTest2Json_no_arg +add_test(NAME ThreeBusBasicJson_no_arg COMMAND ThreeBusBasicJson WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt index a1b3895d6..8e514ab8a 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/CMakeLists.txt @@ -15,9 +15,9 @@ target_link_libraries(TwoBusIeeet1Json install(TARGETS TwoBusIeeet1Json RUNTIME DESTINATION ${_install_path}) gridkit_example_add_file(TwoBusIeeet1.json) -add_test(NAME GenrouTest1_Ieeet1 COMMAND TwoBusIeeet1) -add_test(NAME GenrouTest1_Ieeet1_Json +add_test(NAME TwoBusIeeet1 COMMAND TwoBusIeeet1) +add_test(NAME TwoBusIeeet1Json COMMAND TwoBusIeeet1Json ${CMAKE_CURRENT_BINARY_DIR}/TwoBusIeeet1.json) -add_test(NAME GenrouTest1_Ieeet1_Json_no_arg +add_test(NAME TwoBusIeeet1Json_no_arg COMMAND TwoBusIeeet1Json WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt index c0711d204..c6c8008e4 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/CMakeLists.txt @@ -15,9 +15,9 @@ target_link_libraries(TwoBusTgov1Json install(TARGETS TwoBusTgov1Json RUNTIME DESTINATION ${_install_path}) gridkit_example_add_file(TwoBusTgov1.json) -add_test(NAME GenrouTest1_tgov1 COMMAND TwoBusTgov1) -add_test(NAME GenrouTest1_tgov1_json +add_test(NAME TwoBusTgov1 COMMAND TwoBusTgov1) +add_test(NAME TwoBusTgov1Json COMMAND TwoBusTgov1Json ${CMAKE_CURRENT_BINARY_DIR}/TwoBusTgov1.json) -add_test(NAME GenrouTest1_tgov1_json_no_arg +add_test(NAME TwoBusTgov1Json_no_arg COMMAND TwoBusTgov1Json WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) From 1d3081702fcaf94d394b6abcd49da45bb6997739 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Thu, 8 Jan 2026 17:04:54 -0500 Subject: [PATCH 08/14] Add sstream include to VariableMonitor. --- GridKit/Model/VariableMonitor.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/GridKit/Model/VariableMonitor.hpp b/GridKit/Model/VariableMonitor.hpp index ce7df4f49..c50654aef 100644 --- a/GridKit/Model/VariableMonitor.hpp +++ b/GridKit/Model/VariableMonitor.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include From 59eba69ee9eac0c2268c10cdaed936d0c0b68fd2 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Thu, 8 Jan 2026 17:49:16 -0500 Subject: [PATCH 09/14] Missing LINK_LIBRARIES for Logger. --- GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt | 11 +++++++++-- GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt | 4 ++-- .../Model/PhasorDynamics/BusFault/CMakeLists.txt | 13 ++++++++++--- GridKit/Model/PhasorDynamics/Load/CMakeLists.txt | 11 +++++++++-- .../SynchronousMachine/GenClassical/CMakeLists.txt | 11 +++++++++-- 5 files changed, 39 insertions(+), 11 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt index 10761d6f9..3325ea604 100644 --- a/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Branch/CMakeLists.txt @@ -19,6 +19,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -27,13 +28,19 @@ else() Branch.cpp HEADERS ${_install_headers} + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) endif() gridkit_add_library(phasor_dynamics_branch_dependency_tracking - SOURCES BranchDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + BranchDependencyTracking.cpp + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) # Link to interface target for all components target_link_libraries(phasor_dynamics_components diff --git a/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt index b5086ed76..ca318a3b5 100644 --- a/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Bus/CMakeLists.txt @@ -38,8 +38,8 @@ gridkit_add_library(phasor_dynamics_bus_dependency_tracking SOURCES BusDependencyTracking.cpp BusInfiniteDependencyTracking.cpp - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PUBLIC GridKit::phasor_dynamics_core) diff --git a/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt b/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt index 40f438fb4..53b0de987 100644 --- a/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/BusFault/CMakeLists.txt @@ -13,6 +13,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -22,12 +23,18 @@ else() HEADERS ${_install_headers} INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core) endif() gridkit_add_library(phasor_dynamics_bus_fault_dependency_tracking - SOURCES BusFaultDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + BusFaultDependencyTracking.cpp + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core) # Link to interface target for all components target_link_libraries(phasor_dynamics_components diff --git a/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt index ef86891b4..91fb51641 100644 --- a/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Load/CMakeLists.txt @@ -19,6 +19,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -27,13 +28,19 @@ else() Load.cpp INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core HEADERS ${_install_headers}) endif() gridkit_add_library(phasor_dynamics_load_dependency_tracking - SOURCES LoadDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + LoadDependencyTracking.cpp + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core) # Link to interface target for all components target_link_libraries(phasor_dynamics_components diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt index ef1c782bc..ef03af11a 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/CMakeLists.txt @@ -18,6 +18,7 @@ if(GRIDKIT_ENABLE_ENZYME) PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES PRIVATE ClangEnzymeFlags + PUBLIC GridKit::phasor_dynamics_core COMPILE_OPTIONS PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) else() @@ -26,13 +27,19 @@ else() GenClassical.cpp HEADERS ${_install_headers} + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) endif() gridkit_add_library(phasor_dynamics_gen_classical_dependency_tracking - SOURCES GenClassicalDependencyTracking.cpp - INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) + SOURCES + GenClassicalDependencyTracking.cpp + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include) # Link to interface target for all components target_link_libraries(phasor_dynamics_components From f470a32788b83bcec1c6362ae658307243298c2c Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Mon, 12 Jan 2026 13:29:15 -0500 Subject: [PATCH 10/14] Hide initial Jacobian evaluation workaround in SystemModel a todo item for the sparsity pattern analysis. --- GridKit/Model/PhasorDynamics/SystemModel.hpp | 11 +++++++++++ .../Small/TenGen/Classical/TenGenClassical.cpp | 3 --- .../Small/TenGen/Genrou/TenGenGenrou.cpp | 3 --- .../Tiny/ThreeBus/Basic/ThreeBusBasic.cpp | 3 --- .../Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp | 3 --- .../Tiny/ThreeBus/Classical/ThreeBusClassical.cpp | 3 --- .../Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp | 3 --- .../PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp | 3 --- .../Tiny/TwoBus/Basic/TwoBusBasicJson.cpp | 3 --- .../Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp | 3 --- .../Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp | 3 --- .../PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp | 3 --- .../Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp | 3 --- 13 files changed, 11 insertions(+), 36 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 92cbac041..5db3cb5b3 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -345,6 +345,17 @@ namespace GridKit throw std::runtime_error("SystemModel allocation failed"); } + // Perform an initial Jacobian evaluation for sparse Jacobians, such that + // the dynamic solver can querry the NNZ value when it is configured. + // @todo Replace with a sparsity analysis that sets the NNZ and allocates the Jacobian + // without needing the Jacobian values. + if (hasJacobian()) + { + initialize(); + evaluateResidual(); + evaluateJacobian(); + } + return 0; } diff --git a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp index 374ca81b2..50c9dd2d2 100644 --- a/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Classical/TenGenClassical.cpp @@ -97,9 +97,6 @@ int main() sys.addComponent(&gen10); sys.addComponent(&fault); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); real_type dt = 1.0 / 4.0 / 60.0; diff --git a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp index a3ae0d41f..b3be48eb4 100644 --- a/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp +++ b/examples/PhasorDynamics/Small/TenGen/Genrou/TenGenGenrou.cpp @@ -96,9 +96,6 @@ int main() sys.addComponent(&gen10); sys.addComponent(&fault); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); real_type dt = 1.0 / 4.0 / 60.0; diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp index 467dc7543..baad11072 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasic.cpp @@ -200,9 +200,6 @@ int main() SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp index 952c7f326..225fb919c 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Basic/ThreeBusBasicJson.cpp @@ -124,9 +124,6 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp index 6cfb90456..c7660de6d 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassical.cpp @@ -177,9 +177,6 @@ int main() SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp index 774586b22..98ddc7350 100644 --- a/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp +++ b/examples/PhasorDynamics/Tiny/ThreeBus/Classical/ThreeBusClassicalJson.cpp @@ -124,9 +124,6 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to fault 0 auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp index 86f9ccaf7..117e6e06b 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasic.cpp @@ -96,9 +96,6 @@ int main() SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp index 021dd0650..a2961c704 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp @@ -76,9 +76,6 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp index 5f4d387c1..886aaae1c 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1.cpp @@ -143,9 +143,6 @@ int main() SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp index 61289c97b..4c17ae62b 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Ieeet1/TwoBusIeeet1Json.cpp @@ -76,9 +76,6 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp index 538d93076..b48935189 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1.cpp @@ -117,9 +117,6 @@ int main() SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); diff --git a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp index 046f5a606..2d825a710 100644 --- a/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp +++ b/examples/PhasorDynamics/Tiny/TwoBus/Tgov1/TwoBusTgov1Json.cpp @@ -80,9 +80,6 @@ int main(int argc, const char* argv[]) SystemModel sys(data); sys.allocate(); - sys.initialize(); - sys.evaluateResidual(); - sys.evaluateJacobian(); // Get access to the fault auto* fault = sys.getBusFault(0); From 296316359051f8d2618523793b1a789fdabe5363 Mon Sep 17 00:00:00 2001 From: nkoukpaizan Date: Mon, 12 Jan 2026 18:36:24 +0000 Subject: [PATCH 11/14] Apply pre-commmit fixes --- GridKit/Model/PhasorDynamics/SystemModel.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 5db3cb5b3..a2ce708ef 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -347,9 +347,9 @@ namespace GridKit // Perform an initial Jacobian evaluation for sparse Jacobians, such that // the dynamic solver can querry the NNZ value when it is configured. - // @todo Replace with a sparsity analysis that sets the NNZ and allocates the Jacobian + // @todo Replace with a sparsity analysis that sets the NNZ and allocates the Jacobian // without needing the Jacobian values. - if (hasJacobian()) + if (hasJacobian()) { initialize(); evaluateResidual(); From 30da7cc1e174611cc8e45ca42b9330a2ce404dee Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Tue, 20 Jan 2026 17:42:27 -0500 Subject: [PATCH 12/14] CMake guard for Enzyme -O0 -enzyme-auto-sparsity=1 bug. --- cmake/FindEnzyme.cmake | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/cmake/FindEnzyme.cmake b/cmake/FindEnzyme.cmake index fbfaa440a..ad19c20af 100644 --- a/cmake/FindEnzyme.cmake +++ b/cmake/FindEnzyme.cmake @@ -12,6 +12,15 @@ Author(s): ]] + +# Error if CMAKE_BUILD_TYPE is not Release or RelWithDebInfo due to an Enzyme bug +if (NOT ((CMAKE_BUILD_TYPE STREQUAL "Release") OR (CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo"))) + message(STATUS "CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}") + message(FATAL_ERROR "Enzyme builds currently only support Release and RelWithDebInfo as \ + CMAKE_BUILD_TYPE due to a bug within the Enzyme library.") +endif() + +# Find Enzyme and necessary programs find_package(Enzyme REQUIRED CONFIG PATHS ${ENZYME_DIR} From 6732dc38181250e08840f94bbf7225e6214dd89b Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Tue, 20 Jan 2026 20:29:01 -0500 Subject: [PATCH 13/14] Misc documentation of Enzyme utilities. --- .../AutomaticDifferentiation/Enzyme/DfDwb.hpp | 12 +++++------ .../AutomaticDifferentiation/Enzyme/DfDws.hpp | 6 +++--- .../AutomaticDifferentiation/Enzyme/DfDy.hpp | 20 +++++++++---------- .../AutomaticDifferentiation/Enzyme/DhDwb.hpp | 6 +++--- .../AutomaticDifferentiation/Enzyme/DhDy.hpp | 6 +++--- .../Enzyme/EnzymeDefinitions.hpp | 4 ++++ .../Enzyme/LowerSparseStorage.hpp | 16 ++++++++++++++- .../Enzyme/ModelWrappers.hpp | 8 ++++---- 8 files changed, 48 insertions(+), 30 deletions(-) diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp index f2a6abacf..392a44d40 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp @@ -55,11 +55,11 @@ namespace GridKit { if (n_res > 0 && n_var > 0) { - std::vector> triplets; + std::vector> triplets; //< @todo: Update once sparse storage format changes std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -67,7 +67,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, @@ -122,11 +122,11 @@ namespace GridKit { if (n_res > 0 && n_var > 0) { - std::vector> triplets; + std::vector> triplets; //< @todo: Update once sparse storage format changes std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -134,7 +134,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp index 0155193df..03591b919 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp @@ -58,11 +58,11 @@ namespace GridKit { if (n_res > 0 && n_var > 0) { - std::vector> triplets; + std::vector> triplets; //< @todo: Update once sparse storage format changes std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -70,7 +70,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp index ee2575ef9..3aac93d3f 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp @@ -58,11 +58,11 @@ namespace GridKit if (n_res > 0 && n_var > 0) { // df/dy - std::vector> triplets; + std::vector> triplets; //< @todo: Update once sparse storage format changes std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -70,7 +70,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, @@ -102,7 +102,7 @@ namespace GridKit triplets.clear(); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -110,7 +110,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, @@ -168,11 +168,11 @@ namespace GridKit if (n_res > 0 && n_var > 0) { // df/dy - std::vector> triplets; + std::vector> triplets; //< @todo: Update once sparse storage format changes std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -180,7 +180,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, @@ -214,7 +214,7 @@ namespace GridKit triplets.clear(); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -222,7 +222,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp index 87408de16..234657f4a 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDwb.hpp @@ -55,11 +55,11 @@ namespace GridKit { if (n_res > 0 && n_var > 0) { - std::vector> triplets; + std::vector> triplets; //< @todo: Update once sparse storage format changes std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -67,7 +67,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, diff --git a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp index 9e05f30d2..a0b609725 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/DhDy.hpp @@ -55,11 +55,11 @@ namespace GridKit { if (n_res > 0 && n_var > 0) { - std::vector> triplets; + std::vector> triplets; //< @todo: Update once sparse storage format changes std::vector elementary_v(n_var); for (size_t var_i = 0; var_i < n_var; ++var_i) { - // Sparse storage + // Sparse storage. @see LowerSparseStorage.hpp ScalarT* output = __enzyme_todense((void*) ident_load, (void*) ident_store, var_i); ScalarT* d_output = __enzyme_todense((void*) sparse_load, (void*) sparse_store, var_i, &triplets); @@ -67,7 +67,7 @@ namespace GridKit std::ranges::fill(elementary_v, 0.0); elementary_v[var_i] = 1.0; - // Autodiff + // Core automatic differentiaation intrinsic that will be replaced by a derivative __enzyme_fwddiff((void*) ModelWrapper::eval, enzyme_const, model, diff --git a/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp b/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp index 858385ce7..4edbcd88e 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp @@ -28,6 +28,10 @@ namespace GridKit /** * @brief Enzyme fwddiff template * + * @details This is a core templated intrinsic that the Enzyme pass will + * use to perform automatic differenciation. We define the template here + * so it can later be used in different places. + * * @tparam T - return type * @tparam ModelT - model type */ diff --git a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp index 564584077..646bbf796 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp @@ -1,5 +1,10 @@ /** * @file LowerSparseStorage.hpp + * + * @details This file contains functions used by Enzyme to store sparse Jacobians. + * + * @todo Replace this inner/lower storage by a more CSR-efficient approach. + * * @author Nicholson Koukpaizan (koukpaizannk@ornl.gov) * */ @@ -17,13 +22,16 @@ namespace GridKit /** * @brief Enzyme todense template * + * @details This is used by Enzyme's auto sparsity analysis. It internally maps dense storage + * to sparse ones (where structural zeros are not kept). + * * @tparam T - return type */ template extern T __enzyme_todense(void*...) noexcept; /** - * @brief Enzyme sparse storage in triplet format + * @brief Enzyme inner sparse storage in triplet format * * @tparam ScalarT - scalar data type */ @@ -46,6 +54,8 @@ namespace GridKit /** * @brief Enzyme sparse accumulation for float * + * @todo Separate sparsity analyis from value storage for a known sparsity structure. + * The current implementation always performs sparsity analysis. */ [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(size_t row, size_t col, float val, std::vector>& triplets) { @@ -55,6 +65,8 @@ namespace GridKit /** * @brief Enzyme sparse accumulation for double * + * @todo Separate sparsity analyis from value storage for a known sparsity structure. + * The current implementation always performs sparsity analysis. */ [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(size_t row, size_t col, double val, std::vector>& triplets) { @@ -64,6 +76,8 @@ namespace GridKit /** * @brief Enzyme sparse store * + * @details This takes in a row, column and value and stores them in a vector of triplets + * * @tparam ScalarT - scalar data type */ template diff --git a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp index aff23d92a..09d67f08d 100644 --- a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp +++ b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp @@ -104,7 +104,7 @@ namespace GridKit }; /** - * @brief Residual wrapper partial template specialization for BusResidual11 + * @brief Residual wrapper partial template specialization for BusResidual11 (branch member function) * */ template @@ -124,7 +124,7 @@ namespace GridKit }; /** - * @brief Residual wrapper partial template specialization for BusResidual12 + * @brief Residual wrapper partial template specialization for BusResidual12 (branch member function) * */ template @@ -144,7 +144,7 @@ namespace GridKit }; /** - * @brief Residual wrapper partial template specialization for BusResidual21 + * @brief Residual wrapper partial template specialization for BusResidual21 (branch member function) * */ template @@ -165,7 +165,7 @@ namespace GridKit }; /** - * @brief Residual wrapper partial template specialization for BusResidual22 + * @brief Residual wrapper partial template specialization for BusResidual22 (branch member function) * */ template From 77080fee9505d9096bff30d6ad8cad6f36b44a29 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Tue, 20 Jan 2026 20:37:25 -0500 Subject: [PATCH 14/14] Fix warning. --- GridKit/Model/PowerElectronics/CircuitComponent.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index 55bcf3f17..cbf8a6135 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -40,7 +40,7 @@ namespace GridKit this->alpha_ = a; } - bool hasJacobian() + bool hasJacobian() override { return true; }