diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index f4e8b4b..92855dc 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -67,28 +67,6 @@ void BTCSDiffusion::updateInternals() { bc.resize(cells, {BTCSDiffusion::BC_CLOSED, 0}); } -// BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) { -// this->grid_dim = 1; -// this->dx = 1. / (x - 1); - -// // per default use Neumann condition with gradient of 0 at the end of the -// grid this->bc.resize(2, std::tuple(BTCSDiffusion::BC_CONSTANT, 0.)); -// } -// BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) { - -// // this->grid_dim = 2; - -// // this->bc.reserve(x * 2 + y * 2); -// // // per default use Neumann condition with gradient of 0 at the end of -// the -// // grid std::fill(this->bc.begin(), this->bc.end(), -1); -// } -// BTCSDiffusion::BTCSDiffusion(int x, int y, int z) : n_x(x), n_y(y), n_z(z) { - -// // this->grid_dim = 3; -// // TODO: reserve memory for boundary conditions -// } void BTCSDiffusion::simulate1D(std::vector &c, boundary_condition left, boundary_condition right, @@ -97,11 +75,11 @@ void BTCSDiffusion::simulate1D(std::vector &c, boundary_condition left, bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT); bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT); - int loop_end = size + !right_is_constant; - - // we need 2 more grid cells for ghost cells - // size = size + 2; + //The sizes for matrix and vectors of the equation system is defined by the + //actual size of the input vector and if the system is (partially) closed. + //Then we will need ghost nodes. So this variable will give the count of ghost + //nodes. int bc_offset = !left_is_constant + !right_is_constant; ; @@ -113,11 +91,10 @@ void BTCSDiffusion::simulate1D(std::vector &c, boundary_condition left, * Begin to solve the equation system using LU solver of Eigen. * * But first fill the A matrix and b vector. - * - * At this point there is some debugging output in the code. - * TODO: remove output */ + // Set boundary condition for ghost nodes (for closed or flux system) or outer + // inlet nodes (constant boundary condition) A_matrix.resize(size + bc_offset, size + bc_offset); A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3)); @@ -130,11 +107,13 @@ void BTCSDiffusion::simulate1D(std::vector &c, boundary_condition left, (right_is_constant ? right.value : getBCFromFlux(right, c[size - 1], alpha[size - 1])); - // A_matrix.insert(0, 0) = 1; - // A_matrix.insert(size + 1, size + 1) = 1; - + // Start filling the A matrix + // =i= is used for equation system matrix and vector indexing + // and =j= for indexing of c,alpha and bc for (int i = 1, j = i + !(left_is_constant); i < size - right_is_constant; i++, j++) { + + // if current grid cell is considered as constant boundary conditon if (bc[j].type == BTCSDiffusion::BC_CONSTANT) { A_matrix.insert(i, i) = 1; b_vector[i] = bc[j].value; @@ -150,6 +129,7 @@ void BTCSDiffusion::simulate1D(std::vector &c, boundary_condition left, b_vector[i] = -c[j]; } + // start to solve Eigen::SparseLU, Eigen::COLAMDOrdering> solver; solver.analyzePattern(A_matrix); @@ -158,7 +138,8 @@ void BTCSDiffusion::simulate1D(std::vector &c, boundary_condition left, x_vector = solver.solve(b_vector); - for (int i = 0; i < c.size(); i++) { + //fill solution back in place into =c= vector + for (int i = 0, j = i + !left_is_constant; i < c.size(); i++, j++) { c[i] = x_vector[i + !left_is_constant]; } } @@ -170,10 +151,6 @@ void BTCSDiffusion::setTimestep(double time_step) { void BTCSDiffusion::simulate(std::vector &c, const std::vector &alpha) { if (this->grid_dim == 1) { - // double bc_left = getBCFromTuple(0, c[0], alpha[0]); - // double bc_right = - // getBCFromTuple(1, c[c.size() - 1], alpha[alpha.size() - 1]); - simulate1D(c, bc[0], bc[grid_cells[0] + 1], alpha, this->deltas[0], this->grid_cells[0]); } @@ -197,11 +174,8 @@ inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, return val; } -void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) { +void BTCSDiffusion::setBoundaryCondition(int index, bctype type, double value) { bc[index].type = type; - bc[index].value = val; - - // std::get<0>(bc[index]) = type; - // std::get<1>(bc[index]) = val; + bc[index].value = value; } diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 7954022..ea3473d 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -34,76 +34,76 @@ public: static const int BC_FLUX; /*! - * A boundary condition consists of two features. A type and the according - * value. Here we can differentiate between: + * Creates a diffusion module. * - * - Neumann boundary conditon: type BC_NEUMANN with the value defining the - * gradient - * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of - * the boundary condition - */ - typedef struct boundary_condition { - bctype type; - double value; - } boundary_condition; - - /*! - * A boundary condition consists of two features. A type and the according - * value. Here we can differentiate between: - * - * - Neumann boundary conditon: type BC_NEUMANN with the value defining the - * gradient - * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of - * the boundary condition - */ - // typedef std::vector> boundary_condition; - - /*! - * Datatype to fill the sparse matrix which is used to solve the equation - * system. - */ - typedef Eigen::Triplet T; - - /*! - * Create 1D-diffusion module. - * - * @param x Count of cells in x direction. + * @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(unsigned int 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(unsigned int 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(unsigned int domain_size, unsigned int n_grid_cells); + /*! + * Returns the number of grid cells in x direction. + */ unsigned int getXGridCellsN(); + /*! + * Returns the number of grid cells in y direction. + */ unsigned int getYGridCellsN(); + /*! + * Returns the number of grid cells in z direction. + */ unsigned int getZGridCellsN(); - unsigned int getXDomainSize(); - unsigned int getYDomainSize(); - unsigned int getZDomainSize(); - // /*! - // * Currently not implemented: Create 2D-diffusion module. - // * - // * @param x Count of cells in x direction. - // * @param y Count of cells in y direction. - // */ - // explicit BTCSDiffusion(int x, int y); - // /*! - // * Currently not implemented: Create 3D-diffusion module. - // * - // * @param x Count of cells in x direction. - // * @param y Count of cells in y direction. - // * @param z Count of cells in z direction. - // */ - // explicit BTCSDiffusion(int x, int y, int z); + /*! + * Returns the domain size in x direction. + */ + unsigned int getXDomainSize(); + /*! + * Returns the domain size in y direction. + */ + unsigned int getYDomainSize(); + /*! + * Returns the domain size in z direction. + */ + unsigned int getZDomainSize(); /*! * With given ghost zones simulate diffusion. Only 1D allowed at this moment. * * @param c Vector describing the concentration of one solution of the grid as - * continious memory (Row-wise). - * @param alpha Vector of diffusioncoefficients for each grid element. + * continious memory (row major). + * @param alpha Vector of diffusion coefficients for each grid element. */ void simulate(std::vector &c, const std::vector &alpha); @@ -119,23 +119,29 @@ public: * index (exact order still to be determined), the type of the boundary * condition and the according value. * - * @param index Index of the boundary condition vector. - * @param val Value of the boundary condition (gradient for Neumann, exact - * value for Dirichlet). - * @param Type of the grid cell. + * @param index Index of the grid cell the boundary condition is applied to. + * @param type Type of the boundary condition. Must be constant, closed or + * flux. + * @param value For constant boundary conditions this value is set + * during solving. For flux value refers to a gradient of change for this grid + * cell. For closed this value has no effect since a gradient of 0 is used. */ - void setBoundaryCondition(int index, double val, bctype type); + void setBoundaryCondition(int index, bctype type, double value); private: + typedef struct boundary_condition { + bctype type; + double value; + } boundary_condition; + typedef Eigen::Triplet T; + void simulate1D(std::vector &c, boundary_condition left, boundary_condition right, const std::vector &alpha, double dx, int size); void simulate2D(std::vector &c); void simulate3D(std::vector &c); - inline double getBCFromFlux(boundary_condition bc, double nearest_value, double neighbor_alpha); - void updateInternals(); std::vector bc; @@ -145,8 +151,8 @@ private: Eigen::VectorXd x_vector; double time_step; - int grid_dim; + std::vector grid_cells; std::vector domain_size; std::vector deltas; diff --git a/src/main.cpp b/src/main.cpp index 9802e1f..522c5c6 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,9 +1,8 @@ #include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET -#include -#include // for copy, max +#include // for copy, max #include -#include // for std -#include // for vector +#include // for std +#include // for vector using namespace std; int main(int argc, char *argv[]) { @@ -23,8 +22,8 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); // set the boundary condition for the left ghost cell to dirichlet - diffu.setBoundaryCondition(0, 5. * std::pow(10, -6), - BTCSDiffusion::BC_CONSTANT); + diffu.setBoundaryCondition(0, BTCSDiffusion::BC_CONSTANT, + 5. * std::pow(10, -6)); // set timestep for simulation to 1 second diffu.setTimestep(1.);