diff --git a/include/grid/BoundaryCondition.hpp b/include/BoundaryCondition.hpp similarity index 99% rename from include/grid/BoundaryCondition.hpp rename to include/BoundaryCondition.hpp index adb8e28..0b80d18 100644 --- a/include/grid/BoundaryCondition.hpp +++ b/include/BoundaryCondition.hpp @@ -10,7 +10,7 @@ typedef uint8_t bctype; namespace tug { -namespace boundary_condition { +namespace bc { enum { BC_TYPE_CLOSED, /**< Defines a closed/Neumann boundary condition. */ diff --git a/include/Diffusion.hpp b/include/Diffusion.hpp new file mode 100644 index 0000000..983293d --- /dev/null +++ b/include/Diffusion.hpp @@ -0,0 +1,72 @@ +#ifndef DIFFUSION_H_ +#define DIFFUSION_H_ + +#include "BoundaryCondition.hpp" +#include +#include +#include +#include +#include + +namespace tug { +namespace diffusion { + +/** + * Defines grid dimensions and boundary conditions. + */ +typedef struct { + uint32_t grid_cells[3]; + double domain_size[3]; + bc::BoundaryCondition *bc; +} TugGrid; + +/** + * Besides containing the grid structure it holds also information about the + * desired time step to simulate and which solver to use. + */ +typedef struct { + double time_step; + Eigen::VectorXd (*solver)(Eigen::SparseMatrix, Eigen::VectorXd); + TugGrid grid; +} TugInput; + +/** + * Solving 1D diffusion problems with Backward Time Centred Space scheme. + * + * \param input_param Object with information for the simulation e.g. grid + * dimensions, time step ... + * + * \param field Pointer to continious memory holding the concentrations for each + * grid cell of the grid. It doesn't matter if stored in row (most likely) or + * column major. + * + * \param alpha Pointer to continious memory holding the alpha for each grid + * cell of the grid. (NOTE: only constant alphas are supported currently) + * + * \return Runtime in seconds + */ +auto BTCS_1D(const TugInput &input_param, double *field, const double *alpha) + -> double; + +/** + * Solving 2D diffusion problems with Alternating-direction implicit method. + * + * \param input_param Object with information for the simulation e.g. grid + * dimensions, time step ... + * + * \param field Pointer to continious memory holding the concentrations for each + * grid cell of the grid. It doesn't matter if stored in row (most likely) or + * column major. + * + * \param alpha Pointer to continious memory holding the alpha for each grid + * cell of the grid. (NOTE: only constant alphas are supported currently) + * + * \return Runtime in seconds + */ +auto ADI_2D(const TugInput &input_param, double *field, const double *alpha) + -> double; + +} // namespace diffusion +} // namespace tug + +#endif // DIFFUSION_H_ diff --git a/include/Solver.hpp b/include/Solver.hpp new file mode 100644 index 0000000..6d20493 --- /dev/null +++ b/include/Solver.hpp @@ -0,0 +1,41 @@ +#ifndef SOLVER_H_ +#define SOLVER_H_ + +#include +#include + +namespace tug { +namespace solver { + +/** + * Solving linear equation system with LU-decomposition implemented in Eigen + * library. + * + * \param A_matrix The A matrix represented as a sparse matrix using Eigen + * library. + * + * \param b_vector Right hand side vector of the linear equation system. + * + * \return Solution represented as vector. + */ +auto EigenLU(const Eigen::SparseMatrix A_matrix, + const Eigen::VectorXd b_vector) -> Eigen::VectorXd; + +/** + * Solving linear equation system with brutal implementation of the Thomas + * algorithm. + * + * \param A_matrix The A matrix represented as a sparse matrix using Eigen + * library. + * + * \param b_vector Right hand side vector of the linear equation system. + * + * \return Solution represented as vector. + */ +auto ThomasAlgorithm(const Eigen::SparseMatrix A_matrix, + const Eigen::VectorXd b_vector) -> Eigen::VectorXd; + +} // namespace solver +} // namespace tug + +#endif // SOLVER_H_ diff --git a/include/diffusion/BTCSDiffusion.hpp b/include/diffusion/BTCSDiffusion.hpp deleted file mode 100644 index 3132f85..0000000 --- a/include/diffusion/BTCSDiffusion.hpp +++ /dev/null @@ -1,168 +0,0 @@ -#ifndef BTCSDIFFUSION_H_ -#define BTCSDIFFUSION_H_ - -#include "grid/BoundaryCondition.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace tug { -namespace diffusion { - -/*! - * Class implementing a solution for a 1/2/3D diffusion equation using backward - * euler. - */ -class BTCSDiffusion { - -public: - /*! - * Creates a diffusion module. - * - * \param dim Number of dimensions. Should not be greater than 3 and not less - * than 1. - */ - BTCSDiffusion(unsigned int dim); - - /*! - * Define the grid in x direction. - * - * \param domain_size Size of the domain in x direction. - * \param n_grid_cells Number of grid cells in x direction the domain is - * divided to. - */ - void setXDimensions(double domain_size, unsigned int n_grid_cells); - - /*! - * Define the grid in y direction. - * - * Throws an error if the module wasn't initialized at least as a 2D model. - * - * \param domain_size Size of the domain in y direction. - * \param n_grid_cells Number of grid cells in y direction the domain is - * divided to. - */ - void setYDimensions(double domain_size, unsigned int n_grid_cells); - - /*! - * Define the grid in z direction. - * - * Throws an error if the module wasn't initialized at least as a 3D model. - * - * \param domain_size Size of the domain in z direction. - * \param n_grid_cells Number of grid cells in z direction the domain is - * divided to. - */ - void setZDimensions(double domain_size, unsigned int n_grid_cells); - - /*! - * Returns the number of grid cells in x direction. - */ - auto getXGridCellsN() -> unsigned int; - /*! - * Returns the number of grid cells in y direction. - */ - auto getYGridCellsN() -> unsigned int; - /*! - * Returns the number of grid cells in z direction. - */ - auto getZGridCellsN() -> unsigned int; - - /*! - * Returns the domain size in x direction. - */ - auto getXDomainSize() -> double; - /*! - * Returns the domain size in y direction. - */ - auto getYDomainSize() -> double; - /*! - * Returns the domain size in z direction. - */ - auto getZDomainSize() -> double; - - /*! - * With given ghost zones simulate diffusion. Only 1D allowed at this moment. - * - * \param c Pointer to continious memory describing the current concentration - * state of each grid cell. - * \param alpha Pointer to memory area of diffusion coefficients for each grid - * element. - * \param bc Instance of boundary condition class with. - * - * \return Time in seconds [s] used to simulate one iteration. - */ - auto simulate(double *c, double *alpha, - const tug::boundary_condition::BoundaryCondition &bc) - -> double; - - /*! - * Set the timestep of the simulation - * - * \param time_step Time step (in seconds ???) - */ - void setTimestep(double time_step); - -private: - typedef Eigen::Matrix - DMatrixRowMajor; - typedef Eigen::Matrix - DVectorRowMajor; - - static void simulate_base(DVectorRowMajor &c, - const tug::boundary_condition::bc_tuple &bc_ghosts, - const tug::boundary_condition::bc_vec &bc_inner, - const DVectorRowMajor &alpha, double dx, - double time_step, int size, - const DVectorRowMajor &d_ortho); - - void simulate1D(Eigen::Map &c, - Eigen::Map &alpha, - const tug::boundary_condition::BoundaryCondition &bc); - - void simulate2D(Eigen::Map &c, - Eigen::Map &alpha, - const tug::boundary_condition::BoundaryCondition &bc); - - static auto - calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha, - const tug::boundary_condition::BoundaryCondition &bc, - bool transposed, double time_step, double dx) -> DMatrixRowMajor; - - static void fillMatrixFromRow(Eigen::SparseMatrix &A_matrix, - const DVectorRowMajor &alpha, - const tug::boundary_condition::bc_vec &bc_inner, - int size, double dx, double time_step); - - static void - fillVectorFromRow(Eigen::VectorXd &b_vector, const DVectorRowMajor &c, - const DVectorRowMajor &alpha, - const tug::boundary_condition::bc_tuple &bc_ghosts, - const tug::boundary_condition::bc_vec &bc_inner, - const DVectorRowMajor &d_ortho, int size, double dx, - double time_step); - void simulate3D(std::vector &c); - - inline static auto - getBCFromFlux(tug::boundary_condition::boundary_condition bc, - double neighbor_c, double neighbor_alpha) -> double; - - void updateInternals(); - - double time_step; - unsigned int grid_dim; - - std::vector grid_cells; - std::vector domain_size; - std::vector deltas; -}; -} // namespace diffusion -} // namespace tug -#endif // BTCSDIFFUSION_H_ diff --git a/src/BTCS.cpp b/src/BTCS.cpp new file mode 100644 index 0000000..0a959e4 --- /dev/null +++ b/src/BTCS.cpp @@ -0,0 +1,315 @@ +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "BoundaryCondition.hpp" +#include "TugUtils.hpp" + +#ifdef _OPENMP +#include +#else +#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]; \ + } \ + }) + +namespace { +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; + +inline auto getBCFromFlux(tug::bc::boundary_condition bc, double neighbor_c, + double neighbor_alpha) -> double { + double val = 0; + + if (bc.type == tug::bc::BC_TYPE_CLOSED) { + val = neighbor_c; + } else if (bc.type == tug::bc::BC_TYPE_FLUX) { + // TODO + // val = bc[index].value; + } else { + // TODO: implement error handling here. Type was set to wrong value. + } + + return val; +} + +auto calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha, + const tug::bc::BoundaryCondition &bc, bool transposed, + double time_step, double dx) -> DMatrixRowMajor { + + uint8_t upper = (transposed ? tug::bc::BC_SIDE_LEFT : tug::bc::BC_SIDE_TOP); + uint8_t lower = + (transposed ? tug::bc::BC_SIDE_RIGHT : tug::bc::BC_SIDE_BOTTOM); + + int n_rows = c.rows(); + int n_cols = c.cols(); + + DMatrixRowMajor d_ortho(n_rows, n_cols); + + std::array y_values{}; + + // first, iterate over first row + for (int j = 0; j < n_cols; j++) { + tug::bc::boundary_condition tmp_bc = bc(upper, j); + double sy = (time_step * alpha(0, j)) / (dx * dx); + + y_values[0] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT + ? tmp_bc.value + : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j))); + y_values[1] = c(0, j); + y_values[2] = c(1, j); + + d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]); + } + + // then iterate over inlet +#pragma omp parallel for private(y_values) schedule(dynamic) + for (int i = 1; i < n_rows - 1; i++) { + for (int j = 0; j < n_cols; j++) { + double sy = (time_step * alpha(i, j)) / (dx * dx); + + y_values[0] = c(i - 1, j); + y_values[1] = c(i, j); + y_values[2] = c(i + 1, j); + + d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]); + } + } + + int end = n_rows - 1; + + // and finally over last row + for (int j = 0; j < n_cols; j++) { + tug::bc::boundary_condition tmp_bc = bc(lower, j); + double sy = (time_step * alpha(end, j)) / (dx * dx); + + y_values[0] = c(end - 1, j); + y_values[1] = c(end, j); + y_values[2] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT + ? tmp_bc.value + : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j))); + + d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]); + } + + return d_ortho; +} + +auto fillMatrixFromRow(const DVectorRowMajor &alpha, + const tug::bc::bc_vec &bc_inner, int size, double dx, + double time_step) -> Eigen::SparseMatrix { + + Eigen::SparseMatrix A_matrix(size + 2, size + 2); + + A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL)); + + double sx = 0; + + int A_size = A_matrix.cols(); + + A_matrix.insert(0, 0) = 1; + + if (bc_inner[0].type != tug::bc::BC_UNSET) { + if (bc_inner[0].type != tug::bc::BC_TYPE_CONSTANT) { + throw_invalid_argument("Inner boundary conditions with other type than " + "BC_TYPE_CONSTANT are currently not supported."); + } + A_matrix.insert(1, 1) = 1; + } else { + sx = (alpha[0] * time_step) / (dx * dx); + A_matrix.insert(1, 1) = -1. - 3. * sx; + A_matrix.insert(1, 0) = 2. * sx; + A_matrix.insert(1, 2) = sx; + } + + for (int j = 2, k = j - 1; k < size - 1; j++, k++) { + if (bc_inner[k].type != tug::bc::BC_UNSET) { + if (bc_inner[k].type != tug::bc::BC_TYPE_CONSTANT) { + throw_invalid_argument("Inner boundary conditions with other type than " + "BC_TYPE_CONSTANT are currently not supported."); + } + A_matrix.insert(j, j) = 1; + continue; + } + sx = (alpha[k] * time_step) / (dx * dx); + + A_matrix.insert(j, j) = -1. - 2. * sx; + A_matrix.insert(j, (j - 1)) = sx; + A_matrix.insert(j, (j + 1)) = sx; + } + + if (bc_inner[size - 1].type != tug::bc::BC_UNSET) { + if (bc_inner[size - 1].type != tug::bc::BC_TYPE_CONSTANT) { + throw_invalid_argument("Inner boundary conditions with other type than " + "BC_TYPE_CONSTANT are currently not supported."); + } + A_matrix.insert(A_size - 2, A_size - 2) = 1; + } else { + sx = (alpha[size - 1] * time_step) / (dx * dx); + A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx; + A_matrix.insert(A_size - 2, A_size - 3) = sx; + A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx; + } + + A_matrix.insert(A_size - 1, A_size - 1) = 1; + + return A_matrix; +} + +auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha, + const tug::bc::bc_tuple &bc, + const tug::bc::bc_vec &bc_inner, + const DVectorRowMajor &d_ortho, int size, double dx, + double time_step) -> Eigen::VectorXd { + + Eigen::VectorXd b_vector(size + 2); + + tug::bc::boundary_condition left = bc[0]; + tug::bc::boundary_condition right = bc[1]; + + bool left_constant = (left.type == tug::bc::BC_TYPE_CONSTANT); + bool right_constant = (right.type == tug::bc::BC_TYPE_CONSTANT); + + int b_size = b_vector.size(); + + 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) + 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; + } + b_vector[j + 1] = -c[j] + d_ortho[j]; + } + + // this is not correct currently.We will fix this when we are able to define + // FLUX boundary conditions + b_vector[0] = + (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); + + b_vector[b_size - 1] = + (right_constant ? right.value + : getBCFromFlux(right, c[size - 1], alpha[size - 1])); + + return b_vector; +} + +auto setupBTCSAndSolve( + DVectorRowMajor &c, const tug::bc::bc_tuple bc_ghost, + 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 { + + const Eigen::SparseMatrix A_matrix = + fillMatrixFromRow(alpha, bc_inner, size, dx, time_step); + + const Eigen::VectorXd b_vector = fillVectorFromRow( + c, alpha, bc_ghost, bc_inner, d_ortho, size, dx, time_step); + + // solving of the LEQ + Eigen::VectorXd x_vector = solver(A_matrix, b_vector); + + DVectorRowMajor out_vector = x_vector.segment(1, size); + + return out_vector; +} + +} // namespace + // +auto tug::diffusion::BTCS_1D(const tug::diffusion::TugInput &input_param, + double *field, const double *alpha) -> double { + + 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]; + double dx = deltas[0]; + double time_step = input_param.time_step; + + Eigen::Map c_in(field, size); + Eigen::Map alpha_in(alpha, size); + + DVectorRowMajor input_field = c_in.row(0); + + DVectorRowMajor output = setupBTCSAndSolve( + input_field, bc.row_boundary(0), bc.getInnerRow(0), alpha_in, dx, + time_step, size, Eigen::VectorXd::Constant(size, 0), input_param.solver); + + c_in.row(0) << output; + + auto end = time_marker(); + + return diff_time(start, end); +} + +auto tug::diffusion::ADI_2D(const tug::diffusion::TugInput &input_param, + double *field, const double *alpha) -> double { + + 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]; + double dx = deltas[0]; + double dy = deltas[1]; + double local_dt = input_param.time_step / BTCS_2D_DT_SIZE; + + Eigen::Map c_in(field, n_rows, n_cols); + Eigen::Map alpha_in(alpha, n_rows, n_cols); + + DMatrixRowMajor d_ortho = + calc_d_ortho(c_in, alpha_in, bc, false, local_dt, dx); + +#pragma omp parallel for schedule(dynamic) + for (int i = 0; i < n_rows; i++) { + DVectorRowMajor input_field = c_in.row(i); + DVectorRowMajor output = setupBTCSAndSolve( + input_field, bc.row_boundary(i), bc.getInnerRow(i), alpha_in.row(i), dx, + local_dt, n_cols, d_ortho.row(i), input_param.solver); + c_in.row(i) << output; + } + + d_ortho = calc_d_ortho(c_in.transpose(), alpha_in.transpose(), bc, true, + local_dt, dy); + +#pragma omp parallel for schedule(dynamic) + for (int i = 0; i < n_cols; i++) { + DVectorRowMajor input_field = c_in.col(i); + DVectorRowMajor output = setupBTCSAndSolve( + input_field, bc.col_boundary(i), bc.getInnerCol(i), alpha_in.col(i), dy, + local_dt, n_rows, d_ortho.row(i), input_param.solver); + c_in.col(i) << output.transpose(); + } + + auto end = time_marker(); + + return diff_time(start, end); +} diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp deleted file mode 100644 index c27db1a..0000000 --- a/src/BTCSDiffusion.cpp +++ /dev/null @@ -1,381 +0,0 @@ -#include "diffusion/BTCSDiffusion.hpp" -#include "BTCSUtils.hpp" -#include "grid/BoundaryCondition.hpp" - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#ifdef _OPENMP -#include -#else -#define omp_get_thread_num() 0 -#endif - -constexpr double DOUBLE_MACHINE_EPSILON = 1.93e-34; - -constexpr int BTCS_MAX_DEP_PER_CELL = 3; -constexpr int BTCS_2D_DT_SIZE = 2; -tug::diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { - - grid_cells.resize(dim, 1); - domain_size.resize(dim, 1); - deltas.resize(dim, 1); - - this->time_step = 0; -} - -void tug::diffusion::BTCSDiffusion::setXDimensions(double domain_size, - unsigned int n_grid_cells) { - this->domain_size[0] = domain_size; - this->grid_cells[0] = n_grid_cells; - - updateInternals(); -} - -void tug::diffusion::BTCSDiffusion::setYDimensions(double domain_size, - unsigned int n_grid_cells) { - this->domain_size[1] = domain_size; - this->grid_cells[1] = n_grid_cells; - - updateInternals(); -} - -void tug::diffusion::BTCSDiffusion::setZDimensions(double domain_size, - unsigned int n_grid_cells) { - this->domain_size[2] = domain_size; - this->grid_cells[2] = n_grid_cells; - - updateInternals(); -} - -auto tug::diffusion::BTCSDiffusion::getXGridCellsN() -> unsigned int { - return this->grid_cells[0]; -} -auto tug::diffusion::BTCSDiffusion::getYGridCellsN() -> unsigned int { - return this->grid_cells[1]; -} -auto tug::diffusion::BTCSDiffusion::getZGridCellsN() -> unsigned int { - return this->grid_cells[2]; -} -auto tug::diffusion::BTCSDiffusion::getXDomainSize() -> double { - return this->domain_size[0]; -} -auto tug::diffusion::BTCSDiffusion::getYDomainSize() -> double { - return this->domain_size[1]; -} -auto tug::diffusion::BTCSDiffusion::getZDomainSize() -> double { - return this->domain_size[2]; -} - -void tug::diffusion::BTCSDiffusion::updateInternals() { - for (int i = 0; i < grid_dim; i++) { - deltas[i] = (double)domain_size[i] / grid_cells[i]; - } -} -void tug::diffusion::BTCSDiffusion::simulate_base( - DVectorRowMajor &c, const tug::boundary_condition::bc_tuple &bc_ghosts, - const tug::boundary_condition::bc_vec &bc_inner, - const DVectorRowMajor &alpha, double dx, double time_step, int size, - const DVectorRowMajor &d_ortho) { - - Eigen::SparseMatrix A_matrix; - Eigen::VectorXd b_vector; - Eigen::VectorXd x_vector; - - A_matrix.resize(size + 2, size + 2); - A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL)); - - b_vector.resize(size + 2); - x_vector.resize(size + 2); - - fillMatrixFromRow(A_matrix, alpha, bc_inner, size, dx, time_step); - fillVectorFromRow(b_vector, c, alpha, bc_ghosts, bc_inner, d_ortho, size, dx, - time_step); - - // start to solve - Eigen::SparseLU, Eigen::COLAMDOrdering> - solver; - solver.analyzePattern(A_matrix); - - solver.factorize(A_matrix); - - x_vector = solver.solve(b_vector); - - c = x_vector.segment(1, size); -} - -void tug::diffusion::BTCSDiffusion::simulate1D( - Eigen::Map &c, Eigen::Map &alpha, - const tug::boundary_condition::BoundaryCondition &bc) { - - int size = this->grid_cells[0]; - double dx = this->deltas[0]; - double time_step = this->time_step; - - DVectorRowMajor input_field = c.row(0); - - simulate_base(input_field, bc.row_boundary(0), bc.getInnerRow(0), alpha, dx, - time_step, size, Eigen::VectorXd::Constant(size, 0)); - - c.row(0) << input_field; -} - -void tug::diffusion::BTCSDiffusion::simulate2D( - Eigen::Map &c, Eigen::Map &alpha, - const tug::boundary_condition::BoundaryCondition &bc) { - - int n_rows = this->grid_cells[1]; - int n_cols = this->grid_cells[0]; - double dx = this->deltas[0]; - - double local_dt = this->time_step / BTCS_2D_DT_SIZE; - - DMatrixRowMajor d_ortho = calc_d_ortho(c, alpha, bc, false, local_dt, dx); - -#pragma omp parallel for schedule(dynamic) - for (int i = 0; i < n_rows; i++) { - DVectorRowMajor input_field = c.row(i); - simulate_base(input_field, bc.row_boundary(i), bc.getInnerRow(i), - alpha.row(i), dx, local_dt, n_cols, d_ortho.row(i)); - c.row(i) << input_field; - } - - dx = this->deltas[1]; - - d_ortho = - calc_d_ortho(c.transpose(), alpha.transpose(), bc, true, local_dt, dx); - -#pragma omp parallel for schedule(dynamic) - for (int i = 0; i < n_cols; i++) { - DVectorRowMajor input_field = c.col(i); - simulate_base(input_field, bc.col_boundary(i), bc.getInnerCol(i), - alpha.col(i), dx, local_dt, n_rows, d_ortho.row(i)); - c.col(i) << input_field.transpose(); - } -} - -auto tug::diffusion::BTCSDiffusion::calc_d_ortho( - const DMatrixRowMajor &c, const DMatrixRowMajor &alpha, - const tug::boundary_condition::BoundaryCondition &bc, bool transposed, - double time_step, double dx) -> DMatrixRowMajor { - - uint8_t upper = (transposed ? tug::boundary_condition::BC_SIDE_LEFT - : tug::boundary_condition::BC_SIDE_TOP); - uint8_t lower = (transposed ? tug::boundary_condition::BC_SIDE_RIGHT - : tug::boundary_condition::BC_SIDE_BOTTOM); - - int n_rows = c.rows(); - int n_cols = c.cols(); - - DMatrixRowMajor d_ortho(n_rows, n_cols); - - std::array y_values{}; - - // first, iterate over first row - for (int j = 0; j < n_cols; j++) { - tug::boundary_condition::boundary_condition tmp_bc = bc(upper, j); - double sy = (time_step * alpha(0, j)) / (dx * dx); - - y_values[0] = (tmp_bc.type == tug::boundary_condition::BC_TYPE_CONSTANT - ? tmp_bc.value - : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j))); - y_values[1] = c(0, j); - y_values[2] = c(1, j); - - d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]); - } - -// then iterate over inlet -#pragma omp parallel for private(y_values) schedule(dynamic) - for (int i = 1; i < n_rows - 1; i++) { - for (int j = 0; j < n_cols; j++) { - double sy = (time_step * alpha(i, j)) / (dx * dx); - - y_values[0] = c(i - 1, j); - y_values[1] = c(i, j); - y_values[2] = c(i + 1, j); - - d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]); - } - } - - int end = n_rows - 1; - - // and finally over last row - for (int j = 0; j < n_cols; j++) { - tug::boundary_condition::boundary_condition tmp_bc = bc(lower, j); - double sy = (time_step * alpha(end, j)) / (dx * dx); - - y_values[0] = c(end - 1, j); - y_values[1] = c(end, j); - y_values[2] = (tmp_bc.type == tug::boundary_condition::BC_TYPE_CONSTANT - ? tmp_bc.value - : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j))); - - d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]); - } - - return d_ortho; -} - -void tug::diffusion::BTCSDiffusion::fillMatrixFromRow( - Eigen::SparseMatrix &A_matrix, const DVectorRowMajor &alpha, - const tug::boundary_condition::bc_vec &bc_inner, int size, double dx, - double time_step) { - - double sx = 0; - - int A_size = A_matrix.cols(); - - A_matrix.insert(0, 0) = 1; - - if (bc_inner[0].type != tug::boundary_condition::BC_UNSET) { - if (bc_inner[0].type != tug::boundary_condition::BC_TYPE_CONSTANT) { - throw_invalid_argument("Inner boundary conditions with other type than " - "BC_TYPE_CONSTANT are currently not supported."); - } - A_matrix.insert(1, 1) = 1; - } else { - sx = (alpha[0] * time_step) / (dx * dx); - A_matrix.insert(1, 1) = -1. - 3. * sx; - A_matrix.insert(1, 0) = 2. * sx; - A_matrix.insert(1, 2) = sx; - } - - for (int j = 2, k = j - 1; k < size - 1; j++, k++) { - if (bc_inner[k].type != tug::boundary_condition::BC_UNSET) { - if (bc_inner[k].type != tug::boundary_condition::BC_TYPE_CONSTANT) { - throw_invalid_argument("Inner boundary conditions with other type than " - "BC_TYPE_CONSTANT are currently not supported."); - } - A_matrix.insert(j, j) = 1; - continue; - } - sx = (alpha[k] * time_step) / (dx * dx); - - A_matrix.insert(j, j) = -1. - 2. * sx; - A_matrix.insert(j, (j - 1)) = sx; - A_matrix.insert(j, (j + 1)) = sx; - } - - if (bc_inner[size - 1].type != tug::boundary_condition::BC_UNSET) { - if (bc_inner[size - 1].type != tug::boundary_condition::BC_TYPE_CONSTANT) { - throw_invalid_argument("Inner boundary conditions with other type than " - "BC_TYPE_CONSTANT are currently not supported."); - } - A_matrix.insert(A_size - 2, A_size - 2) = 1; - } else { - sx = (alpha[size - 1] * time_step) / (dx * dx); - A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx; - A_matrix.insert(A_size - 2, A_size - 3) = sx; - A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx; - } - - A_matrix.insert(A_size - 1, A_size - 1) = 1; -} - -void tug::diffusion::BTCSDiffusion::fillVectorFromRow( - Eigen::VectorXd &b_vector, const DVectorRowMajor &c, - const DVectorRowMajor &alpha, const tug::boundary_condition::bc_tuple &bc, - const tug::boundary_condition::bc_vec &bc_inner, - const DVectorRowMajor &d_ortho, int size, double dx, double time_step) { - - tug::boundary_condition::boundary_condition left = bc[0]; - tug::boundary_condition::boundary_condition right = bc[1]; - - bool left_constant = (left.type == tug::boundary_condition::BC_TYPE_CONSTANT); - bool right_constant = - (right.type == tug::boundary_condition::BC_TYPE_CONSTANT); - - int b_size = b_vector.size(); - - for (int j = 0; j < size; j++) { - if (bc_inner[j].type != tug::boundary_condition::BC_UNSET) { - if (bc_inner[j].type != tug::boundary_condition::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; - } - b_vector[j + 1] = -c[j] + d_ortho[j]; - } - - // this is not correct currently.We will fix this when we are able to define - // FLUX boundary conditions - b_vector[0] = - (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); - - b_vector[b_size - 1] = - (right_constant ? right.value - : getBCFromFlux(right, c[size - 1], alpha[size - 1])); -} - -void tug::diffusion::BTCSDiffusion::setTimestep(double time_step) { - this->time_step = time_step; -} - -auto tug::diffusion::BTCSDiffusion::simulate( - double *c, double *alpha, - const tug::boundary_condition::BoundaryCondition &bc) -> double { - - std::chrono::high_resolution_clock::time_point start = - std::chrono::high_resolution_clock::now(); - - if (this->grid_dim == 1) { - Eigen::Map c_in(c, this->grid_cells[0]); - Eigen::Map alpha_in(alpha, this->grid_cells[0]); - - simulate1D(c_in, alpha_in, bc); - } - if (this->grid_dim == 2) { - Eigen::Map c_in(c, this->grid_cells[1], - this->grid_cells[0]); - - Eigen::Map alpha_in(alpha, this->grid_cells[1], - this->grid_cells[0]); - - simulate2D(c_in, alpha_in, bc); - } - - std::chrono::high_resolution_clock::time_point end = - std::chrono::high_resolution_clock::now(); - - std::chrono::duration duration = - std::chrono::duration_cast>(end - start); - - return duration.count(); -} - -inline auto tug::diffusion::BTCSDiffusion::getBCFromFlux( - tug::boundary_condition::boundary_condition bc, double neighbor_c, - double neighbor_alpha) -> double { - - double val = 0; - - if (bc.type == tug::boundary_condition::BC_TYPE_CLOSED) { - val = neighbor_c; - } else if (bc.type == tug::boundary_condition::BC_TYPE_FLUX) { - // TODO - // val = bc[index].value; - } else { - // TODO: implement error handling here. Type was set to wrong value. - } - - return val; -} diff --git a/src/grid/BoundaryCondition.cpp b/src/BoundaryCondition.cpp similarity index 60% rename from src/grid/BoundaryCondition.cpp rename to src/BoundaryCondition.cpp index efb87bf..ba2e403 100644 --- a/src/grid/BoundaryCondition.cpp +++ b/src/BoundaryCondition.cpp @@ -2,13 +2,13 @@ #include #include -#include "../BTCSUtils.hpp" -#include "grid/BoundaryCondition.hpp" +#include "BoundaryCondition.hpp" +#include "TugUtils.hpp" constexpr uint8_t DIM_1D = 2; constexpr uint8_t DIM_2D = 4; -tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x) { +tug::bc::BoundaryCondition::BoundaryCondition(int x) { this->bc_internal.resize(DIM_1D, {0, 0}); this->special_cells.resize(x, {BC_UNSET, 0}); this->dim = 1; @@ -21,7 +21,7 @@ tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x) { this->maxindex = x - 1; } -tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x, int y) { +tug::bc::BoundaryCondition::BoundaryCondition(int x, int y) { this->maxsize = (x >= y ? x : y); this->bc_internal.resize(DIM_2D * maxsize, {0, 0}); this->special_cells.resize(x * y, {BC_UNSET, 0}); @@ -33,8 +33,8 @@ tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x, int y) { this->maxindex = (x * y) - 1; } -void tug::boundary_condition::BoundaryCondition::setSide( - uint8_t side, tug::boundary_condition::boundary_condition &input_bc) { +void tug::bc::BoundaryCondition::setSide( + uint8_t side, tug::bc::boundary_condition &input_bc) { if (this->dim == 1) { throw_invalid_argument("setSide requires at least a 2D grid"); } @@ -42,19 +42,18 @@ void tug::boundary_condition::BoundaryCondition::setSide( throw_out_of_range("Invalid range for 2D grid"); } - uint32_t size = (side == tug::boundary_condition::BC_SIDE_LEFT || - side == tug::boundary_condition::BC_SIDE_RIGHT - ? this->sizes[Y_DIM] - : this->sizes[X_DIM]); + uint32_t size = + (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT + ? this->sizes[Y_DIM] + : this->sizes[X_DIM]); for (uint32_t i = 0; i < size; i++) { this->bc_internal[side * maxsize + i] = input_bc; } } -void tug::boundary_condition::BoundaryCondition::setSide( - uint8_t side, - std::vector &input_bc) { +void tug::bc::BoundaryCondition::setSide( + uint8_t side, std::vector &input_bc) { if (this->dim == 1) { throw_invalid_argument("setSide requires at least a 2D grid"); } @@ -62,10 +61,10 @@ void tug::boundary_condition::BoundaryCondition::setSide( throw_out_of_range("Invalid range for 2D grid"); } - uint32_t size = (side == tug::boundary_condition::BC_SIDE_LEFT || - side == tug::boundary_condition::BC_SIDE_RIGHT - ? this->sizes[Y_DIM] - : this->sizes[X_DIM]); + uint32_t size = + (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT + ? this->sizes[Y_DIM] + : this->sizes[X_DIM]); if (input_bc.size() > size) { throw_out_of_range("Input vector is greater than maximum excpected value"); @@ -76,8 +75,8 @@ void tug::boundary_condition::BoundaryCondition::setSide( } } -auto tug::boundary_condition::BoundaryCondition::getSide(uint8_t side) - -> std::vector { +auto tug::bc::BoundaryCondition::getSide(uint8_t side) + -> std::vector { if (this->dim == 1) { throw_invalid_argument("getSide requires at least a 2D grid"); } @@ -85,12 +84,12 @@ auto tug::boundary_condition::BoundaryCondition::getSide(uint8_t side) throw_out_of_range("Invalid range for 2D grid"); } - uint32_t size = (side == tug::boundary_condition::BC_SIDE_LEFT || - side == tug::boundary_condition::BC_SIDE_RIGHT - ? this->sizes[Y_DIM] - : this->sizes[X_DIM]); + uint32_t size = + (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT + ? this->sizes[Y_DIM] + : this->sizes[X_DIM]); - std::vector out(size); + std::vector out(size); for (int i = 0; i < size; i++) { out[i] = this->bc_internal[this->maxsize * side + i]; @@ -99,8 +98,18 @@ auto tug::boundary_condition::BoundaryCondition::getSide(uint8_t side) return out; } -auto tug::boundary_condition::BoundaryCondition::col_boundary(uint32_t i) const - -> tug::boundary_condition::bc_tuple { +auto tug::bc::BoundaryCondition::row_boundary(uint32_t i) const + -> tug::bc::bc_tuple { + if (i >= this->sizes[Y_DIM]) { + throw_out_of_range("Index out of range"); + } + + return {this->bc_internal[BC_SIDE_LEFT * this->maxsize + i], + this->bc_internal[BC_SIDE_RIGHT * this->maxsize + i]}; +} + +auto tug::bc::BoundaryCondition::col_boundary(uint32_t i) const + -> tug::bc::bc_tuple { if (this->dim == 1) { throw_invalid_argument("Access of column requires at least 2D grid"); } @@ -112,18 +121,7 @@ auto tug::boundary_condition::BoundaryCondition::col_boundary(uint32_t i) const this->bc_internal[BC_SIDE_BOTTOM * this->maxsize + i]}; } -auto tug::boundary_condition::BoundaryCondition::row_boundary(uint32_t i) const - -> tug::boundary_condition::bc_tuple { - if (i >= this->sizes[Y_DIM]) { - throw_out_of_range("Index out of range"); - } - - return {this->bc_internal[BC_SIDE_LEFT * this->maxsize + i], - this->bc_internal[BC_SIDE_RIGHT * this->maxsize + i]}; -} - -auto tug::boundary_condition::BoundaryCondition::getInnerRow(uint32_t i) const - -> bc_vec { +auto tug::bc::BoundaryCondition::getInnerRow(uint32_t i) const -> bc_vec { if (i >= this->sizes[Y_DIM]) { throw_out_of_range("Index is out of range"); } @@ -136,8 +134,7 @@ auto tug::boundary_condition::BoundaryCondition::getInnerRow(uint32_t i) const return row; } -auto tug::boundary_condition::BoundaryCondition::getInnerCol(uint32_t i) const - -> bc_vec { +auto tug::bc::BoundaryCondition::getInnerCol(uint32_t i) const -> bc_vec { if (this->dim != 2) { throw_invalid_argument("getInnerCol is only applicable for 2D grids"); } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 080739a..35932a1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(tug BTCSDiffusion.cpp grid/BoundaryCondition.cpp) +add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp) target_link_libraries(tug Eigen3::Eigen) diff --git a/src/Solver.cpp b/src/Solver.cpp new file mode 100644 index 0000000..4cf6f1c --- /dev/null +++ b/src/Solver.cpp @@ -0,0 +1,62 @@ +#include +#include + +#include +#include + +auto tug::solver::EigenLU(const Eigen::SparseMatrix A_matrix, + const Eigen::VectorXd b_vector) -> Eigen::VectorXd { + + Eigen::SparseLU, Eigen::COLAMDOrdering> + solver; + solver.analyzePattern(A_matrix); + + solver.factorize(A_matrix); + + Eigen::VectorXd x_vec = solver.solve(b_vector); + + return solver.solve(b_vector); +} + +auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix A_matrix, + const Eigen::VectorXd b_vector) + -> Eigen::VectorXd { + uint32_t n = b_vector.size(); + + Eigen::VectorXd a_diag(n); + Eigen::VectorXd b_diag(n); + Eigen::VectorXd c_diag(n); + Eigen::VectorXd x_vec = b_vector; + + // Fill diagonals vectors + b_diag[0] = A_matrix.coeff(0, 0); + c_diag[0] = A_matrix.coeff(0, 1); + + for (int i = 1; i < n - 1; i++) { + a_diag[i] = A_matrix.coeff(i, i - 1); + b_diag[i] = A_matrix.coeff(i, i); + c_diag[i] = A_matrix.coeff(i, i + 1); + } + a_diag[n - 1] = A_matrix.coeff(n - 1, n - 2); + b_diag[n - 1] = A_matrix.coeff(n - 1, n - 1); + + // 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]); + } + + 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]; + } + + return x_vec; +} diff --git a/src/BTCSUtils.hpp b/src/TugUtils.hpp similarity index 52% rename from src/BTCSUtils.hpp rename to src/TugUtils.hpp index f574e7c..502a061 100644 --- a/src/BTCSUtils.hpp +++ b/src/TugUtils.hpp @@ -1,6 +1,7 @@ #ifndef BTCSUTILS_H_ #define BTCSUTILS_H_ +#include #include #include @@ -12,4 +13,14 @@ #define throw_out_of_range(msg) \ throw std::out_of_range(std::string(__FILE__) + ":" + \ std::to_string(__LINE__) + ":" + std::string(msg)) + +#define time_marker() std::chrono::high_resolution_clock::now() + +#define diff_time(start, end) \ + ({ \ + std::chrono::duration duration = \ + std::chrono::duration_cast>(end - \ + start); \ + duration.count(); \ + }) #endif // BTCSUTILS_H_ diff --git a/test/simple.cpp b/test/simple.cpp new file mode 100644 index 0000000..b3b6622 --- /dev/null +++ b/test/simple.cpp @@ -0,0 +1,36 @@ +#include "grid/BoundaryCondition.hpp" +#include +#include +#include + +#define SIZE 5 + +using namespace tug; + +int main() { + boundary_condition::BoundaryCondition bc(SIZE); + + diffusion::TugDiffuInput in = {1., {SIZE, 0, 0}, {1, 0, 0}, 0}; + + std::vector field1(SIZE, 0); + std::vector field2(SIZE, 0); + std::vector alpha(SIZE, 1e-5); + + bc(boundary_condition::BC_SIDE_LEFT) = {boundary_condition::BC_TYPE_CONSTANT, + 1}; + + for (int j = 0; j < 10; j++) { + double time1 = + diffusion::BTCS_Thomas_1D(field1.data(), alpha.data(), bc, in); + double time2 = + diffusion::BTCS_EigenLU_1D(field2.data(), alpha.data(), bc, in); + + for (int i = 0; i < SIZE; i++) { + std::cout << field1[i] << "\t"; + std::cout << field2[i] << "\n"; + } + + std::cout << "Time Thomas = " << time1 << "s \tTime EigenLU = " << time2 + << " s" << std::endl; + } +} diff --git a/test/testBoundaryCondition.cpp b/test/testBoundaryCondition.cpp index 9b6c01e..3eda557 100644 --- a/test/testBoundaryCondition.cpp +++ b/test/testBoundaryCondition.cpp @@ -1,7 +1,7 @@ -#include +#include #include -using namespace tug::boundary_condition; +using namespace tug::bc; #define BC_CONST_VALUE 1e-5 @@ -87,8 +87,8 @@ TEST_CASE("Boundary Condition helpers") { boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; SUBCASE("return boundary condition skeleton") { - boundary_condition bc_test = BoundaryCondition::returnBoundaryCondition( - bc_set.type, bc_set.value); + boundary_condition bc_test = + BoundaryCondition::returnBoundaryCondition(bc_set.type, bc_set.value); CHECK_EQ(bc_test.value, bc_set.value); CHECK_EQ(bc_test.type, bc_set.type); } @@ -98,9 +98,7 @@ TEST_CASE("1D special inner grid cells") { BoundaryCondition bc(5); boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; - SUBCASE("valid set") { - CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set); - } + SUBCASE("valid set") { CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set); } SUBCASE("valid get") { bc(BC_INNER, 0) = bc_set; @@ -112,9 +110,7 @@ TEST_CASE("1D special inner grid cells") { CHECK_THROWS(bc(BC_INNER, 5)); } - SUBCASE("invalid set") { - CHECK_THROWS(bc(BC_INNER, 5) = bc_set); - } + SUBCASE("invalid set") { CHECK_THROWS(bc(BC_INNER, 5) = bc_set); } SUBCASE("valid row getter") { bc(BC_INNER, 1) = bc_set; @@ -125,23 +121,16 @@ TEST_CASE("1D special inner grid cells") { CHECK_EQ(ret[1].type, bc_set.type); } - SUBCASE("invalid row getter") { - CHECK_THROWS(bc.getInnerRow(1)); - } - - SUBCASE("invalid col getter") { - CHECK_THROWS(bc.getInnerCol(0)); - } + SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(1)); } + SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(0)); } } TEST_CASE("2D special inner grid cells") { - BoundaryCondition bc(5,5); + BoundaryCondition bc(5, 5); boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE}; - SUBCASE("valid set") { - CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set); - } + SUBCASE("valid set") { CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set); } SUBCASE("valid get") { bc(BC_INNER, 0) = bc_set; @@ -153,9 +142,7 @@ TEST_CASE("2D special inner grid cells") { CHECK_THROWS(bc(BC_INNER, 25)); } - SUBCASE("invalid set") { - CHECK_THROWS(bc(BC_INNER, 25) = bc_set); - } + SUBCASE("invalid set") { CHECK_THROWS(bc(BC_INNER, 25) = bc_set); } SUBCASE("valid row getter") { bc(BC_INNER, 0) = bc_set; @@ -187,11 +174,7 @@ TEST_CASE("2D special inner grid cells") { CHECK_EQ(ret[1].type, BC_UNSET); } - SUBCASE("invalid row getter") { - CHECK_THROWS(bc.getInnerRow(5)); - } + SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(5)); } - SUBCASE("invalid col getter") { - CHECK_THROWS(bc.getInnerCol(5)); - } + SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(5)); } } diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 9e0796a..2def0c3 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -1,10 +1,11 @@ +#include "Diffusion.hpp" +#include "Solver.hpp" +#include #include -#include #include -#include #include -using namespace tug::boundary_condition; +using namespace tug::bc; using namespace tug::diffusion; #define DIMENSION 2 @@ -14,13 +15,19 @@ using namespace tug::diffusion; static std::vector alpha(N *M, 1e-3); -static BTCSDiffusion setupDiffu(uint32_t n, uint32_t m) { - BTCSDiffusion diffu(DIMENSION); +static TugInput setupDiffu(BoundaryCondition &bc) { + TugInput diffu; - diffu.setXDimensions(n, n); - diffu.setYDimensions(m, m); + diffu.time_step = 1.; + diffu.solver = tug::solver::ThomasAlgorithm; - diffu.setTimestep(1.); + diffu.grid.grid_cells[0] = N; + diffu.grid.grid_cells[1] = M; + + diffu.grid.domain_size[0] = N; + diffu.grid.domain_size[1] = M; + + diffu.grid.bc = &bc; return diffu; } @@ -30,23 +37,22 @@ TEST_CASE("closed boundaries - 1 concentration to 1 - rest 0") { field[MID] = 1; - BTCSDiffusion diffu = setupDiffu(N, M); BoundaryCondition bc(N, M); + TugInput diffu = setupDiffu(bc); + uint32_t iterations = 1000; double sum = 0; for (int t = 0; t < iterations; t++) { - diffu.simulate(field.data(), alpha.data(), bc); + ADI_2D(diffu, field.data(), alpha.data()); + } - if (t == iterations - 1) { - // iterate through rows - for (int i = 0; i < M; i++) { - // iterate through columns - for (int j = 0; j < N; j++) { - sum += field[i * N + j]; - } - } + // iterate through rows + for (int i = 0; i < M; i++) { + // iterate through columns + for (int j = 0; j < N; j++) { + sum += field[i * N + j]; } } CAPTURE(sum); @@ -59,7 +65,6 @@ TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") { field[MID] = 1; - BTCSDiffusion diffu = setupDiffu(N, M); BoundaryCondition bc(N, M); boundary_condition input = {BC_TYPE_CONSTANT, 0}; @@ -69,13 +74,15 @@ 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); + uint32_t max_iterations = 20000; bool reached = false; int t = 0; for (t = 0; t < max_iterations; t++) { - diffu.simulate(field.data(), alpha.data(), bc); + ADI_2D(diffu, field.data(), alpha.data()); if (field[N * M - 1] > 1e-15) { reached = true; @@ -95,7 +102,6 @@ TEST_CASE( "constant top and bottom (1 and 0) - left and right closed - 0 inlet") { std::vector field(N * M, 0); - BTCSDiffusion diffu = setupDiffu(N, M); BoundaryCondition bc(N, M); boundary_condition top = @@ -106,10 +112,12 @@ TEST_CASE( bc.setSide(BC_SIDE_TOP, top); bc.setSide(BC_SIDE_BOTTOM, bottom); + TugInput diffu = setupDiffu(bc); + uint32_t max_iterations = 100; for (int t = 0; t < max_iterations; t++) { - diffu.simulate(field.data(), alpha.data(), bc); + ADI_2D(diffu, field.data(), alpha.data()); } for (int i = 0; i < N; i++) { @@ -129,18 +137,19 @@ TEST_CASE("2D closed boundaries, 1 constant cell in the middle") { std::vector field(N * M, 0); double val = 1e-2; - BTCSDiffusion diffu = setupDiffu(N, M); BoundaryCondition bc(N, M); field[MID] = val; bc(BC_INNER, MID) = {BC_TYPE_CONSTANT, val}; + TugInput diffu = setupDiffu(bc); + uint32_t max_iterations = 100; double sum = val; for (int t = 0; t < max_iterations; t++) { - diffu.simulate(field.data(), alpha.data(), bc); + ADI_2D(diffu, field.data(), alpha.data()); CHECK_EQ(field[MID], val);