Merge branch '19-cmake-requirements-missmatch-with-tug' into 'main'

Resolve "CMake requirements missmatch with TUG."

Closes #19

See merge request naaice/poet!57

CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.14)
+cmake_minimum_required(VERSION 3.20)
 
 project(POET
   LANGUAGES CXX C

README.md

@@ -29,7 +29,7 @@ To compile POET you need following software to be installed:
 
 - C/C++ compiler (tested with GCC)
 - MPI-Implementation (tested with OpenMPI and MVAPICH)
-- CMake 3.9+
+- CMake 3.20+
 - Eigen3 3.4+ (required by `tug`)
 - *optional*: `doxygen` with `dot` bindings for documentation
 - R language and environment including headers or `-dev` packages

ext/tug

@@ -1 +1 @@
-Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab
+Subproject commit 9c4aeee410c71d064f7567143d4f8d6451ade75a

src/CMakeLists.txt

@@ -38,7 +38,7 @@ FetchContent_Declare(
   cli11
   QUIET 
   GIT_REPOSITORY https://github.com/CLIUtils/CLI11.git
-  GIT_TAG v2.4.2
+  GIT_TAG v2.5.0
 )
 
 FetchContent_MakeAvailable(cli11)

src/Init/GridInit.cpp

@@ -1,3 +1,9 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+ * GridInit.cpp - Implementation of grid initialization
+ * Copyright (C) 2024-2025 Max Luebke (University of Potsdam)
+ */
+
 #include "InitialList.hpp"
 
 #include "PhreeqcMatrix.hpp"
@@ -147,7 +153,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 +201,4 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
   // }
 }
 
-} // namespace poet
+} // namespace poet

src/Transport/DiffusionModule.cpp

@@ -1,74 +1,21 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
 /*
-** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
+ * DiffusionModule.cpp - Coupling module betweeen POET and tug
-** Potsdam)
+ * Copyright (C) 2018-2025 Max Luebke (University of Potsdam)
-**
+ */
-** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
-**
-** Copyright (C) 2023-2024 Marco De Lucia (GFZ Potsdam), Max Luebke (University
-** of Potsdam)
-**
-** POET is free software; you can redistribute it and/or modify it under the
-** terms of the GNU General Public License as published by the Free Software
-** Foundation; either version 2 of the License, or (at your option) any later
-** version.
-**
-** POET is distributed in the hope that it will be useful, but WITHOUT ANY
-** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-**
-** You should have received a copy of the GNU General Public License along with
-** this program; if not, write to the Free Software Foundation, Inc., 51
-** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
 
 #include "DiffusionModule.hpp"
 
 #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 +23,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);
+  for (const auto &sol_name : this->param_list.transport_names) {
-  tug::Boundary boundary(grid);
-
-  grid.setDomain(param_list.s_rows, param_list.s_cols);
-
 #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) {
+    diffusion_solver.setAlphaX(alpha_x_map);
-    auto &species_conc = this->transport_field[sol_name];
+    diffusion_solver.setAlphaY(alpha_y_map);
 
-    Eigen::MatrixX<TugType> conc = VecToMatrix(species_conc, n_rows, n_cols);
+    auto &boundary = diffusion_solver.getBoundaryConditions();
-    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);
 
     boundary.deserialize(this->param_list.boundaries[sol_name]);
 
@@ -106,14 +50,8 @@ void DiffusionModule::simulate(double requested_dt) {
       boundary.setInnerBoundaries(this->param_list.inner_boundaries[sol_name]);
     }
 
-    grid.setAlpha(alpha_x, alpha_y);
+    diffusion_solver.setTimestep(requested_dt);
-    grid.setConcentrations(conc);
+    diffusion_solver.run();
-
-    sim.setTimestep(requested_dt);
-
-    sim.run();
-
-    species_conc = MatrixToVec(grid.getConcentrations());
   }
 
   const auto end_diffusion_t = std::chrono::high_resolution_clock::now();