From 3380eb4a8ff1c53dc955ca5c4139ab3348feff4a Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 13 Mar 2024 16:57:44 +0100 Subject: [PATCH] Add iphreeqc and init_r_lib submodules, and make necessary changes to CMakeLists.txt and src/initializer.cpp --- .gitmodules | 9 ++-- CMakeLists.txt | 9 ++-- R_lib/init_r_lib.R | 37 ++++++++++++++ ext/iphreeqc | 1 + ext/phreeqcrm | 1 - src/CMakeLists.txt | 32 ++++++++++-- src/Init/GridInit.cpp | 108 +++++++++++++++++++++++++++++++++++++++ src/Init/InitialList.cpp | 7 +++ src/Init/InitialList.hpp | 44 ++++++++++++++++ src/initializer.cpp | 35 +++++++++++++ src/poet.hpp.in | 2 + 11 files changed, 273 insertions(+), 12 deletions(-) create mode 100644 R_lib/init_r_lib.R create mode 160000 ext/iphreeqc delete mode 160000 ext/phreeqcrm create mode 100644 src/Init/GridInit.cpp create mode 100644 src/Init/InitialList.cpp create mode 100644 src/Init/InitialList.hpp create mode 100644 src/initializer.cpp diff --git a/.gitmodules b/.gitmodules index 744dfdb41..645e48d00 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,9 +2,12 @@ path = ext/tug url = ../tug.git -[submodule "ext/phreeqcrm"] - path = ext/phreeqcrm - url = ../phreeqcrm-gfz.git [submodule "ext/doctest"] path = ext/doctest url = https://github.com/doctest/doctest.git +[submodule "ext/phreeqc"] + path = ext/phreeqc + url = git@git.gfz-potsdam.de:naaice/phreeqcpoet.git +[submodule "ext/iphreeqc"] + path = ext/iphreeqc + url = git@git.gfz-potsdam.de:naaice/iphreeqc.git diff --git a/CMakeLists.txt b/CMakeLists.txt index e6701b564..cef0776ed 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,6 +22,9 @@ find_package(MPI REQUIRED) find_package(RRuntime REQUIRED) +add_compile_options(-fsanitize=address -fno-omit-frame-pointer) +add_link_options(-fsanitize=address) + add_subdirectory(src) add_subdirectory(bench) @@ -29,7 +32,7 @@ add_subdirectory(bench) set(TUG_ENABLE_TESTING OFF CACHE BOOL "" FORCE) add_subdirectory(ext/tug EXCLUDE_FROM_ALL) -add_subdirectory(ext/phreeqcrm EXCLUDE_FROM_ALL) +add_subdirectory(ext/iphreeqc EXCLUDE_FROM_ALL) option(POET_ENABLE_TESTING "Build test suite for POET" OFF) @@ -40,6 +43,4 @@ endif() option(BUILD_DOC "Build documentation with doxygen" OFF) -add_subdirectory(docs) - -add_subdirectory(apps) \ No newline at end of file +add_subdirectory(docs) \ No newline at end of file diff --git a/R_lib/init_r_lib.R b/R_lib/init_r_lib.R new file mode 100644 index 000000000..f2dfbb71c --- /dev/null +++ b/R_lib/init_r_lib.R @@ -0,0 +1,37 @@ +input <- readRDS("/home/max/Storage/poet/build/apps/out.rds") + +grid_def <- matrix(c(2, 3), nrow = 2, ncol = 5) + +library(data.table) + +pqc_to_grid <- function(pqc_in, grid) { + # Convert the input DataFrame to a data.table + dt <- data.table(pqc_in) + + # Flatten the matrix into a vector + id_vector <- as.vector(t(grid)) + + # Initialize an empty data.table to store the results + result_dt <- data.table() + + # Iterate over each ID in the vector + for (id_mat in id_vector) { + # Find the matching row in the data.table + matching_dt <- dt[id == id_mat] + + # Append the matching data.table row to the result data.table + result_dt <- rbind(result_dt, matching_dt) + } + + # Remove all columns which only contain NaN + # result_dt <- result_dt[, colSums(is.na(result_dt)) != nrow(result_dt)] + + res_df <- as.data.frame(result_dt) + + return(res_df[, colSums(is.na(res_df)) != nrow(res_df)]) +} + +pqc_init <- pqc_to_grid(input, grid_def) +test <- pqc_init + +# remove column with all NA diff --git a/ext/iphreeqc b/ext/iphreeqc new file mode 160000 index 000000000..f3f86bb4a --- /dev/null +++ b/ext/iphreeqc @@ -0,0 +1 @@ +Subproject commit f3f86bb4ac1d3b70aec7437781875072f2896ae3 diff --git a/ext/phreeqcrm b/ext/phreeqcrm deleted file mode 160000 index 6ed14c353..000000000 --- a/ext/phreeqcrm +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 6ed14c35322a245e3a9776ef262c0ac0eba3b301 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1494c0849..3f6d4f05c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,3 +1,20 @@ +add_library(poetinitlib + Init/InitialList.cpp + Init/GridInit.cpp +) + +target_link_libraries(poetinitlib PUBLIC + RRuntime + IPhreeqcPOET + tug +) + +target_include_directories(poetinitlib PUBLIC + "${CMAKE_CURRENT_SOURCE_DIR}/Init" +) + +target_compile_definitions(poetinitlib PUBLIC STRICT_R_HEADERS) + add_library(poetlib Base/Grid.cpp Base/SimParams.cpp @@ -17,8 +34,9 @@ target_link_libraries(poetlib PUBLIC MPI::MPI_C ${MATH_LIBRARY} RRuntime - PhreeqcRM + PhreeqcRM tug + poetinitlib ) target_compile_definitions(poetlib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX) @@ -39,11 +57,17 @@ if (POET_PHT_ADDITIONAL_INFO) endif() file(READ "${PROJECT_SOURCE_DIR}/R_lib/kin_r_library.R" R_KIN_LIB ) +file(READ "${PROJECT_SOURCE_DIR}/R_lib/init_r_lib.R" R_INIT_LIB) configure_file(poet.hpp.in poet.hpp @ONLY) -add_executable(poet poet.cpp) -target_link_libraries(poet PRIVATE poetlib MPI::MPI_C RRuntime) -target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") +# add_executable(poet poet.cpp) +# target_link_libraries(poet PRIVATE poetlib MPI::MPI_C RRuntime) +# target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") + +add_executable(poet_init initializer.cpp) +target_link_libraries(poet_init PRIVATE poetinitlib RRuntime) +target_include_directories(poet_init PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") install(TARGETS poet DESTINATION bin) + diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp new file mode 100644 index 000000000..c3cc0b379 --- /dev/null +++ b/src/Init/GridInit.cpp @@ -0,0 +1,108 @@ +#include "InitialList.hpp" + +#include +#include +#include + +namespace poet { +static Rcpp::NumericMatrix +pqcScriptToGrid(RInside &R, const std::string &script, const std::string &db, + std::map &raw_dumps) { + IPhreeqcPOET phreeqc(db, script); + + const auto solution_ids = phreeqc.getSolutionIds(); + auto col_names = phreeqc.getInitNames(); + const auto concentrations = phreeqc.getInitValues(); + + // add "id" to the front of the column names + + const std::size_t ncols = col_names.size(); + const std::size_t nrows = solution_ids.size(); + + col_names.insert(col_names.begin(), "id"); + + Rcpp::NumericMatrix mat(nrows, ncols + 1); + + for (std::size_t i = 0; i < nrows; i++) { + mat(i, 0) = solution_ids[i]; + for (std::size_t j = 0; j < ncols; ++j) { + mat(i, j + 1) = concentrations[i][j]; + } + } + + Rcpp::colnames(mat) = Rcpp::wrap(col_names); + + raw_dumps = phreeqc.raw_dumps(); + + return mat; +} + +static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat, + const Rcpp::NumericMatrix &grid) { + // Rcpp::List res; + // res["init_mat"] = mat; + // res["grid_mat"] = grid; + + Rcpp::Function pqc_to_grid("pqc_to_grid"); + + // R.parseEvalQ("res_df <- pqc_to_grid(init_mat, grid_mat)"); + + return pqc_to_grid(mat, grid); +} + +void InitialList::initGrid(const Rcpp::List &grid_input) { + // parse input values + const std::string script = Rcpp::as(grid_input["pqc_in"]); + const std::string database = Rcpp::as(grid_input["pqc_db"]); + + std::string script_rp(PATH_MAX, '\0'); + std::string database_rp(PATH_MAX, '\0'); + + if (realpath(script.c_str(), script_rp.data()) == nullptr) { + throw std::runtime_error("Failed to resolve the path of the Phreeqc input" + + script); + } + + if (realpath(database.c_str(), database_rp.data()) == nullptr) { + throw std::runtime_error("Failed to resolve the path of the database" + + database); + } + + Rcpp::NumericMatrix grid_def = grid_input["grid_def"]; + Rcpp::NumericVector grid_size = grid_input["grid_size"]; + // Rcpp::NumericVector constant_cells = grid["constant_cells"].; + + this->n_rows = grid_def.nrow(); + this->n_cols = grid_def.ncol(); + + this->s_rows = grid_size[0]; + this->s_cols = grid_size[1]; + + this->dim = n_cols == 1 ? 1 : 2; + + if (this->n_cols > 1 && this->n_rows == 1) { + throw std::runtime_error( + "Dimensions of grid definition does not match the expected format " + "(n_rows > 1, n_cols >= 1)."); + } + + if (this->dim != grid_size.size()) { + throw std::runtime_error( + "Dimensions of grid does not match the dimension of grid definition."); + } + + if (this->s_rows <= 0 || this->s_cols <= 0) { + throw std::runtime_error("Grid size must be positive."); + } + + this->phreeqc_mat = + pqcScriptToGrid(R, script_rp, database_rp, this->pqc_raw_dumps); + this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def); + + // R["pqc_mat"] = this->phreeqc_mat; + // R["grid_def"] = initial_grid; + + // R.parseEval("print(pqc_mat)"); + // R.parseEval("print(grid_def)"); +} +} // namespace poet \ No newline at end of file diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp new file mode 100644 index 000000000..d12c6e314 --- /dev/null +++ b/src/Init/InitialList.cpp @@ -0,0 +1,7 @@ +#include "InitialList.hpp" + +namespace poet { +void InitialList::initializeFromList(const Rcpp::List &setup) { + initGrid(setup["grid"]); +} +} // namespace poet \ No newline at end of file diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp new file mode 100644 index 000000000..e3f86ebf6 --- /dev/null +++ b/src/Init/InitialList.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include "../DataStructures/Field.hpp" + +#include +#include +#include +#include + +namespace poet { + +class InitialList { +public: + InitialList(RInside &R) : R(R){}; + + void initializeFromList(const Rcpp::List &setup); + + void importList(const std::string &file_name); + Rcpp::List exportList(const std::string &file_name); + +private: + RInside &R; + + // Grid members + void initGrid(const Rcpp::List &grid_input); + std::uint8_t dim; + + std::uint32_t n_cols; + std::uint32_t n_rows; + + double s_cols; + double s_rows; + + std::vector constant_cells; + + Rcpp::List initial_grid; + + // No export + Rcpp::NumericMatrix phreeqc_mat; + + // Initialized by grid, modified by chemistry + std::map pqc_raw_dumps; +}; +} // namespace poet \ No newline at end of file diff --git a/src/initializer.cpp b/src/initializer.cpp new file mode 100644 index 000000000..8d1c0221e --- /dev/null +++ b/src/initializer.cpp @@ -0,0 +1,35 @@ +#include "Init/InitialList.hpp" +#include "poet.hpp" + +#include +#include + +#include +#include + +int main(int argc, char **argv) { + if (argc < 2 || argc > 2) { + std::cerr << "Usage: " << argv[0] << " \n"; + return 1; + } + + RInside R(argc, argv); + + R.parseEvalQ(init_r_library); + + const std::string script = argv[1]; + R.parseEvalQ("source('" + script + "')"); + + Rcpp::List setup = R["setup"]; + + // Rcpp::List grid = R.parseEval("setup$grid"); + + // Rcpp::List results; + + poet::InitialList init(R); + + init.initializeFromList(setup); + + // parseGrid(R, grid, results); + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/src/poet.hpp.in b/src/poet.hpp.in index 666fbaa08..115ac7dd7 100644 --- a/src/poet.hpp.in +++ b/src/poet.hpp.in @@ -7,4 +7,6 @@ static const char *poet_version = "@POET_VERSION@"; // using the Raw string literal to avoid escaping the quotes static inline std::string kin_r_library = R"(@R_KIN_LIB@)"; +static inline std::string init_r_library = R"(@R_INIT_LIB@)"; + #endif // POET_H