From 5493cf88a9b5cf626ec79327922e4558a97cd19e Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Thu, 12 Dec 2024 15:24:57 +0100 Subject: [PATCH 01/21] Added qs2 as new default format --- R_lib/kin_r_library.R | 137 ++++++++++++++++++++++-------------------- bench/CMakeLists.txt | 2 +- src/initializer.cpp | 16 ++++- src/poet.cpp | 13 ++-- src/poet.hpp.in | 4 +- 5 files changed, 97 insertions(+), 75 deletions(-) diff --git a/R_lib/kin_r_library.R b/R_lib/kin_r_library.R index cb5552a9b..50497b958 100644 --- a/R_lib/kin_r_library.R +++ b/R_lib/kin_r_library.R @@ -109,69 +109,6 @@ msgm <- function(...) { } -## Function called by master R process to store on disk all relevant -## parameters for the simulation -StoreSetup <- function(setup, filesim, out_dir) { - to_store <- vector(mode = "list", length = 4) - ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT") - names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline") - - ## read the setup R file, which is sourced in kin.cpp - tmpbuff <- file(filesim, "r") - setupfile <- readLines(tmpbuff) - close.connection(tmpbuff) - - to_store$Sim <- setupfile - - ## to_store$Flow <- list( - ## snapshots = setup$snapshots, - ## gridfile = setup$gridfile, - ## phase = setup$phase, - ## density = setup$density, - ## dt_differ = setup$dt_differ, - ## prolong = setup$prolong, - ## maxiter = setup$maxiter, - ## saved_iter = setup$iter_output, - ## out_save = setup$out_save ) - - to_store$Transport <- setup$diffusion - - ## to_store$Chemistry <- list( - ## nprocs = n_procs, - ## wp_size = work_package_size, - ## base = setup$base, - ## first = setup$first, - ## init = setup$initsim, - ## db = db, - ## kin = setup$kin, - ## ann = setup$ann) - - if (dht_enabled) { - to_store$DHT <- list( - enabled = dht_enabled, - log = dht_log - ## signif = dht_final_signif, - ## proptype = dht_final_proptype - ) - } else { - to_store$DHT <- FALSE - } - - if (dht_enabled) { - to_store$DHT <- list( - enabled = dht_enabled, - log = dht_log - # signif = dht_final_signif, - # proptype = dht_final_proptype - ) - } else { - to_store$DHT <- FALSE - } - - saveRDS(to_store, file = paste0(fileout, "/setup.rds")) - msgm("initialization stored in ", paste0(fileout, "/setup.rds")) -} - GetWorkPackageSizesVector <- function(n_packages, package_size, len) { ids <- rep(1:n_packages, times = package_size, each = 1)[1:len] return(as.integer(table(ids))) @@ -179,7 +116,7 @@ GetWorkPackageSizesVector <- function(n_packages, package_size, len) { ## Handler to read R objs from binary files using either builtin -## readRDS() or qs::qread() based on file extension +## readRDS(), qs::qread() or qs2::qs_read() based on file extension ReadRObj <- function(path) { ## code borrowed from tools::file_ext() pos <- regexpr("\\.([[:alnum:]]+)$", path) @@ -187,7 +124,8 @@ ReadRObj <- function(path) { switch(extension, rds = readRDS(path), - qs = qs::qread(path) + qs = qs::qread(path), + qs2 = qs2::qs_read(path) ) } @@ -201,6 +139,73 @@ SaveRObj <- function(x, path) { switch(extension, rds = saveRDS(object = x, file = path), - qs = qs::qsave(x = x, file = path) + qs = qs::qsave(x = x, file = path), + qs2 = qs2::qs_save(object = x, file = path) ) } + + +######## Old relic code + +## ## Function called by master R process to store on disk all relevant +## ## parameters for the simulation +## StoreSetup <- function(setup, filesim, out_dir) { +## to_store <- vector(mode = "list", length = 4) +## ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT") +## names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline") + +## ## read the setup R file, which is sourced in kin.cpp +## tmpbuff <- file(filesim, "r") +## setupfile <- readLines(tmpbuff) +## close.connection(tmpbuff) + +## to_store$Sim <- setupfile + +## ## to_store$Flow <- list( +## ## snapshots = setup$snapshots, +## ## gridfile = setup$gridfile, +## ## phase = setup$phase, +## ## density = setup$density, +## ## dt_differ = setup$dt_differ, +## ## prolong = setup$prolong, +## ## maxiter = setup$maxiter, +## ## saved_iter = setup$iter_output, +## ## out_save = setup$out_save ) + +## to_store$Transport <- setup$diffusion + +## ## to_store$Chemistry <- list( +## ## nprocs = n_procs, +## ## wp_size = work_package_size, +## ## base = setup$base, +## ## first = setup$first, +## ## init = setup$initsim, +## ## db = db, +## ## kin = setup$kin, +## ## ann = setup$ann) + +## if (dht_enabled) { +## to_store$DHT <- list( +## enabled = dht_enabled, +## log = dht_log +## ## signif = dht_final_signif, +## ## proptype = dht_final_proptype +## ) +## } else { +## to_store$DHT <- FALSE +## } + +## if (dht_enabled) { +## to_store$DHT <- list( +## enabled = dht_enabled, +## log = dht_log +## # signif = dht_final_signif, +## # proptype = dht_final_proptype +## ) +## } else { +## to_store$DHT <- FALSE +## } + +## saveRDS(to_store, file = paste0(fileout, "/setup.rds")) +## msgm("initialization stored in ", paste0(fileout, "/setup.rds")) +## } diff --git a/bench/CMakeLists.txt b/bench/CMakeLists.txt index 01dc43caf..794774c4f 100644 --- a/bench/CMakeLists.txt +++ b/bench/CMakeLists.txt @@ -7,7 +7,7 @@ function(ADD_BENCH_TARGET TARGET POET_BENCH_LIST RT_FILES OUT_PATH) foreach(BENCH_FILE ${${POET_BENCH_LIST}}) get_filename_component(BENCH_NAME ${BENCH_FILE} NAME_WE) set(OUT_FILE ${CMAKE_CURRENT_BINARY_DIR}/${BENCH_NAME}) - set(OUT_FILE_EXT ${OUT_FILE}.qs) + set(OUT_FILE_EXT ${OUT_FILE}.qs2) add_custom_command( OUTPUT ${OUT_FILE_EXT} diff --git a/src/initializer.cpp b/src/initializer.cpp index d2e663931..d613b83be 100644 --- a/src/initializer.cpp +++ b/src/initializer.cpp @@ -35,7 +35,11 @@ int main(int argc, char **argv) { ->default_val(false); bool asRDS; - app.add_flag("-r, --rds", asRDS, "Save output as .rds file instead of .qs") + app.add_flag("-r, --rds", asRDS, "Save output as .rds") + ->default_val(false); + + bool asQS; + app.add_flag("-q, --qs", asQS, "Save output as .qs") ->default_val(false); CLI11_PARSE(app, argc, argv); @@ -69,8 +73,14 @@ int main(int argc, char **argv) { } // append the correct file extension - output_file += asRDS ? ".rds" : ".qs"; - + if (asRDS) { + output_file += ".rds"; + } else if (asQS) { + output_file += ".qs"; + } else { + output_file += ".qs2"; + } + // set working directory to the directory of the input script if (setwd) { const std::string dir_path = Rcpp::as( diff --git a/src/poet.cpp b/src/poet.cpp index 9fbf94c18..d4f27525f 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -57,7 +57,7 @@ static std::unique_ptr global_rt_setup; // before the R runtime is initialized static poet::DEFunc master_init_R; static poet::DEFunc master_iteration_end_R; -static poet::DEFunc store_setup_R; +// MDL: unused -> static poet::DEFunc store_setup_R; static poet::DEFunc ReadRObj_R; static poet::DEFunc SaveRObj_R; static poet::DEFunc source_R; @@ -66,7 +66,7 @@ static void init_global_functions(RInside &R) { R.parseEval(kin_r_library); master_init_R = DEFunc("master_init"); master_iteration_end_R = DEFunc("master_iteration_end"); - store_setup_R = DEFunc("StoreSetup"); + // MDL: unused -> store_setup_R = DEFunc("StoreSetup"); source_R = DEFunc("source"); ReadRObj_R = DEFunc("ReadRObj"); SaveRObj_R = DEFunc("SaveRObj"); @@ -146,8 +146,11 @@ int parseInitValues(int argc, char **argv, RuntimeParameters ¶ms) { "Enable AI surrogate for chemistry module"); app.add_flag("--rds", params.as_rds, - "Save output as .rds file instead of .qs"); + "Save output as .rds file instead of default .qs2"); + app.add_flag("--qs", params.as_qs, + "Save output as .qs file instead of default .qs2"); + std::string init_file; std::string runtime_file; @@ -174,7 +177,9 @@ int parseInitValues(int argc, char **argv, RuntimeParameters ¶ms) { } // set the output extension - params.out_ext = params.as_rds ? "rds" : "qs"; + params.out_ext = "qs2"; + if (params.as_rds) params.out_ext = "rds"; + if (params.as_qs) params.out_ext = "qs"; if (MY_RANK == 0) { // MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result)); diff --git a/src/poet.hpp.in b/src/poet.hpp.in index 0e2409f87..c48a0f3df 100644 --- a/src/poet.hpp.in +++ b/src/poet.hpp.in @@ -45,8 +45,10 @@ struct RuntimeParameters { Rcpp::List init_params; + // MDL added to accomodate for qs::qsave/qread bool as_rds = false; - std::string out_ext; // MDL added to accomodate for qs::qsave/qread + bool as_qs = false; + std::string out_ext; bool print_progress = false; From 54d6b6bbf298bf6f44f02953c54d683c622bd057 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Thu, 12 Dec 2024 16:20:09 +0100 Subject: [PATCH 02/21] Less and more informative stdout messages --- R_lib/kin_r_library.R | 4 ++-- src/Transport/DiffusionModule.cpp | 2 +- src/poet.cpp | 30 +++++++++++++++++------------- 3 files changed, 20 insertions(+), 16 deletions(-) diff --git a/R_lib/kin_r_library.R b/R_lib/kin_r_library.R index 50497b958..97fd49e16 100644 --- a/R_lib/kin_r_library.R +++ b/R_lib/kin_r_library.R @@ -94,7 +94,7 @@ master_iteration_end <- function(setup, state_T, state_C) { ## Add last time step to simulation time setup$simulation_time <- setup$simulation_time + setup$timesteps[iter] - msgm("done iteration", iter, "/", length(setup$timesteps)) + ## msgm("done iteration", iter, "/", length(setup$timesteps)) setup$iter <- setup$iter + 1 return(setup) } @@ -132,7 +132,7 @@ ReadRObj <- function(path) { ## Handler to store R objs to binary files using either builtin ## saveRDS() or qs::qsave() based on file extension SaveRObj <- function(x, path) { - msgm("Storing to", path) + ## msgm("Storing to", path) ## code borrowed from tools::file_ext() pos <- regexpr("\\.([[:alnum:]]+)$", path) extension <- ifelse(pos > -1L, substring(path, pos + 1L), "") diff --git a/src/Transport/DiffusionModule.cpp b/src/Transport/DiffusionModule.cpp index 2ea80564a..754d135c9 100644 --- a/src/Transport/DiffusionModule.cpp +++ b/src/Transport/DiffusionModule.cpp @@ -70,7 +70,7 @@ VecToMatrix(const std::vector &vec, std::uint32_t n_rows, // static constexpr double ZERO_MULTIPLICATOR = 10E-14; void DiffusionModule::simulate(double requested_dt) { - MSG("Starting diffusion ..."); + // MSG("Starting diffusion ..."); const auto start_diffusion_t = std::chrono::high_resolution_clock::now(); const auto &n_rows = this->param_list.n_rows; diff --git a/src/poet.cpp b/src/poet.cpp index d4f27525f..93a168496 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -39,13 +39,14 @@ #include #include #include +#include #include #include #include -using namespace std; +// using namespace std; using namespace poet; using namespace Rcpp; @@ -292,18 +293,19 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, const double &dt = params.timesteps[iter - 1]; - // cout << "CPP: Next time step is " << dt << "[s]" << endl; - MSG("Next time step is " + std::to_string(dt) + " [s]"); - + std::cout << std::endl; /* displaying iteration number, with C++ and R iterator */ - MSG("Going through iteration " + std::to_string(iter)); + MSG("Going through iteration " + std::to_string(iter) + "/" + + std::to_string(maxiter)); + + MSG("Current time step is " + std::format("{:.2f}", dt)); /* run transport */ diffusion.simulate(dt); chem.getField().update(diffusion.getField()); - MSG("Chemistry step"); + // MSG("Chemistry start"); if (params.use_ai_surrogate) { double ai_start_t = MPI_Wtime(); // Save current values from the tug field as predictor for the ai step @@ -319,16 +321,16 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, R.parseEval("predictors_scaled <- preprocess(predictors)"); // Predict - MSG("AI Predict"); + MSG("AI Prediction"); R.parseEval( "aipreds_scaled <- prediction_step(model, predictors_scaled)"); // Apply postprocessing - MSG("AI Postprocesing"); + MSG("AI Postprocessing"); R.parseEval("aipreds <- postprocess(aipreds_scaled)"); // Validate prediction and write valid predictions to chem field - MSG("AI Validate"); + MSG("AI Validation"); R.parseEval( "validity_vector <- validate_predictions(predictors, aipreds)"); @@ -338,8 +340,8 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, MSG("AI TempField"); std::vector> RTempField = R.parseEval("set_valid_predictions(predictors,\ - aipreds,\ - validity_vector)"); + aipreds,\ + validity_vector)"); MSG("AI Set Field"); Field predictions_field = @@ -390,9 +392,11 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, MSG("End of *coupling* iteration " + std::to_string(iter) + "/" + std::to_string(maxiter)); - MSG(); + // MSG(); } // END SIMULATION LOOP + std::cout << std::endl; + Rcpp::List chem_profiling; chem_profiling["simtime"] = chem.GetChemistryTime(); chem_profiling["loop"] = chem.GetMasterLoopTime(); @@ -588,7 +592,7 @@ int main(int argc, char *argv[]) { R["setup"] = *global_rt_setup; R["setup$out_ext"] = run_params.out_ext; - string r_vis_code; + std::string r_vis_code; r_vis_code = "SaveRObj(x = profiling, path = paste0(out_dir, " "'/timings.', setup$out_ext));"; R.parseEval(r_vis_code); From 3d0a9041b5ed682a1bbc53aa271ee69d2344d6e4 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Thu, 12 Dec 2024 16:42:44 +0100 Subject: [PATCH 03/21] Update README: qs2 as default output format, gfz.de everywhere --- README.md | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index d71766ecd..cc4e643cf 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ original MPI-based Distributed Hash Table. ## Parsed code documentiation A parsed version of POET's documentation can be found at [Gitlab -pages](https://naaice.git-pages.gfz-potsdam.de/poet). +pages](https://naaice.git-pages.gfz.de/poet). ## External Libraries @@ -18,8 +18,8 @@ The following external libraries are shipped with POET: - **CLI11** - - **IPhreeqc** with patches from GFZ/UP - - - -- **tug** - + +- **tug** - ## Installation @@ -41,7 +41,7 @@ installed: - [Rcpp](https://cran.r-project.org/web/packages/Rcpp/index.html) - [RInside](https://cran.r-project.org/web/packages/RInside/index.html) - [qs](https://cran.r-project.org/web/packages/qs/index.html) - +- [qs2](https://cran.r-project.org/web/packages/qs2/index.html) This can be simply achieved by issuing the following commands: ```sh @@ -49,7 +49,7 @@ This can be simply achieved by issuing the following commands: $ R # install R dependencies (case sensitive!) -> install.packages(c("Rcpp", "RInside","qs")) +> install.packages(c("Rcpp", "RInside","qs","qs2")) > q(save="no") ``` @@ -59,7 +59,7 @@ POET can be anonimously cloned from this repo over https. Make sure to also download the submodules: ```sh -git clone --recurse-submodules https://git.gfz-potsdam.de/naaice/poet.git +git clone --recurse-submodules https://git.gfz.de/naaice/poet.git ``` The `--recurse-submodules` option is a shorthand for: ```sh @@ -110,7 +110,7 @@ follows: $ R # install R dependencies -> install.packages(c("Rcpp", "RInside","qs")) +> install.packages(c("Rcpp", "RInside","qs","qs2")) > q(save="no") # cd into POET project root @@ -138,17 +138,17 @@ poet └── share └── poet ├── barite - │   ├── barite_200.rds + │   ├── barite_200.qs2 │   ├── barite_200_rt.R - │   ├── barite_het.rds + │   ├── barite_het.qs2 │   └── barite_het_rt.R ├── dolo - │   ├── dolo_inner_large.rds + │   ├── dolo_inner_large.qs2 │   ├── dolo_inner_large_rt.R - │   ├── dolo_interp.rds + │   ├── dolo_interp.qs2 │   └── dolo_interp_rt.R └── surfex - ├── PoetEGU_surfex_500.rds + ├── PoetEGU_surfex_500.qs2 └── PoetEGU_surfex_500_rt.R ``` @@ -182,7 +182,8 @@ The following parameters can be set: | **-P, --progress** | | show progress bar | | **--ai-surrogate** | | activates the AI surrogate chemistry model (defaults to _OFF_) | | **--dht** | | enabling DHT usage (defaults to _OFF_) | -| **--qs** | | store results using qs::qsave() (.qs extension) instead of default RDS (.rds) | +| **--qs** | | store results using qs::qsave() (.qs extension) instead of default qs2 (.qs2) | +| **--rds** | | store results using saveRDS() (.rds extension) instead of default qs2 (.qs2) | | **--dht-strategy=** | _0-1_ | change DHT strategy. **NOT IMPLEMENTED YET** (Defaults to _0_) | | **--dht-size=** | _1-n_ | size of DHT per process involved in megabyte (defaults to _1000 MByte_) | | **--dht-snaps=** | _0-2_ | disable or enable storage of DHT snapshots | @@ -284,7 +285,7 @@ produce any valid predictions. In order to provide a model to POET, you need to setup a R script which can then be used by `poet_init` to generate the simulation input. Which parameters are required can be found in the -[Wiki](https://git.gfz-potsdam.de/naaice/poet/-/wikis/Initialization). +[Wiki](https://git.gfz.de/naaice/poet/-/wikis/Initialization). We try to keep the document up-to-date. However, if you encounter missing information or need help, please get in touch with us via the issue tracker or E-Mail. @@ -298,7 +299,7 @@ issue tracker or E-Mail. where: - **output** - name of the output file (defaults to the input file - name with the extension `.rds`) + name with the extension `.qs2`) - **setwd** - set the working directory to the directory of the input file (e.g. to allow relative paths in the input script). However, the output file will be stored in the directory from which From c4780fe9af4fe277c65793713088b4873143bcd2 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Thu, 12 Dec 2024 16:56:05 +0100 Subject: [PATCH 04/21] Update references to .qs2 in README --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index cc4e643cf..9e31ff6d4 100644 --- a/README.md +++ b/README.md @@ -214,7 +214,7 @@ installed POET-dir `/bin` and run: ```sh cp ../share/poet/barite/barite_het* . -mpirun -n 4 ./poet barite_het_rt.R barite_het.rds output +mpirun -n 4 ./poet barite_het_rt.R barite_het.qs2 output ``` After a finished simulation all data generated by POET will be found @@ -227,7 +227,7 @@ DHT snapshot shall be produced. This is done by appending the `--dht-snaps=` option. The resulting call would look like this: ```sh -mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.rds output +mpirun -n 4 ./poet --dht --dht-snaps=2 barite_het_rt.R barite_het.qs2 output ``` ### Example: Preparing Environment and Running with AI surrogate @@ -274,7 +274,7 @@ cp /bench/barite/{barite_50ai*,db_barite.dat,barite.pqi} . ./poet_init barite_50ai.R # run POET with AI surrogate and GPU utilization -srun --gres=gpu -N 1 -n 12 ./poet --ai-surrogate barite_50ai_rt.R barite_50ai.rds output +srun --gres=gpu -N 1 -n 12 ./poet --ai-surrogate barite_50ai_rt.R barite_50ai.qs2 output ``` Keep in mind that the AI surrogate is currently not stable or might also not From 3a164ed71689218bcb58b55a1a18e47e012977e8 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Thu, 12 Dec 2024 18:05:50 +0100 Subject: [PATCH 05/21] reverting since gcc < 13 does not support it --- src/poet.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/poet.cpp b/src/poet.cpp index 93a168496..3da740096 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -39,7 +39,6 @@ #include #include #include -#include #include @@ -294,11 +293,12 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, const double &dt = params.timesteps[iter - 1]; std::cout << std::endl; + /* displaying iteration number, with C++ and R iterator */ MSG("Going through iteration " + std::to_string(iter) + "/" + std::to_string(maxiter)); - MSG("Current time step is " + std::format("{:.2f}", dt)); + MSG("Current time step is " + std::to_string(dt)); /* run transport */ diffusion.simulate(dt); From 07ac2bfd8dbdb3d3aa75f22148c015c547a5eced Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 12 Dec 2024 20:16:02 +0100 Subject: [PATCH 06/21] [wip] fix: add base_totals to SurrogateSetup --- src/Chemistry/ChemistryModule.cpp | 12 ++++++---- src/Chemistry/ChemistryModule.hpp | 7 ++++-- src/poet.cpp | 39 ++++++++++++++++++++++++++----- 3 files changed, 45 insertions(+), 13 deletions(-) diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index dc2d430f1..a6ab03c36 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -309,9 +309,10 @@ void poet::ChemistryModule::initializeInterp( map_copy = this->dht->getKeySpecies(); for (auto i = 0; i < map_copy.size(); i++) { const std::uint32_t signif = - static_cast(map_copy[i]) - (map_copy[i] > InterpolationModule::COARSE_DIFF - ? InterpolationModule::COARSE_DIFF - : 0); + static_cast(map_copy[i]) - + (map_copy[i] > InterpolationModule::COARSE_DIFF + ? InterpolationModule::COARSE_DIFF + : 0); map_copy[i] = signif; } } @@ -368,7 +369,8 @@ void poet::ChemistryModule::unshuffleField(const std::vector &in_buffer, } } } - -void poet::ChemistryModule::set_ai_surrogate_validity_vector(std::vector r_vector) { + +void poet::ChemistryModule::set_ai_surrogate_validity_vector( + std::vector r_vector) { this->ai_surrogate_validity_vector = r_vector; } diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp index d2a9fca0c..f0c164754 100644 --- a/src/Chemistry/ChemistryModule.hpp +++ b/src/Chemistry/ChemistryModule.hpp @@ -76,6 +76,7 @@ public: struct SurrogateSetup { std::vector prop_names; + std::array base_totals; bool dht_enabled; std::uint32_t dht_size_mb; @@ -96,6 +97,8 @@ public: this->interp_enabled = setup.interp_enabled; this->ai_surrogate_enabled = setup.ai_surrogate_enabled; + this->base_totals = setup.base_totals; + if (this->dht_enabled || this->interp_enabled) { this->initializeDHT(setup.dht_size_mb, this->params.dht_species); } @@ -223,8 +226,8 @@ public: }; /** - * **Master only** Set the ai surrogate validity vector from R - */ + * **Master only** Set the ai surrogate validity vector from R + */ void set_ai_surrogate_validity_vector(std::vector r_vector); std::vector GetWorkerInterpolationCalls() const; diff --git a/src/poet.cpp b/src/poet.cpp index 3da740096..df51ce31c 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -34,6 +34,8 @@ #include #include #include +#include +#include #include #include #include @@ -150,7 +152,7 @@ int parseInitValues(int argc, char **argv, RuntimeParameters ¶ms) { app.add_flag("--qs", params.as_qs, "Save output as .qs file instead of default .qs2"); - + std::string init_file; std::string runtime_file; @@ -178,8 +180,10 @@ int parseInitValues(int argc, char **argv, RuntimeParameters ¶ms) { // set the output extension params.out_ext = "qs2"; - if (params.as_rds) params.out_ext = "rds"; - if (params.as_qs) params.out_ext = "qs"; + if (params.as_rds) + params.out_ext = "rds"; + if (params.as_qs) + params.out_ext = "qs"; if (MY_RANK == 0) { // MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result)); @@ -293,7 +297,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, const double &dt = params.timesteps[iter - 1]; std::cout << std::endl; - + /* displaying iteration number, with C++ and R iterator */ MSG("Going through iteration " + std::to_string(iter) + "/" + std::to_string(maxiter)); @@ -396,7 +400,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters ¶ms, } // END SIMULATION LOOP std::cout << std::endl; - + Rcpp::List chem_profiling; chem_profiling["simtime"] = chem.GetChemistryTime(); chem_profiling["loop"] = chem.GetMasterLoopTime(); @@ -483,6 +487,29 @@ std::vector getSpeciesNames(const Field &&field, int root, return species_names_out; } +std::array getBaseTotals(Field &&field, int root, MPI_Comm comm) { + std::array base_totals; + + int rank; + MPI_Comm_rank(comm, &rank); + + const bool is_master = root == rank; + + if (is_master) { + const auto h_col = field["H"]; + const auto o_col = field["O"]; + + base_totals[0] = *std::min_element(h_col.begin(), h_col.end()); + base_totals[1] = *std::min_element(o_col.begin(), o_col.end()); + MPI_Bcast(base_totals.data(), 2, MPI_DOUBLE, root, MPI_COMM_WORLD); + return base_totals; + } + + MPI_Bcast(base_totals.data(), 2, MPI_DOUBLE, root, comm); + + return base_totals; +} + int main(int argc, char *argv[]) { int world_size; @@ -529,8 +556,8 @@ int main(int argc, char *argv[]) { init_list.getChemistryInit(), MPI_COMM_WORLD); const ChemistryModule::SurrogateSetup surr_setup = { - getSpeciesNames(init_list.getInitialGrid(), 0, MPI_COMM_WORLD), + getBaseTotals(init_list.getInitialGrid(), 0, MPI_COMM_WORLD), run_params.use_dht, run_params.dht_size, run_params.use_interp, From 1b4691c1720e24e14694c824ca01d2cf99e1fb01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 13 Dec 2024 08:21:25 +0100 Subject: [PATCH 07/21] feat: implement caching for interpolation calculations and add NaN handling --- src/Chemistry/ChemistryModule.cpp | 1 + .../SurrogateModels/Interpolation.hpp | 2 ++ .../SurrogateModels/InterpolationModule.cpp | 18 +++++++++++++++++- src/Chemistry/SurrogateModels/LookupKey.hpp | 12 ++++++++++++ src/Chemistry/SurrogateModels/Rounding.hpp | 5 +++++ 5 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index a6ab03c36..ed0b17d9a 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -8,6 +8,7 @@ #include #include +#include #include #include #include diff --git a/src/Chemistry/SurrogateModels/Interpolation.hpp b/src/Chemistry/SurrogateModels/Interpolation.hpp index b6b9a69f5..efa79e4f0 100644 --- a/src/Chemistry/SurrogateModels/Interpolation.hpp +++ b/src/Chemistry/SurrogateModels/Interpolation.hpp @@ -261,6 +261,8 @@ private: const InitialList::ChemistryHookFunctions &hooks; const std::vector &out_names; const std::vector dht_names; + + std::unordered_map> to_calc_cache; }; } // namespace poet diff --git a/src/Chemistry/SurrogateModels/InterpolationModule.cpp b/src/Chemistry/SurrogateModels/InterpolationModule.cpp index a2871bb8a..4185e0292 100644 --- a/src/Chemistry/SurrogateModels/InterpolationModule.cpp +++ b/src/Chemistry/SurrogateModels/InterpolationModule.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -116,10 +117,25 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) { this->pht->incrementReadCounter(roundKey(rounded_key)); #endif + const int cell_id = static_cast(work_package.input[wp_i][0]); + + if (!to_calc_cache.contains(cell_id)) { + std::vector to_calc = dht_instance.getKeyElements(); + std::vector keep_indices; + + for (std::size_t i = 0; i < to_calc.size(); i++) { + if (!std::isnan(work_package.input[wp_i][to_calc[i]])) { + keep_indices.push_back(to_calc[i]); + } + } + + to_calc_cache[cell_id] = keep_indices; + } + double start_fc = MPI_Wtime(); work_package.output[wp_i] = - f_interpolate(dht_instance.getKeyElements(), work_package.input[wp_i], + f_interpolate(to_calc_cache[cell_id], work_package.input[wp_i], pht_result.in_values, pht_result.out_values); if (hooks.interp_post.isValid()) { diff --git a/src/Chemistry/SurrogateModels/LookupKey.hpp b/src/Chemistry/SurrogateModels/LookupKey.hpp index e6dc7e697..6e7b0681c 100644 --- a/src/Chemistry/SurrogateModels/LookupKey.hpp +++ b/src/Chemistry/SurrogateModels/LookupKey.hpp @@ -10,9 +10,21 @@ namespace poet { +constexpr std::int8_t SC_NOTATION_EXPONENT_MASK = -128; +constexpr std::int64_t SC_NOTATION_SIGNIFICANT_MASK = 0xFFFFFFFFFFFF; + struct Lookup_SC_notation { std::int8_t exp : 8; std::int64_t significant : 56; + + constexpr static Lookup_SC_notation nan() { + return {SC_NOTATION_EXPONENT_MASK, SC_NOTATION_SIGNIFICANT_MASK}; + } + + constexpr bool isnan() { + return !!(exp == SC_NOTATION_EXPONENT_MASK && + significant == SC_NOTATION_SIGNIFICANT_MASK); + } }; union Lookup_Keyelement { diff --git a/src/Chemistry/SurrogateModels/Rounding.hpp b/src/Chemistry/SurrogateModels/Rounding.hpp index 3f659290e..688ac4707 100644 --- a/src/Chemistry/SurrogateModels/Rounding.hpp +++ b/src/Chemistry/SurrogateModels/Rounding.hpp @@ -20,6 +20,11 @@ class DHT_Rounder { public: Lookup_Keyelement round(const double &value, std::uint32_t signif, bool is_ho) { + + if (std::isnan(value)) { + return {.sc_notation = Lookup_SC_notation::nan()}; + } + std::int8_t exp = static_cast(std::floor(std::log10(std::fabs(value)))); From 1a680d8ca6cb6af1af11302df99e6858723ce853 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 13 Dec 2024 13:11:03 +0100 Subject: [PATCH 08/21] fix: remove double initialization of dht_enabled in ChemistryModule --- src/Chemistry/ChemistryModule.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index ed0b17d9a..cbfe0745c 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -188,8 +188,6 @@ void poet::ChemistryModule::initializeDHT( uint32_t size_mb, const NamedVector &key_species) { constexpr uint32_t MB_FACTOR = 1E6; - this->dht_enabled = true; - MPI_Comm dht_comm; if (this->is_master) { From 5f56ce9e3f594b7e54784bbdb369f9ec901e22c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 13 Dec 2024 13:12:16 +0100 Subject: [PATCH 09/21] fix: use const reference for to_calc in InterpolationModule to improve performance --- src/Chemistry/SurrogateModels/InterpolationModule.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Chemistry/SurrogateModels/InterpolationModule.cpp b/src/Chemistry/SurrogateModels/InterpolationModule.cpp index 4185e0292..8182efe10 100644 --- a/src/Chemistry/SurrogateModels/InterpolationModule.cpp +++ b/src/Chemistry/SurrogateModels/InterpolationModule.cpp @@ -120,7 +120,7 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) { const int cell_id = static_cast(work_package.input[wp_i][0]); if (!to_calc_cache.contains(cell_id)) { - std::vector to_calc = dht_instance.getKeyElements(); + const std::vector &to_calc = dht_instance.getKeyElements(); std::vector keep_indices; for (std::size_t i = 0; i < to_calc.size(); i++) { From cd53f43bd2627c3cbc18395319013f4c0ebf3067 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 13 Dec 2024 14:33:02 +0100 Subject: [PATCH 10/21] feat: add has_het_ids parameter to DHT initialization and related functions --- src/Chemistry/ChemistryModule.cpp | 5 ++-- src/Chemistry/ChemistryModule.hpp | 7 +++-- src/Chemistry/SurrogateModels/DHT_Wrapper.cpp | 15 ++++++++--- src/Chemistry/SurrogateModels/DHT_Wrapper.hpp | 3 ++- src/poet.cpp | 27 +++++++++++++++++++ 5 files changed, 48 insertions(+), 9 deletions(-) diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index cbfe0745c..32c14be90 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -185,7 +185,8 @@ poet::ChemistryModule::~ChemistryModule() { } void poet::ChemistryModule::initializeDHT( - uint32_t size_mb, const NamedVector &key_species) { + uint32_t size_mb, const NamedVector &key_species, + bool has_het_ids) { constexpr uint32_t MB_FACTOR = 1E6; MPI_Comm dht_comm; @@ -217,7 +218,7 @@ void poet::ChemistryModule::initializeDHT( this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices, this->prop_names, params.hooks, - this->prop_count, interp_enabled); + this->prop_count, interp_enabled, has_het_ids); this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1)); } } diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp index f0c164754..697b50744 100644 --- a/src/Chemistry/ChemistryModule.hpp +++ b/src/Chemistry/ChemistryModule.hpp @@ -77,6 +77,7 @@ public: struct SurrogateSetup { std::vector prop_names; std::array base_totals; + bool has_het_ids; bool dht_enabled; std::uint32_t dht_size_mb; @@ -100,7 +101,8 @@ public: this->base_totals = setup.base_totals; if (this->dht_enabled || this->interp_enabled) { - this->initializeDHT(setup.dht_size_mb, this->params.dht_species); + this->initializeDHT(setup.dht_size_mb, this->params.dht_species, + setup.has_het_ids); } if (this->interp_enabled) { @@ -243,7 +245,8 @@ public: protected: void initializeDHT(uint32_t size_mb, - const NamedVector &key_species); + const NamedVector &key_species, + bool has_het_ids); void setDHTSnapshots(int type, const std::string &out_dir); void setDHTReadFile(const std::string &input_file); diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp index 095f83cfa..0d9d2d2f6 100644 --- a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp +++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp @@ -43,11 +43,12 @@ DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size, const std::vector &key_indices, const std::vector &_output_names, const InitialList::ChemistryHookFunctions &_hooks, - uint32_t data_count, bool _with_interp) + uint32_t data_count, bool _with_interp, + bool _has_het_ids) : key_count(key_indices.size()), data_count(data_count), input_key_elements(key_indices), communicator(dht_comm), key_species(key_species), output_names(_output_names), hooks(_hooks), - with_interp(_with_interp) { + with_interp(_with_interp), has_het_ids(_has_het_ids) { // initialize DHT object // key size = count of key elements + timestep uint32_t key_size = (key_count + 1) * sizeof(Lookup_Keyelement); @@ -270,7 +271,7 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector &cell, const std::vector eval_vec = Rcpp::as>(hooks.dht_fuzz(input_nv)); assert(eval_vec.size() == this->key_count); - LookupKey vecFuzz(this->key_count + 1, {.0}); + LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0}); DHT_Rounder rounder; @@ -290,6 +291,9 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector &cell, } // add timestep to the end of the key as double value vecFuzz[this->key_count].fp_element = dt; + if (has_het_ids) { + vecFuzz[this->key_count + 1].fp_element = cell[0]; + } return vecFuzz; } @@ -297,7 +301,7 @@ LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector &cell, LookupKey DHT_Wrapper::fuzzForDHT(const std::vector &cell, double dt) { const auto c_zero_val = std::pow(10, AQUEOUS_EXP); - LookupKey vecFuzz(this->key_count + 1, {.0}); + LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0}); DHT_Rounder rounder; int totals_i = 0; @@ -323,6 +327,9 @@ LookupKey DHT_Wrapper::fuzzForDHT(const std::vector &cell, double dt) { } // add timestep to the end of the key as double value vecFuzz[this->key_count].fp_element = dt; + if (has_het_ids) { + vecFuzz[this->key_count + 1].fp_element = cell[0]; + } return vecFuzz; } diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp index e0b1c022d..9449692b4 100644 --- a/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp +++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.hpp @@ -87,7 +87,7 @@ public: const std::vector &key_indices, const std::vector &output_names, const InitialList::ChemistryHookFunctions &hooks, - uint32_t data_count, bool with_interp); + uint32_t data_count, bool with_interp, bool has_het_ids); /** * @brief Destroy the dht wrapper object * @@ -264,6 +264,7 @@ private: DHT_ResultObject dht_results; std::array base_totals{0}; + bool has_het_ids{false}; }; } // namespace poet diff --git a/src/poet.cpp b/src/poet.cpp index df51ce31c..c7864b4b6 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include @@ -510,6 +511,31 @@ std::array getBaseTotals(Field &&field, int root, MPI_Comm comm) { return base_totals; } +bool getHasID(Field &&field, int root, MPI_Comm comm) { + bool has_id; + + int rank; + MPI_Comm_rank(comm, &rank); + + const bool is_master = root == rank; + + if (is_master) { + const auto ID_field = field["ID"]; + + std::set unique_IDs(ID_field.begin(), ID_field.end()); + + has_id = unique_IDs.size() > 1; + + MPI_Bcast(&has_id, 1, MPI_C_BOOL, root, MPI_COMM_WORLD); + + return has_id; + } + + MPI_Bcast(&has_id, 1, MPI_C_BOOL, root, comm); + + return has_id; +} + int main(int argc, char *argv[]) { int world_size; @@ -558,6 +584,7 @@ int main(int argc, char *argv[]) { const ChemistryModule::SurrogateSetup surr_setup = { getSpeciesNames(init_list.getInitialGrid(), 0, MPI_COMM_WORLD), getBaseTotals(init_list.getInitialGrid(), 0, MPI_COMM_WORLD), + getHasID(init_list.getInitialGrid(), 0, MPI_COMM_WORLD), run_params.use_dht, run_params.dht_size, run_params.use_interp, From 50e856b28bd2b657ef21841cfadca6c864ec8d31 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 16 Dec 2024 10:36:37 +0100 Subject: [PATCH 11/21] feat: make output of H(0) and O(0) optional --- ext/iphreeqc | 2 +- src/Init/DiffusionInit.cpp | 12 +++++------- src/Init/GridInit.cpp | 7 +++---- src/Init/InitialList.cpp | 11 +++++++---- src/Init/InitialList.hpp | 11 +++++------ src/initializer.cpp | 21 ++++++++++++--------- 6 files changed, 33 insertions(+), 31 deletions(-) diff --git a/ext/iphreeqc b/ext/iphreeqc index 38268b4aa..2c8c31193 160000 --- a/ext/iphreeqc +++ b/ext/iphreeqc @@ -1 +1 @@ -Subproject commit 38268b4aad03e6ce4755315f4cd690f007fa2720 +Subproject commit 2c8c311934cd54f0923d8b92dcccd9894825ce1b diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp index 7abd987ec..7b5d9680b 100644 --- a/src/Init/DiffusionInit.cpp +++ b/src/Init/DiffusionInit.cpp @@ -3,7 +3,6 @@ // Rcpp include #include -#include #include "DataStructures/Field.hpp" #include "InitialList.hpp" @@ -15,7 +14,6 @@ #include #include #include -#include #include #include #include @@ -26,9 +24,9 @@ namespace poet { enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC }; const std::map tug_side_mapping = { - {tug::BC_SIDE_RIGHT , "E"}, - {tug::BC_SIDE_LEFT , "W"}, - {tug::BC_SIDE_TOP , "N"}, + {tug::BC_SIDE_RIGHT, "E"}, + {tug::BC_SIDE_LEFT, "W"}, + {tug::BC_SIDE_TOP, "N"}, {tug::BC_SIDE_BOTTOM, "S"}}; static std::vector colMajToRowMaj(const Rcpp::NumericVector &vec, @@ -166,8 +164,8 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list, if (boundaries_list.containsElementNamed(side.second.c_str())) { const Rcpp::List mapping = boundaries_list[side.second]; - const Rcpp::IntegerVector cells = mapping["cell"]; // MDL 2024-06-13 - const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL + const Rcpp::IntegerVector cells = mapping["cell"]; // MDL 2024-06-13 + const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL const Rcpp::CharacterVector type_str = mapping["type"]; if (cells.size() != values.size()) { diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp index 4a0e3fca3..6b0e4ae15 100644 --- a/src/Init/GridInit.cpp +++ b/src/Init/GridInit.cpp @@ -6,10 +6,8 @@ #include #include #include -#include #include #include -#include #include #include #include @@ -123,7 +121,8 @@ static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat, return Rcpp::Function("pqc_to_grid")(pqc_mat_R, grid_def); } -PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) { +PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input, + bool include_h0_o0) { // parse input values std::string script; std::string database; @@ -180,7 +179,7 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) { throw std::runtime_error("Grid size must be positive."); } - PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script); + PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script, include_h0_o0); this->transport_names = pqc_mat.getSolutionNames(); diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp index f4e87af14..88452a8c4 100644 --- a/src/Init/InitialList.cpp +++ b/src/Init/InitialList.cpp @@ -10,8 +10,9 @@ #include namespace poet { -void InitialList::initializeFromList(const Rcpp::List &setup) { - PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]); +void InitialList::initializeFromList(const Rcpp::List &setup, + bool include_h0_o0) { + PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key], include_h0_o0); initDiffusion(setup[diffusion_key], phreeqc); initChemistry(setup[chemistry_key]); } @@ -84,7 +85,8 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) { this->chem_hooks = Rcpp::as(setup[static_cast(ExportList::CHEM_HOOKS)]); - this->ai_surrogate_input_script = Rcpp::as(setup[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)]); + this->ai_surrogate_input_script = Rcpp::as( + setup[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)]); } Rcpp::List InitialList::exportList() { @@ -132,7 +134,8 @@ Rcpp::List InitialList::exportList() { out[static_cast(ExportList::CHEM_INTERP_SPECIES)] = Rcpp::wrap(this->interp_species); out[static_cast(ExportList::CHEM_HOOKS)] = this->chem_hooks; - out[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)] = this->ai_surrogate_input_script; + out[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)] = + this->ai_surrogate_input_script; return out; } diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp index a47d7918f..d3d1bf871 100644 --- a/src/Init/InitialList.hpp +++ b/src/Init/InitialList.hpp @@ -14,7 +14,6 @@ #include #include -#include #include #include #include @@ -27,14 +26,14 @@ using TugType = double; class InitialList { public: - InitialList(RInside &R) : R(R){}; + InitialList(RInside &R) : R(R) {}; - void initializeFromList(const Rcpp::List &setup); + void initializeFromList(const Rcpp::List &setup, bool include_h0_o0 = false); void importList(const Rcpp::List &setup, bool minimal = false); Rcpp::List exportList(); - Field getInitialGrid() const { return Field(this->initial_grid); } + Field getInitialGrid() const { return Field(this->initial_grid); } private: RInside &R; @@ -99,7 +98,7 @@ private: // std::unique_ptr pqc_mat; - PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input); + PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input, bool include_h0_o0); std::uint8_t dim{0}; @@ -196,7 +195,7 @@ private: NamedVector dht_species; NamedVector interp_species; - + // Path to R script that the user defines in the input file std::string ai_surrogate_input_script; diff --git a/src/initializer.cpp b/src/initializer.cpp index d613b83be..e03213c31 100644 --- a/src/initializer.cpp +++ b/src/initializer.cpp @@ -34,12 +34,15 @@ int main(int argc, char **argv) { "input script") ->default_val(false); - bool asRDS; - app.add_flag("-r, --rds", asRDS, "Save output as .rds") - ->default_val(false); + bool asRDS{false}; + app.add_flag("-r, --rds", asRDS, "Save output as .rds")->default_val(false); - bool asQS; - app.add_flag("-q, --qs", asQS, "Save output as .qs") + bool asQS{false}; + app.add_flag("-q, --qs", asQS, "Save output as .qs")->default_val(false); + + bool includeH0O0; + app.add_flag("--include-h0-o0", includeH0O0, + "Include H(0) and O(0) in the output") ->default_val(false); CLI11_PARSE(app, argc, argv); @@ -74,13 +77,13 @@ int main(int argc, char **argv) { // append the correct file extension if (asRDS) { - output_file += ".rds"; + output_file += ".rds"; } else if (asQS) { - output_file += ".qs"; + output_file += ".qs"; } else { - output_file += ".qs2"; + output_file += ".qs2"; } - + // set working directory to the directory of the input script if (setwd) { const std::string dir_path = Rcpp::as( From 034ec5d9f70a87a64d705fd8478ac1480302ef79 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 16 Dec 2024 14:33:50 +0100 Subject: [PATCH 12/21] fix: initialize includeH0O0 to false in command line options --- src/initializer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/initializer.cpp b/src/initializer.cpp index e03213c31..dd99b4ce4 100644 --- a/src/initializer.cpp +++ b/src/initializer.cpp @@ -40,7 +40,7 @@ int main(int argc, char **argv) { bool asQS{false}; app.add_flag("-q, --qs", asQS, "Save output as .qs")->default_val(false); - bool includeH0O0; + bool includeH0O0{false}; app.add_flag("--include-h0-o0", includeH0O0, "Include H(0) and O(0) in the output") ->default_val(false); From 6a889b53647890b09e880c81fda000aa27d99a5b Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 16 Dec 2024 19:53:13 +0100 Subject: [PATCH 13/21] feat: add dht_snaps and dht_out_dir parameters to ChemistryModule and main function --- src/Chemistry/ChemistryModule.hpp | 6 ++++++ src/poet.cpp | 2 ++ 2 files changed, 8 insertions(+) diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp index 697b50744..4bf925400 100644 --- a/src/Chemistry/ChemistryModule.hpp +++ b/src/Chemistry/ChemistryModule.hpp @@ -81,6 +81,8 @@ public: bool dht_enabled; std::uint32_t dht_size_mb; + int dht_snaps; + std::string dht_out_dir; bool interp_enabled; std::uint32_t interp_bucket_size; @@ -103,6 +105,10 @@ public: if (this->dht_enabled || this->interp_enabled) { this->initializeDHT(setup.dht_size_mb, this->params.dht_species, setup.has_het_ids); + + if (setup.dht_snaps != DHT_SNAPS_DISABLED) { + this->setDHTSnapshots(setup.dht_snaps, setup.dht_out_dir); + } } if (this->interp_enabled) { diff --git a/src/poet.cpp b/src/poet.cpp index c7864b4b6..742f4c395 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -587,6 +587,8 @@ int main(int argc, char *argv[]) { getHasID(init_list.getInitialGrid(), 0, MPI_COMM_WORLD), run_params.use_dht, run_params.dht_size, + run_params.dht_snaps, + run_params.out_dir, run_params.use_interp, run_params.interp_bucket_entries, run_params.interp_size, From 7ba9ddfa8360940296530304f573a4c6783ef664 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 17 Dec 2024 10:08:58 +0100 Subject: [PATCH 14/21] refactor: streamline WorkPackage constructor and improve DHT_Wrapper::fillDHT logic --- src/Chemistry/ChemistryDefs.hpp | 7 +- src/Chemistry/SurrogateModels/DHT_Wrapper.cpp | 71 ++++++++++--------- 2 files changed, 38 insertions(+), 40 deletions(-) diff --git a/src/Chemistry/ChemistryDefs.hpp b/src/Chemistry/ChemistryDefs.hpp index cc0aa5232..655eaf2d5 100644 --- a/src/Chemistry/ChemistryDefs.hpp +++ b/src/Chemistry/ChemistryDefs.hpp @@ -14,10 +14,7 @@ struct WorkPackage { std::vector> output; std::vector mapping; - WorkPackage(std::size_t _size) : size(_size) { - input.resize(size); - output.resize(size); - mapping.resize(size, CHEM_PQC); - } + WorkPackage(std::size_t _size) + : size(_size), input(size), output(size), mapping(size, CHEM_PQC) {} }; } // namespace poet \ No newline at end of file diff --git a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp index 0d9d2d2f6..4aa697e3e 100644 --- a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp +++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp @@ -129,42 +129,43 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) { dht_results.filledDHT = std::vector(length, false); for (int i = 0; i < length; i++) { // If true grid cell was simulated, needs to be inserted into dht - if (work_package.mapping[i] == CHEM_PQC) { - - // check if calcite or dolomite is absent and present, resp.n and vice - // versa in input/output. If this is the case -> Do not write to DHT! - // HACK: hardcoded, should be fixed! - if (hooks.dht_fill.isValid()) { - NamedVector old_values(output_names, work_package.input[i]); - NamedVector new_values(output_names, work_package.output[i]); - - if (hooks.dht_fill(old_values, new_values)) { - continue; - } - } - - uint32_t proc, index; - auto &key = dht_results.keys[i]; - const auto data = - (with_interp ? outputToInputAndRates(work_package.input[i], - work_package.output[i]) - : work_package.output[i]); - // void *data = (void *)&(work_package[i * this->data_count]); - // fuzz data (round, logarithm etc.) - - // insert simulated data with fuzzed key into DHT - int res = DHT_write(this->dht_object, key.data(), - const_cast(data.data()), &proc, &index); - - dht_results.locations[i] = {proc, index}; - - // if data was successfully written ... - if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) { - dht_evictions++; - } - - dht_results.filledDHT[i] = true; + if (work_package.mapping[i] != CHEM_PQC) { + continue; } + + // check if calcite or dolomite is absent and present, resp.n and vice + // versa in input/output. If this is the case -> Do not write to DHT! + // HACK: hardcoded, should be fixed! + if (hooks.dht_fill.isValid()) { + NamedVector old_values(output_names, work_package.input[i]); + NamedVector new_values(output_names, work_package.output[i]); + + if (hooks.dht_fill(old_values, new_values)) { + continue; + } + } + + uint32_t proc, index; + auto &key = dht_results.keys[i]; + const auto data = + (with_interp ? outputToInputAndRates(work_package.input[i], + work_package.output[i]) + : work_package.output[i]); + // void *data = (void *)&(work_package[i * this->data_count]); + // fuzz data (round, logarithm etc.) + + // insert simulated data with fuzzed key into DHT + int res = DHT_write(this->dht_object, key.data(), + const_cast(data.data()), &proc, &index); + + dht_results.locations[i] = {proc, index}; + + // if data was successfully written ... + if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) { + dht_evictions++; + } + + dht_results.filledDHT[i] = true; } } From 0e683d98cee9d25ad018cc145d62fbdbf5d280ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 17 Dec 2024 13:29:07 +0100 Subject: [PATCH 15/21] feat: enhance LookupKey with const isnan method and comparison operator --- src/Chemistry/SurrogateModels/LookupKey.hpp | 6 +++++- src/Chemistry/SurrogateModels/Rounding.hpp | 8 ++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/Chemistry/SurrogateModels/LookupKey.hpp b/src/Chemistry/SurrogateModels/LookupKey.hpp index 6e7b0681c..87c626f77 100644 --- a/src/Chemistry/SurrogateModels/LookupKey.hpp +++ b/src/Chemistry/SurrogateModels/LookupKey.hpp @@ -21,7 +21,7 @@ struct Lookup_SC_notation { return {SC_NOTATION_EXPONENT_MASK, SC_NOTATION_SIGNIFICANT_MASK}; } - constexpr bool isnan() { + constexpr bool isnan() const { return !!(exp == SC_NOTATION_EXPONENT_MASK && significant == SC_NOTATION_SIGNIFICANT_MASK); } @@ -35,6 +35,10 @@ union Lookup_Keyelement { return std::memcmp(this, &other, sizeof(Lookup_Keyelement)) == 0 ? true : false; } + + template bool operator>(const T &other) const { + return this->sc_notation.significant > other; + } }; class LookupKey : public std::vector { diff --git a/src/Chemistry/SurrogateModels/Rounding.hpp b/src/Chemistry/SurrogateModels/Rounding.hpp index 688ac4707..6656b8026 100644 --- a/src/Chemistry/SurrogateModels/Rounding.hpp +++ b/src/Chemistry/SurrogateModels/Rounding.hpp @@ -65,6 +65,14 @@ public: std::uint32_t signif) { Lookup_Keyelement new_val = value; + if (value.sc_notation.isnan()) { + return {.sc_notation = Lookup_SC_notation::nan()}; + } + + if (signif == 0) { + return {.sc_notation = {0, value > 0}}; + } + std::uint32_t diff_signif = static_cast( std::ceil(std::log10(std::abs(value.sc_notation.significant)))) - From b4c437e614126ba0ae97767ab12e10423be04157 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 18 Dec 2024 17:34:46 +0100 Subject: [PATCH 16/21] feat: add support for redox in PhreeqcMatrix initialization --- ext/iphreeqc | 2 +- src/Init/GridInit.cpp | 10 +++++++++- src/Init/InitialList.hpp | 7 ++++--- 3 files changed, 14 insertions(+), 5 deletions(-) diff --git a/ext/iphreeqc b/ext/iphreeqc index 2c8c31193..6e727e2f8 160000 --- a/ext/iphreeqc +++ b/ext/iphreeqc @@ -1 +1 @@ -Subproject commit 2c8c311934cd54f0923d8b92dcccd9894825ce1b +Subproject commit 6e727e2f896e853745b4dd123c5772a9b40ad705 diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp index 6b0e4ae15..52ce6943a 100644 --- a/src/Init/GridInit.cpp +++ b/src/Init/GridInit.cpp @@ -179,7 +179,15 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input, throw std::runtime_error("Grid size must be positive."); } - PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script, include_h0_o0); + bool with_redox = + grid_input.containsElementNamed( + GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX)) + ? Rcpp::as( + grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX)]) + : false; + + PhreeqcMatrix pqc_mat = + PhreeqcMatrix(database, script, include_h0_o0, with_redox); this->transport_names = pqc_mat.getSolutionNames(); diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp index d3d1bf871..03132983e 100644 --- a/src/Init/InitialList.hpp +++ b/src/Init/InitialList.hpp @@ -75,6 +75,7 @@ private: enum class GridMembers { PQC_SCRIPT_STRING, PQC_SCRIPT_FILE, + PQC_WITH_REDOX, PQC_DB_STRING, PQC_DB_FILE, GRID_DEF, @@ -88,9 +89,9 @@ private: static_cast(InitialList::GridMembers::ENUM_SIZE); static constexpr std::array - GridMembersString = {"pqc_in_string", "pqc_in_file", "pqc_db_string", - "pqc_db_file", "grid_def", "grid_size", - "constant_cells", "porosity"}; + GridMembersString = {"pqc_in_string", "pqc_in_file", "pqc_with_redox", + "pqc_db_string", "pqc_db_file", "grid_def", + "grid_size", "constant_cells", "porosity"}; constexpr const char *GRID_MEMBER_STR(GridMembers member) const { return GridMembersString[static_cast(member)]; From e1d1a577bcf2c93602a1fcd25c5eab9e9bc25f58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 19 Dec 2024 08:59:56 +0100 Subject: [PATCH 17/21] feat: add with_h0_o0 and with_redox parameters to ChemistryModule and related classes --- src/Chemistry/ChemistryModule.cpp | 3 ++- src/Init/ChemistryInit.cpp | 2 ++ src/Init/GridInit.cpp | 3 +++ src/Init/InitialList.cpp | 8 ++++++++ src/Init/InitialList.hpp | 6 ++++++ 5 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index 32c14be90..c95796b79 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -171,7 +171,8 @@ poet::ChemistryModule::ChemistryModule( if (!is_master) { PhreeqcMatrix pqc_mat = - PhreeqcMatrix(chem_params.database, chem_params.pqc_script); + PhreeqcMatrix(chem_params.database, chem_params.pqc_script, + chem_params.with_h0_o0, chem_params.with_redox); this->pqc_runner = std::make_unique(pqc_mat.subset(chem_params.pqc_ids)); diff --git a/src/Init/ChemistryInit.cpp b/src/Init/ChemistryInit.cpp index b399671ed..db0b098d0 100644 --- a/src/Init/ChemistryInit.cpp +++ b/src/Init/ChemistryInit.cpp @@ -69,6 +69,8 @@ InitialList::ChemistryInit InitialList::getChemistryInit() const { chem_init.database = database; chem_init.pqc_script = pqc_script; chem_init.pqc_ids = pqc_ids; + chem_init.with_h0_o0 = with_h0_o0; + chem_init.with_redox = with_redox; // chem_init.pqc_scripts = pqc_scripts; // chem_init.pqc_ids = pqc_ids; diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp index 52ce6943a..63b33a628 100644 --- a/src/Init/GridInit.cpp +++ b/src/Init/GridInit.cpp @@ -186,6 +186,9 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input, grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX)]) : false; + this->with_h0_o0 = include_h0_o0; + this->with_redox = with_redox; + PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script, include_h0_o0, with_redox); diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp index 88452a8c4..c6282977d 100644 --- a/src/Init/InitialList.cpp +++ b/src/Init/InitialList.cpp @@ -57,6 +57,10 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) { Rcpp::as(setup[static_cast(ExportList::CHEM_DATABASE)]); this->pqc_script = Rcpp::as( setup[static_cast(ExportList::CHEM_PQC_SCRIPT)]); + this->with_h0_o0 = + Rcpp::as(setup[static_cast(ExportList::CHEM_PQC_WITH_H0_O0)]); + this->with_redox = + Rcpp::as(setup[static_cast(ExportList::CHEM_PQC_WITH_REDOX)]); this->field_header = Rcpp::as>( setup[static_cast(ExportList::CHEM_FIELD_HEADER)]); this->pqc_ids = Rcpp::as>( @@ -113,6 +117,10 @@ Rcpp::List InitialList::exportList() { out[static_cast(ExportList::CHEM_DATABASE)] = Rcpp::wrap(this->database); out[static_cast(ExportList::CHEM_PQC_SCRIPT)] = Rcpp::wrap(this->pqc_script); + out[static_cast(ExportList::CHEM_PQC_WITH_H0_O0)] = + Rcpp::wrap(this->with_h0_o0); + out[static_cast(ExportList::CHEM_PQC_WITH_REDOX)] = + Rcpp::wrap(this->with_redox); out[static_cast(ExportList::CHEM_FIELD_HEADER)] = Rcpp::wrap(this->field_header); out[static_cast(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids); diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp index 03132983e..64c2fd08c 100644 --- a/src/Init/InitialList.hpp +++ b/src/Init/InitialList.hpp @@ -52,6 +52,8 @@ private: DIFFU_ALPHA_Y, CHEM_DATABASE, CHEM_PQC_SCRIPT, + CHEM_PQC_WITH_H0_O0, + CHEM_PQC_WITH_REDOX, CHEM_PQC_IDS, CHEM_FIELD_HEADER, // CHEM_PQC_SCRIPTS, @@ -190,6 +192,8 @@ private: std::string database; std::string pqc_script; + bool with_h0_o0{false}; + bool with_redox{false}; // std::vector pqc_scripts; std::vector pqc_ids; @@ -220,6 +224,8 @@ public: std::string database; std::string pqc_script; + bool with_h0_o0; + bool with_redox; std::vector pqc_ids; // std::map pqc_input; From af5d2a7a700719fff4ea42cc3ad8ec8ba4ce153f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 20 Dec 2024 07:33:00 +0100 Subject: [PATCH 18/21] fix: handle zero distance case in inverseDistanceWeighting to prevent division by zero --- src/Chemistry/ChemistryModule.cpp | 59 ++----------------------------- 1 file changed, 3 insertions(+), 56 deletions(-) diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index c95796b79..e6572c092 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -66,7 +66,8 @@ inverseDistanceWeighting(const std::vector &to_calc, distance += std::pow( rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2); } - weights[point_i] = 1 / std::sqrt(distance); + + weights[point_i] = distance != 0 ? 1 / std::sqrt(distance) : 0; assert(!std::isnan(weights[point_i])); inv_sum += weights[point_i]; } @@ -97,63 +98,9 @@ inverseDistanceWeighting(const std::vector &to_calc, key_delta /= inv_sum; results[output_comp_i] = from[output_comp_i] + key_delta; + assert(!std::isnan(results[output_comp_i])); } - // if (!has_h) { - // double new_val = 0; - // for (int j = 0; j < data_set_n; j++) { - // new_val += weights[j] * output[j][0]; - // } - // results[0] = new_val / inv_sum; - // } - - // if (!has_h) { - // double new_val = 0; - // for (int j = 0; j < data_set_n; j++) { - // new_val += weights[j] * output[j][1]; - // } - // results[1] = new_val / inv_sum; - // } - - // for (std::uint32_t i = 0; i < to_calc.size(); i++) { - // const std::uint32_t interp_i = to_calc[i]; - - // // rescale input between 0 and 1 - // for (int j = 0; j < input.size(); j++) { - // buffer[j] = input[j].at(i); - // } - - // buffer[buffer_size - 1] = from[interp_i]; - - // const double min = *std::min_element(buffer, buffer + buffer_size); - // const double max = *std::max_element(buffer, buffer + buffer_size); - - // for (int j = 0; j < input.size(); j++) { - // buffer[j] = ((max - min) != 0 ? (buffer[j] - min) / (max - min) : 1); - // } - // from_rescaled = - // ((max - min) != 0 ? (from[interp_i] - min) / (max - min) : 0); - - // double inv_sum = 0; - - // // calculate distances for each point - // for (int i = 0; i < input.size(); i++) { - // const double distance = std::pow(buffer[i] - from_rescaled, 2); - - // buffer[i] = distance > 0 ? (1 / std::sqrt(distance)) : 0; - // inv_sum += buffer[i]; - // } - // // calculate new values - // double new_val = 0; - // for (int i = 0; i < output.size(); i++) { - // new_val += buffer[i] * output[i][interp_i]; - // } - // results[interp_i] = new_val / inv_sum; - // if (std::isnan(results[interp_i])) { - // std::cout << "nan with new_val = " << output[0][i] << std::endl; - // } - // } - return results; } From f0eacc0f4003f887bb2ad8cf403a7ba6bf853df6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 20 Dec 2024 08:35:00 +0100 Subject: [PATCH 19/21] Revert "feat: make output of H(0) and O(0) optional" This reverts commit bdb5ce944fd001cb36c8369d8c792f2f4a54c7be. --- src/Init/DiffusionInit.cpp | 12 +++--- src/Init/GridInit.cpp | 75 ++++++++------------------------------ src/Init/InitialList.cpp | 11 ++---- src/Init/InitialList.hpp | 22 ++++------- src/initializer.cpp | 5 --- 5 files changed, 34 insertions(+), 91 deletions(-) diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp index 7b5d9680b..7abd987ec 100644 --- a/src/Init/DiffusionInit.cpp +++ b/src/Init/DiffusionInit.cpp @@ -3,6 +3,7 @@ // Rcpp include #include +#include #include "DataStructures/Field.hpp" #include "InitialList.hpp" @@ -14,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -24,9 +26,9 @@ namespace poet { enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC }; const std::map tug_side_mapping = { - {tug::BC_SIDE_RIGHT, "E"}, - {tug::BC_SIDE_LEFT, "W"}, - {tug::BC_SIDE_TOP, "N"}, + {tug::BC_SIDE_RIGHT , "E"}, + {tug::BC_SIDE_LEFT , "W"}, + {tug::BC_SIDE_TOP , "N"}, {tug::BC_SIDE_BOTTOM, "S"}}; static std::vector colMajToRowMaj(const Rcpp::NumericVector &vec, @@ -164,8 +166,8 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list, if (boundaries_list.containsElementNamed(side.second.c_str())) { const Rcpp::List mapping = boundaries_list[side.second]; - const Rcpp::IntegerVector cells = mapping["cell"]; // MDL 2024-06-13 - const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL + const Rcpp::IntegerVector cells = mapping["cell"]; // MDL 2024-06-13 + const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL const Rcpp::CharacterVector type_str = mapping["type"]; if (cells.size() != values.size()) { diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp index 63b33a628..d3ee24d5f 100644 --- a/src/Init/GridInit.cpp +++ b/src/Init/GridInit.cpp @@ -15,42 +15,6 @@ namespace poet { -// static Rcpp::NumericMatrix pqcMatToR(const PhreeqcMatrix &phreeqc, RInside -// &R) { - -// PhreeqcMatrix::STLExport phreeqc_mat = phreeqc.get(); - -// // PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat(); - -// // add "id" to the front of the column names - -// const std::size_t ncols = phreeqc_mat.names.size(); -// const std::size_t nrows = phreeqc_mat.values.size(); - -// phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID"); - -// Rcpp::NumericMatrix mat(nrows, ncols + 1); - -// for (std::size_t i = 0; i < nrows; i++) { -// mat(i, 0) = phreeqc_mat.ids[i]; -// for (std::size_t j = 0; j < ncols; ++j) { -// mat(i, j + 1) = phreeqc_mat.values[i][j]; -// } -// } - -// Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names); - -// return mat; -// } - -// static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix -// &mat, -// const Rcpp::NumericMatrix &grid) { -// Rcpp::Function pqc_to_grid_R("pqc_to_grid"); - -// return pqc_to_grid_R(mat, grid); -// } - static inline std::map replaceRawKeywordIDs(std::map raws) { for (auto &raw : raws) { @@ -64,26 +28,6 @@ replaceRawKeywordIDs(std::map raws) { return raws; } -// static inline uint32_t getSolutionCount(std::unique_ptr -// &phreeqc, -// const Rcpp::List &initial_grid) { -// PhreeqcInit::ModulesArray mod_array; -// Rcpp::Function unique_R("unique"); - -// std::vector row_ids = -// Rcpp::as>(unique_R(initial_grid["ID"])); - -// // std::vector sizes_vec(sizes.begin(), sizes.end()); - -// // Rcpp::Function modify_sizes("modify_module_sizes"); -// // sizes_vec = Rcpp::as>( -// // modify_sizes(sizes_vec, phreeqc_mat, initial_grid)); - -// // std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin()); - -// return phreeqc->getModuleSizes(row_ids)[POET_SOL]; -// } - static std::string readFile(const std::string &path) { std::string string_rpath(PATH_MAX, '\0'); @@ -121,8 +65,7 @@ static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat, return Rcpp::Function("pqc_to_grid")(pqc_mat_R, grid_def); } -PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input, - bool include_h0_o0) { +PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) { // parse input values std::string script; std::string database; @@ -186,11 +129,23 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input, grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_REDOX)]) : false; - this->with_h0_o0 = include_h0_o0; + bool with_h0_o0 = + grid_input.containsElementNamed( + GRID_MEMBER_STR(GridMembers::PQC_WITH_H0_O0)) + ? Rcpp::as( + grid_input[GRID_MEMBER_STR(GridMembers::PQC_WITH_H0_O0)]) + : false; + + if (with_h0_o0 && !with_redox) { + throw std::runtime_error( + "Output of H(0) and O(0) can only be used with redox."); + } + + this->with_h0_o0 = with_h0_o0; this->with_redox = with_redox; PhreeqcMatrix pqc_mat = - PhreeqcMatrix(database, script, include_h0_o0, with_redox); + PhreeqcMatrix(database, script, with_h0_o0, with_redox); this->transport_names = pqc_mat.getSolutionNames(); diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp index c6282977d..25203ebce 100644 --- a/src/Init/InitialList.cpp +++ b/src/Init/InitialList.cpp @@ -10,9 +10,8 @@ #include namespace poet { -void InitialList::initializeFromList(const Rcpp::List &setup, - bool include_h0_o0) { - PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key], include_h0_o0); +void InitialList::initializeFromList(const Rcpp::List &setup) { + PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]); initDiffusion(setup[diffusion_key], phreeqc); initChemistry(setup[chemistry_key]); } @@ -89,8 +88,7 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) { this->chem_hooks = Rcpp::as(setup[static_cast(ExportList::CHEM_HOOKS)]); - this->ai_surrogate_input_script = Rcpp::as( - setup[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)]); + this->ai_surrogate_input_script = Rcpp::as(setup[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)]); } Rcpp::List InitialList::exportList() { @@ -142,8 +140,7 @@ Rcpp::List InitialList::exportList() { out[static_cast(ExportList::CHEM_INTERP_SPECIES)] = Rcpp::wrap(this->interp_species); out[static_cast(ExportList::CHEM_HOOKS)] = this->chem_hooks; - out[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)] = - this->ai_surrogate_input_script; + out[static_cast(ExportList::AI_SURROGATE_INPUT_SCRIPT)] = this->ai_surrogate_input_script; return out; } diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp index 64c2fd08c..576828002 100644 --- a/src/Init/InitialList.hpp +++ b/src/Init/InitialList.hpp @@ -26,9 +26,9 @@ using TugType = double; class InitialList { public: - InitialList(RInside &R) : R(R) {}; + InitialList(RInside &R) : R(R){}; - void initializeFromList(const Rcpp::List &setup, bool include_h0_o0 = false); + void initializeFromList(const Rcpp::List &setup); void importList(const Rcpp::List &setup, bool minimal = false); Rcpp::List exportList(); @@ -56,14 +56,6 @@ private: CHEM_PQC_WITH_REDOX, CHEM_PQC_IDS, CHEM_FIELD_HEADER, - // CHEM_PQC_SCRIPTS, - // CHEM_PQC_SOLUTIONS, - // CHEM_PQC_SOLUTION_PRIMARY, - // CHEM_PQC_EXCHANGER, - // CHEM_PQC_KINETICS, - // CHEM_PQC_EQUILIBRIUM, - // CHEM_PQC_SURFACE_COMPS, - // CHEM_PQC_SURFACE_CHARGES, CHEM_DHT_SPECIES, CHEM_INTERP_SPECIES, CHEM_HOOKS, @@ -78,6 +70,7 @@ private: PQC_SCRIPT_STRING, PQC_SCRIPT_FILE, PQC_WITH_REDOX, + PQC_WITH_H0_O0, PQC_DB_STRING, PQC_DB_FILE, GRID_DEF, @@ -91,9 +84,10 @@ private: static_cast(InitialList::GridMembers::ENUM_SIZE); static constexpr std::array - GridMembersString = {"pqc_in_string", "pqc_in_file", "pqc_with_redox", - "pqc_db_string", "pqc_db_file", "grid_def", - "grid_size", "constant_cells", "porosity"}; + GridMembersString = {"pqc_in_string", "pqc_in_file", "pqc_with_redox", + "pqc_wth_h0_o0", "pqc_db_string", "pqc_db_file", + "grid_def", "grid_size", "constant_cells", + "porosity"}; constexpr const char *GRID_MEMBER_STR(GridMembers member) const { return GridMembersString[static_cast(member)]; @@ -101,7 +95,7 @@ private: // std::unique_ptr pqc_mat; - PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input, bool include_h0_o0); + PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input); std::uint8_t dim{0}; diff --git a/src/initializer.cpp b/src/initializer.cpp index dd99b4ce4..dd7e4ddd0 100644 --- a/src/initializer.cpp +++ b/src/initializer.cpp @@ -40,11 +40,6 @@ int main(int argc, char **argv) { bool asQS{false}; app.add_flag("-q, --qs", asQS, "Save output as .qs")->default_val(false); - bool includeH0O0{false}; - app.add_flag("--include-h0-o0", includeH0O0, - "Include H(0) and O(0) in the output") - ->default_val(false); - CLI11_PARSE(app, argc, argv); // source the input script From cfea048fdbd361663e89405cd3c67b08cf450dd5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 20 Dec 2024 09:18:35 +0100 Subject: [PATCH 20/21] fix: update dht_species values for barite and celestite in benchmark scripts --- bench/barite/barite_200.R | 18 +++++++++--------- bench/barite/barite_50ai.R | 24 ++++++++++++------------ 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/bench/barite/barite_200.R b/bench/barite/barite_200.R index a6c87c749..09483b78b 100644 --- a/bench/barite/barite_200.R +++ b/bench/barite/barite_200.R @@ -35,15 +35,15 @@ diffusion_setup <- list( ) dht_species <- c( - "H" = 7, - "O" = 7, - "Charge" = 4, - "Ba" = 7, - "Cl" = 7, - "S(6)" = 7, - "Sr" = 7, - "Barite" = 4, - "Celestite" = 4 + "H" = 3, + "O" = 3, + "Charge" = 6, + "Ba" = 6, + "Cl" = 6, + "S" = 6, + "Sr" = 6, + "Barite" = 5, + "Celestite" = 5 ) chemistry_setup <- list( diff --git a/bench/barite/barite_50ai.R b/bench/barite/barite_50ai.R index c2a674a85..8dccaa88e 100644 --- a/bench/barite/barite_50ai.R +++ b/bench/barite/barite_50ai.R @@ -11,9 +11,9 @@ grid_def <- matrix(2, nrow = rows, ncol = cols) grid_setup <- list( pqc_in_file = "./barite.pqi", pqc_db_file = "./db_barite.dat", ## Path to the database file for Phreeqc - grid_def = grid_def, ## Definition of the grid, containing IDs according to the Phreeqc input script - grid_size = c(s_rows, s_cols), ## Size of the grid in meters - constant_cells = c() ## IDs of cells with constant concentration + grid_def = grid_def, ## Definition of the grid, containing IDs according to the Phreeqc input script + grid_size = c(s_rows, s_cols), ## Size of the grid in meters + constant_cells = c() ## IDs of cells with constant concentration ) bound_length <- 2 @@ -36,15 +36,15 @@ diffusion_setup <- list( ) dht_species <- c( - "H" = 4, - "O" = 10, - "Charge" = 4, - "Ba" = 7, - "Cl" = 4, - "S(6)" = 7, - "Sr" = 4, - "Barite" = 2, - "Celestite" = 2 + "H" = 3, + "O" = 3, + "Charge" = 3, + "Ba" = 6, + "Cl" = 6, + "S" = 6, + "Sr" = 6, + "Barite" = 5, + "Celestite" = 5 ) chemistry_setup <- list( From b9e18a1059580a0e51b1a519f67cf003b1c780e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 20 Dec 2024 09:18:44 +0100 Subject: [PATCH 21/21] fix: initialize default values for RuntimeParameters in poet.hpp.in --- src/poet.hpp.in | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/poet.hpp.in b/src/poet.hpp.in index c48a0f3df..9462f6d7e 100644 --- a/src/poet.hpp.in +++ b/src/poet.hpp.in @@ -23,7 +23,6 @@ #pragma once #include -#include #include #include @@ -47,27 +46,27 @@ struct RuntimeParameters { // MDL added to accomodate for qs::qsave/qread bool as_rds = false; - bool as_qs = false; - std::string out_ext; + bool as_qs = false; + std::string out_ext; bool print_progress = false; static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32; - std::uint32_t work_package_size; + std::uint32_t work_package_size = WORK_PACKAGE_SIZE_DEFAULT; bool use_dht = false; static constexpr std::uint32_t DHT_SIZE_DEFAULT = 1.5E3; - std::uint32_t dht_size; + std::uint32_t dht_size = DHT_SIZE_DEFAULT; static constexpr std::uint8_t DHT_SNAPS_DEFAULT = 0; - std::uint8_t dht_snaps; + std::uint8_t dht_snaps = DHT_SNAPS_DEFAULT; bool use_interp = false; static constexpr std::uint32_t INTERP_SIZE_DEFAULT = 100; - std::uint32_t interp_size; + std::uint32_t interp_size = INTERP_SIZE_DEFAULT; static constexpr std::uint32_t INTERP_MIN_ENTRIES_DEFAULT = 5; - std::uint32_t interp_min_entries; + std::uint32_t interp_min_entries = INTERP_MIN_ENTRIES_DEFAULT; static constexpr std::uint32_t INTERP_BUCKET_ENTRIES_DEFAULT = 20; - std::uint32_t interp_bucket_entries; + std::uint32_t interp_bucket_entries = INTERP_BUCKET_ENTRIES_DEFAULT; bool use_ai_surrogate = false; };