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

refactor: remove all 'using namespaces' from library

Browse Source
This commit is contained in:
Max Lübke 2023-09-14 12:04:03 +02:00
parent 2483019b89
commit 2096ee5cc3
9 changed files with 164 additions and 169 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
include/tug/Boundary.hpp
View File

@ -10,9 +10,6 @@
#include "Grid.hpp" #include "Grid.hpp"
#include <cstddef> #include <cstddef>
using namespace std;
using namespace Eigen;
/** /**
* @brief Enum defining the two implemented boundary conditions. * @brief Enum defining the two implemented boundary conditions.
* *
@ -156,7 +153,7 @@ public:
* @return vector<BoundaryElement> Contains the boundary conditions as * @return vector<BoundaryElement> Contains the boundary conditions as
* BoundaryElement objects. * BoundaryElement objects.
*/ */
const vector<BoundaryElement> getBoundarySide(BC_SIDE side); const std::vector<BoundaryElement> getBoundarySide(BC_SIDE side);
/** /**
* @brief Get thes Boundary Side Values as a vector. Value is -1 in case some * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some
@ -165,7 +162,7 @@ public:
* @param side Boundary side for which the values are to be returned. * @param side Boundary side for which the values are to be returned.
* @return VectorXd Vector with values as doubles. * @return VectorXd Vector with values as doubles.
*/ */
VectorXd getBoundarySideValues(BC_SIDE side); Eigen::VectorXd getBoundarySideValues(BC_SIDE side);
/** /**
* @brief Returns the boundary condition of a specified element on a given * @brief Returns the boundary condition of a specified element on a given
@ -207,7 +204,7 @@ public:
private: private:
Grid grid; // Boundary is directly dependent on the dimensions of a predefined Grid grid; // Boundary is directly dependent on the dimensions of a predefined
vector<vector<BoundaryElement>> std::vector<std::vector<BoundaryElement>>
boundaries; // Vector with Boundary Element information boundaries; // Vector with Boundary Element information
}; };

36
include/tug/Grid.hpp
View File

@ -11,8 +11,6 @@
#include <Eigen/Core> #include <Eigen/Core>
#include <Eigen/Sparse> #include <Eigen/Sparse>
using namespace Eigen;
class Grid { class Grid {
public: public:
/** /**
@ -45,7 +43,7 @@ public:
* must have correct dimensions as defined in row and col. (Or length, in 1D * must have correct dimensions as defined in row and col. (Or length, in 1D
* case). * case).
*/ */
void setConcentrations(MatrixXd concentrations); void setConcentrations(Eigen::MatrixXd concentrations);
/** /**
* @brief Gets the concentrations matrix for a Grid. * @brief Gets the concentrations matrix for a Grid.
@ -53,7 +51,7 @@ public:
* @return MatrixXd An Eigen3 matrix holding the concentrations and having the * @return MatrixXd An Eigen3 matrix holding the concentrations and having the
* same dimensions as the grid. * same dimensions as the grid.
*/ */
const MatrixXd getConcentrations(); const Eigen::MatrixXd getConcentrations();
/** /**
* @brief Set the alpha coefficients of a 1D-Grid. Grid must be one * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
@ -62,7 +60,7 @@ public:
* @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients.
* Matrix columns must have same size as length of grid. * Matrix columns must have same size as length of grid.
*/ */
void setAlpha(MatrixXd alpha); void setAlpha(Eigen::MatrixXd alpha);
/** /**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
@ -73,7 +71,7 @@ public:
* @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in
* y-direction. Matrix must be of same size as the grid. * y-direction. Matrix must be of same size as the grid.
*/ */
void setAlpha(MatrixXd alphaX, MatrixXd alphaY); void setAlpha(Eigen::MatrixXd alphaX, Eigen::MatrixXd alphaY);
/** /**
* @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
@ -81,7 +79,7 @@ public:
* *
* @return MatrixXd A matrix with 1 row holding the alpha coefficients. * @return MatrixXd A matrix with 1 row holding the alpha coefficients.
*/ */
const MatrixXd getAlpha(); const Eigen::MatrixXd getAlpha();
/** /**
* @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
@ -89,7 +87,7 @@ public:
* *
* @return MatrixXd A matrix holding the alpha coefficients in x-direction. * @return MatrixXd A matrix holding the alpha coefficients in x-direction.
*/ */
const MatrixXd getAlphaX(); const Eigen::MatrixXd getAlphaX();
/** /**
* @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
@ -97,7 +95,7 @@ public:
* *
* @return MatrixXd A matrix holding the alpha coefficients in y-direction. * @return MatrixXd A matrix holding the alpha coefficients in y-direction.
*/ */
const MatrixXd getAlphaY(); const Eigen::MatrixXd getAlphaY();
/** /**
* @brief Gets the dimensions of the grid. * @brief Gets the dimensions of the grid.
@ -166,16 +164,16 @@ public:
double getDeltaRow(); double getDeltaRow();
private: private:
int col; // number of grid columns int col; // number of grid columns
int row; // number of grid rows int row; // number of grid rows
int dim; // 1D or 2D int dim; // 1D or 2D
double domainCol; // number of domain columns double domainCol; // number of domain columns
double domainRow; // number of domain rows double domainRow; // number of domain rows
double deltaCol; // delta in x-direction (between columns) double deltaCol; // delta in x-direction (between columns)
double deltaRow; // delta in y-direction (between rows) double deltaRow; // delta in y-direction (between rows)
MatrixXd concentrations; // Matrix holding grid concentrations Eigen::MatrixXd concentrations; // Matrix holding grid concentrations
MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
}; };
#endif // GRID_H_ #endif // GRID_H_

6
include/tug/Simulation.hpp
View File

@ -12,8 +12,6 @@
#include "Boundary.hpp" #include "Boundary.hpp"
#include "Grid.hpp" #include "Grid.hpp"
using namespace std;
/** /**
* @brief Enum defining the two implemented solution approaches. * @brief Enum defining the two implemented solution approaches.
* *
@ -193,7 +191,7 @@ public:
* *
* @return string Filename with configured simulation parameters. * @return string Filename with configured simulation parameters.
*/ */
string createCSVfile(); std::string createCSVfile();
/** /**
* @brief Writes the currently calculated concentration values of the grid * @brief Writes the currently calculated concentration values of the grid
@ -202,7 +200,7 @@ public:
* @param filename Name of the file to which the concentration values are * @param filename Name of the file to which the concentration values are
* to be written. * to be written.
*/ */
void printConcentrationsCSV(string filename); void printConcentrationsCSV(std::string filename);
/** /**
* @brief Method starts the simulation process with the previously set * @brief Method starts the simulation process with the previously set

132
src/BTCS.cpp
View File

@ -14,13 +14,9 @@
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
#include <tug/Grid.hpp> #include <tug/Grid.hpp>
#define NUM_THREADS_BTCS 10
using namespace Eigen;
// calculates coefficient for left boundary in constant case // calculates coefficient for left boundary in constant case
static tuple<double, double> static std::tuple<double, double>
calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) { calcLeftBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff; double centerCoeff;
double rightCoeff; double rightCoeff;
@ -33,8 +29,8 @@ calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) {
} }
// calculates coefficient for left boundary in closed case // calculates coefficient for left boundary in closed case
static tuple<double, double> static std::tuple<double, double>
calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) { calcLeftBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, double sx) {
double centerCoeff; double centerCoeff;
double rightCoeff; double rightCoeff;
@ -46,9 +42,9 @@ calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) {
} }
// calculates coefficient for right boundary in constant case // calculates coefficient for right boundary in constant case
static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, static std::tuple<double, double>
int rowIndex, int n, calcRightBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, int n,
double sx) { double sx) {
double leftCoeff; double leftCoeff;
double centerCoeff; double centerCoeff;
@ -62,8 +58,9 @@ static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha,
} }
// calculates coefficient for right boundary in closed case // calculates coefficient for right boundary in closed case
static tuple<double, double> static std::tuple<double, double>
calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) { calcRightBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, int n,
double sx) {
double leftCoeff; double leftCoeff;
double centerCoeff; double centerCoeff;
@ -76,15 +73,14 @@ calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) {
} }
// creates coefficient matrix for next time step from alphas in x-direction // creates coefficient matrix for next time step from alphas in x-direction
static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, static Eigen::SparseMatrix<double>
vector<BoundaryElement> &bcLeft, createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
vector<BoundaryElement> &bcRight, std::vector<BoundaryElement> &bcRight, int numCols,
int numCols, int rowIndex, int rowIndex, double sx) {
double sx) {
// square matrix of column^2 dimension for the coefficients // square matrix of column^2 dimension for the coefficients
SparseMatrix<double> cm(numCols, numCols); Eigen::SparseMatrix<double> cm(numCols, numCols);
cm.reserve(VectorXi::Constant(numCols, 3)); cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
// left column // left column
BC_TYPE type = bcLeft[rowIndex].getType(); BC_TYPE type = bcLeft[rowIndex].getType();
@ -140,8 +136,8 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha,
// calculates explicity concentration at top boundary in constant case // calculates explicity concentration at top boundary in constant case
static double calcExplicitConcentrationsTopBoundaryConstant( static double calcExplicitConcentrationsTopBoundaryConstant(
MatrixXd &concentrations, MatrixXd &alpha, vector<BoundaryElement> &bcTop, Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alpha,
int rowIndex, int i, double sy) { std::vector<BoundaryElement> &bcTop, int rowIndex, int i, double sy) {
double c; double c;
c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) * c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
@ -156,8 +152,10 @@ static double calcExplicitConcentrationsTopBoundaryConstant(
} }
// calculates explicit concentration at top boundary in closed case // calculates explicit concentration at top boundary in closed case
static double calcExplicitConcentrationsTopBoundaryClosed( static double
MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) { calcExplicitConcentrationsTopBoundaryClosed(Eigen::MatrixXd &concentrations,
Eigen::MatrixXd &alpha,
int rowIndex, int i, double sy) {
double c; double c;
c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) * c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) *
@ -171,8 +169,8 @@ static double calcExplicitConcentrationsTopBoundaryClosed(
// calculates explicit concentration at bottom boundary in constant case // calculates explicit concentration at bottom boundary in constant case
static double calcExplicitConcentrationsBottomBoundaryConstant( static double calcExplicitConcentrationsBottomBoundaryConstant(
MatrixXd &concentrations, MatrixXd &alpha, Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alpha,
vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) { std::vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) {
double c; double c;
c = sy * alpha(rowIndex, i) * bcBottom[i].getValue() + c = sy * alpha(rowIndex, i) * bcBottom[i].getValue() +
@ -187,8 +185,10 @@ static double calcExplicitConcentrationsBottomBoundaryConstant(
} }
// calculates explicit concentration at bottom boundary in closed case // calculates explicit concentration at bottom boundary in closed case
static double calcExplicitConcentrationsBottomBoundaryClosed( static double
MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) { calcExplicitConcentrationsBottomBoundaryClosed(Eigen::MatrixXd &concentrations,
Eigen::MatrixXd &alpha,
int rowIndex, int i, double sy) {
double c; double c;
c = (1 - c = (1 -
@ -202,13 +202,14 @@ static double calcExplicitConcentrationsBottomBoundaryClosed(
// creates a solution vector for next time step from the current state of // creates a solution vector for next time step from the current state of
// concentrations // concentrations
static VectorXd createSolutionVector( static Eigen::VectorXd createSolutionVector(
MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX,
vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, Eigen::MatrixXd &alphaY, std::vector<BoundaryElement> &bcLeft,
vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom, std::vector<BoundaryElement> &bcRight, std::vector<BoundaryElement> &bcTop,
int length, int rowIndex, double sx, double sy) { std::vector<BoundaryElement> &bcBottom, int length, int rowIndex, double sx,
double sy) {
VectorXd sv(length); Eigen::VectorXd sv(length);
int numRows = concentrations.rows(); int numRows = concentrations.rows();
BC_TYPE type; BC_TYPE type;
@ -283,9 +284,10 @@ static VectorXd createSolutionVector(
// solver for linear equation system; A corresponds to coefficient matrix, // solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector // b to the solution vector
// use of EigenLU solver // use of EigenLU solver
static VectorXd EigenLUAlgorithm(SparseMatrix<double> &A, VectorXd &b) { static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix<double> &A,
Eigen::VectorXd &b) {
SparseLU<SparseMatrix<double>> solver; Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
solver.analyzePattern(A); solver.analyzePattern(A);
solver.factorize(A); solver.factorize(A);
@ -295,7 +297,8 @@ static VectorXd EigenLUAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
// solver for linear equation system; A corresponds to coefficient matrix, // solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector // b to the solution vector
// implementation of Thomas Algorithm // implementation of Thomas Algorithm
static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) { static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix<double> &A,
Eigen::VectorXd &b) {
uint32_t n = b.size(); uint32_t n = b.size();
Eigen::VectorXd a_diag(n); Eigen::VectorXd a_diag(n);
@ -337,22 +340,23 @@ static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {
} }
// BTCS solution for 1D grid // BTCS solution for 1D grid
static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, static void
VectorXd (*solverFunc)(SparseMatrix<double> &A, BTCS_1D(Grid &grid, Boundary &bc, double timestep,
VectorXd &b)) { Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
Eigen::VectorXd &b)) {
int length = grid.getLength(); int length = grid.getLength();
double sx = timestep / (grid.getDelta() * grid.getDelta()); double sx = timestep / (grid.getDelta() * grid.getDelta());
VectorXd concentrations_t1(length); Eigen::VectorXd concentrations_t1(length);
SparseMatrix<double> A; Eigen::SparseMatrix<double> A;
VectorXd b(length); Eigen::VectorXd b(length);
MatrixXd alpha = grid.getAlpha(); Eigen::MatrixXd alpha = grid.getAlpha();
vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
MatrixXd concentrations = grid.getConcentrations(); Eigen::MatrixXd concentrations = grid.getConcentrations();
int rowIndex = 0; int rowIndex = 0;
A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex,
sx); // this is exactly same as in 2D sx); // this is exactly same as in 2D
@ -376,29 +380,31 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double timestep,
} }
// BTCS solution for 2D grid // BTCS solution for 2D grid
static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, static void
VectorXd (*solverFunc)(SparseMatrix<double> &A, BTCS_2D(Grid &grid, Boundary &bc, double timestep,
VectorXd &b), Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
int numThreads) { Eigen::VectorXd &b),
int numThreads) {
int rowMax = grid.getRow(); int rowMax = grid.getRow();
int colMax = grid.getCol(); int colMax = grid.getCol();
double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); Eigen::MatrixXd concentrations_t1 =
VectorXd row_t1(colMax); Eigen::MatrixXd::Constant(rowMax, colMax, 0);
Eigen::VectorXd row_t1(colMax);
SparseMatrix<double> A; Eigen::SparseMatrix<double> A;
VectorXd b; Eigen::VectorXd b;
MatrixXd alphaX = grid.getAlphaX(); Eigen::MatrixXd alphaX = grid.getAlphaX();
MatrixXd alphaY = grid.getAlphaY(); Eigen::MatrixXd alphaY = grid.getAlphaY();
vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP); std::vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); std::vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
MatrixXd concentrations = grid.getConcentrations(); Eigen::MatrixXd 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, row_t1)
for (int i = 0; i < rowMax; i++) { for (int i = 0; i < rowMax; i++) {
@ -407,7 +413,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double timestep,
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
bcTop, bcBottom, colMax, i, sx, sy); bcTop, bcBottom, colMax, i, sx, sy);
SparseLU<SparseMatrix<double>> solver; Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
row_t1 = solverFunc(A, b); row_t1 = solverFunc(A, b);

24
src/Boundary.cpp
View File

@ -5,8 +5,6 @@
#include <stdexcept> #include <stdexcept>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
using namespace std;
BoundaryElement::BoundaryElement() { BoundaryElement::BoundaryElement() {
this->type = BC_TYPE_CLOSED; this->type = BC_TYPE_CLOSED;
@ -38,22 +36,22 @@ double BoundaryElement::getValue() { return this->value; }
Boundary::Boundary(Grid grid) : grid(grid) { Boundary::Boundary(Grid grid) : grid(grid) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
this->boundaries = vector<vector<BoundaryElement>>( this->boundaries = std::vector<std::vector<BoundaryElement>>(
2); // in 1D only left and right boundary 2); // in 1D only left and right boundary
this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement());
this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement());
} else if (grid.getDim() == 2) { } else if (grid.getDim() == 2) {
this->boundaries = vector<vector<BoundaryElement>>(4); this->boundaries = std::vector<std::vector<BoundaryElement>>(4);
this->boundaries[BC_SIDE_LEFT] = this->boundaries[BC_SIDE_LEFT] =
vector<BoundaryElement>(grid.getRow(), BoundaryElement()); std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
this->boundaries[BC_SIDE_RIGHT] = this->boundaries[BC_SIDE_RIGHT] =
vector<BoundaryElement>(grid.getRow(), BoundaryElement()); std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
this->boundaries[BC_SIDE_TOP] = this->boundaries[BC_SIDE_TOP] =
vector<BoundaryElement>(grid.getCol(), BoundaryElement()); std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
this->boundaries[BC_SIDE_BOTTOM] = this->boundaries[BC_SIDE_BOTTOM] =
vector<BoundaryElement>(grid.getCol(), BoundaryElement()); std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
} }
} }
@ -72,7 +70,7 @@ void Boundary::setBoundarySideClosed(BC_SIDE side) {
} else { } else {
n = grid.getCol(); n = grid.getCol();
} }
this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement()); this->boundaries[side] = std::vector<BoundaryElement>(n, BoundaryElement());
} }
void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
@ -90,7 +88,7 @@ void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
} else { } else {
n = grid.getCol(); n = grid.getCol();
} }
this->boundaries[side] = vector<BoundaryElement>(n, BoundaryElement(value)); this->boundaries[side] = std::vector<BoundaryElement>(n, BoundaryElement(value));
} }
void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) {
@ -111,7 +109,7 @@ void Boundary::setBoundaryElementConstant(BC_SIDE side, int index,
this->boundaries[side][index].setValue(value); this->boundaries[side][index].setValue(value);
} }
const vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) { const std::vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
if (grid.getDim() == 1) { if (grid.getDim() == 1) {
if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
throw_invalid_argument( throw_invalid_argument(
@ -122,9 +120,9 @@ const vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
return this->boundaries[side]; return this->boundaries[side];
} }
VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
int length = boundaries[side].size(); int length = boundaries[side].size();
VectorXd values(length); Eigen::VectorXd values(length);
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) {

2
src/CMakeLists.txt
View File

@ -6,4 +6,4 @@ if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND)
target_link_libraries(tug OpenMP::OpenMP_CXX) target_link_libraries(tug OpenMP::OpenMP_CXX)
endif() endif()
target_include_directories(tug PUBLIC ../include) target_include_directories(tug PUBLIC ${PROJECT_SOURCE_DIR}/include)

6
src/FTCS.cpp
View File

@ -13,8 +13,6 @@
#include <omp.h> #include <omp.h>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
using namespace std;
// calculates horizontal change on one cell independent of boundary type // calculates horizontal change on one cell independent of boundary type
static double calcHorizontalChange(Grid &grid, int &row, int &col) { static double calcHorizontalChange(Grid &grid, int &row, int &col) {
@ -222,7 +220,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
double deltaCol = grid.getDeltaCol(); double deltaCol = grid.getDeltaCol();
// matrix for concentrations at time t+1 // matrix for concentrations at time t+1
MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0); Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(1, colMax, 0);
// only one row in 1D case -> row constant at index 0 // only one row in 1D case -> row constant at index 0
int row = 0; int row = 0;
@ -262,7 +260,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
double deltaCol = grid.getDeltaCol(); double deltaCol = grid.getDeltaCol();
// matrix for concentrations at time t+1 // matrix for concentrations at time t+1
MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(rowMax, colMax, 0);
// inner cells // inner cells
// these are independent of the boundary condition type // these are independent of the boundary condition type

24
src/Grid.cpp
View File

@ -15,8 +15,8 @@ Grid::Grid(int length) {
this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
this->dim = 1; this->dim = 1;
this->concentrations = MatrixXd::Constant(1, col, 20); this->concentrations = Eigen::MatrixXd::Constant(1, col, 20);
this->alphaX = MatrixXd::Constant(1, col, 1); this->alphaX = Eigen::MatrixXd::Constant(1, col, 1);
} }
Grid::Grid(int row, int col) { Grid::Grid(int row, int col) {
@ -33,12 +33,12 @@ Grid::Grid(int row, int col) {
this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
this->dim = 2; this->dim = 2;
this->concentrations = MatrixXd::Constant(row, col, 20); this->concentrations = Eigen::MatrixXd::Constant(row, col, 20);
this->alphaX = MatrixXd::Constant(row, col, 1); this->alphaX = Eigen::MatrixXd::Constant(row, col, 1);
this->alphaY = MatrixXd::Constant(row, col, 1); this->alphaY = Eigen::MatrixXd::Constant(row, col, 1);
} }
void Grid::setConcentrations(MatrixXd concentrations) { void Grid::setConcentrations(Eigen::MatrixXd concentrations) {
if (concentrations.rows() != this->row || if (concentrations.rows() != this->row ||
concentrations.cols() != this->col) { concentrations.cols() != this->col) {
throw_invalid_argument( throw_invalid_argument(
@ -48,9 +48,9 @@ void Grid::setConcentrations(MatrixXd concentrations) {
this->concentrations = concentrations; this->concentrations = concentrations;
} }
const MatrixXd Grid::getConcentrations() { return this->concentrations; } const Eigen::MatrixXd Grid::getConcentrations() { return this->concentrations; }
void Grid::setAlpha(MatrixXd alpha) { void Grid::setAlpha(Eigen::MatrixXd alpha) {
if (dim != 1) { if (dim != 1) {
throw_invalid_argument("Grid is not one dimensional, you should probably " throw_invalid_argument("Grid is not one dimensional, you should probably "
"use 2D setter function!"); "use 2D setter function!");
@ -63,7 +63,7 @@ void Grid::setAlpha(MatrixXd alpha) {
this->alphaX = alpha; this->alphaX = alpha;
} }
void Grid::setAlpha(MatrixXd alphaX, MatrixXd alphaY) { void Grid::setAlpha(Eigen::MatrixXd alphaX, Eigen::MatrixXd alphaY) {
if (dim != 2) { if (dim != 2) {
throw_invalid_argument("Grid is not two dimensional, you should probably " throw_invalid_argument("Grid is not two dimensional, you should probably "
"use 1D setter function!"); "use 1D setter function!");
@ -81,7 +81,7 @@ void Grid::setAlpha(MatrixXd alphaX, MatrixXd alphaY) {
this->alphaY = alphaY; this->alphaY = alphaY;
} }
const MatrixXd Grid::getAlpha() { const Eigen::MatrixXd Grid::getAlpha() {
if (dim != 1) { if (dim != 1) {
throw_invalid_argument("Grid is not one dimensional, you should probably " throw_invalid_argument("Grid is not one dimensional, you should probably "
"use either getAlphaX() or getAlphaY()!"); "use either getAlphaX() or getAlphaY()!");
@ -90,7 +90,7 @@ const MatrixXd Grid::getAlpha() {
return this->alphaX; return this->alphaX;
} }
const MatrixXd Grid::getAlphaX() { const Eigen::MatrixXd Grid::getAlphaX() {
if (dim != 2) { if (dim != 2) {
throw_invalid_argument( throw_invalid_argument(
"Grid is not two dimensional, you should probably use getAlpha()!"); "Grid is not two dimensional, you should probably use getAlpha()!");
@ -99,7 +99,7 @@ const MatrixXd Grid::getAlphaX() {
return this->alphaX; return this->alphaX;
} }
const MatrixXd Grid::getAlphaY() { const Eigen::MatrixXd Grid::getAlphaY() {
if (dim != 2) { if (dim != 2) {
throw_invalid_argument( throw_invalid_argument(
"Grid is not two dimensional, you should probably use getAlpha()!"); "Grid is not two dimensional, you should probably use getAlpha()!");

94
src/Simulation.cpp
View File

@ -63,7 +63,7 @@ void Simulation::setTimestep(double timestep) {
double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
double minDeltaSquare; double minDeltaSquare;
double maxAlphaX, maxAlphaY, maxAlpha; double maxAlphaX, maxAlphaY, maxAlpha;
string dim; std::string dim;
if (grid.getDim() == 2) { if (grid.getDim() == 2) {
dim = "2D"; dim = "2D";
@ -91,31 +91,31 @@ void Simulation::setTimestep(double timestep) {
// not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha * // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha *
// ((1/deltaRowSquare) + (1/deltaColSquare))); // ((1/deltaRowSquare) + (1/deltaColSquare)));
string approachPrefix = std::string approachPrefix =
(approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl
<< endl; << std::endl;
cout << approachPrefix << "_" << dim << " :: required dt=" << timestep std::cout << approachPrefix << "_" << dim << " :: required dt=" << timestep
<< endl; << std::endl;
if (timestep > cfl) { if (timestep > cfl) {
this->innerIterations = (int)ceil(timestep / cfl); this->innerIterations = (int)ceil(timestep / cfl);
this->timestep = timestep / (double)innerIterations; this->timestep = timestep / (double)innerIterations;
cerr << "Warning :: Timestep was adjusted, because of stability " std::cerr << "Warning :: Timestep was adjusted, because of stability "
"conditions. Time duration was approximately preserved by " "conditions. Time duration was approximately preserved by "
"adjusting internal number of iterations." "adjusting internal number of iterations."
<< endl; << std::endl;
cout << approachPrefix << "_" << dim << " :: Required " std::cout << approachPrefix << "_" << dim << " :: Required "
<< this->innerIterations << this->innerIterations
<< " inner iterations with dt=" << this->timestep << endl; << " inner iterations with dt=" << this->timestep << std::endl;
} else { } else {
this->timestep = timestep; this->timestep = timestep;
cout << approachPrefix << "_" << dim std::cout << approachPrefix << "_" << dim
<< " :: No inner iterations required, dt=" << timestep << endl; << " :: No inner iterations required, dt=" << timestep << std::endl;
} }
} else { } else {
@ -134,10 +134,10 @@ void Simulation::setIterations(int iterations) {
void Simulation::setSolver(SOLVER solver) { void Simulation::setSolver(SOLVER solver) {
if (this->approach == FTCS_APPROACH) { if (this->approach == FTCS_APPROACH) {
cerr << "Warning: Solver was set, but FTCS approach initialized. Setting " std::cerr << "Warning: Solver was set, but FTCS approach initialized. Setting "
"the solver " "the solver "
"is thus without effect." "is thus without effect."
<< endl; << std::endl;
} }
this->solver = solver; this->solver = solver;
@ -148,9 +148,9 @@ void Simulation::setNumberThreads(int numThreads) {
this->numThreads = numThreads; this->numThreads = numThreads;
} else { } else {
int maxThreadNumber = omp_get_num_procs(); int maxThreadNumber = omp_get_num_procs();
string outputMessage = std::string outputMessage =
"Number of threads exceeds the number of processor cores (" + "Number of threads exceeds the number of processor cores (" +
to_string(maxThreadNumber) + ") or is less than 1."; std::to_string(maxThreadNumber) + ") or is less than 1.";
throw_invalid_argument(outputMessage); throw_invalid_argument(outputMessage);
} }
@ -159,28 +159,28 @@ void Simulation::setNumberThreads(int numThreads) {
int Simulation::getIterations() { return this->iterations; } int Simulation::getIterations() { return this->iterations; }
void Simulation::printConcentrationsConsole() { void Simulation::printConcentrationsConsole() {
cout << grid.getConcentrations() << endl; std::cout << grid.getConcentrations() << std::endl;
cout << endl; std::cout << std::endl;
} }
string Simulation::createCSVfile() { std::string Simulation::createCSVfile() {
ofstream file; std::ofstream file;
int appendIdent = 0; int appendIdent = 0;
string appendIdentString; std::string appendIdentString;
// string approachString = (approach == 0) ? "FTCS" : "BTCS"; // string approachString = (approach == 0) ? "FTCS" : "BTCS";
string approachString = std::string approachString =
(approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
string row = to_string(grid.getRow()); std::string row = std::to_string(grid.getRow());
string col = to_string(grid.getCol()); std::string col = std::to_string(grid.getCol());
string numIterations = to_string(iterations); std::string numIterations = std::to_string(iterations);
string filename = std::string filename =
approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
while (filesystem::exists(filename)) { while (std::filesystem::exists(filename)) {
appendIdent += 1; appendIdent += 1;
appendIdentString = to_string(appendIdent); appendIdentString = std::to_string(appendIdent);
filename = approachString + "_" + row + "_" + col + "_" + numIterations + filename = approachString + "_" + row + "_" + col + "_" + numIterations +
"-" + appendIdentString + ".csv"; "-" + appendIdentString + ".csv";
} }
@ -193,16 +193,16 @@ string Simulation::createCSVfile() {
// adds lines at the beginning of verbose output csv that represent the // adds lines at the beginning of verbose output csv that represent the
// boundary conditions and their values -1 in case of closed boundary // boundary conditions and their values -1 in case of closed boundary
if (csv_output == CSV_OUTPUT_XTREME) { if (csv_output == CSV_OUTPUT_XTREME) {
IOFormat one_row(StreamPrecision, DontAlignCols, "", " "); Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", " ");
file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row)
<< endl; // boundary left << std::endl; // boundary left
file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row)
<< endl; // boundary right << std::endl; // boundary right
file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row)
<< endl; // boundary top << std::endl; // boundary top
file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row)
<< endl; // boundary bottom << std::endl; // boundary bottom
file << endl << endl; file << std::endl << std::endl;
} }
file.close(); file.close();
@ -210,17 +210,17 @@ string Simulation::createCSVfile() {
return filename; return filename;
} }
void Simulation::printConcentrationsCSV(string filename) { void Simulation::printConcentrationsCSV(std::string filename) {
ofstream file; std::ofstream file;
file.open(filename, std::ios_base::app); file.open(filename, std::ios_base::app);
if (!file) { if (!file) {
exit(1); exit(1);
} }
IOFormat do_not_align(StreamPrecision, DontAlignCols); Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols);
file << grid.getConcentrations().format(do_not_align) << endl; file << grid.getConcentrations().format(do_not_align) << std::endl;
file << endl << endl; file << std::endl << std::endl;
file.close(); file.close();
} }
@ -232,7 +232,7 @@ void Simulation::run() {
throw_invalid_argument("Number of iterations are not set!"); throw_invalid_argument("Number of iterations are not set!");
} }
string filename; std::string filename;
if (this->console_output > CONSOLE_OUTPUT_OFF) { if (this->console_output > CONSOLE_OUTPUT_OFF) {
printConcentrationsConsole(); printConcentrationsConsole();
} }
@ -295,9 +295,9 @@ void Simulation::run() {
// TODO this implementation is very inefficient! // TODO this implementation is very inefficient!
// a separate implementation that sets up a specific tridiagonal matrix for // a separate implementation that sets up a specific tridiagonal matrix for
// Crank-Nicolson would be better // Crank-Nicolson would be better
MatrixXd concentrations; Eigen::MatrixXd concentrations;
MatrixXd concentrationsFTCS; Eigen::MatrixXd concentrationsFTCS;
MatrixXd concentrationsResult; Eigen::MatrixXd concentrationsResult;
for (int i = 0; i < iterations * innerIterations; i++) { for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole(); printConcentrationsConsole();
@ -328,10 +328,10 @@ void Simulation::run() {
printConcentrationsCSV(filename); printConcentrationsCSV(filename);
} }
if (this->time_measure > TIME_MEASURE_OFF) { if (this->time_measure > TIME_MEASURE_OFF) {
string approachString = std::string approachString =
(approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI"); (approach == 0) ? "FTCS" : ((approach == 1) ? "BTCS" : "CRNI");
string dimString = (grid.getDim() == 1) ? "-1D" : "-2D"; std::string dimString = (grid.getDim() == 1) ? "-1D" : "-2D";
cout << approachString << dimString << ":: run() finished in " std::cout << approachString << dimString << ":: run() finished in "
<< milliseconds.count() << "ms" << endl; << milliseconds.count() << "ms" << std::endl;
} }
} }