Add support for missing species in init grid, which are injected (by boundary condition)

R_lib/init_r_lib.R

@@ -29,7 +29,7 @@ pqc_to_grid <- function(pqc_in, grid) {
     return(res_df)
 }
 
-resolvePqcBound <- function(pqc_mat, transport_spec, id) {
+resolve_pqc_bound <- function(pqc_mat, transport_spec, id) {
     df <- as.data.frame(pqc_mat, check.names = FALSE)
     value <- df[df$ID == id, transport_spec]
 
@@ -39,3 +39,29 @@ resolvePqcBound <- function(pqc_mat, transport_spec, id) {
 
     return(value)
 }
+
+add_column_after_position <- function(df, new_col, pos, new_col_name) {
+    # Split the data frame into two parts
+    df_left <- df[, 1:(pos)]
+    df_right <- df[, (pos + 1):ncol(df)]
+
+    # Add the new column to the left part
+    df_left[[new_col_name]] <- new_col
+
+    # Combine the left part, new column, and right part
+    df_new <- cbind(df_left, df_right)
+
+    return(df_new)
+}
+
+add_missing_transport_species <- function(init_grid, new_names, old_size) {
+    # skip the ID column
+    column_index <- old_size + 1
+
+    for (name in new_names) {
+        init_grid <- add_column_after_position(init_grid, rep(0, nrow(init_grid)), column_index, name)
+        column_index <- column_index + 1
+    }
+
+    return(init_grid)
+}

ext/iphreeqc

@@ -1 +1 @@
-Subproject commit c0c830eab7c49ec0bf4ef87ef0669f6c15bae58c
+Subproject commit 4ec8b7006c215ad9fc1310505cd56d03ba4b17dd

src/Init/DiffusionInit.cpp

@@ -6,13 +6,16 @@
 #include <set>
 
 #include "DataStructures/Field.hpp"
+#include "IPhreeqcPOET.hpp"
 #include "InitialList.hpp"
 
+#include <Rcpp/Function.h>
 #include <Rcpp/proxy/ProtectedProxy.h>
 #include <Rcpp/vector/instantiation.h>
 #include <Rinternals.h>
 #include <cstddef>
 #include <cstdint>
+#include <memory>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -49,9 +52,81 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
   }
 }
 
+static std::vector<std::string> extend_transport_names(
+    std::unique_ptr<IPhreeqcPOET> &phreeqc, const Rcpp::List &boundaries_list,
+    const std::vector<std::string> &old_trans_names, Rcpp::List &initial_grid) {
+
+  std::vector<std::string> transport_names = old_trans_names;
+  std::set<int> constant_pqc_ids;
+
+  for (const auto &side : tug_side_mapping) {
+    if (!boundaries_list.containsElementNamed(side.second.c_str())) {
+      continue;
+    }
+
+    Rcpp::List mapping = boundaries_list[side.second];
+
+    const Rcpp::NumericVector cells = mapping["cell"];
+    const Rcpp::NumericVector values = mapping["sol_id"];
+    const Rcpp::CharacterVector type_str = mapping["type"];
+
+    if (cells.size() != values.size()) {
+      throw std::runtime_error("Boundary [" + side.second +
+                               "] cells and values are not the same "
+                               "length");
+    }
+
+    for (std::size_t i = 0; i < cells.size(); i++) {
+      if (type_str[i] == "constant") {
+        constant_pqc_ids.insert(values[i]);
+      }
+    }
+  }
+
+  if (!constant_pqc_ids.empty()) {
+    for (const auto &pqc_id : constant_pqc_ids) {
+      const auto solution_names = phreeqc->getSolutionNames(pqc_id);
+
+      // add those strings which are not already in the transport_names
+      for (const auto &name : solution_names) {
+        if (std::find(transport_names.begin(), transport_names.end(), name) ==
+            transport_names.end()) {
+          transport_names.push_back(name);
+        }
+      }
+    }
+  }
+
+  return transport_names;
+}
+
+static Rcpp::List extend_initial_grid(const Rcpp::List &initial_grid,
+                                      std::vector<std::string> transport_names,
+                                      std::size_t old_size) {
+  std::vector<std::string> names_to_add(transport_names.begin() + old_size,
+                                        transport_names.end());
+
+  Rcpp::Function extend_grid_R("add_missing_transport_species");
+
+  return extend_grid_R(initial_grid, Rcpp::wrap(names_to_add), old_size);
+}
+
 Rcpp::List InitialList::resolveBoundaries(const Rcpp::List &boundaries_list) {
   Rcpp::List bound_list;
-  Rcpp::Function resolve_R("resolvePqcBound");
+  Rcpp::Function resolve_R("resolve_pqc_bound");
+
+  const std::size_t old_transport_size = this->transport_names.size();
+
+  this->transport_names =
+      extend_transport_names(this->phreeqc, boundaries_list,
+                             this->transport_names, this->initial_grid);
+
+  const std::size_t new_transport_size = this->transport_names.size();
+
+  if (old_transport_size != new_transport_size) {
+    this->initial_grid = extend_initial_grid(
+        this->initial_grid, this->transport_names, old_transport_size);
+  }
 
   for (const auto &species : this->transport_names) {
     Rcpp::List spec_list;

src/Init/GridInit.cpp

@@ -6,6 +6,7 @@
 #include <Rcpp/vector/instantiation.h>
 #include <cstdint>
 #include <map>
+#include <memory>
 #include <regex>
 #include <sstream>
 #include <string>
@@ -13,8 +14,9 @@
 
 namespace poet {
 
-static Rcpp::NumericMatrix pqcScriptToGrid(IPhreeqcPOET &phreeqc, RInside &R) {
-  IPhreeqcPOET::PhreeqcMat phreeqc_mat = phreeqc.getPhreeqcMat();
+static Rcpp::NumericMatrix
+pqcScriptToGrid(std::unique_ptr<IPhreeqcPOET> &phreeqc, RInside &R) {
+  IPhreeqcPOET::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat();
 
   // add "id" to the front of the column names
 
@@ -57,7 +59,7 @@ replaceRawKeywordIDs(std::map<int, std::string> raws) {
   return raws;
 }
 
-static inline uint32_t getSolutionCount(IPhreeqcPOET &phreeqc,
+static inline uint32_t getSolutionCount(std::unique_ptr<IPhreeqcPOET> &phreeqc,
                                         const Rcpp::List &initial_grid) {
   IPhreeqcPOET::ModulesArray mod_array;
   Rcpp::Function unique_R("unique");
@@ -73,7 +75,7 @@ static inline uint32_t getSolutionCount(IPhreeqcPOET &phreeqc,
 
   // std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin());
 
-  return phreeqc.getModuleSizes(row_ids)[POET_SOL];
+  return phreeqc->getModuleSizes(row_ids)[POET_SOL];
 }
 
 static std::string readFile(const std::string &path) {
@@ -95,7 +97,7 @@ static std::string readFile(const std::string &path) {
   return buffer.str();
 }
 
-void InitialList::initGrid(const Rcpp::List &grid_input) {
+void InitialList::prepareGrid(const Rcpp::List &grid_input) {
   // parse input values
   Rcpp::Function unique_R("unique");
 
@@ -153,7 +155,7 @@ void InitialList::initGrid(const Rcpp::List &grid_input) {
     throw std::runtime_error("Grid size must be positive.");
   }
 
-  IPhreeqcPOET phreeqc(database, script);
+  this->phreeqc = std::make_unique<IPhreeqcPOET>(database, script);
 
   this->phreeqc_mat = pqcScriptToGrid(phreeqc, R);
   this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
@@ -169,7 +171,7 @@ void InitialList::initGrid(const Rcpp::List &grid_input) {
 
   std::map<int, std::string> pqc_raw_dumps;
 
-  pqc_raw_dumps = replaceRawKeywordIDs(phreeqc.raw_dumps());
+  pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps());
 
   this->pqc_ids =
       Rcpp::as<std::vector<int>>(unique_R(this->initial_grid["ID"]));

src/Init/InitialList.cpp

@@ -10,7 +10,7 @@
 
 namespace poet {
 void InitialList::initializeFromList(const Rcpp::List &setup) {
-  initGrid(setup[grid_key]);
+  prepareGrid(setup[grid_key]);
   initDiffusion(setup[diffusion_key]);
   initChemistry(setup[chemistry_key]);
 }

src/Init/InitialList.hpp

@@ -1,9 +1,9 @@
 #pragma once
 
 #include "Base/RInsidePOET.hpp"
-#include "Chemistry/ChemistryDefs.hpp"
 #include "DataStructures/NamedVector.hpp"
 #include <Rcpp/vector/instantiation.h>
+#include <memory>
 #include <set>
 #include <tug/Boundary.hpp>
 
@@ -86,7 +86,10 @@ private:
     return GridMembersString[static_cast<std::size_t>(member)];
   }
 
-  void initGrid(const Rcpp::List &grid_input);
+  std::unique_ptr<IPhreeqcPOET> phreeqc;
+
+  void prepareGrid(const Rcpp::List &grid_input);
+
   std::uint8_t dim;
 
   std::uint32_t n_cols;