diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp new file mode 100644 index 0000000..fcebc92 --- /dev/null +++ b/src/BTCSDiffusion.cpp @@ -0,0 +1,129 @@ +#include "BTCSDiffusion.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +const BCSide BTCSDiffusion::LEFT = 0; +const BCSide BTCSDiffusion::RIGHT = 1; + +BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { + this->grid_dim = 1; + + // per default use Neumann condition with gradient of 0 at the end of the grid + this->bc.resize(2, -1); +} +BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_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) + : dim_x(x), dim_y(y), dim_z(z) { + + this->grid_dim = 3; + // TODO: reserve memory for boundary conditions +} + +void BTCSDiffusion::setBoundaryCondition(std::vector input, + BCSide side) { + if (this->grid_dim == 1) { + bc[side] = input[0]; + } +} +void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, + double timestep) { + // calculate dx + double dx = 1. / (this->dim_x - 1); + + // calculate size needed for A matrix and b,x vectors + int size = this->dim_x + 2; + + Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); + Eigen::VectorXd x_out(size); + + /* + * Initalization of matrix A + * This is done by triplets. See: + * https://eigen.tuxfamily.org/dox/group__TutorialSparse.html + */ + + std::vector tripletList; + tripletList.reserve(c.size() * 3 + bc.size()); + + int A_line = 0; + + // For all concentrations create one row in matrix A + for (int i = 1; i < this->dim_x + 1; i++) { + double sx = (alpha[i - 1] * timestep) / (dx * dx); + + tripletList.push_back(T(A_line, i, (-1. - 2. * sx))); + + tripletList.push_back(T(A_line, i - 1, sx)); + tripletList.push_back(T(A_line, i + 1, sx)); + + b[A_line] = -c[i - 1]; + A_line++; + } + + // append left and right boundary conditions/ghost zones + tripletList.push_back(T(A_line, 0, 1)); + + // if value is -1 apply Neumann condition with given gradient + // TODO: set specific gradient + if (bc[0] == -1) + b[A_line] = c[0]; + // else apply given Dirichlet condition + else + b[A_line] = this->bc[0]; + + A_line++; + tripletList.push_back(T(A_line, size - 1, 1)); + // b[A_line] = bc[1]; + if (bc[1] == -1) + b[A_line] = c[c.size() - 1]; + else + b[A_line] = this->bc[1]; + + /* + * Begin to solve the equation system + * + * At this point there is some debugging output in the code. + * TODO: remove output + */ + + Eigen::SparseMatrix A(size, size); + A.setFromTriplets(tripletList.begin(), tripletList.end()); + + Eigen::SparseLU, Eigen::COLAMDOrdering> + solver; + solver.analyzePattern(A); + + solver.factorize(A); + + std::cout << solver.lastErrorMessage() << std::endl; + + x_out = solver.solve(b); + + std::cout << std::setprecision(10) << x_out << std::endl << std::endl; + + for (int i = 0; i < c.size(); i++) { + c[i] = x_out[i + 1]; + } +} diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp new file mode 100644 index 0000000..cadc1e7 --- /dev/null +++ b/src/BTCSDiffusion.hpp @@ -0,0 +1,89 @@ +#ifndef BTCSDIFFUSION_H_ +#define BTCSDIFFUSION_H_ + +#include +#include + +/*! + * Type defining the side of given boundary condition. + */ +typedef int BCSide; + +/*! + * Datatype to fill the sparse matrix which is used to solve the equation + * system. + */ +typedef Eigen::Triplet T; + +/*! + * Class implementing a solution for a 1/2/3D diffusion equation using backward + * euler. + */ +class BTCSDiffusion { + +public: + /*! + * Set left boundary condition. + */ + static const BCSide LEFT; + + /*! + * Set right boundary condition. + */ + static const BCSide RIGHT; + + /*! + * Create 1D-diffusion module. + * + * @param x Count of cells in x direction. + */ + BTCSDiffusion(int x); + + /*! + * Currently not implemented: Create 2D-diffusion module. + * + * @param x Count of cells in x direction. + * @param y Count of cells in y direction. + */ + 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. + */ + BTCSDiffusion(int x, int y, int z); + + /*! + * Sets internal boundary condition at the end of the grid/ghost zones. + * Currently only implemented for 1D diffusion. + * + * @param input Vector containing all the values to initialize the ghost + * zones. + * @param side Sets the side of the boundary condition. See BCSide for more + * information. + */ + void setBoundaryCondition(std::vector input, BCSide side); + + /*! + * 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. + * @param timestep Time (in seconds ?) to simulate. + */ + void simulate(std::vector &c, std::vector &alpha, + double timestep); + +private: + std::vector bc; + int grid_dim; + int dim_x; + int dim_y; + int dim_z; +}; + +#endif // BTCSDIFFUSION_H_ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 33e2d28..fab0508 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,4 @@ -add_library(diffusion OBJECT diffusion.cpp diffusion.hpp) +add_library(diffusion OBJECT BTCSDiffusion.cpp BTCSDiffusion.hpp) target_link_libraries(diffusion Eigen3::Eigen) add_executable(test main.cpp) diff --git a/src/diffusion.cpp b/src/diffusion.cpp deleted file mode 100644 index 7ecad7f..0000000 --- a/src/diffusion.cpp +++ /dev/null @@ -1,122 +0,0 @@ -#include "diffusion.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -void BTCS2D(int x, int y, std::vector &c, std::vector &alpha, - double timestep) { - - double dx = 1. / x; - double dy = 1. / y; - - int size = (x * y) - 4; - - int local_x = x - 2; - - Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); - Eigen::VectorXd x_out(size); - std::vector tripletList; - tripletList.reserve(size * 5); - - for (int i = x - 1 ; i < 2*x - 3 ; i++) { - double sx = (alpha[i + 2] * timestep) / (dx * dx); - double sy = (alpha[i + 2] * timestep) / (dy * dy); - - - tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy))); - tripletList.push_back(T(i, i - (x - 1), sy)); - tripletList.push_back(T(i, i + x, sy)); - tripletList.push_back(T(i, i + 1, sx)); - tripletList.push_back(T(i, i - 1, sx)); - - b[i] = -c[i+2]; - } - - for (int i = 2*x - 1; i < (y-2)*x-3; i++) { - double sx = (alpha[i + 2] * timestep) / (dx * dx); - double sy = (alpha[i + 2] * timestep) / (dy * dy); - - tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy))); - tripletList.push_back(T(i, i - x, sy)); - tripletList.push_back(T(i, i + x, sy)); - tripletList.push_back(T(i, i + 1, sx)); - tripletList.push_back(T(i, i - 1, sx)); - - b[i] = -c[i+2]; - } - - for (int i = (y-2)*x-1; i < (y-1)*x-3; i++) { - double sx = (alpha[i + 2] * timestep) / (dx * dx); - double sy = (alpha[i + 2] * timestep) / (dy * dy); - - tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy))); - tripletList.push_back(T(i, i - x, sy)); - tripletList.push_back(T(i, i + (x-1), sy)); - tripletList.push_back(T(i, i + 1, sx)); - tripletList.push_back(T(i, i - 1, sx)); - - b[i] = -c[i+2]; - } - - // for (int i = 0; i < (size-local_x); i++) { - // int current = local_x + i; - // double sx = (alpha[current] * timestep) / (dx * dx); - // double sy = (alpha[current] * timestep) / (dy * dy); - - // tripletList.push_back(T(i, current, (1 + 2 * sx + 2 * sy))); - // tripletList.push_back(T(i, current + local_x, sy)); - // tripletList.push_back(T(i, current - local_x, sy)); - // tripletList.push_back(T(i, current + 1, sx)); - // tripletList.push_back(T(i, current - 1, sx)); - - // std::cout << current << std::endl; - - // b[i] = -c[x+i]; - // } - - // for (int i = 0; i < y; i++) { - // for (int j = 0; j < x; j++) { - // double sx = (alpha[i * y + j] * timestep) / (dx * dx); - // double sy = (alpha[i * y + j] * timestep) / (dy * dy); - - // tripletList.push_back(T((i * x) + j, (i * x) + j, (1 + 2 * sx + 2 * - // sy))); tripletList.push_back(T((i * x) + j, ((i * x) + j) + x, sy)); - // tripletList.push_back(T((i * x) + j, ((i * x) + j) - x, sy)); - // tripletList.push_back(T((i * x) + j, ((i * x) + j) + 1, sx)); - // tripletList.push_back(T((i * x) + j, ((i * x) + j) - 1, sx)); - - // b[(i * x) + j] = -c[(i * x) + j]; - // } - // } - - std::cout << b << std::endl; - - Eigen::SparseMatrix A(size, (x * y) - 4); - A.setFromTriplets(tripletList.begin(), tripletList.end()); - -Eigen::SparseQR, Eigen::COLAMDOrdering> solver; - - // Eigen::SparseLU, Eigen::COLAMDOrdering> - // solver; - solver.analyzePattern(A); - - solver.factorize(A); - - std::cout << A << std::endl; - std::cout << solver.lastErrorMessage() << std::endl; - - x_out = solver.solve(b); - - std::cout << x_out << std::endl; -} diff --git a/src/diffusion.hpp b/src/diffusion.hpp deleted file mode 100644 index 6645f8c..0000000 --- a/src/diffusion.hpp +++ /dev/null @@ -1,11 +0,0 @@ -#ifndef DIFFUSION_H_ -#define DIFFUSION_H_ - -#include -#include - -typedef Eigen::Triplet T; - -extern void BTCS2D(int x, int y, std::vector &c, - std::vector &alpha, double timestep); -#endif // DIFFUSION_H_ diff --git a/src/main.cpp b/src/main.cpp index 8269d9b..2836c7c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,4 +1,4 @@ -#include "diffusion.hpp" +#include "BTCSDiffusion.hpp" #include #include #include @@ -7,16 +7,25 @@ using namespace std; int main(int argc, char *argv[]) { - int x = 5; - int y = 5; + int x = 20; - std::vector alpha(x * y, -1.5 * pow(10, -2)); - std::vector input(x * y, 0); - input[x + 1] = 5.55 * std::pow(10, -6); - input[x + 2] = 5.5556554 * std::pow(10, -6); - input[x + 3] = 5.234564213 * std::pow(10, -6); + std::vector alpha(x, 1 * pow(10, -1)); + std::vector input(x, 1 * std::pow(10, -6)); + std::vector bc_left, bc_right; - BTCS2D(x, y, input, alpha, 10.); + bc_left.push_back(5. * std::pow(10, -6)); + bc_right.push_back(-1); + + BTCSDiffusion diffu(x); + + diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT); + // we don't need this since Neumann condition with gradient of 0 is set per + // default + // diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT); + + for (int i = 0; i < 100; i++) { + diffu.simulate(input, alpha, 1.); + } return 0; }