Commit db688d10 authored by Philipp Fensch's avatar Philipp Fensch Committed by Markus Mirz
Browse files

Added option for using sparse matrices in MNA-Solver

parent a1992f5d
......@@ -138,9 +138,10 @@ elseif("${CMAKE_SYSTEM}" MATCHES "Darwin")
set(MacOS_FOUND ON)
endif()
option(BUILD_SHARED_LIBS "Build shared library" OFF)
option(BUILD_EXAMPLES "Build C++ examples" ON)
option(GET_GRID_DATA "Download grid data" ON)
option(BUILD_SHARED_LIBS "Build shared library" OFF)
option(BUILD_EXAMPLES "Build C++ examples" ON )
option(GET_GRID_DATA "Download grid data" ON )
option(WITH_SPARSE "Use sparse matrices in MNA-Solver" ON )
option(PYBIND "Build pybind module" OFF)
if(PYBIND)
......@@ -148,15 +149,15 @@ if(PYBIND)
endif()
include(CMakeDependentOption)
cmake_dependent_option(WITH_GSL "Enable GSL" ON "GSL_FOUND" OFF)
cmake_dependent_option(WITH_SUNDIALS "Enable sundials solver suite" ON "Sundials_FOUND" OFF)
cmake_dependent_option(WITH_SHMEM "Enable shared memory interface" ON "VILLASnode_FOUND" OFF)
cmake_dependent_option(WITH_RT "Enable real-time features" ON "Linux_FOUND" OFF)
cmake_dependent_option(WITH_PYTHON "Enable Python support" ON "Python_FOUND" OFF)
cmake_dependent_option(WITH_CIM "Enable support for parsing CIM" ON "CIMpp_FOUND" OFF)
cmake_dependent_option(WITH_OPENMP "Enable OpenMP-based parallelisation" ON "OPENMP_FOUND" OFF)
cmake_dependent_option(WITH_CUDA "Enable CUDA-based parallelisation" ON "CUDA_FOUND" OFF)
cmake_dependent_option(WITH_GRAPHVIZ "Enable Graphviz Graphs" ON "GRAPHVIZ_FOUND" OFF)
cmake_dependent_option(WITH_GSL "Enable GSL" ON "GSL_FOUND" OFF)
cmake_dependent_option(WITH_SUNDIALS "Enable sundials solver suite" ON "Sundials_FOUND" OFF)
cmake_dependent_option(WITH_SHMEM "Enable shared memory interface" ON "VILLASnode_FOUND" OFF)
cmake_dependent_option(WITH_RT "Enable real-time features" ON "Linux_FOUND" OFF)
cmake_dependent_option(WITH_PYTHON "Enable Python support" ON "Python_FOUND" OFF)
cmake_dependent_option(WITH_CIM "Enable support for parsing CIM" ON "CIMpp_FOUND" OFF)
cmake_dependent_option(WITH_OPENMP "Enable OpenMP-based parallelisation" ON "OPENMP_FOUND" OFF)
cmake_dependent_option(WITH_CUDA "Enable CUDA-based parallelisation" OFF "CUDA_FOUND" OFF)
cmake_dependent_option(WITH_GRAPHVIZ "Enable Graphviz Graphs" ON "GRAPHVIZ_FOUND" OFF)
if(WITH_CUDA)
# BEGIN OF WORKAROUND - enable cuda dynamic linking.
......
......@@ -25,6 +25,7 @@
#cmakedefine WITH_SUNDIALS
#cmakedefine WITH_OPENMP
#cmakedefine WITH_CUDA
#cmakedefine WITH_SPARSE
#cmakedefine CGMES_BUILD
#cmakedefine HAVE_TIMERFD
......
......@@ -20,6 +20,8 @@ namespace DPsim {
using UInt = CPS::UInt;
using Matrix = CPS::Matrix;
using MatrixComp = CPS::MatrixComp;
using SparseMatrix = CPS::SparseMatrixRow;
using SparseMatrixComp = CPS::SparseMatrixCompRow;
template<typename T>
using MatrixVar = CPS::MatrixVar<T>;
......
......@@ -14,6 +14,7 @@
#include <unordered_map>
#include <bitset>
#include <dpsim/Config.h>
#include <dpsim/Solver.h>
#include <dpsim/DataLogger.h>
#include <cps/AttributeList.h>
......@@ -30,6 +31,12 @@
**/
#define SWITCH_NUM sizeof(std::size_t)*8
#ifdef WITH_SPARSE
#define MAT_TYPE SparseMatrix
#else
#define MAT_TYPE Matrix
#endif
namespace DPsim {
/// Solver class using Modified Nodal Analysis (MNA).
template <typename VarType>
......@@ -86,12 +93,16 @@ namespace DPsim {
std::vector< CPS::Attribute<Matrix>::Ptr > mLeftVectorHarmAttributes;
/// Base matrix that includes all static MNA elements to speed up recomputation
Matrix mBaseSystemMatrix;
MAT_TYPE mBaseSystemMatrix;
/// Map of system matrices where the key is the bitset describing the switch states
std::unordered_map< std::bitset<SWITCH_NUM>, Matrix > mSwitchedMatrices;
std::unordered_map< std::bitset<SWITCH_NUM>, MAT_TYPE > mSwitchedMatrices;
std::unordered_map< std::bitset<SWITCH_NUM>, std::vector<Matrix> > mSwitchedMatricesHarm;
#ifdef WITH_SPARSE
/// Map of LU factorizations related to the system matrices
std::unordered_map< std::bitset<SWITCH_NUM>, CPS::LUFactorizedSparse > mLuFactorizations;
#else
std::unordered_map< std::bitset<SWITCH_NUM>, CPS::LUFactorized > mLuFactorizations;
#endif
std::unordered_map< std::bitset<SWITCH_NUM>, std::vector<CPS::LUFactorized> > mLuFactorizationsHarm;
// #### Attributes related to switching ####
......@@ -164,7 +175,9 @@ namespace DPsim {
///
virtual CPS::Task::List getTasks();
///
Matrix& systemMatrix() { return mSwitchedMatrices[mCurrentSwitchStatus]; }
MAT_TYPE& systemMatrix() {
return mSwitchedMatrices[mCurrentSwitchStatus];
}
// #### MNA Solver Tasks ####
///
......
......@@ -214,7 +214,13 @@ void MnaSolver<VarType>::initializeSystemWithPrecomputedMatrices() {
if (mSwitches.size() < 1) {
// Create system matrix if no switches were added
for (auto comp : mMNAComponents) {
#ifdef WITH_SPARSE
Matrix mat = Matrix(mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
comp->mnaApplySystemMatrixStamp(mat);
mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)] = mat.sparseView();
#else
comp->mnaApplySystemMatrixStamp(mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#endif
auto idObj = std::dynamic_pointer_cast<IdentifiedObject>(comp);
mSLog->debug("Stamping {:s} {:s} into system matrix",
idObj->type(), idObj->name());
......@@ -223,17 +229,38 @@ void MnaSolver<VarType>::initializeSystemWithPrecomputedMatrices() {
Logger::matrixToString(mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]));
}
}
#ifdef WITH_SPARSE
mLuFactorizations[std::bitset<SWITCH_NUM>(0)].analyzePattern(mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
mLuFactorizations[std::bitset<SWITCH_NUM>(0)].factorize(mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#else
mLuFactorizations[std::bitset<SWITCH_NUM>(0)] = Eigen::PartialPivLU<Matrix>(mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#endif
}
else {
// Generate switching state dependent system matrices
for (auto& sys : mSwitchedMatrices) {
#ifdef WITH_SPARSE
for (auto comp : mMNAComponents) {
Matrix mat = Matrix(sys.second);
comp->mnaApplySystemMatrixStamp(mat);
sys.second = mat.sparseView();
}
for (UInt i = 0; i < mSwitches.size(); i++) {
Matrix mat = Matrix(sys.second);
mSwitches[i]->mnaApplySwitchSystemMatrixStamp(mat, sys.first[i]);
sys.second = mat.sparseView();
}
// Compute LU-factorization for system matrix
mLuFactorizations[sys.first].analyzePattern(sys.second);
mLuFactorizations[sys.first].factorize(sys.second);
#else
for (auto comp : mMNAComponents)
comp->mnaApplySystemMatrixStamp(sys.second);
for (UInt i = 0; i < mSwitches.size(); i++)
mSwitches[i]->mnaApplySwitchSystemMatrixStamp(sys.second, sys.first[i]);
// Compute LU-factorization for system matrix
mLuFactorizations[sys.first] = Eigen::PartialPivLU<Matrix>(sys.second);
#endif
}
updateSwitchStatus();
}
......@@ -347,10 +374,17 @@ void MnaSolver<Real>::createEmptySystemMatrix() {
throw SystemError("Too many Switches.");
for (std::size_t i = 0; i < (1ULL << mSwitches.size()); i++) {
#ifdef WITH_SPARSE
mSwitchedMatrices[std::bitset<SWITCH_NUM>(i)].resize(mNumMatrixNodeIndices, mNumMatrixNodeIndices);
#else
mSwitchedMatrices[std::bitset<SWITCH_NUM>(i)] = Matrix::Zero(mNumMatrixNodeIndices, mNumMatrixNodeIndices);
#endif
}
#ifdef WITH_SPARSE
mBaseSystemMatrix.resize(mNumMatrixNodeIndices, mNumMatrixNodeIndices);
#else
mBaseSystemMatrix = Matrix::Zero(mNumMatrixNodeIndices, mNumMatrixNodeIndices);
#endif
}
template<>
......@@ -368,11 +402,18 @@ void MnaSolver<Complex>::createEmptySystemMatrix() {
}
else {
for (std::size_t i = 0; i < (1ULL << mSwitches.size()); i++) {
#ifdef WITH_SPARSE
mSwitchedMatrices[std::bitset<SWITCH_NUM>(i)].resize(2*(mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices), 2*(mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices));
#else
mSwitchedMatrices[std::bitset<SWITCH_NUM>(i)] = Matrix::Zero(2*(mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices), 2*(mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices));
#endif
}
}
#ifdef WITH_SPARSE
mBaseSystemMatrix.resize(2 * (mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices), 2 * (mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices));
#else
mBaseSystemMatrix = Matrix::Zero(2 * (mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices), 2 * (mNumMatrixNodeIndices + mNumHarmMatrixNodeIndices));
#endif
}
template <typename VarType>
......@@ -598,8 +639,8 @@ void MnaSolver<VarType>::logSystemMatrices() {
for (UInt i = 0; i < mSwitchedMatricesHarm[std::bitset<SWITCH_NUM>(0)].size(); i++) {
mSLog->info("System matrix for frequency: {:d} \n{:s}", i,
Logger::matrixToString(mSwitchedMatricesHarm[std::bitset<SWITCH_NUM>(0)][i]));
mSLog->info("LU decomposition for frequency: {:d} \n{:s}", i,
Logger::matrixToString(mLuFactorizationsHarm[std::bitset<SWITCH_NUM>(0)][i].matrixLU()));
//mSLog->info("LU decomposition for frequency: {:d} \n{:s}", i,
// Logger::matrixToString(mLuFactorizationsHarm[std::bitset<SWITCH_NUM>(0)][i].matrixLU()));
}
for (UInt i = 0; i < mRightSideVectorHarm.size(); i++)
......@@ -610,7 +651,7 @@ void MnaSolver<VarType>::logSystemMatrices() {
else {
if (mSwitches.size() < 1) {
mSLog->info("System matrix: \n{}", mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
mSLog->info("LU decomposition: \n{}", mLuFactorizations[std::bitset<SWITCH_NUM>(0)].matrixLU());
//mSLog->info("LU decomposition: \n{}", mLuFactorizations[std::bitset<SWITCH_NUM>(0)].matrixLU());
}
else {
mSLog->info("Initial switch status: {:s}", mCurrentSwitchStatus.to_string());
......@@ -618,8 +659,8 @@ void MnaSolver<VarType>::logSystemMatrices() {
for (auto sys : mSwitchedMatrices) {
mSLog->info("Switching System matrix {:s} \n{:s}",
sys.first.to_string(), Logger::matrixToString(sys.second));
mSLog->info("LU Factorization for System Matrix {:s} \n{:s}",
sys.first.to_string(), Logger::matrixToString(mLuFactorizations[sys.first].matrixLU()));
//mSLog->info("LU Factorization for System Matrix {:s} \n{:s}",
// sys.first.to_string(), Logger::matrixToString(mLuFactorizations[sys.first].matrixLU()));
}
}
mSLog->info("Right side vector: \n{}", mRightSideVector);
......
......@@ -96,8 +96,12 @@ void MnaSolverGpu<VarType>::allocateDeviceMemory() {
template <typename VarType>
void MnaSolverGpu<VarType>::copySystemMatrixToDevice() {
auto *mat = &MnaSolver<VarType>::systemMatrix()(0);
CUDA_ERROR_HANDLER(cudaMemcpy(mDeviceCopy.matrix, mat, mDeviceCopy.size * mDeviceCopy.size * sizeof(Real), cudaMemcpyHostToDevice))
#ifdef WITH_SPARSE
auto *data = Matrix(MnaSolver<VarType>::systemMatrix()).data();
#else
auto *data = MnaSolver<VarType>::systemMatrix().data();
#endif
CUDA_ERROR_HANDLER(cudaMemcpy(mDeviceCopy.matrix, data, mDeviceCopy.size * mDeviceCopy.size * sizeof(Real), cudaMemcpyHostToDevice))
}
template <typename VarType>
......
......@@ -34,8 +34,13 @@ void MnaSolverSysRecomp<VarType>::initializeSystem() {
// Do not stamp variable elements yet
auto varcomp = std::dynamic_pointer_cast<MNAVariableCompInterface>(comp);
if (varcomp) continue;
#ifdef WITH_SPARSE
Matrix mat = Matrix(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
comp->mnaApplySystemMatrixStamp(mat);
this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)] = mat.sparseView();
#else
comp->mnaApplySystemMatrixStamp(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#endif
auto idObj = std::dynamic_pointer_cast<IdentifiedObject>(comp);
}
......@@ -46,10 +51,20 @@ void MnaSolverSysRecomp<VarType>::initializeSystem() {
// Now stamp variable elements
this->mSLog->info("Stamping variable elements");
for (auto varElem : this->mMNAIntfVariableComps) {
#ifdef WITH_SPARSE
Matrix mat = Matrix(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
varElem->mnaApplySystemMatrixStamp(mat);
this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)] = mat.sparseView();
#else
varElem->mnaApplySystemMatrixStamp(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#endif
}
#ifdef WITH_SPARSE
this->mLuFactorizations[std::bitset<SWITCH_NUM>(0)].analyzePattern(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
this->mLuFactorizations[std::bitset<SWITCH_NUM>(0)].factorize(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#else
this->mLuFactorizations[std::bitset<SWITCH_NUM>(0)] = Eigen::PartialPivLU<Matrix>(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#endif
// Initialize source vector for debugging
for (auto comp : this->mMNAComponents) {
comp->mnaApplyRightSideVectorStamp(this->mRightSideVector);
......@@ -77,13 +92,23 @@ void MnaSolverSysRecomp<VarType>::updateSystemMatrix(Real time) {
// Create system matrix with changed variable elements
for (auto comp : this->mMNAIntfVariableComps) {
#ifdef WITH_SPARSE
Matrix mat = Matrix(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
comp->mnaApplySystemMatrixStamp(mat);
this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)] = mat.sparseView();
#else
comp->mnaApplySystemMatrixStamp(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#endif
auto idObj = std::dynamic_pointer_cast<IdentifiedObject>(comp);
this->mSLog->debug("Updating {:s} {:s} in system matrix (variabel component)",
idObj->type(), idObj->name());
}
#ifdef WITH_SPARSE
this->mLuFactorizations[std::bitset<SWITCH_NUM>(0)].analyzePattern(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
this->mLuFactorizations[std::bitset<SWITCH_NUM>(0)].factorize(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#else
this->mLuFactorizations[std::bitset<SWITCH_NUM>(0)] = Eigen::PartialPivLU<Matrix>(this->mSwitchedMatrices[std::bitset<SWITCH_NUM>(0)]);
#endif
mUpdateSysMatrix = false;
}
......
......@@ -63,6 +63,8 @@ namespace CPS {
typedef Eigen::Matrix<Real, Eigen::Dynamic, 1> Vector;
/// @brief Sparse matrix for real numbers.
typedef Eigen::SparseMatrix<Real, Eigen::ColMajor> SparseMatrix;
/// @brief Sparse matrix for real numbers (row major).
typedef Eigen::SparseMatrix<Real, Eigen::RowMajor> SparseMatrixRow;
/// @brief Sparse matrix for complex numbers.
typedef Eigen::SparseMatrix<Complex, Eigen::ColMajor> SparseMatrixComp;
/// @brief Sparse matrix for complex numbers (row major).
......@@ -76,6 +78,8 @@ namespace CPS {
///
typedef Eigen::PartialPivLU<Matrix> LUFactorized;
///
typedef Eigen::SparseLU<SparseMatrix> LUFactorizedSparse;
///
typedef Eigen::Matrix<Real, Eigen::Dynamic, 1> Vector;
///
template<typename VarType>
......
......@@ -66,6 +66,14 @@ namespace CPS {
return ss.str();
}
static String sparseMatrixToString(const SparseMatrix& mat) {
return matrixToString(Matrix(mat));
}
static String sparseMatrixCompToString(const SparseMatrixComp& mat) {
return matrixCompToString(MatrixComp(mat));
}
static String phasorMatrixToString(const MatrixComp& mat) {
std::stringstream ss;
ss << std::scientific << Math::abs(mat) << "\n\n" << Math::phase(mat);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment