Merge branch 'boundary' into 'dev'

Implemented 1D diffusion with new data structure

See merge request mluebke/diffusion!3

.gitignore

@@ -7,3 +7,4 @@ build/
 compile_commands.json
 .cache/
 .ccls-cache/
+/iwyu/

src/BTCSDiffusion.cpp

@@ -1,129 +1,127 @@
 #include "BTCSDiffusion.hpp"
 
-#include <Eigen/SparseCholesky>
 #include <Eigen/SparseLU>
-#include <Eigen/SparseQR>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
-#include <Eigen/src/OrderingMethods/Ordering.h>
-#include <Eigen/src/SparseCholesky/SimplicialCholesky.h>
-#include <Eigen/src/SparseCore/SparseMap.h>
-#include <Eigen/src/SparseCore/SparseMatrix.h>
-#include <Eigen/src/SparseCore/SparseMatrixBase.h>
-#include <Eigen/src/SparseLU/SparseLU.h>
-#include <Eigen/src/SparseQR/SparseQR.h>
 
 #include <algorithm>
 #include <iomanip>
 #include <iostream>
+#include <tuple>
+#include <vector>
 
-const BCSide BTCSDiffusion::LEFT = 0;
+const int BTCSDiffusion::BC_NEUMANN = 0;
-const BCSide BTCSDiffusion::RIGHT = 1;
+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<bctype, double>(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);
+  // this->bc.reserve(x * 2 + y * 2);
-  // per default use Neumann condition with gradient of 0 at the end of the grid
+  // // per default use Neumann condition with gradient of 0 at the end of the
-  std::fill(this->bc.begin(), this->bc.end(), -1);
+  // grid std::fill(this->bc.begin(), this->bc.end(), -1);
 }
-BTCSDiffusion::BTCSDiffusion(int x, int y, int z)
+BTCSDiffusion::BTCSDiffusion(int x, int y, int z) : n_x(x), n_y(y), n_z(z) {
-    : dim_x(x), dim_y(y), dim_z(z) {
 
-  this->grid_dim = 3;
+  // this->grid_dim = 3;
   // TODO: reserve memory for boundary conditions
 }
 
-void BTCSDiffusion::setBoundaryCondition(std::vector<double> input,
+void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
-                                         BCSide side) {
+                               double bc_right,
-  if (this->grid_dim == 1) {
+                               const std::vector<double> &alpha, double dx,
-    bc[side] = input[0];
+                               int size) {
-  }
-}
-void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
-                             double timestep) {
-  // calculate dx
-  double dx = 1. / (this->dim_x - 1);
 
-  // calculate size needed for A matrix and b,x vectors
+  // we need 2 more grid cells for ghost cells
-  int size = this->dim_x + 2;
+  size = size + 2;
 
-  Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
+  // set sizes of private and yet allocated vectors
-  Eigen::VectorXd x_out(size);
+  b_vector.resize(size);
+  x_vector.resize(size);
 
   /*
-   * Initalization of matrix A
+   * Begin to solve the equation system using LU solver of Eigen.
-   * This is done by triplets. See:
+   *
-   * https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
+   * But first fill the A matrix and b vector.
-   */
-
-  std::vector<T> 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<double> A(size, size);
+  A_matrix.resize(size, size);
-  A.setFromTriplets(tripletList.begin(), tripletList.end());
+  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::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       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<double> &c,
+                             const std::vector<double> &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;
+}

src/BTCSDiffusion.hpp

@@ -2,19 +2,31 @@
 #define BTCSDIFFUSION_H_
 
 #include <Eigen/Sparse>
+#include <tuple>
 #include <vector>
 
-/*!
- * 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<double> 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<std::tuple<bctype, double>> 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);
+  explicit 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<double> input, BCSide side);
 
   /*!
    * 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<double> &c, std::vector<double> &alpha,
+  void simulate(std::vector<double> &c, const std::vector<double> &alpha);
-                double timestep);
+
+  /*!
+   * 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<double> bc;
+  void simulate1D(std::vector<double> &c, double bc_left, double bc_right,
+                  const std::vector<double> &alpha, double dx, int size);
+  void simulate2D(std::vector<double> &c);
+  void simulate3D(std::vector<double> &c);
+
+  double getBCFromTuple(int index, double nearest_value, double neighbor_alpha);
+
+  boundary_condition bc;
+
+  Eigen::SparseMatrix<double> A_matrix;
+  Eigen::VectorXd b_vector;
+  Eigen::VectorXd x_vector;
+
+  double time_step;
+
   int grid_dim;
-  int dim_x;
+  int n_x;
-  int dim_y;
+  double dx;
-  int dim_z;
+  int n_y;
+  double dy;
+  int n_z;
+  double dz;
 };
 
 #endif // BTCSDIFFUSION_H_

src/main.cpp

@@ -1,30 +1,32 @@
-#include "BTCSDiffusion.hpp"
+#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
-#include <cmath>
+#include <algorithm>         // for copy, max
-#include <iostream>
+#include <iostream>          // for std
-#include <vector>
+#include <vector>            // 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<double> alpha(x, 1 * pow(10, -1));
   std::vector<double> input(x, 1 * std::pow(10, -6));
-  std::vector<double> 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);
+  // set the boundary condition for the left ghost cell to dirichlet
-  // we don't need this since Neumann condition with gradient of 0 is set per
+  diffu.setBoundaryCondition(0, 5. * std::pow(10, -6),
-  // default
+                             BTCSDiffusion::BC_DIRICHLET);
-  // diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT);
 
+  // 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;