From c01d8e8607e002cffccd203d88a315770303411a Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 10 Jun 2024 15:50:35 +0200 Subject: [PATCH] refactor: Use Row-major matrix internally --- include/tug/Core/BTCS.hpp | 17 +++++++------- include/tug/Core/FTCS.hpp | 6 ++--- include/tug/Core/Matrix.hpp | 21 +++++++++++++++++ include/tug/Grid.hpp | 47 ++++++++++++++----------------------- include/tug/Simulation.hpp | 6 ++--- 5 files changed, 51 insertions(+), 46 deletions(-) create mode 100644 include/tug/Core/Matrix.hpp diff --git a/include/tug/Core/BTCS.hpp b/include/tug/Core/BTCS.hpp index 9b952ab..318ddf5 100644 --- a/include/tug/Core/BTCS.hpp +++ b/include/tug/Core/BTCS.hpp @@ -10,6 +10,7 @@ #ifndef BTCS_H_ #define BTCS_H_ +#include "Matrix.hpp" #include "TugUtils.hpp" #include @@ -51,7 +52,7 @@ constexpr std::pair calcBoundaryCoeffClosed(T alpha_center, T alpha_side, // creates coefficient matrix for next time step from alphas in x-direction template static Eigen::SparseMatrix -createCoeffMatrix(const typename Grid::RowMajMat &alpha, +createCoeffMatrix(const RowMajMat &alpha, const std::vector> &bcLeft, const std::vector> &bcRight, const std::vector> &inner_bc, int numCols, @@ -160,9 +161,8 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc, // concentrations template static Eigen::VectorX -createSolutionVector(const typename Grid::RowMajMat &concentrations, - const typename Grid::RowMajMat &alphaX, - const typename Grid::RowMajMat &alphaY, +createSolutionVector(const RowMajMat &concentrations, + const RowMajMat &alphaX, const RowMajMat &alphaY, const std::vector> &bcLeft, const std::vector> &bcRight, const std::vector> &bcTop, @@ -405,22 +405,21 @@ static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); - Eigen::MatrixX concentrations_t1 = - Eigen::MatrixX::Constant(rowMax, colMax, 0); + RowMajMat concentrations_t1 = RowMajMat::Constant(rowMax, colMax, 0); Eigen::VectorX row_t1(colMax); Eigen::SparseMatrix A; Eigen::VectorX b; - auto alphaX = grid.getAlphaX(); - auto alphaY = grid.getAlphaY(); + RowMajMat alphaX = grid.getAlphaX(); + RowMajMat alphaY = grid.getAlphaY(); const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); const auto &bcTop = bc.getBoundarySide(BC_SIDE_TOP); const auto &bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); - Eigen::MatrixX concentrations = grid.getConcentrations(); + RowMajMat &concentrations = grid.getConcentrations(); #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) for (int i = 0; i < rowMax; i++) { diff --git a/include/tug/Core/FTCS.hpp b/include/tug/Core/FTCS.hpp index 76d8d1b..a3012b3 100644 --- a/include/tug/Core/FTCS.hpp +++ b/include/tug/Core/FTCS.hpp @@ -228,8 +228,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixX concentrations_t1 = - Eigen::MatrixX::Constant(1, colMax, 0); + RowMajMat concentrations_t1 = RowMajMat::Constant(1, colMax, 0); // only one row in 1D case -> row constant at index 0 int row = 0; @@ -270,8 +269,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixX concentrations_t1 = - Eigen::MatrixX::Constant(rowMax, colMax, 0); + RowMajMat concentrations_t1 = RowMajMat::Constant(rowMax, colMax, 0); // inner cells // these are independent of the boundary condition type diff --git a/include/tug/Core/Matrix.hpp b/include/tug/Core/Matrix.hpp new file mode 100644 index 0000000..9b1686e --- /dev/null +++ b/include/tug/Core/Matrix.hpp @@ -0,0 +1,21 @@ +#pragma once + +#include + +namespace tug { +/** + * @brief Alias template for a row-major matrix using Eigen library. + * + * This alias template defines a type `RowMajMat` which represents a row-major + * matrix using the Eigen library. It is a template that takes a type `T` as its + * template parameter. The matrix is dynamically sized with `Eigen::Dynamic` for + * both rows and columns. The matrix is stored in row-major order. + * + * @tparam T The type of the matrix elements. + */ +template +using RowMajMat = + Eigen::Matrix; + +template using RowMajMatMap = Eigen::Map>; +} // namespace tug \ No newline at end of file diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 3dbd2bd..b8e471a 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -8,6 +8,7 @@ * */ +#include "Core/Matrix.hpp" #include #include #include @@ -24,9 +25,6 @@ namespace tug { */ template class Grid { public: - using RowMajMat = - Eigen::Matrix; - /** * @brief Constructs a new 1D-Grid object of a given length, which holds a * matrix with concentrations and a respective matrix of alpha coefficients. @@ -45,8 +43,8 @@ public: this->deltaCol = static_cast(this->domainCol) / static_cast(this->col); // -> 1 - this->concentrations = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); + this->concentrations = RowMajMat::Constant(1, col, MAT_INIT_VAL); + this->alphaX = RowMajMat::Constant(1, col, MAT_INIT_VAL); } /** @@ -73,9 +71,9 @@ public: this->deltaCol = static_cast(this->domainCol) / static_cast(this->col); // -> 1 - this->concentrations = RowMajMat::Constant(row, col, MAT_INIT_VAL); - this->alphaX = RowMajMat::Constant(row, col, MAT_INIT_VAL); - this->alphaY = RowMajMat::Constant(row, col, MAT_INIT_VAL); + this->concentrations = RowMajMat::Constant(row, col, MAT_INIT_VAL); + this->alphaX = RowMajMat::Constant(row, col, MAT_INIT_VAL); + this->alphaY = RowMajMat::Constant(row, col, MAT_INIT_VAL); } /** @@ -85,7 +83,7 @@ public: * Matrix must have correct dimensions as defined in row and col. (Or length, * in 1D case). */ - void setConcentrations(const Eigen::MatrixX &concentrations) { + void setConcentrations(const RowMajMat &concentrations) { if (concentrations.rows() != this->row || concentrations.cols() != this->col) { throw std::invalid_argument( @@ -103,9 +101,7 @@ public: * in 1D case). There is no check for correct dimensions, so be careful! */ void setConcentrations(T *concentrations) { - Eigen::Map< - Eigen::Matrix> - map(concentrations, this->row, this->col); + tug::RowMajMatMap map(concentrations, this->row, this->col); this->concentrations = map; } @@ -115,7 +111,7 @@ public: * @return An Eigen3 matrix holding the concentrations and having * the same dimensions as the grid. */ - const auto &getConcentrations() { return this->concentrations; } + auto &getConcentrations() { return this->concentrations; } /** * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one @@ -124,7 +120,7 @@ public: * @param alpha An Eigen3 MatrixX with 1 row holding the alpha * coefficients. Matrix columns must have same size as length of grid. */ - void setAlpha(const Eigen::MatrixX &alpha) { + void setAlpha(const RowMajMat &alpha) { if (dim != 1) { throw std::invalid_argument( "Grid is not one dimensional, you should probably " @@ -152,9 +148,7 @@ public: "Grid is not one dimensional, you should probably " "use 2D setter function!"); } - Eigen::Map< - Eigen::Matrix> - map(alpha, 1, this->col); + RowMajMatMap map(alpha, 1, this->col); this->alphaX = map; } @@ -167,8 +161,7 @@ public: * @param alphaY An Eigen3 MatrixX holding the alpha coefficients in * y-direction. Matrix must be of same size as the grid. */ - void setAlpha(const Eigen::MatrixX &alphaX, - const Eigen::MatrixX &alphaY) { + void setAlpha(const RowMajMat &alphaX, const RowMajMat &alphaY) { if (dim != 2) { throw std::invalid_argument( "Grid is not two dimensional, you should probably " @@ -206,12 +199,8 @@ public: "Grid is not two dimensional, you should probably " "use 1D setter function!"); } - Eigen::Map< - Eigen::Matrix> - mapX(alphaX, this->row, this->col); - Eigen::Map< - Eigen::Matrix> - mapY(alphaY, this->row, this->col); + RowMajMatMap mapX(alphaX, this->row, this->col); + RowMajMatMap mapY(alphaY, this->row, this->col); this->alphaX = mapX; this->alphaY = mapY; } @@ -390,11 +379,9 @@ private: T deltaCol; // delta in x-direction (between columns) T deltaRow{0}; // delta in y-direction (between rows) - RowMajMat concentrations; // Matrix holding grid concentrations - RowMajMat alphaX; // Matrix holding alpha coefficients in x-direction - RowMajMat alphaY; // Matrix holding alpha coefficients in y-direction - - using RowMajMatMap = Eigen::Map; + RowMajMat concentrations; // Matrix holding grid concentrations + RowMajMat alphaX; // Matrix holding alpha coefficients in x-direction + RowMajMat alphaY; // Matrix holding alpha coefficients in y-direction static constexpr T MAT_INIT_VAL = 0; }; diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 469b411..0e14197 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -446,9 +446,9 @@ public: // TODO this implementation is very inefficient! // a separate implementation that sets up a specific tridiagonal matrix // for Crank-Nicolson would be better - Eigen::MatrixX concentrations; - Eigen::MatrixX concentrationsFTCS; - Eigen::MatrixX concentrationsResult; + RowMajMat concentrations; + RowMajMat concentrationsFTCS; + RowMajMat concentrationsResult; for (int i = 0; i < iterations * innerIterations; i++) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole();