From e64e8dfd5e43983ad377f71640b6b68c7d6aaf9d Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 3 Jun 2024 23:48:54 +0200 Subject: [PATCH 1/3] feat: Add support for setting concentrations from a pointer refactor: Use Row-major matrix internally --- include/tug/Core/BTCS.hpp | 8 +-- include/tug/Grid.hpp | 103 +++++++++++++++++++++++++++++++------- test/testGrid.cpp | 15 ++++++ 3 files changed, 105 insertions(+), 21 deletions(-) diff --git a/include/tug/Core/BTCS.hpp b/include/tug/Core/BTCS.hpp index 8b8f7a5..9b952ab 100644 --- a/include/tug/Core/BTCS.hpp +++ b/include/tug/Core/BTCS.hpp @@ -51,7 +51,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 Eigen::MatrixX &alpha, +createCoeffMatrix(const typename Grid::RowMajMat &alpha, const std::vector> &bcLeft, const std::vector> &bcRight, const std::vector> &inner_bc, int numCols, @@ -160,9 +160,9 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc, // concentrations template static Eigen::VectorX -createSolutionVector(const Eigen::MatrixX &concentrations, - const Eigen::MatrixX &alphaX, - const Eigen::MatrixX &alphaY, +createSolutionVector(const typename Grid::RowMajMat &concentrations, + const typename Grid::RowMajMat &alphaX, + const typename Grid::RowMajMat &alphaY, const std::vector> &bcLeft, const std::vector> &bcRight, const std::vector> &bcTop, diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index c4f35ec..3dbd2bd 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -10,6 +10,8 @@ #include #include +#include +#include #include namespace tug { @@ -22,6 +24,9 @@ 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. @@ -68,9 +73,9 @@ public: this->deltaCol = static_cast(this->domainCol) / static_cast(this->col); // -> 1 - this->concentrations = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); - this->alphaY = Eigen::MatrixX::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); } /** @@ -90,13 +95,27 @@ public: this->concentrations = concentrations; } + /** + * @brief Sets the concentrations matrix for a 1D or 2D-Grid. + * + * @param concentrations A pointer to an array holding the concentrations. + * Array must have correct dimensions as defined in row and col. (Or length, + * 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); + this->concentrations = map; + } + /** * @brief Gets the concentrations matrix for a Grid. * * @return An Eigen3 matrix holding the concentrations and having * the same dimensions as the grid. */ - const Eigen::MatrixX &getConcentrations() { return this->concentrations; } + const auto &getConcentrations() { return this->concentrations; } /** * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one @@ -119,6 +138,26 @@ public: this->alphaX = alpha; } + /** + * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one + * dimensional. + * + * @param alpha A pointer to an array holding the alpha coefficients. Array + * must have correct dimensions as defined in length. There is no check for + * correct dimensions, so be careful! + */ + void setAlpha(T *alpha) { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use 2D setter function!"); + } + Eigen::Map< + Eigen::Matrix> + map(alpha, 1, this->col); + this->alphaX = map; + } + /** * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two * dimensional. @@ -150,13 +189,40 @@ public: this->alphaY = alphaY; } + /** + * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two + * dimensional. + * + * @param alphaX A pointer to an array holding the alpha coefficients in + * x-direction. Array must have correct dimensions as defined in row and col. + * There is no check for correct dimensions, so be careful! + * @param alphaY A pointer to an array holding the alpha coefficients in + * y-direction. Array must have correct dimensions as defined in row and col. + * There is no check for correct dimensions, so be careful! + */ + void setAlpha(T *alphaX, T *alphaY) { + if (dim != 2) { + throw std::invalid_argument( + "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); + this->alphaX = mapX; + this->alphaY = mapY; + } + /** * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * * @return A matrix with 1 row holding the alpha coefficients. */ - const Eigen::MatrixX &getAlpha() const { + const auto &getAlpha() const { if (dim != 1) { throw std::invalid_argument( "Grid is not one dimensional, you should probably " @@ -172,7 +238,7 @@ public: * * @return A matrix holding the alpha coefficients in x-direction. */ - const Eigen::MatrixX &getAlphaX() const { + const auto &getAlphaX() const { if (dim != 2) { throw std::invalid_argument( @@ -188,7 +254,7 @@ public: * * @return A matrix holding the alpha coefficients in y-direction. */ - const Eigen::MatrixX &getAlphaY() const { + const auto &getAlphaY() const { if (dim != 2) { throw std::invalid_argument( @@ -316,16 +382,19 @@ public: } private: - int col; // number of grid columns - int row{1}; // number of grid rows - int dim; // 1D or 2D - T domainCol; // number of domain columns - T domainRow{0}; // number of domain rows - T deltaCol; // delta in x-direction (between columns) - T deltaRow{0}; // delta in y-direction (between rows) - Eigen::MatrixX concentrations; // Matrix holding grid concentrations - Eigen::MatrixX alphaX; // Matrix holding alpha coefficients in x-direction - Eigen::MatrixX alphaY; // Matrix holding alpha coefficients in y-direction + int col; // number of grid columns + int row{1}; // number of grid rows + int dim; // 1D or 2D + T domainCol; // number of domain columns + T domainRow{0}; // number of domain rows + 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; static constexpr T MAT_INIT_VAL = 0; }; diff --git a/test/testGrid.cpp b/test/testGrid.cpp index 8e17526..9eee265 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -1,6 +1,7 @@ #include #include #include +#include using namespace Eigen; using namespace std; @@ -250,4 +251,18 @@ TEST_CASE("2D Grid64 non-quadratic") { dr = -2; CHECK_THROWS(grid.setDomain(dr, dc)); } + + SUBCASE("set concentration from pointer") { + std::vector concentrations(r * c); + + for (int i = 0; i < r * c; i++) { + concentrations[i] = i; + } + + grid.setConcentrations(concentrations.data()); + + CHECK_EQ(grid.getConcentrations()(0, 0), 0); + CHECK_EQ(grid.getConcentrations()(0, 1), 1); + CHECK_EQ(grid.getConcentrations()(1, 0), c); + } } From c01d8e8607e002cffccd203d88a315770303411a Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 10 Jun 2024 15:50:35 +0200 Subject: [PATCH 2/3] 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(); From 74b46f111b069c370d46757ab276e831a689e7f2 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 10 Jun 2024 16:02:37 +0200 Subject: [PATCH 3/3] perf: Minimize copy operations --- include/tug/Core/BTCS.hpp | 32 +++++++++++--------------------- 1 file changed, 11 insertions(+), 21 deletions(-) diff --git a/include/tug/Core/BTCS.hpp b/include/tug/Core/BTCS.hpp index 318ddf5..a1081af 100644 --- a/include/tug/Core/BTCS.hpp +++ b/include/tug/Core/BTCS.hpp @@ -369,7 +369,7 @@ static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, const auto inner_bc = bc.getInnerBoundaryRow(0); - Eigen::MatrixX concentrations = grid.getConcentrations(); + RowMajMat &concentrations = grid.getConcentrations(); int rowIndex = 0; A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex, sx); // this is exactly same as in 2D @@ -385,13 +385,13 @@ static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue(); } - concentrations_t1 = solverFunc(A, b); + concentrations = solverFunc(A, b); - for (int j = 0; j < length; j++) { - concentrations(0, j) = concentrations_t1(j); - } + // for (int j = 0; j < length; j++) { + // concentrations(0, j) = concentrations_t1(j); + // } - grid.setConcentrations(concentrations); + // grid.setConcentrations(concentrations); } // BTCS solution for 2D grid @@ -405,8 +405,7 @@ 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()); - RowMajMat concentrations_t1 = RowMajMat::Constant(rowMax, colMax, 0); - Eigen::VectorX row_t1(colMax); + RowMajMat concentrations_t1(rowMax, colMax); Eigen::SparseMatrix A; Eigen::VectorX b; @@ -421,7 +420,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, RowMajMat &concentrations = grid.getConcentrations(); -#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) +#pragma omp parallel for num_threads(numThreads) private(A, b) for (int i = 0; i < rowMax; i++) { auto inner_bc = bc.getInnerBoundaryRow(i); @@ -429,17 +428,14 @@ static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, inner_bc, colMax, i, sx, sy); - row_t1 = solverFunc(A, b); - - concentrations_t1.row(i) = row_t1; + concentrations_t1.row(i) = solverFunc(A, b); } concentrations_t1.transposeInPlace(); - concentrations.transposeInPlace(); alphaX.transposeInPlace(); alphaY.transposeInPlace(); -#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) +#pragma omp parallel for num_threads(numThreads) private(A, b) for (int i = 0; i < colMax; i++) { auto inner_bc = bc.getInnerBoundaryCol(i); // swap alphas, boundary conditions and sx/sy for column-wise calculation @@ -447,14 +443,8 @@ static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, bcLeft, bcRight, inner_bc, rowMax, i, sy, sx); - row_t1 = solverFunc(A, b); - - concentrations.row(i) = row_t1; + concentrations.col(i) = solverFunc(A, b); } - - concentrations.transposeInPlace(); - - grid.setConcentrations(concentrations); } // entry point for EigenLU solver; differentiate between 1D and 2D grid