max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-15 10:28:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

refactor: Use Row-major matrix internally

Browse Source
This commit is contained in:
Max Luebke 2024-06-10 15:50:35 +02:00
parent e64e8dfd5e
commit c01d8e8607
5 changed files with 51 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
include/tug/Core/BTCS.hpp
View File

@ -10,6 +10,7 @@
#ifndef BTCS_H_
#define BTCS_H_
#include "Matrix.hpp"
#include "TugUtils.hpp"
#include <cstddef>
@ -51,7 +52,7 @@ constexpr std::pair<T, T> calcBoundaryCoeffClosed(T alpha_center, T alpha_side,
// creates coefficient matrix for next time step from alphas in x-direction
template <class T>
static Eigen::SparseMatrix<T>
createCoeffMatrix(const typename Grid<T>::RowMajMat &alpha,
createCoeffMatrix(const RowMajMat<T> &alpha,
const std::vector<BoundaryElement<T>> &bcLeft,
const std::vector<BoundaryElement<T>> &bcRight,
const std::vector<std::pair<bool, T>> &inner_bc, int numCols,
@ -160,9 +161,8 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc,
// concentrations
template <class T>
static Eigen::VectorX<T>
createSolutionVector(const typename Grid<T>::RowMajMat &concentrations,
const typename Grid<T>::RowMajMat &alphaX,
const typename Grid<T>::RowMajMat &alphaY,
createSolutionVector(const RowMajMat<T> &concentrations,
const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY,
const std::vector<BoundaryElement<T>> &bcLeft,
const std::vector<BoundaryElement<T>> &bcRight,
const std::vector<BoundaryElement<T>> &bcTop,
@ -405,22 +405,21 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
Eigen::MatrixX<T> concentrations_t1 =
Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
Eigen::VectorX<T> row_t1(colMax);
Eigen::SparseMatrix<T> A;
Eigen::VectorX<T> b;
auto alphaX = grid.getAlphaX();
auto alphaY = grid.getAlphaY();
RowMajMat<T> alphaX = grid.getAlphaX();
RowMajMat<T> 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<T> concentrations = grid.getConcentrations();
RowMajMat<T> &concentrations = grid.getConcentrations();
#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
for (int i = 0; i < rowMax; i++) {

6
include/tug/Core/FTCS.hpp
View File

@ -228,8 +228,7 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
T deltaCol = grid.getDeltaCol();
// matrix for concentrations at time t+1
Eigen::MatrixX<T> concentrations_t1 =
Eigen::MatrixX<T>::Constant(1, colMax, 0);
RowMajMat<T> concentrations_t1 = RowMajMat<T>::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<T> &grid, Boundary<T> &bc, T timestep,
T deltaCol = grid.getDeltaCol();
// matrix for concentrations at time t+1
Eigen::MatrixX<T> concentrations_t1 =
Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
// inner cells
// these are independent of the boundary condition type

21
include/tug/Core/Matrix.hpp Normal file
View File

@ -0,0 +1,21 @@
#pragma once
#include <Eigen/Core>
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 <typename T>
using RowMajMat =
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
template <typename T> using RowMajMatMap = Eigen::Map<RowMajMat<T>>;
} // namespace tug

47
include/tug/Grid.hpp
View File

@ -8,6 +8,7 @@
*
*/
#include "Core/Matrix.hpp"
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <Eigen/src/Core/Matrix.h>
@ -24,9 +25,6 @@ namespace tug {
*/
template <class T> class Grid {
public:
using RowMajMat =
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
/**
* @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<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
this->concentrations = Eigen::MatrixX<T>::Constant(1, col, MAT_INIT_VAL);
this->alphaX = Eigen::MatrixX<T>::Constant(1, col, MAT_INIT_VAL);
this->concentrations = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
this->alphaX = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
}
/**
@ -73,9 +71,9 @@ public:
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(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<T>::Constant(row, col, MAT_INIT_VAL);
this->alphaX = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
this->alphaY = RowMajMat<T>::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<T> &concentrations) {
void setConcentrations(const RowMajMat<T> &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<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
map(concentrations, this->row, this->col);
tug::RowMajMatMap<T> 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<T> with 1 row holding the alpha
* coefficients. Matrix columns must have same size as length of grid.
*/
void setAlpha(const Eigen::MatrixX<T> &alpha) {
void setAlpha(const RowMajMat<T> &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<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
map(alpha, 1, this->col);
RowMajMatMap<T> map(alpha, 1, this->col);
this->alphaX = map;
}
@ -167,8 +161,7 @@ public:
* @param alphaY An Eigen3 MatrixX<T> holding the alpha coefficients in
* y-direction. Matrix must be of same size as the grid.
*/
void setAlpha(const Eigen::MatrixX<T> &alphaX,
const Eigen::MatrixX<T> &alphaY) {
void setAlpha(const RowMajMat<T> &alphaX, const RowMajMat<T> &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<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
mapX(alphaX, this->row, this->col);
Eigen::Map<
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
mapY(alphaY, this->row, this->col);
RowMajMatMap<T> mapX(alphaX, this->row, this->col);
RowMajMatMap<T> 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>;
RowMajMat<T> concentrations; // Matrix holding grid concentrations
RowMajMat<T> alphaX; // Matrix holding alpha coefficients in x-direction
RowMajMat<T> alphaY; // Matrix holding alpha coefficients in y-direction
static constexpr T MAT_INIT_VAL = 0;
};

6
include/tug/Simulation.hpp
View File

@ -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<T> concentrations;
Eigen::MatrixX<T> concentrationsFTCS;
Eigen::MatrixX<T> concentrationsResult;
RowMajMat<T> concentrations;
RowMajMat<T> concentrationsFTCS;
RowMajMat<T> concentrationsResult;
for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();