From 3fdf586e0d7d0707d306d9eafa4d251484811782 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 17 Apr 2023 12:31:39 +0200 Subject: [PATCH 1/2] BREAKING CHANGE: use dump mechanism of PhreeqcRM-GFZ to get and set internal variables feat: enables exchange data: added exchange only benchmark data: applied required changes to benchmarks --- app/poet.cpp | 7 +- app/poet_prm.cpp | 2 - bench/barite/barite.R | 19 +- bench/dolo_diffu_inner/dolo_diffu_inner.R | 23 +- .../dolo_diffu_inner/dolo_diffu_inner_large.R | 13 +- bench/surfex/CMakeLists.txt | 6 +- bench/surfex/ExBase.pqi | 39 +++ bench/surfex/ex.R | 144 +++++++++ ext/phreeqcrm | 2 +- include/poet/ChemistryModule.hpp | 45 +-- include/poet/SimParams.hpp | 3 +- src/ChemistryModule/ChemistryModule.cpp | 273 ++++++++---------- src/ChemistryModule/MasterFunctions.cpp | 7 +- src/ChemistryModule/WorkerFunctions.cpp | 27 +- src/DiffusionModule.cpp | 7 +- src/Grid.cpp | 1 - src/SimParams.cpp | 5 +- 17 files changed, 372 insertions(+), 251 deletions(-) create mode 100644 bench/surfex/ExBase.pqi create mode 100644 bench/surfex/ex.R diff --git a/app/poet.cpp b/app/poet.cpp index 097c9ff33..1af116370 100644 --- a/app/poet.cpp +++ b/app/poet.cpp @@ -21,7 +21,6 @@ #include "poet/ChemistryModule.hpp" #include #include -#include #include #include #include @@ -100,7 +99,8 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, Grid &grid, ChemistryParams &chem_params, const GridParams &g_params, uint32_t nxyz_master) { - DiffusionModule diffusion(poet::DiffusionParams(R), grid); + DiffusionParams d_params{R}; + DiffusionModule diffusion(d_params, grid); /* Iteration Count is dynamic, retrieving value from R (is only needed by * master for the following loop) */ uint32_t maxiter = R.parseEval("mysetup$iterations"); @@ -112,7 +112,8 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, Grid &grid, set_chem_parameters(chem, nxyz_master, chem_params.database_path); chem.RunInitFile(chem_params.input_script); - chem.InitializeField(grid.GetTotalCellCount(), DFToHashMap(g_params.init_df)); + poet::ChemistryModule::SingleCMap init_df = DFToHashMap(d_params.initial_t); + chem.mergeFieldWithModule(init_df, grid.GetTotalCellCount()); if (params.getNumParams().dht_enabled) { chem.SetDHTEnabled(true, params.getNumParams().dht_size_per_process); diff --git a/app/poet_prm.cpp b/app/poet_prm.cpp index a5eb690b9..73734dd3d 100644 --- a/app/poet_prm.cpp +++ b/app/poet_prm.cpp @@ -21,7 +21,6 @@ #include "poet/ChemistryModule.hpp" #include #include -#include #include #include #include @@ -99,7 +98,6 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, Grid &grid, set_chem_parameters(chem, nxyz_master, chem_params.database_path); chem.RunInitFile(chem_params.input_script); - chem.SetSelectedOutputOn(true); chem.SetTimeStep(0); chem.RunCells(); diff --git a/bench/barite/barite.R b/bench/barite/barite.R index 1a215f467..9b0be3971 100644 --- a/bench/barite/barite.R +++ b/bench/barite/barite.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-04-13 17:10:30 mluebke" +## Time-stamp: "Last modified 2023-04-14 11:02:31 mluebke" database <- normalizePath("../share/poet/bench/barite/db_barite.dat") input_script <- normalizePath("../share/poet/bench/barite/barite.pqi") @@ -29,7 +29,7 @@ grid <- list( n_cells = c(n, m), s_cells = c(1, 1), type = types[1], - init_cell = as.data.frame(init_cell), + init_cell = as.data.frame(init_cell, check.names = FALSE), props = names(init_cell), database = database, input_script = input_script @@ -43,13 +43,13 @@ grid <- list( ## initial conditions -init_diffu <- c( +init_diffu <- list( "H" = 110.0124, "O" = 55.5087, "Charge" = -1.217E-09, "Ba" = 1.E-10, "Cl" = 1.E-10, - "S" = 6.205E-4, + "S(6)" = 6.205E-4, "Sr" = 6.205E-4 ) @@ -60,7 +60,7 @@ injection_diff <- list( "Charge" = -3.337E-08, "Ba" = 0.1, "Cl" = 0.2, - "S" = 0, + "S(6)" = 0, "Sr" = 0) ) @@ -71,7 +71,7 @@ alpha_diffu <- c( "Charge" = 1E-06, "Ba" = 1E-06, "Cl" = 1E-06, - "S" = 1E-06, + "S(6)" = 1E-06, "Sr" = 1E-06 ) @@ -91,9 +91,12 @@ boundary <- list( diffu_list <- names(alpha_diffu) +vecinj <- do.call(rbind.data.frame, injection_diff) +names(vecinj) <- names(init_diffu) + diffusion <- list( - init = init_diffu, - vecinj = do.call(rbind.data.frame, injection_diff), + init = as.data.frame(init_diffu, check.names = FALSE), + vecinj = vecinj, # vecinj_inner = vecinj_inner, vecinj_index = boundary, alpha = alpha_diffu diff --git a/bench/dolo_diffu_inner/dolo_diffu_inner.R b/bench/dolo_diffu_inner/dolo_diffu_inner.R index ea9ade839..33c5114a2 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-01-10 13:51:40 delucia" +## Time-stamp: "Last modified 2023-04-17 12:25:25 mluebke" database <- normalizePath("../share/poet/examples/phreeqc_kin.dat") input_script <- normalizePath("../share/poet/bench/dolo_inner.pqi") @@ -17,7 +17,7 @@ init_cell <- list( "H" = 110.683, "O" = 55.3413, "Charge" = -5.0822e-19, - "C" = 1.2279E-4, + "C(4)" = 1.2279E-4, "Ca" = 1.2279E-4, "Cl" = 0, "Mg" = 0, @@ -30,7 +30,7 @@ grid <- list( n_cells = c(n, m), s_cells = c(1, 1), type = types[1], - init_cell = as.data.frame(init_cell), + init_cell = as.data.frame(init_cell, check.names = FALSE), props = names(init_cell), database = database, input_script = input_script @@ -43,11 +43,11 @@ grid <- list( ################################################################## ## initial conditions -init_diffu <- c( +init_diffu <- list( "H" = 110.683, "O" = 55.3413, "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 @@ -70,7 +70,7 @@ vecinj_diffu <- list( "H" = 110.683, "O" = 55.3413, "Charge" = 1.90431e-16, - "C" = 0, + "C(4)" = 0, "Ca" = 0, "Cl" = 0.002, "Mg" = 0.001 @@ -79,7 +79,7 @@ vecinj_diffu <- list( "H" = 110.683, "O" = 55.3413, "Charge" = 1.90431e-16, - "C" = 0, + "C(4)" = 0, "Ca" = 0.0, "Cl" = 0.004, "Mg" = 0.002 @@ -102,9 +102,12 @@ boundary <- list( diffu_list <- names(alpha_diffu) +vecinj <- do.call(rbind.data.frame, vecinj_diffu) +names(vecinj) <- names(init_diffu) + diffusion <- list( - init = init_diffu, - vecinj = do.call(rbind.data.frame, vecinj_diffu), + init = as.data.frame(init_diffu, check.names = FALSE), + vecinj = vecinj, vecinj_inner = vecinj_inner, vecinj_index = boundary, alpha = alpha_diffu diff --git a/bench/dolo_diffu_inner/dolo_diffu_inner_large.R b/bench/dolo_diffu_inner/dolo_diffu_inner_large.R index 57726d925..40f1da3f2 100644 --- a/bench/dolo_diffu_inner/dolo_diffu_inner_large.R +++ b/bench/dolo_diffu_inner/dolo_diffu_inner_large.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-01-10 13:51:40 delucia" +## Time-stamp: "Last modified 2023-04-17 12:27:27 mluebke" database <- normalizePath("../share/poet/examples/phreeqc_kin.dat") input_script <- normalizePath("../share/poet/bench/dolo_inner.pqi") @@ -30,7 +30,7 @@ grid <- list( n_cells = c(n, m), s_cells = c(2, 1), type = types[1], - init_cell = as.data.frame(init_cell), + init_cell = as.data.frame(init_cell, check.names = FALSE), props = names(init_cell), database = database, input_script = input_script @@ -43,7 +43,7 @@ grid <- list( ################################################################## ## initial conditions -init_diffu <- c( +init_diffu <- list( "H" = 110.683, "O" = 55.3413, "Charge" = -5.0822e-19, @@ -102,9 +102,12 @@ boundary <- list( diffu_list <- names(alpha_diffu) +vecinj <- do.call(rbind.data.frame, vecinj_diffu) +names(vecinj) <- names(init_diffu) + diffusion <- list( - init = init_diffu, - vecinj = do.call(rbind.data.frame, vecinj_diffu), + init = as.data.frame(init_diffu, check.names = FALSE), + vecinj = vecinj, vecinj_inner = vecinj_inner, vecinj_index = boundary, alpha = alpha_diffu diff --git a/bench/surfex/CMakeLists.txt b/bench/surfex/CMakeLists.txt index bda3f2353..8b3f44d5b 100644 --- a/bench/surfex/CMakeLists.txt +++ b/bench/surfex/CMakeLists.txt @@ -1,6 +1,8 @@ install(FILES - surfex.R - SurfExBase.pqi + ExBase.pqi + ex.R + #surfex.R + #SurfExBase.pqi SMILE_2021_11_01_TH.dat DESTINATION share/poet/bench/surfex diff --git a/bench/surfex/ExBase.pqi b/bench/surfex/ExBase.pqi new file mode 100644 index 000000000..af9c945a3 --- /dev/null +++ b/bench/surfex/ExBase.pqi @@ -0,0 +1,39 @@ +## Time-stamp: "Last modified 2023-03-21 11:49:43 mluebke" +KNOBS + -logfile false + -iterations 10000 + -convergence_tolerance 1E-12 + -step_size 2 + -pe_step_size 2 +SELECTED_OUTPUT + -reset false + -high_precision true + -solution true + -state true + -step true + -pH true + -pe true + -ionic_strength true + -water true +SOLUTION 1 +temp 13 +units mol/kgw +pH 7.06355 +pe -2.626517 +C(4) 0.001990694 +Ca 0.02172649 +Cl 0.3227673 charge +Fe 0.0001434717 +K 0.001902357 +Mg 0.01739704 +Na 0.2762882 +S(6) 0.01652701 +Sr 0.0004520361 +U(4) 8.147792e-12 +U(6) 2.237946e-09 +-water 0.00133 +EXCHANGE 1 + -equil 1 + Z 0.0012585 + Y 0.0009418 +END diff --git a/bench/surfex/ex.R b/bench/surfex/ex.R new file mode 100644 index 000000000..299c5db0c --- /dev/null +++ b/bench/surfex/ex.R @@ -0,0 +1,144 @@ +## Time-stamp: "Last modified 2023-04-17 12:29:27 mluebke" + +database <- normalizePath("./SMILE_2021_11_01_TH.dat") +input_script <- normalizePath("./ExBase.pqi") + +cat(paste(":: R This is a test 1\n")) + +################################################################# +## Section 1 ## +## Grid initialization ## +################################################################# + +n <- 100 +m <- 100 + +types <- c("scratch", "phreeqc", "rds") + +init_cell <- list(H = 1.476571028625e-01, + O = 7.392297218936e-02, + Charge = -1.765225732724e-18, + `C(-4)` = 2.477908970828e-21, + `C(4)` = 2.647623016916e-06, + Ca = 2.889623169138e-05, + Cl = 4.292806181039e-04, + `Fe(2)` =1.908142472666e-07, + `Fe(3)` =3.173306589931e-12, + `H(0)` =2.675642675119e-15, + K = 2.530134809667e-06, + Mg =2.313806319294e-05, + Na =3.674633059628e-04, + `S(-2)` = 8.589766637180e-15, + `S(2)` = 1.205284362720e-19, + `S(4)` = 9.108958772790e-18, + `S(6)` = 2.198092329098e-05, + Sr = 6.012080128154e-07, + `U(4)` = 1.039668623852e-14, + `U(5)` = 1.208394829796e-15, + `U(6)` = 2.976409147150e-12) + +grid <- list( + n_cells = c(n, m), + s_cells = c(1, 1), + type = "scratch", + init_cell = as.data.frame(init_cell, check.names = FALSE), + props = names(init_cell), + database = database, + input_script = input_script +) + + +################################################################## +## Section 2 ## +## Diffusion parameters and boundary conditions ## +################################################################## + +vecinj_diffu <- list( + list(H = 0.147659686316291, + O = 0.0739242798146046, + Charge = 7.46361643222701e-20, + `C(-4)` = 2.92438561098248e-21, + `C(4)` = 2.65160558871092e-06, + Ca = 2.89001071336443e-05, + Cl = 0.000429291158114428, + `Fe(2)` = 1.90823391198114e-07, + `Fe(3)` = 3.10832423034763e-12, + `H(0)` = 2.7888235127385e-15, + K = 2.5301787e-06, + Mg = 2.31391999937907e-05, + Na = 0.00036746969, + `S(-2)` = 1.01376078438546e-14, + `S(2)` = 1.42247026981542e-19, + `S(4)` = 9.49422092568557e-18, + `S(6)` = 2.19812504654191e-05, + Sr = 6.01218519999999e-07, + `U(4)` = 4.82255946569383e-12, + `U(5)` = 5.49050615347901e-13, + `U(6)` = 1.32462838991902e-09) +) + +vecinj <- do.call(rbind.data.frame, vecinj_diffu) +names(vecinj) <- grid$props + +## diffusion coefficients +alpha_diffu <- c(H = 1E-6, O = 1E-6, Charge = 1E-6, `C(-4)` = 1E-6, + `C(4)` = 1E-6, Ca = 1E-6, Cl = 1E-6, `Fe(2)` = 1E-6, + `Fe(3)` = 1E-6, `H(0)` = 1E-6, K = 1E-6, Mg = 1E-6, + Na = 1E-6, `S(-2)` = 1E-6, `S(2)` = 1E-6, + `S(4)` = 1E-6, `S(6)` = 1E-6, Sr = 1E-6, + `U(4)` = 1E-6, `U(5)` = 1E-6, `U(6)` = 1E-6) + +## list of boundary conditions/inner nodes + +## vecinj_inner <- list( +## list(1,1,1) +## ) + +boundary <- list( + "N" = rep(1, n), + "E" = rep(0, n), + "S" = rep(0, n), + "W" = rep(0, n) +) + +diffu_list <- names(alpha_diffu) + +vecinj <- do.call(rbind.data.frame, vecinj_diffu) +names(vecinj) <- names(init_cell) + +diffusion <- list( + init = as.data.frame(init_cell, check.names = FALSE), + vecinj = vecinj, +# vecinj_inner = vecinj_inner, + vecinj_index = boundary, + alpha = alpha_diffu +) + +################################################################# +## Section 3 ## +## Chemistry module (Phreeqc) ## +################################################################# + + +chemistry <- list( + database = database, + input_script = input_script +) + +################################################################# +## Section 4 ## +## Putting all those things together ## +################################################################# + + +iterations <- 10 +dt <- 200 + +setup <- list( + grid = grid, + diffusion = diffusion, + chemistry = chemistry, + iterations = iterations, + timesteps = rep(dt, iterations), + store_result = TRUE +) diff --git a/ext/phreeqcrm b/ext/phreeqcrm index f862889b6..ade4f5660 160000 --- a/ext/phreeqcrm +++ b/ext/phreeqcrm @@ -1 +1 @@ -Subproject commit f862889b693507458d450c6319abe5e1dce030f7 +Subproject commit ade4f56605b8d0e81b7ccad510f2db21f46b5ff9 diff --git a/include/poet/ChemistryModule.hpp b/include/poet/ChemistryModule.hpp index bab677850..f37f1f2de 100644 --- a/include/poet/ChemistryModule.hpp +++ b/include/poet/ChemistryModule.hpp @@ -5,8 +5,10 @@ #include "PhreeqcRM.h" #include "poet/DHT_Wrapper.hpp" +#include #include #include +#include #include #include #include @@ -127,35 +129,14 @@ public: }; /** - * Initializes field with a "DataFrame" containing a single value for each - * species. + * Merge initial values from existing module with the chemistry module and set + * according internal variables. * - * Each species must have been previously found by - * ChemistryModule::RunInitFile. - * - * \exception std::domain_error Species name of map is not defined. - * - * \param n_cells Count of cells for each species. - * \param mapped_values Unordered map containing one double value for each - * specified species. + * \param input_map Map with name and initial values of another module like + * Diffusion. + * \param n_cells number of cells used to initialize the field with. */ - void InitializeField(uint32_t n_cells, const SingleCMap &mapped_values); - - /** - * Initializes field with a "DataFrame" containing a vector of values for each - * species. - * - * Each species must have been previously found by - * ChemistryModule::RunInitFile. - * - * There is no check if vector length matches count of grid cells defined. - * - * \exception std::domain_error Species name of map is not defined. - * - * \param mapped_values Unordered map containing a vector of multiple values. - * Size of the vectors shall be the count of grid cells defined previously. - */ - void InitializeField(const VectorCMap &mapped_values); + void mergeFieldWithModule(const SingleCMap &input_map, std::uint32_t n_cells); /** * **Only called by workers!** Start the worker listening loop. @@ -298,6 +279,7 @@ protected: enum { CHEM_INIT, CHEM_WP_SIZE, + CHEM_INIT_SPECIES, CHEM_DHT_ENABLE, CHEM_DHT_SIGNIF_VEC, CHEM_DHT_PROP_TYPE_VEC, @@ -365,9 +347,6 @@ protected: std::vector &vecMapping, double dSimTime, double dTimestep); - void GetWPFromInternals(std::vector &vecWP, uint32_t wp_size); - void SetInternalsFromWP(const std::vector &vecWP, uint32_t wp_size); - std::vector CalculateWPSizesVector(uint32_t n_cells, uint32_t wp_size) const; @@ -397,6 +376,7 @@ protected: double idle_t = 0.; double seq_t = 0.; double send_recv_t = 0.; + #endif double chem_t = 0.; @@ -405,8 +385,9 @@ protected: uint32_t prop_count = 0; std::vector prop_names; std::vector field; - std::vector speciesPerModule; - static constexpr uint32_t MODULE_COUNT = 5; + static constexpr uint32_t MODULE_COUNT = 4; + + std::array speciesPerModule{}; }; } // namespace poet diff --git a/include/poet/SimParams.hpp b/include/poet/SimParams.hpp index b68dbce5d..d4f2efb76 100644 --- a/include/poet/SimParams.hpp +++ b/include/poet/SimParams.hpp @@ -21,7 +21,6 @@ #ifndef PARSER_H #define PARSER_H -#include #include #include #include @@ -85,7 +84,7 @@ using GridParams = struct s_GridParams { }; using DiffusionParams = struct s_DiffusionParams { - std::vector prop_names; + Rcpp::DataFrame initial_t; Rcpp::NumericVector alpha; Rcpp::List vecinj_inner; diff --git a/src/ChemistryModule/ChemistryModule.cpp b/src/ChemistryModule/ChemistryModule.cpp index c3d22af83..d46ddf60a 100644 --- a/src/ChemistryModule/ChemistryModule.cpp +++ b/src/ChemistryModule/ChemistryModule.cpp @@ -2,11 +2,15 @@ #include "PhreeqcRM.h" #include "poet/DHT_Wrapper.hpp" +#include +#include #include #include +#include #include #include #include +#include #include #ifndef POET_USE_PRM @@ -58,44 +62,12 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { this->FindComponents(); - std::vector props; - this->speciesPerModule.reserve(this->MODULE_COUNT); + PhreeqcRM::initializePOET(this->speciesPerModule, this->prop_names); + this->prop_count = prop_names.size(); - std::vector curr_prop_names = this->GetComponents(); - props.insert(props.end(), curr_prop_names.begin(), curr_prop_names.end()); - this->speciesPerModule.push_back(curr_prop_names.size()); - - curr_prop_names = this->GetEquilibriumPhases(); - props.insert(props.end(), curr_prop_names.begin(), curr_prop_names.end()); - char equilibrium = (curr_prop_names.empty() ? -1 : 1); - this->speciesPerModule.push_back(curr_prop_names.size()); - - curr_prop_names = this->GetExchangeNames(); - props.insert(props.end(), curr_prop_names.begin(), curr_prop_names.end()); - char exchange = (curr_prop_names.empty() ? -1 : 1); - this->speciesPerModule.push_back(curr_prop_names.size()); - - curr_prop_names = this->GetSurfaceNames(); - props.insert(props.end(), curr_prop_names.begin(), curr_prop_names.end()); - char surface = (curr_prop_names.empty() ? -1 : 1); - this->speciesPerModule.push_back(curr_prop_names.size()); - - // curr_prop_names = this->GetGasComponents(); - // props.insert(props.end(), curr_prop_names.begin(), curr_prop_names.end()); - // char gas = (curr_prop_names.empty() ? -1 : 1); - // this->speciesPerModule.push_back(curr_prop_names.size()); - - // curr_prop_names = this->GetSolidSolutionNames(); - // props.insert(props.end(), curr_prop_names.begin(), curr_prop_names.end()); - // char ssolutions = (curr_prop_names.empty() ? -1 : 1); - // this->speciesPerModule.push_back(curr_prop_names.size()); - - curr_prop_names = this->GetKineticReactions(); - props.insert(props.end(), curr_prop_names.begin(), curr_prop_names.end()); - char kinetics = (curr_prop_names.empty() ? -1 : 1); - this->speciesPerModule.push_back(curr_prop_names.size()); - - this->prop_count = props.size(); + char exchange = (speciesPerModule[1] == 0 ? -1 : 1); + char kinetics = (speciesPerModule[2] == 0 ? -1 : 1); + char equilibrium = (speciesPerModule[3] == 0 ? -1 : 1); #ifdef POET_USE_PRM std::vector ic1; @@ -105,7 +77,7 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { ic1[i] = 1; // Solution 1 ic1[nxyz + i] = equilibrium; // Equilibrium 1 ic1[2 * nxyz + i] = exchange; // Exchange none - ic1[3 * nxyz + i] = surface; // Surface none + ic1[3 * nxyz + i] = -1; // Surface none ic1[4 * nxyz + i] = -1; // Gas phase none ic1[5 * nxyz + i] = -1; // Solid solutions none ic1[6 * nxyz + i] = kinetics; // Kinetics 1 @@ -113,61 +85,105 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { this->InitialPhreeqc2Module(ic1); #else - if (!this->is_master || this->is_sequential) { - std::vector ic1; - ic1.resize(this->nxyz * 7, -1); - // TODO: hardcoded reaction modules - for (int i = 0; i < nxyz; i++) { - ic1[i] = 1; // Solution 1 - ic1[nxyz + i] = equilibrium; // Equilibrium 1 - ic1[2 * nxyz + i] = exchange; // Exchange none - ic1[3 * nxyz + i] = surface; // Surface none - ic1[4 * nxyz + i] = -1; // Gas phase none - ic1[5 * nxyz + i] = -1; // Solid solutions none - ic1[6 * nxyz + i] = kinetics; // Kinetics 1 - } - - this->InitialPhreeqc2Module(ic1); + std::vector ic1; + ic1.resize(this->nxyz * 7, -1); + // TODO: hardcoded reaction modules + for (int i = 0; i < nxyz; i++) { + ic1[i] = 1; // Solution 1 + ic1[nxyz + i] = equilibrium; // Equilibrium 1 + ic1[2 * nxyz + i] = exchange; // Exchange none + ic1[3 * nxyz + i] = -1; // Surface none + ic1[4 * nxyz + i] = -1; // Gas phase none + ic1[5 * nxyz + i] = -1; // Solid solutions none + ic1[6 * nxyz + i] = kinetics; // Kinetics 1 } + + this->InitialPhreeqc2Module(ic1); #endif - this->prop_names = props; } #ifndef POET_USE_PRM -void poet::ChemistryModule::InitializeField( - uint32_t n_cells, const poet::ChemistryModule::SingleCMap &mapped_values) { - if (this->is_master) { - this->field.reserve(this->prop_count * n_cells); - for (const std::string &prop_name : this->prop_names) { - const auto m_it = mapped_values.find(prop_name); - if (m_it == mapped_values.end()) { - throw std::domain_error( - "Prop names vector does not match any key in given map."); - } +void poet::ChemistryModule::mergeFieldWithModule(const SingleCMap &input_map, + std::uint32_t n_cells) { - const std::vector field_row(n_cells, m_it->second); - const auto chem_field_end = this->field.end(); - this->field.insert(chem_field_end, field_row.begin(), field_row.end()); - } - this->n_cells = n_cells; + if (is_master) { + int f_type = CHEM_INIT_SPECIES; + PropagateFunctionType(f_type); } -} -void poet::ChemistryModule::InitializeField( - const poet::ChemistryModule::VectorCMap &mapped_values) { - if (this->is_master) { - for (const std::string &prop_name : this->prop_names) { - const auto m_it = mapped_values.find(prop_name); - if (m_it == mapped_values.end()) { - throw std::domain_error( - "Prop names vector does not match any key in given map."); + std::vector essentials_backup{ + prop_names.begin() + speciesPerModule[0], prop_names.end()}; + + std::vector new_solution_names{ + this->prop_names.begin(), this->prop_names.begin() + speciesPerModule[0]}; + + if (is_master) { + for (const auto &init_val : input_map) { + std::string name = init_val.first; + if (std::find(new_solution_names.begin(), new_solution_names.end(), + name) == new_solution_names.end()) { + int size = name.size(); + ChemBCast(&size, 1, MPI_INT); + ChemBCast(name.data(), name.size(), MPI_CHAR); + new_solution_names.push_back(name); } - - const auto field_row = m_it->second; - const auto chem_field_end = this->field.end(); - this->field.insert(chem_field_end, field_row.begin(), field_row.end()); } - this->n_cells = mapped_values.begin()->second.size(); + int end = 0; + ChemBCast(&end, 1, MPI_INT); + } else { + constexpr int MAXSIZE = 128; + MPI_Status status; + int recv_size; + char recv_buffer[MAXSIZE]; + while (1) { + ChemBCast(&recv_size, 1, MPI_INT); + if (recv_size == 0) { + break; + } + ChemBCast(recv_buffer, recv_size, MPI_CHAR); + recv_buffer[recv_size] = '\0'; + new_solution_names.push_back(std::string(recv_buffer)); + } + } + + // now sort the new values + std::sort(new_solution_names.begin() + 3, new_solution_names.end()); + + // and append other processes than solutions + std::vector new_prop_names = new_solution_names; + new_prop_names.insert(new_prop_names.end(), essentials_backup.begin(), + essentials_backup.end()); + + std::vector old_prop_names{this->prop_names}; + this->prop_names = std::move(new_prop_names); + this->prop_count = prop_names.size(); + + this->SetPOETSolutionNames(new_solution_names); + + if (is_master) { + this->n_cells = n_cells; + + this->field.clear(); + this->field.reserve(this->prop_count * n_cells); + + std::vector init_values; + this->getDumpedField(init_values); + + const std::vector ess_names = old_prop_names; + + for (const auto &name : this->prop_names) { + auto it_find = std::find(ess_names.begin(), ess_names.end(), name); + double value; + if (it_find != ess_names.end()) { + int index = it_find - ess_names.begin(); + value = init_values[index]; + } else { + auto map_it = input_map.find(name); + value = map_it->second; + } + std::vector curr_vec(n_cells, value); + this->field.insert(field.end(), curr_vec.begin(), curr_vec.end()); + } } } @@ -274,89 +290,26 @@ void poet::ChemistryModule::ReadDHTFile(const std::string &input_file) { } } -void poet::ChemistryModule::GetWPFromInternals(std::vector &vecWP, - uint32_t wp_size) { - std::vector vecCurrOutput; - - vecWP.clear(); - vecWP.reserve(this->prop_count * wp_size); - - this->GetConcentrations(vecCurrOutput); - vecWP.insert(vecWP.end(), vecCurrOutput.begin(), vecCurrOutput.end()); - - if (this->speciesPerModule[1] != 0) { - this->GetPPhaseMoles(vecCurrOutput); - vecWP.insert(vecWP.end(), vecCurrOutput.begin(), vecCurrOutput.end()); - } - - // NOTE: Block for 'Surface' and 'Exchange' is missing because of missing - // Getters @ PhreeqcRM - // ... - // BLOCK_END - - if (this->speciesPerModule[4] != 0) { - this->GetKineticsMoles(vecCurrOutput); - vecWP.insert(vecWP.end(), vecCurrOutput.begin(), vecCurrOutput.end()); - } -} -void poet::ChemistryModule::SetInternalsFromWP(const std::vector &vecWP, - uint32_t wp_size) { - uint32_t iCurrElements; - - auto itStart = vecWP.begin(); - auto itEnd = itStart; - - // this->SetMappingForWP(iCurrWPSize); - - int nchem = this->GetChemistryCellCount(); - - iCurrElements = this->speciesPerModule[0]; - - itEnd += iCurrElements * wp_size; - this->SetConcentrations(std::vector(itStart, itEnd)); - itStart = itEnd; - - // Equlibirum Phases - if ((iCurrElements = this->speciesPerModule[1]) != 0) { - itEnd += iCurrElements * wp_size; - this->SetPPhaseMoles(std::vector(itStart, itEnd)); - itStart = itEnd; - } - - // // NOTE: Block for 'Surface' and 'Exchange' is missing because of missing - // // setters @ PhreeqcRM - // // ... - // // BLOCK_END - - if ((iCurrElements = this->speciesPerModule[4]) != 0) { - itEnd += iCurrElements * wp_size; - this->SetKineticsMoles(std::vector(itStart, itEnd)); - itStart = itEnd; - } -} - -#else //POET_USE_PRM +#else // POET_USE_PRM inline void poet::ChemistryModule::PrmToPoetField(std::vector &field) { - GetConcentrations(field); - int col = GetSelectedOutputColumnCount(); - int rows = GetSelectedOutputRowCount(); - - field.reserve(field.size() + 3 * rows); - - std::vector so; - GetSelectedOutput(so); - - for (int j = 0; j < col; j += 2) { - const auto start = so.begin() + (j * rows); - const auto end = start + rows; - field.insert(field.end(), start, end); - } + this->getDumpedField(field); } void poet::ChemistryModule::RunCells() { PhreeqcRM::RunCells(); - PrmToPoetField(this->field); + + std::vector tmp_field; + + PrmToPoetField(tmp_field); + this->field = tmp_field; + + for (uint32_t i = 0; i < field.size(); i++) { + uint32_t back_i = static_cast(i / this->nxyz); + uint32_t mod_i = i % this->nxyz; + + field[i] = tmp_field[back_i + (mod_i * this->prop_count)]; + } } #endif diff --git a/src/ChemistryModule/MasterFunctions.cpp b/src/ChemistryModule/MasterFunctions.cpp index d5667552f..4cd47d6d7 100644 --- a/src/ChemistryModule/MasterFunctions.cpp +++ b/src/ChemistryModule/MasterFunctions.cpp @@ -233,9 +233,12 @@ void poet::ChemistryModule::RunCells() { } void poet::ChemistryModule::MasterRunSequential() { - SetInternalsFromWP(this->field, this->nxyz); + std::vector shuffled_field = + shuffleField(field, n_cells, prop_count, 1); + this->setDumpedField(shuffled_field); PhreeqcRM::RunCells(); - GetWPFromInternals(this->field, this->nxyz); + this->getDumpedField(shuffled_field); + unshuffleField(shuffled_field, n_cells, prop_count, 1, this->field); } void poet::ChemistryModule::MasterRunParallel() { diff --git a/src/ChemistryModule/WorkerFunctions.cpp b/src/ChemistryModule/WorkerFunctions.cpp index f535cc982..5022b2f38 100644 --- a/src/ChemistryModule/WorkerFunctions.cpp +++ b/src/ChemistryModule/WorkerFunctions.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -42,6 +43,12 @@ void poet::ChemistryModule::WorkerLoop() { RunInitFile(get_string(0, this->group_comm)); break; } + case CHEM_INIT_SPECIES: { + SingleCMap dummy_map; + std::uint32_t n_cells_dummy; + mergeFieldWithModule(dummy_map, n_cells_dummy); + break; + } case CHEM_DHT_ENABLE: { bool enable; ChemBCast(&enable, 1, MPI_CXX_BOOL); @@ -339,32 +346,16 @@ poet::ChemistryModule::WorkerRunWorkPackage(std::vector &vecWP, return IRM_OK; } - std::vector vecCopy; - - vecCopy = vecWP; - for (uint32_t i = 0; i < this->prop_count; i++) { - for (uint32_t j = 0; j < this->wp_size; j++) { - vecWP[(i * this->wp_size) + j] = vecCopy[(j * this->prop_count) + i]; - } - } - IRM_RESULT result; this->PhreeqcRM::CreateMapping(vecMapping); - this->SetInternalsFromWP(vecWP, this->wp_size); + this->setDumpedField(vecWP); this->PhreeqcRM::SetTime(dSimTime); this->PhreeqcRM::SetTimeStep(dTimestep); result = this->PhreeqcRM::RunCells(); - this->GetWPFromInternals(vecWP, this->wp_size); - - vecCopy = vecWP; - for (uint32_t i = 0; i < this->prop_count; i++) { - for (uint32_t j = 0; j < this->wp_size; j++) { - vecWP[(j * this->prop_count) + i] = vecCopy[(i * this->wp_size) + j]; - } - } + this->getDumpedField(vecWP); return result; } diff --git a/src/DiffusionModule.cpp b/src/DiffusionModule.cpp index bca4728d1..f2dd35567 100644 --- a/src/DiffusionModule.cpp +++ b/src/DiffusionModule.cpp @@ -74,7 +74,8 @@ void DiffusionModule::initialize(poet::DiffusionParams args) { // const poet::DiffusionParams args(this->R); // name of props - this->prop_names = args.prop_names; + // + this->prop_names = Rcpp::as>(args.initial_t.names()); this->prop_count = this->prop_names.size(); this->state = @@ -91,7 +92,9 @@ void DiffusionModule::initialize(poet::DiffusionParams args) { // initial input this->alpha.push_back(args.alpha[this->prop_names[i]]); - std::vector prop_vec = grid.GetSpeciesByName(this->prop_names[i]); + double val = args.initial_t[prop_names[i]]; + + std::vector prop_vec(n_cells_per_prop, val); std::copy(prop_vec.begin(), prop_vec.end(), field.begin() + (i * this->n_cells_per_prop)); if (this->dim == this->DIM_2D) { diff --git a/src/Grid.cpp b/src/Grid.cpp index ba70e16ee..64561b5ae 100644 --- a/src/Grid.cpp +++ b/src/Grid.cpp @@ -21,7 +21,6 @@ #include "poet/SimParams.hpp" #include #include -#include #include #include #include diff --git a/src/SimParams.cpp b/src/SimParams.cpp index c7bb3b30a..847cf31d2 100644 --- a/src/SimParams.cpp +++ b/src/SimParams.cpp @@ -19,7 +19,6 @@ */ #include "poet/DHT_Types.hpp" -#include #include #include @@ -51,8 +50,8 @@ poet::GridParams::s_GridParams(RInside &R) { } poet::DiffusionParams::s_DiffusionParams(RInside &R) { - this->prop_names = Rcpp::as>( - R.parseEval("names(mysetup$diffusion$init)")); + this->initial_t = + Rcpp::as(R.parseEval("mysetup$diffusion$init")); this->alpha = Rcpp::as(R.parseEval("mysetup$diffusion$alpha")); if (Rcpp::as( From 1716382b840bf018a33a7908e36be4faf5a3a4fe Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 17 Apr 2023 12:34:10 +0200 Subject: [PATCH 2/2] refactor: remove bits/stdint-uintn.h as header --- include/poet/DiffusionModule.hpp | 1 - src/DiffusionModule.cpp | 1 - 2 files changed, 2 deletions(-) diff --git a/include/poet/DiffusionModule.hpp b/include/poet/DiffusionModule.hpp index 63c514189..e2ac32466 100644 --- a/include/poet/DiffusionModule.hpp +++ b/include/poet/DiffusionModule.hpp @@ -23,7 +23,6 @@ #include "poet/SimParams.hpp" #include -#include #include #include #include diff --git a/src/DiffusionModule.cpp b/src/DiffusionModule.cpp index f2dd35567..c3d37ee75 100644 --- a/src/DiffusionModule.cpp +++ b/src/DiffusionModule.cpp @@ -23,7 +23,6 @@ #include "tug/Diffusion.hpp" #include #include -#include #include #include #include