From 00dad33b0f54a1d123043e2be78f07d50b72b4ec Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Max=20L=C3=BCbke?= <mluebke@gfz-potsdam.de>
Date: Fri, 5 Nov 2021 10:07:54 +0100
Subject: [PATCH] use new indexing method for linear system

---
 src/diffusion.cpp | 92 ++++++++++-------------------------------------
 1 file changed, 19 insertions(+), 73 deletions(-)

diff --git a/src/diffusion.cpp b/src/diffusion.cpp
index 7ecad7f..cda6f4b 100644
--- a/src/diffusion.cpp
+++ b/src/diffusion.cpp
@@ -1,6 +1,8 @@
 #include "diffusion.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>
@@ -9,8 +11,6 @@
 #include <Eigen/src/SparseCore/SparseMatrix.h>
 #include <Eigen/src/SparseCore/SparseMatrixBase.h>
 #include <Eigen/src/SparseLU/SparseLU.h>
-#include<Eigen/SparseCholesky>
-#include<Eigen/SparseQR>
 #include <Eigen/src/SparseQR/SparseQR.h>
 #include <iostream>
 
@@ -20,92 +20,38 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
   double dx = 1. / x;
   double dy = 1. / y;
 
-  int size = (x * y) - 4;
-
-  int local_x = x - 2;
+  int size = (x - 2) * (y - 2);
 
   Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
-  Eigen::VectorXd x_out(size);
+  Eigen::VectorXd x_out(x * y);
   std::vector<T> tripletList;
   tripletList.reserve(size * 5);
 
-  for (int i = x - 1 ; i < 2*x - 3 ; i++) {
-    double sx = (alpha[i + 2] * timestep) / (dx * dx);
-    double sy = (alpha[i + 2] * timestep) / (dy * dy);
+  int A_line = 0;
 
+  for (int i = 1; i < y - 1; i++) {
+    for (int j = 1; j < x - 1; j++) {
+      double sx = (alpha[i * x + j] * timestep) / (dx * dx);
+      double sy = (alpha[i * x + j] * timestep) / (dy * dy);
 
-    tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
-    tripletList.push_back(T(i, i - (x - 1), sy));
-    tripletList.push_back(T(i, i + x, sy));
-    tripletList.push_back(T(i, i + 1, sx));
-    tripletList.push_back(T(i, i - 1, sx));
+      tripletList.push_back(T(A_line, i * x + j, (1 + 2 * sx + 2 * sy)));
+      tripletList.push_back(T(A_line, (i - 1) * x + j, sy));
+      tripletList.push_back(T(A_line, (i + 1) * x + j, sy));
+      tripletList.push_back(T(A_line, i * x + (j + 1), sx));
+      tripletList.push_back(T(A_line, i * x + (j - 1), sx));
 
-    b[i] = -c[i+2];
+      b[A_line] = -c[i*x+j];
+      A_line++;
+    }
   }
 
-  for (int i = 2*x - 1; i < (y-2)*x-3; i++) {
-    double sx = (alpha[i + 2] * timestep) / (dx * dx);
-    double sy = (alpha[i + 2] * timestep) / (dy * dy);
-
-    tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
-    tripletList.push_back(T(i, i - x, sy));
-    tripletList.push_back(T(i, i + x, sy));
-    tripletList.push_back(T(i, i + 1, sx));
-    tripletList.push_back(T(i, i - 1, sx));
-
-    b[i] = -c[i+2];
-  }
-
-  for (int i = (y-2)*x-1; i < (y-1)*x-3; i++) {
-    double sx = (alpha[i + 2] * timestep) / (dx * dx);
-    double sy = (alpha[i + 2] * timestep) / (dy * dy);
-
-    tripletList.push_back(T(i, i, (1 + 2 * sx + 2 * sy)));
-    tripletList.push_back(T(i, i - x, sy));
-    tripletList.push_back(T(i, i + (x-1), sy));
-    tripletList.push_back(T(i, i + 1, sx));
-    tripletList.push_back(T(i, i - 1, sx));
-
-    b[i] = -c[i+2];
-  }
-
-  // for (int i = 0; i < (size-local_x); i++) {
-  //   int current = local_x + i;
-  //   double sx = (alpha[current] * timestep) / (dx * dx);
-  //   double sy = (alpha[current] * timestep) / (dy * dy);
-
-  //   tripletList.push_back(T(i, current, (1 + 2 * sx + 2 * sy)));
-  //   tripletList.push_back(T(i, current + local_x, sy));
-  //   tripletList.push_back(T(i, current - local_x, sy));
-  //   tripletList.push_back(T(i, current + 1, sx));
-  //   tripletList.push_back(T(i, current - 1, sx));
-
-  //   std::cout << current << std::endl;
-
-  //   b[i] = -c[x+i];
-  // }
-
-  // for (int i = 0; i < y; i++) {
-  //   for (int j = 0; j < x; j++) {
-  //     double sx = (alpha[i * y + j] * timestep) / (dx * dx);
-  //     double sy = (alpha[i * y + j] * timestep) / (dy * dy);
-
-  //     tripletList.push_back(T((i * x) + j, (i * x) + j, (1 + 2 * sx + 2 *
-  //     sy))); tripletList.push_back(T((i * x) + j, ((i * x) + j) + x, sy));
-  //     tripletList.push_back(T((i * x) + j, ((i * x) + j) - x, sy));
-  //     tripletList.push_back(T((i * x) + j, ((i * x) + j) + 1, sx));
-  //     tripletList.push_back(T((i * x) + j, ((i * x) + j) - 1, sx));
-
-  //     b[(i * x) + j] = -c[(i * x) + j];
-  //   }
-  // }
-
   std::cout << b << std::endl;
 
   Eigen::SparseMatrix<double> A(size, (x * y) - 4);
   A.setFromTriplets(tripletList.begin(), tripletList.end());
 
-Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
+  Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
+      solver;
 
   // Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
   //     solver;