BREAKING CHANGE: substitute R bindings of Phreeqc by PhreeqcRM

app/poet.cpp

@@ -131,9 +131,11 @@ int main(int argc, char *argv[]) {
 
   MPI_Barrier(MPI_COMM_WORLD);
 
+  poet::ChemistryParams chem_params(R);
+
   /* THIS IS EXECUTED BY THE MASTER */
   if (world_rank == 0) {
-    ChemMaster master(params, R, grid);
+    ChemMaster master(params, R, grid, chem_params);
     DiffusionModule diffusion(poet::DiffusionParams(R), grid);
 
     // diffusion.initialize();
@@ -217,7 +219,7 @@ int main(int argc, char *argv[]) {
   }
   /* THIS IS EXECUTED BY THE WORKERS */
   else {
-    ChemWorker worker(params, R, grid, dht_comm);
+    ChemWorker worker(params, R, grid, dht_comm, chem_params);
     worker.loop();
   }
 

data/CMakeLists.txt

@@ -2,5 +2,7 @@ install(
 FILES
     SimDol2D_diffu.R
     SimDol1D_diffu.R
+    phreeqc_kin.dat
+    dol.pqi
 DESTINATION
     data)

data/SimDol2D_diffu.R

@@ -1,10 +1,13 @@
+database <- normalizePath("../data/phreeqc_kin.dat")
+input_script <- normalizePath("../data/dol.pqi")
+
 #################################################################
 ##                          Section 1                          ##
 ##                     Grid initialization                     ##
 #################################################################
 
-n <- 50
-m <- 50
+n <- 5
+m <- 5
 
 types <- c("scratch", "phreeqc", "rds")
 
@@ -28,28 +31,30 @@ types <- c("scratch", "phreeqc", "rds")
 #              "Dolomite 0.0 0.0",
 #              "END")
 
-# needed if init type is set to "scratch" 
+# needed if init type is set to "scratch"
 # prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")
 
 init_cell <- list(
+  "H" = 110.683,
+  "O" = 55.3413,
+  "Charge" = -5.0822e-19,
   "C" = 1.2279E-4,
   "Ca" = 1.2279E-4,
   "Cl" = 0,
   "Mg" = 0,
-  "pH" = 9.91,
-  "pe" = 4,
-  "O2g" = 10,
+  "O2g" = 0.499957,
   "Calcite" = 2.07e-4,
   "Dolomite" = 0
 )
 
 grid <- list(
   n_cells = c(n, m),
-  s_cells = c(n,m),
+  s_cells = c(n, m),
   type = types[1],
   init_cell = as.data.frame(init_cell),
   props = names(init_cell),
-  init_script = NULL
+  database = database,
+  input_script = input_script
 )
 
 
@@ -59,31 +64,34 @@ grid <- list(
 ##################################################################
 
 init_diffu <- c(
+  "H" = 110.683,
+  "O" = 55.3413,
+  "Charge" = -5.0822e-19,
   "C" = 1.2279E-4,
   "Ca" = 1.2279E-4,
   "Cl" = 0,
-  "Mg" = 0,
-  "pe" = 4,
-  "pH" = 7
+  "Mg" = 0
 )
 
 alpha_diffu <- c(
+  "H" = 1E-4,
+  "O" = 1E-4,
+  "Charge" = 1E-4,
   "C" = 1E-4,
   "Ca" = 1E-4,
   "Cl" = 1E-4,
-  "Mg" = 1E-4,
-  "pe" = 1E-4,
-  "pH" = 1E-4
+  "Mg" = 1E-4
 )
 
 vecinj_diffu <- list(
   list(
+    "H" = 110.683,
+    "O" = 55.3413,
+    "Charge" = 1.90431e-16,
     "C" = 0,
     "Ca" = 0,
     "Cl" = 0.002,
-    "Mg" = 0.001,
-    "pe" = 4,
-    "pH" = 9.91
+    "Mg" = 0.001
   )
 )
 
@@ -104,21 +112,21 @@ diffu_list <- names(alpha_diffu)
 diffusion <- list(
   init = init_diffu,
   vecinj = do.call(rbind.data.frame, vecinj_diffu),
-  #vecinj_inner = vecinj_inner,
+  # vecinj_inner = vecinj_inner,
   vecinj_index = boundary,
   alpha = alpha_diffu
 )
 
-##################################################################
-##                          Section 3                           ##
-##                      Phreeqc simulation                      ##
-##################################################################
+#################################################################
+##                          Section 3                          ##
+##                  Chemitry module (Phreeqc)                  ##
+#################################################################
 
-db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
-  package = "RedModRphree"
-), is.db = TRUE)
-
-phreeqc::phrLoadDatabaseString(db)
+# db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
+#  package = "RedModRphree"
+# ), is.db = TRUE)
+#
+# phreeqc::phrLoadDatabaseString(db)
 
 # NOTE: This won't be needed in the future either. Could also be done in a. pqi
 # file
@@ -154,11 +162,23 @@ selout <- c(
   "-kinetic_reactants Calcite Dolomite", "-equilibrium O2g"
 )
 
+
 # Needed when using DHT
 signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5)
 prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act")
 prop <- names(init_cell)
 
+chemistry <- list(
+  database = database,
+  input_script = input_script
+)
+
+#################################################################
+##                          Section 4                          ##
+##              Putting all those things together              ##
+#################################################################
+
+
 iterations <- 10
 
 setup <- list(
@@ -169,6 +189,7 @@ setup <- list(
   # Cf = 1,
   grid = grid,
   diffusion = diffusion,
+  chemistry = chemistry,
   prop = prop,
   immobile = c(7, 8, 9),
   kin = c(8, 9),

data/dol.pqi

@@ -0,0 +1,35 @@
+SELECTED_OUTPUT
+         -high_precision true
+         -reset false
+         -time
+         -soln
+         -temperature true
+         -water true
+         -pH
+         -pe
+         -totals C Ca Cl Mg
+         -kinetic_reactants Calcite Dolomite
+         -equilibrium O2g
+
+SOLUTION 1
+         units mol/kgw
+         temp 25.0
+         water 1
+         pH 9.91 charge
+         pe 4.0
+         C   1.2279E-04
+         Ca  1.2279E-04
+         Cl 1E-12
+         Mg 1E-12
+PURE 1
+         O2g -0.1675 10
+KINETICS 1
+         Calcite
+            -m      0.000207
+            -parms  0.0032
+            -tol    1e-10
+         Dolomite
+            -m      0.0
+            -parms  0.00032
+            -tol    1e-10
+END

ext/phreeqcrm

@@ -1 +1 @@
-Subproject commit 12cc437f6737821f70bcd5d93b7aa871ca2fde07
+Subproject commit 0c160a810aeb64c359531c7f84911b9c1df42a9d

include/poet/ChemSimPar.hpp

@@ -26,6 +26,7 @@
 #include "Grid.hpp"
 #include "RInside.h"
 #include "SimParams.hpp"
+#include "poet/PhreeqcWrapper.hpp"
 #include <PhreeqcRM.h>
 #include <array>
 #include <bits/stdint-uintn.h>
@@ -72,7 +73,8 @@ public:
    * @param R_ R runtime
    * @param grid_ Grid object
    */
-  ChemMaster(SimParams &params, RInside &R_, Grid &grid_);
+  ChemMaster(SimParams &params, RInside &R_, Grid &grid_,
+             poet::ChemistryParams chem_params);
 
   /**
    * @brief Destroy the ChemMaster object
@@ -217,7 +219,8 @@ public:
    * @param grid_ Grid object
    * @param dht_comm Communicator addressing all processes marked as worker
    */
-  ChemWorker(SimParams &params, RInside &R_, Grid &grid_, MPI_Comm dht_comm);
+  ChemWorker(SimParams &params, RInside &R_, Grid &grid_, MPI_Comm dht_comm,
+             poet::ChemistryParams chem_params);
 
   /**
    * @brief Destroy the ChemWorker object
@@ -257,6 +260,8 @@ private:
   int phreeqc_count = 0;
 
   double *mpi_buffer;
+
+  uint32_t ncomps;
 };
 } // namespace poet
 #endif // CHEMSIMPAR_H

include/poet/ChemSimSeq.hpp

@@ -23,6 +23,7 @@
 
 #include "DHT_Wrapper.hpp"
 #include "Grid.hpp"
+#include "PhreeqcWrapper.hpp"
 #include "RInside.h"
 #include "SimParams.hpp"
 #include <bits/stdint-uintn.h>
@@ -69,7 +70,8 @@ public:
    * @param R_ R runtime
    * @param grid_ Initialized grid
    */
-  ChemSim(SimParams &params, RInside &R_, Grid &grid_);
+  ChemSim(SimParams &params, RInside &R_, Grid &grid_,
+          poet::ChemistryParams chem_params);
 
   /**
    * @brief Run iteration of simulation in sequential mode
@@ -119,6 +121,9 @@ protected:
   uint32_t n_cells_per_prop;
   std::vector<std::string> prop_names;
   std::vector<uint32_t> wp_sizes_vector;
+
+  PhreeqcWrapper *phreeqc_rm;
+
 };
 
 } // namespace poet

include/poet/PhreeqcWrapper.hpp

@@ -0,0 +1,40 @@
+#ifndef PHREEQCWRAPPER_H_
+#define PHREEQCWRAPPER_H_
+
+#include "IrmResult.h"
+#include "poet/SimParams.hpp"
+#include <PhreeqcRM.h>
+#include <bits/stdint-uintn.h>
+#include <string>
+#include <vector>
+namespace poet {
+class PhreeqcWrapper : public PhreeqcRM {
+public:
+  PhreeqcWrapper(uint32_t inxyz);
+
+  void SetupAndLoadDB(const poet::ChemistryParams &chemPar);
+  void InitFromFile(const std::string &strInputFile);
+
+  auto GetSpeciesCount() -> uint32_t;
+  auto GetSpeciesCountByModule() -> const std::vector<uint32_t> &;
+  auto GetSpeciesNamesFull() -> const std::vector<std::string> &;
+
+  void SetInternalsFromWP(const std::vector<double> &vecWP,
+                               uint32_t iCurrWPSize);
+  auto GetWPFromInternals(uint32_t iCurrWPSize) -> std::vector<double>;
+
+private:
+  void InitInternals();
+  void SetMappingForWP(uint32_t iCurrWPSize);
+
+  static constexpr int MODULE_COUNT = 5;
+  uint32_t iWPSize;
+  bool bInactiveCells = false;
+  uint32_t iSpeciesCount;
+  std::vector<uint32_t> vecSpeciesPerModule;
+  std::vector<std::string> vecSpeciesNames;
+  std::vector<int> vecDefMapping;
+};
+} // namespace poet
+
+#endif // PHREEQCWRAPPER_H_

include/poet/SimParams.hpp

@@ -75,7 +75,8 @@ using GridParams = struct s_GridParams {
   std::string type;
 
   Rcpp::DataFrame init_df;
-  std::string init_sim;
+  std::string input_script;
+  std::string database_path;
 
   std::vector<std::string> props;
 
@@ -94,6 +95,13 @@ using DiffusionParams = struct s_DiffusionParams {
   s_DiffusionParams(RInside &R);
 };
 
+using ChemistryParams = struct s_ChemistryParams {
+  std::string database_path;
+  std::string input_script;
+
+  s_ChemistryParams(RInside &R);
+};
+
 /**
  * @brief Reads information from program arguments and R runtime
  *

src/ChemMaster.cpp

@@ -19,6 +19,7 @@
 */
 
 #include "poet/DiffusionModule.hpp"
+#include "poet/SimParams.hpp"
 #include <Rcpp.h>
 
 #include <bits/stdint-uintn.h>
@@ -33,8 +34,9 @@ using namespace poet;
 using namespace std;
 using namespace Rcpp;
 
-ChemMaster::ChemMaster(SimParams &params, RInside &R_, Grid &grid_)
-    : ChemSim(params, R_, grid_) {
+ChemMaster::ChemMaster(SimParams &params, RInside &R_, Grid &grid_,
+                       poet::ChemistryParams chem_params)
+    : ChemSim(params, R_, grid_, chem_params) {
   t_simparams tmp = params.getNumParams();
 
   this->wp_size = tmp.wp_size;

src/ChemSim.cpp

@@ -19,6 +19,7 @@
 */
 
 #include "poet/DiffusionModule.hpp"
+#include "poet/SimParams.hpp"
 #include <poet/ChemSimSeq.hpp>
 #include <poet/Grid.hpp>
 
@@ -33,7 +34,8 @@
 using namespace Rcpp;
 using namespace poet;
 
-ChemSim::ChemSim(SimParams &params, RInside &R_, Grid &grid_)
+ChemSim::ChemSim(SimParams &params, RInside &R_, Grid &grid_,
+                 poet::ChemistryParams chem_params)
     : R(R_), grid(grid_) {
   t_simparams tmp = params.getNumParams();
   this->world_rank = tmp.world_rank;
@@ -56,6 +58,10 @@ ChemSim::ChemSim(SimParams &params, RInside &R_, Grid &grid_)
     std::copy(prop_vec.begin(), prop_vec.end(),
               field.begin() + (i * this->n_cells_per_prop));
   }
+  this->phreeqc_rm = new PhreeqcWrapper(this->n_cells_per_prop);
+
+  this->phreeqc_rm->SetupAndLoadDB(chem_params);
+  this->phreeqc_rm->InitFromFile(chem_params.input_script);
 }
 
 void ChemSim::simulate(double dt) {
@@ -86,20 +92,23 @@ void ChemSim::simulate(double dt) {
                std::to_string(this->n_cells_per_prop) +
                ")), mysetup$grid$props)");
 
-  R.parseEvalQ(
-      "result <- slave_chemistry(setup=mysetup, data=mysetup$state_T)");
-  R.parseEvalQ("mysetup <- master_chemistry(setup=mysetup, data=result)");
+  this->phreeqc_rm->SetInternalsFromWP(field, this->n_cells_per_prop);
+  this->phreeqc_rm->SetTime(0);
+  this->phreeqc_rm->SetTimeStep(dt);
+  this->phreeqc_rm->RunCells();
 
-  // HACK: copy R data structure back to C++ field
-  Rcpp::DataFrame result = R.parseEval("mysetup$state_C");
+  // HACK: we will copy resulting field into global grid field. Maybe we will
+  // find a more performant way here ...
+  std::vector<double> vecSimResult =
+      this->phreeqc_rm->GetWPFromInternals(this->n_cells_per_prop);
+  std::copy(vecSimResult.begin(), vecSimResult.end(), field.begin());
 
-  for (uint32_t i = 0; i < this->prop_names.size(); i++) {
-    std::vector<double> c_prop =
-        Rcpp::as<std::vector<double>>(result[this->prop_names[i].c_str()]);
+  R["TMP_C"] = field;
 
-    std::copy(c_prop.begin(), c_prop.end(),
-              field.begin() + (i * this->n_cells_per_prop));
-  }
+  R.parseEvalQ("mysetup$state_C <- setNames(data.frame(matrix(TMP_C, "
+               "ncol=length(mysetup$grid$props), nrow=" +
+               std::to_string(this->n_cells_per_prop) +
+               ")), mysetup$grid$props)");
 
   /* end time measuring */
   chem_b = MPI_Wtime();

src/ChemWorker.cpp

@@ -33,8 +33,8 @@ using namespace std;
 using namespace Rcpp;
 
 ChemWorker::ChemWorker(SimParams &params, RInside &R_, Grid &grid_,
-                       MPI_Comm dht_comm)
-    : ChemSim(params, R_, grid_) {
+                       MPI_Comm dht_comm, poet::ChemistryParams chem_args)
+    : ChemSim(params, R_, grid_, chem_args) {
   t_simparams tmp = params.getNumParams();
 
   this->dt_differ = tmp.dt_differ;

src/PhreeqcWrapper.cpp

@@ -0,0 +1,215 @@
+#include <algorithm>
+#include <bits/stdint-uintn.h>
+#include <iterator>
+#include <poet/PhreeqcWrapper.hpp>
+#include <string>
+#include <vector>
+
+poet::PhreeqcWrapper::PhreeqcWrapper(uint32_t inxyz)
+    : PhreeqcRM(inxyz, 1), iWPSize(inxyz) {
+  this->vecDefMapping.reserve(inxyz);
+
+  for (uint32_t i = 0; i < inxyz; i++) {
+    this->vecDefMapping[i] = i;
+  }
+}
+
+void poet::PhreeqcWrapper::SetupAndLoadDB(
+    const poet::ChemistryParams &chemPar) {
+  // TODO: hardcoded options ...
+  this->SetErrorHandlerMode(1);
+  this->SetComponentH2O(false);
+  this->SetRebalanceFraction(0.5);
+  this->SetRebalanceByCell(true);
+  this->UseSolutionDensityVolume(false);
+  this->SetPartitionUZSolids(false);
+
+  // Set concentration units
+  // 1, mg/L; 2, mol/L; 3, kg/kgs
+  this->SetUnitsSolution(2);
+  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+  this->SetUnitsPPassemblage(1);
+  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+  this->SetUnitsExchange(1);
+  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+  this->SetUnitsSurface(1);
+  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+  this->SetUnitsGasPhase(1);
+  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+  this->SetUnitsSSassemblage(1);
+  // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
+  this->SetUnitsKinetics(1);
+
+  // Set representative volume
+  std::vector<double> rv;
+  rv.resize(this->iWPSize, 1.0);
+  this->SetRepresentativeVolume(rv);
+
+  // Set initial porosity
+  std::vector<double> por;
+  por.resize(this->iWPSize, 0.05);
+  this->SetPorosity(por);
+
+  // Set initial saturation
+  std::vector<double> sat;
+  sat.resize(this->iWPSize, 1.0);
+  this->SetSaturation(sat);
+
+  // Load database
+  this->LoadDatabase(chemPar.database_path);
+}
+
+void poet::PhreeqcWrapper::InitFromFile(const std::string &strInputFile) {
+  this->RunFile(true, true, true, strInputFile);
+  this->RunString(true, false, true, "DELETE; -all");
+
+  this->FindComponents();
+
+  std::vector<int> ic1;
+  ic1.resize(this->iWPSize * 7, -1);
+  // TODO: hardcoded reaction modules
+  for (int i = 0; i < nxyz; i++) {
+    ic1[i] = 1;             // Solution 1
+    ic1[nxyz + i] = 1;      // Equilibrium 1
+    ic1[2 * nxyz + i] = -1; // Exchange none
+    ic1[3 * nxyz + i] = -1; // Surface none
+    ic1[4 * nxyz + i] = -1; // Gas phase none
+    ic1[5 * nxyz + i] = -1; // Solid solutions none
+    ic1[6 * nxyz + i] = 1;  // Kinetics 1
+  }
+
+  this->InitialPhreeqc2Module(ic1);
+
+  // Initial equilibration of cells
+  double time = 0.0;
+  double time_step = 0.0;
+  this->SetTime(time);
+  this->SetTimeStep(time_step);
+  this->RunCells();
+
+  this->InitInternals();
+}
+
+void poet::PhreeqcWrapper::InitInternals() {
+  uint32_t sum = 0, curr_count;
+  std::vector<std::string> vecCurrSpecies;
+
+  this->vecSpeciesPerModule.reserve(this->MODULE_COUNT);
+
+  vecCurrSpecies = this->GetComponents();
+  this->vecSpeciesNames.insert(this->vecSpeciesNames.end(),
+                               vecCurrSpecies.begin(), vecCurrSpecies.end());
+  this->vecSpeciesPerModule.push_back(vecCurrSpecies.size());
+
+  vecCurrSpecies = this->GetEquilibriumPhases();
+  this->vecSpeciesNames.insert(this->vecSpeciesNames.end(),
+                               vecCurrSpecies.begin(), vecCurrSpecies.end());
+  this->vecSpeciesPerModule.push_back(vecCurrSpecies.size());
+
+  vecCurrSpecies = this->GetExchangeSpecies();
+  this->vecSpeciesNames.insert(this->vecSpeciesNames.end(),
+                               vecCurrSpecies.begin(), vecCurrSpecies.end());
+  this->vecSpeciesPerModule.push_back(vecCurrSpecies.size());
+
+  vecCurrSpecies = this->GetSurfaceSpecies();
+  this->vecSpeciesNames.insert(this->vecSpeciesNames.end(),
+                               vecCurrSpecies.begin(), vecCurrSpecies.end());
+  this->vecSpeciesPerModule.push_back(vecCurrSpecies.size());
+
+  vecCurrSpecies = this->GetKineticReactions();
+  this->vecSpeciesNames.insert(this->vecSpeciesNames.end(),
+                               vecCurrSpecies.begin(), vecCurrSpecies.end());
+  this->vecSpeciesPerModule.push_back(vecCurrSpecies.size());
+
+  this->iSpeciesCount = this->vecSpeciesNames.size();
+}
+
+auto poet::PhreeqcWrapper::GetSpeciesCount() -> uint32_t {
+  return this->iSpeciesCount;
+}
+
+auto poet::PhreeqcWrapper::GetSpeciesCountByModule()
+    -> const std::vector<uint32_t> & {
+  return this->vecSpeciesPerModule;
+}
+
+auto poet::PhreeqcWrapper::GetSpeciesNamesFull()
+    -> const std::vector<std::string> & {
+  return this->vecSpeciesNames;
+}
+
+void poet::PhreeqcWrapper::SetInternalsFromWP(const std::vector<double> &vecWP,
+                                              uint32_t iCurrWPSize) {
+  uint32_t iCurrElements;
+
+  if (iCurrWPSize != this->iWPSize) {
+    this->SetMappingForWP(iCurrWPSize);
+  }
+
+  auto itStart = vecWP.begin();
+  auto itEnd = itStart;
+
+  iCurrElements = this->vecSpeciesPerModule[0];
+
+  itEnd += iCurrElements * iCurrWPSize;
+  this->SetConcentrations(std::vector<double>(itStart, itEnd));
+  itStart = itEnd;
+
+  // Equlibirum Phases
+  if ((iCurrElements = this->vecSpeciesPerModule[1]) != 0) {
+    itEnd += iCurrElements * iCurrWPSize;
+    this->SetPPhaseMoles(std::vector<double>(itStart, itEnd));
+    itStart = itEnd;
+  }
+
+  // // NOTE: Block for 'Surface' and 'Exchange' is missing because of missing
+  // // setters @ PhreeqcRM
+  // // ...
+  // // BLOCK_END
+
+  if ((iCurrElements = this->vecSpeciesPerModule[4]) != 0) {
+    itEnd += iCurrElements * iCurrWPSize;
+    this->SetKineticsMoles(std::vector<double>(itStart, itEnd));
+    itStart = itEnd;
+  }
+}
+auto poet::PhreeqcWrapper::GetWPFromInternals(uint32_t iCurrWPSize)
+    -> std::vector<double> {
+  std::vector<double> vecOutput, vecCurrOutput;
+
+  vecOutput.reserve(this->iSpeciesCount * iCurrWPSize);
+
+  this->GetConcentrations(vecCurrOutput);
+  vecOutput.insert(vecOutput.end(), vecCurrOutput.begin(), vecCurrOutput.end());
+
+  if (this->vecSpeciesPerModule[1] != 0) {
+    this->GetPPhaseMoles(vecCurrOutput);
+    vecOutput.insert(vecOutput.end(), vecCurrOutput.begin(),
+                     vecCurrOutput.end());
+  }
+
+  // NOTE: Block for 'Surface' and 'Exchange' is missing because of missing
+  // Getters @ PhreeqcRM
+  // ...
+  // BLOCK_END
+
+  if (this->vecSpeciesPerModule[4] != 0) {
+    this->GetKineticsMoles(vecCurrOutput);
+    vecOutput.insert(vecOutput.end(), vecCurrOutput.begin(),
+                     vecCurrOutput.end());
+  }
+
+  return vecOutput;
+}
+
+void poet::PhreeqcWrapper::SetMappingForWP(uint32_t iCurrWPSize) {
+  if (iCurrWPSize == this->iWPSize) {
+    this->CreateMapping(this->vecDefMapping);
+  } else {
+    std::vector<int> vecCurrMapping = this->vecDefMapping;
+    for (uint32_t i = iCurrWPSize; i < vecCurrMapping.size(); i++) {
+      vecCurrMapping[i] = -1;
+    }
+    this->CreateMapping(vecCurrMapping);
+  }
+}

src/SimParams.cpp

@@ -42,8 +42,10 @@ poet::GridParams::s_GridParams(RInside &R) {
       Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$grid$init_cell"));
   this->props =
       Rcpp::as<std::vector<std::string>>(R.parseEval("mysetup$grid$props"));
-  // this->init_sim =
-  // Rcpp::as<std::string>(R.parseEval("mysetup$grid$init_sim"));
+  this->input_script =
+  Rcpp::as<std::string>(R.parseEval("mysetup$grid$input_script"));
+  this->database_path =
+  Rcpp::as<std::string>(R.parseEval("mysetup$grid$database"));
 }
 
 poet::DiffusionParams::s_DiffusionParams(RInside &R) {
@@ -61,6 +63,13 @@ poet::DiffusionParams::s_DiffusionParams(RInside &R) {
       Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj_index"));
 }
 
+poet::ChemistryParams::s_ChemistryParams(RInside &R) {
+  this->database_path = Rcpp::as<std::string>(
+      R.parseEval("mysetup$chemistry$database"));
+  this->input_script = Rcpp::as<std::string>(
+      R.parseEval("mysetup$chemistry$input_script"));
+}
+
 SimParams::SimParams(int world_rank_, int world_size_) {
   this->simparams.world_rank = world_rank_;
   this->simparams.world_size = world_size_;