Merge branch '6-dht-support-with-poet-dump' into 'main'

Resolve "DHT support with POET dump"

Closes #6

See merge request naaice/poet!11

app/poet.cpp

@@ -70,7 +70,7 @@ void writeFieldsToR(RInside &R, const Field &trans, const Field &chem) {
 void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
                          const std::string &database_path) {
   chem.SetErrorHandlerMode(1);
-  chem.SetComponentH2O(false);
+  chem.SetComponentH2O(true);
   chem.SetRebalanceFraction(0.5);
   chem.SetRebalanceByCell(true);
   chem.UseSolutionDensityVolume(false);
@@ -132,9 +132,10 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
   chem.initializeField(diffusion.getField());
 
   if (params.getNumParams().dht_enabled) {
-    chem.SetDHTEnabled(true, params.getNumParams().dht_size_per_process);
+    chem.SetDHTEnabled(true, params.getNumParams().dht_size_per_process,
-    if (!params.getDHTSignifVector().empty()) {
+                       chem_params.dht_species);
-      chem.SetDHTSignifVector(params.getDHTSignifVector());
+    if (!chem_params.dht_signif.empty()) {
+      chem.SetDHTSignifVector(chem_params.dht_signif);
     }
     if (!params.getDHTPropTypeVector().empty()) {
       chem.SetDHTPropTypeVector(params.getDHTPropTypeVector());

bench/barite/barite.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-04-14 11:02:31 mluebke"
+## Time-stamp: "Last modified 2023-04-24 16:51:23 mluebke"
 
 database <- normalizePath("../share/poet/bench/barite/db_barite.dat")
 input_script <- normalizePath("../share/poet/bench/barite/barite.pqi")
@@ -44,9 +44,12 @@ grid <- list(
 ## initial conditions
 
 init_diffu <- list(
-  "H" = 110.0124,
+  #"H" = 110.0124,
-  "O" = 55.5087,
+  "H" = 0.00000028904,
-  "Charge" = -1.217E-09,
+  #"O" = 55.5087,
+  "O" = 0.000000165205,
+  #"Charge" = -1.217E-09,
+  "Charge" = -3.337E-08,
   "Ba" = 1.E-10,
   "Cl" = 1.E-10,
   "S(6)" = 6.205E-4,
@@ -55,9 +58,12 @@ init_diffu <- list(
 
 injection_diff <- list(
     list(
-        "H" = 111.0124,
+        #"H" = 111.0124,
-        "O" = 55.50622,
+        "H" = 0.0000002890408,
-        "Charge" = -3.337E-08,
+        #"O" = 55.50622,
+        "O" = 0.00002014464,
+        #"Charge" = -3.337E-08,
+        "Charge" = -3.337000004885E-08,
         "Ba" = 0.1,
         "Cl" = 0.2,
         "S(6)"  = 0,
@@ -109,13 +115,16 @@ diffusion <- list(
 
 
 ## # Needed when using DHT
-signif_vector <- c(9, 9, 10, 5, 5, 5, 5, 5, 5)
+
-prop_type <- c("", "", "", "act", "act", "act", "act", "", "")
+dht_species <- names(init_diffu)
-prop <- names(init_cell)
+#dht_signif <- rep(15, length(dht_species))
+dht_signif <- c(10, 10, 3, 5, 5, 5, 5)
 
 chemistry <- list(
   database = database,
-  input_script = input_script
+  input_script = input_script,
+  dht_species = dht_species,
+  dht_signif = dht_signif
 )
 
 #################################################################

bench/dolo_diffu_inner/dolo_diffu_inner.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-04-24 14:11:01 mluebke"
+## Time-stamp: "Last modified 2023-04-24 15:12:02 mluebke"
 
 database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
 input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
@@ -44,8 +44,8 @@ grid <- list(
 
 ## initial conditions
 init_diffu <- list(
-  "H" = 110.683,
+  "H" = 0.000211313883539788,
-  "O" = 55.3413,
+  "O" = 0.00398302904424952,
   "Charge" = -5.0822e-19,
   "C(4)" = 1.2279E-4,
   "Ca" = 1.2279E-4,
@@ -67,8 +67,8 @@ alpha_diffu <- c(
 ## list of boundary conditions/inner nodes
 vecinj_diffu <- list(
     list(
-        "H" = 110.683,
+        "H" = 0.0001540445,
-        "O" = 55.3413,
+        "O" = 0.002148006,
         "Charge" = 1.90431e-16,
         "C(4)" = 0,
         "Ca" = 0,
@@ -76,8 +76,8 @@ vecinj_diffu <- list(
         "Mg" = 0.001
     ),
     list(
-        "H" = 110.683,
+        "H" = 0.0001610193,
-        "O" = 55.3413,
+        "O" = 0.002386934,
         "Charge" = 1.90431e-16,
         "C(4)" = 0,
         "Ca" = 0.0,
@@ -120,13 +120,16 @@ diffusion <- list(
 
 
 ## # Needed when using DHT
-signif_vector <- c(10, 10, 2, 5, 5, 5, 5, 0, 5, 5)
+dht_species <- c("H", "O", "Charge", "C(4)", "Ca", "Cl", "Mg", "Calcite",
-prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
+                 "Dolomite")
-prop <- names(init_cell)
+#dht_signif <- rep(15, length(dht_species))
+dht_signif <- c(10, 10, 3, 5, 5, 5, 5, 5, 5)
 
 chemistry <- list(
   database = database,
-  input_script = input_script
+  input_script = input_script,
+  dht_species = dht_species,
+  dht_signif = dht_signif
 )
 
 #################################################################
@@ -135,7 +138,7 @@ chemistry <- list(
 #################################################################
 
 
-iterations <- 1000
+iterations <- 10
 dt <- 200
 
 setup <- list(
@@ -144,6 +147,5 @@ setup <- list(
   chemistry = chemistry,
   iterations = iterations,
   timesteps = rep(dt, iterations),
-  store_result = TRUE,
+  store_result = TRUE
-  out_save = c(5, iterations)
 )

bench/dolo_diffu_inner/dolo_diffu_inner_large.R

@@ -1,7 +1,7 @@
-## Time-stamp: "Last modified 2023-04-17 12:27:27 mluebke"
+## Time-stamp: "Last modified 2023-04-24 15:43:26 mluebke"
 
 database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
-input_script <- normalizePath("../share/poet/bench/dolo_inner.pqi")
+input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
 
 #################################################################
 ##                          Section 1                          ##
@@ -44,10 +44,10 @@ grid <- list(
 
 ## initial conditions
 init_diffu <- list(
-  "H" = 110.683,
+  "H" = 0.000211313883539788,
-  "O" = 55.3413,
+  "O" = 0.00398302904424952,
   "Charge" = -5.0822e-19,
-  "C" = 1.2279E-4,
+  "C(4)" = 1.2279E-4,
   "Ca" = 1.2279E-4,
   "Cl" = 0,
   "Mg" = 0
@@ -58,7 +58,7 @@ alpha_diffu <- c(
   "H" = 1E-6,
   "O" = 1E-6,
   "Charge" = 1E-6,
-  "C" = 1E-6,
+  "C(4)" = 1E-6,
   "Ca" = 1E-6,
   "Cl" = 1E-6,
   "Mg" = 1E-6
@@ -67,19 +67,19 @@ alpha_diffu <- c(
 ## list of boundary conditions/inner nodes
 vecinj_diffu <- list(
     list(
-        "H" = 110.683,
+        "H" = 0.0001540445,
-        "O" = 55.3413,
+        "O" = 0.002148006,
         "Charge" = 1.90431e-16,
-        "C" = 0,
+        "C(4)" = 0,
         "Ca" = 0,
         "Cl" = 0.002,
         "Mg" = 0.001
     ),
     list(
-        "H" = 110.683,
+        "H" = 0.0001610193,
-        "O" = 55.3413,
+        "O" = 0.002386934,
         "Charge" = 1.90431e-16,
-        "C" = 0,
+        "C(4)" = 0,
         "Ca" = 0.0,
         "Cl" = 0.004,
         "Mg" = 0.002
@@ -120,13 +120,16 @@ diffusion <- list(
 
 
 ## # Needed when using DHT
-signif_vector <- c(10, 10, 2, 5, 5, 5, 5, 0, 5, 5)
+dht_species <- c("H", "O", "Charge", "C(4)", "Ca", "Cl", "Mg", "Calcite",
-prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
+                 "Dolomite")
-prop <- names(init_cell)
+#dht_signif <- rep(15, length(dht_species))
+dht_signif <- c(10, 10, 3, 5, 5, 5, 5, 5, 5)
 
 chemistry <- list(
   database = database,
-  input_script = input_script
+  input_script = input_script,
+  dht_species = dht_species,
+  dht_signif = dht_signif
 )
 
 #################################################################
@@ -145,5 +148,5 @@ setup <- list(
   iterations = iterations,
   timesteps = rep(dt, iterations),
   store_result = TRUE,
-  out_save = c(5, iterations)
+  out_save = seq(5, iterations, by=5)
 )

ext/phreeqcrm

@@ -1 +1 @@
-Subproject commit 6a4c14c2857f0ad3558af3f1e1aa7555a805e8e2
+Subproject commit 6db8938592a487a701b3bca19858ab946b42f1a0

include/poet/ChemistryModule.hpp

@@ -1,3 +1,5 @@
+//  Time-stamp: "Last modified 2023-04-24 14:30:06 mluebke"
+
 #ifndef CHEMISTRYMODULE_H_
 #define CHEMISTRYMODULE_H_
 
@@ -158,8 +160,10 @@ public:
    *
    * \param enable Enables or disables the usage of the DHT.
    * \param size_mb Size in megabyte to allocate for the DHT if enabled.
+   * \param key_species Name of species to be used for key creation.
    */
-  void SetDHTEnabled(bool enable, uint32_t size_mb);
+  void SetDHTEnabled(bool enable, uint32_t size_mb,
+                     const std::vector<std::string> &key_species);
   /**
    * **Master only** Set DHT snapshots to specific mode.
    *
@@ -174,7 +178,7 @@ public:
    * \param signif_vec Vector defining significant digit for each species. Order
    * is defined by prop_type vector (ChemistryModule::GetPropNames).
    */
-  void SetDHTSignifVector(std::vector<uint32_t> signif_vec);
+  void SetDHTSignifVector(std::vector<uint32_t> &signif_vec);
   /**
    * **Master only** Set the DHT rounding type of each species. See
    * DHT_PROP_TYPES enumemartion for explanation.
@@ -343,7 +347,7 @@ protected:
   void WorkerMetricsToMaster(int type);
 
   IRM_RESULT WorkerRunWorkPackage(std::vector<double> &vecWP,
-                                  std::vector<int32_t> &vecMapping,
+                                  std::vector<std::uint32_t> &vecMapping,
                                   double dSimTime, double dTimestep);
 
   std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
@@ -355,6 +359,8 @@ protected:
   void unshuffleField(const std::vector<double> &in_buffer,
                       uint32_t size_per_prop, uint32_t prop_count,
                       uint32_t wp_count, std::vector<double> &out_field);
+  std::vector<std::uint32_t>
+  parseDHTSpeciesVec(const std::vector<std::string> &species_vec) const;
 
   int comm_size, comm_rank;
   MPI_Comm group_comm;
@@ -362,12 +368,12 @@ protected:
   bool is_sequential;
   bool is_master;
 
-  uint32_t wp_size;
+  uint32_t wp_size{CHEM_DEFAULT_WP_SIZE};
-  bool dht_enabled = false;
+  bool dht_enabled{false};
-  int dht_snaps_type = DHT_FILES_DISABLED;
+  int dht_snaps_type{DHT_FILES_DISABLED};
   std::string dht_file_out_dir;
 
-  poet::DHT_Wrapper *dht = std::nullptr_t();
+  poet::DHT_Wrapper *dht = nullptr;
 
   static constexpr uint32_t BUFFER_OFFSET = 5;
 
@@ -383,6 +389,8 @@ protected:
   double seq_t = 0.;
   double send_recv_t = 0.;
 
+  std::array<double, 2> base_totals{0};
+
 #endif
 
   double chem_t = 0.;

include/poet/DHT_Types.hpp

@@ -2,7 +2,7 @@
 #define DHT_TYPES_H_
 
 namespace poet {
-enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_ACT, DHT_TYPE_IGNORE };
+enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_IGNORE, DHT_TYPE_TOTAL };
 }
 
 #endif // DHT_TYPES_H_

include/poet/DHT_Wrapper.hpp

@@ -1,3 +1,5 @@
+//  Time-stamp: "Last modified 2023-04-24 16:23:42 mluebke"
+
 /*
 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
 ** Potsdam)
@@ -21,6 +23,8 @@
 #ifndef DHT_WRAPPER_H
 #define DHT_WRAPPER_H
 
+#include "poet/DHT_Types.hpp"
+#include <array>
 #include <cstdint>
 #include <string>
 #include <vector>
@@ -33,19 +37,21 @@ extern "C" {
 
 namespace poet {
 
-using DHT_Keyelement = struct keyelem {
+struct DHT_SCNotation {
   std::int8_t exp : 8;
   std::int64_t significant : 56;
 };
 
+union DHT_Keyelement {
+  double fp_elemet;
+  DHT_SCNotation sc_notation;
+};
+
 using DHT_ResultObject = struct DHTResobj {
   uint32_t length;
   std::vector<std::vector<DHT_Keyelement>> keys;
   std::vector<std::vector<double>> results;
   std::vector<bool> needPhreeqc;
-
-  void ResultsToMapping(std::vector<int32_t> &curr_mapping);
-  void ResultsToWP(std::vector<double> &curr_wp);
 };
 
 /**
@@ -79,11 +85,14 @@ public:
    *
    * @param dht_comm Communicator which addresses all participating DHT
    * processes
-   * @param buckets_per_process Count of buckets to allocate for each process
+   * @param buckets_per_process Count of buckets to allocate for each
-   * @param key_count Count of key entries
+   * process
+   * @param key_indices Vector indexing elements of one grid cell used
+   * for key creation.
    * @param data_count Count of data entries
    */
-  DHT_Wrapper(MPI_Comm dht_comm, uint32_t dht_size, uint32_t key_count,
+  DHT_Wrapper(MPI_Comm dht_comm, uint32_t dht_size,
+              const std::vector<std::uint32_t> &key_indices,
               uint32_t data_count);
   /**
    * @brief Destroy the dht wrapper object
@@ -112,8 +121,9 @@ public:
    * @param[in,out] work_package Pointer to current work package
    * @param dt Current timestep of simulation
    */
-  auto checkDHT(int length, double dt, const std::vector<double> &work_package)
+  auto checkDHT(int length, double dt, const std::vector<double> &work_package,
-      -> poet::DHT_ResultObject;
+                std::vector<std::uint32_t> &curr_mapping)
+      -> const poet::DHT_ResultObject &;
 
   /**
    * @brief Write simulated values into DHT
@@ -131,8 +141,9 @@ public:
    * outputs of the PHREEQC simulation
    * @param dt Current timestep of simulation
    */
-  void fillDHT(int length, const DHT_ResultObject &dht_check_data,
+  void fillDHT(int length, const std::vector<double> &work_package);
-               const std::vector<double> &results);
+
+  void resultsToWP(std::vector<double> &work_package);
 
   /**
    * @brief Dump current DHT state into file.
@@ -189,6 +200,10 @@ public:
   void SetSignifVector(std::vector<uint32_t> signif_vec);
   void SetPropTypeVector(std::vector<uint32_t> prop_type_vec);
 
+  void setBaseTotals(const std::array<double, 2> &bt) {
+    this->base_totals = bt;
+  }
+
 private:
   uint32_t key_count;
   uint32_t data_count;
@@ -202,11 +217,16 @@ private:
   uint32_t dht_evictions = 0;
 
   std::vector<uint32_t> dht_signif_vector;
-  std::vector<uint32_t> dht_prop_type_vector;
+  std::vector<std::uint32_t> dht_prop_type_vector;
+  std::vector<std::uint32_t> input_key_elements;
 
   static constexpr int DHT_KEY_SIGNIF_DEFAULT = 5;
-  static constexpr int DHT_KEY_SIGNIF_TOTALS = 12;
+  static constexpr int DHT_KEY_SIGNIF_TOTALS = 10;
   static constexpr int DHT_KEY_SIGNIF_CHARGE = 3;
+
+  DHT_ResultObject dht_results;
+
+  std::array<double, 2> base_totals{0};
 };
 } // namespace poet
 

include/poet/SimParams.hpp

@@ -101,6 +101,8 @@ using DiffusionParams = struct s_DiffusionParams {
 using ChemistryParams = struct s_ChemistryParams {
   std::string database_path;
   std::string input_script;
+  std::vector<std::string> dht_species;
+  std::vector<std::uint32_t> dht_signif;
 
   s_ChemistryParams(RInside &R);
 };

src/ChemistryModule/ChemistryModule.cpp

@@ -3,7 +3,6 @@
 #include "poet/DHT_Wrapper.hpp"
 
 #include <algorithm>
-#include <bits/stdint-uintn.h>
 #include <cassert>
 #include <cstdint>
 #include <map>
@@ -146,7 +145,7 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) {
   }
 
   // now sort the new values
-  std::sort(new_solution_names.begin() + 3, new_solution_names.end());
+  std::sort(new_solution_names.begin() + 4, new_solution_names.end());
 
   // and append other processes than solutions
   std::vector<std::string> new_prop_names = new_solution_names;
@@ -184,14 +183,28 @@ void poet::ChemistryModule::initializeField(const Field &trans_field) {
   }
 }
 
-void poet::ChemistryModule::SetDHTEnabled(bool enable, uint32_t size_mb) {
+void poet::ChemistryModule::SetDHTEnabled(
+    bool enable, uint32_t size_mb,
+    const std::vector<std::string> &key_species) {
+
   constexpr uint32_t MB_FACTOR = 1E6;
+  std::vector<std::uint32_t> key_inidices;
 
   if (this->is_master) {
     int ftype = CHEM_DHT_ENABLE;
     PropagateFunctionType(ftype);
     ChemBCast(&enable, 1, MPI_CXX_BOOL);
     ChemBCast(&size_mb, 1, MPI_UINT32_T);
+
+    key_inidices = parseDHTSpeciesVec(key_species);
+    int vec_size = key_inidices.size();
+    ChemBCast(&vec_size, 1, MPI_INT);
+    ChemBCast(key_inidices.data(), key_inidices.size(), MPI_UINT32_T);
+  } else {
+    int vec_size;
+    ChemBCast(&vec_size, 1, MPI_INT);
+    key_inidices.resize(vec_size);
+    ChemBCast(key_inidices.data(), vec_size, MPI_UINT32_T);
   }
 
   this->dht_enabled = enable;
@@ -214,10 +227,39 @@ void poet::ChemistryModule::SetDHTEnabled(bool enable, uint32_t size_mb) {
     const uint32_t dht_size = size_mb * MB_FACTOR;
 
     this->dht =
-        new DHT_Wrapper(dht_comm, dht_size, this->prop_count, this->prop_count);
+        new DHT_Wrapper(dht_comm, dht_size, key_inidices, this->prop_count);
+    // this->dht->setBaseTotals(this->base_totals);
   }
 }
 
+inline std::vector<std::uint32_t> poet::ChemistryModule::parseDHTSpeciesVec(
+    const std::vector<std::string> &species_vec) const {
+  std::vector<uint32_t> species_indices;
+
+  if (species_vec.empty()) {
+    species_indices.resize(this->prop_count);
+
+    int i = 0;
+    std::generate(species_indices.begin(), species_indices.end(),
+                  [&] { return i++; });
+    return species_indices;
+  }
+
+  species_indices.reserve(species_vec.size());
+
+  for (const auto &name : species_vec) {
+    auto it = std::find(this->prop_names.begin(), this->prop_names.end(), name);
+    if (it == prop_names.end()) {
+      throw std::runtime_error(
+          "DHT species name was not found in prop name vector!");
+    }
+    const std::uint32_t index = it - prop_names.begin();
+    species_indices.push_back(index);
+  }
+
+  return species_indices;
+}
+
 void poet::ChemistryModule::SetDHTSnaps(int type, const std::string &out_dir) {
   if (this->is_master) {
     int ftype = CHEM_DHT_SNAPS;
@@ -230,24 +272,28 @@ void poet::ChemistryModule::SetDHTSnaps(int type, const std::string &out_dir) {
     ChemBCast(const_cast<char *>(out_dir.data()), str_size, MPI_CHAR);
   }
 
+  this->dht_snaps_type = type;
   this->dht_file_out_dir = out_dir;
 }
 
 void poet::ChemistryModule::SetDHTSignifVector(
-    std::vector<uint32_t> signif_vec) {
+    std::vector<uint32_t> &signif_vec) {
   if (this->is_master) {
-    if (signif_vec.size() != this->prop_count) {
-      throw std::runtime_error(
-          "Significant vector sizes mismatches prop count.");
-    }
-
     int ftype = CHEM_DHT_SIGNIF_VEC;
     PropagateFunctionType(ftype);
+
+    int data_count = signif_vec.size();
+    ChemBCast(&data_count, 1, MPI_INT);
     ChemBCast(signif_vec.data(), signif_vec.size(), MPI_UINT32_T);
 
     return;
   }
 
+  int data_count;
+  ChemBCast(&data_count, 1, MPI_INT);
+  signif_vec.resize(data_count);
+  ChemBCast(signif_vec.data(), data_count, MPI_UINT32_T);
+
   this->dht->SetSignifVector(signif_vec);
 }
 

src/ChemistryModule/DHT_Wrapper.cpp

@@ -35,41 +35,24 @@
 using namespace poet;
 using namespace std;
 
-inline DHT_Keyelement round_key_element(double value, std::uint32_t signif) {
+inline DHT_SCNotation round_key_element(double value, std::uint32_t signif) {
-  std::int8_t exp = (std::int8_t)std::floor(std::log10(std::fabs(value)));
+  std::int8_t exp =
-  std::int64_t significant = value * std::pow(10, signif - exp - 1);
+      static_cast<std::int8_t>(std::floor(std::log10(std::fabs(value))));
+  std::int64_t significant =
+      static_cast<std::int64_t>(value * std::pow(10, signif - exp - 1));
 
-  DHT_Keyelement elem;
+  DHT_SCNotation elem;
   elem.exp = exp;
   elem.significant = significant;
 
   return elem;
 }
 
-void poet::DHT_ResultObject::ResultsToMapping(
-    std::vector<int32_t> &curr_mapping) {
-  uint32_t iMappingIndex = 0;
-  for (uint32_t i = 0; i < this->length; i++) {
-    if (curr_mapping[i] == -1) {
-      continue;
-    }
-    curr_mapping[i] = (this->needPhreeqc[i] ? iMappingIndex++ : -1);
-  }
-}
-
-void poet::DHT_ResultObject::ResultsToWP(std::vector<double> &curr_wp) {
-  for (uint32_t i = 0; i < this->length; i++) {
-    if (!this->needPhreeqc[i]) {
-      const uint32_t length = this->results[i].end() - this->results[i].begin();
-      std::copy(this->results[i].begin(), this->results[i].end(),
-                curr_wp.begin() + (length * i));
-    }
-  }
-}
-
 DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, uint32_t dht_size,
-                         uint32_t key_count, uint32_t data_count)
+                         const std::vector<std::uint32_t> &key_indices,
-    : key_count(key_count), data_count(data_count) {
+                         uint32_t data_count)
+    : key_count(key_indices.size()), data_count(data_count),
+      input_key_elements(key_indices) {
   // initialize DHT object
   uint32_t key_size = (key_count + 1) * sizeof(DHT_Keyelement);
   uint32_t data_size = data_count * sizeof(double);
@@ -77,45 +60,39 @@ DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, uint32_t dht_size,
   dht_object = DHT_create(dht_comm, buckets_per_process, data_size, key_size,
                           &poet::Murmur2_64A);
 
-  // extract needed values from sim_param struct
-  // t_simparams tmp = params.getNumParams();
-
-  // this->dt_differ = tmp.dt_differ;
-  // this->dht_log = tmp.dht_log;
-
-  // this->dht_signif_vector = params.getDHTSignifVector();
-  // this->dht_prop_type_vector = params.getDHTPropTypeVector();
-
   this->dht_signif_vector.resize(key_size, DHT_KEY_SIGNIF_DEFAULT);
   this->dht_signif_vector[0] = DHT_KEY_SIGNIF_TOTALS;
   this->dht_signif_vector[1] = DHT_KEY_SIGNIF_TOTALS;
   this->dht_signif_vector[2] = DHT_KEY_SIGNIF_CHARGE;
 
-  this->dht_prop_type_vector.resize(key_size, DHT_TYPE_DEFAULT);
+  this->dht_prop_type_vector.resize(key_count, DHT_TYPE_DEFAULT);
+  this->dht_prop_type_vector[0] = DHT_TYPE_TOTAL;
+  this->dht_prop_type_vector[1] = DHT_TYPE_TOTAL;
+  this->dht_prop_type_vector[2] = DHT_TYPE_CHARGE;
 }
 
 DHT_Wrapper::~DHT_Wrapper() {
   // free DHT
   DHT_free(dht_object, NULL, NULL);
-
 }
 auto DHT_Wrapper::checkDHT(int length, double dt,
-                           const std::vector<double> &work_package)
+                           const std::vector<double> &work_package,
-    -> poet::DHT_ResultObject {
+                           std::vector<std::uint32_t> &curr_mapping)
+    -> const poet::DHT_ResultObject & {
 
-  DHT_ResultObject check_data;
+  dht_results.length = length;
+  dht_results.keys.resize(length);
+  dht_results.results.resize(length);
+  dht_results.needPhreeqc.resize(length);
 
-  check_data.length = length;
+  std::vector<std::uint32_t> new_mapping;
-  check_data.keys.resize(length);
-  check_data.results.resize(length);
-  check_data.needPhreeqc.resize(length);
 
   // loop over every grid cell contained in work package
   for (int i = 0; i < length; i++) {
     // point to current grid cell
-    void *key = (void *)&(work_package[i * this->key_count]);
+    void *key = (void *)&(work_package[i * this->data_count]);
-    auto &data = check_data.results[i];
+    auto &data = dht_results.results[i];
-    auto &key_vector = check_data.keys[i];
+    auto &key_vector = dht_results.keys[i];
 
     data.resize(this->data_count);
     key_vector = fuzzForDHT(this->key_count, key, dt);
@@ -125,27 +102,29 @@ auto DHT_Wrapper::checkDHT(int length, double dt,
 
     switch (res) {
     case DHT_SUCCESS:
-      check_data.needPhreeqc[i] = false;
+      dht_results.needPhreeqc[i] = false;
       this->dht_hits++;
       break;
     case DHT_READ_MISS:
-      check_data.needPhreeqc[i] = true;
+      dht_results.needPhreeqc[i] = true;
+      new_mapping.push_back(curr_mapping[i]);
       this->dht_miss++;
       break;
     }
   }
 
-  return check_data;
+  curr_mapping = std::move(new_mapping);
+
+  return dht_results;
 }
 
-void DHT_Wrapper::fillDHT(int length, const DHT_ResultObject &dht_check_data,
+void DHT_Wrapper::fillDHT(int length, const std::vector<double> &work_package) {
-                          const std::vector<double> &results) {
   // loop over every grid cell contained in work package
   for (int i = 0; i < length; i++) {
     // If true grid cell was simulated, needs to be inserted into dht
-    if (dht_check_data.needPhreeqc[i]) {
+    if (dht_results.needPhreeqc[i]) {
-      const auto &key = dht_check_data.keys[i];
+      const auto &key = dht_results.keys[i];
-      void *data = (void *)&(results[i * this->data_count]);
+      void *data = (void *)&(work_package[i * this->data_count]);
       // fuzz data (round, logarithm etc.)
 
       // insert simulated data with fuzzed key into DHT
@@ -159,6 +138,15 @@ void DHT_Wrapper::fillDHT(int length, const DHT_ResultObject &dht_check_data,
   }
 }
 
+void DHT_Wrapper::resultsToWP(std::vector<double> &work_package) {
+  for (int i = 0; i < dht_results.length; i++) {
+    if (!dht_results.needPhreeqc[i]) {
+      std::copy(dht_results.results[i].begin(), dht_results.results[i].end(),
+                work_package.begin() + (data_count * i));
+    }
+  }
+}
+
 int DHT_Wrapper::tableToFile(const char *filename) {
   int res = DHT_to_file(dht_object, filename);
   return res;
@@ -194,25 +182,28 @@ std::vector<DHT_Keyelement> DHT_Wrapper::fuzzForDHT(int var_count, void *key,
   std::vector<DHT_Keyelement> vecFuzz(var_count + 1);
   std::memset(&vecFuzz[0], 0, sizeof(DHT_Keyelement) * var_count);
 
-  unsigned int i = 0;
+  int totals_i = 0;
   // introduce fuzzing to allow more hits in DHT
   // loop over every variable of grid cell
-  for (i = 0; i < (unsigned int)var_count; i++) {
+  for (std::uint32_t i = 0; i < input_key_elements.size(); i++) {
-    double &curr_key = ((double *)key)[i];
+    double curr_key = ((double *)key)[input_key_elements[i]];
     if (curr_key != 0) {
       if (curr_key < zero_val &&
-          this->dht_prop_type_vector[i] == DHT_TYPE_ACT) {
+          this->dht_prop_type_vector[i] == DHT_TYPE_DEFAULT) {
         continue;
       }
       if (this->dht_prop_type_vector[i] == DHT_TYPE_IGNORE) {
         continue;
       }
-      vecFuzz[i] = round_key_element(curr_key, dht_signif_vector[i]);
+      if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL) {
+        curr_key -= base_totals[totals_i++];
+      }
+      vecFuzz[i].sc_notation = round_key_element(curr_key, dht_signif_vector[i]);
     }
   }
   // if timestep differs over iterations set current current time step at the
   // end of fuzzing buffer
-  vecFuzz[var_count] = round_key_element(dt, 55);
+  vecFuzz[var_count].fp_elemet = dt;
 
   return vecFuzz;
 }

src/ChemistryModule/HashFunctions.cpp

@@ -1,4 +1,4 @@
-//  Time-stamp: "Last modified 2023-03-31 15:00:11 mluebke"
+//  Time-stamp: "Last modified 2023-04-24 16:56:18 mluebke"
 /*
 **-----------------------------------------------------------------------------
 ** MurmurHash2 was written by Austin Appleby, and is placed in the public
@@ -26,9 +26,6 @@
 
 #include "poet/HashFunctions.hpp"
 
-#include <cstddef>
-#include <stdexcept>
-
 #if defined(_MSC_VER)
 
 #define BIG_CONSTANT(x) (x)

src/ChemistryModule/WorkerFunctions.cpp

@@ -1,6 +1,10 @@
+//  Time-stamp: "Last modified 2023-04-24 16:55:34 mluebke"
+
 #include "IrmResult.h"
 #include "poet/ChemistryModule.hpp"
 
+#include <algorithm>
+#include <bits/stdint-uintn.h>
 #include <cstddef>
 #include <cstdint>
 #include <iomanip>
@@ -55,12 +59,13 @@ void poet::ChemistryModule::WorkerLoop() {
       uint32_t size_mb;
       ChemBCast(&size_mb, 1, MPI_UINT32_T);
 
-      SetDHTEnabled(enable, size_mb);
+      std::vector<std::string> name_dummy;
+
+      SetDHTEnabled(enable, size_mb, name_dummy);
       break;
     }
     case CHEM_DHT_SIGNIF_VEC: {
-      std::vector<uint32_t> input_vec(this->prop_count);
+      std::vector<uint32_t> input_vec;
-      ChemBCast(input_vec.data(), this->prop_count, MPI_UINT32_T);
 
       SetDHTSignifVector(input_vec);
       break;
@@ -185,46 +190,35 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
   /* 4th double value is currently a placeholder */
   // placeholder = mpi_buffer[count+4];
 
-  // std::vector<double> vecCurrWP(
-  //     mpi_buffer,
-  //     mpi_buffer + (local_work_package_size * this->prop_names.size()));
   vecCurrWP.resize(n_cells_times_props);
-  std::vector<int32_t> vecMappingWP(this->wp_size);
+  std::vector<std::uint32_t> vecMappingWP(local_work_package_size);
 
-  DHT_ResultObject DHT_Results;
+  {
-
+    std::uint32_t i = 0;
-  for (uint32_t i = 0; i < local_work_package_size; i++) {
+    std::generate(vecMappingWP.begin(), vecMappingWP.end(),
-    vecMappingWP[i] = i;
+                  [&] { return i++; });
-  }
-
-  if (local_work_package_size != this->wp_size) {
-    // std::vector<double> vecFiller(
-    //     (this->wp_size - local_work_package_size) * prop_count, 0);
-    // vecCurrWP.insert(vecCurrWP.end(), vecFiller.begin(), vecFiller.end());
-
-    // set all remaining cells to inactive
-    for (int i = local_work_package_size; i < this->wp_size; i++) {
-      vecMappingWP[i] = -1;
-    }
   }
 
   if (dht_enabled) {
     /* check for values in DHT */
     dht_get_start = MPI_Wtime();
-    DHT_Results = dht->checkDHT(local_work_package_size, dt, vecCurrWP);
+    dht->checkDHT(local_work_package_size, dt, vecCurrWP, vecMappingWP);
     dht_get_end = MPI_Wtime();
 
-    DHT_Results.ResultsToMapping(vecMappingWP);
+    // DHT_Results.ResultsToMapping(vecMappingWP);
   }
 
   phreeqc_time_start = MPI_Wtime();
 
-  WorkerRunWorkPackage(vecCurrWP, vecMappingWP, current_sim_time, dt);
+  if (WorkerRunWorkPackage(vecCurrWP, vecMappingWP, current_sim_time, dt) !=
+      IRM_OK) {
+    throw std::runtime_error("Phreeqc threw an error!");
+  };
 
   phreeqc_time_end = MPI_Wtime();
 
   if (dht_enabled) {
-    DHT_Results.ResultsToWP(vecCurrWP);
+    dht->resultsToWP(vecCurrWP);
   }
 
   /* send results to master */
@@ -235,7 +229,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
   if (dht_enabled) {
     /* write results to DHT */
     dht_fill_start = MPI_Wtime();
-    dht->fillDHT(local_work_package_size, DHT_Results, vecCurrWP);
+    dht->fillDHT(local_work_package_size, vecCurrWP);
     dht_fill_end = MPI_Wtime();
 
     timings.dht_get += dht_get_end - dht_get_start;
@@ -260,37 +254,14 @@ void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status,
   }
 }
 void poet::ChemistryModule::WorkerPostSim(uint32_t iteration) {
-  /* before death, submit profiling/timings to master*/
+  if (this->dht_enabled && this->dht_snaps_type == DHT_FILES_SIMEND) {
-
-  // double timings_serialized[4];
-  // timings_serialized[0] = timings.phreeqc_t;
-  // timings_serialized[1] = timings.dht_get;
-  // timings_serialized[2] = timings.dht_fill;
-  // timings_serialized[3] = timings.idle_t;
-
-  // // timings
-  // MPI_Send(timings_serialized, 4, MPI_DOUBLE, 0, 0, this->group_comm);
-
-  // // MPI_Send(&phreeqc_count, 1, MPI_INT, 0, TAG_TIMING, MPI_COMM_WORLD);
-  // // MPI_Send(&idle_t, 1, MPI_DOUBLE, 0, TAG_TIMING, MPI_COMM_WORLD);
-
-  // if (this->dht_enabled) {
-  //   // dht_perf
-  //   int dht_perf[3];
-  //   dht_perf[0] = dht->getHits();
-  //   dht_perf[1] = dht->getMisses();
-  //   dht_perf[2] = dht->getEvictions();
-  //   MPI_Send(dht_perf, 3, MPI_INT, 0, 0, this->group_comm);
-  // }
-
-  if (this->dht_enabled && this->dht_snaps_type > DHT_FILES_SIMEND) {
     WorkerWriteDHTDump(iteration);
   }
 }
 
 void poet::ChemistryModule::WorkerWriteDHTDump(uint32_t iteration) {
   std::stringstream out;
-  out << this->dht_file_out_dir << "/iter_" << std::setfill('0') << std::setw(4)
+  out << this->dht_file_out_dir << "/iter_" << std::setfill('0') << std::setw(3)
       << iteration << ".dht";
   int res = dht->tableToFile(out.str().c_str());
   if (res != DHT_SUCCESS && this->comm_rank == 2)
@@ -321,13 +292,9 @@ void poet::ChemistryModule::WorkerReadDHTDump(
 }
 
 IRM_RESULT
-poet::ChemistryModule::WorkerRunWorkPackage(std::vector<double> &vecWP,
+poet::ChemistryModule::WorkerRunWorkPackage(
-                                            std::vector<int32_t> &vecMapping,
+    std::vector<double> &vecWP, std::vector<std::uint32_t> &vecMapping,
-                                            double dSimTime, double dTimestep) {
+    double dSimTime, double dTimestep) {
-  if (this->wp_size != vecMapping.size()) {
-    return IRM_INVALIDARG;
-  }
-
   if ((this->wp_size * this->prop_count) != vecWP.size()) {
     return IRM_INVALIDARG;
   }
@@ -346,7 +313,7 @@ poet::ChemistryModule::WorkerRunWorkPackage(std::vector<double> &vecWP,
   }
 
   IRM_RESULT result;
-  this->PhreeqcRM::CreateMapping(vecMapping);
+  this->PhreeqcRM::setPOETMapping(vecMapping);
   this->setDumpedField(vecWP);
 
   this->PhreeqcRM::SetTime(dSimTime);

src/SimParams.cpp

@@ -21,6 +21,7 @@
 #include "poet/DHT_Types.hpp"
 #include <algorithm>
 #include <cassert>
+#include <cstdint>
 #include <poet/SimParams.hpp>
 
 #include <RInside.h>
@@ -86,6 +87,16 @@ poet::ChemistryParams::s_ChemistryParams(RInside &R) {
       Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$database"));
   this->input_script =
       Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$input_script"));
+  if (Rcpp::as<bool>(
+          R.parseEval("'dht_species' %in% names(mysetup$chemistry)"))) {
+    this->dht_species = Rcpp::as<std::vector<std::string>>(
+        R.parseEval("mysetup$chemistry$dht_species"));
+  }
+  if (Rcpp::as<bool>(
+          R.parseEval("'dht_signif' %in% names(mysetup$chemistry)"))) {
+    this->dht_signif = Rcpp::as<std::vector<std::uint32_t>>(
+        R.parseEval("mysetup$chemistry$dht_signif"));
+  }
 }
 
 SimParams::SimParams(int world_rank_, int world_size_) {
@@ -225,7 +236,7 @@ void SimParams::initVectorParams(RInside &R) {
 
       for (const auto &type : prop_type_R) {
         if (type == "act") {
-          this->dht_prop_type_vector.push_back(DHT_TYPE_ACT);
+          this->dht_prop_type_vector.push_back(DHT_TYPE_CHARGE);
           continue;
         }
 

util/data_evaluation/RFun_Eval.R

@@ -1,6 +1,6 @@
 ## Simple library of functions to assess and visualize the results of the coupled simulations
 
-## Time-stamp: "Last modified 2023-01-18 16:02:58 delucia"
+## Time-stamp: "Last modified 2023-04-24 16:09:55 mluebke"
 
 require(RedModRphree)
 require(Rmufits)  ## essentially for PlotCartCellData 
@@ -59,10 +59,10 @@ ReadDHT <- function(file, new_scheme = TRUE) {
 
     if (new_scheme) {
       nkeys <- dims[1] / 8
-      keys <- res[, 1:nkeys]
+      keys <- res[, 1:nkeys - 1]
 
       conv <- apply(keys, 2, ConvertDHTKey)
-      res[, 1:nkeys] <- conv
+      res[, 1:nkeys - 1] <- conv
     }
 
     return(res)