Skip to content

Commit bc14bf9

Browse files
nkoukpaizanSteven-Roberts
nkoukpaizan
and
Steven-Roberts
authored
System-level Jacobian assembly (#229)
* Functional system level Jacobian for internal variables/residuals and bus variables/residuals. * Bug fix after rebasing. * Expand documentation of Enzyme::Sparse::Wrapper/Jacobian * Use struct for Enzyme::Sparse::*Jacobian. * Documentation updates. * Revert unnecessary, distracting changes. * Add headers for *Enzyme.cpp files. * Swap res_i and var_i in Enzyme::Sparse::Wrapper. * Fix elementary_v. * Use ranges::fill instead of manual loop Co-authored-by: Steven Roberts <[email protected]> * Specify inline in function definition as well Co-authored-by: Steven Roberts <[email protected]> * Document variable_indices and residual_indices. --------- Co-authored-by: nkoukpaizan <[email protected]> Co-authored-by: Steven Roberts <[email protected]>
1 parent 5de1437 commit bc14bf9

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

48 files changed

+1660
-532
lines changed

‎src/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp‎

Lines changed: 286 additions & 43 deletions
Large diffs are not rendered by default.

‎src/LinearAlgebra/SparseMatrix/COO_Matrix.hpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -54,8 +54,10 @@ namespace GridKit
5454
std::tuple<std::vector<IdxT>, std::vector<IdxT>, std::vector<ScalarT>> setDataToCSR();
5555
std::vector<IdxT> getCSRRowData();
5656

57+
// Set values from vector storage. Will sort before storing
58+
void setValues(std::vector<IdxT> r, std::vector<IdxT> c, std::vector<ScalarT> v);
59+
5760
// BLAS. Will sort before running
58-
void setValues(std::vector<IdxT> r, std::vector<IdxT> c, std::vector<ScalarT> v);
5961
void axpy(ScalarT alpha, COO_Matrix<ScalarT, IdxT>& a);
6062
void axpy(ScalarT alpha, std::vector<IdxT> r, std::vector<IdxT> c, std::vector<ScalarT> v);
6163
void scal(ScalarT alpha);
@@ -66,8 +68,8 @@ namespace GridKit
6668
void permutation(std::vector<IdxT> row_perm, std::vector<IdxT> col_perm);
6769
void permutationSizeMap(std::vector<IdxT> row_perm, std::vector<IdxT> col_perm, IdxT m, IdxT n);
6870

71+
// Special matrices (zero and identity)
6972
void zeroMatrix();
70-
7173
void identityMatrix(IdxT n);
7274

7375
// Resort values_

‎src/Model/PhasorDynamics/Branch/Branch.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,22 @@ namespace GridKit
1010
{
1111
namespace PhasorDynamics
1212
{
13+
/**
14+
* @brief Jacobian evaluation not implemented
15+
*
16+
* @tparam ScalarT - scalar data type
17+
* @tparam IdxT - matrix index data type
18+
* @return int - error code, 0 = success
19+
*/
20+
template <class ScalarT, typename IdxT>
21+
int Branch<ScalarT, IdxT>::evaluateJacobian()
22+
{
23+
std::cout << "Evaluate Jacobian for Branch..." << std::endl;
24+
std::cout << "Jacobian evaluation is not implemented!" << std::endl;
25+
26+
return 0;
27+
}
28+
1329
// Available template instantiations
1430
template class Branch<double, long int>;
1531
template class Branch<double, size_t>;

‎src/Model/PhasorDynamics/Branch/Branch.hpp‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,9 @@ namespace GridKit
4747
using Component<ScalarT, IdxT>::yp_;
4848
using Component<ScalarT, IdxT>::tag_;
4949
using Component<ScalarT, IdxT>::f_;
50+
using Component<ScalarT, IdxT>::w_;
51+
using Component<ScalarT, IdxT>::h_;
52+
using Component<ScalarT, IdxT>::J_;
5053

5154
using real_type = typename Component<ScalarT, IdxT>::real_type;
5255
using bus_type = BusBase<ScalarT, IdxT>;
@@ -73,25 +76,31 @@ namespace GridKit
7376
void setR(real_type R)
7477
{
7578
R_ = R;
79+
setDerivedParams();
7680
}
7781

7882
void setX(real_type X)
7983
{
8084
// std::cout << "Setting X ...\n";
8185
X_ = X;
86+
setDerivedParams();
8287
}
8388

8489
void setG(real_type G)
8590
{
8691
G_ = G;
92+
setDerivedParams();
8793
}
8894

8995
void setB(real_type B)
9096
{
9197
B_ = B;
98+
setDerivedParams();
9299
}
93100

94101
private:
102+
void setDerivedParams();
103+
95104
ScalarT& Vr1()
96105
{
97106
return bus1_->Vr();
@@ -132,6 +141,12 @@ namespace GridKit
132141
return bus2_->Ii();
133142
}
134143

144+
public:
145+
__attribute__((always_inline)) inline int evaluateBusResidual11(ScalarT*, ScalarT*, ScalarT*, ScalarT*);
146+
__attribute__((always_inline)) inline int evaluateBusResidual12(ScalarT*, ScalarT*, ScalarT*, ScalarT*);
147+
__attribute__((always_inline)) inline int evaluateBusResidual21(ScalarT*, ScalarT*, ScalarT*, ScalarT*);
148+
__attribute__((always_inline)) inline int evaluateBusResidual22(ScalarT*, ScalarT*, ScalarT*, ScalarT*);
149+
135150
private:
136151
bus_type* bus1_;
137152
bus_type* bus2_;
@@ -141,6 +156,10 @@ namespace GridKit
141156
real_type B_{0.0};
142157
IdxT bus1_id_{0};
143158
IdxT bus2_id_{0};
159+
160+
/* Derivied parameters */
161+
real_type b_;
162+
real_type g_;
144163
};
145164

146165
} // namespace PhasorDynamics

‎src/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,22 @@ namespace GridKit
1010
{
1111
namespace PhasorDynamics
1212
{
13+
/**
14+
* @brief Jacobian evaluation not implemented
15+
*
16+
* @tparam ScalarT - scalar data type
17+
* @tparam IdxT - matrix index data type
18+
* @return int - error code, 0 = success
19+
*/
20+
template <class ScalarT, typename IdxT>
21+
int Branch<ScalarT, IdxT>::evaluateJacobian()
22+
{
23+
std::cout << "Evaluate Jacobian for Branch..." << std::endl;
24+
std::cout << "Jacobian evaluation is not implemented!" << std::endl;
25+
26+
return 0;
27+
}
28+
1329
// Available template instantiations
1430
template class Branch<DependencyTracking::Variable, long int>;
1531
template class Branch<DependencyTracking::Variable, size_t>;

‎src/Model/PhasorDynamics/Branch/BranchEnzyme.cpp‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
/**
2+
* @file BranchEnzyme.cpp
3+
* @author Nicholson Koukpaizan ([email protected])
4+
*
5+
*/
6+
7+
#include "BranchImpl.hpp"
8+
#include <AutomaticDifferentiation/Enzyme/SparseWrapper.hpp>
9+
10+
namespace GridKit
11+
{
12+
namespace PhasorDynamics
13+
{
14+
/**
15+
* @brief Jacobian evaluation experimental
16+
*
17+
* @tparam ScalarT - scalar data type
18+
* @tparam IdxT - matrix index data type
19+
* @return int - error code, 0 = success
20+
*/
21+
template <class ScalarT, typename IdxT>
22+
int Branch<ScalarT, IdxT>::evaluateJacobian()
23+
{
24+
std::cout << "Evaluate Jacobian for Branch..." << std::endl;
25+
std::cout << "Jacobian evaluation is experimental!" << std::endl;
26+
27+
// Bus 1 diagonal Jacobian block owned by the bus
28+
GridKit::Enzyme::Sparse::BusJacobian<GridKit::PhasorDynamics::Branch<ScalarT, IdxT>,
29+
GridKit::Enzyme::Sparse::MemberFunctions::BusResidual11,
30+
ScalarT,
31+
IdxT>::eval(this,
32+
static_cast<size_t>(bus1_->size()),
33+
static_cast<size_t>((bus1_->y()).size()),
34+
bus1_->getResidualIndices(),
35+
bus1_->getVariableIndices(),
36+
y_.data(),
37+
yp_.data(),
38+
(bus1_->y()).data(),
39+
bus1_->getJacobian());
40+
41+
// Bus 2 diagonal Jacobian block owned by the bus
42+
GridKit::Enzyme::Sparse::BusJacobian<GridKit::PhasorDynamics::Branch<ScalarT, IdxT>,
43+
GridKit::Enzyme::Sparse::MemberFunctions::BusResidual22,
44+
ScalarT,
45+
IdxT>::eval(this,
46+
static_cast<size_t>(bus2_->size()),
47+
static_cast<size_t>((bus2_->y()).size()),
48+
bus2_->getResidualIndices(),
49+
bus2_->getVariableIndices(),
50+
y_.data(),
51+
yp_.data(),
52+
(bus2_->y()).data(),
53+
bus2_->getJacobian());
54+
55+
// Off-diagonal Jacobian block (Bus2 variables) owned by the branch
56+
GridKit::Enzyme::Sparse::BranchJacobian<GridKit::PhasorDynamics::Branch<ScalarT, IdxT>,
57+
GridKit::Enzyme::Sparse::MemberFunctions::BusResidual12,
58+
ScalarT,
59+
IdxT>::eval(this,
60+
static_cast<size_t>(bus1_->size()),
61+
static_cast<size_t>((bus2_->y()).size()),
62+
bus1_->getResidualIndices(),
63+
bus2_->getVariableIndices(),
64+
y_.data(),
65+
yp_.data(),
66+
(bus2_->y()).data(),
67+
J_);
68+
69+
J_.printMatrix("Branch after BusJacobian12 call");
70+
71+
// Off-diagonal Jacobian block (Bus1 variables) owned by the branch
72+
GridKit::Enzyme::Sparse::BranchJacobian<GridKit::PhasorDynamics::Branch<ScalarT, IdxT>,
73+
GridKit::Enzyme::Sparse::MemberFunctions::BusResidual21,
74+
ScalarT,
75+
IdxT>::eval(this,
76+
static_cast<size_t>(bus2_->size()),
77+
static_cast<size_t>((bus1_->y()).size()),
78+
bus2_->getResidualIndices(),
79+
bus1_->getVariableIndices(),
80+
y_.data(),
81+
yp_.data(),
82+
(bus1_->y()).data(),
83+
J_);
84+
85+
J_.printMatrix("Branch after BusJacobian21 call");
86+
87+
return 0;
88+
}
89+
90+
// Available template instantiations
91+
template class Branch<double, long int>;
92+
template class Branch<double, size_t>;
93+
94+
} // namespace PhasorDynamics
95+
} // namespace GridKit

0 commit comments

Comments
 (0)