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

BREAKING CHANGE: Rework Grid definition

Browse Source 
Now the API does not rely on `Grid` structure but lay a wrapper around
an existing memory region, which defines for example a diffusion grid.

All simulation steps are done in place.

The user has to make sure the memory existing throughout the lifetime of
a simulation grid.
This commit is contained in:
Max Luebke 2025-01-31 15:43:37 +01:00
parent 5c68f8b6b2
commit 3612dcf034
12 changed files with 1099 additions and 1156 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
include/tug/Boundary.hpp
View File

@ -7,12 +7,13 @@
#ifndef BOUNDARY_H_
#define BOUNDARY_H_
#include "Grid.hpp"
#include "tug/Core/TugUtils.hpp"
#include <Eigen/Dense>
#include <cstddef>
#include <cstdint>
#include <map>
#include <stdexcept>
#include <utility>
#include <vector>
@ -45,7 +46,7 @@ public:
* BC_TYPE_CLOSED, where the value takes -1 and does not hold any
* physical meaning.
*/
BoundaryElement(){};
BoundaryElement() {};
/**
* @brief Construct a new Boundary Element object for the constant case.
@ -115,7 +116,7 @@ public:
*
* @param length Length of the grid
*/
Boundary(std::uint32_t length) : Boundary(Grid<T>(length)){};
Boundary(std::uint32_t length) : Boundary(1, length) {};
/**
* @brief Creates a boundary object for a 2D grid
@ -124,17 +125,7 @@ public:
* @param cols Number of columns of the grid
*/
Boundary(std::uint32_t rows, std::uint32_t cols)
: Boundary(Grid<T>(rows, cols)){};
/**
* @brief Creates a boundary object based on the passed grid object and
* initializes the boundaries as closed.
*
* @param grid Grid object on the basis of which the simulation takes place
* and from which the dimensions (in 2D case) are taken.
*/
Boundary(const Grid<T> &grid)
: dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) {
: dim(rows == 1 ? 1 : 2), cols(cols), rows(rows) {
if (this->dim == 1) {
this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(
2); // in 1D only left and right boundary
@ -153,8 +144,37 @@ public:
this->boundaries[BC_SIDE_BOTTOM] =
std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
}
}
};
/**
* @brief Creates a boundary object based on the passed grid object and
* initializes the boundaries as closed.
*
* @param grid Grid object on the basis of which the simulation takes place
* and from which the dimensions (in 2D case) are taken.
*/
// Boundary(const Grid<T> &grid)
// : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) {
// if (this->dim == 1) {
// this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(
// 2); // in 1D only left and right boundary
//
// this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement<T>());
// this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement<T>());
// } else if (this->dim == 2) {
// this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(4);
//
// this->boundaries[BC_SIDE_LEFT] =
// std::vector<BoundaryElement<T>>(this->rows, BoundaryElement<T>());
// this->boundaries[BC_SIDE_RIGHT] =
// std::vector<BoundaryElement<T>>(this->rows, BoundaryElement<T>());
// this->boundaries[BC_SIDE_TOP] =
// std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
// this->boundaries[BC_SIDE_BOTTOM] =
// std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
// }
// }
//
/**
* @brief Sets all elements of the specified boundary side to the boundary
* condition closed.

243
include/tug/Core/BaseSimulation.hpp
View File

@ -1,5 +1,6 @@
#pragma once
#include "tug/Boundary.hpp"
#include <cstddef>
#include <cstdint>
#include <tug/Core/Matrix.hpp>
@ -38,14 +39,25 @@ enum class TIME_MEASURE {
ON /*!< print time measure after last iteration */
};
/**
* @brief A base class for simulation grids.
*
* This class provides a base implementation for simulation grids, including
* methods for setting and getting grid dimensions, domain sizes, and output
* options. It also includes methods for running simulations, which must be
* implemented by derived classes.
*
* @tparam T The type of the elements in the grid.
*/
template <typename T> class BaseSimulationGrid {
protected:
private:
CSV_OUTPUT csv_output{CSV_OUTPUT::OFF};
CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT::OFF};
TIME_MEASURE time_measure{TIME_MEASURE::OFF};
int iterations{1};
RowMajMatMap<T> concentration_matrix;
RowMajMatMap<T> concentrationMatrix;
Boundary<T> boundaryConditions;
const std::uint8_t dim;
@ -53,29 +65,196 @@ protected:
T delta_row;
public:
/**
* @brief Constructs a BaseSimulationGrid from a given RowMajMat object.
*
* This constructor initializes a BaseSimulationGrid using the data, number of
* rows, and number of columns from the provided RowMajMat object.
*
* @tparam T The type of the elements in the RowMajMat.
* @param origin The RowMajMat object from which to initialize the
* BaseSimulationGrid.
*/
BaseSimulationGrid(RowMajMat<T> &origin)
: BaseSimulationGrid(origin.data(), origin.rows(), origin.cols()) {}
/**
* @brief Constructs a BaseSimulationGrid object.
*
* @tparam T The type of the data elements.
* @param data Pointer to the data array.
* @param rows Number of rows in the grid.
* @param cols Number of columns in the grid.
*
* Initializes the concentration_matrix with the provided data, rows, and
* columns. Sets delta_col and delta_row to 1. Determines the dimension (dim)
* based on the number of rows: if rows == 1, dim is set to 1; otherwise, dim
* is set to 2.
*/
BaseSimulationGrid(T *data, std::size_t rows, std::size_t cols)
: concentrationMatrix(data, rows, cols), boundaryConditions(rows, cols),
delta_col(1), delta_row(1), dim(rows == 1 ? 1 : 2) {}
/**
* @brief Constructs a BaseSimulationGrid with a single dimension.
*
* This constructor initializes a BaseSimulationGrid object with the provided
* data and length. It assumes the grid has only one dimension.
*
* @param data Pointer to the data array.
* @param length The length of the data array.
*/
BaseSimulationGrid(T *data, std::size_t length)
: BaseSimulationGrid(data, 1, length) {}
template <typename EigenType>
BaseSimulationGrid(const EigenType &origin)
: BaseSimulationGrid(origin.data(), origin.rows(), origin.cols()) {}
/**
* @brief Overloaded function call operator to access elements in a
* one-dimensional grid.
*
* This operator provides access to elements in the concentration matrix using
* a single index. It asserts that the grid is one-dimensional before
* accessing the element.
*
* @tparam T The type of elements in the concentration matrix.
* @param index The index of the element to access.
* @return A reference to the element at the specified index in the
* concentration matrix.
*/
constexpr T &operator()(std::size_t index) {
tug_assert(dim == 1, "Grid is not one dimensional, use 2D index operator!");
BaseSimulationGrid(T *data, std::size_t rows, std::size_t cols)
: concentration_matrix(data, rows, cols), delta_col(1), delta_row(1),
dim(rows == 1 ? 1 : 2) {}
return concentrationMatrix(index);
}
std::size_t rows() const { return concentration_matrix.rows(); }
std::size_t cols() const { return concentration_matrix.cols(); }
/**
* @brief Overloaded function call operator to access elements in a 2D
* concentration matrix.
*
* This operator allows accessing elements in the concentration matrix using
* row and column indices. It asserts that the grid is two-dimensional before
* accessing the element.
*
* @param row The row index of the element to access.
* @param col The column index of the element to access.
* @return A reference to the element at the specified row and column in the
* concentration matrix.
*/
constexpr T &operator()(std::size_t row, std::size_t col) {
tug_assert(dim == 2, "Grid is not two dimensional, use 1D index operator!");
return concentrationMatrix(row, col);
}
/**
* @brief Retrieves the concentration matrix.
*
* @tparam T The data type of the elements in the concentration matrix.
* @return RowMajMat<T>& Reference to the concentration matrix.
*/
RowMajMatMap<T> &getConcentrationMatrix() { return concentrationMatrix; }
const RowMajMatMap<T> &getConcentrationMatrix() const {
return concentrationMatrix;
}
/**
* @brief Retrieves the boundary conditions for the simulation.
*
* @tparam T The type parameter for the Boundary class.
* @return Boundary<T>& A reference to the boundary conditions.
*/
Boundary<T> &getBoundaryConditions() { return boundaryConditions; }
const Boundary<T> &getBoundaryConditions() const {
return boundaryConditions;
}
/**
* @brief Retrieves the dimension value.
*
* @return The dimension value as an 8-bit unsigned integer.
*/
std::uint8_t getDim() const { return dim; }
/**
* @brief Returns the number of rows in the concentration matrix.
*
* @return std::size_t The number of rows in the concentration matrix.
*/
std::size_t rows() const { return concentrationMatrix.rows(); }
/**
* @brief Get the number of columns in the concentration matrix.
*
* @return std::size_t The number of columns in the concentration matrix.
*/
std::size_t cols() const { return concentrationMatrix.cols(); }
/**
* @brief Returns the cell size in meter of the x-direction.
*
* This function returns the value of the delta column, which is used
* to represent the difference or change in the column value.
*
* @return T The cell size in meter of the x-direction.
*/
T deltaCol() const { return delta_col; }
/**
* @brief Returns the cell size in meter of the y-direction.
*
* This function asserts that the grid is two-dimensional. If the grid is not
* two-dimensional, an assertion error is raised with the message "Grid is not
* two dimensional, there is no delta in y-direction!".
*
* @return The cell size in meter of the y-direction.
*/
T deltaRow() const {
tug_assert(
dim == 1,
dim == 2,
"Grid is not two dimensional, there is no delta in y-direction!");
return delta_row;
}
/**
* @brief Computes the domain size in the X direction.
*
* This function calculates the size of the domain in the X direction by
* multiplying the column spacing (delta_col) by the number of columns (cols).
*
* @return The size of the domain in the X direction.
*/
T domainX() const { return delta_col * cols(); }
/**
* @brief Returns the size of the domain in the y-direction.
*
* This function calculates the size of the domain in the y-direction
* by multiplying the row spacing (delta_row) by the number of rows.
* It asserts that the grid is two-dimensional before performing the
* calculation.
*
* @return The size of the domain in the y-direction.
*/
T domainY() const {
tug_assert(
dim == 2,
"Grid is not two dimensional, there is no domain in y-direction!");
return delta_row * rows();
}
/**
* @brief Sets the domain length for a one-dimensional grid.
*
* This function sets the domain length for a one-dimensional grid and
* calculates the column width (delta_col) based on the given domain length
* and the number of columns. It asserts that the grid is one-dimensional and
* that the given domain length is positive.
*
* @param domain_length The length of the domain. Must be positive.
*/
void setDomain(T domain_length) {
tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!");
tug_assert(domain_length > 0, "Given domain length is not positive!");
@ -83,6 +262,17 @@ public:
delta_col = domain_length / cols();
}
/**
* @brief Sets the domain size for a 2D grid simulation.
*
* This function sets the domain size in the x and y directions for a
* two-dimensional grid simulation. It asserts that the grid is indeed
* two-dimensional and that the provided domain sizes are positive.
*
* @tparam T The type of the domain size parameters.
* @param domain_row The size of the domain in the y-direction.
* @param domain_col The size of the domain in the x-direction.
*/
void setDomain(T domain_row, T domain_col) {
tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!");
tug_assert(domain_col > 0,
@ -110,6 +300,15 @@ public:
*/
void setOutputCSV(CSV_OUTPUT csv_output) { this->csv_output = csv_output; }
/**
* @brief Retrieves the CSV output.
*
* This function returns the CSV output associated with the simulation.
*
* @return CSV_OUTPUT The CSV output of the simulation.
*/
constexpr CSV_OUTPUT getOutputCSV() const { return this->csv_output; }
/**
* @brief Set the options for outputting information to the console. Off by
* default.
@ -127,6 +326,17 @@ public:
this->console_output = console_output;
}
/**
* @brief Retrieves the console output.
*
* This function returns the current state of the console output.
*
* @return CONSOLE_OUTPUT The current console output.
*/
constexpr CONSOLE_OUTPUT getOutputConsole() const {
return this->console_output;
}
/**
* @brief Set the Time Measure option. Off by default.
*
@ -140,6 +350,13 @@ public:
this->time_measure = time_measure;
}
/**
* @brief Retrieves the current time measurement.
*
* @return TIME_MEASURE The current time measurement.
*/
constexpr TIME_MEASURE getTimeMeasure() const { return this->time_measure; }
/**
* @brief Set the desired iterations to be calculated. A value greater
* than zero must be specified here. Setting iterations is required.
@ -165,5 +382,7 @@ public:
* parameters.
*/
virtual void run() = 0;
virtual void setTimestep(T timestep) = 0;
};
} // namespace tug
} // namespace tug

84
include/tug/Core/Numeric/BTCS.hpp
View File

@ -10,15 +10,16 @@
#ifndef BTCS_H_
#define BTCS_H_
#include "../Matrix.hpp"
#include "../TugUtils.hpp"
#include <cstddef>
#include <tug/Boundary.hpp>
#include <tug/Grid.hpp>
#include <tug/Core/Matrix.hpp>
#include <tug/Core/Numeric/SimulationInput.hpp>
#include <utility>
#include <vector>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#ifdef _OPENMP
#include <omp.h>
#else
@ -159,9 +160,9 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc,
// creates a solution vector for next time step from the current state of
// concentrations
template <class T>
template <class T, class EigenType>
static Eigen::VectorX<T>
createSolutionVector(const RowMajMat<T> &concentrations,
createSolutionVector(const EigenType &concentrations,
const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY,
const std::vector<BoundaryElement<T>> &bcLeft,
const std::vector<BoundaryElement<T>> &bcRight,
@ -351,25 +352,27 @@ static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
// BTCS solution for 1D grid
template <class T>
static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
static void BTCS_1D(SimulationInput<T> &input,
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b)) {
int length = grid.getCol();
T sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol());
const std::size_t &length = input.colMax;
T sx = input.timestep / (input.deltaCol * input.deltaCol);
Eigen::VectorX<T> concentrations_t1(length);
Eigen::SparseMatrix<T> A;
Eigen::VectorX<T> b(length);
const auto &alpha = grid.getAlphaX();
const auto &alpha = input.alphaX;
const auto &bc = input.boundaries;
const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
const auto inner_bc = bc.getInnerBoundaryRow(0);
RowMajMat<T> &concentrations = grid.getConcentrations();
RowMajMatMap<T> &concentrations = input.concentrations;
int rowIndex = 0;
A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex,
sx); // this is exactly same as in 2D
@ -396,29 +399,31 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
// BTCS solution for 2D grid
template <class T>
static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
static void BTCS_2D(SimulationInput<T> &input,
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b),
int numThreads) {
int rowMax = grid.getRow();
int colMax = grid.getCol();
T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
const std::size_t &rowMax = input.rowMax;
const std::size_t &colMax = input.colMax;
const T sx = input.timestep / (2 * input.deltaCol * input.deltaCol);
const T sy = input.timestep / (2 * input.deltaRow * input.deltaRow);
RowMajMat<T> concentrations_t1(rowMax, colMax);
Eigen::SparseMatrix<T> A;
Eigen::VectorX<T> b;
RowMajMat<T> alphaX = grid.getAlphaX();
RowMajMat<T> alphaY = grid.getAlphaY();
const RowMajMat<T> &alphaX = input.alphaX;
const RowMajMat<T> &alphaY = input.alphaY;
const auto &bc = input.boundaries;
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);
RowMajMat<T> &concentrations = grid.getConcentrations();
RowMajMatMap<T> &concentrations = input.concentrations;
#pragma omp parallel for num_threads(numThreads) private(A, b)
for (int i = 0; i < rowMax; i++) {
@ -432,44 +437,43 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
}
concentrations_t1.transposeInPlace();
alphaX.transposeInPlace();
alphaY.transposeInPlace();
const RowMajMat<T> alphaX_t = alphaX.transpose();
const RowMajMat<T> alphaY_t = alphaY.transpose();
#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
A = createCoeffMatrix(alphaY, bcTop, bcBottom, inner_bc, rowMax, i, sy);
b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
bcLeft, bcRight, inner_bc, rowMax, i, sy, sx);
A = createCoeffMatrix(alphaY_t, bcTop, bcBottom, inner_bc, rowMax, i, sy);
b = createSolutionVector(concentrations_t1, alphaY_t, alphaX_t, bcTop,
bcBottom, bcLeft, bcRight, inner_bc, rowMax, i, sy,
sx);
concentrations.col(i) = solverFunc(A, b);
}
}
// entry point for EigenLU solver; differentiate between 1D and 2D grid
template <class T>
void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
} else if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads);
template <class T> void BTCS_LU(SimulationInput<T> &input, int numThreads) {
tug_assert(input.dim <= 2,
"Error: Only 1- and 2-dimensional grids are defined!");
if (input.dim == 1) {
BTCS_1D(input, EigenLUAlgorithm);
} else {
throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
BTCS_2D(input.dim, EigenLUAlgorithm, numThreads);
}
}
// entry point for Thomas algorithm solver; differentiate 1D and 2D grid
template <class T>
void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
} else if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads);
template <class T> void BTCS_Thomas(SimulationInput<T> &input, int numThreads) {
tug_assert(input.dim <= 2,
"Error: Only 1- and 2-dimensional grids are defined!");
if (input.dim == 1) {
BTCS_1D(input, ThomasAlgorithm);
} else {
throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
BTCS_2D(input, ThomasAlgorithm, numThreads);
}
}
} // namespace tug

863
include/tug/Core/Numeric/FTCS.hpp
View File

@ -8,11 +8,12 @@
#ifndef FTCS_H_
#define FTCS_H_
#include "../TugUtils.hpp"
#include "tug/Core/Matrix.hpp"
#include "tug/Core/TugUtils.hpp"
#include <cstddef>
#include <cstring>
#include <tug/Boundary.hpp>
#include <tug/Core/Matrix.hpp>
#include <tug/Core/Numeric/SimulationInput.hpp>
#ifdef _OPENMP
#include <omp.h>
@ -22,214 +23,258 @@
namespace tug {
template <class T>
constexpr T calcChangeInner(T conc_c, T conc_left, T conc_right, T alpha_c,
T alpha_left, T alpha_right) {
const T alpha_center_left = calcAlphaIntercell(alpha_left, alpha_c);
const T alpha_center_right = calcAlphaIntercell(alpha_right, alpha_c);
return alpha_center_left * conc_left -
(alpha_center_left + alpha_center_right) * conc_c +
alpha_center_right * conc_right;
}
// calculates horizontal change on one cell independent of boundary type
template <class T>
static inline T calcHorizontalChange(Grid<T> &grid, int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col + 1) -
(calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) +
calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col))) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col - 1);
}
// calculates vertical change on one cell independent of boundary type
template <class T>
static inline T calcVerticalChange(Grid<T> &grid, int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) +
calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
grid.getAlphaY()(row, col))) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row - 1, col);
}
// template <class T>
// static inline T calcHorizontalChange(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaX, int &row,
// int &col) {
// return calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) *
// concentrations(row, col + 1) -
// (calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) +
// calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col))) *
// concentrations(row, col) +
// calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) *
// concentrations(row, col - 1);
// }
//
// // calculates vertical change on one cell independent of boundary type
// template <class T>
// static inline T calcVerticalChange(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaY, int &row,
// int &col) {
// return calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) *
// concentrations(row + 1, col) -
// (calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) +
// calcAlphaIntercell(alphaY(row - 1, col), alphaY(row, col))) *
// concentrations(row, col) +
// calcAlphaIntercell(alphaY(row - 1, col), alphaY(row, col)) *
// concentrations(row - 1, col);
// }
// calculates horizontal change on one cell with a constant left boundary
// template <class T>
// static inline T calcHorizontalChangeLeftBoundaryConstant(
// const RowMajMatMap<T> &concentrations, const RowMajMatMap<T> &alphaX,
// const Boundary<T> &bc, int &row, int &col) {
// return calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) *
// concentrations(row, col + 1) -
// (calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) +
// 2 * alphaX(row, col)) *
// concentrations(row, col) +
// 2 * alphaX(row, col) * bc.getBoundaryElementValue(BC_SIDE_LEFT,
// row);
// }
//
// // calculates horizontal change on one cell with a closed left boundary
// template <class T>
// static inline T
// calcHorizontalChangeLeftBoundaryClosed(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaX, int
// &row, int &col) {
// return calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) *
// (concentrations(row, col + 1) - concentrations(row, col));
// }
//
// // checks boundary condition type for a cell on the left edge of grid
// template <class T>
// static inline T
// calcHorizontalChangeLeftBoundary(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaX, Boundary<T>
// &bc, int &row, int &col) {
// if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) {
// return calcHorizontalChangeLeftBoundaryConstant(concentrations, alphaX,
// bc,
// row, col);
// } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED)
// {
// return calcHorizontalChangeLeftBoundaryClosed(concentrations, alphaX,
// row,
// col);
// } else {
// throw_invalid_argument("Undefined Boundary Condition Type!");
// }
// }
//
// // calculates horizontal change on one cell with a constant right boundary
// template <class T>
// static inline T
// calcHorizontalChangeRightBoundaryConstant(const RowMajMatMap<T>
// &concentrations,
// const RowMajMatMap<T> &alphaX,
// Boundary<T> &bc, int &row, int
// &col) {
// return 2 * alphaX(row, col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT,
// row) -
// (calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) +
// 2 * alphaX(row, col)) *
// concentrations(row, col) +
// calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) *
// concentrations(row, col - 1);
// }
//
// // calculates horizontal change on one cell with a closed right boundary
// template <class T>
// static inline T
// calcHorizontalChangeRightBoundaryClosed(const RowMajMatMap<T>
// &concentrations,
// const RowMajMatMap<T> &alphaX, int
// &row, int &col) {
// return -(calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) *
// (concentrations(row, col) - concentrations(row, col - 1)));
// }
//
// // checks boundary condition type for a cell on the right edge of grid
// template <class T>
// static inline T
// calcHorizontalChangeRightBoundary(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaX,
// Boundary<T> &bc, int &row, int &col) {
// if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) {
// return calcHorizontalChangeRightBoundaryConstant(concentrations, alphaX,
// bc,
// row, col);
// } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CLOSED)
// {
// return calcHorizontalChangeRightBoundaryClosed(concentrations, alphaX,
// row,
// col);
// } else {
// throw_invalid_argument("Undefined Boundary Condition Type!");
// }
// }
//
// // calculates vertical change on one cell with a constant top boundary
// template <class T>
// static inline T
// calcVerticalChangeTopBoundaryConstant(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaY,
// Boundary<T> &bc, int &row, int &col) {
// return calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) *
// concentrations(row + 1, col) -
// (calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) +
// 2 * alphaY(row, col)) *
// concentrations(row, col) +
// 2 * alphaY(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col);
// }
//
// // calculates vertical change on one cell with a closed top boundary
// template <class T>
// static inline T
// calcVerticalChangeTopBoundaryClosed(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaY, int &row,
// int &col) {
// return calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) *
// (concentrations(row + 1, col) - concentrations(row, col));
// }
//
// // checks boundary condition type for a cell on the top edge of grid
// template <class T>
// static inline T
// calcVerticalChangeTopBoundary(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaY, Boundary<T> &bc,
// int &row, int &col) {
// if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
// return calcVerticalChangeTopBoundaryConstant(concentrations, alphaY, bc,
// row, col);
// } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
// return calcVerticalChangeTopBoundaryClosed(concentrations, alphaY, row,
// col);
// } else {
// throw_invalid_argument("Undefined Boundary Condition Type!");
// }
// }
//
// // calculates vertical change on one cell with a constant bottom boundary
// template <class T>
// static inline T
// calcVerticalChangeBottomBoundaryConstant(const RowMajMatMap<T>
// &concentrations,
// const RowMajMatMap<T> &alphaY,
// Boundary<T> &bc, int &row, int &col)
// {
// return 2 * alphaY(row, col) *
// bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) -
// (calcAlphaIntercell(alphaY(row, col), alphaY(row - 1, col)) +
// 2 * alphaY(row, col)) *
// concentrations(row, col) +
// calcAlphaIntercell(alphaY(row, col), alphaY(row - 1, col)) *
// concentrations(row - 1, col);
// }
//
// // calculates vertical change on one cell with a closed bottom boundary
// template <class T>
// static inline T
// calcVerticalChangeBottomBoundaryClosed(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaY, int
// &row, int &col) {
// return -(calcAlphaIntercell(alphaY(row, col), alphaY(row - 1, col)) *
// (concentrations(row, col) - concentrations(row - 1, col)));
// }
//
// // checks boundary condition type for a cell on the bottom edge of grid
// template <class T>
// static inline T
// calcVerticalChangeBottomBoundary(const RowMajMatMap<T> &concentrations,
// const RowMajMatMap<T> &alphaY, Boundary<T>
// &bc, int &row, int &col) {
// if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
// return calcVerticalChangeBottomBoundaryConstant(concentrations, alphaY,
// bc,
// row, col);
// } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) ==
// BC_TYPE_CLOSED) {
// return calcVerticalChangeBottomBoundaryClosed(concentrations, alphaY,
// row,
// col);
// } else {
// throw_invalid_argument("Undefined Boundary Condition Type!");
// }
// }
template <class T>
static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid,
Boundary<T> &bc,
int &row, int &col) {
constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center,
T alpha_neighbor, const BoundaryElement<T> &bc) {
const T alpha_center_neighbor =
calcAlphaIntercell(alpha_center, alpha_neighbor);
const T &conc_boundary = bc.getValue();
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col + 1) -
(calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) +
2 * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col) +
2 * grid.getAlphaX()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
}
// calculates horizontal change on one cell with a closed left boundary
template <class T>
static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col + 1) -
grid.getConcentrations()(row, col));
}
// checks boundary condition type for a cell on the left edge of grid
template <class T>
static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &bc,
int &row, int &col) {
if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED) {
return calcHorizontalChangeLeftBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
switch (bc.getType()) {
case BC_TYPE_CONSTANT: {
return 2 * alpha_center * conc_boundary -
(alpha_center_neighbor + 2 * alpha_center) * conc_c +
alpha_center_neighbor * conc_neighbor;
}
}
// calculates horizontal change on one cell with a constant right boundary
template <class T>
static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid,
Boundary<T> &bc,
int &row, int &col) {
return 2 * grid.getAlphaX()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) -
(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) +
2 * grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
grid.getConcentrations()(row, col - 1);
}
// calculates horizontal change on one cell with a closed right boundary
template <class T>
static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
(grid.getConcentrations()(row, col) -
grid.getConcentrations()(row, col - 1)));
}
// checks boundary condition type for a cell on the right edge of grid
template <class T>
static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid,
Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CLOSED) {
return calcHorizontalChangeRightBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
case BC_TYPE_CLOSED: {
return (alpha_center_neighbor * (conc_neighbor - conc_c));
}
}
// calculates vertical change on one cell with a constant top boundary
template <class T>
static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid,
Boundary<T> &bc, int &row,
int &col) {
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row + 1, col) -
(calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) +
2 * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row, col) +
2 * grid.getAlphaY()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_TOP, col);
}
// calculates vertical change on one cell with a closed top boundary
template <class T>
static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
(grid.getConcentrations()(row + 1, col) -
grid.getConcentrations()(row, col));
}
// checks boundary condition type for a cell on the top edge of grid
template <class T>
static inline T calcVerticalChangeTopBoundary(Grid<T> &grid, Boundary<T> &bc,
int &row, int &col) {
if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
return calcVerticalChangeTopBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
}
}
// calculates vertical change on one cell with a constant bottom boundary
template <class T>
static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid,
Boundary<T> &bc,
int &row, int &col) {
return 2 * grid.getAlphaY()(row, col) *
bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) -
(calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) +
2 * grid.getAlphaY()(row, col)) *
grid.getConcentrations()(row, col) +
calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) *
grid.getConcentrations()(row - 1, col);
}
// calculates vertical change on one cell with a closed bottom boundary
template <class T>
static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
return -(calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) *
(grid.getConcentrations()(row, col) -
grid.getConcentrations()(row - 1, col)));
}
// checks boundary condition type for a cell on the bottom edge of grid
template <class T>
static inline T calcVerticalChangeBottomBoundary(Grid<T> &grid, Boundary<T> &bc,
int &row, int &col) {
if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
return calcVerticalChangeBottomBoundaryClosed(grid, row, col);
} else {
throw_invalid_argument("Undefined Boundary Condition Type!");
}
tug_assert(false, "Undefined Boundary Condition Type!");
}
// FTCS solution for 1D grid
template <class T>
static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
int colMax = grid.getCol();
T deltaCol = grid.getDeltaCol();
template <class T> static void FTCS_1D(SimulationInput<T> &input) {
const std::size_t &colMax = input.colMax;
const T &deltaCol = input.deltaCol;
const T &timestep = input.timestep;
RowMajMat<T> &concentrations_grid = grid.getConcentrations();
RowMajMatMap<T> &concentrations_grid = input.concentrations;
// matrix for concentrations at time t+1
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(1, colMax, 0);
RowMajMat<T> concentrations_t1 = concentrations_grid;
const auto &alphaX = input.alphaX;
const auto &bc = input.boundaries;
// only one row in 1D case -> row constant at index 0
int row = 0;
@ -237,169 +282,325 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
// inner cells
// independent of boundary condition type
for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) = concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChange(grid, row, col));
const T &conc_c = concentrations_grid(row, col);
const T &conc_left = concentrations_grid(row, col - 1);
const T &conc_right = concentrations_grid(row, col + 1);
const T &alpha_c = alphaX(row, col);
const T &alpha_left = alphaX(row, col - 1);
const T &alpha_right = alphaX(row, col + 1);
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
calcChangeInner(conc_c, conc_left, conc_right, alpha_c, alpha_left,
alpha_right);
}
// left boundary; hold column constant at index 0
int col = 0;
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col));
{
int col = 0;
const T &conc_c = concentrations_grid(row, col);
const T &conc_right = concentrations_grid(row, col + 1);
const T &alpha_c = alphaX(row, col);
const T &alpha_right = alphaX(row, col + 1);
const BoundaryElement<T> &bc_element =
input.boundaries.getBoundaryElement(BC_SIDE_LEFT, row);
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
calcChangeBoundary(conc_c, conc_right, alpha_c, alpha_right,
bc_element);
}
// right boundary; hold column constant at max index
col = colMax - 1;
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col));
{
int col = colMax - 1;
const T &conc_c = concentrations_grid(row, col);
const T &conc_left = concentrations_grid(row, col - 1);
const T &alpha_c = alphaX(row, col);
const T &alpha_left = alphaX(row, col - 1);
const BoundaryElement<T> &bc_element =
bc.getBoundaryElement(BC_SIDE_RIGHT, row);
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
calcChangeBoundary(conc_c, conc_left, alpha_c, alpha_left,
bc_element);
}
// overwrite obsolete concentrations
std::memcpy(concentrations_grid.data(), concentrations_t1.data(),
colMax * sizeof(T));
concentrations_grid = concentrations_t1;
}
// FTCS solution for 2D grid
template <class T>
static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
int numThreads) {
int rowMax = grid.getRow();
int colMax = grid.getCol();
T deltaRow = grid.getDeltaRow();
T deltaCol = grid.getDeltaCol();
static void FTCS_2D(SimulationInput<T> &input, int numThreads) {
const std::size_t &rowMax = input.rowMax;
const std::size_t &colMax = input.colMax;
const T &deltaRow = input.deltaRow;
const T &deltaCol = input.deltaCol;
const T &timestep = input.timestep;
RowMajMat<T> &concentrations_grid = grid.getConcentrations();
RowMajMatMap<T> &concentrations_grid = input.concentrations;
// matrix for concentrations at time t+1
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
RowMajMat<T> concentrations_t1 = concentrations_grid;
const auto &alphaX = input.alphaX;
const auto &alphaY = input.alphaY;
const auto &bc = input.boundaries;
const T sx = timestep / (deltaCol * deltaCol);
const T sy = timestep / (deltaRow * deltaRow);
#pragma omp parallel for num_threads(numThreads)
for (std::size_t row_i = 0; row_i < rowMax; row_i++) {
for (std::size_t col_i = 0; col_i < colMax; col_i++) {
// horizontal change
T horizontal_change;
{
const T &conc_c = concentrations_grid(row_i, col_i);
const T &alpha_c = alphaX(row_i, col_i);
if (col_i == 0 || col_i == colMax - 1) {
// left or right boundary
const T &conc_neigbor =
concentrations_grid(row_i, col_i == 0 ? col_i + 1 : col_i - 1);
const T &alpha_neigbor =
alphaX(row_i, col_i == 0 ? col_i + 1 : col_i - 1);
const BoundaryElement<T> &bc_element = bc.getBoundaryElement(
col_i == 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i);
horizontal_change = calcChangeBoundary(conc_c, conc_neigbor, alpha_c,
alpha_neigbor, bc_element);
} else {
// inner cell
const T &conc_left = concentrations_grid(row_i, col_i - 1);
const T &conc_right = concentrations_grid(row_i, col_i + 1);
const T &alpha_left = alphaX(row_i, col_i - 1);
const T &alpha_right = alphaX(row_i, col_i + 1);
horizontal_change = calcChangeInner(conc_c, conc_left, conc_right,
alpha_c, alpha_left, alpha_right);
}
}
// vertical change
T vertical_change;
{
const T &conc_c = concentrations_grid(row_i, col_i);
const T &alpha_c = alphaY(row_i, col_i);
if (row_i == 0 || row_i == rowMax - 1) {
// top or bottom boundary
const T &conc_neigbor =
concentrations_grid(row_i == 0 ? row_i + 1 : row_i - 1, col_i);
const T &alpha_neigbor =
alphaY(row_i == 0 ? row_i + 1 : row_i - 1, col_i);
const BoundaryElement<T> &bc_element = bc.getBoundaryElement(
row_i == 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i);
vertical_change = calcChangeBoundary(conc_c, conc_neigbor, alpha_c,
alpha_neigbor, bc_element);
} else {
// inner cell
const T &conc_bottom = concentrations_grid(row_i - 1, col_i);
const T &conc_top = concentrations_grid(row_i + 1, col_i);
const T &alpha_bottom = alphaY(row_i - 1, col_i);
const T &alpha_top = alphaY(row_i + 1, col_i);
vertical_change = calcChangeInner(conc_c, conc_bottom, conc_top,
alpha_c, alpha_bottom, alpha_top);
}
}
concentrations_t1(row_i, col_i) = concentrations_grid(row_i, col_i) +
sx * horizontal_change +
sy * vertical_change;
}
}
// inner cells
// these are independent of the boundary condition type
// omp_set_num_threads(10);
#pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax - 1; row++) {
for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) = concentrations_grid(row, col) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChange(grid, row, col)) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChange(grid, row, col));
}
}
// #pragma omp parallel for num_threads(numThreads)
// for (int row = 1; row < rowMax - 1; row++) {
// for (int col = 1; col < colMax - 1; col++) {
// const T &conc_c = concentrations_grid(row, col);
// const T &conc_left = concentrations_grid(row, col - 1);
// const T &conc_right = concentrations_grid(row, col + 1);
// const T &conc_top = concentrations_grid(row + 1, col);
// const T &conc_bottom = concentrations_grid(row - 1, col);
//
// const T &alpha_c = alphaX(row, col);
// const T &alpha_left = alphaX(row, col - 1);
// const T &alpha_right = alphaX(row, col + 1);
// const T &alpha_top = alphaY(row + 1, col);
// const T &alpha_bottom = alphaY(row - 1, col);
//
// const T horizontal_change = calcChangesInner(
// conc_c, conc_left, conc_right, alpha_c, alpha_left, alpha_right);
//
// const T vertical_change = calcChangesInner(
// conc_c, conc_bottom, conc_top, alpha_c, alpha_bottom, alpha_top);
//
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaRow * deltaRow) * vertical_change +
// timestep / (deltaCol * deltaCol) * horizontal_change;
// }
// }
// boundary conditions
// left without corners / looping over rows
// hold column constant at index 0
int col = 0;
#pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax - 1; row++) {
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
}
// right without corners / looping over rows
// hold column constant at max index
col = colMax - 1;
#pragma omp parallel for num_threads(numThreads)
for (int row = 1; row < rowMax - 1; row++) {
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
}
// top without corners / looping over columns
// hold row constant at index 0
int row = 0;
#pragma omp parallel for num_threads(numThreads)
for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeTopBoundary(grid, bc, row, col)) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChange(grid, row, col));
}
// bottom without corners / looping over columns
// hold row constant at max index
row = rowMax - 1;
#pragma omp parallel for num_threads(numThreads)
for (int col = 1; col < colMax - 1; col++) {
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col)) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChange(grid, row, col));
}
// corner top left
// hold row and column constant at 0
row = 0;
col = 0;
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeTopBoundary(grid, bc, row, col));
// corner top right
// hold row constant at 0 and column constant at max index
row = 0;
col = colMax - 1;
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeTopBoundary(grid, bc, row, col));
// corner bottom left
// hold row constant at max index and column constant at 0
row = rowMax - 1;
col = 0;
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col));
// corner bottom right
// hold row and column constant at max index
row = rowMax - 1;
col = colMax - 1;
concentrations_t1(row, col) =
concentrations_grid(row, col) +
timestep / (deltaCol * deltaCol) *
(calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
timestep / (deltaRow * deltaRow) *
(calcVerticalChangeBottomBoundary(grid, bc, row, col));
// int col = 0;
// #pragma omp parallel for num_threads(numThreads)
// for (int row = 1; row < rowMax - 1; row++) {
// const T horizontal_change = calcChangeBoundary(
// concentrations_grid(row, col), concentrations_grid(row, col + 1),
// alphaX(row, col), alphaX(row, col + 1),
// bc.getBoundaryElement(BC_SIDE_LEFT, row));
//
// const T vertical_change = calcChangesInner(
// concentrations_grid(row, col), concentrations_grid(row - 1, col),
// concentrations_grid(row + 1, col), alphaX(row, col),
// alphaY(row - 1, col), alphaY(row + 1, col));
//
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaCol * deltaCol) * horizontal_change +
// timestep / (deltaRow * deltaRow) * vertical_change;
// }
//
// // right without corners / looping over rows
// // hold column constant at max index
// col = colMax - 1;
// #pragma omp parallel for num_threads(numThreads)
// for (int row = 1; row < rowMax - 1; row++) {
//
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaCol * deltaCol) *
// (calcHorizontalChangeRightBoundary(concentrations_grid, alphaX,
// bc,
// row, col)) +
// timestep / (deltaRow * deltaRow) *
// (calcVerticalChange(concentrations_grid, alphaY, row, col));
// }
//
// // top without corners / looping over columns
// // hold row constant at index 0
// int row = 0;
// #pragma omp parallel for num_threads(numThreads)
// for (int col = 1; col < colMax - 1; col++) {
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaRow * deltaRow) *
// (calcVerticalChangeTopBoundary(concentrations_grid, alphaY, bc,
// row,
// col)) +
// timestep / (deltaCol * deltaCol) *
// (calcHorizontalChange(concentrations_grid, alphaX, row, col));
// }
//
// // bottom without corners / looping over columns
// // hold row constant at max index
// row = rowMax - 1;
// #pragma omp parallel for num_threads(numThreads)
// for (int col = 1; col < colMax - 1; col++) {
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaRow * deltaRow) *
// (calcVerticalChangeBottomBoundary(concentrations_grid, alphaY,
// bc,
// row, col)) +
// timestep / (deltaCol * deltaCol) *
// (calcHorizontalChange(concentrations_grid, alphaX, row, col));
// }
//
// // corner top left
// // hold row and column constant at 0
// row = 0;
// col = 0;
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaCol * deltaCol) *
// (calcHorizontalChangeLeftBoundary(concentrations_grid, alphaX,
// bc,
// row, col)) +
// timestep / (deltaRow * deltaRow) *
// (calcVerticalChangeTopBoundary(concentrations_grid, alphaY, bc,
// row,
// col));
//
// // corner top right
// // hold row constant at 0 and column constant at max index
// row = 0;
// col = colMax - 1;
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaCol * deltaCol) *
// (calcHorizontalChangeRightBoundary(concentrations_grid, alphaX,
// bc,
// row, col)) +
// timestep / (deltaRow * deltaRow) *
// (calcVerticalChangeTopBoundary(concentrations_grid, alphaY, bc,
// row,
// col));
//
// // corner bottom left
// // hold row constant at max index and column constant at 0
// row = rowMax - 1;
// col = 0;
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaCol * deltaCol) *
// (calcHorizontalChangeLeftBoundary(concentrations_grid, alphaX,
// bc,
// row, col)) +
// timestep / (deltaRow * deltaRow) *
// (calcVerticalChangeBottomBoundary(concentrations_grid, alphaY,
// bc,
// row, col));
//
// // corner bottom right
// // hold row and column constant at max index
// row = rowMax - 1;
// col = colMax - 1;
// concentrations_t1(row, col) =
// concentrations_grid(row, col) +
// timestep / (deltaCol * deltaCol) *
// (calcHorizontalChangeRightBoundary(concentrations_grid, alphaX,
// bc,
// row, col)) +
// timestep / (deltaRow * deltaRow) *
// (calcVerticalChangeBottomBoundary(concentrations_grid, alphaY,
// bc,
// row, col));
// overwrite obsolete concentrations
std::memcpy(concentrations_grid.data(), concentrations_t1.data(),
rowMax * colMax * sizeof(T));
concentrations_grid = concentrations_t1;
// }
}
// entry point; differentiate between 1D and 2D grid
template <class T>
void FTCS(Grid<T> &grid, Boundary<T> &bc, T timestep, int &numThreads) {
if (grid.getDim() == 1) {
FTCS_1D(grid, bc, timestep);
} else if (grid.getDim() == 2) {
FTCS_2D(grid, bc, timestep, numThreads);
template <class T> void FTCS(SimulationInput<T> &input, int &numThreads) {
tug_assert(input.dim <= 2,
"Error: Only 1- and 2-dimensional grids are defined!");
if (input.dim == 1) {
FTCS_1D(input);
} else {
throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!");
FTCS_2D(input, numThreads);
}
}
} // namespace tug

21
include/tug/Core/Numeric/SimulationInput.hpp Normal file
View File

@ -0,0 +1,21 @@
#pragma once
#include <tug/Boundary.hpp>
#include <tug/Core/Matrix.hpp>
namespace tug {
template <typename T> struct SimulationInput {
RowMajMatMap<T> &concentrations;
const RowMajMat<T> &alphaX;
const RowMajMat<T> &alphaY;
const Boundary<T> boundaries;
const std::uint8_t dim;
const T timestep;
const std::size_t rowMax;
const std::size_t colMax;
const T deltaRow;
const T deltaCol;
};
} // namespace tug

199
include/tug/Diffusion.hpp
View File

@ -8,8 +8,7 @@
#pragma once
#include "Boundary.hpp"
#include "Grid.hpp"
#include "tug/Core/Matrix.hpp"
#include <algorithm>
#include <filesystem>
#include <fstream>
@ -18,10 +17,9 @@
#include <string>
#include <vector>
#include "Core/Numeric/BTCS.hpp"
#include "Core/Numeric/FTCS.hpp"
#include "Core/TugUtils.hpp"
#include "tug/Core/BaseSimulation.hpp"
#include <tug/Core/BaseSimulation.hpp>
#include <tug/Core/Numeric/BTCS.hpp>
#include <tug/Core/Numeric/FTCS.hpp>
#ifdef _OPENMP
#include <omp.h>
@ -63,31 +61,77 @@ enum SOLVER {
*/
template <class T, APPROACH approach = BTCS_APPROACH,
SOLVER solver = THOMAS_ALGORITHM_SOLVER>
class Diffusion : public BaseSimulation {
class Diffusion : public BaseSimulationGrid<T> {
private:
T timestep{-1};
int innerIterations{1};
int numThreads{omp_get_num_procs()};
Grid<T> &grid;
Boundary<T> &bc;
RowMajMat<T> alphaX;
RowMajMat<T> alphaY;
const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
static constexpr T DEFAULT_ALPHA = 1E-8;
void init_alpha() {
this->alphaX =
RowMajMat<T>::Constant(this->rows(), this->cols(), DEFAULT_ALPHA);
if (this->getDim() == 2) {
this->alphaY =
RowMajMat<T>::Constant(this->rows(), this->cols(), DEFAULT_ALPHA);
}
}
public:
/**
* @brief Set up a simulation environment. The timestep and number of
* iterations must be set. For the BTCS approach, the Thomas algorithm is used
* as the default linear equation solver as this is faster for tridiagonal
* matrices. CSV output, console output and time measure are off by
* default. Also, the number of cores is set to the maximum number of cores -1
* by default.
*
* @param grid Valid grid object
* @param bc Valid boundary condition object
* @param approach Approach to solving the problem. Either FTCS or BTCS.
*/
Diffusion(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
Diffusion(RowMajMat<T> &origin) : BaseSimulationGrid<T>(origin) {
init_alpha();
}
Diffusion(T *data, int rows, int cols)
: BaseSimulationGrid<T>(data, rows, cols) {
init_alpha();
}
Diffusion(T *data, std::size_t length) : BaseSimulationGrid<T>(data, length) {
init_alpha();
}
RowMajMat<T> &getAlphaX() { return alphaX; }
RowMajMat<T> &getAlphaY() {
tug_assert(
this->getDim(),
"Grid is not two dimensional, there is no domain in y-direction!");
return alphaY;
}
void setAlphaX(const RowMajMat<T> &alphaX) { this->alphaX = alphaX; }
void setAlphaY(const RowMajMat<T> &alphaY) {
tug_assert(
this->getDim(),
"Grid is not two dimensional, there is no domain in y-direction!");
this->alphaY = alphaY;
}
// /**
// * @brief Set up a simulation environment. The timestep and number of
// * iterations must be set. For the BTCS approach, the Thomas algorithm is
// used
// * as the default linear equation solver as this is faster for tridiagonal
// * matrices. CSV output, console output and time measure are off by
// * default. Also, the number of cores is set to the maximum number of cores
// -1
// * by default.
// *
// * @param grid Valid grid object
// * @param bc Valid boundary condition object
// * @param approach Approach to solving the problem. Either FTCS or BTCS.
// */
// Diffusion(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc) {};
/**
* @brief Setting the time step for each iteration step. Time step must be
@ -95,31 +139,31 @@ public:
*
* @param timestep Valid timestep greater than zero.
*/
void setTimestep(T timestep) {
tug_assert(timestep > 0, "Timestep has to be greater than zero.");
void setTimestep(T timestep) override {
tug_assert(timestep > 0, "Timestep has to be greater than zero.");
if constexpr (approach == FTCS_APPROACH ||
approach == CRANK_NICOLSON_APPROACH) {
T cfl;
if (grid.getDim() == 1) {
if (this->getDim() == 1) {
const T deltaSquare = grid.getDeltaCol();
const T maxAlpha = grid.getAlphaX().maxCoeff();
const T deltaSquare = this->deltaCol();
const T maxAlpha = this->alphaX.maxCoeff();
// Courant-Friedrichs-Lewy condition
cfl = deltaSquare / (4 * maxAlpha);
} else if (grid.getDim() == 2) {
const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
} else if (this->getDim() == 2) {
const T deltaColSquare = this->deltaCol() * this->deltaCol();
// will be 0 if 1D, else ...
const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
const T deltaRowSquare = this->deltaRow() * this->deltaRow();
const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
const T maxAlpha =
std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff());
std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff());
cfl = minDeltaSquare / (4 * maxAlpha);
}
const std::string dim = std::to_string(grid.getDim()) + "D";
const std::string dim = std::to_string(this->getDim()) + "D";
const std::string &approachPrefix = this->approach_names[approach];
std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl
@ -183,8 +227,8 @@ public:
* @brief Outputs the current concentrations of the grid on the console.
*
*/
inline void printConcentrationsConsole() const {
std::cout << grid.getConcentrations() << std::endl;
void printConcentrationsConsole() const {
std::cout << this->getConcentrationMatrix() << std::endl;
std::cout << std::endl;
}
@ -204,9 +248,9 @@ public:
// string approachString = (approach == 0) ? "FTCS" : "BTCS";
const std::string &approachString = this->approach_names[approach];
std::string row = std::to_string(grid.getRow());
std::string col = std::to_string(grid.getCol());
std::string numIterations = std::to_string(iterations);
std::string row = std::to_string(this->rows());
std::string col = std::to_string(this->cols());
std::string numIterations = std::to_string(this->getIterations());
std::string filename =
approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
@ -225,7 +269,9 @@ public:
// adds lines at the beginning of verbose output csv that represent the
// boundary conditions and their values -1 in case of closed boundary
if (csv_output == CSV_OUTPUT_XTREME) {
if (this->getOutputCSV() == CSV_OUTPUT::XTREME) {
const auto &bc = this->getBoundaryConditions();
Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "",
" ");
file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row)
@ -260,7 +306,7 @@ public:
}
Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols);
file << grid.getConcentrations().format(do_not_align) << std::endl;
file << this->getConcentrationMatrix().format(do_not_align) << std::endl;
file << std::endl << std::endl;
file.close();
}
@ -269,30 +315,42 @@ public:
* @brief Method starts the simulation process with the previously set
* parameters.
*/
void run() {
tug_assert(this->timestep > 0, "Timestep is not set!");
tug_assert(this->iterations > 0, "Number of iterations are not set!");
void run() override {
tug_assert(this->getTimestep() > 0, "Timestep is not set!");
tug_assert(this->getIterations() > 0, "Number of iterations are not set!");
std::string filename;
if (this->console_output > CONSOLE_OUTPUT_OFF) {
if (this->getOutputConsole() > CONSOLE_OUTPUT::OFF) {
printConcentrationsConsole();
}
if (this->csv_output > CSV_OUTPUT_OFF) {
if (this->getOutputCSV() > CSV_OUTPUT::OFF) {
filename = createCSVfile();
}
auto begin = std::chrono::high_resolution_clock::now();
SimulationInput<T> sim_input = {.concentrations =
this->getConcentrationMatrix(),
.alphaX = this->getAlphaX(),
.alphaY = this->getAlphaY(),
.boundaries = this->getBoundaryConditions(),
.dim = this->getDim(),
.timestep = this->getTimestep(),
.rowMax = this->rows(),
.colMax = this->cols(),
.deltaRow = this->deltaRow(),
.deltaCol = this->deltaCol()};
if constexpr (approach == FTCS_APPROACH) { // FTCS case
for (int i = 0; i < iterations * innerIterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
for (int i = 0; i < this->getIterations() * innerIterations; i++) {
if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
printConcentrationsCSV(filename);
}
FTCS(this->grid, this->bc, this->timestep, this->numThreads);
FTCS(sim_input, this->numThreads);
// if (i % (iterations * innerIterations / 100) == 0) {
// double percentage = (double)i / ((double)iterations *
@ -305,29 +363,28 @@ public:
} else if constexpr (approach == BTCS_APPROACH) { // BTCS case
if constexpr (solver == EIGEN_LU_SOLVER) {
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
for (int i = 0; i < this->getIterations(); i++) {
if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
printConcentrationsCSV(filename);
}
BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads);
BTCS_LU(sim_input, this->numThreads);
}
} else if constexpr (solver == THOMAS_ALGORITHM_SOLVER) {
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
for (int i = 0; i < this->getIterations(); i++) {
if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
printConcentrationsCSV(filename);
}
BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
BTCS_Thomas(sim_input, this->numThreads);
}
}
} else if constexpr (approach ==
CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case
@ -339,22 +396,22 @@ public:
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) {
for (int i = 0; i < this->getIterations() * innerIterations; i++) {
if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
printConcentrationsCSV(filename);
}
concentrations = grid.getConcentrations();
concentrations = this->getConcentrationMatrix();
FTCS(this->grid, this->bc, this->timestep, this->numThreads);
concentrationsFTCS = grid.getConcentrations();
grid.setConcentrations(concentrations);
BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
concentrationsResult =
beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations();
grid.setConcentrations(concentrationsResult);
concentrationsFTCS = this->getConcentrationMatrix();
this->getConcentrationMatrix() = concentrations;
BTCS_Thomas(sim_input, this->numThreads);
concentrationsResult = beta * concentrationsFTCS +
(1 - beta) * this->getConcentrationMatrix();
this->getConcentrationMatrix() = concentrationsResult;
}
}
@ -362,15 +419,15 @@ public:
auto milliseconds =
std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
if (this->console_output > CONSOLE_OUTPUT_OFF) {
if (this->getOutputConsole() > CONSOLE_OUTPUT::OFF) {
printConcentrationsConsole();
}
if (this->csv_output > CSV_OUTPUT_OFF) {
if (this->getOutputCSV() > CSV_OUTPUT::OFF) {
printConcentrationsCSV(filename);
}
if (this->time_measure > TIME_MEASURE_OFF) {
if (this->getTimeMeasure() > TIME_MEASURE::OFF) {
const std::string &approachString = this->approach_names[approach];
const std::string dimString = std::to_string(grid.getDim()) + "D";
const std::string dimString = std::to_string(this->getDim()) + "D";
std::cout << approachString << dimString << ":: run() finished in "
<< milliseconds.count() << "ms" << std::endl;
}

342
include/tug/Grid.hpp
View File

@ -1,342 +0,0 @@
#pragma once
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients.
*
*/
#include "Core/Matrix.hpp"
#include "tug/Core/TugUtils.hpp"
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
#include <cstddef>
namespace tug {
/**
* @brief Holds a matrix with concenctration and respective matrix/matrices of
* alpha coefficients.
*
* @tparam T Type to be used for matrices, e.g. double or float
*/
template <class T> class Grid {
public:
/**
* @brief Construct a new Grid object.
*
* Constructs a new Grid object with given concentrations, defined by an
* Eigen::Matrix. The domain length is per default the same as the length. The
* alpha coefficients are set to 1. The dimensions of the grid are determined
* by the given matrix, which can also be an Eigen::Vector for a 1D-Grid.
*
* @param concentrations An Eigen3 MatrixX<T> holding the concentrations.
*/
Grid(const RowMajMat<T> &concentrations) {
if (concentrations.rows() == 1) {
this->dim = 1;
this->domainCol = static_cast<T>(concentrations.cols());
this->deltaCol = static_cast<T>(this->domainCol) /
static_cast<T>(concentrations.cols()); // -> 1
this->concentrations = concentrations;
return;
}
if (concentrations.cols() == 1) {
this->dim = 1;
this->domainCol = static_cast<T>(concentrations.rows());
this->deltaCol = static_cast<T>(this->domainCol) /
static_cast<T>(concentrations.rows()); // -> 1
this->concentrations = concentrations.transpose();
return;
}
this->dim = 2;
this->domainRow = static_cast<T>(concentrations.rows());
this->domainCol = static_cast<T>(concentrations.cols());
this->deltaRow = static_cast<T>(this->domainRow) /
static_cast<T>(concentrations.rows()); // -> 1
this->deltaCol = static_cast<T>(this->domainCol) /
static_cast<T>(concentrations.cols()); // -> 1
this->concentrations = concentrations;
// this->alphaX = RowMajMat<T>::Constant(concentrations.rows(),
// concentrations.cols(),
// MAT_INIT_VAL);
// this->alphaY = RowMajMat<T>::Constant(concentrations.rows(),
// concentrations.cols(),
// MAT_INIT_VAL);
}
/**
* @brief Construct a new Grid object.
*
* Constructs a new 1D Grid object with given concentrations, defined by a
* pointer to consecutive memory and the length of the array. The domain
* length is per default the same as the count of grid cells (length of
* array). The memory region is mapped internally, changes will affect the
* original array and the memory shall not be freed. There is no check for
* correct memory size!
*
* @param concentrations Pointer to consecutive memory holding concentrations.
* @param length Length of the array/the 1D grid.
*/
Grid(T *concentrations, std::size_t length) : dim(1) {
this->domainCol = static_cast<T>(length); // -> 1
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(length); // -> 1
this->concentrations = RowMajMatMap<T>(concentrations, 1, length);
}
/**
* @brief Construct a new Grid object.
*
* Constructs a new 2D Grid object with given concentrations, defined by a
* pointer to consecutive memory and the number of rows and columns. The
* domain size is per default the same as the number of rows and columns. The
* memory region is mapped internally, changes will affect the original array
* and the memory shall not be freed. There is no check for correct memory
* size!
*
* @param concentrations Pointer to consecutive memory holding concentrations.
* @param row Number of rows.
* @param col Number of columns.
*/
Grid(T *concentrations, std::size_t row, std::size_t col) : dim(2) {
this->domainRow = static_cast<T>(row); // -> 1
this->domainCol = static_cast<T>(col); // -> 1
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(col); // -> 1
this->deltaRow =
static_cast<T>(this->domainRow) / static_cast<T>(row); // -> 1
this->concentrations = RowMajMatMap<T>(concentrations, row, col);
}
/**
* @brief Gets the concentrations matrix for a Grid.
*
* @return An Eigen3 matrix holding the concentrations and having
* the same dimensions as the grid.
*/
auto &getConcentrations() { return this->concentrations; }
void initAlpha() {
this->alphaX = RowMajMat<T>::Constant(
this->concentrations.rows(), this->concentrations.cols(), MAT_INIT_VAL);
if (dim > 1) {
this->alphaY =
RowMajMat<T>::Constant(this->concentrations.rows(),
this->concentrations.cols(), MAT_INIT_VAL);
}
}
/**
* @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
* dimensional.
*
* @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 RowMajMat<T> &alpha) {
tug_assert(dim == 1,
"Grid is not one dimensional, use 2D setter function!");
tug_assert(
alpha.rows() == 1 && alpha.cols() == this->concentrations.cols(),
"Given matrix of alpha coefficients mismatch with Grid dimensions!");
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) {
tug_assert(dim == 1,
"Grid is not one dimensional, use 2D setter function!");
RowMajMatMap<T> map(alpha, 1, this->col);
this->alphaX = map;
}
/**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
* dimensional.
*
* @param alphaX An Eigen3 MatrixX<T> holding the alpha coefficients in
* x-direction. Matrix must be of same size as the grid.
* @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 RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY) {
tug_assert(dim == 2,
"Grid is not two dimensional, use 1D setter function!");
tug_assert(alphaX.rows() == this->concentrations.rows(),
"Alpha in x-direction "
"has wrong number of rows!");
tug_assert(alphaX.cols() == this->concentrations.cols(),
"Alpha in x-direction has wrong number of columns!");
tug_assert(alphaY.rows() == this->concentrations.rows(),
"Alpha in y-direction "
"has wrong number of rows!");
tug_assert(alphaY.cols() == this->concentrations.cols(),
"Alpha in y-direction has wrong number of columns!");
this->alphaX = alphaX;
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) {
tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!");
RowMajMatMap<T> mapX(alphaX, this->row, this->col);
RowMajMatMap<T> mapY(alphaY, this->row, this->col);
this->alphaX = mapX;
this->alphaY = mapY;
}
/**
* @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return A matrix holding the alpha coefficients in x-direction.
*/
const auto &getAlphaX() const {
tug_assert(this->alphaX.size() > 0, "AlphaX is empty!");
return this->alphaX;
}
/**
* @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return A matrix holding the alpha coefficients in y-direction.
*/
const auto &getAlphaY() const {
tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!");
tug_assert(this->alphaY.size() > 0, "AlphaY is empty!");
return this->alphaY;
}
/**
* @brief Gets the dimensions of the grid.
*
* @return Dimensions, either 1 or 2.
*/
int getDim() const { return this->dim; }
/**
* @brief Gets the number of rows of the grid.
*
* @return Number of rows.
*/
int getRow() const { return this->concentrations.rows(); }
/**
* @brief Gets the number of columns of the grid.
*
* @return Number of columns.
*/
int getCol() const { return this->concentrations.cols(); }
/**
* @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional.
*
* @param domainLength A double value of the domain length. Must be positive.
*/
void setDomain(double domainLength) {
tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!");
tug_assert(domainLength > 0, "Given domain length is not positive!");
this->domainCol = domainLength;
this->deltaCol =
double(this->domainCol) / double(this->concentrations.cols());
}
/**
* @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional.
*
* @param domainRow A double value of the domain size in y-direction. Must
* be positive.
* @param domainCol A double value of the domain size in x-direction. Must
* be positive.
*/
void setDomain(double domainRow, double domainCol) {
tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!");
tug_assert(domainCol > 0,
"Given domain size in x-direction is not positive!");
tug_assert(domainRow > 0,
"Given domain size in y-direction is not positive!");
this->domainRow = domainRow;
this->domainCol = domainCol;
this->deltaRow =
double(this->domainRow) / double(this->concentrations.rows());
this->deltaCol =
double(this->domainCol) / double(this->concentrations.cols());
}
/**
* @brief Gets the delta value in x-direction.
*
* @return Delta value in x-direction.
*/
T getDeltaCol() const { return this->deltaCol; }
/**
* @brief Gets the delta value in y-direction. Must be two dimensional grid.
*
* @return Delta value in y-direction.
*/
T getDeltaRow() const {
tug_assert(dim == 2, "Grid is not two dimensional, there is no delta in "
"y-direction!");
return this->deltaRow;
}
private:
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; // delta in y-direction (between rows)
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;
};
using Grid64 = Grid<double>;
using Grid32 = Grid<float>;
} // namespace tug

3
test/CMakeLists.txt
View File

@ -12,7 +12,6 @@ FetchContent_MakeAvailable(googletest)
add_executable(testTug
setup.cpp
testDiffusion.cpp
testGrid.cpp
testFTCS.cpp
testBoundary.cpp
)
@ -26,4 +25,4 @@ get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH)
# set relative path in header file
configure_file(testSimulation.hpp.in testSimulation.hpp)
# include test directory with generated header file from above
target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src")
target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src")

18
test/TestUtils.hpp
View File

@ -1,3 +1,4 @@
#include "tug/Core/Matrix.hpp"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <fstream>
@ -8,7 +9,7 @@
#define TUG_TEST(x) TEST(Tug, x)
inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) {
inline tug::RowMajMat<double> CSV2Eigen(std::string file2Convert) {
std::vector<double> matrixEntries;
@ -31,21 +32,20 @@ inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) {
}
}
return Eigen::Map<
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
matrixEntries.data(), matrixRowNumber,
matrixEntries.size() / matrixRowNumber);
return tug::RowMajMatMap<double>(matrixEntries.data(), matrixRowNumber,
matrixEntries.size() / matrixRowNumber);
}
inline bool checkSimilarity(Eigen::MatrixXd a, Eigen::MatrixXd b,
inline bool checkSimilarity(tug::RowMajMat<double> &a,
tug::RowMajMatMap<double> &b,
double precision = 1e-5) {
return a.isApprox(b, precision);
}
inline bool checkSimilarityV2(Eigen::MatrixXd a, Eigen::MatrixXd b,
double maxDiff) {
inline bool checkSimilarityV2(tug::RowMajMat<double> &a,
tug::RowMajMatMap<double> &b, double maxDiff) {
Eigen::MatrixXd diff = a - b;
tug::RowMajMat<double> diff = a - b;
double maxCoeff = diff.maxCoeff();
return abs(maxCoeff) < maxDiff;
}

10
test/testBoundary.cpp
View File

@ -28,12 +28,8 @@ BOUNDARY_TEST(Element) {
}
BOUNDARY_TEST(Class) {
Eigen::VectorXd conc(10);
Grid grid1D = Grid64(conc);
Eigen::MatrixXd conc2D(10, 12);
Grid grid2D = Grid64(conc2D);
Boundary boundary1D = Boundary(grid1D);
Boundary boundary2D = Boundary(grid2D);
Boundary<double> boundary1D(10);
Boundary<double> boundary2D(10, 12);
vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));
constexpr double inner_condition_value = -5;
@ -47,7 +43,6 @@ BOUNDARY_TEST(Class) {
col_ibc[0] = innerBoundary;
{
EXPECT_NO_THROW(Boundary boundary(grid1D));
EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1);
EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1);
EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
@ -76,7 +71,6 @@ BOUNDARY_TEST(Class) {
}
{
EXPECT_NO_THROW(Boundary boundary(grid1D));
EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10);
EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10);
EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12);

158
test/testDiffusion.cpp
View File

@ -1,4 +1,5 @@
#include "TestUtils.hpp"
#include "tug/Core/Matrix.hpp"
#include "gtest/gtest.h"
#include <gtest/gtest.h>
#include <stdexcept>
@ -16,18 +17,27 @@ using namespace Eigen;
using namespace std;
using namespace tug;
Grid64 setupSimulation(double timestep, int iterations) {
int row = 11;
int col = 11;
constexpr int row = 11;
constexpr int col = 11;
template <tug::APPROACH approach>
Diffusion<double, approach> setupSimulation(RowMajMat<double> &concentrations,
double timestep, int iterations) {
int domain_row = 10;
int domain_col = 10;
// Grid
MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
// RowMajMat<double> concentrations = MatrixXd::Constant(row, col, 0);
concentrations(5, 5) = 1;
Grid grid = Grid64(concentrations);
grid.setDomain(domain_row, domain_col);
Diffusion<double, approach> diffusiongrid(concentrations);
diffusiongrid.getConcentrationMatrix() = concentrations;
diffusiongrid.setDomain(domain_row, domain_col);
diffusiongrid.setTimestep(timestep);
diffusiongrid.setIterations(iterations);
diffusiongrid.setDomain(domain_row, domain_col);
MatrixXd alpha = MatrixXd::Constant(row, col, 1);
for (int i = 0; i < 5; i++) {
@ -45,9 +55,10 @@ Grid64 setupSimulation(double timestep, int iterations) {
alpha(i, j) = 0.1;
}
}
grid.setAlpha(alpha, alpha);
diffusiongrid.setAlphaX(alpha);
diffusiongrid.setAlphaY(alpha);
return grid;
return diffusiongrid;
}
constexpr double timestep = 0.001;
@ -56,127 +67,150 @@ constexpr double iterations = 7000;
DIFFUSION_TEST(EqualityFTCS) {
// set string from the header file
string test_path = testSimulationCSVDir;
MatrixXd reference = CSV2Eigen(test_path);
RowMajMat<double> reference = CSV2Eigen(test_path);
cout << "FTCS Test: " << endl;
Grid grid = setupSimulation(timestep, iterations); // Boundary
Boundary bc = Boundary(grid);
RowMajMat<double> concentrations = MatrixXd::Constant(row, col, 0);
Diffusion<double, tug::FTCS_APPROACH> sim =
setupSimulation<tug::FTCS_APPROACH>(concentrations, timestep, iterations);
// Boundary bc = Boundary(grid);
// Simulation
Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
// Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
// sim.setOutputConsole(CONSOLE_OUTPUT_ON);
sim.setTimestep(timestep);
sim.setIterations(iterations);
// sim.setTimestep(timestep);
// sim.setIterations(iterations);
sim.run();
cout << endl;
EXPECT_TRUE(checkSimilarity(reference, grid.getConcentrations(), 0.1));
EXPECT_TRUE(checkSimilarity(reference, sim.getConcentrationMatrix(), 0.1));
}
DIFFUSION_TEST(EqualityBTCS) {
// set string from the header file
string test_path = testSimulationCSVDir;
MatrixXd reference = CSV2Eigen(test_path);
RowMajMat<double> reference = CSV2Eigen(test_path);
cout << "BTCS Test: " << endl;
Grid grid = setupSimulation(timestep, iterations); // Boundary
Boundary bc = Boundary(grid);
RowMajMat<double> concentrations = MatrixXd::Constant(row, col, 0);
Diffusion<double, tug::BTCS_APPROACH> sim =
setupSimulation<tug::BTCS_APPROACH>(concentrations, timestep,
iterations); // Boundary
// Boundary bc = Boundary(grid);
// Simulation
Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
// Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
// sim.setOutputConsole(CONSOLE_OUTPUT_ON);
sim.setTimestep(timestep);
sim.setIterations(iterations);
// sim.setTimestep(timestep);
// sim.setIterations(iterations);
sim.run();
cout << endl;
EXPECT_TRUE(checkSimilarityV2(reference, grid.getConcentrations(), 0.01));
EXPECT_TRUE(checkSimilarityV2(reference, sim.getConcentrationMatrix(), 0.01));
}
DIFFUSION_TEST(InitializeEnvironment) {
int rc = 12;
Eigen::MatrixXd concentrations(rc, rc);
Grid64 grid(concentrations);
Boundary boundary(grid);
RowMajMat<double> concentrations(rc, rc);
// Grid64 grid(concentrations);
// Boundary boundary(grid);
EXPECT_NO_THROW(Diffusion sim(grid, boundary));
EXPECT_NO_FATAL_FAILURE(Diffusion<double> sim(concentrations));
}
DIFFUSION_TEST(SimulationEnvironment) {
int rc = 12;
Eigen::MatrixXd concentrations(rc, rc);
Grid64 grid(concentrations);
grid.initAlpha();
Boundary boundary(grid);
Diffusion<double, tug::FTCS_APPROACH> sim(grid, boundary);
// DIFFUSION_TEST(SimulationEnvironment) {
// int rc = 12;
// Eigen::MatrixXd concentrations(rc, rc);
// Grid64 grid(concentrations);
// grid.initAlpha();
// Boundary boundary(grid);
// Diffusion<double, tug::FTCS_APPROACH> sim(grid, boundary);
EXPECT_EQ(sim.getIterations(), 1);
// EXPECT_EQ(sim.getIterations(), 1);
EXPECT_NO_THROW(sim.setIterations(2000));
EXPECT_EQ(sim.getIterations(), 2000);
EXPECT_THROW(sim.setIterations(-300), std::invalid_argument);
// EXPECT_NO_THROW(sim.setIterations(2000));
// EXPECT_EQ(sim.getIterations(), 2000);
// EXPECT_THROW(sim.setIterations(-300), std::invalid_argument);
EXPECT_NO_THROW(sim.setTimestep(0.1));
EXPECT_DOUBLE_EQ(sim.getTimestep(), 0.1);
EXPECT_DEATH(sim.setTimestep(-0.3), ".* greater than zero.*");
}
// EXPECT_NO_THROW(sim.setTimestep(0.1));
// EXPECT_DOUBLE_EQ(sim.getTimestep(), 0.1);
// EXPECT_DEATH(sim.setTimestep(-0.3), ".* greater than zero.*");
// }
DIFFUSION_TEST(ClosedBoundaries) {
constexpr std::uint32_t nrows = 5;
constexpr std::uint32_t ncols = 5;
auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
RowMajMat<double> concentrations =
RowMajMat<double>::Constant(nrows, ncols, 1.0);
RowMajMat<double> alphax = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
RowMajMat<double> alphay = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
tug::Grid64 grid(concentrations);
Diffusion<double> sim(concentrations);
sim.getAlphaX() = alphax;
sim.getAlphaY() = alphay;
grid.setAlpha(alphax, alphay);
// tug::Grid64 grid(concentrations);
tug::Boundary bc(grid);
// grid.setAlpha(alphax, alphay);
// tug::Boundary bc(grid);
auto &bc = sim.getBoundaryConditions();
bc.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1.0);
bc.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1.0);
bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0);
bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0);
tug::Diffusion<double> sim(grid, bc);
// tug::Diffusion<double> sim(grid, bc);
sim.setTimestep(1);
sim.setIterations(1);
MatrixXd input_values(concentrations);
RowMajMat<double> input_values(concentrations);
sim.run();
EXPECT_TRUE(checkSimilarityV2(input_values, grid.getConcentrations(), 1E-12));
EXPECT_TRUE(
checkSimilarityV2(input_values, sim.getConcentrationMatrix(), 1E-12));
}
DIFFUSION_TEST(ConstantInnerCell) {
constexpr std::uint32_t nrows = 5;
constexpr std::uint32_t ncols = 5;
auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
RowMajMat<double> concentrations =
RowMajMat<double>::Constant(nrows, ncols, 1.0);
RowMajMat<double> alphax = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
RowMajMat<double> alphay = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
tug::Grid64 grid(concentrations);
grid.setAlpha(alphax, alphay);
Diffusion<double> sim(concentrations);
sim.getAlphaX() = alphax;
sim.getAlphaY() = alphay;
tug::Boundary bc(grid);
// tug::Grid64 grid(concentrations);
// grid.setAlpha(alphax, alphay);
// tug::Boundary bc(grid);
auto &bc = sim.getBoundaryConditions();
// inner
bc.setInnerBoundary(2, 2, 0);
tug::Diffusion<double> sim(grid, bc);
// tug::Diffusion<double> sim(grid, bc);
sim.setTimestep(1);
sim.setIterations(1);
MatrixXd input_values(concentrations);
sim.run();
EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 2), 0);
EXPECT_LT(grid.getConcentrations().sum(), input_values.sum());
const auto &concentrations_result = sim.getConcentrationMatrix();
EXPECT_FALSE((grid.getConcentrations().array() > 1.0).any());
EXPECT_DOUBLE_EQ(concentrations_result(2, 2), 0);
EXPECT_LT(concentrations_result.sum(), input_values.sum());
EXPECT_FALSE((grid.getConcentrations().array() < 0.0).any());
EXPECT_FALSE((concentrations_result.array() > 1.0).any());
EXPECT_FALSE((concentrations_result.array() < 0.0).any());
}

264
test/testGrid.cpp
View File

@ -1,264 +0,0 @@
#include "gtest/gtest.h"
#include <Eigen/Core>
#include <tug/Grid.hpp>
#include <gtest/gtest.h>
using namespace Eigen;
using namespace std;
using namespace tug;
#define GRID_TEST(x) TEST(Grid, x)
GRID_TEST(Grid64OneDimensional) {
int l = 12;
Eigen::VectorXd conc(l);
Grid64 grid(conc);
grid.initAlpha();
{
EXPECT_EQ(grid.getDim(), 1);
EXPECT_EQ(grid.getCol(), l);
EXPECT_EQ(grid.getCol(), l);
EXPECT_EQ(grid.getRow(), 1);
EXPECT_EQ(grid.getConcentrations().rows(), 1);
EXPECT_EQ(grid.getConcentrations().cols(), l);
EXPECT_EQ(grid.getAlphaX().rows(), 1);
EXPECT_EQ(grid.getAlphaX().cols(), l);
EXPECT_EQ(grid.getDeltaCol(), 1);
EXPECT_DEATH(grid.getAlphaY(), ".* no alphaY!.*");
EXPECT_DEATH(grid.getDeltaRow(), ".* not two dimensional, .*");
}
{
// correct alpha matrix
MatrixXd alpha = MatrixXd::Constant(1, l, 3);
EXPECT_NO_THROW(grid.setAlpha(alpha));
EXPECT_DEATH(grid.setAlpha(alpha, alpha), ".* is not two dimensional, .*");
grid.setAlpha(alpha);
EXPECT_EQ(grid.getAlphaX(), alpha);
EXPECT_NO_THROW(grid.getAlphaX());
EXPECT_DEATH(grid.getAlphaY(), ".* no alphaY!.*");
// false alpha matrix
MatrixXd wAlpha = MatrixXd::Constant(3, l, 2);
EXPECT_DEATH(grid.setAlpha(wAlpha), ".* mismatch with Grid dimensions!.*");
}
{
int d = 8;
// set 1D domain
EXPECT_NO_THROW(grid.setDomain(d));
// set 2D domain
EXPECT_DEATH(grid.setDomain(d, d), ".* not two dimensional, .*");
grid.setDomain(d);
EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(d) / double(l));
EXPECT_DEATH(grid.getDeltaRow(), ".* not two dimensional, .*");
// set too small domain
EXPECT_DEATH(grid.setDomain(-2), "Given domain length .*");
}
}
GRID_TEST(Grid64Quadratic) {
int rc = 12;
Eigen::MatrixXd conc(rc, rc);
Grid64 grid(conc);
grid.initAlpha();
{
EXPECT_EQ(grid.getDim(), 2);
EXPECT_EQ(grid.getCol(), rc);
EXPECT_EQ(grid.getRow(), rc);
EXPECT_EQ(grid.getConcentrations().rows(), rc);
EXPECT_EQ(grid.getConcentrations().cols(), rc);
EXPECT_EQ(grid.getAlphaX().rows(), rc);
EXPECT_EQ(grid.getAlphaX().cols(), rc);
EXPECT_EQ(grid.getAlphaY().rows(), rc);
EXPECT_EQ(grid.getAlphaY().cols(), rc);
EXPECT_EQ(grid.getDeltaRow(), 1);
EXPECT_EQ(grid.getDeltaCol(), 1);
}
{
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
EXPECT_DEATH(grid.setAlpha(alphax), ".* 2D setter function!.*");
grid.setAlpha(alphax, alphay);
EXPECT_EQ(grid.getAlphaX(), alphax);
EXPECT_EQ(grid.getAlphaY(), alphay);
// false alpha matrices
alphax = MatrixXd::Constant(rc + 3, rc + 1, 3);
EXPECT_DEATH(grid.setAlpha(alphax, alphay),
".*has wrong number of rows!.*");
alphay = MatrixXd::Constant(rc + 2, rc + 1, 3);
EXPECT_DEATH(grid.setAlpha(alphax, alphay),
".*has wrong number of rows!.*");
}
{
int dr = 8;
int dc = 9;
// set 1D domain
EXPECT_DEATH(grid.setDomain(dr), ".* 2D domain setter!.*");
// set 2D domain
EXPECT_NO_THROW(grid.setDomain(dr, dc));
grid.setDomain(dr, dc);
EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(rc));
EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(rc));
// set too small domain
dr = 0;
EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*");
dr = 8;
dc = 0;
EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*");
dr = -2;
EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*");
}
}
GRID_TEST(Grid64NonQuadratic) {
int r = 12;
int c = 15;
Eigen::MatrixXd conc(r, c);
Grid64 grid(conc);
grid.initAlpha();
{
EXPECT_EQ(grid.getDim(), 2);
EXPECT_EQ(grid.getCol(), c);
EXPECT_EQ(grid.getRow(), r);
EXPECT_EQ(grid.getConcentrations().rows(), r);
EXPECT_EQ(grid.getConcentrations().cols(), c);
EXPECT_EQ(grid.getAlphaX().rows(), r);
EXPECT_EQ(grid.getAlphaX().cols(), c);
EXPECT_EQ(grid.getAlphaY().rows(), r);
EXPECT_EQ(grid.getAlphaY().cols(), c);
EXPECT_EQ(grid.getDeltaRow(), 1);
EXPECT_EQ(grid.getDeltaCol(), 1);
}
{
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(r, c, 2);
MatrixXd alphay = MatrixXd::Constant(r, c, 4);
EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
grid.setAlpha(alphax, alphay);
EXPECT_EQ(grid.getAlphaX(), alphax);
EXPECT_EQ(grid.getAlphaY(), alphay);
}
{
int dr = 8;
int dc = 9;
// set 2D domain
EXPECT_NO_THROW(grid.setDomain(dr, dc));
grid.setDomain(dr, dc);
EXPECT_EQ(grid.getDeltaCol(), double(dc) / double(c));
EXPECT_EQ(grid.getDeltaRow(), double(dr) / double(r));
}
{
int r = 4;
int c = 5;
std::vector<double> concentrations(r * c);
for (int i = 0; i < r * c; i++) {
concentrations[i] = i;
}
Grid64 grid(concentrations.data(), r, c);
grid.initAlpha();
{
EXPECT_EQ(grid.getDim(), 2);
EXPECT_EQ(grid.getCol(), c);
EXPECT_EQ(grid.getRow(), r);
EXPECT_EQ(grid.getConcentrations().rows(), r);
EXPECT_EQ(grid.getConcentrations().cols(), c);
EXPECT_EQ(grid.getAlphaX().rows(), r);
EXPECT_EQ(grid.getAlphaX().cols(), c);
EXPECT_EQ(grid.getAlphaY().rows(), r);
EXPECT_EQ(grid.getAlphaY().cols(), c);
EXPECT_EQ(grid.getDeltaRow(), 1);
EXPECT_EQ(grid.getDeltaCol(), 1);
}
{
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(r, c, 2);
MatrixXd alphay = MatrixXd::Constant(r, c, 4);
EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
EXPECT_DEATH(grid.setAlpha(alphax), ".* 2D setter function!.*");
grid.setAlpha(alphax, alphay);
EXPECT_EQ(grid.getAlphaX(), alphax);
EXPECT_EQ(grid.getAlphaY(), alphay);
// false alpha matrices
alphax = MatrixXd::Constant(r + 3, c + 1, 3);
EXPECT_DEATH(grid.setAlpha(alphax, alphay),
".*has wrong number of rows!.*");
alphay = MatrixXd::Constant(r + 2, c + 1, 5);
EXPECT_DEATH(grid.setAlpha(alphax, alphay),
".*has wrong number of rows!.*");
{
int dr = 8;
int dc = 9;
// set 1D domain
EXPECT_DEATH(grid.setDomain(dr), ".* 2D domain setter!.*");
// set 2D domain
EXPECT_NO_THROW(grid.setDomain(dr, dc));
grid.setDomain(dr, dc);
EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(c));
EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(r));
}
{
auto &concentrations = grid.getConcentrations();
for (int i = 0; i < r; i++) {
for (int j = 0; j < c; j++) {
concentrations(i, j) = i * c + j;
}
}
EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0);
EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1);
EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c);
EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0);
EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1);
EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c);
EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 1), 2 * c + 1);
}
}
}
}