Add iphreeqc and init_r_lib submodules, and make necessary changes to CMakeLists.txt and src/initializer.cpp

.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

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(docs)
-
-add_subdirectory(apps)

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

ext/iphreeqc

@@ -0,0 +1 @@
+Subproject commit f3f86bb4ac1d3b70aec7437781875072f2896ae3

ext/phreeqcrm

@@ -1 +0,0 @@
-Subproject commit 6ed14c35322a245e3a9776ef262c0ac0eba3b301

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)
+# add_executable(poet poet.cpp)
-target_link_libraries(poet PRIVATE poetlib MPI::MPI_C RRuntime)
+# target_link_libraries(poet PRIVATE poetlib MPI::MPI_C RRuntime)
-target_include_directories(poet PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
+# 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)
+

src/Init/GridInit.cpp

@@ -0,0 +1,108 @@
+#include "InitialList.hpp"
+
+#include <IPhreeqcPOET.hpp>
+#include <Rcpp/Function.h>
+#include <Rcpp/vector/instantiation.h>
+
+namespace poet {
+static Rcpp::NumericMatrix
+pqcScriptToGrid(RInside &R, const std::string &script, const std::string &db,
+                std::map<int, std::string> &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<std::string>(grid_input["pqc_in"]);
+  const std::string database = Rcpp::as<std::string>(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

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

src/Init/InitialList.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include "../DataStructures/Field.hpp"
+
+#include <RInside.h>
+#include <Rcpp.h>
+#include <Rcpp/vector/instantiation.h>
+#include <cstdint>
+
+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<std::uint32_t> constant_cells;
+
+  Rcpp::List initial_grid;
+
+  // No export
+  Rcpp::NumericMatrix phreeqc_mat;
+
+  // Initialized by grid, modified by chemistry
+  std::map<int, std::string> pqc_raw_dumps;
+};
+} // namespace poet

src/initializer.cpp

@@ -0,0 +1,35 @@
+#include "Init/InitialList.hpp"
+#include "poet.hpp"
+
+#include <Rcpp/vector/instantiation.h>
+#include <cstdlib>
+
+#include <RInside.h>
+#include <Rcpp.h>
+
+int main(int argc, char **argv) {
+  if (argc < 2 || argc > 2) {
+    std::cerr << "Usage: " << argv[0] << " <script.R>\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;
+}

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