From 6e388f3d990e2dc9848e94b0942cea5700faae8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 5 Sep 2023 16:38:54 +0200 Subject: [PATCH 01/17] write input of thomas algortithm to file --- src/BTCSv2.cpp | 27 +++++++++++++++++++++++++++ src/CMakeLists.txt | 6 ++++++ 2 files changed, 33 insertions(+) diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index 49e4d1c..705b152 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -11,6 +11,11 @@ #include #include +#ifdef WRITE_THOMAS_CSV +#include +#include +#endif + #define NUM_THREADS_BTCS 10 using namespace Eigen; @@ -295,6 +300,28 @@ static VectorXd ThomasAlgorithm(SparseMatrix &A, VectorXd &b) { a_diag[n - 1] = A.coeff(n - 1, n - 2); b_diag[n - 1] = A.coeff(n - 1, n - 1); + // HACK: write CSV to file +#ifdef WRITE_THOMAS_CSV + #include + #include + static std::uint32_t file_index = 0; + std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; + + std::ofstream out_file; + + out_file.open(file_name, std::ofstream::trunc | std::ofstream::out); + + // print header + out_file << "Aa, Ab, Ac, b\n"; + + // iterate through all elements + for (std::size_t i = 0; i < n; i++) { + out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " << b[i] << "\n"; + } + + out_file.close(); +#endif + // start solving - c_diag and x_vec are overwritten n--; c_diag[0] /= b_diag[0]; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b7aceb2..3cc79d4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,11 @@ add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCSv2.cpp) +option(TUG_WRITE_CSV "Write CSV during Thomas algorithm with consecutive numbers for each call" OFF) + +IF(TUG_WRITE_CSV) + target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV) +endif() + target_link_libraries(tug Eigen3::Eigen) if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND) From 7d05320f24de39abbafdbf00fa331340c2eb3c06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 5 Sep 2023 16:42:05 +0200 Subject: [PATCH 02/17] apply format changes (LLVM) --- src/BTCSv2.cpp | 693 +++++++++++++++++++++++++------------------------ 1 file changed, 358 insertions(+), 335 deletions(-) diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index 705b152..12b9844 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -1,464 +1,487 @@ /** * @file BTCSv2.cpp - * @brief Implementation of heterogenous BTCS (backward time-centered space) solution - * of diffusion equation in 1D and 2D space. Internally the alternating-direction - * implicit (ADI) method is used. Version 2, because Version 1 was an - * implementation for the homogeneous BTCS solution. - * + * @brief Implementation of heterogenous BTCS (backward time-centered space) + * solution of diffusion equation in 1D and 2D space. Internally the + * alternating-direction implicit (ADI) method is used. Version 2, because + * Version 1 was an implementation for the homogeneous BTCS solution. + * */ #include "FTCS.cpp" -#include #include +#include #ifdef WRITE_THOMAS_CSV -#include -#include +#include +#include #endif #define NUM_THREADS_BTCS 10 using namespace Eigen; - // calculates coefficient for left boundary in constant case -static tuple calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) { - double centerCoeff; - double rightCoeff; +static tuple +calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) { + double centerCoeff; + double rightCoeff; - centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)) - + 2 * alpha(rowIndex,0)); - rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); + centerCoeff = + 1 + sx * (calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)) + + 2 * alpha(rowIndex, 0)); + rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); - return {centerCoeff, rightCoeff}; + return {centerCoeff, rightCoeff}; } - // calculates coefficient for left boundary in closed case -static tuple calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) { - double centerCoeff; - double rightCoeff; +static tuple +calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) { + double centerCoeff; + double rightCoeff; - centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); - rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1)); + centerCoeff = + 1 + sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); + rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); - return {centerCoeff, rightCoeff}; + return {centerCoeff, rightCoeff}; } - // calculates coefficient for right boundary in constant case -static tuple calcRightBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, int n, double sx) { - double leftCoeff; - double centerCoeff; +static tuple calcRightBoundaryCoeffConstant(MatrixXd &alpha, + int rowIndex, int n, + double sx) { + double leftCoeff; + double centerCoeff; - leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); - centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)) - + 2 * alpha(rowIndex,n)); + leftCoeff = + -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); + centerCoeff = + 1 + sx * (calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)) + + 2 * alpha(rowIndex, n)); - return {leftCoeff, centerCoeff}; + return {leftCoeff, centerCoeff}; } - // calculates coefficient for right boundary in closed case -static tuple calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) { - double leftCoeff; - double centerCoeff; +static tuple +calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) { + double leftCoeff; + double centerCoeff; - leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); - centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n)); + leftCoeff = + -sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); + centerCoeff = + 1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)); - return {leftCoeff, centerCoeff}; + return {leftCoeff, centerCoeff}; } - // creates coefficient matrix for next time step from alphas in x-direction -static SparseMatrix createCoeffMatrix(MatrixXd &alpha, vector &bcLeft, vector &bcRight, int numCols, int rowIndex, double sx) { +static SparseMatrix createCoeffMatrix(MatrixXd &alpha, + vector &bcLeft, + vector &bcRight, + int numCols, int rowIndex, + double sx) { - // square matrix of column^2 dimension for the coefficients - SparseMatrix cm(numCols, numCols); - cm.reserve(VectorXi::Constant(numCols, 3)); + // square matrix of column^2 dimension for the coefficients + SparseMatrix cm(numCols, numCols); + cm.reserve(VectorXi::Constant(numCols, 3)); - // left column - BC_TYPE type = bcLeft[rowIndex].getType(); - if (type == BC_TYPE_CONSTANT) { - auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx); - cm.insert(0,0) = centerCoeffTop; - cm.insert(0,1) = rightCoeffTop; - } else if (type == BC_TYPE_CLOSED) { - auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx); - cm.insert(0,0) = centerCoeffTop; - cm.insert(0,1) = rightCoeffTop; - } else { - throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!"); - } + // left column + BC_TYPE type = bcLeft[rowIndex].getType(); + if (type == BC_TYPE_CONSTANT) { + auto [centerCoeffTop, rightCoeffTop] = + calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx); + cm.insert(0, 0) = centerCoeffTop; + cm.insert(0, 1) = rightCoeffTop; + } else if (type == BC_TYPE_CLOSED) { + auto [centerCoeffTop, rightCoeffTop] = + calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx); + cm.insert(0, 0) = centerCoeffTop; + cm.insert(0, 1) = rightCoeffTop; + } else { + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Left or Top!"); + } - // inner columns - int n = numCols-1; - for (int i = 1; i < n; i++) { - cm.insert(i,i-1) = -sx * calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i)); - cm.insert(i,i) = 1 + sx * ( - calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1)) - + calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i)) - ) - ; - cm.insert(i,i+1) = -sx * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1)); - } + // inner columns + int n = numCols - 1; + for (int i = 1; i < n; i++) { + cm.insert(i, i - 1) = + -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)); + cm.insert(i, i) = + 1 + + sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) + + calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i))); + cm.insert(i, i + 1) = + -sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)); + } - // right column - type = bcRight[rowIndex].getType(); - if (type == BC_TYPE_CONSTANT) { - auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx); - cm.insert(n,n-1) = leftCoeffBottom; - cm.insert(n,n) = centerCoeffBottom; - } else if (type == BC_TYPE_CLOSED) { - auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx); - cm.insert(n,n-1) = leftCoeffBottom; - cm.insert(n,n) = centerCoeffBottom; - } else { - throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!"); - } + // right column + type = bcRight[rowIndex].getType(); + if (type == BC_TYPE_CONSTANT) { + auto [leftCoeffBottom, centerCoeffBottom] = + calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx); + cm.insert(n, n - 1) = leftCoeffBottom; + cm.insert(n, n) = centerCoeffBottom; + } else if (type == BC_TYPE_CLOSED) { + auto [leftCoeffBottom, centerCoeffBottom] = + calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx); + cm.insert(n, n - 1) = leftCoeffBottom; + cm.insert(n, n) = centerCoeffBottom; + } else { + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Right or Bottom!"); + } - cm.makeCompressed(); // important for Eigen solver + cm.makeCompressed(); // important for Eigen solver - return cm; + return cm; } - // calculates explicity concentration at top boundary in constant case -static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations, - MatrixXd &alpha, vector &bcTop, int rowIndex, int i, double sy) { - double c; +static double calcExplicitConcentrationsTopBoundaryConstant( + MatrixXd &concentrations, MatrixXd &alpha, vector &bcTop, + int rowIndex, int i, double sy) { + double c; - c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) - * concentrations(rowIndex,i) - + ( - 1 - sy * ( - calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) - + 2 * alpha(rowIndex,i) - ) - ) * concentrations(rowIndex,i) - + sy * alpha(rowIndex,i) * bcTop[i].getValue(); + c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) * + concentrations(rowIndex, i) + + (1 - + sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) + + 2 * alpha(rowIndex, i))) * + concentrations(rowIndex, i) + + sy * alpha(rowIndex, i) * bcTop[i].getValue(); - return c; + return c; } - // calculates explicit concentration at top boundary in closed case -static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations, - MatrixXd &alpha, int rowIndex, int i, double sy) { - double c; +static double calcExplicitConcentrationsTopBoundaryClosed( + MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) { + double c; - c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) - * concentrations(rowIndex,i) - + ( - 1 - sy * ( - calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i)) - ) - ) * concentrations(rowIndex,i); + c = sy * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)) * + concentrations(rowIndex, i) + + (1 - + sy * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex + 1, i)))) * + concentrations(rowIndex, i); - return c; + return c; } - // calculates explicit concentration at bottom boundary in constant case -static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations, - MatrixXd &alpha, vector &bcBottom, int rowIndex, int i, double sy) { - double c; +static double calcExplicitConcentrationsBottomBoundaryConstant( + MatrixXd &concentrations, MatrixXd &alpha, + vector &bcBottom, int rowIndex, int i, double sy) { + double c; - c = sy * alpha(rowIndex,i) * bcBottom[i].getValue() - + ( - 1 - sy * ( - 2 * alpha(rowIndex,i) - + calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) - ) - ) * concentrations(rowIndex,i) - + sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) - * concentrations(rowIndex-1,i); + c = sy * alpha(rowIndex, i) * bcBottom[i].getValue() + + (1 - + sy * (2 * alpha(rowIndex, i) + + calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) * + concentrations(rowIndex, i) + + sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) * + concentrations(rowIndex - 1, i); - return c; + return c; } - // calculates explicit concentration at bottom boundary in closed case -static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations, - MatrixXd &alpha, int rowIndex, int i, double sy) { - double c; +static double calcExplicitConcentrationsBottomBoundaryClosed( + MatrixXd &concentrations, MatrixXd &alpha, int rowIndex, int i, double sy) { + double c; - c = ( - 1 - sy * ( - + calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) - ) - ) * concentrations(rowIndex,i) - + sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) - * concentrations(rowIndex-1,i); + c = (1 - + sy * (+calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)))) * + concentrations(rowIndex, i) + + sy * calcAlphaIntercell(alpha(rowIndex - 1, i), alpha(rowIndex, i)) * + concentrations(rowIndex - 1, i); - return c; + return c; } +// creates a solution vector for next time step from the current state of +// concentrations +static VectorXd createSolutionVector( + MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, + vector &bcLeft, vector &bcRight, + vector &bcTop, vector &bcBottom, + int length, int rowIndex, double sx, double sy) { -// creates a solution vector for next time step from the current state of concentrations -static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, - vector &bcLeft, vector &bcRight, - vector &bcTop, vector &bcBottom, - int length, int rowIndex, double sx, double sy) { + VectorXd sv(length); + int numRows = concentrations.rows(); + BC_TYPE type; - VectorXd sv(length); - int numRows = concentrations.rows(); - BC_TYPE type; - - // inner rows - if (rowIndex > 0 && rowIndex < numRows-1) { - for (int i = 0; i < length; i++) { - sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) - * concentrations(rowIndex+1,i) - + ( - 1 - sy * ( - calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) - + calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i)) - ) - ) * concentrations(rowIndex,i) - + sy * calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i)) - * concentrations(rowIndex-1,i) - ; - } + // inner rows + if (rowIndex > 0 && rowIndex < numRows - 1) { + for (int i = 0; i < length; i++) { + sv(i) = + sy * + calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) * + concentrations(rowIndex + 1, i) + + (1 - sy * (calcAlphaIntercell(alphaY(rowIndex, i), + alphaY(rowIndex + 1, i)) + + calcAlphaIntercell(alphaY(rowIndex - 1, i), + alphaY(rowIndex, i)))) * + concentrations(rowIndex, i) + + sy * + calcAlphaIntercell(alphaY(rowIndex - 1, i), alphaY(rowIndex, i)) * + concentrations(rowIndex - 1, i); } + } - // first row - else if (rowIndex == 0) { - for (int i = 0; i < length; i++) { - type = bcTop[i].getType(); - if (type == BC_TYPE_CONSTANT) { - sv(i) = calcExplicitConcentrationsTopBoundaryConstant(concentrations, alphaY, bcTop, rowIndex, i, sy); - } else if (type == BC_TYPE_CLOSED) { - sv(i) = calcExplicitConcentrationsTopBoundaryClosed(concentrations, alphaY, rowIndex, i, sy); - } else { - throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!"); - } - } + // first row + else if (rowIndex == 0) { + for (int i = 0; i < length; i++) { + type = bcTop[i].getType(); + if (type == BC_TYPE_CONSTANT) { + sv(i) = calcExplicitConcentrationsTopBoundaryConstant( + concentrations, alphaY, bcTop, rowIndex, i, sy); + } else if (type == BC_TYPE_CLOSED) { + sv(i) = calcExplicitConcentrationsTopBoundaryClosed( + concentrations, alphaY, rowIndex, i, sy); + } else { + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Left or Top!"); + } } + } - // last row - else if (rowIndex == numRows-1) { - for (int i = 0; i < length; i++) { - type = bcBottom[i].getType(); - if (type == BC_TYPE_CONSTANT) { - sv(i) = calcExplicitConcentrationsBottomBoundaryConstant(concentrations, alphaY, bcBottom, rowIndex, i, sy); - } else if (type == BC_TYPE_CLOSED) { - sv(i) = calcExplicitConcentrationsBottomBoundaryClosed(concentrations, alphaY, rowIndex, i, sy); - } else { - throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!"); - } - } + // last row + else if (rowIndex == numRows - 1) { + for (int i = 0; i < length; i++) { + type = bcBottom[i].getType(); + if (type == BC_TYPE_CONSTANT) { + sv(i) = calcExplicitConcentrationsBottomBoundaryConstant( + concentrations, alphaY, bcBottom, rowIndex, i, sy); + } else if (type == BC_TYPE_CLOSED) { + sv(i) = calcExplicitConcentrationsBottomBoundaryClosed( + concentrations, alphaY, rowIndex, i, sy); + } else { + throw_invalid_argument( + "Undefined Boundary Condition Type somewhere on Right or Bottom!"); + } } + } - // first column -> additional fixed concentration change from perpendicular dimension in constant bc case - if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) { - sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft[rowIndex].getValue(); - } + // first column -> additional fixed concentration change from perpendicular + // dimension in constant bc case + if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) { + sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue(); + } - // last column -> additional fixed concentration change from perpendicular dimension in constant bc case - if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) { - sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight[rowIndex].getValue(); - } + // last column -> additional fixed concentration change from perpendicular + // dimension in constant bc case + if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) { + sv(length - 1) += + 2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue(); + } - return sv; + return sv; } - // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // use of EigenLU solver static VectorXd EigenLUAlgorithm(SparseMatrix &A, VectorXd &b) { - SparseLU> solver; - solver.analyzePattern(A); - solver.factorize(A); + SparseLU> solver; + solver.analyzePattern(A); + solver.factorize(A); - return solver.solve(b); + return solver.solve(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 // implementation of Thomas Algorithm static VectorXd ThomasAlgorithm(SparseMatrix &A, VectorXd &b) { - uint32_t n = b.size(); + uint32_t n = b.size(); - Eigen::VectorXd a_diag(n); - Eigen::VectorXd b_diag(n); - Eigen::VectorXd c_diag(n); - Eigen::VectorXd x_vec = b; + Eigen::VectorXd a_diag(n); + Eigen::VectorXd b_diag(n); + Eigen::VectorXd c_diag(n); + Eigen::VectorXd x_vec = b; - // Fill diagonals vectors - b_diag[0] = A.coeff(0, 0); - c_diag[0] = A.coeff(0, 1); + // Fill diagonals vectors + b_diag[0] = A.coeff(0, 0); + c_diag[0] = A.coeff(0, 1); - for (int i = 1; i < n - 1; i++) { - a_diag[i] = A.coeff(i, i - 1); - b_diag[i] = A.coeff(i, i); - c_diag[i] = A.coeff(i, i + 1); - } - a_diag[n - 1] = A.coeff(n - 1, n - 2); - b_diag[n - 1] = A.coeff(n - 1, n - 1); + for (int i = 1; i < n - 1; i++) { + a_diag[i] = A.coeff(i, i - 1); + b_diag[i] = A.coeff(i, i); + c_diag[i] = A.coeff(i, i + 1); + } + a_diag[n - 1] = A.coeff(n - 1, n - 2); + b_diag[n - 1] = A.coeff(n - 1, n - 1); - // HACK: write CSV to file + // HACK: write CSV to file #ifdef WRITE_THOMAS_CSV - #include - #include - static std::uint32_t file_index = 0; - std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; +#include +#include + static std::uint32_t file_index = 0; + std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; - std::ofstream out_file; + std::ofstream out_file; - out_file.open(file_name, std::ofstream::trunc | std::ofstream::out); + out_file.open(file_name, std::ofstream::trunc | std::ofstream::out); - // print header - out_file << "Aa, Ab, Ac, b\n"; + // print header + out_file << "Aa, Ab, Ac, b\n"; - // iterate through all elements - for (std::size_t i = 0; i < n; i++) { - out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " << b[i] << "\n"; - } + // iterate through all elements + for (std::size_t i = 0; i < n; i++) { + out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " + << b[i] << "\n"; + } - out_file.close(); + out_file.close(); #endif - // start solving - c_diag and x_vec are overwritten - n--; - c_diag[0] /= b_diag[0]; - x_vec[0] /= b_diag[0]; + // start solving - c_diag and x_vec are overwritten + n--; + c_diag[0] /= b_diag[0]; + x_vec[0] /= b_diag[0]; - for (int i = 1; i < n; i++) { - c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1]; - x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) / - (b_diag[i] - a_diag[i] * c_diag[i - 1]); - } + for (int i = 1; i < n; i++) { + c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1]; + x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) / + (b_diag[i] - a_diag[i] * c_diag[i - 1]); + } - x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) / - (b_diag[n] - a_diag[n] * c_diag[n - 1]); + x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) / + (b_diag[n] - a_diag[n] * c_diag[n - 1]); - for (int i = n; i-- > 0;) { - x_vec[i] -= c_diag[i] * x_vec[i + 1]; - } + for (int i = n; i-- > 0;) { + x_vec[i] -= c_diag[i] * x_vec[i + 1]; + } - return x_vec; + return x_vec; } +// BTCS solution for 1D grid +static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, + VectorXd (*solverFunc)(SparseMatrix &A, + VectorXd &b)) { + int length = grid.getLength(); + double sx = timestep / (grid.getDelta() * grid.getDelta()); -// BTCS solution for 1D grid -static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix &A, VectorXd &b)) { - int length = grid.getLength(); - double sx = timestep / (grid.getDelta() * grid.getDelta()); + VectorXd concentrations_t1(length); - VectorXd concentrations_t1(length); + SparseMatrix A; + VectorXd b(length); - SparseMatrix A; - VectorXd b(length); + MatrixXd alpha = grid.getAlpha(); + vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - MatrixXd alpha = grid.getAlpha(); - vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + MatrixXd concentrations = grid.getConcentrations(); + int rowIndex = 0; + A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, + sx); // this is exactly same as in 2D + for (int i = 0; i < length; i++) { + b(i) = concentrations(0, i); + } + if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) { + b(0) += 2 * sx * alpha(0, 0) * bcLeft[0].getValue(); + } + if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) { + b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue(); + } - MatrixXd concentrations = grid.getConcentrations(); - int rowIndex = 0; - A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D - for (int i = 0; i < length; i++) { - b(i) = concentrations(0,i); - } - if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) { - b(0) += 2 * sx * alpha(0,0) * bcLeft[0].getValue(); - } - if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) { - b(length-1) += 2 * sx * alpha(0,length-1) * bcRight[0].getValue(); - } + concentrations_t1 = solverFunc(A, b); - concentrations_t1 = solverFunc(A, b); + for (int j = 0; j < length; j++) { + concentrations(0, j) = concentrations_t1(j); + } - for (int j = 0; j < length; j++) { - concentrations(0,j) = concentrations_t1(j); - } - - grid.setConcentrations(concentrations); + grid.setConcentrations(concentrations); } - // BTCS solution for 2D grid -static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix &A, VectorXd &b), int numThreads) { - int rowMax = grid.getRow(); - int colMax = grid.getCol(); - double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); - double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); +static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, + VectorXd (*solverFunc)(SparseMatrix &A, + VectorXd &b), + int numThreads) { + int rowMax = grid.getRow(); + int colMax = grid.getCol(); + double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); + double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); - MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); - VectorXd row_t1(colMax); + MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0); + VectorXd row_t1(colMax); - SparseMatrix A; - VectorXd b; + SparseMatrix A; + VectorXd b; - MatrixXd alphaX = grid.getAlphaX(); - MatrixXd alphaY = grid.getAlphaY(); - vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); - vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + MatrixXd alphaX = grid.getAlphaX(); + MatrixXd alphaY = grid.getAlphaY(); + vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); + vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); - MatrixXd concentrations = grid.getConcentrations(); + MatrixXd concentrations = grid.getConcentrations(); - #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) - for (int i = 0; i < rowMax; i++) { - - - A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); - b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, - bcTop, bcBottom, colMax, i, sx, sy); - - SparseLU> solver; +#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) + for (int i = 0; i < rowMax; i++) { - row_t1 = solverFunc(A, b); - - concentrations_t1.row(i) = row_t1; - } - - concentrations_t1.transposeInPlace(); - concentrations.transposeInPlace(); - alphaX.transposeInPlace(); - alphaY.transposeInPlace(); - - #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) - for (int i = 0; i < colMax; i++) { - // swap alphas, boundary conditions and sx/sy for column-wise calculation - A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); - b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, - bcLeft, bcRight, rowMax, i, sy, sx); - - row_t1 = solverFunc(A, b); + A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); + b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, + bcTop, bcBottom, colMax, i, sx, sy); - concentrations.row(i) = row_t1; - } - - concentrations.transposeInPlace(); + SparseLU> solver; - grid.setConcentrations(concentrations); + row_t1 = solverFunc(A, b); + + concentrations_t1.row(i) = row_t1; + } + + concentrations_t1.transposeInPlace(); + concentrations.transposeInPlace(); + alphaX.transposeInPlace(); + alphaY.transposeInPlace(); + +#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) + for (int i = 0; i < colMax; i++) { + // swap alphas, boundary conditions and sx/sy for column-wise calculation + A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); + b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, + bcLeft, bcRight, rowMax, i, sy, sx); + + row_t1 = solverFunc(A, b); + + concentrations.row(i) = row_t1; + } + + concentrations.transposeInPlace(); + + grid.setConcentrations(concentrations); } - // entry point for EigenLU solver; differentiate between 1D and 2D grid static void BTCS_LU(Grid &grid, Boundary &bc, double 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); - } else { - throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); - } + if (grid.getDim() == 1) { + BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); + } else if (grid.getDim() == 2) { + BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads); + } else { + throw_invalid_argument( + "Error: Only 1- and 2-dimensional grids are defined!"); + } } // entry point for Thomas algorithm solver; differentiate 1D and 2D grid -static void BTCS_Thomas(Grid &grid, Boundary &bc, double 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); - } else { - throw_invalid_argument("Error: Only 1- and 2-dimensional grids are defined!"); - } -} \ No newline at end of file +static void BTCS_Thomas(Grid &grid, Boundary &bc, double 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); + } else { + throw_invalid_argument( + "Error: Only 1- and 2-dimensional grids are defined!"); + } +} From fab0f35ed0e401f919534a4ac8776266efaf31b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 6 Sep 2023 13:20:40 +0200 Subject: [PATCH 03/17] Add benchmark --- CMakeLists.txt | 6 ++ naaice/BTCS_2D_NAAICE.cpp | 150 ++++++++++++++++++++++++++++++++++++++ naaice/CMakeLists.txt | 9 +++ naaice/alphax.csv | 5 ++ naaice/files.hpp.in | 7 ++ naaice/init_conc.csv | 5 ++ src/CMakeLists.txt | 4 +- 7 files changed, 183 insertions(+), 3 deletions(-) create mode 100644 naaice/BTCS_2D_NAAICE.cpp create mode 100644 naaice/CMakeLists.txt create mode 100644 naaice/alphax.csv create mode 100644 naaice/files.hpp.in create mode 100644 naaice/init_conc.csv diff --git a/CMakeLists.txt b/CMakeLists.txt index a76e995..ec9da4c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,3 +42,9 @@ if(TUG_ENABLE_TESTING) endif() add_subdirectory(examples) + +option(TUG_NAAICE_EXAMPLE "Enables NAAICE examples with export of CSV files" OFF) + +if (TUG_NAAICE_EXAMPLE) + add_subdirectory(naaice) +endif() diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp new file mode 100644 index 0000000..1b99849 --- /dev/null +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -0,0 +1,150 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "files.hpp" + +template inline T parseString(const std::string &str) { + T result; + std::istringstream iss(str); + + if (!(iss >> result)) { + throw std::invalid_argument("Invalid input for parsing."); + } + + return result; +} + +template +std::vector tokenize(const std::string &input, char delimiter) { + std::vector tokens; + std::istringstream tokenStream(input); + std::string token; + + while (std::getline(tokenStream, token, delimiter)) { + tokens.push_back(parseString(token)); + } + + return tokens; +} + +template +std::vector> CSVToVector(const char *filename) { + std::ifstream in_file(filename); + if (!in_file.is_open()) { + throw std::runtime_error("Error opening file \'" + std::string(filename) + + "\'."); + } + + std::vector> csv_data; + + std::string line; + while (std::getline(in_file, line)) { + csv_data.push_back(tokenize(line, ',')); + } + + in_file.close(); + + return csv_data; +} + +template +Eigen::MatrixXd CMVecToRMMatrix(const std::vector> &vec, + std::uint32_t exp_rows, + std::uint32_t exp_cols) { + if (exp_rows != vec.size()) { + throw std::runtime_error( + "Mismatch in y dimension while converting to Eigen::Matrix."); + } + + Eigen::MatrixXd out_mat(exp_rows, exp_cols); + + for (std::uint32_t ri = 0; ri < exp_rows; ri++) { + const auto &vec_row = vec[ri]; + if (vec[ri].size() != exp_cols) { + throw std::runtime_error( + "Mismatch in x dimension while converting to Eigen::Matrix."); + } + for (std::uint32_t cj = 0; cj < exp_cols; cj++) { + out_mat(ri, cj) = vec_row[cj]; + } + } + + return out_mat; +} + +int main(int argc, char *argv[]) { + // EASY_PROFILER_ENABLE; + // profiler::startListen(); + // ************** + // **** GRID **** + // ************** + // profiler::startListen(); + // create a grid with a 5 x 10 field + constexpr int row = 5; + constexpr int col = 10; + Grid grid = Grid(row, col); + + // (optional) set the domain, e.g.: + grid.setDomain(0.005, 0.01); + + const auto init_values_vec = CSVToVector(INPUT_CONC_FILE); + MatrixXd concentrations = CMVecToRMMatrix(init_values_vec, row, col); + grid.setConcentrations(concentrations); + + // // (optional) set alphas of the grid, e.g.: + const auto alphax_vec = CSVToVector(INPUT_ALPHAX_FILE); + MatrixXd alphax = CMVecToRMMatrix(alphax_vec, row, col); + + constexpr double alphay_val = 5e-10; + MatrixXd alphay = MatrixXd::Constant(row, col, alphay_val); // row,col,value + grid.setAlpha(alphax, alphay); + + // // ****************** + // // **** BOUNDARY **** + // // ****************** + + // create a boundary with constant values + Boundary bc = Boundary(grid); + bc.setBoundarySideClosed(BC_SIDE_LEFT); + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); + + // // ************************ + // // **** SIMULATION ENV **** + // // ************************ + + // set up a simulation environment + Simulation simulation = Simulation( + grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + + // set the timestep of the simulation + simulation.setTimestep(360); // timestep + + // set the number of iterations + simulation.setIterations(1); + + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + simulation.setOutputCSV(CSV_OUTPUT_ON); + + simulation.setOutputConsole(CONSOLE_OUTPUT_ON); + + // // **** RUN SIMULATION **** + + // // run the simulation + + // // EASY_BLOCK("SIMULATION") + simulation.run(); + // // EASY_END_BLOCK; + // // profiler::dumpBlocksToFile("test_profile.prof"); + // // profiler::stopListen(); +} diff --git a/naaice/CMakeLists.txt b/naaice/CMakeLists.txt new file mode 100644 index 0000000..947e720 --- /dev/null +++ b/naaice/CMakeLists.txt @@ -0,0 +1,9 @@ +add_executable(naaice BTCS_2D_NAAICE.cpp) + +get_filename_component(IN_CONC_FILE "init_conc.csv" REALPATH) +get_filename_component(IN_ALPHAX_FILE "alphax.csv" REALPATH) + +configure_file(files.hpp.in files.hpp) +target_include_directories(naaice PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) + +target_link_libraries(naaice PUBLIC tug) diff --git a/naaice/alphax.csv b/naaice/alphax.csv new file mode 100644 index 0000000..17db4c6 --- /dev/null +++ b/naaice/alphax.csv @@ -0,0 +1,5 @@ +9.35863547772169e-10,4.1475358910393e-10,7.40997499064542e-10,4.63898729439825e-10,5.50232032872736e-10,1.83670640387572e-10,5.84096802584827e-10,4.72976535512134e-10,2.92526249657385e-10,5.8247926668264e-10 +8.08492869278416e-10,3.5943602763582e-10,6.87254608655348e-10,1.52116895909421e-10,5.95404988038354e-10,8.90064929216169e-10,6.13143724552356e-10,1.23507722397335e-10,2.898759260308e-10,9.90528965741396e-10 +3.12359050870873e-10,6.34084914484993e-10,8.28234328492545e-10,1.60584925650619e-10,1.31777092232369e-10,7.34155903453939e-10,9.86383097898215e-10,5.42474469379522e-10,1.23030534153804e-10,2.33146838657558e-10 +7.76231796317734e-10,1.69479642040096e-10,6.74477553134784e-10,2.4903867142275e-10,2.70820381003432e-10,1.67315319390036e-10,7.1631961625535e-10,4.51548034301959e-10,2.41610987577587e-10,4.98075650655665e-10 +7.33895266707987e-10,5.82150935800746e-10,1.49088777275756e-10,5.00345975253731e-10,8.26257051993161e-10,5.28838745504618e-10,9.94136832957156e-10,7.44971914449707e-10,7.53557282453403e-10,7.54089470859617e-10 diff --git a/naaice/files.hpp.in b/naaice/files.hpp.in new file mode 100644 index 0000000..c7dd813 --- /dev/null +++ b/naaice/files.hpp.in @@ -0,0 +1,7 @@ +#ifndef FILES_H_ +#define FILES_H_ + +const char *INPUT_CONC_FILE = "@IN_CONC_FILE@"; +const char *INPUT_ALPHAX_FILE = "@IN_ALPHAX_FILE@"; + +#endif // FILES_H_ diff --git a/naaice/init_conc.csv b/naaice/init_conc.csv new file mode 100644 index 0000000..ad55817 --- /dev/null +++ b/naaice/init_conc.csv @@ -0,0 +1,5 @@ +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 +6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,6.92023e-07,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08,2.02396e-08 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3cc79d4..04db461 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,8 +1,6 @@ add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCSv2.cpp) -option(TUG_WRITE_CSV "Write CSV during Thomas algorithm with consecutive numbers for each call" OFF) - -IF(TUG_WRITE_CSV) +IF(TUG_NAAICE_EXAMPLE) target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV) endif() From 72107c944d452eb284971234dbdf524c1f789c85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 6 Sep 2023 14:21:22 +0200 Subject: [PATCH 04/17] Add readme --- naaice/README.org | 72 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 naaice/README.org diff --git a/naaice/README.org b/naaice/README.org new file mode 100644 index 0000000..658eb12 --- /dev/null +++ b/naaice/README.org @@ -0,0 +1,72 @@ +#+title: NAAICE Benchmark + +This directory contains a concise benchmark designed for validating FPGA +offloading of the Thomas algorithm, primarily employed for solving linear +equation systems structured within a tridiagonal matrix. + +* Benchmark Setup + +The benchmark involves a domain measuring $0.5 \text{cm} \times 1 \text{cm}$, +divided into a grid of dimensions $10 \times 5$. Each grid cell initially +contains a specific concentration. The concentration in the first half along the +x-dimension is set at $6.92023 \times 10^{-7}$, while in the second half, it's +$2.02396 \times 10^{-8}$, creating a concentration gradient along the y-axis at +the center of the grid. + +To achieve concentration equilibrium, we employ a simulation based on a +heterogeneous 2D-ADI BTCS diffusion approach, detailed in the [[file:../doc/ADI_scheme.pdf][ADI_scheme.pdf]] +file. In the x-direction, diffusion coefficients range from $\alpha = 10^{-9}$ +to $10^{-10}$, while in the y-direction, a constant value of $5 \times 10^{-10}$ +is applied. A closed boundary condition is implemented, meaning concentrations +cannot enter or exit the system. The diffusion process is simulated for a single +iteration with a time step ($\Delta t$) of 360 seconds. + +* Setup + +To generate new makefiles using the =-DTUG_NAAICE_EXAMPLE=ON= option in CMake, +compile the executable, and run it to generate the benchmark output, follow +these steps: + +1. Navigate to your project's build directory. +2. Run the following CMake command with the =-DTUG_NAAICE_EXAMPLE=ON= option to + generate the makefiles: + + #+BEGIN_SRC sh + cmake -DTUG_NAAICE_EXAMPLE=ON .. + #+END_SRC + +3. After CMake configuration is complete, build the =naaice= executable by running =make=: + + #+BEGIN_SRC sh + make naaice + #+END_SRC + +4. Once the compilation is successful, navigate to the build directory by =cd + /naaice= + +5. Finally, run the =naaice= executable to generate the benchmark output: + + #+BEGIN_SRC sh + ./naaice + #+END_SRC + +** Output Files + +*** =Thomas_.csv= + +These files contain the values of the tridiagonal coefficient matrix $A$, where: + +- $Aa$ represents the leftmost value, +- $Ab$ represents the middle value, and +- $Ac$ represents the rightmost value of one row of the matrix. + +Additionally, the corresponding values of the right-hand-side vector $b$ are +provided. + +Since the 2D-ADI BTCS scheme processes each row first and then proceeds +column-wise through the grid, each iteration is saved separately in +consecutively numbered files. + +*** =BTCS_5_10_1.csv= + +The result of the simulation, *separated by whitespaces*! From f8bdfe39ea2b50cb9bb3bb4e81c9e11b2144cf26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 6 Sep 2023 14:23:45 +0200 Subject: [PATCH 05/17] Use Markdown over Org-Mode --- naaice/README.md | 72 +++++++++++++++++++++++++++++++++++++++++++++++ naaice/README.org | 72 ----------------------------------------------- 2 files changed, 72 insertions(+), 72 deletions(-) create mode 100644 naaice/README.md delete mode 100644 naaice/README.org diff --git a/naaice/README.md b/naaice/README.md new file mode 100644 index 0000000..eff0ede --- /dev/null +++ b/naaice/README.md @@ -0,0 +1,72 @@ + +This directory contains a concise benchmark designed for validating FPGA +offloading of the Thomas algorithm, primarily employed for solving linear +equation systems structured within a tridiagonal matrix. + + +# Benchmark Setup + +The benchmark involves a domain measuring $0.5 \text{cm} \times 1 \text{cm}$, +divided into a grid of dimensions $10 \times 5$. Each grid cell initially +contains a specific concentration. The concentration in the first half along the +x-dimension is set at $6.92023 \times 10^{-7}$, while in the second half, it’s +$2.02396 \times 10^{-8}$, creating a concentration gradient along the y-axis at +the center of the grid. + +To achieve concentration equilibrium, we employ a simulation based on a +heterogeneous 2D-ADI BTCS diffusion approach, detailed in the +[ADIscheme.pdf](../doc/ADI_scheme.pdf) file. In the x-direction, +diffusion coefficients range from $\alpha = 10^{-9}$ to $10^{-10}$, while in the +y-direction, a constant value of $5 \times 10^{-10}$ is applied. A closed +boundary condition is implemented, meaning concentrations cannot enter or exit +the system. The diffusion process is simulated for a single iteration with a +time step ($\Delta t$) of 360 seconds. + + +# Setup + +To generate new makefiles using the `-DTUG_NAAICE_EXAMPLE=ON` option in CMake, +compile the executable, and run it to generate the benchmark output, follow +these steps: + +1. Navigate to your project's build directory. +2. Run the following CMake command with the `-DTUG_NAAICE_EXAMPLE=ON` option to + generate the makefiles: + + cmake -DTUG_NAAICE_EXAMPLE=ON .. + +3. After CMake configuration is complete, build the `naaice` executable by running `make`: + + make naaice + +4. Once the compilation is successful, navigate to the build directory by `cd + /naaice` + +5. Finally, run the `naaice` executable to generate the benchmark output: + + ./naaice + + +## Output Files + + +### `Thomas_.csv` + +These files contain the values of the tridiagonal coefficient matrix $A$, where: + +- $Aa$ represents the leftmost value, +- $Ab$ represents the middle value, and +- $Ac$ represents the rightmost value of one row of the matrix. + +Additionally, the corresponding values of the right-hand-side vector $b$ are +provided. + +Since the 2D-ADI BTCS scheme processes each row first and then proceeds +column-wise through the grid, each iteration is saved separately in +consecutively numbered files. + + +### `BTCS_5_10_1.csv` + +The result of the simulation, **separated by whitespaces**! + diff --git a/naaice/README.org b/naaice/README.org deleted file mode 100644 index 658eb12..0000000 --- a/naaice/README.org +++ /dev/null @@ -1,72 +0,0 @@ -#+title: NAAICE Benchmark - -This directory contains a concise benchmark designed for validating FPGA -offloading of the Thomas algorithm, primarily employed for solving linear -equation systems structured within a tridiagonal matrix. - -* Benchmark Setup - -The benchmark involves a domain measuring $0.5 \text{cm} \times 1 \text{cm}$, -divided into a grid of dimensions $10 \times 5$. Each grid cell initially -contains a specific concentration. The concentration in the first half along the -x-dimension is set at $6.92023 \times 10^{-7}$, while in the second half, it's -$2.02396 \times 10^{-8}$, creating a concentration gradient along the y-axis at -the center of the grid. - -To achieve concentration equilibrium, we employ a simulation based on a -heterogeneous 2D-ADI BTCS diffusion approach, detailed in the [[file:../doc/ADI_scheme.pdf][ADI_scheme.pdf]] -file. In the x-direction, diffusion coefficients range from $\alpha = 10^{-9}$ -to $10^{-10}$, while in the y-direction, a constant value of $5 \times 10^{-10}$ -is applied. A closed boundary condition is implemented, meaning concentrations -cannot enter or exit the system. The diffusion process is simulated for a single -iteration with a time step ($\Delta t$) of 360 seconds. - -* Setup - -To generate new makefiles using the =-DTUG_NAAICE_EXAMPLE=ON= option in CMake, -compile the executable, and run it to generate the benchmark output, follow -these steps: - -1. Navigate to your project's build directory. -2. Run the following CMake command with the =-DTUG_NAAICE_EXAMPLE=ON= option to - generate the makefiles: - - #+BEGIN_SRC sh - cmake -DTUG_NAAICE_EXAMPLE=ON .. - #+END_SRC - -3. After CMake configuration is complete, build the =naaice= executable by running =make=: - - #+BEGIN_SRC sh - make naaice - #+END_SRC - -4. Once the compilation is successful, navigate to the build directory by =cd - /naaice= - -5. Finally, run the =naaice= executable to generate the benchmark output: - - #+BEGIN_SRC sh - ./naaice - #+END_SRC - -** Output Files - -*** =Thomas_.csv= - -These files contain the values of the tridiagonal coefficient matrix $A$, where: - -- $Aa$ represents the leftmost value, -- $Ab$ represents the middle value, and -- $Ac$ represents the rightmost value of one row of the matrix. - -Additionally, the corresponding values of the right-hand-side vector $b$ are -provided. - -Since the 2D-ADI BTCS scheme processes each row first and then proceeds -column-wise through the grid, each iteration is saved separately in -consecutively numbered files. - -*** =BTCS_5_10_1.csv= - -The result of the simulation, *separated by whitespaces*! From 40710a0b39a37a325495b14da14c3736a17eba1d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 6 Sep 2023 14:25:02 +0200 Subject: [PATCH 06/17] Usage instead of setup --- naaice/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/naaice/README.md b/naaice/README.md index eff0ede..54a2c95 100644 --- a/naaice/README.md +++ b/naaice/README.md @@ -23,7 +23,7 @@ the system. The diffusion process is simulated for a single iteration with a time step ($\Delta t$) of 360 seconds. -# Setup +# Usage To generate new makefiles using the `-DTUG_NAAICE_EXAMPLE=ON` option in CMake, compile the executable, and run it to generate the benchmark output, follow From bf4444fc8419385bdc5d6c0d654ecf1de421d183 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 6 Sep 2023 14:39:27 +0200 Subject: [PATCH 07/17] add comments and rename function --- naaice/BTCS_2D_NAAICE.cpp | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index 1b99849..a3937b4 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -11,6 +11,9 @@ #include "files.hpp" +/** + * Try to parse an input string into a given template type. + */ template inline T parseString(const std::string &str) { T result; std::istringstream iss(str); @@ -22,6 +25,9 @@ template inline T parseString(const std::string &str) { return result; } +/** + * Splits a given string into a vector by using a delimiter character. + */ template std::vector tokenize(const std::string &input, char delimiter) { std::vector tokens; @@ -35,6 +41,9 @@ std::vector tokenize(const std::string &input, char delimiter) { return tokens; } +/** + * Opens a file containing CSV and transform it into row-major 2D STL vector. + */ template std::vector> CSVToVector(const char *filename) { std::ifstream in_file(filename); @@ -55,8 +64,12 @@ std::vector> CSVToVector(const char *filename) { return csv_data; } +/** + * Converts a 2D STL vector, where values are stored row-major into a + * column-major Eigen::Matrix. + */ template -Eigen::MatrixXd CMVecToRMMatrix(const std::vector> &vec, +Eigen::MatrixXd rmVecTocmMatrix(const std::vector> &vec, std::uint32_t exp_rows, std::uint32_t exp_cols) { if (exp_rows != vec.size()) { @@ -96,12 +109,12 @@ int main(int argc, char *argv[]) { grid.setDomain(0.005, 0.01); const auto init_values_vec = CSVToVector(INPUT_CONC_FILE); - MatrixXd concentrations = CMVecToRMMatrix(init_values_vec, row, col); + MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); grid.setConcentrations(concentrations); // // (optional) set alphas of the grid, e.g.: const auto alphax_vec = CSVToVector(INPUT_ALPHAX_FILE); - MatrixXd alphax = CMVecToRMMatrix(alphax_vec, row, col); + MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); constexpr double alphay_val = 5e-10; MatrixXd alphay = MatrixXd::Constant(row, col, alphay_val); // row,col,value @@ -123,8 +136,8 @@ int main(int argc, char *argv[]) { // // ************************ // set up a simulation environment - Simulation simulation = Simulation( - grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + Simulation simulation = + Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach // set the timestep of the simulation simulation.setTimestep(360); // timestep @@ -136,15 +149,9 @@ int main(int argc, char *argv[]) { // CSV_OUTPUT_VERBOSE] simulation.setOutputCSV(CSV_OUTPUT_ON); + // set output to the console to 'ON' simulation.setOutputConsole(CONSOLE_OUTPUT_ON); // // **** RUN SIMULATION **** - - // // run the simulation - - // // EASY_BLOCK("SIMULATION") simulation.run(); - // // EASY_END_BLOCK; - // // profiler::dumpBlocksToFile("test_profile.prof"); - // // profiler::stopListen(); } From 42ad07c252d59803fc9176c0e61203a65a7d5358 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Thu, 7 Sep 2023 09:24:03 +0200 Subject: [PATCH 08/17] Fix: some updates to naaice/README.md --- naaice/README.md | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/naaice/README.md b/naaice/README.md index 54a2c95..b4f1f9d 100644 --- a/naaice/README.md +++ b/naaice/README.md @@ -1,4 +1,3 @@ - This directory contains a concise benchmark designed for validating FPGA offloading of the Thomas algorithm, primarily employed for solving linear equation systems structured within a tridiagonal matrix. @@ -6,20 +5,19 @@ equation systems structured within a tridiagonal matrix. # Benchmark Setup -The benchmark involves a domain measuring $0.5 \text{cm} \times 1 \text{cm}$, -divided into a grid of dimensions $10 \times 5$. Each grid cell initially -contains a specific concentration. The concentration in the first half along the -x-dimension is set at $6.92023 \times 10^{-7}$, while in the second half, it’s -$2.02396 \times 10^{-8}$, creating a concentration gradient along the y-axis at -the center of the grid. +The benchmark defines a domain measuring $1 \text{cm} \times 0.5 \text{cm}$ (easting $\times$ northing), +discretized in a $10 \times 5$ grid. Each grid cell initially +contains a specific concentration. The concentration in the left domain half is set to $6.92023 \times 10^{-7}$, while in the right half to +$2.02396 \times 10^{-8}$, creating an horizontal concentration discontinuity at +the center of the grid. These initial concentrations are read from headerless csv file [init_conc.csv](./init_conc.csv). -To achieve concentration equilibrium, we employ a simulation based on a -heterogeneous 2D-ADI BTCS diffusion approach, detailed in the -[ADIscheme.pdf](../doc/ADI_scheme.pdf) file. In the x-direction, -diffusion coefficients range from $\alpha = 10^{-9}$ to $10^{-10}$, while in the -y-direction, a constant value of $5 \times 10^{-10}$ is applied. A closed -boundary condition is implemented, meaning concentrations cannot enter or exit -the system. The diffusion process is simulated for a single iteration with a +A diffusion time step is simulated with the +heterogeneous 2D-ADI approach detailed in the +[ADI_scheme.pdf](../doc/ADI_scheme.pdf) file. The x component of the +diffusion coefficients, read from headerless csv file [alphax.csv](./alphax.csv) ranges from $\alpha = 10^{-9}$ to $10^{-10}$ (distributed randomly), while the +y-component is held constant at $5 \times 10^{-10}$. Closed +boundary conditions are enforced at all domain boundaries, meaning that concentration cannot enter or exit +the system, or in other terms, that the sum of concentrations over the domain must stay constant. The benchmark simulates a single iteration with a time step ($\Delta t$) of 360 seconds. From 5490216bc084ded28a841751f617706d8fd49b31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 7 Sep 2023 11:10:56 +0200 Subject: [PATCH 09/17] feat: print sums of input and output field --- naaice/BTCS_2D_NAAICE.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index a3937b4..fe15455 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -112,6 +112,8 @@ int main(int argc, char *argv[]) { MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); grid.setConcentrations(concentrations); + const double sum_init = concentrations.sum(); + // // (optional) set alphas of the grid, e.g.: const auto alphax_vec = CSVToVector(INPUT_ALPHAX_FILE); MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); @@ -154,4 +156,12 @@ int main(int argc, char *argv[]) { // // **** RUN SIMULATION **** simulation.run(); + + const double sum_after = grid.getConcentrations().sum(); + + std::cout << "Sum of init field: " << std::to_string(sum_init) + << "\nSum after iteration: " << std::to_string(sum_after) + << std::endl; + + return 0; } From ba627b66245e37ba0681a9412c5f102eebcfe426 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 15 Sep 2023 10:39:51 +0200 Subject: [PATCH 10/17] feat: rewrite library as template library --- examples/BTCS_1D_proto_example.cpp | 2 +- examples/BTCS_2D_proto_example.cpp | 2 +- examples/CRNI_2D_proto_example.cpp | 2 +- examples/FTCS_1D_proto_example.cpp | 2 +- examples/FTCS_2D_proto_closed_mdl.cpp | 2 +- examples/FTCS_2D_proto_example.cpp | 2 +- examples/FTCS_2D_proto_example_mdl.cpp | 2 +- examples/profiling_openmp.cpp | 2 +- examples/profiling_speedup.cpp | 2 +- examples/reference-FTCS_2D_closed.cpp | 2 +- include/tug/Boundary.hpp | 172 +++++++++++++++--- include/tug/Grid.hpp | 236 ++++++++++++++++++++----- include/tug/Simulation.hpp | 145 +++++++++++++-- naaice/BTCS_2D_NAAICE.cpp | 8 +- src/BTCS.cpp | 132 ++++++++------ src/Boundary.cpp | 149 ---------------- src/CMakeLists.txt | 2 +- src/FTCS.cpp | 105 ++++++----- src/Grid.cpp | 167 ----------------- src/Schemes.hpp | 9 +- src/Simulation.cpp | 166 +++-------------- test/testBoundary.cpp | 11 +- test/testGrid.cpp | 20 +-- test/testSimulation.cpp | 8 +- 24 files changed, 664 insertions(+), 686 deletions(-) delete mode 100644 src/Boundary.cpp delete mode 100644 src/Grid.cpp diff --git a/examples/BTCS_1D_proto_example.cpp b/examples/BTCS_1D_proto_example.cpp index 1a53738..4bfd1b4 100644 --- a/examples/BTCS_1D_proto_example.cpp +++ b/examples/BTCS_1D_proto_example.cpp @@ -10,7 +10,7 @@ int main(int argc, char *argv[]) { // create a linear grid with 20 cells int cells = 20; - Grid grid = Grid(cells); + Grid64 grid(cells); MatrixXd concentrations = MatrixXd::Constant(1, 20, 0); concentrations(0, 0) = 2000; diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index 9a103cb..f14e0f5 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -13,7 +13,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 40; int col = 50; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/CRNI_2D_proto_example.cpp b/examples/CRNI_2D_proto_example.cpp index 39262e7..d605038 100644 --- a/examples/CRNI_2D_proto_example.cpp +++ b/examples/CRNI_2D_proto_example.cpp @@ -6,7 +6,7 @@ using namespace Eigen; int main(int argc, char *argv[]) { int row = 20; int col = 20; - Grid grid(row, col); + Grid64 grid(row, col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); concentrations(10, 10) = 2000; diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp index e97f0be..74b6f3e 100644 --- a/examples/FTCS_1D_proto_example.cpp +++ b/examples/FTCS_1D_proto_example.cpp @@ -11,7 +11,7 @@ int main(int argc, char *argv[]) { // create a linear grid with 20 cells int cells = 20; - Grid grid = Grid(cells); + Grid64 grid(cells); MatrixXd concentrations = MatrixXd::Constant(1, 20, 20); grid.setConcentrations(concentrations); diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 845d36f..e17f92b 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -30,7 +30,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int n2 = row / 2 - 1; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 5f5e520..89f821d 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -23,7 +23,7 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 20; int col = 20; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 2c63c0c..d77710d 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -21,7 +21,7 @@ int main(int argc, char *argv[]) { int row = 64; int col = 64; int n2 = row / 2 - 1; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); diff --git a/examples/profiling_openmp.cpp b/examples/profiling_openmp.cpp index 5dc13ef..c781900 100644 --- a/examples/profiling_openmp.cpp +++ b/examples/profiling_openmp.cpp @@ -29,7 +29,7 @@ int main(int argc, char *argv[]) { // myfile << "Iterations: " << iterations[j] << endl; for (int k = 0; k < repetition; k++) { cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); + Grid64 grid(n[i], n[i]); grid.setDomain(1, 1); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); diff --git a/examples/profiling_speedup.cpp b/examples/profiling_speedup.cpp index 5dc13ef..c781900 100644 --- a/examples/profiling_speedup.cpp +++ b/examples/profiling_speedup.cpp @@ -29,7 +29,7 @@ int main(int argc, char *argv[]) { // myfile << "Iterations: " << iterations[j] << endl; for (int k = 0; k < repetition; k++) { cout << "Wiederholung: " << k << endl; - Grid grid = Grid(n[i], n[i]); + Grid64 grid(n[i], n[i]); grid.setDomain(1, 1); MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp index 340a6fd..fa16063 100644 --- a/examples/reference-FTCS_2D_closed.cpp +++ b/examples/reference-FTCS_2D_closed.cpp @@ -12,7 +12,7 @@ int main(int argc, char *argv[]) { int domain_col = 10; // Grid - Grid grid = Grid(row, col); + Grid64 grid(row, col); grid.setDomain(domain_row, domain_col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index b803e75..b51df24 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -9,6 +9,8 @@ #include "Grid.hpp" #include +#include +#include /** * @brief Enum defining the two implemented boundary conditions. @@ -27,7 +29,7 @@ enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM }; * These can be flexibly used and combined later in other classes. * The class serves as an auxiliary class for structuring the Boundary class. */ -class BoundaryElement { +template class BoundaryElement { public: /** * @brief Construct a new Boundary Element object for the closed case. @@ -35,7 +37,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. @@ -45,7 +47,7 @@ public: * @param value Value of the constant concentration to be assumed at the * corresponding boundary element. */ - BoundaryElement(double value); + BoundaryElement(T _value) : value(_value), type(BC_TYPE_CONSTANT) {} /** * @brief Allows changing the boundary type of a corresponding @@ -54,7 +56,7 @@ public: * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or BC_TYPE_CLOSED. */ - void setType(BC_TYPE type); + void setType(BC_TYPE type) { this->type = type; }; /** * @brief Sets the value of a boundary condition for the constant case. @@ -62,7 +64,17 @@ public: * @param value Concentration to be considered constant for the * corresponding boundary element. */ - void setValue(double value); + void setValue(double value) { + if (value < 0) { + throw std::invalid_argument("No negative concentration allowed."); + } + if (type == BC_TYPE_CLOSED) { + throw std::invalid_argument( + "No constant boundary concentrations can be set for closed " + "boundaries. Please change type first."); + } + this->value = value; + } /** * @brief Return the type of the boundary condition, i.e. whether the @@ -71,25 +83,25 @@ public: * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or BC_TYPE_CONSTANT. */ - BC_TYPE getType(); + BC_TYPE getType() const { return this->type; } /** * @brief Return the concentration value for the constant boundary condition. * * @return double Value of the concentration. */ - double getValue(); + T getValue() const { return this->value; } private: BC_TYPE type{BC_TYPE_CLOSED}; - double value{-1}; + T value{-1}; }; /** * This class implements the functionality and management of the boundary * conditions in the grid to be simulated. */ -class Boundary { +template class Boundary { public: /** * @brief Creates a boundary object based on the passed grid object and @@ -98,7 +110,27 @@ public: * @param grid Grid object on the basis of which the simulation takes place * and from which the dimensions (in 2D case) are taken. */ - Boundary(Grid grid); + Boundary(const Grid &grid) + : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) { + if (this->dim == 1) { + this->boundaries = std::vector>>( + 2); // in 1D only left and right boundary + + this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); + } else if (this->dim == 2) { + this->boundaries = std::vector>>(4); + + this->boundaries[BC_SIDE_LEFT] = + std::vector>(this->rows, BoundaryElement()); + this->boundaries[BC_SIDE_RIGHT] = + std::vector>(this->rows, BoundaryElement()); + this->boundaries[BC_SIDE_TOP] = + std::vector>(this->cols, BoundaryElement()); + this->boundaries[BC_SIDE_BOTTOM] = + std::vector>(this->cols, BoundaryElement()); + } + } /** * @brief Sets all elements of the specified boundary side to the boundary @@ -106,7 +138,21 @@ public: * * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. */ - void setBoundarySideClosed(BC_SIDE side); + void setBoundarySideClosed(BC_SIDE side) { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? this->rows : this->cols; + + this->boundaries[side] = + std::vector>(n, BoundaryElement()); + } /** * @brief Sets all elements of the specified boundary side to the boundary @@ -117,7 +163,21 @@ public: * @param value Concentration to be set for all elements of the specified * page. */ - void setBoundarySideConstant(BC_SIDE side, double value); + void setBoundarySideConstant(BC_SIDE side, double value) { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional case, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + + const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; + const int n = is_vertical ? this->rows : this->cols; + + this->boundaries[side] = + std::vector>(n, BoundaryElement(value)); + } /** * @brief Specifically sets the boundary element of the specified side @@ -128,7 +188,13 @@ public: * boundary side. Must index an element of the corresponding * side. */ - void setBoundaryElementClosed(BC_SIDE side, int index); + void setBoundaryElemenClosed(BC_SIDE side, int index) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CLOSED); + } /** * @brief Specifically sets the boundary element of the specified side @@ -142,7 +208,14 @@ public: * @param value Concentration value to which the boundary element should be set. */ - void setBoundaryElementConstant(BC_SIDE side, int index, double value); + void setBoundaryElementConstant(BC_SIDE side, int index, double value) { + // tests whether the index really points to an element of the boundary side. + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + this->boundaries[side][index].setType(BC_TYPE_CONSTANT); + this->boundaries[side][index].setValue(value); + } /** * @brief Returns the boundary condition of a specified side as a vector @@ -150,19 +223,41 @@ public: * * @param side Boundary side from which the boundary conditions are to be * returned. - * @return vector Contains the boundary conditions as - * BoundaryElement objects. + * @return vector> Contains the boundary conditions as + * BoundaryElement objects. */ - const std::vector getBoundarySide(BC_SIDE side); + const std::vector> &getBoundarySide(BC_SIDE side) const { + if (this->dim == 1) { + if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { + throw std::invalid_argument( + "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "BC_SIDE_RIGHT borders exist."); + } + } + return this->boundaries[side]; + } /** * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some specific boundary is closed. * * @param side Boundary side for which the values are to be returned. - * @return VectorXd Vector with values as doubles. + * @return VectorX Vector with values as T. */ - Eigen::VectorXd getBoundarySideValues(BC_SIDE side); + Eigen::VectorX getBoundarySideValues(BC_SIDE side) const { + const std::size_t length = boundaries[side].size(); + Eigen::VectorX values(length); + + for (int i = 0; i < length; i++) { + if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { + values(i) = -1; + continue; + } + values(i) = getBoundaryElementValue(side, i); + } + + return values; + } /** * @brief Returns the boundary condition of a specified element on a given @@ -172,9 +267,15 @@ public: * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return BoundaryElement Boundary condition as a BoundaryElement object. + * @return BoundaryElement Boundary condition as a BoundaryElement + * object. */ - BoundaryElement getBoundaryElement(BC_SIDE side, int index); + BoundaryElement getBoundaryElement(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index]; + } /** * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED @@ -186,25 +287,42 @@ public: side. * @return BC_TYPE Boundary Type of the corresponding boundary condition. */ - BC_TYPE getBoundaryElementType(BC_SIDE side, int index); + BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + return this->boundaries[side][index].getType(); + } /** * @brief Returns the concentration value of a corresponding - * BoundaryElement object if it is a constant boundary condition. + * BoundaryElement object if it is a constant boundary condition. * * @param side Boundary side in which the boundary condition value is * located. * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return double Concentration of the corresponding BoundaryElement object. + * @return double Concentration of the corresponding BoundaryElement + * object. */ - double getBoundaryElementValue(BC_SIDE side, int index); + T getBoundaryElementValue(BC_SIDE side, int index) const { + if ((boundaries[side].size() < index) || index < 0) { + throw std::invalid_argument("Index is selected either too large or too small."); + } + if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { + throw std::invalid_argument( + "A value can only be output if it is a constant boundary condition."); + } + return this->boundaries[side][index].getValue(); + } private: - Grid grid; // Boundary is directly dependent on the dimensions of a predefined + const std::uint8_t dim; + const std::uint32_t cols; + const std::uint32_t rows; - std::vector> + std::vector>> boundaries; // Vector with Boundary Element information }; diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 6a6f973..3d3638d 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -10,8 +10,9 @@ #include #include +#include -class Grid { +template class Grid { public: /** * @brief Constructs a new 1D-Grid object of a given length, which holds a @@ -21,7 +22,18 @@ public: * * @param length Length of the 1D-Grid. Must be greater than 3. */ - Grid(int length); + Grid(int length) : col(length), domainCol(length) { + if (length <= 3) { + throw std::invalid_argument( + "Given grid length too small. Must be greater than 3."); + } + + this->dim = 1; + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + + this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); + } /** * @brief Constructs a new 2D-Grid object of given dimensions, which holds a @@ -34,146 +46,280 @@ public: * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. */ - Grid(int row, int col); + Grid(int _row, int _col) + : row(_row), col(_col), domainRow(_row), domainCol(_col) { + if (row <= 3 || col <= 3) { + throw std::invalid_argument( + "Given grid dimensions too small. Must each be greater than 3."); + } + + this->dim = 2; + this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 + this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + + this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + } /** * @brief Sets the concentrations matrix for a 1D or 2D-Grid. * - * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix - * must have correct dimensions as defined in row and col. (Or length, in 1D - * case). + * @param concentrations An Eigen3 MatrixX holding the concentrations. + * Matrix must have correct dimensions as defined in row and col. (Or length, + * in 1D case). */ - void setConcentrations(const Eigen::MatrixXd &concentrations); + void setConcentrations(const Eigen::MatrixX &concentrations) { + if (concentrations.rows() != this->row || + concentrations.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of concentrations mismatch with Grid dimensions!"); + } + + this->concentrations = concentrations; + } /** * @brief Gets the concentrations matrix for a Grid. * - * @return MatrixXd An Eigen3 matrix holding the concentrations and having the - * same dimensions as the grid. + * @return MatrixX An Eigen3 matrix holding the concentrations and having + * the same dimensions as the grid. */ - const Eigen::MatrixXd &getConcentrations() { return this->concentrations; } + const Eigen::MatrixX &getConcentrations() { return this->concentrations; } /** * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients. - * Matrix columns must have same size as length of grid. + * @param alpha An Eigen3 MatrixX with 1 row holding the alpha + * coefficients. Matrix columns must have same size as length of grid. */ - void setAlpha(const Eigen::MatrixXd &alpha); + void setAlpha(const Eigen::MatrixX &alpha) { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use 2D setter function!"); + } + if (alpha.rows() != 1 || alpha.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients mismatch with Grid dimensions!"); + } + + this->alphaX = alpha; + } /** * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two * dimensional. * - * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in + * @param alphaX An Eigen3 MatrixX holding the alpha coefficients in * x-direction. Matrix must be of same size as the grid. - * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in + * @param alphaY An Eigen3 MatrixX holding the alpha coefficients in * y-direction. Matrix must be of same size as the grid. */ - void setAlpha(const Eigen::MatrixXd &alphaX, const Eigen::MatrixXd &alphaY); + void setAlpha(const Eigen::MatrixX &alphaX, + const Eigen::MatrixX &alphaY) { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably " + "use 1D setter function!"); + } + if (alphaX.rows() != this->row || alphaX.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients in x-direction " + "mismatch with GRid dimensions!"); + } + if (alphaY.rows() != this->row || alphaY.cols() != this->col) { + throw std::invalid_argument( + "Given matrix of alpha coefficients in y-direction " + "mismatch with GRid dimensions!"); + } + + this->alphaX = alphaX; + this->alphaY = alphaY; + } /** * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @return MatrixXd A matrix with 1 row holding the alpha coefficients. + * @return MatrixX A matrix with 1 row holding the alpha coefficients. */ - const Eigen::MatrixXd &getAlpha(); + const Eigen::MatrixX &getAlpha() const { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use either getAlphaX() or getAlphaY()!"); + } + + return this->alphaX; + } /** * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixXd A matrix holding the alpha coefficients in x-direction. + * @return MatrixX A matrix holding the alpha coefficients in x-direction. */ - const Eigen::MatrixXd &getAlphaX(); + const Eigen::MatrixX &getAlphaX() const { + + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } + + return this->alphaX; + } /** * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixXd A matrix holding the alpha coefficients in y-direction. + * @return MatrixX A matrix holding the alpha coefficients in y-direction. */ - const Eigen::MatrixXd &getAlphaY(); + const Eigen::MatrixX &getAlphaY() const { + + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably use getAlpha()!"); + } + + return this->alphaY; + } /** * @brief Gets the dimensions of the grid. * * @return int Dimensions, either 1 or 2. */ - int getDim(); + int getDim() const { return this->dim; } /** * @brief Gets length of 1D grid. Must be one dimensional grid. * * @return int Length of 1D grid. */ - int getLength(); + int getLength() const { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use getRow() or getCol()!"); + } + + return col; + } /** * @brief Gets the number of rows of the grid. * * @return int Number of rows. */ - int getRow(); + int getRow() const { return this->row; } /** * @brief Gets the number of columns of the grid. * * @return int Number of columns. */ - int getCol(); + int getCol() const { return this->col; } /** * @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); + void setDomain(double domainLength) { + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probaly " + "use the 2D domain setter!"); + } + if (domainLength <= 0) { + throw std::invalid_argument("Given domain length is not positive!"); + } + + this->domainCol = domainLength; + this->deltaCol = double(this->domainCol) / double(this->col); + } /** * @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. + * @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); + void setDomain(double domainRow, double domainCol) { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, you should probably " + "use the 1D domain setter!"); + } + if (domainRow <= 0 || domainCol <= 0) { + throw std::invalid_argument("Given domain size is not positive!"); + } + + this->domainRow = domainRow; + this->domainCol = domainCol; + this->deltaRow = double(this->domainRow) / double(this->row); + this->deltaCol = double(this->domainCol) / double(this->col); + } /** * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. * * @return double Delta value. */ - double getDelta(); + T getDelta() const { + + if (dim != 1) { + throw std::invalid_argument( + "Grid is not one dimensional, you should probably " + "use the 2D delta getters"); + } + + return this->deltaCol; + } /** * @brief Gets the delta value in x-direction. * * @return double Delta value in x-direction. */ - double getDeltaCol(); + T getDeltaCol() const { return this->deltaCol; } /** * @brief Gets the delta value in y-direction. Must be two dimensional grid. * * @return double Delta value in y-direction. */ - double getDeltaRow(); + T getDeltaRow() const { + if (dim != 2) { + throw std::invalid_argument( + "Grid is not two dimensional, meaning there is no " + "delta in y-direction!"); + } + + return this->deltaRow; + } private: - int col; // number of grid columns - int row{1}; // number of grid rows - int dim; // 1D or 2D - double domainCol; // number of domain columns - double domainRow{0}; // number of domain rows - double deltaCol; // delta in x-direction (between columns) - double deltaRow{0}; // delta in y-direction (between rows) - Eigen::MatrixXd concentrations; // Matrix holding grid concentrations - Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction - Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction + int col; // number of grid columns + int row{1}; // number of grid rows + int dim; // 1D or 2D + T domainCol; // number of domain columns + T domainRow{0}; // number of domain rows + T deltaCol; // delta in x-direction (between columns) + T deltaRow{0}; // delta in y-direction (between rows) + Eigen::MatrixX concentrations; // Matrix holding grid concentrations + Eigen::MatrixX alphaX; // Matrix holding alpha coefficients in x-direction + Eigen::MatrixX alphaY; // Matrix holding alpha coefficients in y-direction + + static constexpr double MAT_INIT_VAL = 0; }; +using Grid64 = Grid; +using Grid32 = Grid; + #endif // GRID_H_ diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 0c09a7e..7ed50c4 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -11,6 +11,11 @@ #include "Boundary.hpp" #include "Grid.hpp" +#include +#include +#include +#include +#include #include #include @@ -77,7 +82,8 @@ enum TIME_MEASURE { * time step, number of simulations, etc. * */ -class Simulation { + +template class Simulation { public: /** * @brief Set up a simulation environment. The timestep and number of @@ -91,7 +97,8 @@ public: * @param bc Valid boundary condition object * @param approach Approach to solving the problem. Either FTCS or BTCS. */ - Simulation(Grid &grid, Boundary &bc, APPROACH approach); + Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach) + : grid(_grid), bc(_bc), approach(_approach){}; /** * @brief Set the option to output the results to a CSV file. Off by default. @@ -107,7 +114,13 @@ public: * - CSV_OUTPUT_XTREME: produce csv output with all * concentration matrices and simulation environment */ - void setOutputCSV(CSV_OUTPUT csv_output); + void setOutputCSV(CSV_OUTPUT csv_output) { + if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid CSV output option given!"); + } + + this->csv_output = csv_output; + } /** * @brief Set the options for outputting information to the console. Off by @@ -122,7 +135,14 @@ public: * - CONSOLE_OUTPUT_VERBOSE: print all concentration * matrices to console */ - void setOutputConsole(CONSOLE_OUTPUT console_output); + void setOutputConsole(CONSOLE_OUTPUT console_output) { + if (console_output < CONSOLE_OUTPUT_OFF && + console_output > CONSOLE_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid console output option given!"); + } + + this->console_output = console_output; + } /** * @brief Set the Time Measure option. Off by default. @@ -133,7 +153,13 @@ public: * - TIME_MEASURE_ON: Time of simulation run is printed to * console */ - void setTimeMeasure(TIME_MEASURE time_measure); + void setTimeMeasure(TIME_MEASURE time_measure) { + if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { + throw std::invalid_argument("Invalid time measure option given!"); + } + + this->time_measure = time_measure; + } /** * @brief Setting the time step for each iteration step. Time step must be @@ -141,14 +167,14 @@ public: * * @param timestep Valid timestep greater than zero. */ - void setTimestep(double timestep); + void setTimestep(T timestep); /** * @brief Currently set time step is returned. * * @return double timestep */ - double getTimestep(); + T getTimestep() const { return this->timestep; } /** * @brief Set the desired iterations to be calculated. A value greater @@ -156,7 +182,13 @@ public: * * @param iterations Number of iterations to be simulated. */ - void setIterations(int iterations); + void setIterations(int iterations) { + if (iterations <= 0) { + throw std::invalid_argument( + "Number of iterations must be greater than zero."); + } + this->iterations = iterations; + } /** * @brief Set the desired linear equation solver to be used for BTCS approach. @@ -165,7 +197,17 @@ public: * @param solver Solver to be used. Default is Thomas Algorithm as it is more * efficient for tridiagonal Matrices. */ - void setSolver(SOLVER solver); + void setSolver(SOLVER solver) { + if (this->approach == FTCS_APPROACH) { + std::cerr + << "Warning: Solver was set, but FTCS approach initialized. Setting " + "the solver " + "is thus without effect." + << std::endl; + } + + this->solver = solver; + } /** * @brief Set the number of desired openMP Threads. @@ -175,20 +217,32 @@ public: * maximum number of processors is set as the default case during Simulation * construction. */ - void setNumberThreads(int num_threads); + void setNumberThreads(int num_threads) { + if (numThreads > 0 && numThreads <= omp_get_num_procs()) { + this->numThreads = numThreads; + } else { + int maxThreadNumber = omp_get_num_procs(); + throw std::invalid_argument( + "Number of threads exceeds the number of processor cores (" + + std::to_string(maxThreadNumber) + ") or is less than 1."); + } + } /** * @brief Return the currently set iterations to be calculated. * * @return int Number of iterations. */ - int getIterations(); + int getIterations() const { return this->iterations; } /** * @brief Outputs the current concentrations of the grid on the console. * */ - void printConcentrationsConsole(); + inline void printConcentrationsConsole() const { + std::cout << grid.getConcentrations() << std::endl; + std::cout << std::endl; + } /** * @brief Creates a CSV file with a name containing the current simulation @@ -199,7 +253,52 @@ public: * * @return string Filename with configured simulation parameters. */ - std::string createCSVfile(); + std::string createCSVfile() const { + std::ofstream file; + int appendIdent = 0; + std::string appendIdentString; + + // 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 filename = + approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; + + while (std::filesystem::exists(filename)) { + appendIdent += 1; + appendIdentString = std::to_string(appendIdent); + filename = approachString + "_" + row + "_" + col + "_" + numIterations + + "-" + appendIdentString + ".csv"; + } + + file.open(filename); + if (!file) { + exit(1); + } + + // 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) { + Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", + " "); + file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) + << std::endl; // boundary left + file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) + << std::endl; // boundary right + file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) + << std::endl; // boundary top + file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) + << std::endl; // boundary bottom + file << std::endl << std::endl; + } + + file.close(); + + return filename; + } /** * @brief Writes the currently calculated concentration values of the grid @@ -208,7 +307,19 @@ public: * @param filename Name of the file to which the concentration values are * to be written. */ - void printConcentrationsCSV(const std::string &filename); + void printConcentrationsCSV(const std::string &filename) const { + std::ofstream file; + + file.open(filename, std::ios_base::app); + if (!file) { + exit(1); + } + + Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << std::endl; + file << std::endl << std::endl; + file.close(); + } /** * @brief Method starts the simulation process with the previously set @@ -217,7 +328,7 @@ public: void run(); private: - double timestep{-1}; + T timestep{-1}; int iterations{-1}; int innerIterations{1}; int numThreads{omp_get_num_procs()}; @@ -225,8 +336,8 @@ private: CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF}; TIME_MEASURE time_measure{TIME_MEASURE_OFF}; - Grid &grid; - Boundary &bc; + Grid &grid; + Boundary &bc; APPROACH approach; SOLVER solver{THOMAS_ALGORITHM_SOLVER}; diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index fe15455..bc80a8b 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -103,23 +103,23 @@ int main(int argc, char *argv[]) { // create a grid with a 5 x 10 field constexpr int row = 5; constexpr int col = 10; - Grid grid = Grid(row, col); + Grid64 grid(row, col); // (optional) set the domain, e.g.: grid.setDomain(0.005, 0.01); const auto init_values_vec = CSVToVector(INPUT_CONC_FILE); - MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); + Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); grid.setConcentrations(concentrations); const double sum_init = concentrations.sum(); // // (optional) set alphas of the grid, e.g.: const auto alphax_vec = CSVToVector(INPUT_ALPHAX_FILE); - MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); + Eigen::MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); constexpr double alphay_val = 5e-10; - MatrixXd alphay = MatrixXd::Constant(row, col, alphay_val); // row,col,value + Eigen::MatrixXd alphay = Eigen::MatrixXd::Constant(row, col, alphay_val); // row,col,value grid.setAlpha(alphax, alphay); // // ****************** diff --git a/src/BTCS.cpp b/src/BTCS.cpp index 3b314e7..7a9dc7e 100644 --- a/src/BTCS.cpp +++ b/src/BTCS.cpp @@ -45,13 +45,15 @@ constexpr std::pair calcBoundaryCoeffClosed(T alpha_center, T alpha_side, } // creates coefficient matrix for next time step from alphas in x-direction -static Eigen::SparseMatrix -createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector &bcLeft, - std::vector &bcRight, int numCols, - int rowIndex, double sx) { +template +static Eigen::SparseMatrix +createCoeffMatrix(const Eigen::MatrixX &alpha, + const std::vector> &bcLeft, + const std::vector> &bcRight, int numCols, + int rowIndex, T sx) { // square matrix of column^2 dimension for the coefficients - Eigen::SparseMatrix cm(numCols, numCols); + Eigen::SparseMatrix cm(numCols, numCols); cm.reserve(Eigen::VectorXi::Constant(numCols, 3)); // left column @@ -142,14 +144,18 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc, // creates a solution vector for next time step from the current state of // concentrations -static Eigen::VectorXd createSolutionVector( - Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX, - Eigen::MatrixXd &alphaY, std::vector &bcLeft, - std::vector &bcRight, std::vector &bcTop, - std::vector &bcBottom, int length, int rowIndex, double sx, - double sy) { +template +static Eigen::VectorX +createSolutionVector(const Eigen::MatrixX &concentrations, + const Eigen::MatrixX &alphaX, + const Eigen::MatrixX &alphaY, + const std::vector> &bcLeft, + const std::vector> &bcRight, + const std::vector> &bcTop, + const std::vector> &bcBottom, + int length, int rowIndex, T sx, T sy) { - Eigen::VectorXd sv(length); + Eigen::VectorX sv(length); const std::size_t numRows = concentrations.rows(); // inner rows @@ -235,10 +241,11 @@ static Eigen::VectorXd createSolutionVector( // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // use of EigenLU solver -static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix &A, - Eigen::VectorXd &b) { +template +static Eigen::VectorX EigenLUAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorX &b) { - Eigen::SparseLU> solver; + Eigen::SparseLU> solver; solver.analyzePattern(A); solver.factorize(A); @@ -248,14 +255,15 @@ static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix &A, // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector // implementation of Thomas Algorithm -static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, - Eigen::VectorXd &b) { +template +static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, + Eigen::VectorX &b) { Eigen::Index n = b.size(); - Eigen::VectorXd a_diag(n); - Eigen::VectorXd b_diag(n); - Eigen::VectorXd c_diag(n); - Eigen::VectorXd x_vec = b; + Eigen::VectorX a_diag(n); + Eigen::VectorX b_diag(n); + Eigen::VectorX c_diag(n); + Eigen::VectorX x_vec = b; // Fill diagonals vectors b_diag[0] = A.coeff(0, 0); @@ -291,23 +299,24 @@ static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix &A, } // BTCS solution for 1D grid -static void -BTCS_1D(Grid &grid, Boundary &bc, double timestep, - Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, - Eigen::VectorXd &b)) { +template +static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, + Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX &b)) { int length = grid.getLength(); - double sx = timestep / (grid.getDelta() * grid.getDelta()); + T sx = timestep / (grid.getDelta() * grid.getDelta()); - Eigen::VectorXd concentrations_t1(length); + Eigen::VectorX concentrations_t1(length); - Eigen::SparseMatrix A; - Eigen::VectorXd b(length); + Eigen::SparseMatrix A; + Eigen::VectorX b(length); - Eigen::MatrixXd alpha = grid.getAlpha(); - std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + const auto &alpha = grid.getAlpha(); - Eigen::MatrixXd concentrations = grid.getConcentrations(); + const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + + Eigen::MatrixX concentrations = grid.getConcentrations(); int rowIndex = 0; A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D @@ -331,31 +340,32 @@ BTCS_1D(Grid &grid, Boundary &bc, double timestep, } // BTCS solution for 2D grid -static void -BTCS_2D(Grid &grid, Boundary &bc, double timestep, - Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix &A, - Eigen::VectorXd &b), - int numThreads) { +template +static void BTCS_2D(Grid &grid, Boundary &bc, T timestep, + Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX &b), + int numThreads) { int rowMax = grid.getRow(); int colMax = grid.getCol(); - double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); - double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); + T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); + T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow()); - Eigen::MatrixXd concentrations_t1 = - Eigen::MatrixXd::Constant(rowMax, colMax, 0); - Eigen::VectorXd row_t1(colMax); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(rowMax, colMax, 0); + Eigen::VectorX row_t1(colMax); - Eigen::SparseMatrix A; - Eigen::VectorXd b; + Eigen::SparseMatrix A; + Eigen::VectorX b; - Eigen::MatrixXd alphaX = grid.getAlphaX(); - Eigen::MatrixXd alphaY = grid.getAlphaY(); - std::vector bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); - std::vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); - std::vector bcTop = bc.getBoundarySide(BC_SIDE_TOP); - std::vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + auto alphaX = grid.getAlphaX(); + auto alphaY = grid.getAlphaY(); - Eigen::MatrixXd concentrations = grid.getConcentrations(); + const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); + const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); + const auto &bcTop = bc.getBoundarySide(BC_SIDE_TOP); + const auto &bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); + + Eigen::MatrixX concentrations = grid.getConcentrations(); #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1) for (int i = 0; i < rowMax; i++) { @@ -364,8 +374,6 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy); - Eigen::SparseLU> solver; - row_t1 = solverFunc(A, b); concentrations_t1.row(i) = row_t1; @@ -394,7 +402,8 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep, } // entry point for EigenLU solver; differentiate between 1D and 2D grid -void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { +template +void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, EigenLUAlgorithm); } else if (grid.getDim() == 2) { @@ -406,7 +415,8 @@ void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) { } // entry point for Thomas algorithm solver; differentiate 1D and 2D grid -void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { +template +void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads) { if (grid.getDim() == 1) { BTCS_1D(grid, bc, timestep, ThomasAlgorithm); } else if (grid.getDim() == 2) { @@ -416,3 +426,13 @@ void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } + +template void BTCS_Thomas(Grid &grid, Boundary &bc, + double timestep, int numThreads); +template void BTCS_Thomas(Grid &grid, Boundary &bc, + float timestep, int numThreads); + +template void BTCS_LU(Grid &grid, Boundary &bc, double timestep, + int numThreads); +template void BTCS_LU(Grid &grid, Boundary &bc, float timestep, + int numThreads); diff --git a/src/Boundary.cpp b/src/Boundary.cpp deleted file mode 100644 index def9bbd..0000000 --- a/src/Boundary.cpp +++ /dev/null @@ -1,149 +0,0 @@ -#include "TugUtils.hpp" - -#include -#include -#include - -BoundaryElement::BoundaryElement() {} - -BoundaryElement::BoundaryElement(double _value) - : value(_value), type(BC_TYPE_CONSTANT) {} - -void BoundaryElement::setType(BC_TYPE type) { this->type = type; } - -void BoundaryElement::setValue(double value) { - if (value < 0) { - throw_invalid_argument("No negative concentration allowed."); - } - if (type == BC_TYPE_CLOSED) { - throw_invalid_argument( - "No constant boundary concentrations can be set for closed " - "boundaries. Please change type first."); - } - this->value = value; -} - -BC_TYPE BoundaryElement::getType() { return this->type; } - -double BoundaryElement::getValue() { return this->value; } - -Boundary::Boundary(Grid grid) : grid(grid) { - if (grid.getDim() == 1) { - this->boundaries = std::vector>( - 2); // in 1D only left and right boundary - - this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement()); - } else if (grid.getDim() == 2) { - this->boundaries = std::vector>(4); - - this->boundaries[BC_SIDE_LEFT] = - std::vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_RIGHT] = - std::vector(grid.getRow(), BoundaryElement()); - this->boundaries[BC_SIDE_TOP] = - std::vector(grid.getCol(), BoundaryElement()); - this->boundaries[BC_SIDE_BOTTOM] = - std::vector(grid.getCol(), BoundaryElement()); - } -} - -void Boundary::setBoundarySideClosed(BC_SIDE side) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - - const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; - const int n = is_vertical ? grid.getRow() : grid.getCol(); - - this->boundaries[side] = std::vector(n, BoundaryElement()); -} - -void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - - const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; - const int n = is_vertical ? grid.getRow() : grid.getCol(); - - this->boundaries[side] = - std::vector(n, BoundaryElement(value)); -} - -void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) { - // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CLOSED); -} - -void Boundary::setBoundaryElementConstant(BC_SIDE side, int index, - double value) { - // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - this->boundaries[side][index].setType(BC_TYPE_CONSTANT); - this->boundaries[side][index].setValue(value); -} - -const std::vector Boundary::getBoundarySide(BC_SIDE side) { - if (grid.getDim() == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw_invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } - return this->boundaries[side]; -} - -Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) { - const std::size_t length = boundaries[side].size(); - Eigen::VectorXd values(length); - - for (int i = 0; i < length; i++) { - if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) { - values(i) = -1; - continue; - } - values(i) = getBoundaryElementValue(side, i); - } - - return values; -} - -BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index]; -} - -BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - return this->boundaries[side][index].getType(); -} - -double Boundary::getBoundaryElementValue(BC_SIDE side, int index) { - if ((boundaries[side].size() < index) || index < 0) { - throw_invalid_argument("Index is selected either too large or too small."); - } - if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { - throw_invalid_argument( - "A value can only be output if it is a constant boundary condition."); - } - return this->boundaries[side][index].getValue(); -} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5527d09..91d141d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCS.cpp) +add_library(tug Simulation.cpp FTCS.cpp BTCS.cpp) IF(TUG_NAAICE_EXAMPLE) target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV) diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 375fe50..5ee4657 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -19,7 +19,8 @@ #endif // calculates horizontal change on one cell independent of boundary type -static inline double calcHorizontalChange(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChange(Grid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -35,7 +36,8 @@ static inline double calcHorizontalChange(Grid &grid, int &row, int &col) { } // calculates vertical change on one cell independent of boundary type -static inline double calcVerticalChange(Grid &grid, int &row, int &col) { +template +static inline T calcVerticalChange(Grid &grid, int &row, int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -51,10 +53,10 @@ static inline double calcVerticalChange(Grid &grid, int &row, int &col) { } // calculates horizontal change on one cell with a constant left boundary -static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcHorizontalChangeLeftBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -68,8 +70,9 @@ static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid, } // calculates horizontal change on one cell with a closed left boundary -static inline double -calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaX()(row, col + 1), grid.getAlphaX()(row, col)) * @@ -78,8 +81,9 @@ calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the left edge of grid -static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, + int &row, int &col) { if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) { return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) { @@ -90,10 +94,10 @@ static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc, } // calculates horizontal change on one cell with a constant right boundary -static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcHorizontalChangeRightBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return 2 * grid.getAlphaX()(row, col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) - @@ -107,8 +111,9 @@ static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid, } // calculates horizontal change on one cell with a closed right boundary -static inline double -calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, + int &col) { return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1), grid.getAlphaX()(row, col)) * @@ -117,8 +122,10 @@ calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the right edge of grid -static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcHorizontalChangeRightBoundary(Grid &grid, + Boundary &bc, int &row, + int &col) { if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) { return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col); } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) { @@ -129,9 +136,10 @@ static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc, } // calculates vertical change on one cell with a constant top boundary -static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, int &col) { +template +static inline T calcVerticalChangeTopBoundaryConstant(Grid &grid, + Boundary &bc, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getAlphaY()(row, col)) * @@ -145,8 +153,9 @@ static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid, } // calculates vertical change on one cell with a closed top boundary -static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, - int &col) { +template +static inline T calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, + int &col) { return calcAlphaIntercell(grid.getAlphaY()(row + 1, col), grid.getConcentrations()(row, col)) * @@ -155,8 +164,9 @@ static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row, } // checks boundary condition type for a cell on the top edge of grid -static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcVerticalChangeTopBoundary(Grid &grid, Boundary &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) { @@ -167,10 +177,10 @@ static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc, } // calculates vertical change on one cell with a constant bottom boundary -static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid, - Boundary &bc, - int &row, - int &col) { +template +static inline T calcVerticalChangeBottomBoundaryConstant(Grid &grid, + Boundary &bc, + int &row, int &col) { return 2 * grid.getAlphaY()(row, col) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) - @@ -184,8 +194,9 @@ static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid, } // calculates vertical change on one cell with a closed bottom boundary -static inline double -calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { +template +static inline T calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, + int &col) { return -(calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row - 1, col)) * @@ -194,8 +205,9 @@ calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) { } // checks boundary condition type for a cell on the bottom edge of grid -static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, - int &row, int &col) { +template +static inline T calcVerticalChangeBottomBoundary(Grid &grid, Boundary &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) { @@ -206,12 +218,14 @@ static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc, } // FTCS solution for 1D grid -static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { +template +static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { int colMax = grid.getCol(); - double deltaCol = grid.getDeltaCol(); + T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(1, colMax, 0); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(1, colMax, 0); // only one row in 1D case -> row constant at index 0 int row = 0; @@ -243,16 +257,17 @@ static void FTCS_1D(Grid &grid, Boundary &bc, double ×tep) { } // FTCS solution for 2D grid -static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, +template +static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, int numThreads) { int rowMax = grid.getRow(); int colMax = grid.getCol(); - double deltaRow = grid.getDeltaRow(); - double deltaCol = grid.getDeltaCol(); + T deltaRow = grid.getDeltaRow(); + T deltaCol = grid.getDeltaCol(); // matrix for concentrations at time t+1 - Eigen::MatrixXd concentrations_t1 = - Eigen::MatrixXd::Constant(rowMax, colMax, 0); + Eigen::MatrixX concentrations_t1 = + Eigen::MatrixX::Constant(rowMax, colMax, 0); // inner cells // these are independent of the boundary condition type @@ -369,7 +384,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double ×tep, } // entry point; differentiate between 1D and 2D grid -void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads) { +template +void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads) { if (grid.getDim() == 1) { FTCS_1D(grid, bc, timestep); } else if (grid.getDim() == 2) { @@ -379,3 +395,8 @@ void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } + +template void FTCS(Grid &grid, Boundary &bc, double timestep, + int &numThreads); +template void FTCS(Grid &grid, Boundary &bc, float timestep, + int &numThreads); diff --git a/src/Grid.cpp b/src/Grid.cpp deleted file mode 100644 index c6573e4..0000000 --- a/src/Grid.cpp +++ /dev/null @@ -1,167 +0,0 @@ -#include "TugUtils.hpp" - -#include -#include - -constexpr double MAT_INIT_VAL = 0; - -Grid::Grid(int length) : col(length), domainCol(length) { - if (length <= 3) { - throw_invalid_argument( - "Given grid length too small. Must be greater than 3."); - } - - this->dim = 1; - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - - this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); -} - -Grid::Grid(int _row, int _col) - : row(_row), col(_col), domainRow(_row), domainCol(_col) { - if (row <= 3 || col <= 3) { - throw_invalid_argument( - "Given grid dimensions too small. Must each be greater than 3."); - } - - this->dim = 2; - this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - - this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); -} - -void Grid::setConcentrations(const Eigen::MatrixXd &concentrations) { - if (concentrations.rows() != this->row || - concentrations.cols() != this->col) { - throw_invalid_argument( - "Given matrix of concentrations mismatch with Grid dimensions!"); - } - - this->concentrations = concentrations; -} - -void Grid::setAlpha(const Eigen::MatrixXd &alpha) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use 2D setter function!"); - } - if (alpha.rows() != 1 || alpha.cols() != this->col) { - throw_invalid_argument( - "Given matrix of alpha coefficients mismatch with Grid dimensions!"); - } - - this->alphaX = alpha; -} - -void Grid::setAlpha(const Eigen::MatrixXd &alphaX, - const Eigen::MatrixXd &alphaY) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably " - "use 1D setter function!"); - } - if (alphaX.rows() != this->row || alphaX.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in x-direction " - "mismatch with GRid dimensions!"); - } - if (alphaY.rows() != this->row || alphaY.cols() != this->col) { - throw_invalid_argument("Given matrix of alpha coefficients in y-direction " - "mismatch with GRid dimensions!"); - } - - this->alphaX = alphaX; - this->alphaY = alphaY; -} - -const Eigen::MatrixXd &Grid::getAlpha() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use either getAlphaX() or getAlphaY()!"); - } - - return this->alphaX; -} - -const Eigen::MatrixXd &Grid::getAlphaX() { - if (dim != 2) { - throw_invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } - - return this->alphaX; -} - -const Eigen::MatrixXd &Grid::getAlphaY() { - if (dim != 2) { - throw_invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } - - return this->alphaY; -} - -int Grid::getDim() { return dim; } - -int Grid::getLength() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use getRow() or getCol()!"); - } - - return col; -} - -int Grid::getRow() { return row; } - -int Grid::getCol() { return col; } - -void Grid::setDomain(double domainLength) { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probaly " - "use the 2D domain setter!"); - } - if (domainLength <= 0) { - throw_invalid_argument("Given domain length is not positive!"); - } - - this->domainCol = domainLength; - this->deltaCol = double(this->domainCol) / double(this->col); -} - -void Grid::setDomain(double domainRow, double domainCol) { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, you should probably " - "use the 1D domain setter!"); - } - if (domainRow <= 0 || domainCol <= 0) { - throw_invalid_argument("Given domain size is not positive!"); - } - - this->domainRow = domainRow; - this->domainCol = domainCol; - this->deltaRow = double(this->domainRow) / double(this->row); - this->deltaCol = double(this->domainCol) / double(this->col); -} - -double Grid::getDelta() { - if (dim != 1) { - throw_invalid_argument("Grid is not one dimensional, you should probably " - "use the 2D delta getters"); - } - - return this->deltaCol; -} - -double Grid::getDeltaCol() { return this->deltaCol; } - -double Grid::getDeltaRow() { - if (dim != 2) { - throw_invalid_argument("Grid is not two dimensional, meaning there is no " - "delta in y-direction!"); - } - - return this->deltaRow; -} diff --git a/src/Schemes.hpp b/src/Schemes.hpp index d10e0c5..771a2a4 100644 --- a/src/Schemes.hpp +++ b/src/Schemes.hpp @@ -16,13 +16,16 @@ #include // entry point; differentiate between 1D and 2D grid -extern void FTCS(Grid &grid, Boundary &bc, double ×tep, int &numThreads); +template +extern void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads); // entry point for EigenLU solver; differentiate between 1D and 2D grid -extern void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads); +template +extern void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads); // entry point for Thomas algorithm solver; differentiate 1D and 2D grid -extern void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, +template +extern void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads); #endif // SCHEMES_H_ diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 0a5a961..90a1459 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -12,56 +12,28 @@ #include "Schemes.hpp" #include "TugUtils.hpp" -Simulation::Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach) - : grid(_grid), bc(_bc), approach(_approach) {} - -void Simulation::setOutputCSV(CSV_OUTPUT csv_output) { - if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid CSV output option given!"); - } - - this->csv_output = csv_output; -} - -void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) { - if (console_output < CONSOLE_OUTPUT_OFF && - console_output > CONSOLE_OUTPUT_VERBOSE) { - throw_invalid_argument("Invalid console output option given!"); - } - - this->console_output = console_output; -} - -void Simulation::setTimeMeasure(TIME_MEASURE time_measure) { - if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) { - throw_invalid_argument("Invalid time measure option given!"); - } - - this->time_measure = time_measure; -} - -void Simulation::setTimestep(double timestep) { +template void Simulation::setTimestep(T timestep) { if (timestep <= 0) { throw_invalid_argument("Timestep has to be greater than zero."); } if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { - const double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); // will be 0 if 1D, else ... - const double deltaRowSquare = grid.getDim() != 1 - ? grid.getDeltaRow() * grid.getDeltaRow() - : deltaColSquare; - const double minDeltaSquare = + const T deltaRowSquare = grid.getDim() != 1 + ? grid.getDeltaRow() * grid.getDeltaRow() + : deltaColSquare; + const T minDeltaSquare = (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - double maxAlpha = std::numeric_limits::quiet_NaN(); + T maxAlpha = std::numeric_limits::quiet_NaN(); // determine maximum alpha if (grid.getDim() == 2) { - const double maxAlphaX = grid.getAlphaX().maxCoeff(); - const double maxAlphaY = grid.getAlphaY().maxCoeff(); + const T maxAlphaX = grid.getAlphaX().maxCoeff(); + const T maxAlphaY = grid.getAlphaY().maxCoeff(); maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; } else if (grid.getDim() == 1) { @@ -74,7 +46,7 @@ void Simulation::setTimestep(double timestep) { const std::string dim = std::to_string(grid.getDim()) + "D"; // Courant-Friedrichs-Lewy condition - double cfl = minDeltaSquare / (4 * maxAlpha); + T cfl = minDeltaSquare / (4 * maxAlpha); // stability equation from Wikipedia; might be useful if applied cfl does // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha * @@ -112,109 +84,7 @@ void Simulation::setTimestep(double timestep) { } } -double Simulation::getTimestep() { return this->timestep; } - -void Simulation::setIterations(int iterations) { - if (iterations <= 0) { - throw_invalid_argument("Number of iterations must be greater than zero."); - } - this->iterations = iterations; -} - -void Simulation::setSolver(SOLVER solver) { - if (this->approach == FTCS_APPROACH) { - std::cerr - << "Warning: Solver was set, but FTCS approach initialized. Setting " - "the solver " - "is thus without effect." - << std::endl; - } - - this->solver = solver; -} - -void Simulation::setNumberThreads(int numThreads) { - if (numThreads > 0 && numThreads <= omp_get_num_procs()) { - this->numThreads = numThreads; - } else { - int maxThreadNumber = omp_get_num_procs(); - std::string outputMessage = - "Number of threads exceeds the number of processor cores (" + - std::to_string(maxThreadNumber) + ") or is less than 1."; - - throw_invalid_argument(outputMessage); - } -} - -int Simulation::getIterations() { return this->iterations; } - -void Simulation::printConcentrationsConsole() { - std::cout << grid.getConcentrations() << std::endl; - std::cout << std::endl; -} - -std::string Simulation::createCSVfile() { - std::ofstream file; - int appendIdent = 0; - std::string appendIdentString; - - // 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 filename = - approachString + "_" + row + "_" + col + "_" + numIterations + ".csv"; - - while (std::filesystem::exists(filename)) { - appendIdent += 1; - appendIdentString = std::to_string(appendIdent); - filename = approachString + "_" + row + "_" + col + "_" + numIterations + - "-" + appendIdentString + ".csv"; - } - - file.open(filename); - if (!file) { - exit(1); - } - - // 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) { - Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "", - " "); - file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row) - << std::endl; // boundary left - file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row) - << std::endl; // boundary right - file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row) - << std::endl; // boundary top - file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row) - << std::endl; // boundary bottom - file << std::endl << std::endl; - } - - file.close(); - - return filename; -} - -void Simulation::printConcentrationsCSV(const std::string &filename) { - std::ofstream file; - - file.open(filename, std::ios_base::app); - if (!file) { - exit(1); - } - - Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << std::endl; - file << std::endl << std::endl; - file.close(); -} - -void Simulation::run() { +template void Simulation::run() { if (this->timestep == -1) { throw_invalid_argument("Timestep is not set!"); } @@ -280,14 +150,14 @@ void Simulation::run() { } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case - double beta = 0.5; + constexpr T beta = 0.5; // TODO this implementation is very inefficient! // a separate implementation that sets up a specific tridiagonal matrix for // Crank-Nicolson would be better - Eigen::MatrixXd concentrations; - Eigen::MatrixXd concentrationsFTCS; - Eigen::MatrixXd concentrationsResult; + Eigen::MatrixX concentrations; + Eigen::MatrixX concentrationsFTCS; + Eigen::MatrixX concentrationsResult; for (int i = 0; i < iterations * innerIterations; i++) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { printConcentrationsConsole(); @@ -324,3 +194,9 @@ void Simulation::run() { << milliseconds.count() << "ms" << std::endl; } } + +template void Simulation::setTimestep(double timestep); +template void Simulation::setTimestep(float timestep); + +template void Simulation::run(); +template void Simulation::run(); diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 1d2bff2..be91b1e 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -6,12 +6,11 @@ #include using namespace std; - TEST_CASE("BoundaryElement") { SUBCASE("Closed case") { - BoundaryElement boundaryElementClosed = BoundaryElement(); - CHECK_NOTHROW(BoundaryElement()); + BoundaryElement boundaryElementClosed = BoundaryElement(); + CHECK_NOTHROW(BoundaryElement()); CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); CHECK_EQ(boundaryElementClosed.getValue(), -1); CHECK_THROWS(boundaryElementClosed.setValue(0.2)); @@ -28,11 +27,11 @@ TEST_CASE("BoundaryElement") { } TEST_CASE("Boundary Class") { - Grid grid1D = Grid(10); - Grid grid2D = Grid(10, 12); + Grid grid1D = Grid64(10); + Grid grid2D = Grid64(10, 12); Boundary boundary1D = Boundary(grid1D); Boundary boundary2D = Boundary(grid2D); - vector boundary1DVector(1, BoundaryElement(1.0)); + vector> boundary1DVector(1, BoundaryElement(1.0)); SUBCASE("Boundaries 1D case") { CHECK_NOTHROW(Boundary boundary(grid1D)); diff --git a/test/testGrid.cpp b/test/testGrid.cpp index a1f2886..ba831d3 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -8,22 +8,22 @@ using namespace std; TEST_CASE("1D Grid, too small length") { int l = 2; - CHECK_THROWS(Grid(l)); + CHECK_THROWS(Grid64(l)); } -TEST_CASE("2D Grid, too small side") { +TEST_CASE("2D Grid64, too small side") { int r = 2; int c = 4; - CHECK_THROWS(Grid(r, c)); + CHECK_THROWS(Grid64(r, c)); r = 4; c = 2; - CHECK_THROWS(Grid(r, c)); + CHECK_THROWS(Grid64(r, c)); } -TEST_CASE("1D Grid") { +TEST_CASE("1D Grid64") { int l = 12; - Grid grid(l); + Grid64 grid(l); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 1); @@ -89,9 +89,9 @@ TEST_CASE("1D Grid") { } } -TEST_CASE("2D Grid quadratic") { +TEST_CASE("2D Grid64 quadratic") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); @@ -170,10 +170,10 @@ TEST_CASE("2D Grid quadratic") { } } -TEST_CASE("2D Grid non-quadratic") { +TEST_CASE("2D Grid64 non-quadratic") { int r = 12; int c = 15; - Grid grid(r, c); + Grid64 grid(r, c); SUBCASE("correct construction") { CHECK_EQ(grid.getDim(), 2); diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index 004799a..e1b0cf3 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -11,7 +11,7 @@ using namespace Eigen; using namespace std; -static Grid setupSimulation(APPROACH approach, double timestep, +static Grid64 setupSimulation(APPROACH approach, double timestep, int iterations) { int row = 11; int col = 11; @@ -19,7 +19,7 @@ static Grid setupSimulation(APPROACH approach, double timestep, int domain_col = 10; // Grid - Grid grid = Grid(row, col); + Grid grid = Grid64(row, col); grid.setDomain(domain_row, domain_col); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); @@ -80,7 +80,7 @@ TEST_CASE("equality to reference matrix with BTCS") { TEST_CASE("Initialize environment") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); Boundary boundary(grid); CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); @@ -88,7 +88,7 @@ TEST_CASE("Initialize environment") { TEST_CASE("Simulation environment") { int rc = 12; - Grid grid(rc, rc); + Grid64 grid(rc, rc); Boundary boundary(grid); Simulation sim(grid, boundary, FTCS_APPROACH); From 5196c36ec58ec61def4084b19065d4c1face0b7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 15 Sep 2023 11:37:42 +0200 Subject: [PATCH 11/17] fix: reintroduce tug namespace --- include/tug/Boundary.hpp | 19 +++++++++++++------ include/tug/Grid.hpp | 4 +++- include/tug/Simulation.hpp | 4 +++- src/BTCS.cpp | 3 +++ src/FTCS.cpp | 3 +++ src/Schemes.hpp | 11 ++++++++--- src/Simulation.cpp | 3 +++ test/testBoundary.cpp | 2 ++ test/testGrid.cpp | 2 +- test/testSimulation.cpp | 1 + 10 files changed, 40 insertions(+), 12 deletions(-) diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index b51df24..4870ae2 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -12,6 +12,8 @@ #include #include +namespace tug { + /** * @brief Enum defining the two implemented boundary conditions. * @@ -191,7 +193,8 @@ public: void setBoundaryElemenClosed(BC_SIDE side, int index) { // tests whether the index really points to an element of the boundary side. if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument("Index is selected either too large or too small."); + throw std::invalid_argument( + "Index is selected either too large or too small."); } this->boundaries[side][index].setType(BC_TYPE_CLOSED); } @@ -211,7 +214,8 @@ public: void setBoundaryElementConstant(BC_SIDE side, int index, double value) { // tests whether the index really points to an element of the boundary side. if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument("Index is selected either too large or too small."); + throw std::invalid_argument( + "Index is selected either too large or too small."); } this->boundaries[side][index].setType(BC_TYPE_CONSTANT); this->boundaries[side][index].setValue(value); @@ -272,7 +276,8 @@ public: */ BoundaryElement getBoundaryElement(BC_SIDE side, int index) const { if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument("Index is selected either too large or too small."); + throw std::invalid_argument( + "Index is selected either too large or too small."); } return this->boundaries[side][index]; } @@ -289,7 +294,8 @@ public: */ BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const { if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument("Index is selected either too large or too small."); + throw std::invalid_argument( + "Index is selected either too large or too small."); } return this->boundaries[side][index].getType(); } @@ -308,7 +314,8 @@ public: */ T getBoundaryElementValue(BC_SIDE side, int index) const { if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument("Index is selected either too large or too small."); + throw std::invalid_argument( + "Index is selected either too large or too small."); } if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { throw std::invalid_argument( @@ -325,5 +332,5 @@ private: std::vector>> boundaries; // Vector with Boundary Element information }; - +} // namespace tug #endif // BOUNDARY_H_ diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 3d3638d..9698184 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -12,6 +12,8 @@ #include #include +namespace tug { + template class Grid { public: /** @@ -321,5 +323,5 @@ private: using Grid64 = Grid; using Grid32 = Grid; - +} // namespace tug #endif // GRID_H_ diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index 7ed50c4..e83d8f4 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -25,6 +25,8 @@ #define omp_get_num_procs() 1 #endif +namespace tug { + /** * @brief Enum defining the two implemented solution approaches. * @@ -343,5 +345,5 @@ private: const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; }; - +} // namespace tug #endif // SIMULATION_H_ diff --git a/src/BTCS.cpp b/src/BTCS.cpp index 7a9dc7e..e166519 100644 --- a/src/BTCS.cpp +++ b/src/BTCS.cpp @@ -22,6 +22,8 @@ #define omp_get_thread_num() 0 #endif +namespace tug { + // calculates coefficient for boundary in constant case template constexpr std::pair calcBoundaryCoeffConstant(T alpha_center, @@ -436,3 +438,4 @@ template void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads); template void BTCS_LU(Grid &grid, Boundary &bc, float timestep, int numThreads); +} // namespace tug diff --git a/src/FTCS.cpp b/src/FTCS.cpp index 5ee4657..914cc41 100644 --- a/src/FTCS.cpp +++ b/src/FTCS.cpp @@ -18,6 +18,8 @@ #define omp_get_thread_num() 0 #endif +namespace tug { + // calculates horizontal change on one cell independent of boundary type template static inline T calcHorizontalChange(Grid &grid, int &row, int &col) { @@ -400,3 +402,4 @@ template void FTCS(Grid &grid, Boundary &bc, double timestep, int &numThreads); template void FTCS(Grid &grid, Boundary &bc, float timestep, int &numThreads); +} // namespace tug diff --git a/src/Schemes.hpp b/src/Schemes.hpp index 771a2a4..09e6d9e 100644 --- a/src/Schemes.hpp +++ b/src/Schemes.hpp @@ -15,17 +15,22 @@ #include #include +namespace tug { + // entry point; differentiate between 1D and 2D grid template -extern void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads); +extern void FTCS(tug::Grid &grid, tug::Boundary &bc, T timestep, + int &numThreads); // entry point for EigenLU solver; differentiate between 1D and 2D grid template -extern void BTCS_LU(Grid &grid, Boundary &bc, T timestep, int numThreads); +extern void BTCS_LU(tug::Grid &grid, tug::Boundary &bc, T timestep, + int numThreads); // entry point for Thomas algorithm solver; differentiate 1D and 2D grid template -extern void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, +extern void BTCS_Thomas(tug::Grid &grid, tug::Boundary &bc, T timestep, int numThreads); +} // namespace tug #endif // SCHEMES_H_ diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 90a1459..ca623d4 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -12,6 +12,8 @@ #include "Schemes.hpp" #include "TugUtils.hpp" +namespace tug { + template void Simulation::setTimestep(T timestep) { if (timestep <= 0) { throw_invalid_argument("Timestep has to be greater than zero."); @@ -200,3 +202,4 @@ template void Simulation::setTimestep(float timestep); template void Simulation::run(); template void Simulation::run(); +} // namespace tug diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index be91b1e..a7b2b5e 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -6,6 +6,8 @@ #include using namespace std; +using namespace tug; + TEST_CASE("BoundaryElement") { SUBCASE("Closed case") { diff --git a/test/testGrid.cpp b/test/testGrid.cpp index ba831d3..d0179dd 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -4,7 +4,7 @@ using namespace Eigen; using namespace std; - +using namespace tug; TEST_CASE("1D Grid, too small length") { int l = 2; diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index e1b0cf3..a7600c5 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -10,6 +10,7 @@ using namespace Eigen; using namespace std; +using namespace tug; static Grid64 setupSimulation(APPROACH approach, double timestep, int iterations) { From 9a3fc678858b2b53db83ad69030741fad8780721 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Mon, 16 Oct 2023 12:11:52 +0200 Subject: [PATCH 12/17] Fix: Eigen::MatrixX instead of Eigen::MatrixXd in Grid.hpp --- include/tug/Grid.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 9698184..0aa9f58 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -59,9 +59,9 @@ public: this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 - this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); - this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL); + this->concentrations = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); + this->alphaY = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); } /** From 8cfb61587d7d877da89c98918f4879c1a738ca77 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Mon, 16 Oct 2023 17:31:43 +0200 Subject: [PATCH 13/17] feat: add naaice/NAAICE_dble_vs_float.cpp --- naaice/CMakeLists.txt | 3 + naaice/NAAICE_dble_vs_float.cpp | 247 ++++++++++++++++++++++++++++++++ 2 files changed, 250 insertions(+) create mode 100644 naaice/NAAICE_dble_vs_float.cpp diff --git a/naaice/CMakeLists.txt b/naaice/CMakeLists.txt index 947e720..b0a45d0 100644 --- a/naaice/CMakeLists.txt +++ b/naaice/CMakeLists.txt @@ -1,9 +1,12 @@ add_executable(naaice BTCS_2D_NAAICE.cpp) +add_executable(NAAICE_dble_vs_float NAAICE_dble_vs_float.cpp) get_filename_component(IN_CONC_FILE "init_conc.csv" REALPATH) get_filename_component(IN_ALPHAX_FILE "alphax.csv" REALPATH) configure_file(files.hpp.in files.hpp) target_include_directories(naaice PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) +target_include_directories(NAAICE_dble_vs_float PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) target_link_libraries(naaice PUBLIC tug) +target_link_libraries(NAAICE_dble_vs_float PUBLIC tug) diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp new file mode 100644 index 0000000..5be74a4 --- /dev/null +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -0,0 +1,247 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "files.hpp" + +using namespace tug; + +/** + * Try to parse an input string into a given template type. + */ +template inline T parseString(const std::string &str) { + T result; + std::istringstream iss(str); + + if (!(iss >> result)) { + throw std::invalid_argument("Invalid input for parsing."); + } + + return result; +} + +/** + * Splits a given string into a vector by using a delimiter character. + */ +template +std::vector tokenize(const std::string &input, char delimiter) { + std::vector tokens; + std::istringstream tokenStream(input); + std::string token; + + while (std::getline(tokenStream, token, delimiter)) { + tokens.push_back(parseString(token)); + } + + return tokens; +} + +/** + * Opens a file containing CSV and transform it into row-major 2D STL vector. + */ +template +std::vector> CSVToVector(const char *filename) { + std::ifstream in_file(filename); + if (!in_file.is_open()) { + throw std::runtime_error("Error opening file \'" + std::string(filename) + + "\'."); + } + + std::vector> csv_data; + + std::string line; + while (std::getline(in_file, line)) { + csv_data.push_back(tokenize(line, ',')); + } + + in_file.close(); + + return csv_data; +} + +/** + * Converts a 2D STL vector, where values are stored row-major into a + * column-major Eigen::Matrix. + */ +template +Eigen::MatrixXd rmVecTocmMatrix(const std::vector> &vec, + std::uint32_t exp_rows, + std::uint32_t exp_cols) { + if (exp_rows != vec.size()) { + throw std::runtime_error( + "Mismatch in y dimension while converting to Eigen::Matrix."); + } + + Eigen::MatrixXd out_mat(exp_rows, exp_cols); + + for (std::uint32_t ri = 0; ri < exp_rows; ri++) { + const auto &vec_row = vec[ri]; + if (vec[ri].size() != exp_cols) { + throw std::runtime_error( + "Mismatch in x dimension while converting to Eigen::Matrix."); + } + for (std::uint32_t cj = 0; cj < exp_cols; cj++) { + out_mat(ri, cj) = vec_row[cj]; + } + } + + return out_mat; +} + + +int DoDble(int ngrid, APPROACH approach) { + + constexpr double dt = 10; + + // create a grid + std::cout << "DOUBLE grid: " << ngrid << std::endl; + + Grid64 grid64(ngrid, ngrid); + + // (optional) set the domain, e.g.: + grid64.setDomain(0.1, 0.1); + + Eigen::MatrixXd initConc64 = Eigen::MatrixXd::Constant(ngrid, ngrid, 0); + initConc64(50, 50) = 1E-6; + grid64.setConcentrations(initConc64); + + const double sum_init64 = initConc64.sum(); + + constexpr double alphax_val = 5e-10; + Eigen::MatrixXd alphax = Eigen::MatrixXd::Constant(ngrid, ngrid, alphax_val); // row,col,value + constexpr double alphay_val = 1e-10; + Eigen::MatrixXd alphay = Eigen::MatrixXd::Constant(ngrid, ngrid, alphay_val); // row,col,value + grid64.setAlpha(alphax, alphay); + + // create a boundary with constant values + Boundary bc64 = Boundary(grid64); + bc64.setBoundarySideClosed(BC_SIDE_LEFT); + bc64.setBoundarySideClosed(BC_SIDE_RIGHT); + bc64.setBoundarySideClosed(BC_SIDE_TOP); + bc64.setBoundarySideClosed(BC_SIDE_BOTTOM); + + // set up a simulation environment + Simulation Sim64 = + Simulation(grid64, bc64, approach); // grid_64,boundary,simulation-approach + + //Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); + + // set the timestep of the simulation + Sim64.setTimestep(dt); // timestep + + // set the number of iterations + Sim64.setIterations(2); + + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + Sim64.setOutputCSV(CSV_OUTPUT_ON); + + // set output to the console to 'ON' + Sim64.setOutputConsole(CONSOLE_OUTPUT_OFF); + + // // **** RUN SIM64 **** + auto begin64 = std::chrono::high_resolution_clock::now(); + + Sim64.run(); + + auto end64 = std::chrono::high_resolution_clock::now(); + auto ms64 = std::chrono::duration_cast(end64 - begin64); + + const double sum_after64 = grid64.getConcentrations().sum(); + + std::cout << "Sum of init field: " << std::setprecision(15) << sum_init64 + << "\nSum after 2 iterations: " << sum_after64 + << "\nMilliseconds: " << ms64.count() << std::endl << std::endl; + return 0; + +} + +int DoFloat(int ngrid, APPROACH approach) { + + constexpr double dt = 10; + + // create a grid + std::cout << "FLOAT grid: " << ngrid << std::endl; + + Grid32 grid32(ngrid, ngrid); + + // (optional) set the domain, e.g.: + grid32.setDomain(0.1, 0.1); + + Eigen::MatrixXf initConc32 = Eigen::MatrixXf::Constant(ngrid, ngrid, 0); + initConc32(50, 50) = 1E-6; + grid32.setConcentrations(initConc32); + + const float sum_init32 = initConc32.sum(); + + constexpr float alphax_val = 5e-10; + Eigen::MatrixXf alphax = Eigen::MatrixXf::Constant(ngrid, ngrid, alphax_val); // row,col,value + constexpr float alphay_val = 1e-10; + Eigen::MatrixXf alphay = Eigen::MatrixXf::Constant(ngrid, ngrid, alphay_val); // row,col,value + grid32.setAlpha(alphax, alphay); + + // create a boundary with constant values + Boundary bc32 = Boundary(grid32); + bc32.setBoundarySideClosed(BC_SIDE_LEFT); + bc32.setBoundarySideClosed(BC_SIDE_RIGHT); + bc32.setBoundarySideClosed(BC_SIDE_TOP); + bc32.setBoundarySideClosed(BC_SIDE_BOTTOM); + + // set up a simulation environment + Simulation Sim32 = + Simulation(grid32, bc32, approach); // grid_32,boundary,simulation-approach + + // Sim32.setSolver(THOMAS_ALGORITHM_SOLVER); + + // set the timestep of the simulation + Sim32.setTimestep(dt); // timestep + + // set the number of iterations + Sim32.setIterations(2); + + // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, + // CSV_OUTPUT_VERBOSE] + Sim32.setOutputCSV(CSV_OUTPUT_ON); + + // set output to the console to 'ON' + Sim32.setOutputConsole(CONSOLE_OUTPUT_OFF); + + // // **** RUN SIM32 **** + auto begin32 = std::chrono::high_resolution_clock::now(); + + Sim32.run(); + + auto end32 = std::chrono::high_resolution_clock::now(); + auto ms32 = std::chrono::duration_cast(end32 - begin32); + + const float sum_after32 = grid32.getConcentrations().sum(); + + std::cout << "Sum of init field: " << std::setprecision(15) << sum_init32 + << "\nSum after 2 iteration: " << sum_after32 + << "\nMilliseconds: " << ms32.count() << std::endl << std::endl; + return 0; + +} + + +int main(int argc, char *argv[]) { + + int n[] = {101, 201, 501, 1001, 2001}; + + for (int i = 0; i < std::size(n); i++) { + DoFloat(n[i], BTCS_APPROACH); + DoDble(n[i], BTCS_APPROACH); + DoFloat(n[i], FTCS_APPROACH); + DoDble(n[i], FTCS_APPROACH); + } + + return 0; +} From 0471f3d8f9f12ae43d755a48033c4828bd5a2523 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Mon, 16 Oct 2023 17:34:52 +0200 Subject: [PATCH 14/17] fix: readded "using namespace tug;" in naaice/BTCS_2D_NAAICE.cpp --- naaice/BTCS_2D_NAAICE.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index bc80a8b..c356aea 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -11,6 +11,8 @@ #include "files.hpp" +using namespace tug; + /** * Try to parse an input string into a given template type. */ From 8456f2212d058eeec118f2d65a034fcb8bc55a3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 17 Oct 2023 12:00:35 +0200 Subject: [PATCH 15/17] BREAKING CHANGE: tug as header-only library build: installation of library is now possible --- CMakeLists.txt | 102 ++++++--- cmake/tugConfig.cmake.in | 4 + docs_sphinx/Boundary.rst | 6 +- docs_sphinx/Grid.rst | 2 +- docs_sphinx/Simulation.rst | 8 +- docs_sphinx/developper.rst | 2 - docs_sphinx/images/class_diagram.svg | 2 +- docs_sphinx/theory.rst | 2 +- docs_sphinx/user.rst | 3 +- examples/BTCS_1D_proto_example.cpp | 4 +- examples/BTCS_2D_proto_example.cpp | 3 +- examples/CRNI_2D_proto_example.cpp | 4 +- examples/FTCS_1D_proto_example.cpp | 3 +- examples/FTCS_2D_proto_closed_mdl.cpp | 3 +- examples/FTCS_2D_proto_example.cpp | 3 +- examples/FTCS_2D_proto_example_mdl.cpp | 3 +- examples/profiling_openmp.cpp | 4 +- examples/profiling_speedup.cpp | 4 +- examples/reference-FTCS_2D_closed.cpp | 3 +- include/tug/Boundary.hpp | 18 +- src/BTCS.cpp => include/tug/Core/BTCS.hpp | 19 +- src/FTCS.cpp => include/tug/Core/FTCS.hpp | 13 +- {src => include/tug/Core}/TugUtils.hpp | 0 include/tug/Grid.hpp | 43 ++-- include/tug/Simulation.hpp | 245 ++++++++++++++++++---- naaice/BTCS_2D_NAAICE.cpp | 6 +- src/CMakeLists.txt | 13 -- src/README.md | 36 ---- src/Schemes.hpp | 36 ---- src/Simulation.cpp | 205 ------------------ test/testFTCS.cpp | 2 +- test/testSimulation.cpp | 46 ++-- 32 files changed, 396 insertions(+), 451 deletions(-) create mode 100644 cmake/tugConfig.cmake.in delete mode 100644 docs_sphinx/developper.rst rename src/BTCS.cpp => include/tug/Core/BTCS.hpp (96%) rename src/FTCS.cpp => include/tug/Core/FTCS.hpp (98%) rename {src => include/tug/Core}/TugUtils.hpp (100%) delete mode 100644 src/CMakeLists.txt delete mode 100644 src/README.md delete mode 100644 src/Schemes.hpp delete mode 100644 src/Simulation.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index afed992..78e45c9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,49 +1,93 @@ -#debian stable (currently bullseye) +# debian stable (currently bullseye) cmake_minimum_required(VERSION 3.18) -project(tug CXX) +project( + tug + VERSION 0.4 + LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) find_package(Eigen3 REQUIRED NO_MODULE) find_package(OpenMP) -# find_package(easy_profiler) -# option(EASY_OPTION_LOG "Verbose easy_profiler" 1) -## SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -mfma") -option(TUG_USE_OPENMP "Compile with OpenMP support" ON) +include(GNUInstallDirs) -set(CMAKE_CXX_FLAGS_GENERICOPT "-O3 -march=native" CACHE STRING - "Flags used by the C++ compiler during opt builds." - FORCE) +# find_package(easy_profiler) option(EASY_OPTION_LOG "Verbose easy_profiler" 1) -set(CMAKE_BUILD_TYPE "${CMAKE_BUILD_TYPE}" CACHE STRING - "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel GenericOpt." - FORCE) +# SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -mfma") +option(TUG_USE_OPENMP "Compile tug with OpenMP support" ON) -option(TUG_USE_UNSAFE_MATH_OPT - "Use compiler options to break IEEE compliances by +set(CMAKE_CXX_FLAGS_GENERICOPT + "-O3 -march=native" + CACHE STRING "Flags used by the C++ compiler during opt builds." FORCE) + +set(CMAKE_BUILD_TYPE + "${CMAKE_BUILD_TYPE}" + CACHE + STRING + "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel GenericOpt." + FORCE) + +option( + TUG_USE_UNSAFE_MATH_OPT "Use compiler options to break IEEE compliances by oenabling reordering of instructions when adding/multiplying of floating - points." - OFF) + points." OFF) -if(TUG_USE_UNSAFE_MATH_OPT) - add_compile_options(-ffast-math) + +option(TUG_ENABLE_TESTING "Run tests after succesfull compilation" OFF) + +option(TUG_HANNESPHILIPP_EXAMPLES "Compile example applications" OFF) + +option(TUG_NAAICE_EXAMPLE "Enables NAAICE examples with export of CSV files" + OFF) + +add_library(tug INTERFACE) +target_include_directories( + tug INTERFACE $ + $) + +target_link_libraries(tug INTERFACE Eigen3::Eigen) + +target_compile_features(tug INTERFACE cxx_std_17) + +if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND) + target_link_libraries(tug INTERFACE OpenMP::OpenMP_CXX) endif() -option(TUG_ENABLE_TESTING - "Run tests after succesfull compilation" - OFF) +if(TUG_USE_UNSAFE_MATH_OPT) + target_compile_options(tug INTERFACE -ffast-math) +endif() -option(TUG_HANNESPHILIPP_EXAMPLES - "Compile example applications" - OFF) +install( + TARGETS tug + EXPORT ${PROJECT_NAME}_Targets + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) -option(TUG_NAAICE_EXAMPLE - "Enables NAAICE examples with export of CSV files" - OFF) +include(CMakePackageConfigHelpers) +write_basic_package_version_file( + "tugConfigVersion.cmake" + VERSION ${PROJECT_VERSION} + COMPATIBILITY SameMajorVersion) -add_subdirectory(src) +configure_package_config_file( + "${PROJECT_SOURCE_DIR}/cmake/${PROJECT_NAME}Config.cmake.in" + "${PROJECT_BINARY_DIR}/${PROJECT_NAME}Config.cmake" + INSTALL_DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/${PROJECT_NAME}/cmake) + +install( + EXPORT ${PROJECT_NAME}_Targets + FILE ${PROJECT_NAME}Targets.cmake + NAMESPACE ${PROJECT_NAME}:: + DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/${PROJECT_NAME}/cmake) + +install(FILES "${PROJECT_BINARY_DIR}/${PROJECT_NAME}Config.cmake" + "${PROJECT_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake" + DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/${PROJECT_NAME}/cmake) + +install(DIRECTORY ${PROJECT_SOURCE_DIR}/include/tug DESTINATION include) if(TUG_ENABLE_TESTING) add_subdirectory(test) @@ -53,6 +97,6 @@ if(TUG_HANNESPHILIPP_EXAMPLES) add_subdirectory(examples) endif() -if (TUG_NAAICE_EXAMPLE) +if(TUG_NAAICE_EXAMPLE) add_subdirectory(naaice) endif() diff --git a/cmake/tugConfig.cmake.in b/cmake/tugConfig.cmake.in new file mode 100644 index 0000000..9c15f36 --- /dev/null +++ b/cmake/tugConfig.cmake.in @@ -0,0 +1,4 @@ +@PACKAGE_INIT@ + +include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake") +check_required_components("@PROJECT_NAME@") diff --git a/docs_sphinx/Boundary.rst b/docs_sphinx/Boundary.rst index d22e133..7f3ce75 100644 --- a/docs_sphinx/Boundary.rst +++ b/docs_sphinx/Boundary.rst @@ -1,4 +1,8 @@ Boundary ======== -.. doxygenclass:: Boundary \ No newline at end of file +.. doxygenenum:: tug::BC_TYPE +.. doxygenenum:: tug::BC_SIDE + +.. doxygenclass:: tug::Boundary +.. doxygenclass:: tug::BoundaryElement diff --git a/docs_sphinx/Grid.rst b/docs_sphinx/Grid.rst index c680ab4..e4dc497 100644 --- a/docs_sphinx/Grid.rst +++ b/docs_sphinx/Grid.rst @@ -1,4 +1,4 @@ Grid ==== -.. doxygenclass:: Grid \ No newline at end of file +.. doxygenclass:: tug::Grid diff --git a/docs_sphinx/Simulation.rst b/docs_sphinx/Simulation.rst index a31ec9b..a461073 100644 --- a/docs_sphinx/Simulation.rst +++ b/docs_sphinx/Simulation.rst @@ -1,4 +1,10 @@ Simulation ========== -.. doxygenclass:: Simulation \ No newline at end of file +.. doxygenenum:: tug::APPROACH +.. doxygenenum:: tug::SOLVER +.. doxygenenum:: tug::CSV_OUTPUT +.. doxygenenum:: tug::CONSOLE_OUTPUT +.. doxygenenum:: tug::TIME_MEASURE + +.. doxygenclass:: tug::Simulation diff --git a/docs_sphinx/developper.rst b/docs_sphinx/developper.rst deleted file mode 100644 index 6b97cc5..0000000 --- a/docs_sphinx/developper.rst +++ /dev/null @@ -1,2 +0,0 @@ -Developper API -============== \ No newline at end of file diff --git a/docs_sphinx/images/class_diagram.svg b/docs_sphinx/images/class_diagram.svg index 3392a01..d33462e 100644 --- a/docs_sphinx/images/class_diagram.svg +++ b/docs_sphinx/images/class_diagram.svg @@ -1,4 +1,4 @@ -
Grid
Grid
- dim:int
- col:int
- row:int
- domain_col: int
- domain_row: int
- delta_col: double
- delta_row: double
- dim:int...
+ Grid(col: int)
+ Grid(row: int, col: int)
+ setConcentrations(concentrations: MatrixXd)
+ getConcentrations(): MatrixXd
+ setAlpha(alpha: MatrixXd)
+ setAlpha(alpha_x: MatrixXd, alpha_y: MatrixXd)
+ getAlphaX(): MatrixXd
+ getAlphaY(): MatrixXd
+ getDim(): int
+ getRow(): int
+ getCol(): int
+ setDomain(domain_col: int)
+ setDomain(domain_row: int, domain_col: int)
+ getDeltaCol(): double
+ getDeltaRow(): double
+ Grid(col: int)...
Boundary
Boundary
- grid: Grid
- type: BC_TYPE
- left: VectorXd
- right: VectorXd
- top: VectorXd
- bottom: VectorXd
- grid: Grid...
+ Boundary(grid: Grid, type: BC_TYPE)
+ getBoundaryConditionType(): BC_TYPE
+ setBoundaryConditionValue(side: BC_SIDE, value: VectorXd)
+ getBoundaryConditionValue(): VectorXd
+ Boundary(grid: Grid, type: BC_TYPE)...
Simulation
Simulation
- timestep: double
- iterations: int
- csv_output: CSV_OUTPUT
- grid: GRID
- bc: Boundary
- approach: APPROACH
- timestep: double...
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH
+ setOutputCSV(csv_output: CSV_OUTPUT)
+ setTimestep(timestep: double)
+ getTimestep(): double
+ setIterations(iterations: int)
+ getIterations(): int
+ printConcentrationsConsole()
+ printConcentrationsCSV(ident: string, appendMode: bool)
+ run()
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH...
FTCS
FTCS
BTCS
BTCS
functionalities
functionalities
TUG API
TUG API
BoundaryElement
BoundaryElement
- type: BC_TYPE
- value: double
- type: BC_TYPE...
+ BoundaryElement()
+ BoundaryElement(value: double)
+ setType(type: BC_TYPE)
+ setValue(value: double)
+ getType(): BC_TYPE
+ getValue(): double
+ BoundaryElement()...
Text is not SVG - cannot display
\ No newline at end of file +
Grid<T>
Grid<T>
- dim:int
- col:int
- row:int
- domain_col: int
- domain_row: int
- delta_col: T
- delta_row: T
- dim:int...
+ Grid(col: int)
+ Grid(row: int, col: int)
+ setConcentrations(concentrations: MatrixX<T>)
+ getConcentrations(): MatrixX<T>
+ setAlpha(alpha: MatrixX<T>)
+ setAlpha(alpha_x: MatrixX<T>, alpha_y:   MatrixX<T>)
+ getAlphaX(): MatrixX<T>
+ getAlphaY(): MatrixX<T>
+ getDim(): int
+ getRow(): int
+ getCol(): int
+ setDomain(domain_col: int)
+ setDomain(domain_row: int, domain_col: int)
+ getDeltaCol(): T
+ getDeltaRow(): T
+ Grid(col: int)...
Boundary<T>
Boundary<T>
- grid: Grid
- type: BC_TYPE
- left: std::vector<T>
- right: std::vector<T>
- top: std::vector<T>
- bottom: std::vector<T>
- grid: Grid...
+ Boundary(grid: Grid, type: BC_TYPE)
+ getBoundaryConditionType(): BC_TYPE
+ setBoundaryConditionValue(side: BC_SIDE, 
value: std::vector<T>)
+ getBoundaryConditionValue(): std::vector<T>
+ Boundary(grid: Grid, type: BC_TYPE)...
Simulation<T, approach, solver>
Simulation<T, approach, solver>
- timestep: double
- iterations: int
- csv_output: CSV_OUTPUT
- grid: GRID
- bc: Boundary
- timestep: double...
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH
+ setOutputCSV(csv_output: CSV_OUTPUT)
+ setTimestep(timestep: T)
+ getTimestep(): T
+ setIterations(iterations: int)
+ getIterations(): int
+ printConcentrationsConsole()
+ printConcentrationsCSV(ident: string, appendMode: bool)
+ run()
+ Simulation(grid: Grid, bc: Boundary, approach: APPROACH...
FTCS
FTCS
BTCS
BTCS
functionalities
functionalities
TUG API
TUG API
BoundaryElement<T>
BoundaryElement<T>
- type: BC_TYPE
- value: T
- type: BC_TYPE...
+ BoundaryElement()
+ BoundaryElement(value: T)
+ setType(type: BC_TYPE)
+ setValue(value: T)
+ getType(): BC_TYPE
+ getValue(): T
+ BoundaryElement()...
Text is not SVG - cannot display
\ No newline at end of file diff --git a/docs_sphinx/theory.rst b/docs_sphinx/theory.rst index 27c8146..4d3468f 100644 --- a/docs_sphinx/theory.rst +++ b/docs_sphinx/theory.rst @@ -12,4 +12,4 @@ Numerical Solver **Backward Time-Centered Space (BTCS) Method** -**Forward Time-Centered Space (BTCS) Method** \ No newline at end of file +**Forward Time-Centered Space (FTCS) Method** diff --git a/docs_sphinx/user.rst b/docs_sphinx/user.rst index 8d73236..ddd4941 100644 --- a/docs_sphinx/user.rst +++ b/docs_sphinx/user.rst @@ -3,7 +3,6 @@ User API .. toctree:: - :maxdepth: 2 Grid Boundary - Simulation \ No newline at end of file + Simulation diff --git a/examples/BTCS_1D_proto_example.cpp b/examples/BTCS_1D_proto_example.cpp index 4bfd1b4..af9888e 100644 --- a/examples/BTCS_1D_proto_example.cpp +++ b/examples/BTCS_1D_proto_example.cpp @@ -2,6 +2,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { // ************** @@ -31,8 +32,7 @@ int main(int argc, char *argv[]) { // ************************ // set up a simulation environment - Simulation simulation = - Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + Simulation simulation = Simulation(grid, bc); // grid,boundary // set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index f14e0f5..98acdb9 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -2,6 +2,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { // EASY_PROFILER_ENABLE; @@ -61,7 +62,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary // set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/CRNI_2D_proto_example.cpp b/examples/CRNI_2D_proto_example.cpp index d605038..fc3f046 100644 --- a/examples/CRNI_2D_proto_example.cpp +++ b/examples/CRNI_2D_proto_example.cpp @@ -2,6 +2,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { int row = 20; @@ -18,7 +19,8 @@ int main(int argc, char *argv[]) { bc.setBoundarySideClosed(BC_SIDE_TOP); bc.setBoundarySideClosed(BC_SIDE_BOTTOM); - Simulation simulation = Simulation(grid, bc, CRANK_NICOLSON_APPROACH); + Simulation simulation = + Simulation(grid, bc); simulation.setTimestep(0.1); simulation.setIterations(50); simulation.setOutputCSV(CSV_OUTPUT_XTREME); diff --git a/examples/FTCS_1D_proto_example.cpp b/examples/FTCS_1D_proto_example.cpp index 74b6f3e..b3dc905 100644 --- a/examples/FTCS_1D_proto_example.cpp +++ b/examples/FTCS_1D_proto_example.cpp @@ -3,6 +3,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { // ************** @@ -31,7 +32,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary,simulation-approach // (optional) set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index e17f92b..c6ed02f 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -12,6 +12,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { @@ -69,7 +70,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary,simulation-approach // set the timestep of the simulation simulation.setTimestep(10000); // timestep diff --git a/examples/FTCS_2D_proto_example.cpp b/examples/FTCS_2D_proto_example.cpp index 89f821d..9b8987a 100644 --- a/examples/FTCS_2D_proto_example.cpp +++ b/examples/FTCS_2D_proto_example.cpp @@ -10,6 +10,7 @@ #include using namespace Eigen; +using namespace tug; // #include // #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true); @@ -67,7 +68,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary,simulation-approach // set the timestep of the simulation simulation.setTimestep(0.1); // timestep diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index d77710d..fdf40d1 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -10,6 +10,7 @@ #include using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { @@ -60,7 +61,7 @@ int main(int argc, char *argv[]) { // set up a simulation environment Simulation simulation = - Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach + Simulation(grid, bc); // grid,boundary,simulation-approach // (optional) set the timestep of the simulation simulation.setTimestep(1000); // timestep diff --git a/examples/profiling_openmp.cpp b/examples/profiling_openmp.cpp index c781900..8ad62e0 100644 --- a/examples/profiling_openmp.cpp +++ b/examples/profiling_openmp.cpp @@ -7,6 +7,7 @@ using namespace Eigen; using namespace std; +using namespace tug; int main(int argc, char *argv[]) { @@ -39,8 +40,7 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid); - Simulation sim = Simulation(grid, bc, BTCS_APPROACH); - sim.setSolver(THOMAS_ALGORITHM_SOLVER); + Simulation sim = Simulation(grid, bc); if (argc == 2) { int numThreads = atoi(argv[1]); diff --git a/examples/profiling_speedup.cpp b/examples/profiling_speedup.cpp index c781900..8ad62e0 100644 --- a/examples/profiling_speedup.cpp +++ b/examples/profiling_speedup.cpp @@ -7,6 +7,7 @@ using namespace Eigen; using namespace std; +using namespace tug; int main(int argc, char *argv[]) { @@ -39,8 +40,7 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid); - Simulation sim = Simulation(grid, bc, BTCS_APPROACH); - sim.setSolver(THOMAS_ALGORITHM_SOLVER); + Simulation sim = Simulation(grid, bc); if (argc == 2) { int numThreads = atoi(argv[1]); diff --git a/examples/reference-FTCS_2D_closed.cpp b/examples/reference-FTCS_2D_closed.cpp index fa16063..f08c689 100644 --- a/examples/reference-FTCS_2D_closed.cpp +++ b/examples/reference-FTCS_2D_closed.cpp @@ -4,6 +4,7 @@ using namespace std; using namespace Eigen; +using namespace tug; int main(int argc, char *argv[]) { int row = 50; @@ -41,7 +42,7 @@ int main(int argc, char *argv[]) { Boundary bc = Boundary(grid); // Simulation - Simulation sim = Simulation(grid, bc, FTCS_APPROACH); + Simulation sim = Simulation(grid, bc); sim.setTimestep(0.001); sim.setIterations(10000); sim.setOutputCSV(CSV_OUTPUT_OFF); diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index 4870ae2..3a7f9b8 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -30,6 +30,8 @@ enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM }; * This class defines the boundary conditions of individual boundary elements. * These can be flexibly used and combined later in other classes. * The class serves as an auxiliary class for structuring the Boundary class. + * + * @tparam T Data type of the boundary condition element */ template class BoundaryElement { public: @@ -82,7 +84,7 @@ public: * @brief Return the type of the boundary condition, i.e. whether the * boundary is considered closed or constant. * - * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or + * @return Type of boundary condition, either BC_TYPE_CLOSED or BC_TYPE_CONSTANT. */ BC_TYPE getType() const { return this->type; } @@ -90,7 +92,7 @@ public: /** * @brief Return the concentration value for the constant boundary condition. * - * @return double Value of the concentration. + * @return Value of the concentration. */ T getValue() const { return this->value; } @@ -102,6 +104,8 @@ private: /** * This class implements the functionality and management of the boundary * conditions in the grid to be simulated. + * + * @tparam Data type of the boundary condition value */ template class Boundary { public: @@ -227,7 +231,7 @@ public: * * @param side Boundary side from which the boundary conditions are to be * returned. - * @return vector> Contains the boundary conditions as + * @return Contains the boundary conditions as * BoundaryElement objects. */ const std::vector> &getBoundarySide(BC_SIDE side) const { @@ -246,7 +250,7 @@ public: specific boundary is closed. * * @param side Boundary side for which the values are to be returned. - * @return VectorX Vector with values as T. + * @return Vector with values as T. */ Eigen::VectorX getBoundarySideValues(BC_SIDE side) const { const std::size_t length = boundaries[side].size(); @@ -271,7 +275,7 @@ public: * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return BoundaryElement Boundary condition as a BoundaryElement + * @return Boundary condition as a BoundaryElement * object. */ BoundaryElement getBoundaryElement(BC_SIDE side, int index) const { @@ -290,7 +294,7 @@ public: * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding side. - * @return BC_TYPE Boundary Type of the corresponding boundary condition. + * @return Boundary Type of the corresponding boundary condition. */ BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const { if ((boundaries[side].size() < index) || index < 0) { @@ -309,7 +313,7 @@ public: * @param index Index of the boundary element on the corresponding * boundary side. Must index an element of the corresponding * side. - * @return double Concentration of the corresponding BoundaryElement + * @return Concentration of the corresponding BoundaryElement * object. */ T getBoundaryElementValue(BC_SIDE side, int index) const { diff --git a/src/BTCS.cpp b/include/tug/Core/BTCS.hpp similarity index 96% rename from src/BTCS.cpp rename to include/tug/Core/BTCS.hpp index e166519..cea830a 100644 --- a/src/BTCS.cpp +++ b/include/tug/Core/BTCS.hpp @@ -1,5 +1,5 @@ /** - * @file BTCSv2.cpp + * @file BTCS.hpp * @brief Implementation of heterogenous BTCS (backward time-centered space) * solution of diffusion equation in 1D and 2D space. Internally the * alternating-direction implicit (ADI) method is used. Version 2, because @@ -7,10 +7,11 @@ * */ -#include "Schemes.hpp" +#ifndef BTCS_H_ +#define BTCS_H_ + #include "TugUtils.hpp" -#include #include #include #include @@ -428,14 +429,6 @@ void BTCS_Thomas(Grid &grid, Boundary &bc, T timestep, int numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } - -template void BTCS_Thomas(Grid &grid, Boundary &bc, - double timestep, int numThreads); -template void BTCS_Thomas(Grid &grid, Boundary &bc, - float timestep, int numThreads); - -template void BTCS_LU(Grid &grid, Boundary &bc, double timestep, - int numThreads); -template void BTCS_LU(Grid &grid, Boundary &bc, float timestep, - int numThreads); } // namespace tug + +#endif // BTCS_H_ diff --git a/src/FTCS.cpp b/include/tug/Core/FTCS.hpp similarity index 98% rename from src/FTCS.cpp rename to include/tug/Core/FTCS.hpp index 914cc41..7f12780 100644 --- a/src/FTCS.cpp +++ b/include/tug/Core/FTCS.hpp @@ -1,11 +1,13 @@ /** - * @file FTCS.cpp + * @file FTCS.hpp * @brief Implementation of heterogenous FTCS (forward time-centered space) * solution of diffusion equation in 1D and 2D space. * */ -#include "Schemes.hpp" +#ifndef FTCS_H_ +#define FTCS_H_ + #include "TugUtils.hpp" #include @@ -397,9 +399,6 @@ void FTCS(Grid &grid, Boundary &bc, T timestep, int &numThreads) { "Error: Only 1- and 2-dimensional grids are defined!"); } } - -template void FTCS(Grid &grid, Boundary &bc, double timestep, - int &numThreads); -template void FTCS(Grid &grid, Boundary &bc, float timestep, - int &numThreads); } // namespace tug + +#endif // FTCS_H_ diff --git a/src/TugUtils.hpp b/include/tug/Core/TugUtils.hpp similarity index 100% rename from src/TugUtils.hpp rename to include/tug/Core/TugUtils.hpp diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index 0aa9f58..dee004a 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -14,6 +14,12 @@ 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 Grid { public: /** @@ -31,10 +37,11 @@ public: } this->dim = 1; - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(this->col); // -> 1 - this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); - this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL); + this->concentrations = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); + this->alphaX = Eigen::MatrixX::Constant(1, col, MAT_INIT_VAL); } /** @@ -56,8 +63,10 @@ public: } this->dim = 2; - this->deltaRow = double(this->domainRow) / double(this->row); // -> 1 - this->deltaCol = double(this->domainCol) / double(this->col); // -> 1 + this->deltaRow = + static_cast(this->domainRow) / static_cast(this->row); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(this->col); // -> 1 this->concentrations = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); this->alphaX = Eigen::MatrixX::Constant(row, col, MAT_INIT_VAL); @@ -84,7 +93,7 @@ public: /** * @brief Gets the concentrations matrix for a Grid. * - * @return MatrixX An Eigen3 matrix holding the concentrations and having + * @return An Eigen3 matrix holding the concentrations and having * the same dimensions as the grid. */ const Eigen::MatrixX &getConcentrations() { return this->concentrations; } @@ -145,7 +154,7 @@ public: * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one * dimensional. * - * @return MatrixX A matrix with 1 row holding the alpha coefficients. + * @return A matrix with 1 row holding the alpha coefficients. */ const Eigen::MatrixX &getAlpha() const { if (dim != 1) { @@ -161,7 +170,7 @@ public: * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixX A matrix holding the alpha coefficients in x-direction. + * @return A matrix holding the alpha coefficients in x-direction. */ const Eigen::MatrixX &getAlphaX() const { @@ -177,7 +186,7 @@ public: * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid. * Grid must be two dimensional. * - * @return MatrixX A matrix holding the alpha coefficients in y-direction. + * @return A matrix holding the alpha coefficients in y-direction. */ const Eigen::MatrixX &getAlphaY() const { @@ -192,14 +201,14 @@ public: /** * @brief Gets the dimensions of the grid. * - * @return int Dimensions, either 1 or 2. + * @return Dimensions, either 1 or 2. */ int getDim() const { return this->dim; } /** * @brief Gets length of 1D grid. Must be one dimensional grid. * - * @return int Length of 1D grid. + * @return Length of 1D grid. */ int getLength() const { if (dim != 1) { @@ -214,14 +223,14 @@ public: /** * @brief Gets the number of rows of the grid. * - * @return int Number of rows. + * @return Number of rows. */ int getRow() const { return this->row; } /** * @brief Gets the number of columns of the grid. * - * @return int Number of columns. + * @return Number of columns. */ int getCol() const { return this->col; } @@ -271,7 +280,7 @@ public: /** * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. * - * @return double Delta value. + * @return Delta value. */ T getDelta() const { @@ -287,14 +296,14 @@ public: /** * @brief Gets the delta value in x-direction. * - * @return double 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 double Delta value in y-direction. + * @return Delta value in y-direction. */ T getDeltaRow() const { if (dim != 2) { @@ -318,7 +327,7 @@ private: Eigen::MatrixX alphaX; // Matrix holding alpha coefficients in x-direction Eigen::MatrixX alphaY; // Matrix holding alpha coefficients in y-direction - static constexpr double MAT_INIT_VAL = 0; + static constexpr T MAT_INIT_VAL = 0; }; using Grid64 = Grid; diff --git a/include/tug/Simulation.hpp b/include/tug/Simulation.hpp index e83d8f4..4c0f5c8 100644 --- a/include/tug/Simulation.hpp +++ b/include/tug/Simulation.hpp @@ -11,7 +11,7 @@ #include "Boundary.hpp" #include "Grid.hpp" -#include +#include #include #include #include @@ -19,6 +19,10 @@ #include #include +#include "Core/BTCS.hpp" +#include "Core/FTCS.hpp" +#include "Core/TugUtils.hpp" + #ifdef _OPENMP #include #else @@ -28,13 +32,13 @@ namespace tug { /** - * @brief Enum defining the two implemented solution approaches. + * @brief Enum defining the implemented solution approaches. * */ enum APPROACH { - FTCS_APPROACH, // Forward Time-Centered Space - BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver - CRANK_NICOLSON_APPROACH + FTCS_APPROACH, /*!< Forward Time-Centered Space */ + BTCS_APPROACH, /*!< Backward Time-Centered Space */ + CRANK_NICOLSON_APPROACH /*!< Crank-Nicolson method */ }; /** @@ -42,9 +46,9 @@ enum APPROACH { * */ enum SOLVER { - EIGEN_LU_SOLVER, // EigenLU solver - THOMAS_ALGORITHM_SOLVER // Thomas Algorithm solver; more efficient for - // tridiagonal matrices + EIGEN_LU_SOLVER, /*!< EigenLU solver */ + THOMAS_ALGORITHM_SOLVER /*!< Thomas Algorithm solver; more efficient for + tridiagonal matrices */ }; /** @@ -52,11 +56,11 @@ enum SOLVER { * */ enum CSV_OUTPUT { - CSV_OUTPUT_OFF, // do not produce csv output - CSV_OUTPUT_ON, // produce csv output with last concentration matrix - CSV_OUTPUT_VERBOSE, // produce csv output with all concentration matrices - CSV_OUTPUT_XTREME // csv output like VERBOSE but additional boundary - // conditions at beginning + CSV_OUTPUT_OFF, /*!< do not produce csv output */ + CSV_OUTPUT_ON, /*!< produce csv output with last concentration matrix */ + CSV_OUTPUT_VERBOSE, /*!< produce csv output with all concentration matrices */ + CSV_OUTPUT_XTREME /*!< csv output like VERBOSE but additional boundary + conditions at beginning */ }; /** @@ -64,9 +68,9 @@ enum CSV_OUTPUT { * */ enum CONSOLE_OUTPUT { - CONSOLE_OUTPUT_OFF, // do not print any output to console - CONSOLE_OUTPUT_ON, // print before and after concentrations to console - CONSOLE_OUTPUT_VERBOSE // print all concentration matrices to console + CONSOLE_OUTPUT_OFF, /*!< do not print any output to console */ + CONSOLE_OUTPUT_ON, /*!< print before and after concentrations to console */ + CONSOLE_OUTPUT_VERBOSE /*!< print all concentration matrices to console */ }; /** @@ -74,8 +78,8 @@ enum CONSOLE_OUTPUT { * */ enum TIME_MEASURE { - TIME_MEASURE_OFF, // do not print any time measures - TIME_MEASURE_ON // print time measure after last iteration + TIME_MEASURE_OFF, /*!< do not print any time measures */ + TIME_MEASURE_ON /*!< print time measure after last iteration */ }; /** @@ -83,9 +87,14 @@ enum TIME_MEASURE { * and contains all the methods for controlling the desired parameters, such as * time step, number of simulations, etc. * + * @tparam T the type of the internal data structures for grid, boundary + * condition and timestep + * @tparam approach Set the SLE scheme to be used + * @tparam solver Set the solver to be used */ - -template class Simulation { +template +class Simulation { public: /** * @brief Set up a simulation environment. The timestep and number of @@ -99,8 +108,7 @@ public: * @param bc Valid boundary condition object * @param approach Approach to solving the problem. Either FTCS or BTCS. */ - Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach) - : grid(_grid), bc(_bc), approach(_approach){}; + Simulation(Grid &_grid, Boundary &_bc) : grid(_grid), bc(_bc){}; /** * @brief Set the option to output the results to a CSV file. Off by default. @@ -169,7 +177,65 @@ public: * * @param timestep Valid timestep greater than zero. */ - void setTimestep(T timestep); + void setTimestep(T timestep) { + if (timestep <= 0) { + throw_invalid_argument("Timestep has to be greater than zero."); + } + + if constexpr (approach == FTCS_APPROACH || + approach == CRANK_NICOLSON_APPROACH) { + T cfl; + if (grid.getDim() == 1) { + + const T deltaSquare = grid.getDelta(); + const T maxAlpha = grid.getAlpha().maxCoeff(); + + // Courant-Friedrichs-Lewy condition + cfl = deltaSquare / (4 * maxAlpha); + } else if (grid.getDim() == 2) { + const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + // will be 0 if 1D, else ... + const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); + + const T maxAlpha = + std::min(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); + + cfl = minDeltaSquare / (4 * maxAlpha); + } + const std::string dim = std::to_string(grid.getDim()) + "D"; + + const std::string &approachPrefix = this->approach_names[approach]; + std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl + << std::endl; + std::cout << approachPrefix << "_" << dim + << " :: required dt=" << timestep << std::endl; + + if (timestep > cfl) { + + this->innerIterations = (int)ceil(timestep / cfl); + this->timestep = timestep / (double)innerIterations; + + std::cerr << "Warning :: Timestep was adjusted, because of stability " + "conditions. Time duration was approximately preserved by " + "adjusting internal number of iterations." + << std::endl; + std::cout << approachPrefix << "_" << dim << " :: Required " + << this->innerIterations + << " inner iterations with dt=" << this->timestep + << std::endl; + + } else { + + this->timestep = timestep; + std::cout << approachPrefix << "_" << dim + << " :: No inner iterations required, dt=" << timestep + << std::endl; + } + } else { + this->timestep = timestep; + } + } /** * @brief Currently set time step is returned. @@ -192,25 +258,6 @@ public: this->iterations = iterations; } - /** - * @brief Set the desired linear equation solver to be used for BTCS approach. - * Without effect in case of FTCS approach. - * - * @param solver Solver to be used. Default is Thomas Algorithm as it is more - * efficient for tridiagonal Matrices. - */ - void setSolver(SOLVER solver) { - if (this->approach == FTCS_APPROACH) { - std::cerr - << "Warning: Solver was set, but FTCS approach initialized. Setting " - "the solver " - "is thus without effect." - << std::endl; - } - - this->solver = solver; - } - /** * @brief Set the number of desired openMP Threads. * @@ -327,7 +374,117 @@ public: * @brief Method starts the simulation process with the previously set * parameters. */ - void run(); + void run() { + if (this->timestep == -1) { + throw_invalid_argument("Timestep is not set!"); + } + if (this->iterations == -1) { + throw_invalid_argument("Number of iterations are not set!"); + } + + std::string filename; + if (this->console_output > CONSOLE_OUTPUT_OFF) { + printConcentrationsConsole(); + } + if (this->csv_output > CSV_OUTPUT_OFF) { + filename = createCSVfile(); + } + + auto begin = std::chrono::high_resolution_clock::now(); + + if constexpr (approach == FTCS_APPROACH) { // FTCS case + + for (int i = 0; i < iterations * innerIterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + FTCS(this->grid, this->bc, this->timestep, this->numThreads); + + // if (i % (iterations * innerIterations / 100) == 0) { + // double percentage = (double)i / ((double)iterations * + // (double)innerIterations) * 100; if ((int)percentage % 10 == 0) { + // cout << "Progress: " << percentage << "%" << endl; + // } + // } + } + + } 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) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads); + } + } else if constexpr (solver == THOMAS_ALGORITHM_SOLVER) { + for (int i = 0; i < iterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); + } + } + + } else if constexpr (approach == + CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case + + constexpr T beta = 0.5; + + // TODO this implementation is very inefficient! + // a separate implementation that sets up a specific tridiagonal matrix + // for Crank-Nicolson would be better + Eigen::MatrixX concentrations; + Eigen::MatrixX concentrationsFTCS; + Eigen::MatrixX concentrationsResult; + for (int i = 0; i < iterations * innerIterations; i++) { + if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { + printConcentrationsConsole(); + } + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationsCSV(filename); + } + + concentrations = grid.getConcentrations(); + 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); + } + } + + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = + std::chrono::duration_cast(end - begin); + + if (this->console_output > CONSOLE_OUTPUT_OFF) { + printConcentrationsConsole(); + } + if (this->csv_output > CSV_OUTPUT_OFF) { + printConcentrationsCSV(filename); + } + if (this->time_measure > TIME_MEASURE_OFF) { + const std::string &approachString = this->approach_names[approach]; + const std::string dimString = std::to_string(grid.getDim()) + "D"; + std::cout << approachString << dimString << ":: run() finished in " + << milliseconds.count() << "ms" << std::endl; + } + } private: T timestep{-1}; @@ -340,8 +497,6 @@ private: Grid &grid; Boundary &bc; - APPROACH approach; - SOLVER solver{THOMAS_ALGORITHM_SOLVER}; const std::vector approach_names = {"FTCS", "BTCS", "CRNI"}; }; diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index c356aea..9ed7ccb 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -121,7 +121,8 @@ int main(int argc, char *argv[]) { Eigen::MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col); constexpr double alphay_val = 5e-10; - Eigen::MatrixXd alphay = Eigen::MatrixXd::Constant(row, col, alphay_val); // row,col,value + Eigen::MatrixXd alphay = + Eigen::MatrixXd::Constant(row, col, alphay_val); // row,col,value grid.setAlpha(alphax, alphay); // // ****************** @@ -140,8 +141,7 @@ int main(int argc, char *argv[]) { // // ************************ // set up a simulation environment - Simulation simulation = - Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach + Simulation simulation = Simulation(grid, bc); // grid,boundary // set the timestep of the simulation simulation.setTimestep(360); // timestep diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt deleted file mode 100644 index 91d141d..0000000 --- a/src/CMakeLists.txt +++ /dev/null @@ -1,13 +0,0 @@ -add_library(tug Simulation.cpp FTCS.cpp BTCS.cpp) - -IF(TUG_NAAICE_EXAMPLE) - target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV) -endif() - -target_link_libraries(tug Eigen3::Eigen) - -if(TUG_USE_OPENMP AND OpenMP_CXX_FOUND) - target_link_libraries(tug OpenMP::OpenMP_CXX) -endif() - -target_include_directories(tug PUBLIC ${PROJECT_SOURCE_DIR}/include) diff --git a/src/README.md b/src/README.md deleted file mode 100644 index 82cd165..0000000 --- a/src/README.md +++ /dev/null @@ -1,36 +0,0 @@ -

src-Directory

-This is the src-directory that holds the source code to the implementation of the TUG framework. - -
-src/  
-  ├── CMakeFiles/
-  ├── Boundary.cpp  
-  ├── BTCSv2.cpp
-  ├── CMakeLists.txt  
-  ├── FTCS.cpp
-  ├── Grid.cpp 
-  ├── README.md
-  ├── Simulation.cpp
-  └── TugUtils.hpp 
-
- -**src/** Directory with the source code. - -**CMakeFiles/** Automatically generated directory by CMake. - -**Boundary.cpp** Implementation of Boundary class, that holds the boundary conditions. - -**BTCSv2.cpp** Implementation of BTCS solution to heterogeneous diffusion in 1D and 2D. - -**CMakeLists.txt** CMakeLists for this directory. - -**FTCS.cpp** Implementation of FTCS solution to heterogeneous diffusion in 1D and 2D. - -**Grid.cpp** Implementation of Grid class, that holds all of the concentrations alpha coefficient in x- and y-direction. - -**README.md** This file. - -**Simulation.cpp** Implementation of Simulation class, that holds all of the information for a specific simulation run, as well as the Boundary and Grid objects. - -**TugUtils.hpp** Helper functions for other source files. - diff --git a/src/Schemes.hpp b/src/Schemes.hpp deleted file mode 100644 index 09e6d9e..0000000 --- a/src/Schemes.hpp +++ /dev/null @@ -1,36 +0,0 @@ -/** - * @file BTCSv2.cpp - * @brief Implementation of heterogenous BTCS (backward time-centered space) - * solution of diffusion equation in 1D and 2D space. Internally the - * alternating-direction implicit (ADI) method is used. Version 2, because - * Version 1 was an implementation for the homogeneous BTCS solution. - * - */ - -#ifndef SCHEMES_H_ -#define SCHEMES_H_ - -#include "TugUtils.hpp" - -#include -#include - -namespace tug { - -// entry point; differentiate between 1D and 2D grid -template -extern void FTCS(tug::Grid &grid, tug::Boundary &bc, T timestep, - int &numThreads); - -// entry point for EigenLU solver; differentiate between 1D and 2D grid -template -extern void BTCS_LU(tug::Grid &grid, tug::Boundary &bc, T timestep, - int numThreads); - -// entry point for Thomas algorithm solver; differentiate 1D and 2D grid -template -extern void BTCS_Thomas(tug::Grid &grid, tug::Boundary &bc, T timestep, - int numThreads); -} // namespace tug - -#endif // SCHEMES_H_ diff --git a/src/Simulation.cpp b/src/Simulation.cpp deleted file mode 100644 index ca623d4..0000000 --- a/src/Simulation.cpp +++ /dev/null @@ -1,205 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include "Schemes.hpp" -#include "TugUtils.hpp" - -namespace tug { - -template void Simulation::setTimestep(T timestep) { - if (timestep <= 0) { - throw_invalid_argument("Timestep has to be greater than zero."); - } - - if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { - - const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - // will be 0 if 1D, else ... - const T deltaRowSquare = grid.getDim() != 1 - ? grid.getDeltaRow() * grid.getDeltaRow() - : deltaColSquare; - const T minDeltaSquare = - (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare; - - T maxAlpha = std::numeric_limits::quiet_NaN(); - - // determine maximum alpha - if (grid.getDim() == 2) { - - const T maxAlphaX = grid.getAlphaX().maxCoeff(); - const T maxAlphaY = grid.getAlphaY().maxCoeff(); - maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY; - - } else if (grid.getDim() == 1) { - maxAlpha = grid.getAlpha().maxCoeff(); - - } else { - throw_invalid_argument("Critical error: Undefined number of dimensions!"); - } - - const std::string dim = std::to_string(grid.getDim()) + "D"; - - // Courant-Friedrichs-Lewy condition - T cfl = minDeltaSquare / (4 * maxAlpha); - - // stability equation from Wikipedia; might be useful if applied cfl does - // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha * - // ((1/deltaRowSquare) + (1/deltaColSquare))); - - const std::string &approachPrefix = this->approach_names[approach]; - std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl - << std::endl; - std::cout << approachPrefix << "_" << dim << " :: required dt=" << timestep - << std::endl; - - if (timestep > cfl) { - - this->innerIterations = (int)ceil(timestep / cfl); - this->timestep = timestep / (double)innerIterations; - - std::cerr << "Warning :: Timestep was adjusted, because of stability " - "conditions. Time duration was approximately preserved by " - "adjusting internal number of iterations." - << std::endl; - std::cout << approachPrefix << "_" << dim << " :: Required " - << this->innerIterations - << " inner iterations with dt=" << this->timestep << std::endl; - - } else { - - this->timestep = timestep; - std::cout << approachPrefix << "_" << dim - << " :: No inner iterations required, dt=" << timestep - << std::endl; - } - - } else { - this->timestep = timestep; - } -} - -template void Simulation::run() { - if (this->timestep == -1) { - throw_invalid_argument("Timestep is not set!"); - } - if (this->iterations == -1) { - throw_invalid_argument("Number of iterations are not set!"); - } - - std::string filename; - if (this->console_output > CONSOLE_OUTPUT_OFF) { - printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { - filename = createCSVfile(); - } - - auto begin = std::chrono::high_resolution_clock::now(); - - if (approach == FTCS_APPROACH) { // FTCS case - - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - FTCS(this->grid, this->bc, this->timestep, this->numThreads); - - // if (i % (iterations * innerIterations / 100) == 0) { - // double percentage = (double)i / ((double)iterations * - // (double)innerIterations) * 100; if ((int)percentage % 10 == 0) { - // cout << "Progress: " << percentage << "%" << endl; - // } - // } - } - - } else if (approach == BTCS_APPROACH) { // BTCS case - - if (solver == EIGEN_LU_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads); - } - } else if (solver == THOMAS_ALGORITHM_SOLVER) { - for (int i = 0; i < iterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads); - } - } - - } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case - - constexpr T beta = 0.5; - - // TODO this implementation is very inefficient! - // a separate implementation that sets up a specific tridiagonal matrix for - // Crank-Nicolson would be better - Eigen::MatrixX concentrations; - Eigen::MatrixX concentrationsFTCS; - Eigen::MatrixX concentrationsResult; - for (int i = 0; i < iterations * innerIterations; i++) { - if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { - printConcentrationsConsole(); - } - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationsCSV(filename); - } - - concentrations = grid.getConcentrations(); - 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); - } - } - - auto end = std::chrono::high_resolution_clock::now(); - auto milliseconds = - std::chrono::duration_cast(end - begin); - - if (this->console_output > CONSOLE_OUTPUT_OFF) { - printConcentrationsConsole(); - } - if (this->csv_output > CSV_OUTPUT_OFF) { - printConcentrationsCSV(filename); - } - if (this->time_measure > TIME_MEASURE_OFF) { - const std::string &approachString = this->approach_names[approach]; - const std::string dimString = std::to_string(grid.getDim()) + "D"; - std::cout << approachString << dimString << ":: run() finished in " - << milliseconds.count() << "ms" << std::endl; - } -} - -template void Simulation::setTimestep(double timestep); -template void Simulation::setTimestep(float timestep); - -template void Simulation::run(); -template void Simulation::run(); -} // namespace tug diff --git a/test/testFTCS.cpp b/test/testFTCS.cpp index 209ac8a..8e577ad 100644 --- a/test/testFTCS.cpp +++ b/test/testFTCS.cpp @@ -1,4 +1,4 @@ -#include +#include #include #include diff --git a/test/testSimulation.cpp b/test/testSimulation.cpp index a7600c5..ca6888d 100644 --- a/test/testSimulation.cpp +++ b/test/testSimulation.cpp @@ -12,8 +12,7 @@ using namespace Eigen; using namespace std; using namespace tug; -static Grid64 setupSimulation(APPROACH approach, double timestep, - int iterations) { +Grid64 setupSimulation(double timestep, int iterations) { int row = 11; int col = 11; int domain_row = 10; @@ -45,26 +44,29 @@ static Grid64 setupSimulation(APPROACH approach, double timestep, } grid.setAlpha(alpha, alpha); - // Boundary - Boundary bc = Boundary(grid); - - // Simulation - Simulation sim = Simulation(grid, bc, approach); - // sim.setOutputConsole(CONSOLE_OUTPUT_ON); - sim.setTimestep(timestep); - sim.setIterations(iterations); - sim.run(); - - // RUN return grid; } +constexpr double timestep = 0.001; +constexpr double iterations = 7000; + TEST_CASE("equality to reference matrix with FTCS") { // set string from the header file string test_path = testSimulationCSVDir; MatrixXd reference = CSV2Eigen(test_path); cout << "FTCS Test: " << endl; - Grid grid = setupSimulation(FTCS_APPROACH, 0.001, 7000); + + Grid grid = setupSimulation(timestep, iterations); // Boundary + Boundary bc = Boundary(grid); + + // Simulation + + Simulation sim = Simulation(grid, bc); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + sim.setTimestep(timestep); + sim.setIterations(iterations); + sim.run(); + cout << endl; CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true); } @@ -74,7 +76,17 @@ TEST_CASE("equality to reference matrix with BTCS") { string test_path = testSimulationCSVDir; MatrixXd reference = CSV2Eigen(test_path); cout << "BTCS Test: " << endl; - Grid grid = setupSimulation(BTCS_APPROACH, 1, 7); + + Grid grid = setupSimulation(timestep, iterations); // Boundary + Boundary bc = Boundary(grid); + + // Simulation + Simulation sim = Simulation(grid, bc); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + sim.setTimestep(timestep); + sim.setIterations(iterations); + sim.run(); + cout << endl; CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true); } @@ -84,14 +96,14 @@ TEST_CASE("Initialize environment") { Grid64 grid(rc, rc); Boundary boundary(grid); - CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH)); + CHECK_NOTHROW(Simulation sim(grid, boundary)); } TEST_CASE("Simulation environment") { int rc = 12; Grid64 grid(rc, rc); Boundary boundary(grid); - Simulation sim(grid, boundary, FTCS_APPROACH); + Simulation sim(grid, boundary); SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); } From 77914ea69faffe8cddcdcfd6a8e7c7172566d748 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 17 Oct 2023 12:07:33 +0200 Subject: [PATCH 16/17] fix: include optional output of csv during thomas algorithm fix: marco's benchmark --- include/tug/Core/BTCS.hpp | 23 +++++ naaice/CMakeLists.txt | 2 + naaice/NAAICE_dble_vs_float.cpp | 168 +++++++++++--------------------- 3 files changed, 82 insertions(+), 111 deletions(-) diff --git a/include/tug/Core/BTCS.hpp b/include/tug/Core/BTCS.hpp index cea830a..b844ba4 100644 --- a/include/tug/Core/BTCS.hpp +++ b/include/tug/Core/BTCS.hpp @@ -280,6 +280,29 @@ static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, a_diag[n - 1] = A.coeff(n - 1, n - 2); b_diag[n - 1] = A.coeff(n - 1, n - 1); + // HACK: write CSV to file +#ifdef WRITE_THOMAS_CSV +#include +#include + static std::uint32_t file_index = 0; + std::string file_name = "Thomas_" + std::to_string(file_index++) + ".csv"; + + std::ofstream out_file; + + out_file.open(file_name, std::ofstream::trunc | std::ofstream::out); + + // print header + out_file << "Aa, Ab, Ac, b\n"; + + // iterate through all elements + for (std::size_t i = 0; i < n; i++) { + out_file << a_diag[i] << ", " << b_diag[i] << ", " << c_diag[i] << ", " + << b[i] << "\n"; + } + + out_file.close(); +#endif + // start solving - c_diag and x_vec are overwritten n--; c_diag[0] /= b_diag[0]; diff --git a/naaice/CMakeLists.txt b/naaice/CMakeLists.txt index b0a45d0..209dc91 100644 --- a/naaice/CMakeLists.txt +++ b/naaice/CMakeLists.txt @@ -8,5 +8,7 @@ configure_file(files.hpp.in files.hpp) target_include_directories(naaice PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) target_include_directories(NAAICE_dble_vs_float PRIVATE ${CMAKE_CURRENT_BINARY_DIR}) +target_compile_definitions(naaice PRIVATE WRITE_THOMAS_CSV) + target_link_libraries(naaice PUBLIC tug) target_link_libraries(NAAICE_dble_vs_float PUBLIC tug) diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp index 5be74a4..ac3718c 100644 --- a/naaice/NAAICE_dble_vs_float.cpp +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -7,8 +8,8 @@ #include #include #include +#include #include -#include #include "files.hpp" @@ -96,151 +97,96 @@ Eigen::MatrixXd rmVecTocmMatrix(const std::vector> &vec, return out_mat; } +template int doWork(int ngrid) { + + constexpr T dt = 10; -int DoDble(int ngrid, APPROACH approach) { - - constexpr double dt = 10; - // create a grid - std::cout << "DOUBLE grid: " << ngrid << std::endl; + std::string name; - Grid64 grid64(ngrid, ngrid); + if constexpr (std::is_same_v) { + name = "DOUBLE"; + } else if constexpr (std::is_same_v) { + name = "FLOAT"; + } else { + name = "unknown"; + } + + std::cout << name << " grid: " << ngrid << std::endl; + + Grid grid(ngrid, ngrid); + // Grid64 grid(ngrid, ngrid); // (optional) set the domain, e.g.: - grid64.setDomain(0.1, 0.1); + grid.setDomain(0.1, 0.1); - Eigen::MatrixXd initConc64 = Eigen::MatrixXd::Constant(ngrid, ngrid, 0); + Eigen::MatrixX initConc64 = Eigen::MatrixX::Constant(ngrid, ngrid, 0); initConc64(50, 50) = 1E-6; - grid64.setConcentrations(initConc64); - - const double sum_init64 = initConc64.sum(); + grid.setConcentrations(initConc64); - constexpr double alphax_val = 5e-10; - Eigen::MatrixXd alphax = Eigen::MatrixXd::Constant(ngrid, ngrid, alphax_val); // row,col,value - constexpr double alphay_val = 1e-10; - Eigen::MatrixXd alphay = Eigen::MatrixXd::Constant(ngrid, ngrid, alphay_val); // row,col,value - grid64.setAlpha(alphax, alphay); + const T sum_init64 = initConc64.sum(); + + constexpr T alphax_val = 5e-10; + Eigen::MatrixX alphax = + Eigen::MatrixX::Constant(ngrid, ngrid, alphax_val); // row,col,value + constexpr T alphay_val = 1e-10; + Eigen::MatrixX alphay = + Eigen::MatrixX::Constant(ngrid, ngrid, alphay_val); // row,col,value + grid.setAlpha(alphax, alphay); // create a boundary with constant values - Boundary bc64 = Boundary(grid64); - bc64.setBoundarySideClosed(BC_SIDE_LEFT); - bc64.setBoundarySideClosed(BC_SIDE_RIGHT); - bc64.setBoundarySideClosed(BC_SIDE_TOP); - bc64.setBoundarySideClosed(BC_SIDE_BOTTOM); + Boundary bc = Boundary(grid); + bc.setBoundarySideClosed(BC_SIDE_LEFT); + bc.setBoundarySideClosed(BC_SIDE_RIGHT); + bc.setBoundarySideClosed(BC_SIDE_TOP); + bc.setBoundarySideClosed(BC_SIDE_BOTTOM); // set up a simulation environment - Simulation Sim64 = - Simulation(grid64, bc64, approach); // grid_64,boundary,simulation-approach + Simulation Sim = + Simulation(grid, bc); // grid_64,boundary,simulation-approach + + // Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); - //Sim64.setSolver(THOMAS_ALGORITHM_SOLVER); - // set the timestep of the simulation - Sim64.setTimestep(dt); // timestep + Sim.setTimestep(dt); // timestep // set the number of iterations - Sim64.setIterations(2); + Sim.setIterations(2); // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, // CSV_OUTPUT_VERBOSE] - Sim64.setOutputCSV(CSV_OUTPUT_ON); + Sim.setOutputCSV(CSV_OUTPUT_ON); // set output to the console to 'ON' - Sim64.setOutputConsole(CONSOLE_OUTPUT_OFF); + Sim.setOutputConsole(CONSOLE_OUTPUT_OFF); // // **** RUN SIM64 **** - auto begin64 = std::chrono::high_resolution_clock::now(); + auto begin_t = std::chrono::high_resolution_clock::now(); - Sim64.run(); + Sim.run(); - auto end64 = std::chrono::high_resolution_clock::now(); - auto ms64 = std::chrono::duration_cast(end64 - begin64); + auto end_t = std::chrono::high_resolution_clock::now(); + auto ms_t = + std::chrono::duration_cast(end_t - begin_t); - const double sum_after64 = grid64.getConcentrations().sum(); + const double sum_after64 = grid.getConcentrations().sum(); - std::cout << "Sum of init field: " << std::setprecision(15) << sum_init64 - << "\nSum after 2 iterations: " << sum_after64 - << "\nMilliseconds: " << ms64.count() << std::endl << std::endl; + std::cout << "Sum of init field: " << std::setprecision(15) << sum_init64 + << "\nSum after 2 iterations: " << sum_after64 + << "\nMilliseconds: " << ms_t.count() << std::endl + << std::endl; return 0; - } -int DoFloat(int ngrid, APPROACH approach) { - - constexpr double dt = 10; - - // create a grid - std::cout << "FLOAT grid: " << ngrid << std::endl; - - Grid32 grid32(ngrid, ngrid); - - // (optional) set the domain, e.g.: - grid32.setDomain(0.1, 0.1); - - Eigen::MatrixXf initConc32 = Eigen::MatrixXf::Constant(ngrid, ngrid, 0); - initConc32(50, 50) = 1E-6; - grid32.setConcentrations(initConc32); - - const float sum_init32 = initConc32.sum(); - - constexpr float alphax_val = 5e-10; - Eigen::MatrixXf alphax = Eigen::MatrixXf::Constant(ngrid, ngrid, alphax_val); // row,col,value - constexpr float alphay_val = 1e-10; - Eigen::MatrixXf alphay = Eigen::MatrixXf::Constant(ngrid, ngrid, alphay_val); // row,col,value - grid32.setAlpha(alphax, alphay); - - // create a boundary with constant values - Boundary bc32 = Boundary(grid32); - bc32.setBoundarySideClosed(BC_SIDE_LEFT); - bc32.setBoundarySideClosed(BC_SIDE_RIGHT); - bc32.setBoundarySideClosed(BC_SIDE_TOP); - bc32.setBoundarySideClosed(BC_SIDE_BOTTOM); - - // set up a simulation environment - Simulation Sim32 = - Simulation(grid32, bc32, approach); // grid_32,boundary,simulation-approach - - // Sim32.setSolver(THOMAS_ALGORITHM_SOLVER); - - // set the timestep of the simulation - Sim32.setTimestep(dt); // timestep - - // set the number of iterations - Sim32.setIterations(2); - - // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, - // CSV_OUTPUT_VERBOSE] - Sim32.setOutputCSV(CSV_OUTPUT_ON); - - // set output to the console to 'ON' - Sim32.setOutputConsole(CONSOLE_OUTPUT_OFF); - - // // **** RUN SIM32 **** - auto begin32 = std::chrono::high_resolution_clock::now(); - - Sim32.run(); - - auto end32 = std::chrono::high_resolution_clock::now(); - auto ms32 = std::chrono::duration_cast(end32 - begin32); - - const float sum_after32 = grid32.getConcentrations().sum(); - - std::cout << "Sum of init field: " << std::setprecision(15) << sum_init32 - << "\nSum after 2 iteration: " << sum_after32 - << "\nMilliseconds: " << ms32.count() << std::endl << std::endl; - return 0; - -} - - int main(int argc, char *argv[]) { int n[] = {101, 201, 501, 1001, 2001}; for (int i = 0; i < std::size(n); i++) { - DoFloat(n[i], BTCS_APPROACH); - DoDble(n[i], BTCS_APPROACH); - DoFloat(n[i], FTCS_APPROACH); - DoDble(n[i], FTCS_APPROACH); + doWork(n[i]); + doWork(n[i]); + doWork(n[i]); + doWork(n[i]); } return 0; From 5a39f5377ea45f2e1d75006f05440e24a75b6085 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 19 Oct 2023 13:09:22 +0200 Subject: [PATCH 17/17] doc: update example pages --- docs_sphinx/examples.rst | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/docs_sphinx/examples.rst b/docs_sphinx/examples.rst index e9736f2..9de5a8f 100644 --- a/docs_sphinx/examples.rst +++ b/docs_sphinx/examples.rst @@ -11,18 +11,21 @@ Two dimensional grid with constant boundaries and FTCS method ------------------------------------------------------------- **Initialization of the grid** -For example, the initalization of a grid with 20 by 20 cells and a domain size (physical extent of the grid) of -also 20 by 20 length units can be done as follows. The setting of the domain is optional here and is set to the -same size as the number of cells in the standard case. As seen in the code, the cells of the grid are set to an -initial value of 0 and only in the upper left corner (0,0) the starting concentration is set to the value 20. +For example, the initalization of a grid with 20 by 20 cells using double values +and a domain size (physical extent of the grid) of also 20 by 20 length units +can be done as follows. The setting of the domain is optional here and is set to +the same size as the number of cells in the standard case. As seen in the code, +the cells of the grid are set to an initial value of 0 and only in the upper +left corner (0,0) the starting concentration is set to the value 20. .. code-block:: cpp int row = 20 int col = 20; - Grid grid = Grid(row,col); + Grid grid(row,col); grid.setDomain(row, col); MatrixXd concentrations = MatrixXd::Constant(row,col,0); + // or MatrixX concentrations = MatrixX::Constant(row,col,0); concentrations(0,0) = 20; grid.setConcentrations(concentrations); @@ -40,14 +43,17 @@ of the grid are set as constant edges with a concentration of 0. bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); **Setting of the simulation parameters and simulation start** -In the last block, a simulation class is created and the objects of the grid and the boundary conditions are passed. The solution -method is also specified (either FCTS or BTCS). Furthermore, the desired time step and the number of iterations are set. The penultimate -parameter specifies the output of the simulated results in a CSV file. In the present case, the result of each iteration step is written -one below the other into the corresponding CSV file. + +In the last block, a simulation class is created and the objects of the grid and +the boundary conditions are passed. The solution method is also specified +(either FCTS or BTCS). Furthermore, the desired time step and the number of +iterations are set. The penultimate parameter specifies the output of the +simulated results in a CSV file. In the present case, the result of each +iteration step is written one below the other into the corresponding CSV file. .. code-block:: cpp - Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); + Simulation simulation(grid, bc); simulation.setTimestep(0.1); simulation.setIterations(1000); simulation.setOutputCSV(CSV_OUTPUT_VERBOSE); @@ -59,4 +65,4 @@ one below the other into the corresponding CSV file. Setting special boundary conditions on individual cells -------------------------------------------------------- \ No newline at end of file +-------------------------------------------------------