Skip to content

PhasorDynamics simulations with sparse Jacobians #312

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pelesh merged 14 commits into from
Jan 21, 2026
Merged

PhasorDynamics simulations with sparse Jacobians #312

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
cf71e3b
Functional PhasorDynamics simulations with sparse Jacobians.
nkoukpaizan Nov 24, 2025
e49d7b3
Cleaner Enzyme sparse wrapper implementation.
nkoukpaizan Dec 21, 2025
d70afce
Replace int-int maps with vector of integers for residual/variable in…
nkoukpaizan Dec 22, 2025
d0763c9
Make use of INVALID_INDEX constant.
nkoukpaizan Dec 22, 2025
0429a10
Suppress warnings.
nkoukpaizan Dec 22, 2025
ffe01f5
Apply pre-commmit fixes
nkoukpaizan Dec 22, 2025
02410e0
More intuitive PhasorDynamics example test names.
nkoukpaizan Dec 22, 2025
1d30817
Add sstream include to VariableMonitor.
nkoukpaizan Jan 8, 2026
59eba69
Missing LINK_LIBRARIES for Logger.
nkoukpaizan Jan 8, 2026
f470a32
Hide initial Jacobian evaluation workaround in SystemModel a todo ite…
nkoukpaizan Jan 12, 2026
2963163
Apply pre-commmit fixes
nkoukpaizan Jan 12, 2026
30da7cc
CMake guard for Enzyme -O0 -enzyme-auto-sparsity=1 bug.
nkoukpaizan Jan 20, 2026
6732dc3
Misc documentation of Enzyme utilities.
nkoukpaizan Jan 21, 2026
77080fe
Fix warning.
nkoukpaizan Jan 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
3 changes: 2 additions & 1 deletion GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include <string>
#include <vector>

#include <GridKit/Constants.hpp>
#include <GridKit/ScalarTraits.hpp>

namespace GridKit
Expand Down Expand Up @@ -266,7 +267,7 @@ namespace GridKit
bool is_fixed_; ///< Constant parameter flag.

mutable DependencyMap* dependencies_;
static const size_t INVALID_VAR_NUMBER = static_cast<size_t>(-1);
static const size_t INVALID_VAR_NUMBER = INVALID_INDEX<size_t>;
};

//------------------------------------
Expand Down
171 changes: 171 additions & 0 deletions GridKit/AutomaticDifferentiation/Enzyme/DfDwb.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
/**
* @file DfDwb.hpp
* @author Nicholson Koukpaizan (koukpaizannk@ornl.gov)
*
*/

#pragma once

#include <GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp>
#include <GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp>
#include <GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp>
#include <GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp>
#include <GridKit/ScalarTraits.hpp>

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 <typename ModelT, MemberFunctions function, class ScalarT, typename IdxT>
struct DfDwb
{
using RealT = typename GridKit::ScalarTraits<ScalarT>::RealT;
using MatrixT = GridKit::LinearAlgebra::COO_Matrix<RealT, IdxT>;

/**
* @param[in] model - Pointer to the model to be differentiated
* @param[in] n_res - Number of residual functions
* @param[in] n_var - Number of independent variables
* @param[in] res_indices - Global residual indices
* @param[in] var_indices - Global variable indices
* @param[in] y - Internal variables
* @param[in] yp - Internal variable derivatives
* @param[in] wb - Bus variables
* @param[in,out] jac - Jacobian
*/
static void eval(ModelT* model,
size_t n_res,
size_t n_var,
const std::vector<IdxT>& res_indices,
const std::vector<IdxT>& var_indices,
ScalarT* y,
ScalarT* yp,
ScalarT* wb,
MatrixT& jac)
{
if (n_res > 0 && n_var > 0)
{
std::vector<Triple<ScalarT>> triplets; //< @todo: Update once sparse storage format changes
std::vector<ScalarT> elementary_v(n_var);
for (size_t var_i = 0; var_i < n_var; ++var_i)
{
// Sparse storage. @see LowerSparseStorage.hpp
ScalarT* output = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT>, (void*) ident_store<ScalarT>, var_i);
ScalarT* d_output = __enzyme_todense<ScalarT*>((void*) sparse_load<ScalarT>, (void*) sparse_store<ScalarT>, var_i, &triplets);

// Elementary vector for Jacobian-vector product
std::ranges::fill(elementary_v, 0.0);
elementary_v[var_i] = 1.0;

// Core automatic differentiaation intrinsic that will be replaced by a derivative
__enzyme_fwddiff<void>((void*) ModelWrapper<ModelT, function, ScalarT>::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<IdxT> ctemp{};
std::vector<IdxT> rtemp{};
std::vector<ScalarT> valtemp{};
for (auto& tup : triplets)
{
rtemp.push_back(res_indices[static_cast<size_t>(tup.row)]);
ctemp.push_back(var_indices[static_cast<size_t>(tup.col)]);
valtemp.push_back(tup.val);
}
jac.setValues(rtemp, ctemp, valtemp); //< @todo: Update once sparse storage format changes
}
}

/**
* @param[in] model - Pointer to the model to be differentiated
* @param[in] n_res - Number of residual functions
* @param[in] n_var - Number of independent variables
* @param[in] res_indices - Global residual indices
* @param[in] var_indices - Global variable indices
* @param[in] y - Internal variables
* @param[in] yp - Internal variable derivatives
* @param[in] wb - Bus variables
* @param[in] ws - Signal variables
* @param[in,out] jac - Jacobian
*/
static void eval(ModelT* model,
size_t n_res,
size_t n_var,
const std::vector<IdxT>& res_indices,
const std::vector<IdxT>& var_indices,
ScalarT* y,
ScalarT* yp,
ScalarT* wb,
ScalarT* ws,
MatrixT& jac)
{
if (n_res > 0 && n_var > 0)
{
std::vector<Triple<ScalarT>> triplets; //< @todo: Update once sparse storage format changes
std::vector<ScalarT> elementary_v(n_var);
for (size_t var_i = 0; var_i < n_var; ++var_i)
{
// Sparse storage. @see LowerSparseStorage.hpp
ScalarT* output = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT>, (void*) ident_store<ScalarT>, var_i);
ScalarT* d_output = __enzyme_todense<ScalarT*>((void*) sparse_load<ScalarT>, (void*) sparse_store<ScalarT>, var_i, &triplets);

// Elementary vector for Jacobian-vector product
std::ranges::fill(elementary_v, 0.0);
elementary_v[var_i] = 1.0;

// Core automatic differentiaation intrinsic that will be replaced by a derivative
__enzyme_fwddiff<void>((void*) ModelWrapper<ModelT, function, ScalarT>::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<IdxT> ctemp{};
std::vector<IdxT> rtemp{};
std::vector<ScalarT> valtemp{};
for (auto& tup : triplets)
{
rtemp.push_back(res_indices[static_cast<size_t>(tup.row)]);
ctemp.push_back(var_indices[static_cast<size_t>(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
112 changes: 112 additions & 0 deletions GridKit/AutomaticDifferentiation/Enzyme/DfDws.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,112 @@
/**
* @file DfDws.hpp
* @author Nicholson Koukpaizan (koukpaizannk@ornl.gov)
*
*/

#pragma once

#include <GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp>
#include <GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp>
#include <GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp>
#include <GridKit/Constants.hpp>
#include <GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp>
#include <GridKit/ScalarTraits.hpp>

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 <typename ModelT, MemberFunctions function, class ScalarT, typename IdxT>
struct DfDws
{
using RealT = typename GridKit::ScalarTraits<ScalarT>::RealT;
using MatrixT = GridKit::LinearAlgebra::COO_Matrix<RealT, IdxT>;

/**
* @param[in] model - Pointer to the model to be differentiated
* @param[in] n_res - Number of residual functions
* @param[in] n_var - Number of independent variables
* @param[in] res_indices - Global residual indices
* @param[in] var_indices - Global variable indices
* @param[in] y - Internal variables
* @param[in] yp - Internal variable derivatives
* @param[in] wb - Bus variables
* @param[in] ws - Signal variables
* @param[in,out] jac - Jacobian
*/
static void eval(ModelT* model,
size_t n_res,
size_t n_var,
const std::vector<IdxT>& res_indices,
const std::vector<IdxT>& var_indices,
ScalarT* y,
ScalarT* yp,
ScalarT* wb,
ScalarT* ws,
MatrixT& jac)
{
if (n_res > 0 && n_var > 0)
{
std::vector<Triple<ScalarT>> triplets; //< @todo: Update once sparse storage format changes
std::vector<ScalarT> elementary_v(n_var);
for (size_t var_i = 0; var_i < n_var; ++var_i)
{
// Sparse storage. @see LowerSparseStorage.hpp
ScalarT* output = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT>, (void*) ident_store<ScalarT>, var_i);
ScalarT* d_output = __enzyme_todense<ScalarT*>((void*) sparse_load<ScalarT>, (void*) sparse_store<ScalarT>, var_i, &triplets);

// Elementary vector for Jacobian-vector product
std::ranges::fill(elementary_v, 0.0);
elementary_v[var_i] = 1.0;

// Core automatic differentiaation intrinsic that will be replaced by a derivative
__enzyme_fwddiff<void>((void*) ModelWrapper<ModelT, function, ScalarT>::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<IdxT> ctemp{};
std::vector<IdxT> rtemp{};
std::vector<ScalarT> valtemp{};
for (auto& tup : triplets)
{
IdxT row = res_indices[static_cast<size_t>(tup.row)];
IdxT col = var_indices[static_cast<size_t>(tup.col)];
if (col != INVALID_INDEX<IdxT>)
{
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
Loading
Loading