diff --git a/app/poet.cpp b/app/poet.cpp index 086948bc7..4fb2ef86e 100644 --- a/app/poet.cpp +++ b/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 ¶ms, RInside &R, chem.initializeField(diffusion.getField()); if (params.getNumParams().dht_enabled) { - chem.SetDHTEnabled(true, params.getNumParams().dht_size_per_process); - if (!params.getDHTSignifVector().empty()) { - chem.SetDHTSignifVector(params.getDHTSignifVector()); + chem.SetDHTEnabled(true, params.getNumParams().dht_size_per_process, + chem_params.dht_species); + if (!chem_params.dht_signif.empty()) { + chem.SetDHTSignifVector(chem_params.dht_signif); } if (!params.getDHTPropTypeVector().empty()) { chem.SetDHTPropTypeVector(params.getDHTPropTypeVector()); diff --git a/bench/barite/barite.R b/bench/barite/barite.R index 9b0be3971..dad3b79ea 100644 --- a/bench/barite/barite.R +++ b/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, - "O" = 55.5087, - "Charge" = -1.217E-09, + #"H" = 110.0124, + "H" = 0.00000028904, + #"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, - "O" = 55.50622, - "Charge" = -3.337E-08, + #"H" = 111.0124, + "H" = 0.0000002890408, + #"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", "", "") -prop <- names(init_cell) + +dht_species <- names(init_diffu) +#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 ) ################################################################# diff --git a/bench/dolo_diffu_inner/dolo_diffu_inner.R b/bench/dolo_diffu_inner/dolo_diffu_inner.R index 4821160e8..f243889d4 100644 --- a/bench/dolo_diffu_inner/dolo_diffu_inner.R +++ b/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, - "O" = 55.3413, + "H" = 0.000211313883539788, + "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, - "O" = 55.3413, + "H" = 0.0001540445, + "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, - "O" = 55.3413, + "H" = 0.0001610193, + "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) -prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "") -prop <- names(init_cell) +dht_species <- c("H", "O", "Charge", "C(4)", "Ca", "Cl", "Mg", "Calcite", + "Dolomite") +#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, - out_save = c(5, iterations) + store_result = TRUE ) diff --git a/bench/dolo_diffu_inner/dolo_diffu_inner_large.R b/bench/dolo_diffu_inner/dolo_diffu_inner_large.R index 40f1da3f2..69b528ce3 100644 --- a/bench/dolo_diffu_inner/dolo_diffu_inner_large.R +++ b/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, - "O" = 55.3413, + "H" = 0.000211313883539788, + "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, - "O" = 55.3413, + "H" = 0.0001540445, + "O" = 0.002148006, "Charge" = 1.90431e-16, - "C" = 0, + "C(4)" = 0, "Ca" = 0, "Cl" = 0.002, "Mg" = 0.001 ), list( - "H" = 110.683, - "O" = 55.3413, + "H" = 0.0001610193, + "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) -prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "") -prop <- names(init_cell) +dht_species <- c("H", "O", "Charge", "C(4)", "Ca", "Cl", "Mg", "Calcite", + "Dolomite") +#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) ) diff --git a/ext/phreeqcrm b/ext/phreeqcrm index 6a4c14c28..6db893859 160000 --- a/ext/phreeqcrm +++ b/ext/phreeqcrm @@ -1 +1 @@ -Subproject commit 6a4c14c2857f0ad3558af3f1e1aa7555a805e8e2 +Subproject commit 6db8938592a487a701b3bca19858ab946b42f1a0 diff --git a/include/poet/ChemistryModule.hpp b/include/poet/ChemistryModule.hpp index 07b854b1b..af63f87eb 100644 --- a/include/poet/ChemistryModule.hpp +++ b/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 &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 signif_vec); + void SetDHTSignifVector(std::vector &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 &vecWP, - std::vector &vecMapping, + std::vector &vecMapping, double dSimTime, double dTimestep); std::vector CalculateWPSizesVector(uint32_t n_cells, @@ -355,6 +359,8 @@ protected: void unshuffleField(const std::vector &in_buffer, uint32_t size_per_prop, uint32_t prop_count, uint32_t wp_count, std::vector &out_field); + std::vector + parseDHTSpeciesVec(const std::vector &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; - bool dht_enabled = false; - int dht_snaps_type = DHT_FILES_DISABLED; + uint32_t wp_size{CHEM_DEFAULT_WP_SIZE}; + bool dht_enabled{false}; + 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 base_totals{0}; + #endif double chem_t = 0.; diff --git a/include/poet/DHT_Types.hpp b/include/poet/DHT_Types.hpp index 3a57a4dc1..be0ed295b 100644 --- a/include/poet/DHT_Types.hpp +++ b/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_ diff --git a/include/poet/DHT_Wrapper.hpp b/include/poet/DHT_Wrapper.hpp index 081a23cb8..c1ebffcb1 100644 --- a/include/poet/DHT_Wrapper.hpp +++ b/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 #include #include #include @@ -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> keys; std::vector> results; std::vector needPhreeqc; - - void ResultsToMapping(std::vector &curr_mapping); - void ResultsToWP(std::vector &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 key_count Count of key entries + * @param buckets_per_process Count of buckets to allocate for each + * 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 &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 &work_package) - -> poet::DHT_ResultObject; + auto checkDHT(int length, double dt, const std::vector &work_package, + std::vector &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, - const std::vector &results); + void fillDHT(int length, const std::vector &work_package); + + void resultsToWP(std::vector &work_package); /** * @brief Dump current DHT state into file. @@ -189,6 +200,10 @@ public: void SetSignifVector(std::vector signif_vec); void SetPropTypeVector(std::vector prop_type_vec); + void setBaseTotals(const std::array &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 dht_signif_vector; - std::vector dht_prop_type_vector; + std::vector dht_prop_type_vector; + std::vector 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 base_totals{0}; }; } // namespace poet diff --git a/include/poet/SimParams.hpp b/include/poet/SimParams.hpp index 44ffb1263..0d67977c9 100644 --- a/include/poet/SimParams.hpp +++ b/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 dht_species; + std::vector dht_signif; s_ChemistryParams(RInside &R); }; diff --git a/src/ChemistryModule/ChemistryModule.cpp b/src/ChemistryModule/ChemistryModule.cpp index 497fdb517..ad58881ef 100644 --- a/src/ChemistryModule/ChemistryModule.cpp +++ b/src/ChemistryModule/ChemistryModule.cpp @@ -3,7 +3,6 @@ #include "poet/DHT_Wrapper.hpp" #include -#include #include #include #include @@ -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 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 &key_species) { + constexpr uint32_t MB_FACTOR = 1E6; + std::vector 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 poet::ChemistryModule::parseDHTSpeciesVec( + const std::vector &species_vec) const { + std::vector 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(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 signif_vec) { + std::vector &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); } diff --git a/src/ChemistryModule/DHT_Wrapper.cpp b/src/ChemistryModule/DHT_Wrapper.cpp index 320785ccd..129b45340 100644 --- a/src/ChemistryModule/DHT_Wrapper.cpp +++ b/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) { - std::int8_t exp = (std::int8_t)std::floor(std::log10(std::fabs(value))); - std::int64_t significant = value * std::pow(10, signif - exp - 1); +inline DHT_SCNotation round_key_element(double value, std::uint32_t signif) { + std::int8_t exp = + static_cast(std::floor(std::log10(std::fabs(value)))); + std::int64_t significant = + static_cast(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 &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 &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) - : key_count(key_count), data_count(data_count) { + const std::vector &key_indices, + 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 &work_package) - -> poet::DHT_ResultObject { + const std::vector &work_package, + std::vector &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; - check_data.keys.resize(length); - check_data.results.resize(length); - check_data.needPhreeqc.resize(length); + std::vector new_mapping; // 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]); - auto &data = check_data.results[i]; - auto &key_vector = check_data.keys[i]; + void *key = (void *)&(work_package[i * this->data_count]); + auto &data = dht_results.results[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, - const std::vector &results) { +void DHT_Wrapper::fillDHT(int length, const std::vector &work_package) { // 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]) { - const auto &key = dht_check_data.keys[i]; - void *data = (void *)&(results[i * this->data_count]); + if (dht_results.needPhreeqc[i]) { + const auto &key = dht_results.keys[i]; + 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 &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_Wrapper::fuzzForDHT(int var_count, void *key, std::vector 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++) { - double &curr_key = ((double *)key)[i]; + for (std::uint32_t i = 0; i < input_key_elements.size(); 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; } diff --git a/src/ChemistryModule/HashFunctions.cpp b/src/ChemistryModule/HashFunctions.cpp index da6e10138..b2739aea6 100644 --- a/src/ChemistryModule/HashFunctions.cpp +++ b/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 -#include - #if defined(_MSC_VER) #define BIG_CONSTANT(x) (x) diff --git a/src/ChemistryModule/WorkerFunctions.cpp b/src/ChemistryModule/WorkerFunctions.cpp index c4f940767..a9ed3d3bf 100644 --- a/src/ChemistryModule/WorkerFunctions.cpp +++ b/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 +#include #include #include #include @@ -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 name_dummy; + + SetDHTEnabled(enable, size_mb, name_dummy); break; } case CHEM_DHT_SIGNIF_VEC: { - std::vector input_vec(this->prop_count); - ChemBCast(input_vec.data(), this->prop_count, MPI_UINT32_T); + std::vector input_vec; 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 vecCurrWP( - // mpi_buffer, - // mpi_buffer + (local_work_package_size * this->prop_names.size())); vecCurrWP.resize(n_cells_times_props); - std::vector vecMappingWP(this->wp_size); + std::vector vecMappingWP(local_work_package_size); - DHT_ResultObject DHT_Results; - - for (uint32_t i = 0; i < local_work_package_size; i++) { - vecMappingWP[i] = i; - } - - if (local_work_package_size != this->wp_size) { - // std::vector 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; - } + { + std::uint32_t i = 0; + std::generate(vecMappingWP.begin(), vecMappingWP.end(), + [&] { return i++; }); } 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*/ - - // 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) { + 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 &vecWP, - std::vector &vecMapping, - double dSimTime, double dTimestep) { - if (this->wp_size != vecMapping.size()) { - return IRM_INVALIDARG; - } - +poet::ChemistryModule::WorkerRunWorkPackage( + std::vector &vecWP, std::vector &vecMapping, + double dSimTime, double dTimestep) { if ((this->wp_size * this->prop_count) != vecWP.size()) { return IRM_INVALIDARG; } @@ -346,7 +313,7 @@ poet::ChemistryModule::WorkerRunWorkPackage(std::vector &vecWP, } IRM_RESULT result; - this->PhreeqcRM::CreateMapping(vecMapping); + this->PhreeqcRM::setPOETMapping(vecMapping); this->setDumpedField(vecWP); this->PhreeqcRM::SetTime(dSimTime); diff --git a/src/SimParams.cpp b/src/SimParams.cpp index 7a536134c..8e25868d1 100644 --- a/src/SimParams.cpp +++ b/src/SimParams.cpp @@ -21,6 +21,7 @@ #include "poet/DHT_Types.hpp" #include #include +#include #include #include @@ -86,6 +87,16 @@ poet::ChemistryParams::s_ChemistryParams(RInside &R) { Rcpp::as(R.parseEval("mysetup$chemistry$database")); this->input_script = Rcpp::as(R.parseEval("mysetup$chemistry$input_script")); + if (Rcpp::as( + R.parseEval("'dht_species' %in% names(mysetup$chemistry)"))) { + this->dht_species = Rcpp::as>( + R.parseEval("mysetup$chemistry$dht_species")); + } + if (Rcpp::as( + R.parseEval("'dht_signif' %in% names(mysetup$chemistry)"))) { + this->dht_signif = Rcpp::as>( + 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; } diff --git a/util/data_evaluation/RFun_Eval.R b/util/data_evaluation/RFun_Eval.R index 4175d9395..ed0068689 100644 --- a/util/data_evaluation/RFun_Eval.R +++ b/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)