diff --git a/include/BoundaryCondition.hpp b/include/BoundaryCondition.hpp index 0b80d18..6dc8dae 100644 --- a/include/BoundaryCondition.hpp +++ b/include/BoundaryCondition.hpp @@ -166,7 +166,7 @@ private: uint8_t dim; - uint32_t sizes[2]; + std::array sizes; uint32_t maxsize; uint32_t maxindex; @@ -260,6 +260,6 @@ public: } }; -} // namespace boundary_condition +} // namespace bc } // namespace tug #endif // BOUNDARYCONDITION_H_ diff --git a/include/Diffusion.hpp b/include/Diffusion.hpp index 0f03088..3a93722 100644 --- a/include/Diffusion.hpp +++ b/include/Diffusion.hpp @@ -12,14 +12,17 @@ namespace tug { namespace diffusion { +constexpr uint8_t MAX_ARR_SIZE = 3; + /** * Defines grid dimensions and boundary conditions. */ -typedef struct { - uint32_t - grid_cells[3]; /**< Count of grid cells in each of the 3 directions.*/ - double domain_size[3]; /**< Domain sizes in each of the 3 directions.*/ - bc::BoundaryCondition *bc; /**< Boundary conditions for the grid.*/ +typedef struct tug_grid_s { + std::array + grid_cells; /**< Count of grid cells in each of the 3 directions.*/ + std::array + domain_size; /**< Domain sizes in each of the 3 directions.*/ + bc::BoundaryCondition *bc = NULL; /**< Boundary conditions for the grid.*/ } TugGrid; /** @@ -28,7 +31,8 @@ typedef struct { */ typedef struct tug_input_s { double time_step; /**< Time step which should be simulated by diffusion.*/ - Eigen::VectorXd (*solver)(Eigen::SparseMatrix, Eigen::VectorXd) = + Eigen::VectorXd (*solver)(const Eigen::SparseMatrix &, + const Eigen::VectorXd &) = tug::solver::ThomasAlgorithm; /**< Solver function to use.*/ TugGrid grid; /**< Grid specification.*/ @@ -86,8 +90,9 @@ typedef struct tug_input_s { * \param f_in Pointer to function which takes a sparse matrix and a vector as * input and returns another vector. */ - void setSolverFunction(Eigen::VectorXd (*f_in)(Eigen::SparseMatrix, - Eigen::VectorXd)) { + void + setSolverFunction(Eigen::VectorXd (*f_in)(const Eigen::SparseMatrix &, + const Eigen::VectorXd &)) { solver = f_in; } } TugInput; diff --git a/include/Solver.hpp b/include/Solver.hpp index 6d20493..96f8d55 100644 --- a/include/Solver.hpp +++ b/include/Solver.hpp @@ -18,8 +18,8 @@ namespace solver { * * \return Solution represented as vector. */ -auto EigenLU(const Eigen::SparseMatrix A_matrix, - const Eigen::VectorXd b_vector) -> Eigen::VectorXd; +auto EigenLU(const Eigen::SparseMatrix &A_matrix, + const Eigen::VectorXd &b_vector) -> Eigen::VectorXd; /** * Solving linear equation system with brutal implementation of the Thomas @@ -32,8 +32,8 @@ auto EigenLU(const Eigen::SparseMatrix A_matrix, * * \return Solution represented as vector. */ -auto ThomasAlgorithm(const Eigen::SparseMatrix A_matrix, - const Eigen::VectorXd b_vector) -> Eigen::VectorXd; +auto ThomasAlgorithm(const Eigen::SparseMatrix &A_matrix, + const Eigen::VectorXd &b_vector) -> Eigen::VectorXd; } // namespace solver } // namespace tug diff --git a/src/BTCS.cpp b/src/BTCS.cpp index 0a959e4..525ba4b 100644 --- a/src/BTCS.cpp +++ b/src/BTCS.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -7,6 +8,7 @@ #include #include #include +#include #include "BoundaryCondition.hpp" #include "TugUtils.hpp" @@ -17,12 +19,16 @@ #define omp_get_thread_num() 0 #endif -#define init_delta(in, out, dim) \ - ({ \ - for (uint8_t i = 0; i < dim; i++) { \ - out[i] = (double)in.domain_size[i] / in.grid_cells[i]; \ - } \ - }) +inline auto +init_delta(const std::array &domain_size, + const std::array &grid_cells, + const uint8_t dim) -> std::vector { + std::vector out(dim); + for (uint8_t i = 0; i < dim; i++) { + out[i] = (double)(domain_size.at(i) / grid_cells.at(i)); + } + return out; +} namespace { enum { GRID_1D = 1, GRID_2D, GRID_3D }; @@ -30,10 +36,10 @@ enum { GRID_1D = 1, GRID_2D, GRID_3D }; constexpr int BTCS_MAX_DEP_PER_CELL = 3; constexpr int BTCS_2D_DT_SIZE = 2; -typedef Eigen::Matrix - DMatrixRowMajor; -typedef Eigen::Matrix - DVectorRowMajor; +using DMatrixRowMajor = + Eigen::Matrix; +using DVectorRowMajor = + Eigen::Matrix; inline auto getBCFromFlux(tug::bc::boundary_condition bc, double neighbor_c, double neighbor_alpha) -> double { @@ -192,9 +198,10 @@ auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha, for (int j = 0; j < size; j++) { if (bc_inner[j].type != tug::bc::BC_UNSET) { - if (bc_inner[j].type != tug::bc::BC_TYPE_CONSTANT) + if (bc_inner[j].type != tug::bc::BC_TYPE_CONSTANT) { throw_invalid_argument("Inner boundary conditions with other type than " "BC_TYPE_CONSTANT are currently not supported."); + } b_vector[j + 1] = bc_inner[j].value; continue; } @@ -215,10 +222,10 @@ auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha, auto setupBTCSAndSolve( DVectorRowMajor &c, const tug::bc::bc_tuple bc_ghost, - const tug::bc::bc_vec bc_inner, const DVectorRowMajor &alpha, double dx, + const tug::bc::bc_vec &bc_inner, const DVectorRowMajor &alpha, double dx, double time_step, int size, const DVectorRowMajor &d_ortho, - Eigen::VectorXd (*solver)(Eigen::SparseMatrix, Eigen::VectorXd)) - -> DVectorRowMajor { + Eigen::VectorXd (*solver)(const Eigen::SparseMatrix &, + const Eigen::VectorXd &)) -> DVectorRowMajor { const Eigen::SparseMatrix A_matrix = fillMatrixFromRow(alpha, bc_inner, size, dx, time_step); @@ -241,15 +248,18 @@ auto tug::diffusion::BTCS_1D(const tug::diffusion::TugInput &input_param, auto start = time_marker(); - const tug::bc::BoundaryCondition bc = *input_param.grid.bc; - double deltas[GRID_1D]; - - init_delta(input_param.grid, deltas, GRID_1D); - uint32_t size = input_param.grid.grid_cells[0]; + + auto deltas = init_delta(input_param.grid.domain_size, + input_param.grid.grid_cells, GRID_1D); double dx = deltas[0]; + double time_step = input_param.time_step; + const tug::bc::BoundaryCondition bc = + (input_param.grid.bc != nullptr ? *input_param.grid.bc + : tug::bc::BoundaryCondition(size)); + Eigen::Map c_in(field, size); Eigen::Map alpha_in(alpha, size); @@ -271,17 +281,21 @@ auto tug::diffusion::ADI_2D(const tug::diffusion::TugInput &input_param, auto start = time_marker(); - tug::bc::BoundaryCondition bc = *input_param.grid.bc; - double deltas[GRID_2D]; - - init_delta(input_param.grid, deltas, GRID_2D); - uint32_t n_cols = input_param.grid.grid_cells[0]; uint32_t n_rows = input_param.grid.grid_cells[1]; + + auto deltas = init_delta(input_param.grid.domain_size, + input_param.grid.grid_cells, GRID_2D); double dx = deltas[0]; double dy = deltas[1]; + double local_dt = input_param.time_step / BTCS_2D_DT_SIZE; + tug::bc::BoundaryCondition bc = + (input_param.grid.bc != nullptr + ? *input_param.grid.bc + : tug::bc::BoundaryCondition(n_cols, n_rows)); + Eigen::Map c_in(field, n_rows, n_cols); Eigen::Map alpha_in(alpha, n_rows, n_cols); diff --git a/src/Solver.cpp b/src/Solver.cpp index 4cf6f1c..a3514bf 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -4,8 +4,8 @@ #include #include -auto tug::solver::EigenLU(const Eigen::SparseMatrix A_matrix, - const Eigen::VectorXd b_vector) -> Eigen::VectorXd { +auto tug::solver::EigenLU(const Eigen::SparseMatrix &A_matrix, + const Eigen::VectorXd &b_vector) -> Eigen::VectorXd { Eigen::SparseLU, Eigen::COLAMDOrdering> solver; @@ -18,8 +18,8 @@ auto tug::solver::EigenLU(const Eigen::SparseMatrix A_matrix, return solver.solve(b_vector); } -auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix A_matrix, - const Eigen::VectorXd b_vector) +auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix &A_matrix, + const Eigen::VectorXd &b_vector) -> Eigen::VectorXd { uint32_t n = b_vector.size(); diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 1ecb321..a940d57 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -15,13 +15,12 @@ using namespace tug::diffusion; static std::vector alpha(N *M, 1e-3); -static TugInput setupDiffu(BoundaryCondition &bc) { +static TugInput setupDiffu() { TugInput diffu; diffu.setTimestep(1); diffu.setGridCellN(N, M); diffu.setDomainSize(N, M); - diffu.setBoundaryCondition(bc); return diffu; } @@ -33,7 +32,7 @@ TEST_CASE("closed boundaries - 1 concentration to 1 - rest 0") { BoundaryCondition bc(N, M); - TugInput diffu = setupDiffu(bc); + TugInput diffu = setupDiffu(); uint32_t iterations = 1000; double sum = 0; @@ -68,7 +67,8 @@ TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") { bc.setSide(BC_SIDE_TOP, input); bc.setSide(BC_SIDE_BOTTOM, input); - TugInput diffu = setupDiffu(bc); + TugInput diffu = setupDiffu(); + diffu.setBoundaryCondition(bc); uint32_t max_iterations = 20000; bool reached = false; @@ -106,7 +106,8 @@ TEST_CASE( bc.setSide(BC_SIDE_TOP, top); bc.setSide(BC_SIDE_BOTTOM, bottom); - TugInput diffu = setupDiffu(bc); + TugInput diffu = setupDiffu(); + diffu.setBoundaryCondition(bc); uint32_t max_iterations = 100; @@ -136,7 +137,8 @@ TEST_CASE("2D closed boundaries, 1 constant cell in the middle") { field[MID] = val; bc(BC_INNER, MID) = {BC_TYPE_CONSTANT, val}; - TugInput diffu = setupDiffu(bc); + TugInput diffu = setupDiffu(); + diffu.setBoundaryCondition(bc); uint32_t max_iterations = 100;