Add DHT functionality, still need to be validated

bench/dolo/het/dolo_200_rt.R

@@ -1,4 +1,4 @@
-iterations <- 500
+iterations <- 100
 dt <- 50
 
 list(

src/Base/RInsidePOET.hpp

@@ -3,6 +3,7 @@
 
 #include <RInside.h>
 #include <Rcpp.h>
+#include <Rinternals.h>
 #include <cstddef>
 #include <exception>
 #include <optional>
@@ -42,6 +43,13 @@ public:
     }
   }
 
+  RHookFunction(SEXP f) {
+    try {
+      this->func = Rcpp::Function(f);
+    } catch (const std::exception &e) {
+    }
+  }
+
   template <typename... Args> T operator()(Args... args) const {
     if (func.has_value()) {
       return (Rcpp::as<T>(this->func.value()(args...)));
@@ -52,6 +60,8 @@ public:
 
   bool isValid() const { return this->func.has_value(); }
 
+  SEXP asSEXP() const { return Rcpp::as<SEXP>(this->func.value()); }
+
 private:
   std::optional<Rcpp::Function> func;
 };

src/Chemistry/ChemistryModule.cpp

@@ -157,7 +157,9 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
 poet::ChemistryModule::ChemistryModule(
     uint32_t wp_size_, const InitialList::ChemistryInit &chem_params,
     MPI_Comm communicator)
-    : params(chem_params), wp_size(wp_size_), group_comm(communicator) {
+    : params(chem_params), wp_size(wp_size_), group_comm(communicator),
+      dht_species(chem_params.dht_species),
+      interp_species(chem_params.interp_species) {
   MPI_Comm_rank(communicator, &comm_rank);
   MPI_Comm_size(communicator, &comm_size);
 
@@ -184,43 +186,42 @@ poet::ChemistryModule::~ChemistryModule() {
 
 void poet::ChemistryModule::initializeDHT(
     uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
-  // constexpr uint32_t MB_FACTOR = 1E6;
+  constexpr uint32_t MB_FACTOR = 1E6;
 
-  // this->dht_enabled = true;
+  this->dht_enabled = true;
 
-  // MPI_Comm dht_comm;
+  MPI_Comm dht_comm;
 
-  // if (this->is_master) {
+  if (this->is_master) {
-  //   MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank,
+    MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank, &dht_comm);
-  //   &dht_comm); return;
+    return;
-  // }
+  }
 
-  // if (!this->is_master) {
+  if (!this->is_master) {
 
-  //   MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
+    MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
 
-  //   auto map_copy = key_species;
+    auto map_copy = key_species;
 
-  //   if (key_species.empty()) {
+    if (key_species.empty()) {
-  //     std::vector<std::uint32_t> default_signif(
+      std::vector<std::uint32_t> default_signif(
-  //         this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
+          this->prop_count, DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
-  //     map_copy = NamedVector<std::uint32_t>(this->prop_names,
+      map_copy = NamedVector<std::uint32_t>(this->prop_names, default_signif);
-  //     default_signif);
+    }
-  //   }
 
-  //   auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
+    auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
 
-  //   if (this->dht) {
+    if (this->dht) {
-  //     delete this->dht;
+      delete this->dht;
-  //   }
+    }
 
-  //   const std::uint64_t dht_size = size_mb * MB_FACTOR;
+    const std::uint64_t dht_size = size_mb * MB_FACTOR;
 
-  //   this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
+    this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
-  //                               this->prop_names, params.hooks,
+                                this->prop_names, params.hooks,
-  //                               this->prop_count, params.use_interp);
+                                this->prop_count, interp_enabled);
-  //   this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
+    this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
-  // }
+  }
 }
 
 inline std::vector<std::int32_t> poet::ChemistryModule::parseDHTSpeciesVec(
@@ -229,6 +230,8 @@ inline std::vector<std::int32_t> poet::ChemistryModule::parseDHTSpeciesVec(
   std::vector<int32_t> species_indices;
   species_indices.reserve(key_species.size());
 
+  const auto test = key_species.getNames();
+
   for (const auto &species : key_species.getNames()) {
     auto it = std::find(to_compare.begin(), to_compare.end(), species);
     if (it == to_compare.end()) {
@@ -293,42 +296,42 @@ void poet::ChemistryModule::initializeInterp(
     std::uint32_t bucket_size, std::uint32_t size_mb, std::uint32_t min_entries,
     const NamedVector<std::uint32_t> &key_species) {
 
-  // if (!this->is_master) {
+  if (!this->is_master) {
 
-  //   constexpr uint32_t MB_FACTOR = 1E6;
+    constexpr uint32_t MB_FACTOR = 1E6;
 
-  //   assert(this->dht);
+    assert(this->dht);
 
-  //   this->interp_enabled = true;
+    this->interp_enabled = true;
 
-  //   auto map_copy = key_species;
+    auto map_copy = key_species;
 
-  //   if (key_species.empty()) {
+    if (key_species.empty()) {
-  //     map_copy = this->dht->getKeySpecies();
+      map_copy = this->dht->getKeySpecies();
-  //     for (std::size_t i = 0; i < map_copy.size(); i++) {
+      for (std::size_t i = 0; i < map_copy.size(); i++) {
-  //       const std::uint32_t signif =
+        const std::uint32_t signif =
-  //           map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
+            map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
-  //                              ? InterpolationModule::COARSE_DIFF
+                               ? InterpolationModule::COARSE_DIFF
-  //                              : 0);
+                               : 0);
-  //       map_copy[i] = signif;
+        map_copy[i] = signif;
-  //     }
+      }
-  //   }
+    }
 
-  //   auto key_indices =
+    auto key_indices =
-  //       parseDHTSpeciesVec(map_copy, dht->getKeySpecies().getNames());
+        parseDHTSpeciesVec(map_copy, dht->getKeySpecies().getNames());
 
-  //   if (this->interp) {
+    if (this->interp) {
-  //     this->interp.reset();
+      this->interp.reset();
-  //   }
+    }
 
-  //   const uint64_t pht_size = size_mb * MB_FACTOR;
+    const uint64_t pht_size = size_mb * MB_FACTOR;
 
-  //   interp = std::make_unique<poet::InterpolationModule>(
+    interp = std::make_unique<poet::InterpolationModule>(
-  //       bucket_size, pht_size, min_entries, *(this->dht), map_copy,
+        bucket_size, pht_size, min_entries, *(this->dht), map_copy, key_indices,
-  //       key_indices, this->prop_names, this->params.hooks);
+        this->prop_names, this->params.hooks);
 
-  //   interp->setInterpolationFunction(inverseDistanceWeighting);
+    interp->setInterpolationFunction(inverseDistanceWeighting);
-  // }
+  }
 }
 
 std::vector<double>

src/Chemistry/ChemistryModule.hpp

@@ -76,22 +76,35 @@ public:
     this->file_pad = std::ceil(std::log10(maxiter + 1));
   }
 
-  /**
+  struct SurrogateSetup {
-   * Create a new worker instance inside given MPI communicator.
+    std::vector<std::string> prop_names;
-   *
-   * Wraps communication needed before instanciation can take place.
-   *
-   * \param communicator MPI communicator to distribute work in.
-   *
-   * \returns A worker instance with fixed work package size.
-   */
-  static ChemistryModule createWorker(MPI_Comm communicator,
-                                      const ChemistryParams &chem_param);
 
-  /**
+    bool dht_enabled;
-   * Default work package size.
+    std::uint32_t dht_size_mb;
-   */
+
-  static constexpr uint32_t CHEM_DEFAULT_WP_SIZE = 5;
+    bool interp_enabled;
+    std::uint32_t interp_bucket_size;
+    std::uint32_t interp_size_mb;
+    std::uint32_t interp_min_entries;
+  };
+
+  void masterEnableSurrogates(const SurrogateSetup &setup) {
+    // FIXME: This is a hack to get the prop_names and prop_count from the setup
+    this->prop_names = setup.prop_names;
+    this->prop_count = setup.prop_names.size();
+
+    this->dht_enabled = setup.dht_enabled;
+    this->interp_enabled = setup.interp_enabled;
+
+    if (this->dht_enabled || this->interp_enabled) {
+      this->initializeDHT(setup.dht_size_mb, this->dht_species);
+    }
+
+    if (this->interp_enabled) {
+      this->initializeInterp(setup.interp_bucket_size, setup.interp_size_mb,
+                             setup.interp_min_entries, this->interp_species);
+    }
+  }
 
   /**
    * Intended to alias input parameters for grid initialization with a single
@@ -327,7 +340,7 @@ protected:
   bool is_sequential;
   bool is_master;
 
-  uint32_t wp_size{CHEM_DEFAULT_WP_SIZE};
+  uint32_t wp_size;
   bool dht_enabled{false};
   int dht_snaps_type{DHT_SNAPS_DISABLED};
   std::string dht_file_out_dir;
@@ -361,6 +374,10 @@ protected:
 
   uint32_t n_cells = 0;
   uint32_t prop_count = 0;
+  std::vector<std::string> prop_names;
+
+  NamedVector<std::uint32_t> dht_species;
+  NamedVector<std::uint32_t> interp_species;
 
   Field chem_field;
 

src/Chemistry/SurrogateModels/DHT_Wrapper.cpp

@@ -22,6 +22,7 @@
 
 #include "DHT_Wrapper.hpp"
 
+#include "Init/InitialList.hpp"
 #include "Rounding.hpp"
 
 #include <algorithm>
@@ -41,7 +42,7 @@ DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
                          const NamedVector<std::uint32_t> &key_species,
                          const std::vector<std::int32_t> &key_indices,
                          const std::vector<std::string> &_output_names,
-                         const ChemistryParams::Chem_Hook_Functions &_hooks,
+                         const InitialList::ChemistryHookFunctions &_hooks,
                          uint32_t data_count, bool _with_interp)
     : key_count(key_indices.size()), data_count(data_count),
       input_key_elements(key_indices), communicator(dht_comm),

src/Chemistry/SurrogateModels/DHT_Wrapper.hpp

@@ -29,6 +29,7 @@
 #include "../../Base/SimParams.hpp"
 #include "Chemistry/ChemistryDefs.hpp"
 
+#include "Init/InitialList.hpp"
 #include "LookupKey.hpp"
 
 #include <array>
@@ -87,7 +88,7 @@ public:
               const NamedVector<std::uint32_t> &key_species,
               const std::vector<std::int32_t> &key_indices,
               const std::vector<std::string> &output_names,
-              const ChemistryParams::Chem_Hook_Functions &hooks,
+              const InitialList::ChemistryHookFunctions &hooks,
               uint32_t data_count, bool with_interp);
   /**
    * @brief Destroy the dht wrapper object
@@ -259,7 +260,7 @@ private:
 
   const std::vector<std::string> &output_names;
 
-  const ChemistryParams::Chem_Hook_Functions &hooks;
+  const InitialList::ChemistryHookFunctions &hooks;
   const bool with_interp;
 
   DHT_ResultObject dht_results;

src/Chemistry/SurrogateModels/Interpolation.hpp

@@ -6,6 +6,7 @@
 #include "DataStructures/NamedVector.hpp"
 
 #include "DHT_Wrapper.hpp"
+#include "Init/InitialList.hpp"
 #include "LookupKey.hpp"
 #include "Rounding.hpp"
 
@@ -162,7 +163,7 @@ public:
                       const NamedVector<std::uint32_t> &interp_key_signifs,
                       const std::vector<std::int32_t> &dht_key_indices,
                       const std::vector<std::string> &out_names,
-                      const ChemistryParams::Chem_Hook_Functions &hooks);
+                      const InitialList::ChemistryHookFunctions &hooks);
 
   enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED };
 
@@ -258,7 +259,7 @@ private:
     return out_key;
   }
 
-  const ChemistryParams::Chem_Hook_Functions &hooks;
+  const InitialList::ChemistryHookFunctions &hooks;
   const std::vector<std::string> &out_names;
   const std::vector<std::string> dht_names;
 };

src/Chemistry/SurrogateModels/InterpolationModule.cpp

@@ -1,4 +1,5 @@
 //  Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
+#include "Init/InitialList.hpp"
 #include "Interpolation.hpp"
 
 #include "DHT_Wrapper.hpp"
@@ -35,7 +36,7 @@ InterpolationModule::InterpolationModule(
     const NamedVector<std::uint32_t> &interp_key_signifs,
     const std::vector<std::int32_t> &dht_key_indices,
     const std::vector<std::string> &_out_names,
-    const ChemistryParams::Chem_Hook_Functions &_hooks)
+    const InitialList::ChemistryHookFunctions &_hooks)
     : dht_instance(dht), key_signifs(interp_key_signifs),
       key_indices(dht_key_indices), min_entries_needed(min_entries_needed),
       dht_names(dht.getKeySpecies().getNames()), out_names(_out_names),

src/Init/ChemistryInit.cpp

@@ -1,13 +1,30 @@
 #include "InitialList.hpp"
 
+#include <Rcpp.h>
+#include <vector>
+
 namespace poet {
 
 void InitialList::initChemistry(const Rcpp::List &chem) {
-  this->dht_defined = chem.containsElementNamed("dht_species");
+  if (chem.containsElementNamed("dht_species")) {
-
-  if (this->dht_defined) {
     this->dht_species = Rcpp::as<NamedVector<uint32_t>>(chem["dht_species"]);
   }
+
+  if (chem.containsElementNamed("pht_species")) {
+    this->interp_species = Rcpp::as<NamedVector<uint32_t>>(chem["pht_species"]);
+  }
+
+  if (chem.containsElementNamed("hooks")) {
+    this->chem_hooks = Rcpp::List(chem["hooks"]);
+
+    std::vector<std::string> hook_names = this->chem_hooks.names();
+
+    for (const auto &name : hook_names) {
+      if (this->hook_name_list.find(name) == this->hook_name_list.end()) {
+        Rcpp::Rcerr << "Unknown chemistry hook: " << name << std::endl;
+      }
+    }
+  }
 }
 
 InitialList::ChemistryInit InitialList::getChemistryInit() const {
@@ -21,8 +38,15 @@ InitialList::ChemistryInit InitialList::getChemistryInit() const {
 
   chem_init.pqc_sol_order = pqc_sol_order;
 
-  chem_init.dht_defined = dht_defined;
   chem_init.dht_species = dht_species;
+  chem_init.interp_species = interp_species;
+
+  if (this->chem_hooks.size() > 0) {
+    chem_init.hooks.dht_fill = this->chem_hooks["dht_fill"];
+    chem_init.hooks.dht_fuzz = this->chem_hooks["dht_fuzz"];
+    chem_init.hooks.interp_pre = this->chem_hooks["interp_pre"];
+    chem_init.hooks.interp_post = this->chem_hooks["interp_post"];
+  }
 
   return chem_init;
 }

src/Init/InitialList.cpp

@@ -1,6 +1,7 @@
 #include "InitialList.hpp"
 #include "DataStructures/NamedVector.hpp"
 #include <Rcpp/internal/wrap.h>
+#include <Rcpp/iostream/Rstreambuf.h>
 #include <Rcpp/proxy/ProtectedProxy.h>
 #include <Rcpp/vector/instantiation.h>
 #include <cstdint>
@@ -54,10 +55,14 @@ void InitialList::importList(const Rcpp::List &setup) {
   this->pqc_sol_order = Rcpp::as<std::vector<std::string>>(
       setup[static_cast<int>(ExportList::CHEM_PQC_SOL_ORDER)]);
 
-  this->dht_defined =
+  this->dht_species = NamedVector<uint32_t>(
-      Rcpp::as<bool>(setup[static_cast<int>(ExportList::CHEM_DHT_DEFINED)]);
-  this->dht_species = Rcpp::as<NamedVector<uint32_t>>(
       setup[static_cast<int>(ExportList::CHEM_DHT_SPECIES)]);
+
+  this->interp_species = Rcpp::as<NamedVector<uint32_t>>(
+      setup[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)]);
+
+  this->chem_hooks =
+      Rcpp::as<Rcpp::List>(setup[static_cast<int>(ExportList::CHEM_HOOKS)]);
 }
 
 Rcpp::List InitialList::exportList() {
@@ -85,10 +90,10 @@ Rcpp::List InitialList::exportList() {
   out[static_cast<int>(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids);
   out[static_cast<int>(ExportList::CHEM_PQC_SOL_ORDER)] =
       Rcpp::wrap(this->pqc_sol_order);
-  out[static_cast<int>(ExportList::CHEM_DHT_DEFINED)] =
+  out[static_cast<int>(ExportList::CHEM_DHT_SPECIES)] = this->dht_species;
-      Rcpp::wrap(this->dht_defined);
+  out[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)] =
-  out[static_cast<int>(ExportList::CHEM_DHT_SPECIES)] =
+      Rcpp::wrap(this->interp_species);
-      Rcpp::wrap(this->dht_species);
+  out[static_cast<int>(ExportList::CHEM_HOOKS)] = this->chem_hooks;
 
   return out;
 }

src/Init/InitialList.hpp

@@ -1,8 +1,10 @@
 #pragma once
 
+#include "Base/RInsidePOET.hpp"
 #include "Chemistry/ChemistryDefs.hpp"
 #include "DataStructures/NamedVector.hpp"
 #include <Rcpp/vector/instantiation.h>
+#include <set>
 #include <tug/Boundary.hpp>
 
 #include <RInside.h>
@@ -51,8 +53,9 @@ private:
     CHEM_PQC_SCRIPTS,
     CHEM_PQC_IDS,
     CHEM_PQC_SOL_ORDER,
-    CHEM_DHT_DEFINED,
     CHEM_DHT_SPECIES,
+    CHEM_INTERP_SPECIES,
+    CHEM_HOOKS,
     ENUM_SIZE
   };
 
@@ -161,10 +164,23 @@ private:
 
   std::vector<std::string> pqc_sol_order;
 
-  bool dht_defined;
   NamedVector<std::uint32_t> dht_species;
 
+  NamedVector<std::uint32_t> interp_species;
+
+  Rcpp::List chem_hooks;
+
+  const std::set<std::string> hook_name_list{"dht_fill", "dht_fuzz",
+                                             "interp_pre", "interp_post"};
+
 public:
+  struct ChemistryHookFunctions {
+    RHookFunction<bool> dht_fill;
+    RHookFunction<std::vector<double>> dht_fuzz;
+    RHookFunction<std::vector<std::size_t>> interp_pre;
+    RHookFunction<bool> interp_post;
+  };
+
   struct ChemistryInit {
     uint32_t total_grid_cells;
 
@@ -173,8 +189,9 @@ public:
     std::vector<int> pqc_ids;
     std::vector<std::string> pqc_sol_order;
 
-    bool dht_defined;
     NamedVector<std::uint32_t> dht_species;
+    NamedVector<std::uint32_t> interp_species;
+    ChemistryHookFunctions hooks;
   };
 
   ChemistryInit getChemistryInit() const;

src/poet.cpp

@@ -110,64 +110,55 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
   cmdl("work-package-size", CHEM_DEFAULT_WORK_PACKAGE_SIZE) >>
       params.work_package_size;
 
-  // /* Parse DHT arguments */
+  /* Parse DHT arguments */
-  // chem_params.use_dht = cmdl["dht"];
+  params.use_dht = cmdl["dht"];
-  // chem_params.use_interp = cmdl["interp"];
+  params.use_interp = cmdl["interp"];
-  // // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
+  // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';
 
-  // cmdl("dht-size", DHT_SIZE_PER_PROCESS_MB) >> chem_params.dht_size;
+  cmdl("dht-size", CHEM_DHT_SIZE_PER_PROCESS_MB) >> params.dht_size;
-  // // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process <<
+  // cout << "CPP: DHT size per process (Byte) = " << dht_size_per_process <<
-  // // endl;
+  // endl;
 
-  // cmdl("dht-snaps", 0) >> chem_params.dht_snaps;
+  cmdl("dht-snaps", 0) >> params.dht_snaps;
 
-  // cmdl("dht-file") >> chem_params.dht_file;
+  params.use_interp = cmdl["interp"];
-  // /*Parse work package size*/
+  cmdl("interp-size", 100) >> params.interp_size;
-  // cmdl("work-package-size", WORK_PACKAGE_SIZE_DEFAULT) >> simparams.wp_size;
+  cmdl("interp-min", 5) >> params.interp_min_entries;
+  cmdl("interp-bucket-entries", 20) >> params.interp_bucket_entries;
 
-  // cmdl("interp-size", 100) >> chem_params.pht_size;
+  /*Parse output options*/
-  // cmdl("interp-min", 5) >> chem_params.interp_min_entries;
-  // cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries;
-
-  // /*Parse output options*/
   // simparams.store_result = !cmdl["ignore-result"];
 
-  // chem_params.use_interp = cmdl["interp"];
+  /*Parse output options*/
-  // cmdl("interp-size", 100) >> chem_params.pht_size;
-  // cmdl("interp-min", 5) >> chem_params.interp_min_entries;
-  // cmdl("interp-bucket-entries", 20) >> chem_params.pht_max_entries;
-
-  // /*Parse output options*/
   // simparams.store_result = !cmdl["ignore-result"];
 
-  // if (simparams.world_rank == 0) {
+  if (MY_RANK == 0) {
-  //   MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
+    // MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
-  //   MSG("Work Package Size: " + std::to_string(simparams.wp_size));
+    MSG("Work Package Size: " + std::to_string(params.work_package_size));
-  //   MSG("DHT is " + BOOL_PRINT(chem_params.use_dht));
+    MSG("DHT is " + BOOL_PRINT(params.use_dht));
 
-  //   if (chem_params.use_dht) {
+    if (params.use_dht) {
-  //     MSG("DHT strategy is " + std::to_string(simparams.dht_strategy));
+      // MSG("DHT strategy is " + std::to_string(simparams.dht_strategy));
-  //     // MDL: these should be outdated (?)
+      // MDL: these should be outdated (?)
-  //     // MSG("DHT key default digits (ignored if 'signif_vector' is "
+      // MSG("DHT key default digits (ignored if 'signif_vector' is "
-  //     // 	"defined) = "
+      // 	"defined) = "
-  //     // 	 << simparams.dht_significant_digits);
+      // 	 << simparams.dht_significant_digits);
-  //     // MSG("DHT logarithm before rounding: "
+      // MSG("DHT logarithm before rounding: "
-  //     // 	 << (simparams.dht_log ? "ON" : "OFF"));
+      // 	 << (simparams.dht_log ? "ON" : "OFF"));
-  //     MSG("DHT size per process (Megabyte) = " +
+      MSG("DHT size per process (Megabyte) = " +
-  //         std::to_string(chem_params.dht_size));
+          std::to_string(params.dht_size));
-  //     MSG("DHT save snapshots is " + BOOL_PRINT(chem_params.dht_snaps));
+      MSG("DHT save snapshots is " + BOOL_PRINT(params.dht_snaps));
-  //     MSG("DHT load file is " + chem_params.dht_file);
+      // MSG("DHT load file is " + chem_params.dht_file);
-  //   }
+    }
 
-  //   if (chem_params.use_interp) {
+    if (params.use_interp) {
-  //     MSG("PHT interpolation enabled: " +
+      MSG("PHT interpolation enabled: " + BOOL_PRINT(params.use_interp));
-  //     BOOL_PRINT(chem_params.use_interp)); MSG("PHT interp-size = " +
+      MSG("PHT interp-size = " + std::to_string(params.interp_size));
-  //     std::to_string(chem_params.pht_size)); MSG("PHT interp-min  = " +
+      MSG("PHT interp-min  = " + std::to_string(params.interp_min_entries));
-  //         std::to_string(chem_params.interp_min_entries));
+      MSG("PHT interp-bucket-entries = " +
-  //     MSG("PHT interp-bucket-entries = " +
+          std::to_string(params.interp_bucket_entries));
-  //         std::to_string(chem_params.pht_max_entries));
+    }
-  //   }
+  }
-  // }
 
   std::string init_file;
   std::string runtime_file;
@@ -225,50 +216,6 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
   return ParseRet::PARSER_OK;
 }
 
-// void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
-//                          const std::string &database_path) {
-//   chem.SetErrorHandlerMode(1);
-//   chem.SetComponentH2O(false);
-//   chem.SetRebalanceFraction(0.5);
-//   chem.SetRebalanceByCell(true);
-//   chem.UseSolutionDensityVolume(false);
-//   chem.SetPartitionUZSolids(false);
-
-//   // Set concentration units
-//   // 1, mg/L; 2, mol/L; 3, kg/kgs
-//   chem.SetUnitsSolution(2);
-//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-//   chem.SetUnitsPPassemblage(1);
-//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-//   chem.SetUnitsExchange(1);
-//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-//   chem.SetUnitsSurface(1);
-//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-//   chem.SetUnitsGasPhase(1);
-//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-//   chem.SetUnitsSSassemblage(1);
-//   // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
-//   chem.SetUnitsKinetics(1);
-
-//   // Set representative volume
-//   std::vector<double> rv;
-//   rv.resize(wp_size, 1.0);
-//   chem.SetRepresentativeVolume(rv);
-
-//   // Set initial porosity
-//   std::vector<double> por;
-//   por.resize(wp_size, 1);
-//   chem.SetPorosity(por);
-
-//   // Set initial saturation
-//   std::vector<double> sat;
-//   sat.resize(wp_size, 1.0);
-//   chem.SetSaturation(sat);
-
-//   // Load database
-//   chem.LoadDatabase(database_path);
-// }
-
 static Rcpp::List RunMasterLoop(RInside &R, const RuntimeParameters &params,
                                 DiffusionModule &diffusion,
                                 ChemistryModule &chem) {
@@ -279,20 +226,9 @@ static Rcpp::List RunMasterLoop(RInside &R, const RuntimeParameters &params,
 
   double sim_time = .0;
 
-  // ChemistryModule chem(nxyz_master, params.getNumParams().wp_size, maxiter,
+  if (params.print_progressbar) {
-  //                      params.getChemParams(), MPI_COMM_WORLD);
+    chem.setProgressBarPrintout(true);
-
+  }
-  // set_chem_parameters(chem, nxyz_master,
-  // params.getChemParams().database_path);
-  // chem.RunInitFile(params.getChemParams().input_script);
-
-  // poet::ChemistryModule::SingleCMap init_df =
-  // DFToHashMap(d_params.initial_t);
-  // chem.initializeField(diffusion.getField());
-
-  // if (params.getNumParams().print_progressbar) {
-  chem.setProgressBarPrintout(true);
-  // }
 
   /* SIMULATION LOOP */
 
@@ -405,17 +341,16 @@ static Rcpp::List RunMasterLoop(RInside &R, const RuntimeParameters &params,
 int main(int argc, char *argv[]) {
   double dSimTime, sim_end;
 
-  int world_size, world_rank;
+  int world_size;
 
   MPI_Init(&argc, &argv);
 
   MPI_Comm_size(MPI_COMM_WORLD, &world_size);
-
+  MPI_Comm_rank(MPI_COMM_WORLD, &MY_RANK);
-  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
 
   RInsidePOET &R = RInsidePOET::getInstance();
 
-  if (world_rank == 0) {
+  if (MY_RANK == 0) {
     MSG("Running POET version " + std::string(poet_version));
   }
 
@@ -437,18 +372,29 @@ int main(int argc, char *argv[]) {
   InitialList init_list(R);
   init_list.importList(run_params.init_params);
 
-  MSG("RInside initialized on process " + std::to_string(world_rank));
+  MSG("RInside initialized on process " + std::to_string(MY_RANK));
 
   ChemistryModule chemistry(run_params.work_package_size,
                             init_list.getChemistryInit(), MPI_COMM_WORLD);
 
-  if (world_rank > 0) {
+  const ChemistryModule::SurrogateSetup surr_setup = {
+      init_list.getInitialGrid().GetProps(),
+      run_params.use_dht,
+      run_params.dht_size,
+      run_params.use_interp,
+      run_params.interp_bucket_entries,
+      run_params.interp_size,
+      run_params.interp_min_entries};
+
+  chemistry.masterEnableSurrogates(surr_setup);
+
+  if (MY_RANK > 0) {
 
     chemistry.WorkerLoop();
 
     MPI_Barrier(MPI_COMM_WORLD);
 
-    MSG("finished, cleanup of process " + std::to_string(world_rank));
+    MSG("finished, cleanup of process " + std::to_string(MY_RANK));
 
     MPI_Finalize();
 
@@ -456,7 +402,7 @@ int main(int argc, char *argv[]) {
   }
 
   // R.parseEvalQ("mysetup <- setup");
-  // // if (world_rank == 0) { // get timestep vector from
+  // // if (MY_RANK == 0) { // get timestep vector from
   // // grid_init function ... //
   std::string master_init_code = "mysetup <- master_init(setup=mysetup)";
   R.parseEval(master_init_code);
@@ -464,13 +410,13 @@ int main(int argc, char *argv[]) {
   // run_params.initVectorParams(R);
 
   // MDL: store all parameters
-  if (world_rank == 0) {
+  if (MY_RANK == 0) {
     MSG("Calling R Function to store calling parameters");
     // R.parseEvalQ("StoreSetup(setup=mysetup)");
   }
 
-  if (world_rank == 0) {
+  if (MY_RANK == 0) {
-    MSG("Init done on process with rank " + std::to_string(world_rank));
+    MSG("Init done on process with rank " + std::to_string(MY_RANK));
   }
 
   // MPI_Barrier(MPI_COMM_WORLD);
@@ -501,10 +447,10 @@ int main(int argc, char *argv[]) {
 
   MPI_Barrier(MPI_COMM_WORLD);
 
-  MSG("finished, cleanup of process " + std::to_string(world_rank));
+  MSG("finished, cleanup of process " + std::to_string(MY_RANK));
   MPI_Finalize();
 
-  if (world_rank == 0) {
+  if (MY_RANK == 0) {
     MSG("done, bye!");
   }
 

src/poet.hpp.in

@@ -46,6 +46,8 @@ const std::set<std::string> paramlist{
 
 constexpr uint32_t CHEM_DEFAULT_WORK_PACKAGE_SIZE = 32;
 
+constexpr uint32_t CHEM_DHT_SIZE_PER_PROCESS_MB = 1.5E3;
+
 struct RuntimeParameters {
   std::vector<double> timesteps;
 
@@ -54,6 +56,15 @@ struct RuntimeParameters {
 
   Rcpp::List init_params;
 
+  bool use_dht;
+  std::uint32_t dht_size;
+  std::uint8_t dht_snaps;
+
+  bool use_interp;
+  std::uint32_t interp_size;
+  std::uint32_t interp_min_entries;
+  std::uint32_t interp_bucket_entries;
+
   struct ChemistryParams {
     // std::string database_path;
     // std::string input_script;