From 9c1afe8e2d46313d22f8085e0c2155253322d89e Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Tue, 1 Mar 2022 13:56:32 +0100
Subject: [PATCH] Implement 2D-row-wise in one direction

---
 src/BTCSDiffusion.cpp | 86 +++++++++++++++++++++++++++++++++----------
 src/BTCSDiffusion.hpp |  8 +++-
 2 files changed, 72 insertions(+), 22 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 2f6b477..a015eb2 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -7,6 +7,7 @@
 #include <Eigen/src/Core/Matrix.h>
 #include <algorithm>
 #include <cassert>
+#include <cstddef>
 #include <iomanip>
 #include <iterator>
 #include <ostream>
@@ -125,32 +126,28 @@ void Diffusion::BTCSDiffusion::simulate2D(
     Eigen::Map<DMatrixRowMajor> &c, Eigen::Map<const DMatrixRowMajor> &alpha,
     Eigen::Map<const BCMatrixRowMajor> &bc) {
 
-  // double local_dt = this->time_step / 2.;
-  // DMatrixRowMajor tmp_vector;
+  int n_rows, n_cols;
+
+  double local_dt = this->time_step / 2.;
+  double dx = this->deltas[0];
+  DMatrixRowMajor tmp_vector;
+
+  n_rows = this->grid_cells[1];
+  n_cols = this->grid_cells[0];
 
-  // int n_cols = c.cols();
   // unsigned int size = (this->grid_cells[0] + 2) * (this->grid_cells[1]);
 
-  // A_matrix.resize(size, size);
-  // A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
+  DMatrixRowMajor t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
 
-  // b_vector.resize(size);
-  // x_vector.resize(size);
+  std::cout << t0_c.row(0) << std::endl;
 
-  // for (int i = 0; i < c.rows(); i++) {
-  //   boundary_condition left = bc(i, 0);
-  //   bool left_constant = left.type == Diffusion::BC_CONSTANT;
-  //   boundary_condition right = bc(i, n_cols - 1);
-  //   bool right_constant = right.type == Diffusion::BC_CONSTANT;
+  for (int i = 0; i < n_rows; i++) {
+    DVectorRowMajor input_field = c.row(i);
+    simulate_base(input_field, bc.row(i), alpha.row(i), dx, local_dt, n_cols, t0_c.row(i));
+    c.row(i) << input_field;
+  }
 
-  //   fillMatrixFromRow(alpha.row(i), n_cols, i, left_constant, right_constant,
-  //                     deltas[0], this->time_step / 2, bc.row(i));
-  //   fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right,
-  //   local_dt,
-  //                        bc.row(i));
-  // }
-
-  // solveLES();
+  solveLES();
 
   // tmp_vector = x_vector;
   // tmp_vector.transposeInPlace();
@@ -195,6 +192,55 @@ void Diffusion::BTCSDiffusion::simulate2D(
   // c.transposeInPlace();
 }
 
+Diffusion::BTCSDiffusion::DMatrixRowMajor Diffusion::BTCSDiffusion::calc_t0_c(
+    const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
+    const BCMatrixRowMajor &bc, double time_step, double dx) {
+
+  int n_rows = this->grid_cells[1];
+  int n_cols = this->grid_cells[0];
+
+  DMatrixRowMajor t0_c(n_rows, n_cols);
+
+  double y_values[3];
+
+  // first, iterate over first row
+  for (int j = 0; j < n_cols; j++) {
+    y_values[0] = getBCFromFlux(bc(0, j), c(0, j), alpha(0, j));
+    y_values[1] = c(0, j);
+    y_values[2] = c(1, j);
+
+    t0_c(0, j) = time_step * alpha(0, j) *
+                 (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
+  }
+
+  // then iterate over inlet
+  for (int i = 1; i < n_rows - 1; i++) {
+    for (int j = 0; j < n_cols; j++) {
+
+      y_values[0] = c(i - 1, j);
+      y_values[1] = c(i, j);
+      y_values[2] = c(i + 1, j);
+
+    t0_c(i, j) = time_step * alpha(i, j) *
+                 (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
+    }
+  }
+
+  int end = n_rows - 1;
+
+  // and finally over last row
+  for (int j = 0; j < n_cols; j++) {
+    y_values[0] = c(end-1,j);
+    y_values[1] = c(end, j);
+    y_values[2] = getBCFromFlux(bc(end, j), c(end, j), alpha(end, j));
+
+    t0_c(end, j) = time_step * alpha(end, j) *
+                 (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
+  }
+
+  return t0_c;
+}
+
 inline void Diffusion::BTCSDiffusion::fillMatrixFromRow(
     const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size,
     double dx, double time_step) {
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 6657f97..1c9556a 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -116,8 +116,8 @@ private:
       BCVectorRowMajor;
 
   void simulate_base(DVectorRowMajor &c, const BCVectorRowMajor &bc,
-                     const DVectorRowMajor &alpha, double dx,
-                     double time_step, int size, const DVectorRowMajor &t0_c);
+                     const DVectorRowMajor &alpha, double dx, double time_step,
+                     int size, const DVectorRowMajor &t0_c);
 
   void simulate1D(Eigen::Map<DVectorRowMajor> &c,
                   Eigen::Map<const DVectorRowMajor> &alpha,
@@ -127,6 +127,10 @@ private:
                   Eigen::Map<const DMatrixRowMajor> &alpha,
                   Eigen::Map<const BCMatrixRowMajor> &bc);
 
+  DMatrixRowMajor calc_t0_c(const DMatrixRowMajor &c,
+                            const DMatrixRowMajor &alpha,
+                            const BCMatrixRowMajor &bc, double time_step, double dx);
+
   inline void fillMatrixFromRow(const DVectorRowMajor &alpha,
                                 const BCVectorRowMajor &bc, int size, double dx,
                                 double time_step);