From 816a1cc67ab8b8bf88332a70d4f2d7f011547ad1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Mon, 13 Dec 2021 19:39:22 +0100 Subject: [PATCH] Implemented 1D diffusion with new data structure --- .gitignore | 1 + src/BTCSDiffusion.cpp | 174 +++++++++++++++++++++--------------------- src/BTCSDiffusion.hpp | 96 +++++++++++++++-------- src/main.cpp | 30 ++++---- 4 files changed, 168 insertions(+), 133 deletions(-) diff --git a/.gitignore b/.gitignore index a072bc9..1e841cf 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,4 @@ build/ compile_commands.json .cache/ .ccls-cache/ +/iwyu/ diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index fcebc92..d5d7ac1 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -1,129 +1,127 @@ #include "BTCSDiffusion.hpp" -#include #include -#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; +const int BTCSDiffusion::BC_NEUMANN = 0; +const int BTCSDiffusion::BC_DIRICHLET = 1; -BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { +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, -1); + this->bc.resize(2, std::tuple(BTCSDiffusion::BC_NEUMANN, 0.)); } -BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) { +BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) { - this->grid_dim = 2; + // 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); + // 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) { +BTCSDiffusion::BTCSDiffusion(int x, int y, int z) : n_x(x), n_y(y), n_z(z) { - this->grid_dim = 3; + // 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); +void BTCSDiffusion::simulate1D(std::vector &c, double bc_left, + double bc_right, + const std::vector &alpha, double dx, + int size) { - // calculate size needed for A matrix and b,x vectors - int size = this->dim_x + 2; + // we need 2 more grid cells for ghost cells + size = size + 2; - Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0); - Eigen::VectorXd x_out(size); + // set sizes of private and yet allocated vectors + b_vector.resize(size); + x_vector.resize(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 + * 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 */ - Eigen::SparseMatrix A(size, size); - A.setFromTriplets(tripletList.begin(), tripletList.end()); + A_matrix.resize(size, size); + A_matrix.reserve(Eigen::VectorXi::Constant(size, 3)); + + A_matrix.insert(0, 0) = 1; + A_matrix.insert(size - 1, size - 1) = 1; + + b_vector[0] = bc_left; + b_vector[size - 1] = bc_right; + + for (int i = 1; i < this->n_x + 1; i++) { + double sx = (alpha[i - 1] * 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]; + } Eigen::SparseLU, Eigen::COLAMDOrdering> solver; - solver.analyzePattern(A); + solver.analyzePattern(A_matrix); - solver.factorize(A); + solver.factorize(A_matrix); std::cout << solver.lastErrorMessage() << std::endl; - x_out = solver.solve(b); + x_vector = solver.solve(b_vector); - std::cout << std::setprecision(10) << x_out << std::endl << std::endl; + std::cout << std::setprecision(10) << x_vector << std::endl << std::endl; for (int i = 0; i < c.size(); i++) { - c[i] = x_out[i + 1]; + c[i] = x_vector[i + 1]; } } + +void BTCSDiffusion::setTimestep(double time_step) { + this->time_step = 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_left, bc_right, alpha, this->dx, this->n_x); + } +} + +double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c, + double neighbor_alpha) { + double val = -1; + int type = std::get<0>(bc[index]); + + 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]); + } else { + // TODO: implement error handling here. Type was set to wrong value. + } + + return val; +} + +void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) { + std::get<0>(bc[index]) = type; + std::get<1>(bc[index]) = val; +} diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index cadc1e7..355c646 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -2,19 +2,31 @@ #define BTCSDIFFUSION_H_ #include +#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; +/*! + * 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. @@ -23,21 +35,20 @@ class BTCSDiffusion { public: /*! - * Set left boundary condition. + * Defines a Neumann boundary condition. */ - static const BCSide LEFT; - + static const int BC_NEUMANN; /*! - * Set right boundary condition. + * Defines a Dirichlet boundary condition. */ - static const BCSide RIGHT; + static const int BC_DIRICHLET; /*! * Create 1D-diffusion module. * * @param x Count of cells in x direction. */ - BTCSDiffusion(int x); + explicit BTCSDiffusion(int x); /*! * Currently not implemented: Create 2D-diffusion module. @@ -45,7 +56,7 @@ public: * @param x Count of cells in x direction. * @param y Count of cells in y direction. */ - BTCSDiffusion(int x, int y); + explicit BTCSDiffusion(int x, int y); /*! * Currently not implemented: Create 3D-diffusion module. @@ -54,18 +65,7 @@ public: * @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); + explicit BTCSDiffusion(int x, int y, int z); /*! * With given ghost zones simulate diffusion. Only 1D allowed at this moment. @@ -73,17 +73,51 @@ public: * @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); + void simulate(std::vector &c, const std::vector &alpha); + + /*! + * Set the timestep of the simulation + * + * @param time_step Time step (in seconds ???) + */ + void setTimestep(double time_step); + + /*! + * Set the boundary condition of the given grid. This is done by defining an + * 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. + */ + void setBoundaryCondition(int index, double val, bctype type); private: - std::vector bc; + void simulate1D(std::vector &c, double bc_left, double bc_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); + + boundary_condition bc; + + Eigen::SparseMatrix A_matrix; + Eigen::VectorXd b_vector; + Eigen::VectorXd x_vector; + + double time_step; + int grid_dim; - int dim_x; - int dim_y; - int dim_z; + int n_x; + double dx; + int n_y; + double dy; + int n_z; + double dz; }; #endif // BTCSDIFFUSION_H_ diff --git a/src/main.cpp b/src/main.cpp index 2836c7c..9735f11 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,30 +1,32 @@ -#include "BTCSDiffusion.hpp" -#include -#include -#include - +#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET +#include // for copy, max +#include // for std +#include // for vector using namespace std; int main(int argc, char *argv[]) { + // count of grid cells int x = 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 bc_left, bc_right; - - bc_left.push_back(5. * std::pow(10, -6)); - bc_right.push_back(-1); + // create instance of diffusion module 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); + // set the boundary condition for the left ghost cell to dirichlet + diffu.setBoundaryCondition(0, 5. * std::pow(10, -6), + BTCSDiffusion::BC_DIRICHLET); + // set timestep for simulation to 1 second + diffu.setTimestep(1.); + + // loop 100 times + // output is currently generated by the method itself for (int i = 0; i < 100; i++) { - diffu.simulate(input, alpha, 1.); + diffu.simulate(input, alpha); } return 0;