Merge branch '2D' of git.gfz-potsdam.de:mluebke/diffusion into 2D

.gitlab-ci.yml

@@ -0,0 +1,22 @@
+image: sobc/gitlab-ci
+
+stages:
+    - build
+    - test
+
+before_script:
+    - apt-get update && apt-get install -y libeigen3-dev
+
+build:
+    stage: build
+    script:
+        - mkdir build && cd build
+        - cmake ..
+        - make
+
+lint:
+    stage: test
+    script:
+        - mkdir build && cd build
+        - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*,readability-*, modernize-*" ..
+        - make

doc/ADI_scheme.org

@@ -96,7 +96,7 @@ c(i,N-1) & \text{else}
 
 *** Inlet
 
-$p(i,j) = \alpha(i,j)\frac{c(i-1,j) - 2\cdot c(i,j) + c(i+1,j)}{\Delta x^2}$[fn:1]
+$p(i,j) = \frac{\Delta t}{2}\alpha(i,j)\frac{c(i-1,j) - 2\cdot c(i,j) + c(i+1,j)}{\Delta x^2}$[fn:1]
 
 $b(i,j) = \begin{cases}
 bc(i,j).\text{value} & \text{if } bc(i,N-1) = \text{constant} \\

src/BTCSDiffusion.cpp

@@ -139,6 +139,7 @@ void BTCSDiffusion::simulate1D(Eigen::Map<DVectorRowMajor> &c,
 void BTCSDiffusion::simulate2D(Eigen::Map<DMatrixRowMajor> &c,
                                Eigen::Map<const DMatrixRowMajor> &alpha) {
 
+  double local_dt = this->time_step/ 2.;
   DMatrixRowMajor tmp_vector;
 
   int n_cols = c.cols();
@@ -158,7 +159,7 @@ void BTCSDiffusion::simulate2D(Eigen::Map<DMatrixRowMajor> &c,
 
     fillMatrixFromRow(alpha.row(i), n_cols, i, left_constant, right_constant,
                       deltas[0], this->time_step / 2);
-    fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right);
+    fillVectorFromRowADI(c, alpha.row(i), i, deltas[0], left, right, local_dt);
   }
 
   solveLES();
@@ -190,7 +191,7 @@ void BTCSDiffusion::simulate2D(Eigen::Map<DMatrixRowMajor> &c,
 
     fillMatrixFromRow(alpha.col(i), n_cols, i, left_constant, right_constant,
                       deltas[1], this->time_step / 2);
-    fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right);
+    fillVectorFromRowADI(c, alpha.row(i), i, deltas[1], left, right, local_dt);
   }
 
   solveLES();
@@ -238,9 +239,10 @@ void BTCSDiffusion::fillMatrixFromRow(const DVectorRowMajor &alpha, int n_cols,
 }
 
 void BTCSDiffusion::fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
-                                        const Eigen::VectorXd alpha, int row,
+                                         const Eigen::VectorXd alpha, int row,
-                                        double delta, boundary_condition left,
+                                         double delta, boundary_condition left,
-                                        boundary_condition right) {
+                                         boundary_condition right,
+                                         double time_step) {
 
   int ncol = c.cols();
   int nrow = c.rows();
@@ -257,27 +259,26 @@ void BTCSDiffusion::fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
         getBCFromFlux(right, c(row, ncol - 1), alpha[ncol - 1]);
   }
 
-  for (int j = 1; j < offset - 1; j++) {
+  for (int j = 0; j < ncol; j++) {
-    boundary_condition tmp_bc = this->bc(row, j-1);
+    boundary_condition tmp_bc = this->bc(row, j);
 
     if (tmp_bc.type == BTCSDiffusion::BC_CONSTANT) {
-      b_vector[offset * row + j] = tmp_bc.value;
+      b_vector[offset * row + (j+1)] = tmp_bc.value;
       continue;
     }
 
     double y_values[3];
     y_values[0] =
-        (row != 0 ? c(row - 1, j - 1)
+        (row != 0 ? c(row - 1, j) : getBCFromFlux(tmp_bc, c(row, j), alpha[j]));
-                  : getBCFromFlux(tmp_bc, c(row, j - 1), alpha[j - 1]));
+    y_values[1] = c(row, j);
-    y_values[1] = c(row, j - 1);
     y_values[2] =
-        (row != nrow - 1 ? c(row + 1, j - 1)
+        (row != nrow - 1 ? c(row + 1, j)
-                         : getBCFromFlux(tmp_bc, c(row, j - 1), alpha[j - 1]));
+                         : getBCFromFlux(tmp_bc, c(row, j), alpha[j]));
 
     double t0_c =
-        alpha[j - 1] *
+        time_step * alpha[j] *
         ((y_values[0] - 2 * y_values[1] + y_values[2]) / (delta * delta));
-    b_vector[offset * row + j] = -c(row, j - 1) - t0_c;
+    b_vector[offset * row + (j + 1)] = -c(row, j) - (t0_c);
   }
 }
 

src/BTCSDiffusion.hpp

@@ -153,8 +153,9 @@ private:
                          bool left_constant, bool right_constant, double delta,
                          double time_step);
   void fillVectorFromRowADI(Eigen::Map<DMatrixRowMajor> &c,
-                           const Eigen::VectorXd alpha, int row, double delta,
+                            const Eigen::VectorXd alpha, int row, double delta,
-                           boundary_condition left, boundary_condition right);
+                            boundary_condition left, boundary_condition right,
+                            double time_step);
   void simulate3D(std::vector<double> &c);
   inline double getBCFromFlux(boundary_condition bc, double nearest_value,
                               double neighbor_alpha);