From c8bedaac28422a4ddf2478e3aa0e425caa891f59 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Max=20L=C3=BCbke?= <mluebke@uni-potsdam.de>
Date: Fri, 24 Oct 2025 09:41:03 +0200
Subject: [PATCH] refactor(diffusion): integrate tug::Diffusion solver

---
 src/Init/GridInit.cpp             |  4 +-
 src/Transport/DiffusionModule.cpp | 77 +++++++------------------------
 2 files changed, 18 insertions(+), 63 deletions(-)

diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp
index d3ee24d5f..75adc7061 100644
--- a/src/Init/GridInit.cpp
+++ b/src/Init/GridInit.cpp
@@ -147,7 +147,7 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
   PhreeqcMatrix pqc_mat =
       PhreeqcMatrix(database, script, with_h0_o0, with_redox);
 
-  this->transport_names = pqc_mat.getSolutionNames();
+  this->transport_names = pqc_mat.getMatrixTransported();
 
   Rcpp::Function unique_R("unique");
   Rcpp::Function sort_R("sort");
@@ -195,4 +195,4 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
   // }
 }
 
-} // namespace poet
\ No newline at end of file
+} // namespace poet
diff --git a/src/Transport/DiffusionModule.cpp b/src/Transport/DiffusionModule.cpp
index d7f1df5f3..0538662b9 100644
--- a/src/Transport/DiffusionModule.cpp
+++ b/src/Transport/DiffusionModule.cpp
@@ -26,49 +26,13 @@
 #include "Base/Macros.hpp"
 #include "Init/InitialList.hpp"
 
-#include <Eigen/src/Core/Map.h>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
 #include <chrono>
-
-#include "tug/Boundary.hpp"
-#include "tug/Grid.hpp"
-#include "tug/Simulation.hpp"
-
 #include <string>
-#include <vector>
+#include <tug/Boundary.hpp>
+#include <tug/Diffusion.hpp>
 
 using namespace poet;
 
-static inline std::vector<TugType>
-MatrixToVec(const Eigen::MatrixX<TugType> &mat) {
-  std::vector<TugType> vec(mat.rows() * mat.cols());
-
-  for (std::uint32_t i = 0; i < mat.cols(); i++) {
-    for (std::uint32_t j = 0; j < mat.rows(); j++) {
-      vec[j * mat.cols() + i] = mat(j, i);
-    }
-  }
-
-  return vec;
-}
-
-static inline Eigen::MatrixX<TugType>
-VecToMatrix(const std::vector<TugType> &vec, std::uint32_t n_rows,
-            std::uint32_t n_cols) {
-  Eigen::MatrixX<TugType> mat(n_rows, n_cols);
-
-  for (std::uint32_t i = 0; i < n_cols; i++) {
-    for (std::uint32_t j = 0; j < n_rows; j++) {
-      mat(j, i) = vec[j * n_cols + i];
-    }
-  }
-
-  return mat;
-}
-
-// static constexpr double ZERO_MULTIPLICATOR = 10E-14;
-
 void DiffusionModule::simulate(double requested_dt) {
   // MSG("Starting diffusion ...");
   const auto start_diffusion_t = std::chrono::high_resolution_clock::now();
@@ -76,29 +40,26 @@ void DiffusionModule::simulate(double requested_dt) {
   const auto &n_rows = this->param_list.n_rows;
   const auto &n_cols = this->param_list.n_cols;
 
-  tug::Grid<TugType> grid(param_list.n_rows, param_list.n_cols);
-  tug::Boundary boundary(grid);
-
-  grid.setDomain(param_list.s_rows, param_list.s_cols);
-
+  for (const auto &sol_name : this->param_list.transport_names) {
 #if defined(POET_TUG_BTCS)
-  tug::Simulation<TugType> sim(grid, boundary);
+    tug::Diffusion<TugType> diffusion_solver(
+        this->transport_field[sol_name].data(), n_rows, n_cols);
 #elif defined(POET_TUG_FTCS)
-  tug::Simulation<TugType, tug::FTCS_APPROACH> sim(grid, boundary);
+    tug::Diffusion<TugType, tug::FTCS_APPROACH> diffusion_solver(
+        this->transport_field[sol_name].data(), n_rows, n_cols);
 #else
 #error "No valid approach defined"
 #endif
 
-  sim.setIterations(1);
+    tug::RowMajMatMap<TugType> alpha_x_map(
+        this->param_list.alpha_x[sol_name].data(), n_rows, n_cols);
+    tug::RowMajMatMap<TugType> alpha_y_map(
+        this->param_list.alpha_y[sol_name].data(), n_rows, n_cols);
 
-  for (const auto &sol_name : this->param_list.transport_names) {
-    auto &species_conc = this->transport_field[sol_name];
+    diffusion_solver.setAlphaX(alpha_x_map);
+    diffusion_solver.setAlphaY(alpha_y_map);
 
-    Eigen::MatrixX<TugType> conc = VecToMatrix(species_conc, n_rows, n_cols);
-    Eigen::MatrixX<TugType> alpha_x =
-        VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols);
-    Eigen::MatrixX<TugType> alpha_y =
-        VecToMatrix(this->param_list.alpha_y[sol_name], n_rows, n_cols);
+    auto &boundary = diffusion_solver.getBoundaryConditions();
 
     boundary.deserialize(this->param_list.boundaries[sol_name]);
 
@@ -106,14 +67,8 @@ void DiffusionModule::simulate(double requested_dt) {
       boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]);
     }
 
-    grid.setAlpha(alpha_x, alpha_y);
-    grid.setConcentrations(conc);
-
-    sim.setTimestep(requested_dt);
-
-    sim.run();
-
-    species_conc = MatrixToVec(grid.getConcentrations());
+    diffusion_solver.setTimestep(requested_dt);
+    diffusion_solver.run();
   }
 
   const auto end_diffusion_t = std::chrono::high_resolution_clock::now();