From 45ea77ae0ff5b798c224f4305d5df5f9d73bf2de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 15 Oct 2024 12:47:48 +0200 Subject: [PATCH 1/2] Update subproject commit in ext/iphreeqc --- R_lib/init_r_lib.R | 30 +++++- ext/iphreeqc | 2 +- src/Chemistry/ChemistryModule.cpp | 9 +- src/Init/ChemistryInit.cpp | 25 ++--- src/Init/DiffusionInit.cpp | 160 ++++++++++++++-------------- src/Init/GridInit.cpp | 167 +++++++++++++++++++----------- src/Init/InitialList.cpp | 69 ++++++------ src/Init/InitialList.hpp | 54 ++++------ 8 files changed, 292 insertions(+), 224 deletions(-) diff --git a/R_lib/init_r_lib.R b/R_lib/init_r_lib.R index 468b485bc..f3f1adeed 100644 --- a/R_lib/init_r_lib.R +++ b/R_lib/init_r_lib.R @@ -13,6 +13,34 @@ ### this program; if not, write to the Free Software Foundation, Inc., 51 ### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +##' @param pqc_mat matrix, containing IDs and PHREEQC outputs +##' @param grid matrix, zonation referring to pqc_mat$ID +##' @return a data.frame +# pqc_to_grid <- function(pqc_mat, grid) { +# # Convert the input DataFrame to a matrix +# pqc_mat <- as.matrix(pqc_mat) + +# # Flatten the matrix into a vector +# id_vector <- as.integer(t(grid)) + +# # Find the matching rows in the matrix +# row_indices <- match(id_vector, pqc_mat[, "ID"]) + +# # Extract the matching rows from pqc_mat to size of grid matrix +# result_mat <- pqc_mat[row_indices, ] + +# # Convert the result matrix to a data frame +# res_df <- as.data.frame(result_mat) + +# # Remove all columns which only contain NaN +# res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)] + +# # Remove row names +# rownames(res_df) <- NULL + +# return(res_df) +# } + ##' @param pqc_mat matrix, containing IDs and PHREEQC outputs ##' @param grid matrix, zonation referring to pqc_mat$ID ##' @return a data.frame @@ -33,7 +61,7 @@ pqc_to_grid <- function(pqc_mat, grid) { res_df <- as.data.frame(result_mat) # Remove all columns which only contain NaN - res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)] + # res_df <- res_df[, colSums(is.na(res_df)) != nrow(res_df)] # Remove row names rownames(res_df) <- NULL diff --git a/ext/iphreeqc b/ext/iphreeqc index e6e5e0d51..61214a01a 160000 --- a/ext/iphreeqc +++ b/ext/iphreeqc @@ -1 +1 @@ -Subproject commit e6e5e0d5156c093241a53e6ce074ef346d64ae26 +Subproject commit 61214a01ad4cf99527f657e6a41c68282e4886c4 diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp index f2edf5c28..a7462b0ec 100644 --- a/src/Chemistry/ChemistryModule.cpp +++ b/src/Chemistry/ChemistryModule.cpp @@ -1,6 +1,7 @@ #include "ChemistryModule.hpp" #include "PhreeqcEngine.hpp" +#include "PhreeqcMatrix.hpp" #include "SurrogateModels/DHT_Wrapper.hpp" #include "SurrogateModels/Interpolation.hpp" @@ -167,9 +168,11 @@ poet::ChemistryModule::ChemistryModule( this->n_cells = chem_params.total_grid_cells; if (!is_master) { - for (const auto &[pqc_id, pqc_config] : chem_params.pqc_config) { - this->phreeqc_instances[pqc_id] = - std::make_unique(pqc_config); + PhreeqcMatrix pqc_mat = + PhreeqcMatrix(chem_params.database, chem_params.pqc_script); + for (const auto &cell_id : chem_params.pqc_ids) { + this->phreeqc_instances[cell_id] = + std::make_unique(pqc_mat, cell_id); } // for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) { // this->phreeqc_instances[chem_params.pqc_ids[i]] = diff --git a/src/Init/ChemistryInit.cpp b/src/Init/ChemistryInit.cpp index 3c7a9871c..b399671ed 100644 --- a/src/Init/ChemistryInit.cpp +++ b/src/Init/ChemistryInit.cpp @@ -2,16 +2,14 @@ #include #include +#include +#include +#include #include namespace poet { void InitialList::initChemistry(const Rcpp::List &chem) { - this->pqc_solutions = std::vector( - this->transport_names.begin() + 3, this->transport_names.end()); - - this->pqc_solution_primaries = this->phreeqc->getSolutionPrimaries(); - if (chem.containsElementNamed("dht_species")) { this->dht_species = Rcpp::as>(chem["dht_species"]); } @@ -69,21 +67,14 @@ InitialList::ChemistryInit InitialList::getChemistryInit() const { // chem_init.field_header = this->field_header; chem_init.database = database; + chem_init.pqc_script = pqc_script; + chem_init.pqc_ids = pqc_ids; // chem_init.pqc_scripts = pqc_scripts; // chem_init.pqc_ids = pqc_ids; - for (std::size_t i = 0; i < pqc_scripts.size(); i++) { - POETInitCell cell = { - pqc_solutions, - pqc_solution_primaries, - Rcpp::as>(pqc_exchanger[i]), - Rcpp::as>(pqc_kinetics[i]), - Rcpp::as>(pqc_equilibrium[i]), - Rcpp::as>(pqc_surface_comps[i]), - Rcpp::as>(pqc_surface_charges[i])}; - - chem_init.pqc_config[pqc_ids[i]] = {database, pqc_scripts[i], cell}; - } + // for (std::size_t i = 0; i < pqc_ids.size(); ++i) { + // chem_init.pqc_input[pqc_ids[i]] = pqc_scripts[i]; + // } // chem_init.pqc_sol_order = pqc_solutions; diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp index 242929b47..7abd987ec 100644 --- a/src/Init/DiffusionInit.cpp +++ b/src/Init/DiffusionInit.cpp @@ -1,4 +1,3 @@ -#include #include // leave above Rcpp includes, as eigen seem to have problems with a preceding // Rcpp include @@ -8,7 +7,7 @@ #include "DataStructures/Field.hpp" #include "InitialList.hpp" -#include "PhreeqcInit.hpp" +#include "PhreeqcMatrix.hpp" #include #include @@ -53,98 +52,102 @@ static std::vector colMajToRowMaj(const Rcpp::NumericVector &vec, } } -static std::vector -extend_transport_names(std::unique_ptr &phreeqc, - const Rcpp::List &boundaries_list, - const Rcpp::List &inner_boundaries, - const std::vector &old_trans_names) { +// static std::vector +// extend_transport_names(std::unique_ptr &phreeqc, +// const Rcpp::List &boundaries_list, +// const Rcpp::List &inner_boundaries, +// const std::vector &old_trans_names) { - std::vector transport_names = old_trans_names; - std::set constant_pqc_ids; +// std::vector transport_names = old_trans_names; +// std::set constant_pqc_ids; - for (const auto &side : tug_side_mapping) { - if (!boundaries_list.containsElementNamed(side.second.c_str())) { - continue; - } +// for (const auto &side : tug_side_mapping) { +// if (!boundaries_list.containsElementNamed(side.second.c_str())) { +// continue; +// } - Rcpp::List mapping = boundaries_list[side.second]; +// Rcpp::List mapping = boundaries_list[side.second]; - const Rcpp::NumericVector cells = mapping["cell"]; - const Rcpp::NumericVector values = mapping["sol_id"]; - const Rcpp::CharacterVector type_str = mapping["type"]; +// const Rcpp::NumericVector cells = mapping["cell"]; +// const Rcpp::NumericVector values = mapping["sol_id"]; +// const Rcpp::CharacterVector type_str = mapping["type"]; - if (cells.size() != values.size()) { - throw std::runtime_error("Boundary [" + side.second + - "] cells and values are not the same " - "length"); - } +// if (cells.size() != values.size()) { +// throw std::runtime_error("Boundary [" + side.second + +// "] cells and values are not the same " +// "length"); +// } - for (auto i = 0; i < cells.size(); i++) { - if (type_str[i] == "constant") { - constant_pqc_ids.insert(static_cast(values[i])); - } - } - } +// for (auto i = 0; i < cells.size(); i++) { +// if (type_str[i] == "constant") { +// constant_pqc_ids.insert(static_cast(values[i])); +// } +// } +// } - if (inner_boundaries.size() > 0) { - const Rcpp::NumericVector values = inner_boundaries["sol_id"]; - for (auto i = 0; i < values.size(); i++) { - constant_pqc_ids.insert(static_cast(values[i])); - } - } +// if (inner_boundaries.size() > 0) { +// const Rcpp::NumericVector values = inner_boundaries["sol_id"]; +// for (auto i = 0; i < values.size(); i++) { +// constant_pqc_ids.insert(static_cast(values[i])); +// } +// } - if (!constant_pqc_ids.empty()) { - constexpr std::size_t keep_h_o_charge = 3; +// if (!constant_pqc_ids.empty()) { +// constexpr std::size_t keep_h_o_charge = 3; - for (const auto &pqc_id : constant_pqc_ids) { - const auto solution_names = phreeqc->getSolutionNames(pqc_id); +// for (const auto &pqc_id : constant_pqc_ids) { +// const auto solution_names = phreeqc->getSolutionNames(pqc_id); - // add those strings which are not already in the transport_names - for (const auto &name : solution_names) { - if (std::find(transport_names.begin(), transport_names.end(), name) == - transport_names.end()) { - auto position = - std::lower_bound(transport_names.begin() + keep_h_o_charge, - transport_names.end(), name); +// // add those strings which are not already in the transport_names +// for (const auto &name : solution_names) { +// if (std::find(transport_names.begin(), transport_names.end(), name) +// == +// transport_names.end()) { +// auto position = +// std::lower_bound(transport_names.begin() + keep_h_o_charge, +// transport_names.end(), name); - transport_names.insert(position, name); - } - } - } - } +// transport_names.insert(position, name); +// } +// } +// } +// } - return transport_names; -} +// return transport_names; +// } -static Rcpp::List -extend_initial_grid(const Rcpp::List &initial_grid, - const std::vector &transport_names) { - // std::vector names_to_add(transport_names.begin() + old_size, - // transport_names.end()); +// static Rcpp::List +// extend_initial_grid(const Rcpp::List &initial_grid, +// const std::vector &transport_names) { +// // std::vector names_to_add(transport_names.begin() + +// old_size, +// // transport_names.end()); - Rcpp::Function extend_grid_R("add_missing_transport_species"); +// Rcpp::Function extend_grid_R("add_missing_transport_species"); - return extend_grid_R(initial_grid, Rcpp::wrap(transport_names)); -} +// return extend_grid_R(initial_grid, Rcpp::wrap(transport_names)); +// } std::pair InitialList::resolveBoundaries(const Rcpp::List &boundaries_list, - const Rcpp::List &inner_boundaries) { + const Rcpp::List &inner_boundaries, + const PhreeqcMatrix &phreeqc_mat) { Rcpp::List bound_list; Rcpp::List inner_bound; Rcpp::Function resolve_R("resolve_pqc_bound"); - const std::size_t old_transport_size = this->transport_names.size(); + // const std::size_t old_transport_size = this->transport_names.size(); - this->transport_names = extend_transport_names( - this->phreeqc, boundaries_list, inner_boundaries, this->transport_names); + // this->transport_names = extend_transport_names( + // this->phreeqc, boundaries_list, inner_boundaries, + // this->transport_names); - const std::size_t new_transport_size = this->transport_names.size(); + // const std::size_t new_transport_size = this->transport_names.size(); - if (old_transport_size != new_transport_size) { - this->initial_grid = - extend_initial_grid(this->initial_grid, this->transport_names); - } + // if (old_transport_size != new_transport_size) { + // this->initial_grid = + // extend_initial_grid(this->initial_grid, this->transport_names); + // } for (const auto &species_name : this->transport_names) { Rcpp::List spec_list; @@ -179,8 +182,11 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list, c_type[c_id] = tug::BC_TYPE_CLOSED; } else if (type_str[i] == "constant") { c_type[c_id] = tug::BC_TYPE_CONSTANT; - c_value[c_id] = Rcpp::as( - resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name), values[i])); + c_value[c_id] = + static_cast(phreeqc_mat(values[i], species_name)); + // c_value[c_id] = Rcpp::as( + // resolve_R(this->phreeqc_mat, Rcpp::wrap(species_name), + // values[i])); } else { throw std::runtime_error("Unknown boundary type"); } @@ -211,8 +217,11 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list, for (std::size_t i = 0; i < inner_row_vec.size(); i++) { rows.push_back(inner_row_vec[i] - 1); cols.push_back(inner_col_vec[i] - 1); - c_value.push_back(Rcpp::as(resolve_R( - this->phreeqc_mat, Rcpp::wrap(species_name), inner_pqc_id_vec[i]))); + c_value.push_back(static_cast( + phreeqc_mat(inner_pqc_id_vec[i], species_name))); + // c_value.push_back(Rcpp::as(resolve_R( + // this->phreeqc_mat, Rcpp::wrap(species_name), + // inner_pqc_id_vec[i]))); } inner_bound[species_name] = @@ -276,7 +285,8 @@ static Rcpp::List parseAlphas(const SEXP &input, return out_list; } -void InitialList::initDiffusion(const Rcpp::List &diffusion_input) { +void InitialList::initDiffusion(const Rcpp::List &diffusion_input, + const PhreeqcMatrix &phreeqc) { Rcpp::List boundaries; Rcpp::List inner_boundaries; @@ -298,7 +308,7 @@ void InitialList::initDiffusion(const Rcpp::List &diffusion_input) { diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_Y)]; const auto resolved_boundaries = - resolveBoundaries(boundaries, inner_boundaries); + resolveBoundaries(boundaries, inner_boundaries, phreeqc); this->boundaries = resolved_boundaries.first; this->inner_boundaries = resolved_boundaries.second; diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp index 39305dbb2..1b8f24c7f 100644 --- a/src/Init/GridInit.cpp +++ b/src/Init/GridInit.cpp @@ -1,11 +1,13 @@ #include "InitialList.hpp" -#include "PhreeqcInit.hpp" +#include "PhreeqcMatrix.hpp" #include #include +#include #include #include +#include #include #include #include @@ -15,37 +17,41 @@ namespace poet { -static Rcpp::NumericMatrix -pqcScriptToGrid(std::unique_ptr &phreeqc, RInside &R) { - PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat(); +// static Rcpp::NumericMatrix pqcMatToR(const PhreeqcMatrix &phreeqc, RInside +// &R) { - // add "id" to the front of the column names +// PhreeqcMatrix::STLExport phreeqc_mat = phreeqc.get(); - const std::size_t ncols = phreeqc_mat.names.size(); - const std::size_t nrows = phreeqc_mat.values.size(); +// // PhreeqcInit::PhreeqcMat phreeqc_mat = phreeqc->getPhreeqcMat(); - phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID"); +// // add "id" to the front of the column names - Rcpp::NumericMatrix mat(nrows, ncols + 1); +// const std::size_t ncols = phreeqc_mat.names.size(); +// const std::size_t nrows = phreeqc_mat.values.size(); - 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]; - } - } +// phreeqc_mat.names.insert(phreeqc_mat.names.begin(), "ID"); - Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names); +// Rcpp::NumericMatrix mat(nrows, ncols + 1); - return mat; -} +// 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]; +// } +// } -static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat, - const Rcpp::NumericMatrix &grid) { - Rcpp::Function pqc_to_grid_R("pqc_to_grid"); +// Rcpp::colnames(mat) = Rcpp::wrap(phreeqc_mat.names); - return pqc_to_grid_R(mat, grid); -} +// 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) { @@ -60,24 +66,25 @@ 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"); +// 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 row_ids = +// Rcpp::as>(unique_R(initial_grid["ID"])); - // std::vector sizes_vec(sizes.begin(), sizes.end()); +// // 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)); +// // 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()); +// // std::copy(sizes_vec.begin(), sizes_vec.end(), sizes.begin()); - return phreeqc->getModuleSizes(row_ids)[POET_SOL]; -} +// return phreeqc->getModuleSizes(row_ids)[POET_SOL]; +// } static std::string readFile(const std::string &path) { std::string string_rpath(PATH_MAX, '\0'); @@ -98,10 +105,26 @@ static std::string readFile(const std::string &path) { return buffer.str(); } -void InitialList::prepareGrid(const Rcpp::List &grid_input) { - // parse input values - Rcpp::Function unique_R("unique"); +static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat, + const std::vector unique_ids, + const Rcpp::IntegerMatrix &grid_def) { + PhreeqcMatrix subset_pqc_mat = pqc_mat.subset(unique_ids); + + PhreeqcMatrix::STLExport phreeqc_mat = + subset_pqc_mat.get(PhreeqcMatrix::VectorExportType::COLUMN_MAJOR); + + const std::size_t ncols = phreeqc_mat.names.size(); + const std::size_t nrows = phreeqc_mat.values.size() / ncols; + + Rcpp::NumericMatrix pqc_mat_R(nrows, ncols, phreeqc_mat.values.begin()); + Rcpp::colnames(pqc_mat_R) = Rcpp::wrap(phreeqc_mat.names); + + return Rcpp::Function("pqc_to_grid")(pqc_mat_R, grid_def); +} + +PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) { + // parse input values std::string script; std::string database; @@ -126,8 +149,9 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) { } this->database = database; + this->pqc_script = script; - Rcpp::NumericMatrix grid_def = + Rcpp::IntegerMatrix grid_def = grid_input[GRID_MEMBER_STR(GridMembers::GRID_DEF)]; Rcpp::NumericVector grid_size = grid_input[GRID_MEMBER_STR(GridMembers::GRID_SIZE)]; @@ -156,35 +180,54 @@ void InitialList::prepareGrid(const Rcpp::List &grid_input) { throw std::runtime_error("Grid size must be positive."); } - this->phreeqc = std::make_unique(database, script); + PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script); - this->phreeqc_mat = pqcScriptToGrid(phreeqc, R); - this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def); + this->transport_names = pqc_mat.getSolutionNames(true); - const uint32_t solution_count = getSolutionCount(phreeqc, this->initial_grid); + Rcpp::Function unique_R("unique"); + Rcpp::Function sort_R("sort"); - std::vector colnames = - Rcpp::as>(this->initial_grid.names()); + std::vector unique_ids = Rcpp::as>( + unique_R(Rcpp::IntegerVector(grid_def.begin(), grid_def.end()))); - this->transport_names = std::vector( - colnames.begin() + 1, - colnames.begin() + 1 + solution_count); // skip ID + this->initial_grid = expandGrid(pqc_mat, unique_ids, grid_def); - std::map pqc_raw_dumps; + const auto pqc_raw_dumps = replaceRawKeywordIDs(pqc_mat.getDumpStringsPQI()); - pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps()); - - this->pqc_ids = - Rcpp::as>(unique_R(this->initial_grid["ID"])); - - for (const auto &id : this->pqc_ids) { - this->pqc_scripts.push_back(pqc_raw_dumps[id]); - this->pqc_exchanger.push_back(phreeqc->getExchanger(id)); - this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id)); - this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id)); - this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id)); - this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id)); + for (const auto &id : unique_ids) { + this->pqc_ids.push_back(id); } + + return pqc_mat; + // this->phreeqc_mat = pqcScriptToGrid(phreeqc, R); + // this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def); + + // const uint32_t solution_count = getSolutionCount(phreeqc, + // this->initial_grid); + + // std::vector colnames = + // Rcpp::as>(this->initial_grid.names()); + + // this->transport_names = std::vector( + // colnames.begin() + 1, + // colnames.begin() + 1 + solution_count); // skip ID + + // std::map pqc_raw_dumps; + + // pqc_raw_dumps = replaceRawKeywordIDs(phreeqc->raw_dumps()); + + // this->pqc_ids = + // Rcpp::as>(unique_R(this->initial_grid["ID"])); + + // for (const auto &id : this->pqc_ids) { + // this->pqc_scripts.push_back(pqc_raw_dumps[id]); + // // this->pqc_exchanger.push_back(phreeqc->getExchanger(id)); + // // this->pqc_kinetics.push_back(phreeqc->getKineticsNames(id)); + // // this->pqc_equilibrium.push_back(phreeqc->getEquilibriumNames(id)); + // // this->pqc_surface_comps.push_back(phreeqc->getSurfaceCompNames(id)); + // // + // this->pqc_surface_charges.push_back(phreeqc->getSurfaceChargeNames(id)); + // } } } // namespace poet \ No newline at end of file diff --git a/src/Init/InitialList.cpp b/src/Init/InitialList.cpp index 1aced099e..f4e87af14 100644 --- a/src/Init/InitialList.cpp +++ b/src/Init/InitialList.cpp @@ -1,5 +1,6 @@ #include "InitialList.hpp" #include "DataStructures/NamedVector.hpp" +#include "PhreeqcMatrix.hpp" #include #include #include @@ -10,8 +11,8 @@ namespace poet { void InitialList::initializeFromList(const Rcpp::List &setup) { - prepareGrid(setup[grid_key]); - initDiffusion(setup[diffusion_key]); + PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]); + initDiffusion(setup[diffusion_key], phreeqc); initChemistry(setup[chemistry_key]); } @@ -53,26 +54,26 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) { } this->database = Rcpp::as(setup[static_cast(ExportList::CHEM_DATABASE)]); + this->pqc_script = Rcpp::as( + setup[static_cast(ExportList::CHEM_PQC_SCRIPT)]); this->field_header = Rcpp::as>( setup[static_cast(ExportList::CHEM_FIELD_HEADER)]); - this->pqc_scripts = Rcpp::as>( - setup[static_cast(ExportList::CHEM_PQC_SCRIPTS)]); this->pqc_ids = Rcpp::as>( setup[static_cast(ExportList::CHEM_PQC_IDS)]); - this->pqc_solutions = Rcpp::as>( - setup[static_cast(ExportList::CHEM_PQC_SOLUTIONS)]); - this->pqc_solution_primaries = Rcpp::as>( - setup[static_cast(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]); - this->pqc_exchanger = - Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_EXCHANGER)]); - this->pqc_kinetics = - Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_KINETICS)]); - this->pqc_equilibrium = - Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_EQUILIBRIUM)]); - this->pqc_surface_comps = - Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_SURFACE_COMPS)]); - this->pqc_surface_charges = - Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_SURFACE_CHARGES)]); + // this->pqc_solutions = Rcpp::as>( + // setup[static_cast(ExportList::CHEM_PQC_SOLUTIONS)]); + // this->pqc_solution_primaries = Rcpp::as>( + // setup[static_cast(ExportList::CHEM_PQC_SOLUTION_PRIMARY)]); + // this->pqc_exchanger = + // Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_EXCHANGER)]); + // this->pqc_kinetics = + // Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_KINETICS)]); + // this->pqc_equilibrium = + // Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_EQUILIBRIUM)]); + // this->pqc_surface_comps = + // Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_SURFACE_COMPS)]); + // this->pqc_surface_charges = + // Rcpp::List(setup[static_cast(ExportList::CHEM_PQC_SURFACE_CHARGES)]); this->dht_species = NamedVector( setup[static_cast(ExportList::CHEM_DHT_SPECIES)]); @@ -108,25 +109,25 @@ Rcpp::List InitialList::exportList() { out[static_cast(ExportList::DIFFU_ALPHA_Y)] = this->alpha_y; 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_FIELD_HEADER)] = Rcpp::wrap(this->field_header); - out[static_cast(ExportList::CHEM_PQC_SCRIPTS)] = - Rcpp::wrap(this->pqc_scripts); out[static_cast(ExportList::CHEM_PQC_IDS)] = Rcpp::wrap(this->pqc_ids); - out[static_cast(ExportList::CHEM_PQC_SOLUTIONS)] = - Rcpp::wrap(this->pqc_solutions); - out[static_cast(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] = - Rcpp::wrap(this->pqc_solution_primaries); - out[static_cast(ExportList::CHEM_PQC_EXCHANGER)] = - Rcpp::wrap(this->pqc_exchanger); - out[static_cast(ExportList::CHEM_PQC_KINETICS)] = - Rcpp::wrap(this->pqc_kinetics); - out[static_cast(ExportList::CHEM_PQC_EQUILIBRIUM)] = - Rcpp::wrap(this->pqc_equilibrium); - out[static_cast(ExportList::CHEM_PQC_SURFACE_COMPS)] = - Rcpp::wrap(this->pqc_surface_comps); - out[static_cast(ExportList::CHEM_PQC_SURFACE_CHARGES)] = - Rcpp::wrap(this->pqc_surface_charges); + // out[static_cast(ExportList::CHEM_PQC_SOLUTIONS)] = + // Rcpp::wrap(this->pqc_solutions); + // out[static_cast(ExportList::CHEM_PQC_SOLUTION_PRIMARY)] = + // Rcpp::wrap(this->pqc_solution_primaries); + // out[static_cast(ExportList::CHEM_PQC_EXCHANGER)] = + // Rcpp::wrap(this->pqc_exchanger); + // out[static_cast(ExportList::CHEM_PQC_KINETICS)] = + // Rcpp::wrap(this->pqc_kinetics); + // out[static_cast(ExportList::CHEM_PQC_EQUILIBRIUM)] = + // Rcpp::wrap(this->pqc_equilibrium); + // out[static_cast(ExportList::CHEM_PQC_SURFACE_COMPS)] = + // Rcpp::wrap(this->pqc_surface_comps); + // out[static_cast(ExportList::CHEM_PQC_SURFACE_CHARGES)] = + // Rcpp::wrap(this->pqc_surface_charges); out[static_cast(ExportList::CHEM_DHT_SPECIES)] = this->dht_species; out[static_cast(ExportList::CHEM_INTERP_SPECIES)] = Rcpp::wrap(this->interp_species); diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp index 3a3c5ea23..a47d7918f 100644 --- a/src/Init/InitialList.hpp +++ b/src/Init/InitialList.hpp @@ -2,7 +2,7 @@ #include "Base/RInsidePOET.hpp" #include "DataStructures/NamedVector.hpp" -#include "POETInit.hpp" +#include "PhreeqcMatrix.hpp" #include #include #include @@ -13,7 +13,6 @@ #include #include -#include #include #include #include @@ -53,22 +52,23 @@ private: DIFFU_ALPHA_X, DIFFU_ALPHA_Y, CHEM_DATABASE, - CHEM_FIELD_HEADER, - CHEM_PQC_SCRIPTS, + CHEM_PQC_SCRIPT, CHEM_PQC_IDS, - 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_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, AI_SURROGATE_INPUT_SCRIPT, ENUM_SIZE // Hack: Last element of the enum to show enum size - }; + }; // Grid members static constexpr const char *grid_key = "Grid"; @@ -97,9 +97,9 @@ private: return GridMembersString[static_cast(member)]; } - std::unique_ptr phreeqc; + // std::unique_ptr pqc_mat; - void prepareGrid(const Rcpp::List &grid_input); + PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input); std::uint8_t dim{0}; @@ -114,9 +114,6 @@ private: Rcpp::List initial_grid; - // No export - Rcpp::NumericMatrix phreeqc_mat; - public: struct DiffusionInit { using BoundaryMap = @@ -169,10 +166,12 @@ private: return DiffusionMembersString[static_cast(member)]; } - void initDiffusion(const Rcpp::List &diffusion_input); + void initDiffusion(const Rcpp::List &diffusion_input, + const PhreeqcMatrix &phreeqc); std::pair resolveBoundaries(const Rcpp::List &boundaries_list, - const Rcpp::List &inner_boundaries); + const Rcpp::List &inner_boundaries, + const PhreeqcMatrix &phreeqc_mat); Rcpp::List boundaries; Rcpp::List inner_boundaries; @@ -190,17 +189,10 @@ private: std::vector field_header; std::string database; - std::vector pqc_scripts; + std::string pqc_script; + // std::vector pqc_scripts; std::vector pqc_ids; - std::vector pqc_solutions; - std::vector pqc_solution_primaries; - Rcpp::List pqc_exchanger; - Rcpp::List pqc_kinetics; - Rcpp::List pqc_equilibrium; - Rcpp::List pqc_surface_comps; - Rcpp::List pqc_surface_charges; - NamedVector dht_species; NamedVector interp_species; @@ -227,10 +219,10 @@ public: // std::vector field_header; std::string database; - // std::vector pqc_scripts; - // std::vector pqc_ids; + std::string pqc_script; + std::vector pqc_ids; - std::map pqc_config; + // std::map pqc_input; // std::vector pqc_sol_order; From cdbf344329c2e01291130ae5a12a41e0fbdefddc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 7 Nov 2024 14:16:53 +0100 Subject: [PATCH 2/2] fix: Update getSolutionNames call to remove unnecessary argument --- src/Init/GridInit.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Init/GridInit.cpp b/src/Init/GridInit.cpp index 1b8f24c7f..4a0e3fca3 100644 --- a/src/Init/GridInit.cpp +++ b/src/Init/GridInit.cpp @@ -182,7 +182,7 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) { PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script); - this->transport_names = pqc_mat.getSolutionNames(true); + this->transport_names = pqc_mat.getSolutionNames(); Rcpp::Function unique_R("unique"); Rcpp::Function sort_R("sort");