feat: Remove class BTCSDiffusion

BREAKING CHANGE: Functionality is now provided by function calls and
scheme generation is decoupled from LEqS solving.

include/BoundaryCondition.hpp

@@ -10,7 +10,7 @@
 typedef uint8_t bctype;
 
 namespace tug {
-namespace boundary_condition {
+namespace bc {
 
 enum {
   BC_TYPE_CLOSED,   /**< Defines a closed/Neumann boundary condition. */

include/Diffusion.hpp

@@ -0,0 +1,72 @@
+#ifndef DIFFUSION_H_
+#define DIFFUSION_H_
+
+#include "BoundaryCondition.hpp"
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <array>
+#include <bits/stdint-uintn.h>
+#include <vector>
+
+namespace tug {
+namespace diffusion {
+
+/**
+ * Defines grid dimensions and boundary conditions.
+ */
+typedef struct {
+  uint32_t grid_cells[3];
+  double domain_size[3];
+  bc::BoundaryCondition *bc;
+} TugGrid;
+
+/**
+ * Besides containing the grid structure it holds also information about the
+ * desired time step to simulate and which solver to use.
+ */
+typedef struct {
+  double time_step;
+  Eigen::VectorXd (*solver)(Eigen::SparseMatrix<double>, Eigen::VectorXd);
+  TugGrid grid;
+} TugInput;
+
+/**
+ * Solving 1D diffusion problems with Backward Time Centred Space scheme.
+ *
+ * \param input_param Object with information for the simulation e.g. grid
+ * dimensions, time step ...
+ *
+ * \param field Pointer to continious memory holding the concentrations for each
+ * grid cell of the grid. It doesn't matter if stored in row (most likely) or
+ * column major.
+ *
+ * \param alpha Pointer to continious memory holding the alpha for each grid
+ * cell of the grid. (NOTE: only constant alphas are supported currently)
+ *
+ * \return Runtime in seconds
+ */
+auto BTCS_1D(const TugInput &input_param, double *field, const double *alpha)
+    -> double;
+
+/**
+ * Solving 2D diffusion problems with Alternating-direction implicit method.
+ *
+ * \param input_param Object with information for the simulation e.g. grid
+ * dimensions, time step ...
+ *
+ * \param field Pointer to continious memory holding the concentrations for each
+ * grid cell of the grid. It doesn't matter if stored in row (most likely) or
+ * column major.
+ *
+ * \param alpha Pointer to continious memory holding the alpha for each grid
+ * cell of the grid. (NOTE: only constant alphas are supported currently)
+ *
+ * \return Runtime in seconds
+ */
+auto ADI_2D(const TugInput &input_param, double *field, const double *alpha)
+    -> double;
+
+} // namespace diffusion
+} // namespace tug
+
+#endif // DIFFUSION_H_

include/Solver.hpp

@@ -0,0 +1,41 @@
+#ifndef SOLVER_H_
+#define SOLVER_H_
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace tug {
+namespace solver {
+
+/**
+ * Solving linear equation system with LU-decomposition implemented in Eigen
+ * library.
+ *
+ * \param A_matrix The A matrix represented as a sparse matrix using Eigen
+ * library.
+ *
+ * \param b_vector Right hand side vector of the linear equation system.
+ *
+ * \return Solution represented as vector.
+ */
+auto EigenLU(const Eigen::SparseMatrix<double> A_matrix,
+             const Eigen::VectorXd b_vector) -> Eigen::VectorXd;
+
+/**
+ * Solving linear equation system with brutal implementation of the Thomas
+ * algorithm.
+ *
+ * \param A_matrix The A matrix represented as a sparse matrix using Eigen
+ * library.
+ *
+ * \param b_vector Right hand side vector of the linear equation system.
+ *
+ * \return Solution represented as vector.
+ */
+auto ThomasAlgorithm(const Eigen::SparseMatrix<double> A_matrix,
+                     const Eigen::VectorXd b_vector) -> Eigen::VectorXd;
+
+} // namespace solver
+} // namespace tug
+
+#endif // SOLVER_H_

include/diffusion/BTCSDiffusion.hpp

@@ -1,168 +0,0 @@
-#ifndef BTCSDIFFUSION_H_
-#define BTCSDIFFUSION_H_
-
-#include "grid/BoundaryCondition.hpp"
-
-#include <Eigen/Sparse>
-#include <Eigen/src/Core/Map.h>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
-#include <Eigen/src/SparseCore/SparseMatrix.h>
-#include <cstddef>
-#include <tuple>
-#include <type_traits>
-#include <vector>
-
-namespace tug {
-namespace diffusion {
-
-/*!
- * Class implementing a solution for a 1/2/3D diffusion equation using backward
- * euler.
- */
-class BTCSDiffusion {
-
-public:
-  /*!
-   * Creates a diffusion module.
-   *
-   * \param dim Number of dimensions. Should not be greater than 3 and not less
-   * than 1.
-   */
-  BTCSDiffusion(unsigned int dim);
-
-  /*!
-   * Define the grid in x direction.
-   *
-   * \param domain_size Size of the domain in x direction.
-   * \param n_grid_cells Number of grid cells in x direction the domain is
-   * divided to.
-   */
-  void setXDimensions(double domain_size, unsigned int n_grid_cells);
-
-  /*!
-   * Define the grid in y direction.
-   *
-   * Throws an error if the module wasn't initialized at least as a 2D model.
-   *
-   * \param domain_size Size of the domain in y direction.
-   * \param n_grid_cells Number of grid cells in y direction the domain is
-   * divided to.
-   */
-  void setYDimensions(double domain_size, unsigned int n_grid_cells);
-
-  /*!
-   * Define the grid in z direction.
-   *
-   * Throws an error if the module wasn't initialized at least as a 3D model.
-   *
-   * \param domain_size Size of the domain in z direction.
-   * \param n_grid_cells Number of grid cells in z direction the domain is
-   * divided to.
-   */
-  void setZDimensions(double domain_size, unsigned int n_grid_cells);
-
-  /*!
-   * Returns the number of grid cells in x direction.
-   */
-  auto getXGridCellsN() -> unsigned int;
-  /*!
-   * Returns the number of grid cells in y direction.
-   */
-  auto getYGridCellsN() -> unsigned int;
-  /*!
-   * Returns the number of grid cells in z direction.
-   */
-  auto getZGridCellsN() -> unsigned int;
-
-  /*!
-   * Returns the domain size in x direction.
-   */
-  auto getXDomainSize() -> double;
-  /*!
-   * Returns the domain size in y direction.
-   */
-  auto getYDomainSize() -> double;
-  /*!
-   * Returns the domain size in z direction.
-   */
-  auto getZDomainSize() -> double;
-
-  /*!
-   * With given ghost zones simulate diffusion. Only 1D allowed at this moment.
-   *
-   * \param c Pointer to continious memory describing the current concentration
-   * state of each grid cell.
-   * \param alpha Pointer to memory area of diffusion coefficients for each grid
-   * element.
-   * \param bc Instance of boundary condition class with.
-   *
-   * \return Time in seconds [s] used to simulate one iteration.
-   */
-  auto simulate(double *c, double *alpha,
-                const tug::boundary_condition::BoundaryCondition &bc)
-      -> double;
-
-  /*!
-   * Set the timestep of the simulation
-   *
-   * \param time_step Time step (in seconds ???)
-   */
-  void setTimestep(double time_step);
-
-private:
-  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
-      DMatrixRowMajor;
-  typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>
-      DVectorRowMajor;
-
-  static void simulate_base(DVectorRowMajor &c,
-                            const tug::boundary_condition::bc_tuple &bc_ghosts,
-                            const tug::boundary_condition::bc_vec &bc_inner,
-                            const DVectorRowMajor &alpha, double dx,
-                            double time_step, int size,
-                            const DVectorRowMajor &d_ortho);
-
-  void simulate1D(Eigen::Map<DVectorRowMajor> &c,
-                  Eigen::Map<const DVectorRowMajor> &alpha,
-                  const tug::boundary_condition::BoundaryCondition &bc);
-
-  void simulate2D(Eigen::Map<DMatrixRowMajor> &c,
-                  Eigen::Map<const DMatrixRowMajor> &alpha,
-                  const tug::boundary_condition::BoundaryCondition &bc);
-
-  static auto
-  calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
-               const tug::boundary_condition::BoundaryCondition &bc,
-               bool transposed, double time_step, double dx) -> DMatrixRowMajor;
-
-  static void fillMatrixFromRow(Eigen::SparseMatrix<double> &A_matrix,
-                                const DVectorRowMajor &alpha,
-                                const tug::boundary_condition::bc_vec &bc_inner,
-                                int size, double dx, double time_step);
-
-  static void
-  fillVectorFromRow(Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
-                    const DVectorRowMajor &alpha,
-                    const tug::boundary_condition::bc_tuple &bc_ghosts,
-                    const tug::boundary_condition::bc_vec &bc_inner,
-                    const DVectorRowMajor &d_ortho, int size, double dx,
-                    double time_step);
-  void simulate3D(std::vector<double> &c);
-
-  inline static auto
-  getBCFromFlux(tug::boundary_condition::boundary_condition bc,
-                double neighbor_c, double neighbor_alpha) -> double;
-
-  void updateInternals();
-
-  double time_step;
-  unsigned int grid_dim;
-
-  std::vector<unsigned int> grid_cells;
-  std::vector<double> domain_size;
-  std::vector<double> deltas;
-};
-} // namespace diffusion
-} // namespace tug
-#endif // BTCSDIFFUSION_H_

src/BTCS.cpp

@@ -0,0 +1,315 @@
+#include <Diffusion.hpp>
+#include <Solver.hpp>
+#include <iostream>
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/src/Core/Matrix.h>
+#include <bits/stdint-uintn.h>
+#include <chrono>
+
+#include "BoundaryCondition.hpp"
+#include "TugUtils.hpp"
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_thread_num() 0
+#endif
+
+#define init_delta(in, out, dim)                                               \
+  ({                                                                           \
+    for (uint8_t i = 0; i < dim; i++) {                                        \
+      out[i] = (double)in.domain_size[i] / in.grid_cells[i];                   \
+    }                                                                          \
+  })
+
+namespace {
+enum { GRID_1D = 1, GRID_2D, GRID_3D };
+
+constexpr int BTCS_MAX_DEP_PER_CELL = 3;
+constexpr int BTCS_2D_DT_SIZE = 2;
+
+typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
+    DMatrixRowMajor;
+typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>
+    DVectorRowMajor;
+
+inline auto getBCFromFlux(tug::bc::boundary_condition bc, double neighbor_c,
+                          double neighbor_alpha) -> double {
+  double val = 0;
+
+  if (bc.type == tug::bc::BC_TYPE_CLOSED) {
+    val = neighbor_c;
+  } else if (bc.type == tug::bc::BC_TYPE_FLUX) {
+    // TODO
+    //  val = bc[index].value;
+  } else {
+    // TODO: implement error handling here. Type was set to wrong value.
+  }
+
+  return val;
+}
+
+auto calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
+                  const tug::bc::BoundaryCondition &bc, bool transposed,
+                  double time_step, double dx) -> DMatrixRowMajor {
+
+  uint8_t upper = (transposed ? tug::bc::BC_SIDE_LEFT : tug::bc::BC_SIDE_TOP);
+  uint8_t lower =
+      (transposed ? tug::bc::BC_SIDE_RIGHT : tug::bc::BC_SIDE_BOTTOM);
+
+  int n_rows = c.rows();
+  int n_cols = c.cols();
+
+  DMatrixRowMajor d_ortho(n_rows, n_cols);
+
+  std::array<double, 3> y_values{};
+
+  // first, iterate over first row
+  for (int j = 0; j < n_cols; j++) {
+    tug::bc::boundary_condition tmp_bc = bc(upper, j);
+    double sy = (time_step * alpha(0, j)) / (dx * dx);
+
+    y_values[0] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT
+                       ? tmp_bc.value
+                       : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j)));
+    y_values[1] = c(0, j);
+    y_values[2] = c(1, j);
+
+    d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]);
+  }
+
+  // then iterate over inlet
+#pragma omp parallel for private(y_values) schedule(dynamic)
+  for (int i = 1; i < n_rows - 1; i++) {
+    for (int j = 0; j < n_cols; j++) {
+      double sy = (time_step * alpha(i, j)) / (dx * dx);
+
+      y_values[0] = c(i - 1, j);
+      y_values[1] = c(i, j);
+      y_values[2] = c(i + 1, j);
+
+      d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]);
+    }
+  }
+
+  int end = n_rows - 1;
+
+  // and finally over last row
+  for (int j = 0; j < n_cols; j++) {
+    tug::bc::boundary_condition tmp_bc = bc(lower, j);
+    double sy = (time_step * alpha(end, j)) / (dx * dx);
+
+    y_values[0] = c(end - 1, j);
+    y_values[1] = c(end, j);
+    y_values[2] = (tmp_bc.type == tug::bc::BC_TYPE_CONSTANT
+                       ? tmp_bc.value
+                       : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j)));
+
+    d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]);
+  }
+
+  return d_ortho;
+}
+
+auto fillMatrixFromRow(const DVectorRowMajor &alpha,
+                       const tug::bc::bc_vec &bc_inner, int size, double dx,
+                       double time_step) -> Eigen::SparseMatrix<double> {
+
+  Eigen::SparseMatrix<double> A_matrix(size + 2, size + 2);
+
+  A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL));
+
+  double sx = 0;
+
+  int A_size = A_matrix.cols();
+
+  A_matrix.insert(0, 0) = 1;
+
+  if (bc_inner[0].type != tug::bc::BC_UNSET) {
+    if (bc_inner[0].type != tug::bc::BC_TYPE_CONSTANT) {
+      throw_invalid_argument("Inner boundary conditions with other type than "
+                             "BC_TYPE_CONSTANT are currently not supported.");
+    }
+    A_matrix.insert(1, 1) = 1;
+  } else {
+    sx = (alpha[0] * time_step) / (dx * dx);
+    A_matrix.insert(1, 1) = -1. - 3. * sx;
+    A_matrix.insert(1, 0) = 2. * sx;
+    A_matrix.insert(1, 2) = sx;
+  }
+
+  for (int j = 2, k = j - 1; k < size - 1; j++, k++) {
+    if (bc_inner[k].type != tug::bc::BC_UNSET) {
+      if (bc_inner[k].type != tug::bc::BC_TYPE_CONSTANT) {
+        throw_invalid_argument("Inner boundary conditions with other type than "
+                               "BC_TYPE_CONSTANT are currently not supported.");
+      }
+      A_matrix.insert(j, j) = 1;
+      continue;
+    }
+    sx = (alpha[k] * time_step) / (dx * dx);
+
+    A_matrix.insert(j, j) = -1. - 2. * sx;
+    A_matrix.insert(j, (j - 1)) = sx;
+    A_matrix.insert(j, (j + 1)) = sx;
+  }
+
+  if (bc_inner[size - 1].type != tug::bc::BC_UNSET) {
+    if (bc_inner[size - 1].type != tug::bc::BC_TYPE_CONSTANT) {
+      throw_invalid_argument("Inner boundary conditions with other type than "
+                             "BC_TYPE_CONSTANT are currently not supported.");
+    }
+    A_matrix.insert(A_size - 2, A_size - 2) = 1;
+  } else {
+    sx = (alpha[size - 1] * time_step) / (dx * dx);
+    A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
+    A_matrix.insert(A_size - 2, A_size - 3) = sx;
+    A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
+  }
+
+  A_matrix.insert(A_size - 1, A_size - 1) = 1;
+
+  return A_matrix;
+}
+
+auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha,
+                       const tug::bc::bc_tuple &bc,
+                       const tug::bc::bc_vec &bc_inner,
+                       const DVectorRowMajor &d_ortho, int size, double dx,
+                       double time_step) -> Eigen::VectorXd {
+
+  Eigen::VectorXd b_vector(size + 2);
+
+  tug::bc::boundary_condition left = bc[0];
+  tug::bc::boundary_condition right = bc[1];
+
+  bool left_constant = (left.type == tug::bc::BC_TYPE_CONSTANT);
+  bool right_constant = (right.type == tug::bc::BC_TYPE_CONSTANT);
+
+  int b_size = b_vector.size();
+
+  for (int j = 0; j < size; j++) {
+    if (bc_inner[j].type != tug::bc::BC_UNSET) {
+      if (bc_inner[j].type != tug::bc::BC_TYPE_CONSTANT)
+        throw_invalid_argument("Inner boundary conditions with other type than "
+                               "BC_TYPE_CONSTANT are currently not supported.");
+      b_vector[j + 1] = bc_inner[j].value;
+      continue;
+    }
+    b_vector[j + 1] = -c[j] + d_ortho[j];
+  }
+
+  // this is not correct currently.We will fix this when we are able to define
+  // FLUX boundary conditions
+  b_vector[0] =
+      (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
+
+  b_vector[b_size - 1] =
+      (right_constant ? right.value
+                      : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
+
+  return b_vector;
+}
+
+auto setupBTCSAndSolve(
+    DVectorRowMajor &c, const tug::bc::bc_tuple bc_ghost,
+    const tug::bc::bc_vec bc_inner, const DVectorRowMajor &alpha, double dx,
+    double time_step, int size, const DVectorRowMajor &d_ortho,
+    Eigen::VectorXd (*solver)(Eigen::SparseMatrix<double>, Eigen::VectorXd))
+    -> DVectorRowMajor {
+
+  const Eigen::SparseMatrix<double> A_matrix =
+      fillMatrixFromRow(alpha, bc_inner, size, dx, time_step);
+
+  const Eigen::VectorXd b_vector = fillVectorFromRow(
+      c, alpha, bc_ghost, bc_inner, d_ortho, size, dx, time_step);
+
+  // solving of the LEQ
+  Eigen::VectorXd x_vector = solver(A_matrix, b_vector);
+
+  DVectorRowMajor out_vector = x_vector.segment(1, size);
+
+  return out_vector;
+}
+
+} // namespace
+  //
+auto tug::diffusion::BTCS_1D(const tug::diffusion::TugInput &input_param,
+                             double *field, const double *alpha) -> double {
+
+  auto start = time_marker();
+
+  const tug::bc::BoundaryCondition bc = *input_param.grid.bc;
+  double deltas[GRID_1D];
+
+  init_delta(input_param.grid, deltas, GRID_1D);
+
+  uint32_t size = input_param.grid.grid_cells[0];
+  double dx = deltas[0];
+  double time_step = input_param.time_step;
+
+  Eigen::Map<DVectorRowMajor> c_in(field, size);
+  Eigen::Map<const DVectorRowMajor> alpha_in(alpha, size);
+
+  DVectorRowMajor input_field = c_in.row(0);
+
+  DVectorRowMajor output = setupBTCSAndSolve(
+      input_field, bc.row_boundary(0), bc.getInnerRow(0), alpha_in, dx,
+      time_step, size, Eigen::VectorXd::Constant(size, 0), input_param.solver);
+
+  c_in.row(0) << output;
+
+  auto end = time_marker();
+
+  return diff_time(start, end);
+}
+
+auto tug::diffusion::ADI_2D(const tug::diffusion::TugInput &input_param,
+                            double *field, const double *alpha) -> double {
+
+  auto start = time_marker();
+
+  tug::bc::BoundaryCondition bc = *input_param.grid.bc;
+  double deltas[GRID_2D];
+
+  init_delta(input_param.grid, deltas, GRID_2D);
+
+  uint32_t n_cols = input_param.grid.grid_cells[0];
+  uint32_t n_rows = input_param.grid.grid_cells[1];
+  double dx = deltas[0];
+  double dy = deltas[1];
+  double local_dt = input_param.time_step / BTCS_2D_DT_SIZE;
+
+  Eigen::Map<DMatrixRowMajor> c_in(field, n_rows, n_cols);
+  Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, n_rows, n_cols);
+
+  DMatrixRowMajor d_ortho =
+      calc_d_ortho(c_in, alpha_in, bc, false, local_dt, dx);
+
+#pragma omp parallel for schedule(dynamic)
+  for (int i = 0; i < n_rows; i++) {
+    DVectorRowMajor input_field = c_in.row(i);
+    DVectorRowMajor output = setupBTCSAndSolve(
+        input_field, bc.row_boundary(i), bc.getInnerRow(i), alpha_in.row(i), dx,
+        local_dt, n_cols, d_ortho.row(i), input_param.solver);
+    c_in.row(i) << output;
+  }
+
+  d_ortho = calc_d_ortho(c_in.transpose(), alpha_in.transpose(), bc, true,
+                         local_dt, dy);
+
+#pragma omp parallel for schedule(dynamic)
+  for (int i = 0; i < n_cols; i++) {
+    DVectorRowMajor input_field = c_in.col(i);
+    DVectorRowMajor output = setupBTCSAndSolve(
+        input_field, bc.col_boundary(i), bc.getInnerCol(i), alpha_in.col(i), dy,
+        local_dt, n_rows, d_ortho.row(i), input_param.solver);
+    c_in.col(i) << output.transpose();
+  }
+
+  auto end = time_marker();
+
+  return diff_time(start, end);
+}

src/BTCSDiffusion.cpp

@@ -1,381 +0,0 @@
-#include "diffusion/BTCSDiffusion.hpp"
-#include "BTCSUtils.hpp"
-#include "grid/BoundaryCondition.hpp"
-
-#include <Eigen/SparseLU>
-
-#include <Eigen/src/Core/Map.h>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/SparseCore/SparseMatrix.h>
-#include <Eigen/src/SparseCore/SparseMatrixBase.h>
-#include <algorithm>
-#include <array>
-#include <bits/stdint-uintn.h>
-#include <cassert>
-#include <chrono>
-#include <cstddef>
-#include <iomanip>
-#include <iostream>
-#include <iterator>
-#include <ostream>
-#include <tuple>
-#include <vector>
-
-#ifdef _OPENMP
-#include <omp.h>
-#else
-#define omp_get_thread_num() 0
-#endif
-
-constexpr double DOUBLE_MACHINE_EPSILON = 1.93e-34;
-
-constexpr int BTCS_MAX_DEP_PER_CELL = 3;
-constexpr int BTCS_2D_DT_SIZE = 2;
-tug::diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
-
-  grid_cells.resize(dim, 1);
-  domain_size.resize(dim, 1);
-  deltas.resize(dim, 1);
-
-  this->time_step = 0;
-}
-
-void tug::diffusion::BTCSDiffusion::setXDimensions(double domain_size,
-                                                   unsigned int n_grid_cells) {
-  this->domain_size[0] = domain_size;
-  this->grid_cells[0] = n_grid_cells;
-
-  updateInternals();
-}
-
-void tug::diffusion::BTCSDiffusion::setYDimensions(double domain_size,
-                                                   unsigned int n_grid_cells) {
-  this->domain_size[1] = domain_size;
-  this->grid_cells[1] = n_grid_cells;
-
-  updateInternals();
-}
-
-void tug::diffusion::BTCSDiffusion::setZDimensions(double domain_size,
-                                                   unsigned int n_grid_cells) {
-  this->domain_size[2] = domain_size;
-  this->grid_cells[2] = n_grid_cells;
-
-  updateInternals();
-}
-
-auto tug::diffusion::BTCSDiffusion::getXGridCellsN() -> unsigned int {
-  return this->grid_cells[0];
-}
-auto tug::diffusion::BTCSDiffusion::getYGridCellsN() -> unsigned int {
-  return this->grid_cells[1];
-}
-auto tug::diffusion::BTCSDiffusion::getZGridCellsN() -> unsigned int {
-  return this->grid_cells[2];
-}
-auto tug::diffusion::BTCSDiffusion::getXDomainSize() -> double {
-  return this->domain_size[0];
-}
-auto tug::diffusion::BTCSDiffusion::getYDomainSize() -> double {
-  return this->domain_size[1];
-}
-auto tug::diffusion::BTCSDiffusion::getZDomainSize() -> double {
-  return this->domain_size[2];
-}
-
-void tug::diffusion::BTCSDiffusion::updateInternals() {
-  for (int i = 0; i < grid_dim; i++) {
-    deltas[i] = (double)domain_size[i] / grid_cells[i];
-  }
-}
-void tug::diffusion::BTCSDiffusion::simulate_base(
-    DVectorRowMajor &c, const tug::boundary_condition::bc_tuple &bc_ghosts,
-    const tug::boundary_condition::bc_vec &bc_inner,
-    const DVectorRowMajor &alpha, double dx, double time_step, int size,
-    const DVectorRowMajor &d_ortho) {
-
-  Eigen::SparseMatrix<double> A_matrix;
-  Eigen::VectorXd b_vector;
-  Eigen::VectorXd x_vector;
-
-  A_matrix.resize(size + 2, size + 2);
-  A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, BTCS_MAX_DEP_PER_CELL));
-
-  b_vector.resize(size + 2);
-  x_vector.resize(size + 2);
-
-  fillMatrixFromRow(A_matrix, alpha, bc_inner, size, dx, time_step);
-  fillVectorFromRow(b_vector, c, alpha, bc_ghosts, bc_inner, d_ortho, size, dx,
-                    time_step);
-
-  // start to solve
-  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
-      solver;
-  solver.analyzePattern(A_matrix);
-
-  solver.factorize(A_matrix);
-
-  x_vector = solver.solve(b_vector);
-
-  c = x_vector.segment(1, size);
-}
-
-void tug::diffusion::BTCSDiffusion::simulate1D(
-    Eigen::Map<DVectorRowMajor> &c, Eigen::Map<const DVectorRowMajor> &alpha,
-    const tug::boundary_condition::BoundaryCondition &bc) {
-
-  int size = this->grid_cells[0];
-  double dx = this->deltas[0];
-  double time_step = this->time_step;
-
-  DVectorRowMajor input_field = c.row(0);
-
-  simulate_base(input_field, bc.row_boundary(0), bc.getInnerRow(0), alpha, dx,
-                time_step, size, Eigen::VectorXd::Constant(size, 0));
-
-  c.row(0) << input_field;
-}
-
-void tug::diffusion::BTCSDiffusion::simulate2D(
-    Eigen::Map<DMatrixRowMajor> &c, Eigen::Map<const DMatrixRowMajor> &alpha,
-    const tug::boundary_condition::BoundaryCondition &bc) {
-
-  int n_rows = this->grid_cells[1];
-  int n_cols = this->grid_cells[0];
-  double dx = this->deltas[0];
-
-  double local_dt = this->time_step / BTCS_2D_DT_SIZE;
-
-  DMatrixRowMajor d_ortho = calc_d_ortho(c, alpha, bc, false, local_dt, dx);
-
-#pragma omp parallel for schedule(dynamic)
-  for (int i = 0; i < n_rows; i++) {
-    DVectorRowMajor input_field = c.row(i);
-    simulate_base(input_field, bc.row_boundary(i), bc.getInnerRow(i),
-                  alpha.row(i), dx, local_dt, n_cols, d_ortho.row(i));
-    c.row(i) << input_field;
-  }
-
-  dx = this->deltas[1];
-
-  d_ortho =
-      calc_d_ortho(c.transpose(), alpha.transpose(), bc, true, local_dt, dx);
-
-#pragma omp parallel for schedule(dynamic)
-  for (int i = 0; i < n_cols; i++) {
-    DVectorRowMajor input_field = c.col(i);
-    simulate_base(input_field, bc.col_boundary(i), bc.getInnerCol(i),
-                  alpha.col(i), dx, local_dt, n_rows, d_ortho.row(i));
-    c.col(i) << input_field.transpose();
-  }
-}
-
-auto tug::diffusion::BTCSDiffusion::calc_d_ortho(
-    const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
-    const tug::boundary_condition::BoundaryCondition &bc, bool transposed,
-    double time_step, double dx) -> DMatrixRowMajor {
-
-  uint8_t upper = (transposed ? tug::boundary_condition::BC_SIDE_LEFT
-                              : tug::boundary_condition::BC_SIDE_TOP);
-  uint8_t lower = (transposed ? tug::boundary_condition::BC_SIDE_RIGHT
-                              : tug::boundary_condition::BC_SIDE_BOTTOM);
-
-  int n_rows = c.rows();
-  int n_cols = c.cols();
-
-  DMatrixRowMajor d_ortho(n_rows, n_cols);
-
-  std::array<double, 3> y_values{};
-
-  // first, iterate over first row
-  for (int j = 0; j < n_cols; j++) {
-    tug::boundary_condition::boundary_condition tmp_bc = bc(upper, j);
-    double sy = (time_step * alpha(0, j)) / (dx * dx);
-
-    y_values[0] = (tmp_bc.type == tug::boundary_condition::BC_TYPE_CONSTANT
-                       ? tmp_bc.value
-                       : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j)));
-    y_values[1] = c(0, j);
-    y_values[2] = c(1, j);
-
-    d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]);
-  }
-
-// then iterate over inlet
-#pragma omp parallel for private(y_values) schedule(dynamic)
-  for (int i = 1; i < n_rows - 1; i++) {
-    for (int j = 0; j < n_cols; j++) {
-      double sy = (time_step * alpha(i, j)) / (dx * dx);
-
-      y_values[0] = c(i - 1, j);
-      y_values[1] = c(i, j);
-      y_values[2] = c(i + 1, j);
-
-      d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]);
-    }
-  }
-
-  int end = n_rows - 1;
-
-  // and finally over last row
-  for (int j = 0; j < n_cols; j++) {
-    tug::boundary_condition::boundary_condition tmp_bc = bc(lower, j);
-    double sy = (time_step * alpha(end, j)) / (dx * dx);
-
-    y_values[0] = c(end - 1, j);
-    y_values[1] = c(end, j);
-    y_values[2] = (tmp_bc.type == tug::boundary_condition::BC_TYPE_CONSTANT
-                       ? tmp_bc.value
-                       : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j)));
-
-    d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]);
-  }
-
-  return d_ortho;
-}
-
-void tug::diffusion::BTCSDiffusion::fillMatrixFromRow(
-    Eigen::SparseMatrix<double> &A_matrix, const DVectorRowMajor &alpha,
-    const tug::boundary_condition::bc_vec &bc_inner, int size, double dx,
-    double time_step) {
-
-  double sx = 0;
-
-  int A_size = A_matrix.cols();
-
-  A_matrix.insert(0, 0) = 1;
-
-  if (bc_inner[0].type != tug::boundary_condition::BC_UNSET) {
-    if (bc_inner[0].type != tug::boundary_condition::BC_TYPE_CONSTANT) {
-      throw_invalid_argument("Inner boundary conditions with other type than "
-                             "BC_TYPE_CONSTANT are currently not supported.");
-    }
-    A_matrix.insert(1, 1) = 1;
-  } else {
-    sx = (alpha[0] * time_step) / (dx * dx);
-    A_matrix.insert(1, 1) = -1. - 3. * sx;
-    A_matrix.insert(1, 0) = 2. * sx;
-    A_matrix.insert(1, 2) = sx;
-  }
-
-  for (int j = 2, k = j - 1; k < size - 1; j++, k++) {
-    if (bc_inner[k].type != tug::boundary_condition::BC_UNSET) {
-      if (bc_inner[k].type != tug::boundary_condition::BC_TYPE_CONSTANT) {
-        throw_invalid_argument("Inner boundary conditions with other type than "
-                               "BC_TYPE_CONSTANT are currently not supported.");
-      }
-      A_matrix.insert(j, j) = 1;
-      continue;
-    }
-    sx = (alpha[k] * time_step) / (dx * dx);
-
-    A_matrix.insert(j, j) = -1. - 2. * sx;
-    A_matrix.insert(j, (j - 1)) = sx;
-    A_matrix.insert(j, (j + 1)) = sx;
-  }
-
-  if (bc_inner[size - 1].type != tug::boundary_condition::BC_UNSET) {
-    if (bc_inner[size - 1].type != tug::boundary_condition::BC_TYPE_CONSTANT) {
-      throw_invalid_argument("Inner boundary conditions with other type than "
-                             "BC_TYPE_CONSTANT are currently not supported.");
-    }
-    A_matrix.insert(A_size - 2, A_size - 2) = 1;
-  } else {
-    sx = (alpha[size - 1] * time_step) / (dx * dx);
-    A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
-    A_matrix.insert(A_size - 2, A_size - 3) = sx;
-    A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
-  }
-
-  A_matrix.insert(A_size - 1, A_size - 1) = 1;
-}
-
-void tug::diffusion::BTCSDiffusion::fillVectorFromRow(
-    Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
-    const DVectorRowMajor &alpha, const tug::boundary_condition::bc_tuple &bc,
-    const tug::boundary_condition::bc_vec &bc_inner,
-    const DVectorRowMajor &d_ortho, int size, double dx, double time_step) {
-
-  tug::boundary_condition::boundary_condition left = bc[0];
-  tug::boundary_condition::boundary_condition right = bc[1];
-
-  bool left_constant = (left.type == tug::boundary_condition::BC_TYPE_CONSTANT);
-  bool right_constant =
-      (right.type == tug::boundary_condition::BC_TYPE_CONSTANT);
-
-  int b_size = b_vector.size();
-
-  for (int j = 0; j < size; j++) {
-    if (bc_inner[j].type != tug::boundary_condition::BC_UNSET) {
-      if (bc_inner[j].type != tug::boundary_condition::BC_TYPE_CONSTANT)
-        throw_invalid_argument("Inner boundary conditions with other type than "
-                               "BC_TYPE_CONSTANT are currently not supported.");
-      b_vector[j + 1] = bc_inner[j].value;
-      continue;
-    }
-    b_vector[j + 1] = -c[j] + d_ortho[j];
-  }
-
-  // this is not correct currently.We will fix this when we are able to define
-  // FLUX boundary conditions
-  b_vector[0] =
-      (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
-
-  b_vector[b_size - 1] =
-      (right_constant ? right.value
-                      : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
-}
-
-void tug::diffusion::BTCSDiffusion::setTimestep(double time_step) {
-  this->time_step = time_step;
-}
-
-auto tug::diffusion::BTCSDiffusion::simulate(
-    double *c, double *alpha,
-    const tug::boundary_condition::BoundaryCondition &bc) -> double {
-
-  std::chrono::high_resolution_clock::time_point start =
-      std::chrono::high_resolution_clock::now();
-
-  if (this->grid_dim == 1) {
-    Eigen::Map<DVectorRowMajor> c_in(c, this->grid_cells[0]);
-    Eigen::Map<const DVectorRowMajor> alpha_in(alpha, this->grid_cells[0]);
-
-    simulate1D(c_in, alpha_in, bc);
-  }
-  if (this->grid_dim == 2) {
-    Eigen::Map<DMatrixRowMajor> c_in(c, this->grid_cells[1],
-                                     this->grid_cells[0]);
-
-    Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, this->grid_cells[1],
-                                               this->grid_cells[0]);
-
-    simulate2D(c_in, alpha_in, bc);
-  }
-
-  std::chrono::high_resolution_clock::time_point end =
-      std::chrono::high_resolution_clock::now();
-
-  std::chrono::duration<double> duration =
-      std::chrono::duration_cast<std::chrono::duration<double>>(end - start);
-
-  return duration.count();
-}
-
-inline auto tug::diffusion::BTCSDiffusion::getBCFromFlux(
-    tug::boundary_condition::boundary_condition bc, double neighbor_c,
-    double neighbor_alpha) -> double {
-
-  double val = 0;
-
-  if (bc.type == tug::boundary_condition::BC_TYPE_CLOSED) {
-    val = neighbor_c;
-  } else if (bc.type == tug::boundary_condition::BC_TYPE_FLUX) {
-    // TODO
-    //  val = bc[index].value;
-  } else {
-    // TODO: implement error handling here. Type was set to wrong value.
-  }
-
-  return val;
-}

src/BoundaryCondition.cpp

@@ -2,13 +2,13 @@
 #include <bits/stdint-uintn.h>
 #include <vector>
 
-#include "../BTCSUtils.hpp"
+#include "BoundaryCondition.hpp"
-#include "grid/BoundaryCondition.hpp"
+#include "TugUtils.hpp"
 
 constexpr uint8_t DIM_1D = 2;
 constexpr uint8_t DIM_2D = 4;
 
-tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x) {
+tug::bc::BoundaryCondition::BoundaryCondition(int x) {
   this->bc_internal.resize(DIM_1D, {0, 0});
   this->special_cells.resize(x, {BC_UNSET, 0});
   this->dim = 1;
@@ -21,7 +21,7 @@ tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x) {
   this->maxindex = x - 1;
 }
 
-tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x, int y) {
+tug::bc::BoundaryCondition::BoundaryCondition(int x, int y) {
   this->maxsize = (x >= y ? x : y);
   this->bc_internal.resize(DIM_2D * maxsize, {0, 0});
   this->special_cells.resize(x * y, {BC_UNSET, 0});
@@ -33,8 +33,8 @@ tug::boundary_condition::BoundaryCondition::BoundaryCondition(int x, int y) {
   this->maxindex = (x * y) - 1;
 }
 
-void tug::boundary_condition::BoundaryCondition::setSide(
+void tug::bc::BoundaryCondition::setSide(
-    uint8_t side, tug::boundary_condition::boundary_condition &input_bc) {
+    uint8_t side, tug::bc::boundary_condition &input_bc) {
   if (this->dim == 1) {
     throw_invalid_argument("setSide requires at least a 2D grid");
   }
@@ -42,19 +42,18 @@ void tug::boundary_condition::BoundaryCondition::setSide(
     throw_out_of_range("Invalid range for 2D grid");
   }
 
-  uint32_t size = (side == tug::boundary_condition::BC_SIDE_LEFT ||
+  uint32_t size =
-                           side == tug::boundary_condition::BC_SIDE_RIGHT
+      (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT
-                       ? this->sizes[Y_DIM]
+           ? this->sizes[Y_DIM]
-                       : this->sizes[X_DIM]);
+           : this->sizes[X_DIM]);
 
   for (uint32_t i = 0; i < size; i++) {
     this->bc_internal[side * maxsize + i] = input_bc;
   }
 }
 
-void tug::boundary_condition::BoundaryCondition::setSide(
+void tug::bc::BoundaryCondition::setSide(
-    uint8_t side,
+    uint8_t side, std::vector<tug::bc::boundary_condition> &input_bc) {
-    std::vector<tug::boundary_condition::boundary_condition> &input_bc) {
   if (this->dim == 1) {
     throw_invalid_argument("setSide requires at least a 2D grid");
   }
@@ -62,10 +61,10 @@ void tug::boundary_condition::BoundaryCondition::setSide(
     throw_out_of_range("Invalid range for 2D grid");
   }
 
-  uint32_t size = (side == tug::boundary_condition::BC_SIDE_LEFT ||
+  uint32_t size =
-                           side == tug::boundary_condition::BC_SIDE_RIGHT
+      (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT
-                       ? this->sizes[Y_DIM]
+           ? this->sizes[Y_DIM]
-                       : this->sizes[X_DIM]);
+           : this->sizes[X_DIM]);
 
   if (input_bc.size() > size) {
     throw_out_of_range("Input vector is greater than maximum excpected value");
@@ -76,8 +75,8 @@ void tug::boundary_condition::BoundaryCondition::setSide(
   }
 }
 
-auto tug::boundary_condition::BoundaryCondition::getSide(uint8_t side)
+auto tug::bc::BoundaryCondition::getSide(uint8_t side)
-    -> std::vector<tug::boundary_condition::boundary_condition> {
+    -> std::vector<tug::bc::boundary_condition> {
   if (this->dim == 1) {
     throw_invalid_argument("getSide requires at least a 2D grid");
   }
@@ -85,12 +84,12 @@ auto tug::boundary_condition::BoundaryCondition::getSide(uint8_t side)
     throw_out_of_range("Invalid range for 2D grid");
   }
 
-  uint32_t size = (side == tug::boundary_condition::BC_SIDE_LEFT ||
+  uint32_t size =
-                           side == tug::boundary_condition::BC_SIDE_RIGHT
+      (side == tug::bc::BC_SIDE_LEFT || side == tug::bc::BC_SIDE_RIGHT
-                       ? this->sizes[Y_DIM]
+           ? this->sizes[Y_DIM]
-                       : this->sizes[X_DIM]);
+           : this->sizes[X_DIM]);
 
-  std::vector<tug::boundary_condition::boundary_condition> out(size);
+  std::vector<tug::bc::boundary_condition> out(size);
 
   for (int i = 0; i < size; i++) {
     out[i] = this->bc_internal[this->maxsize * side + i];
@@ -99,8 +98,18 @@ auto tug::boundary_condition::BoundaryCondition::getSide(uint8_t side)
   return out;
 }
 
-auto tug::boundary_condition::BoundaryCondition::col_boundary(uint32_t i) const
+auto tug::bc::BoundaryCondition::row_boundary(uint32_t i) const
-    -> tug::boundary_condition::bc_tuple {
+    -> tug::bc::bc_tuple {
+  if (i >= this->sizes[Y_DIM]) {
+    throw_out_of_range("Index out of range");
+  }
+
+  return {this->bc_internal[BC_SIDE_LEFT * this->maxsize + i],
+          this->bc_internal[BC_SIDE_RIGHT * this->maxsize + i]};
+}
+
+auto tug::bc::BoundaryCondition::col_boundary(uint32_t i) const
+    -> tug::bc::bc_tuple {
   if (this->dim == 1) {
     throw_invalid_argument("Access of column requires at least 2D grid");
   }
@@ -112,18 +121,7 @@ auto tug::boundary_condition::BoundaryCondition::col_boundary(uint32_t i) const
           this->bc_internal[BC_SIDE_BOTTOM * this->maxsize + i]};
 }
 
-auto tug::boundary_condition::BoundaryCondition::row_boundary(uint32_t i) const
+auto tug::bc::BoundaryCondition::getInnerRow(uint32_t i) const -> bc_vec {
-    -> tug::boundary_condition::bc_tuple {
-  if (i >= this->sizes[Y_DIM]) {
-    throw_out_of_range("Index out of range");
-  }
-
-  return {this->bc_internal[BC_SIDE_LEFT * this->maxsize + i],
-          this->bc_internal[BC_SIDE_RIGHT * this->maxsize + i]};
-}
-
-auto tug::boundary_condition::BoundaryCondition::getInnerRow(uint32_t i) const
-    -> bc_vec {
   if (i >= this->sizes[Y_DIM]) {
     throw_out_of_range("Index is out of range");
   }
@@ -136,8 +134,7 @@ auto tug::boundary_condition::BoundaryCondition::getInnerRow(uint32_t i) const
   return row;
 }
 
-auto tug::boundary_condition::BoundaryCondition::getInnerCol(uint32_t i) const
+auto tug::bc::BoundaryCondition::getInnerCol(uint32_t i) const -> bc_vec {
-    -> bc_vec {
   if (this->dim != 2) {
     throw_invalid_argument("getInnerCol is only applicable for 2D grids");
   }

src/CMakeLists.txt

@@ -1,4 +1,4 @@
-add_library(tug BTCSDiffusion.cpp grid/BoundaryCondition.cpp)
+add_library(tug BTCS.cpp BoundaryCondition.cpp Solver.cpp)
 
 target_link_libraries(tug Eigen3::Eigen)
 

src/Solver.cpp

@@ -0,0 +1,62 @@
+#include <bits/stdint-uintn.h>
+#include <Solver.hpp>
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+auto tug::solver::EigenLU(const Eigen::SparseMatrix<double> A_matrix,
+                          const Eigen::VectorXd b_vector) -> Eigen::VectorXd {
+
+  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
+      solver;
+  solver.analyzePattern(A_matrix);
+
+  solver.factorize(A_matrix);
+
+  Eigen::VectorXd x_vec = solver.solve(b_vector);
+
+  return solver.solve(b_vector);
+}
+
+auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix<double> A_matrix,
+                                  const Eigen::VectorXd b_vector)
+    -> Eigen::VectorXd {
+  uint32_t n = b_vector.size();
+
+  Eigen::VectorXd a_diag(n);
+  Eigen::VectorXd b_diag(n);
+  Eigen::VectorXd c_diag(n);
+  Eigen::VectorXd x_vec = b_vector;
+
+  // Fill diagonals vectors
+  b_diag[0] = A_matrix.coeff(0, 0);
+  c_diag[0] = A_matrix.coeff(0, 1);
+
+  for (int i = 1; i < n - 1; i++) {
+    a_diag[i] = A_matrix.coeff(i, i - 1);
+    b_diag[i] = A_matrix.coeff(i, i);
+    c_diag[i] = A_matrix.coeff(i, i + 1);
+  }
+  a_diag[n - 1] = A_matrix.coeff(n - 1, n - 2);
+  b_diag[n - 1] = A_matrix.coeff(n - 1, n - 1);
+
+  // start solving - c_diag and x_vec are overwritten
+  n--;
+  c_diag[0] /= b_diag[0];
+  x_vec[0] /= b_diag[0];
+
+  for (int i = 1; i < n; i++) {
+    c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1];
+    x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) /
+               (b_diag[i] - a_diag[i] * c_diag[i - 1]);
+  }
+
+  x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) /
+             (b_diag[n] - a_diag[n] * c_diag[n - 1]);
+
+  for (int i = n; i-- > 0;) {
+    x_vec[i] -= c_diag[i] * x_vec[i + 1];
+  }
+
+  return x_vec;
+}

src/TugUtils.hpp

@@ -1,6 +1,7 @@
 #ifndef BTCSUTILS_H_
 #define BTCSUTILS_H_
 
+#include <chrono>
 #include <stdexcept>
 #include <string>
 
@@ -12,4 +13,14 @@
 #define throw_out_of_range(msg)                                                \
   throw std::out_of_range(std::string(__FILE__) + ":" +                        \
                           std::to_string(__LINE__) + ":" + std::string(msg))
+
+#define time_marker() std::chrono::high_resolution_clock::now()
+
+#define diff_time(start, end)                                                  \
+  ({                                                                           \
+    std::chrono::duration<double> duration =                                   \
+        std::chrono::duration_cast<std::chrono::duration<double>>(end -        \
+                                                                  start);      \
+    duration.count();                                                          \
+  })
 #endif // BTCSUTILS_H_

test/simple.cpp

@@ -0,0 +1,36 @@
+#include "grid/BoundaryCondition.hpp"
+#include <Diffusion.hpp>
+#include <iostream>
+#include <vector>
+
+#define SIZE 5
+
+using namespace tug;
+
+int main() {
+  boundary_condition::BoundaryCondition bc(SIZE);
+
+  diffusion::TugDiffuInput in = {1., {SIZE, 0, 0}, {1, 0, 0}, 0};
+
+  std::vector<double> field1(SIZE, 0);
+  std::vector<double> field2(SIZE, 0);
+  std::vector<double> alpha(SIZE, 1e-5);
+
+  bc(boundary_condition::BC_SIDE_LEFT) = {boundary_condition::BC_TYPE_CONSTANT,
+                                          1};
+
+  for (int j = 0; j < 10; j++) {
+    double time1 =
+        diffusion::BTCS_Thomas_1D(field1.data(), alpha.data(), bc, in);
+    double time2 =
+        diffusion::BTCS_EigenLU_1D(field2.data(), alpha.data(), bc, in);
+
+    for (int i = 0; i < SIZE; i++) {
+      std::cout << field1[i] << "\t";
+      std::cout << field2[i] << "\n";
+    }
+
+    std::cout << "Time Thomas = " << time1 << "s \tTime EigenLU = " << time2
+              << " s" << std::endl;
+  }
+}

test/testBoundaryCondition.cpp

@@ -1,7 +1,7 @@
-#include <grid/BoundaryCondition.hpp>
+#include <BoundaryCondition.hpp>
 #include <doctest/doctest.h>
 
-using namespace tug::boundary_condition;
+using namespace tug::bc;
 
 #define BC_CONST_VALUE 1e-5
 
@@ -87,8 +87,8 @@ TEST_CASE("Boundary Condition helpers") {
   boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
 
   SUBCASE("return boundary condition skeleton") {
-    boundary_condition bc_test = BoundaryCondition::returnBoundaryCondition(
+    boundary_condition bc_test =
-        bc_set.type, bc_set.value);
+        BoundaryCondition::returnBoundaryCondition(bc_set.type, bc_set.value);
     CHECK_EQ(bc_test.value, bc_set.value);
     CHECK_EQ(bc_test.type, bc_set.type);
   }
@@ -98,9 +98,7 @@ TEST_CASE("1D special inner grid cells") {
   BoundaryCondition bc(5);
   boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
 
-  SUBCASE("valid set") {
+  SUBCASE("valid set") { CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set); }
-    CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set);
-  }
 
   SUBCASE("valid get") {
     bc(BC_INNER, 0) = bc_set;
@@ -112,9 +110,7 @@ TEST_CASE("1D special inner grid cells") {
     CHECK_THROWS(bc(BC_INNER, 5));
   }
 
-  SUBCASE("invalid set") {
+  SUBCASE("invalid set") { CHECK_THROWS(bc(BC_INNER, 5) = bc_set); }
-    CHECK_THROWS(bc(BC_INNER, 5) = bc_set);
-  }
 
   SUBCASE("valid row getter") {
     bc(BC_INNER, 1) = bc_set;
@@ -125,23 +121,16 @@ TEST_CASE("1D special inner grid cells") {
     CHECK_EQ(ret[1].type, bc_set.type);
   }
 
-  SUBCASE("invalid row getter") {
+  SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(1)); }
-    CHECK_THROWS(bc.getInnerRow(1));
-  }
-
-  SUBCASE("invalid col getter") {
-    CHECK_THROWS(bc.getInnerCol(0));
-  }
 
+  SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(0)); }
 }
 
 TEST_CASE("2D special inner grid cells") {
-  BoundaryCondition bc(5,5);
+  BoundaryCondition bc(5, 5);
   boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
 
-  SUBCASE("valid set") {
+  SUBCASE("valid set") { CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set); }
-    CHECK_NOTHROW(bc(BC_INNER, 0) = bc_set);
-  }
 
   SUBCASE("valid get") {
     bc(BC_INNER, 0) = bc_set;
@@ -153,9 +142,7 @@ TEST_CASE("2D special inner grid cells") {
     CHECK_THROWS(bc(BC_INNER, 25));
   }
 
-  SUBCASE("invalid set") {
+  SUBCASE("invalid set") { CHECK_THROWS(bc(BC_INNER, 25) = bc_set); }
-    CHECK_THROWS(bc(BC_INNER, 25) = bc_set);
-  }
 
   SUBCASE("valid row getter") {
     bc(BC_INNER, 0) = bc_set;
@@ -187,11 +174,7 @@ TEST_CASE("2D special inner grid cells") {
     CHECK_EQ(ret[1].type, BC_UNSET);
   }
 
-  SUBCASE("invalid row getter") {
+  SUBCASE("invalid row getter") { CHECK_THROWS(bc.getInnerRow(5)); }
-    CHECK_THROWS(bc.getInnerRow(5));
-  }
 
-  SUBCASE("invalid col getter") {
+  SUBCASE("invalid col getter") { CHECK_THROWS(bc.getInnerCol(5)); }
-    CHECK_THROWS(bc.getInnerCol(5));
-  }
 }

test/testDiffusion.cpp

@@ -1,10 +1,11 @@
+#include "Diffusion.hpp"
+#include "Solver.hpp"
+#include <BoundaryCondition.hpp>
 #include <bits/stdint-uintn.h>
-#include <diffusion/BTCSDiffusion.hpp>
 #include <doctest/doctest.h>
-#include <grid/BoundaryCondition.hpp>
 #include <vector>
 
-using namespace tug::boundary_condition;
+using namespace tug::bc;
 using namespace tug::diffusion;
 
 #define DIMENSION 2
@@ -14,13 +15,19 @@ using namespace tug::diffusion;
 
 static std::vector<double> alpha(N *M, 1e-3);
 
-static BTCSDiffusion setupDiffu(uint32_t n, uint32_t m) {
+static TugInput setupDiffu(BoundaryCondition &bc) {
-  BTCSDiffusion diffu(DIMENSION);
+  TugInput diffu;
 
-  diffu.setXDimensions(n, n);
+  diffu.time_step = 1.;
-  diffu.setYDimensions(m, m);
+  diffu.solver = tug::solver::ThomasAlgorithm;
 
-  diffu.setTimestep(1.);
+  diffu.grid.grid_cells[0] = N;
+  diffu.grid.grid_cells[1] = M;
+
+  diffu.grid.domain_size[0] = N;
+  diffu.grid.domain_size[1] = M;
+
+  diffu.grid.bc = &bc;
 
   return diffu;
 }
@@ -30,23 +37,22 @@ TEST_CASE("closed boundaries - 1 concentration to 1 - rest 0") {
 
   field[MID] = 1;
 
-  BTCSDiffusion diffu = setupDiffu(N, M);
   BoundaryCondition bc(N, M);
 
+  TugInput diffu = setupDiffu(bc);
+
   uint32_t iterations = 1000;
   double sum = 0;
 
   for (int t = 0; t < iterations; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc);
+    ADI_2D(diffu, field.data(), alpha.data());
+  }
 
-    if (t == iterations - 1) {
+  // iterate through rows
-      // iterate through rows
+  for (int i = 0; i < M; i++) {
-      for (int i = 0; i < M; i++) {
+    // iterate through columns
-        // iterate through columns
+    for (int j = 0; j < N; j++) {
-        for (int j = 0; j < N; j++) {
+      sum += field[i * N + j];
-          sum += field[i * N + j];
-        }
-      }
     }
   }
   CAPTURE(sum);
@@ -59,7 +65,6 @@ TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") {
 
   field[MID] = 1;
 
-  BTCSDiffusion diffu = setupDiffu(N, M);
   BoundaryCondition bc(N, M);
 
   boundary_condition input = {BC_TYPE_CONSTANT, 0};
@@ -69,13 +74,15 @@ TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") {
   bc.setSide(BC_SIDE_TOP, input);
   bc.setSide(BC_SIDE_BOTTOM, input);
 
+  TugInput diffu = setupDiffu(bc);
+
   uint32_t max_iterations = 20000;
   bool reached = false;
 
   int t = 0;
 
   for (t = 0; t < max_iterations; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc);
+    ADI_2D(diffu, field.data(), alpha.data());
 
     if (field[N * M - 1] > 1e-15) {
       reached = true;
@@ -95,7 +102,6 @@ TEST_CASE(
     "constant top and bottom (1 and 0) - left and right closed - 0 inlet") {
   std::vector<double> field(N * M, 0);
 
-  BTCSDiffusion diffu = setupDiffu(N, M);
   BoundaryCondition bc(N, M);
 
   boundary_condition top =
@@ -106,10 +112,12 @@ TEST_CASE(
   bc.setSide(BC_SIDE_TOP, top);
   bc.setSide(BC_SIDE_BOTTOM, bottom);
 
+  TugInput diffu = setupDiffu(bc);
+
   uint32_t max_iterations = 100;
 
   for (int t = 0; t < max_iterations; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc);
+    ADI_2D(diffu, field.data(), alpha.data());
   }
 
   for (int i = 0; i < N; i++) {
@@ -129,18 +137,19 @@ TEST_CASE("2D closed boundaries, 1 constant cell in the middle") {
   std::vector<double> field(N * M, 0);
   double val = 1e-2;
 
-  BTCSDiffusion diffu = setupDiffu(N, M);
   BoundaryCondition bc(N, M);
 
   field[MID] = val;
   bc(BC_INNER, MID) = {BC_TYPE_CONSTANT, val};
 
+  TugInput diffu = setupDiffu(bc);
+
   uint32_t max_iterations = 100;
 
   double sum = val;
 
   for (int t = 0; t < max_iterations; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc);
+    ADI_2D(diffu, field.data(), alpha.data());
 
     CHECK_EQ(field[MID], val);