diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index d5d7ac1..88c816f 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -3,46 +3,101 @@ #include #include +#include #include #include +#include #include #include -const int BTCSDiffusion::BC_NEUMANN = 0; -const int BTCSDiffusion::BC_DIRICHLET = 1; +const int BTCSDiffusion::BC_CONSTANT = 0; +const int BTCSDiffusion::BC_CLOSED = 1; +const int BTCSDiffusion::BC_FLUX = 2; -BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) { - this->grid_dim = 1; - this->dx = 1. / (x - 1); +BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { + assert(dim <= 3); - // per default use Neumann condition with gradient of 0 at the end of the grid - this->bc.resize(2, std::tuple(BTCSDiffusion::BC_NEUMANN, 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 + grid_cells.resize(dim, 1); + domain_size.resize(dim, 1); + deltas.resize(dim, 1); } -void BTCSDiffusion::simulate1D(std::vector &c, double bc_left, - double bc_right, +std::vector BTCSDiffusion::getNumberOfGridCells() { + return this->grid_cells; +} +std::vector BTCSDiffusion::getSpatialDiscretization() { + return this->domain_size; +} +void BTCSDiffusion::setNumberOfGridCells(std::vector &n_grid) { + grid_cells = n_grid; + assert(grid_cells.size() == grid_dim); + updateInternals(); +} +void BTCSDiffusion::setSpatialDiscretization(std::vector &s_grid) { + domain_size = s_grid; + assert(domain_size.size() == grid_dim); + updateInternals(); +} + +void BTCSDiffusion::updateInternals() { + for (int i = 0; i < grid_dim; i++) { + deltas[i] = (double)domain_size[i] / grid_cells[i]; + } + + switch (grid_dim) { + case 1: + bc.resize(2, {BTCSDiffusion::BC_CLOSED, 0}); + break; + case 2: + bc.resize(2 * grid_cells[0] + 2 * grid_cells[1], + {BTCSDiffusion::BC_CLOSED, 0}); + break; + case 3: + // TODO + break; + } +} +// 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, const std::vector &alpha, double dx, int size) { + 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; + // size = size + 2; + + int bc_offset = !left_is_constant + !right_is_constant; + ; // set sizes of private and yet allocated vectors - b_vector.resize(size); - x_vector.resize(size); + b_vector.resize(size + bc_offset); + x_vector.resize(size + bc_offset); /* * Begin to solve the equation system using LU solver of Eigen. @@ -53,25 +108,33 @@ void BTCSDiffusion::simulate1D(std::vector &c, double bc_left, * TODO: remove output */ - A_matrix.resize(size, size); - A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); + A_matrix.resize(size + bc_offset, size + bc_offset); + A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3)); A_matrix.insert(0, 0) = 1; - A_matrix.insert(size - 1, size - 1) = 1; + b_vector[0] = + (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); - b_vector[0] = bc_left; - b_vector[size - 1] = bc_right; + A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1; + b_vector[size + bc_offset - 1] = + (right_is_constant ? right.value + : getBCFromFlux(right, c[size - 1], alpha[size - 1])); - for (int i = 1; i < this->n_x + 1; i++) { - double sx = (alpha[i - 1] * time_step) / (dx * dx); + // A_matrix.insert(0, 0) = 1; + // A_matrix.insert(size + 1, size + 1) = 1; + + for (int i = 1; i < size - right_is_constant; i++) { + double sx = (alpha[i + !(left_is_constant)] * time_step) / (dx * dx); A_matrix.insert(i, i) = -1. - 2. * sx; A_matrix.insert(i, i - 1) = sx; A_matrix.insert(i, i + 1) = sx; - b_vector[i] = -c[i - 1]; + b_vector[i] = -c[i + !(left_is_constant)]; } + std::cout << b_vector << "\n" << A_matrix << std::endl; + Eigen::SparseLU, Eigen::COLAMDOrdering> solver; solver.analyzePattern(A_matrix); @@ -85,7 +148,7 @@ void BTCSDiffusion::simulate1D(std::vector &c, double bc_left, std::cout << std::setprecision(10) << x_vector << std::endl << std::endl; for (int i = 0; i < c.size(); i++) { - c[i] = x_vector[i + 1]; + c[i] = x_vector[i + !left_is_constant]; } } @@ -96,24 +159,25 @@ 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]); + // 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_left, bc_right, alpha, this->dx, this->n_x); + simulate1D(c, bc[0], bc[1], alpha, this->deltas[0], this->grid_cells[0]); } } -double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c, - double neighbor_alpha) { - double val = -1; - int type = std::get<0>(bc[index]); +inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, + double neighbor_c, + double neighbor_alpha) { - if (type == BTCSDiffusion::BC_NEUMANN) { - val = neighbor_c + (this->time_step / (dx * dx)) * neighbor_alpha * - std::get<1>(bc[index]); - } else if (type == BTCSDiffusion::BC_DIRICHLET) { - val = std::get<1>(bc[index]); + double val; + + if (bc.type == BTCSDiffusion::BC_CLOSED) { + val = neighbor_c; + } else if (bc.type == BTCSDiffusion::BC_FLUX) { + // TODO + // val = bc[index].value; } else { // TODO: implement error handling here. Type was set to wrong value. } @@ -122,6 +186,10 @@ double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c, } void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) { - std::get<0>(bc[index]) = type; - std::get<1>(bc[index]) = val; + + bc[index].type = type; + bc[index].value = val; + + // std::get<0>(bc[index]) = type; + // std::get<1>(bc[index]) = val; } diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 355c646..af08e0a 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -5,28 +5,11 @@ #include #include -/*! - * Datatype to fill the sparse matrix which is used to solve the equation - * system. - */ -typedef Eigen::Triplet T; - /*! * Defines both types of boundary condition as a datatype. */ typedef int bctype; -/*! - * 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; - /*! * Class implementing a solution for a 1/2/3D diffusion equation using backward * euler. @@ -35,37 +18,79 @@ class BTCSDiffusion { public: /*! - * Defines a Neumann boundary condition. + * Defines a constant/Dirichlet boundary condition. */ - static const int BC_NEUMANN; + static const int BC_CONSTANT; + /*! - * Defines a Dirichlet boundary condition. + * Defines a closed/Neumann boundary condition. */ - static const int BC_DIRICHLET; + static const int BC_CLOSED; + + /*! + * Defines a flux/Cauchy boundary condition. + */ + static const int BC_FLUX; + + /*! + * 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 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. */ - explicit BTCSDiffusion(int x); + BTCSDiffusion(unsigned int dim); - /*! - * 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); + std::vector getNumberOfGridCells(); + std::vector getSpatialDiscretization(); + void setNumberOfGridCells(std::vector &n_grid); + void setSpatialDiscretization(std::vector &s_grid); - /*! - * 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); + // /*! + // * 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); /*! * With given ghost zones simulate diffusion. Only 1D allowed at this moment. @@ -96,14 +121,17 @@ public: void setBoundaryCondition(int index, double val, bctype type); private: - void simulate1D(std::vector &c, double bc_left, double bc_right, - const std::vector &alpha, double dx, int size); + 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); - double getBCFromTuple(int index, double nearest_value, double neighbor_alpha); + inline double getBCFromFlux(boundary_condition bc, double nearest_value, double neighbor_alpha); - boundary_condition bc; + void updateInternals(); + + std::vector bc; Eigen::SparseMatrix A_matrix; Eigen::VectorXd b_vector; @@ -112,12 +140,9 @@ private: double time_step; int grid_dim; - int n_x; - double dx; - int n_y; - double dy; - int n_z; - double dz; + std::vector grid_cells; + std::vector domain_size; + std::vector deltas; }; #endif // BTCSDIFFUSION_H_ diff --git a/src/main.cpp b/src/main.cpp index 9735f11..c4549ee 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -6,19 +6,27 @@ using namespace std; int main(int argc, char *argv[]) { - // count of grid cells - int x = 20; + // dimension of grid + int dim = 1; + + int n = 20; // create input + diffusion coefficients for each grid cell - std::vector alpha(x, 1 * pow(10, -1)); - std::vector input(x, 1 * std::pow(10, -6)); + std::vector alpha(n, 1 * pow(10, -1)); + std::vector field(n, 1 * std::pow(10, -6)); // create instance of diffusion module - BTCSDiffusion diffu(x); + BTCSDiffusion diffu(dim); + + std::vector vec_n = diffu.getNumberOfGridCells(); + + vec_n[0] = n; + + diffu.setNumberOfGridCells(vec_n); // set the boundary condition for the left ghost cell to dirichlet diffu.setBoundaryCondition(0, 5. * std::pow(10, -6), - BTCSDiffusion::BC_DIRICHLET); + BTCSDiffusion::BC_CONSTANT); // set timestep for simulation to 1 second diffu.setTimestep(1.); @@ -26,7 +34,7 @@ int main(int argc, char *argv[]) { // loop 100 times // output is currently generated by the method itself for (int i = 0; i < 100; i++) { - diffu.simulate(input, alpha); + diffu.simulate(field, alpha); } return 0;