diff --git a/CMakeLists.txt b/CMakeLists.txt index fa7f009a1..249953424 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.14) +cmake_minimum_required(VERSION 3.20) project(POET LANGUAGES CXX C diff --git a/README.md b/README.md index 4e01e1c59..f92d70c7f 100644 --- a/README.md +++ b/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 diff --git a/ext/tug b/ext/tug index 449647010..9c4aeee41 160000 --- a/ext/tug +++ b/ext/tug @@ -1 +1 @@ -Subproject commit 449647010ab9cdf9e405139f360424a2b21ab3ab +Subproject commit 9c4aeee410c71d064f7567143d4f8d6451ade75a diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ca39b6106..1727685f1 100644 --- a/src/CMakeLists.txt +++ b/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) diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp index d3ee24d5f..3e2fa9dee 100644 --- a/src/Init/GridInit.cpp +++ b/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 \ No newline at end of file +} // namespace poet diff --git a/src/Transport/DiffusionModule.cpp b/src/Transport/DiffusionModule.cpp index d7f1df5f3..85fc07aa4 100644 --- a/src/Transport/DiffusionModule.cpp +++ b/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 -** 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. -*/ + * DiffusionModule.cpp - Coupling module betweeen POET and tug + * Copyright (C) 2018-2025 Max Luebke (University of Potsdam) + */ #include "DiffusionModule.hpp" #include "Base/Macros.hpp" #include "Init/InitialList.hpp" -#include -#include -#include #include - -#include "tug/Boundary.hpp" -#include "tug/Grid.hpp" -#include "tug/Simulation.hpp" - #include -#include +#include +#include using namespace poet; -static inline std::vector -MatrixToVec(const Eigen::MatrixX &mat) { - std::vector 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 -VecToMatrix(const std::vector &vec, std::uint32_t n_rows, - std::uint32_t n_cols) { - Eigen::MatrixX 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 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 sim(grid, boundary); + tug::Diffusion diffusion_solver( + this->transport_field[sol_name].data(), n_rows, n_cols); #elif defined(POET_TUG_FTCS) - tug::Simulation sim(grid, boundary); + tug::Diffusion diffusion_solver( + this->transport_field[sol_name].data(), n_rows, n_cols); #else #error "No valid approach defined" #endif - sim.setIterations(1); + tug::RowMajMatMap alpha_x_map( + this->param_list.alpha_x[sol_name].data(), n_rows, n_cols); + tug::RowMajMatMap 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 conc = VecToMatrix(species_conc, n_rows, n_cols); - Eigen::MatrixX alpha_x = - VecToMatrix(this->param_list.alpha_x[sol_name], n_rows, n_cols); - Eigen::MatrixX 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 +50,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();