From 971f8212afaba568a09fa0da92066e9306c4dbc2 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 2 Dec 2021 09:25:34 +0100 Subject: [PATCH] Added comments --- src/BTCSDiffusion.cpp | 32 +++++++++++++++++++++-- src/BTCSDiffusion.hpp | 59 ++++++++++++++++++++++++++++++++++++++++++- src/main.cpp | 4 ++- 3 files changed, 91 insertions(+), 4 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index d2baef1..287fcb7 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -13,6 +13,7 @@ #include #include +#include #include #include @@ -21,18 +22,23 @@ const BCSide BTCSDiffusion::RIGHT = 1; BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) { this->grid_dim = 1; - this->bc.reserve(2); + + // 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 + // TODO: reserve memory for boundary conditions } void BTCSDiffusion::setBoundaryCondition(std::vector input, @@ -43,17 +49,27 @@ void BTCSDiffusion::setBoundaryCondition(std::vector input, } void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, double timestep) { + // calculate dx double dx = 1. / this->dim_x; + // 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); @@ -66,9 +82,14 @@ void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, 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]; @@ -80,6 +101,13 @@ void BTCSDiffusion::simulate(std::vector &c, std::vector &alpha, 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()); diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp index 693b8f7..cadc1e7 100644 --- a/src/BTCSDiffusion.hpp +++ b/src/BTCSDiffusion.hpp @@ -1,23 +1,80 @@ #ifndef BTCSDIFFUSION_H_ #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; +/*! + * 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); diff --git a/src/main.cpp b/src/main.cpp index 03b5152..2836c7c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -19,7 +19,9 @@ int main(int argc, char *argv[]) { BTCSDiffusion diffu(x); diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT); - diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT); + // 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.);