Add a separate EigenvalueUtils class

Signed-off-by: Georgii Tishenin <[email protected]>

dpsim-models/include/dpsim-models/MathUtils.h

@@ -63,10 +63,6 @@ class Math {
   // | Interharmonics harm1-harm2                  | Re(row,col)_harm2 | Re(row,col)_harm2 |
   // | Interharmonics harm1-harm2                  | Im(row,col)_harm2 | Im(row,col)_harm2 |
 
-  static MatrixComp returnNonZeroElements(const MatrixComp &mat);
-
-  static MatrixComp convertRealEquivalentToComplexMatrix(const Matrix &realEquivalentMatrix);
-
   static void setMatrixElement(SparseMatrixRow &mat, Matrix::Index row,
                                Matrix::Index column, Complex value,
                                Int maxFreq = 1, Int freqIdx = 0);

dpsim-models/src/MathUtils.cpp

@@ -97,47 +97,6 @@ Real Math::realFromVectorElement(const Matrix &mat, Matrix::Index row) {
   return mat(row, 0);
 }
 
-// TODO: consider resizing existing matrix instead of creating a new one (performance, memory)
-MatrixComp Math::returnNonZeroElements(const MatrixComp &mat) {
-  Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> mask =
-      (mat.array().abs() > 1e-14);
-  int nonZeroCount = mask.cast<int>().sum();
-
-  Eigen::MatrixXcd nonZeroMatrix(nonZeroCount, 1);
-  int index = 0;
-  for (int i = 0; i < mask.rows(); ++i) {
-    for (int j = 0; j < mask.cols(); ++j) {
-      if (mask(i, j)) {
-        nonZeroMatrix(index++) = mat(i, j);
-      }
-    }
-  }
-  return nonZeroMatrix;
-}
-
-MatrixComp Math::convertRealEquivalentToComplexMatrix(const Matrix &realEquivalentMatrix) {
-  // The size of the complex matrix is half the size of the real matrix
-  int size = realEquivalentMatrix.rows() / 2;
-
-  // Create a complex matrix of the appropriate size
-  MatrixComp complexMatrix(size, size);
-
-  // Iterate over the complex matrix
-  for (int i = 0; i < size; ++i) {
-    for (int j = 0; j < size; ++j) {
-      // The real part is in the upper left quadrant of the real matrix
-      double realPart = realEquivalentMatrix(i, j);
-
-      // The imaginary part is in the lower left quadrant of the real matrix
-      double imagPart = realEquivalentMatrix(i + size, j);
-
-      // Assign the complex number to the complex matrix
-      complexMatrix(i, j) = std::complex<double>(realPart, imagPart);
-    }
-  }
-  return complexMatrix;
-}
-
 void Math::setMatrixElement(SparseMatrixRow &mat, Matrix::Index row,
                             Matrix::Index column, Complex value, Int maxFreq,
                             Int freqIdx) {

dpsim/include/dpsim/EigenvalueUtils.h

@@ -0,0 +1,14 @@
+#pragma once
+
+#include <dpsim/Definitions.h>
+
+namespace DPsim {
+class EigenvalueUtils {
+public:
+  static MatrixComp returnNonZeroElements(const MatrixComp &mat);
+
+  static MatrixComp
+  convertRealEquivalentToComplexMatrix(const Matrix &realEquivalentMatrix);
+};
+
+} // namespace DPsim

dpsim/include/dpsim/MNAEigenvalueExtractor.h

@@ -2,6 +2,7 @@
 
 #include <dpsim-models/Solver/EigenvalueDynamicCompInterface.h>
 #include <dpsim-models/SystemTopology.h>
+#include <dpsim/EigenvalueUtils.h>
 #include <dpsim/MNAEigenvalueLogger.h>
 
 namespace DPsim {

dpsim/src/CMakeLists.txt

@@ -5,6 +5,7 @@ set(DPSIM_SOURCES
 	MNASolverDirect.cpp
 	MNAEigenvalueExtractor.cpp
 	MNAEigenvalueLogger.cpp
+	EigenvalueUtils.cpp
 	DenseLUAdapter.cpp
 	SparseLUAdapter.cpp
 	DirectLinearSolverConfiguration.cpp

dpsim/src/EigenvalueUtils.cpp

@@ -0,0 +1,45 @@
+#include <dpsim/EigenvalueUtils.h>
+
+using namespace DPsim;
+
+// TODO: consider resizing existing matrix instead of creating a new one (performance, memory)
+MatrixComp EigenvalueUtils::returnNonZeroElements(const MatrixComp &mat) {
+  Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> mask =
+      (mat.array().abs() > 1e-14);
+  int nonZeroCount = mask.cast<int>().sum();
+
+  Eigen::MatrixXcd nonZeroMatrix(nonZeroCount, 1);
+  int index = 0;
+  for (int i = 0; i < mask.rows(); ++i) {
+    for (int j = 0; j < mask.cols(); ++j) {
+      if (mask(i, j)) {
+        nonZeroMatrix(index++) = mat(i, j);
+      }
+    }
+  }
+  return nonZeroMatrix;
+}
+
+MatrixComp EigenvalueUtils::convertRealEquivalentToComplexMatrix(
+    const Matrix &realEquivalentMatrix) {
+  // The size of the complex matrix is half the size of the real matrix
+  int size = realEquivalentMatrix.rows() / 2;
+
+  // Create a complex matrix of the appropriate size
+  MatrixComp complexMatrix(size, size);
+
+  // Iterate over the complex matrix
+  for (int i = 0; i < size; ++i) {
+    for (int j = 0; j < size; ++j) {
+      // The real part is in the upper left quadrant of the real matrix
+      double realPart = realEquivalentMatrix(i, j);
+
+      // The imaginary part is in the lower left quadrant of the real matrix
+      double imagPart = realEquivalentMatrix(i + size, j);
+
+      // Assign the complex number to the complex matrix
+      complexMatrix(i, j) = std::complex<double>(realPart, imagPart);
+    }
+  }
+  return complexMatrix;
+}

dpsim/src/MNAEigenvalueExtractor.cpp

@@ -110,7 +110,7 @@ void MNAEigenvalueExtractor<Complex>::calculateStateMatrix(
     const Matrix &powerSystemMatrix) {
   // TODO: use right hand side solving of factorized power system matrix instead of inversion (performance).
   MatrixComp compPowerSystemMatrix =
-      CPS::Math::convertRealEquivalentToComplexMatrix(powerSystemMatrix);
+      EigenvalueUtils::convertRealEquivalentToComplexMatrix(powerSystemMatrix);
   MatrixComp intermediateResult =
       compPowerSystemMatrix.inverse() * mNodeBranchIncidenceMatrix;
   mStateMatrix = mSignMatrix - mDiscretizationMatrix *
@@ -122,7 +122,7 @@ template <typename VarType>
 void MNAEigenvalueExtractor<VarType>::computeDiscreteEigenvalues() {
   auto discreteEigenvaluesIncludingZeros = mStateMatrix.eigenvalues();
   **mDiscreteEigenvalues =
-      CPS::Math::returnNonZeroElements(discreteEigenvaluesIncludingZeros);
+      EigenvalueUtils::returnNonZeroElements(discreteEigenvaluesIncludingZeros);
   // TODO: filter out eigenvalues = -1 + 0i to avoid division by zero in recoverEigenvalues().
 }